DualTraverse.cc 5.86 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
14
15
16
17
18
#include "DualTraverse.h"
#include "Mesh.h"
#include "ElInfo.h"

namespace AMDiS {

19
20
21
22
23
24
  bool DualTraverse::traverseFirst(Mesh *mesh1, 
				   Mesh *mesh2, 
				   int level1, 
				   int level2, 
				   Flag flag1,
				   Flag flag2,
25
26
27
28
29
30
				   ElInfo **elInfo1,
				   ElInfo **elInfo2,
				   ElInfo **elInfoSmall,
				   ElInfo **elInfoLarge)
  {
    // replace CALL_EL_LEVEL by CALL_MG_LEVEL (covers whole domain)
31
    if (flag1.isSet(Mesh::CALL_EL_LEVEL)) {
32
33
34
35
36
37
38
      flag1 &= ~Mesh::CALL_EL_LEVEL;
      flag1 |= Mesh::CALL_MG_LEVEL;
      level1_ = level1;
    } else {
      level1_ = -1;
    }

39
    if (flag2.isSet(Mesh::CALL_EL_LEVEL)) {
40
41
42
43
44
45
46
47
      flag2 &= ~Mesh::CALL_EL_LEVEL;
      flag2 |= Mesh::CALL_MG_LEVEL;
      level2_ = level2;
    } else {
      level2_ = -1;
    }

    // replace CALL_LEAF_EL_LEVEL by CALL_MG_LEVEL (covers whole domain)
48
    if (flag1.isSet(Mesh::CALL_LEAF_EL_LEVEL)) {
49
50
51
52
53
54
55
56
57
      flag1 &= ~Mesh::CALL_LEAF_EL_LEVEL;
      flag1 |= Mesh::CALL_MG_LEVEL;
      level1_ = level1;
      callLeafElLevel1_ = true;
    } else {
      level1_ = -1;
      callLeafElLevel1_ = false;
    }

58
    if (flag2.isSet(Mesh::CALL_LEAF_EL_LEVEL)) {
59
60
61
62
63
64
65
66
67
68
69
      flag2 &= ~Mesh::CALL_LEAF_EL_LEVEL;
      flag2 |= Mesh::CALL_MG_LEVEL;
      level2_ = level2;
      callLeafElLevel2_ = true;
    } else {
      level2_ = -1;
      callLeafElLevel2_ = false;
    }

    // call standard traverse
    *elInfo1 = stack1.traverseFirst(mesh1, level1, flag1);
Naumann, Andreas's avatar
Naumann, Andreas committed
70
    while (*elInfo1 != NULL && skipEl1(*elInfo1)) {
71
      *elInfo1 = stack1.traverseNext(*elInfo1);
Naumann, Andreas's avatar
Naumann, Andreas committed
72
    }
73

74
    *elInfo2 = stack2.traverseFirst(mesh2, level2, flag2);
Naumann, Andreas's avatar
Naumann, Andreas committed
75
    while (*elInfo2 != NULL && skipEl2(*elInfo2)) {
76
      *elInfo2 = stack2.traverseNext(*elInfo2);
Naumann, Andreas's avatar
Naumann, Andreas committed
77
    }
78
 
79
    // finished ?
80
    if (*elInfo1 == NULL || *elInfo2 == NULL) {
81
      TEST_EXIT(*elInfo1 == *elInfo2)("invalid dual traverse\n");
82
83
84
      return false;
    }

85
    rest = 1.0;
86

87
    bool accepted = check(elInfo1, elInfo2, elInfoSmall, elInfoLarge);
88
89

    // check for non domain covering level traverse
90
91
92
93
94
95
96
    if (!accepted ||
	(level1_ != -1 && (*elInfo1)->getLevel() != level1_) ||
	(callLeafElLevel1_ && !(*elInfo1)->getElement()->isLeaf()) ||
	(level2_ != -1 && (*elInfo2)->getLevel() != level2_) ||
	(callLeafElLevel2_ && !(*elInfo2)->getElement()->isLeaf())) {
      return traverseNext(elInfo1, elInfo2, elInfoSmall, elInfoLarge);
    }
97

98
99
    fillSubElInfo(*elInfo1, *elInfo2, *elInfoSmall, *elInfoLarge);

100
101
102
    return true;
  }

103

104
105
106
107
108
109
  bool DualTraverse::traverseNext(ElInfo **elInfo1,
				  ElInfo **elInfo2,
				  ElInfo **elInfoSmall,
				  ElInfo **elInfoLarge)
  {
    // call standard traverse
110
111
112
113
114
115
116
117
118
119
    if (inc1) {
      do {
	*elInfo1 = stack1.traverseNext(*elInfo1);
      } while(*elInfo1 != NULL && skipEl1(*elInfo1));
    }
    if (inc2) {
      do {
	*elInfo2 = stack2.traverseNext(*elInfo2);
      } while (*elInfo2 != NULL && skipEl2(*elInfo2));
    }
Naumann, Andreas's avatar
Naumann, Andreas committed
120
121

    // finished ?
122
123
124
125
    if (*elInfo1 == NULL || *elInfo2 == NULL) {
      TEST_EXIT(*elInfo1 == *elInfo2)("invalid dual traverse\n");
      return false;
    }
126
127

    // finished ?
128
    if (*elInfo1 == NULL || *elInfo2 == NULL) {
129
      TEST_EXIT(*elInfo1 == *elInfo2)("invalid dual traverse\n");
130
131
132
      return false;
    }

133
    bool accepted = check(elInfo1, elInfo2, elInfoSmall, elInfoLarge);
134
135

    // check for non domain covering level traverse
136
137
138
139
140
141
142
    if (!accepted ||
	(level1_ != -1 && (*elInfo1)->getLevel() != level1_) ||
	(callLeafElLevel1_ && !(*elInfo1)->getElement()->isLeaf()) ||
	(level2_ != -1 && (*elInfo2)->getLevel() != level2_) ||
	(callLeafElLevel2_ && !(*elInfo2)->getElement()->isLeaf())) {
      return traverseNext(elInfo1, elInfo2, elInfoSmall, elInfoLarge);
    }
143
144
145

    fillSubElInfo(*elInfo1, *elInfo2, *elInfoSmall, *elInfoLarge);

146
147
148
    return true;
  }

149

150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
  void DualTraverse::prepareNextStep(ElInfo **elInfo1,
				     ElInfo **elInfo2,
				     ElInfo **elInfoSmall,
				     ElInfo **elInfoLarge)
  {
    // which is the smaller element ?
    *elInfoSmall = 
      (*elInfo1)->getLevel() > (*elInfo2)->getLevel() ?
      *elInfo1 :
      *elInfo2;
    *elInfoLarge = 
      (*elInfo1)->getLevel() <= (*elInfo2)->getLevel() ?
      *elInfo1 :
      *elInfo2;

    // update rest
166
    rest -= 1.0 / (1 << ((*elInfoSmall)->getLevel() - (*elInfoLarge)->getLevel()));
167

168
    if (rest < 1e-32) {
169
      // large element finished -> increment both elements
170
171
172
      rest = 1.0; 
      inc1 = true;
      inc2 = true;
173
174
    } else {
      // increment only small element
175
176
      inc1 = (*elInfo1 == *elInfoSmall) ? true : false;
      inc2 = (*elInfo2 == *elInfoSmall) ? true : false;
177
178
179
    }
  }

180

181
182
183
184
185
  void DualTraverse::fillSubElInfo(ElInfo *elInfo1, 
				   ElInfo *elInfo2,
				   ElInfo *elInfoSmall,
				   ElInfo *elInfoLarge)
  {
186
    if (!fillSubElemMat)
187
188
      return;

189
    if (elInfo1 == elInfoSmall)
190
      stack1.fillRefinementPath(*elInfoSmall, *elInfo2);
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
    else
      stack2.fillRefinementPath(*elInfoSmall, *elInfo1);    
  }


  int DualTraverse::getFace(DualElInfo *dualElInfo, int largeFace)
  {
    FUNCNAME("DualTraverse::getFace()");

    TEST_EXIT_DBG(dualElInfo)("No dual element info object!\n");

    ElInfo *largeElInfo = dualElInfo->largeElInfo;
    ElInfo *smallElInfo = dualElInfo->smallElInfo;

    TEST_EXIT_DBG(largeElInfo->getLevel() <= smallElInfo->getLevel())
      ("Should not happen!\n");

    if (largeElInfo->getLevel() == smallElInfo->getLevel())
      return largeFace;

    int refinementPathLength = smallElInfo->getRefinementPathLength();
    unsigned long refinementPath = smallElInfo->getRefinementPath();

    TEST_EXIT_DBG(refinementPathLength != 0)
      ("Refinement path is empty. This should not happen!\n");

    std::cout << "REFINEMENTPATHLENGTH = " << refinementPathLength << std::endl;

    return -1;
220
  }
221
}