DualTraverse.cc 8.62 KB
Newer Older
1
2
3
4
5
6
7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9
10
11
12
13
14
15
16
17
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
18
 *
19
 ******************************************************************************/
20
21


22
23
24
25
26
27
#include "DualTraverse.h"
#include "Mesh.h"
#include "ElInfo.h"

namespace AMDiS {

28
29
30
31
  bool DualTraverse::traverseFirst(Mesh *mesh1,
				   Mesh *mesh2,
				   int level1,
				   int level2,
32
33
				   Flag flag1,
				   Flag flag2,
34
35
36
37
38
				   ElInfo **elInfo1,
				   ElInfo **elInfo2,
				   ElInfo **elInfoSmall,
				   ElInfo **elInfoLarge)
  {
39
    FUNCNAME("DualTraverse::traverseFirst()");
40
    // replace CALL_EL_LEVEL by CALL_MG_LEVEL (covers whole domain)
41
    if (flag1.isSet(Mesh::CALL_EL_LEVEL)) {
42
43
44
45
46
47
48
      flag1 &= ~Mesh::CALL_EL_LEVEL;
      flag1 |= Mesh::CALL_MG_LEVEL;
      level1_ = level1;
    } else {
      level1_ = -1;
    }

49
    if (flag2.isSet(Mesh::CALL_EL_LEVEL)) {
50
51
52
53
54
55
56
57
      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)
58
    if (flag1.isSet(Mesh::CALL_LEAF_EL_LEVEL)) {
59
60
61
62
63
64
65
66
67
      flag1 &= ~Mesh::CALL_LEAF_EL_LEVEL;
      flag1 |= Mesh::CALL_MG_LEVEL;
      level1_ = level1;
      callLeafElLevel1_ = true;
    } else {
      level1_ = -1;
      callLeafElLevel1_ = false;
    }

68
    if (flag2.isSet(Mesh::CALL_LEAF_EL_LEVEL)) {
69
70
71
72
73
74
75
76
77
78
79
      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);
80
    while (*elInfo1 != NULL && skipEl1(*elInfo1)) {
81
      *elInfo1 = stack1.traverseNext(*elInfo1);
Naumann, Andreas's avatar
Naumann, Andreas committed
82
    }
83

84
    *elInfo2 = stack2.traverseFirst(mesh2, level2, flag2);
85
    while (*elInfo2 != NULL && skipEl2(*elInfo2)) {
86
      *elInfo2 = stack2.traverseNext(*elInfo2);
Naumann, Andreas's avatar
Naumann, Andreas committed
87
    }
88

89
    // finished ?
90
    if (*elInfo1 == NULL || *elInfo2 == NULL) {
91
      TEST_EXIT(*elInfo1 == *elInfo2)("invalid dual traverse\n");
92
93
94
      return false;
    }

95
    rest = 1.0;
96

97
    bool accepted = check(elInfo1, elInfo2, elInfoSmall, elInfoLarge);
98
99

    // check for non domain covering level traverse
100
101
102
103
104
105
106
    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);
    }
107

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

110
111
112
    return true;
  }

113

Praetorius, Simon's avatar
Praetorius, Simon committed
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
  bool DualTraverse::traverseFirstOneMacro(Mesh *mesh1,
				   Mesh *mesh2,
                                   int macroIndex,
				   int level1,
				   int level2,
				   Flag flag1,
				   Flag flag2,
				   ElInfo **elInfo1,
				   ElInfo **elInfo2,
				   ElInfo **elInfoSmall,
				   ElInfo **elInfoLarge)
  {
    FUNCNAME("DualTraverse::traverseFirst()");
    // replace CALL_EL_LEVEL by CALL_MG_LEVEL (covers whole domain)
    if (flag1.isSet(Mesh::CALL_EL_LEVEL)) {
      flag1 &= ~Mesh::CALL_EL_LEVEL;
      flag1 |= Mesh::CALL_MG_LEVEL;
      level1_ = level1;
    } else {
      level1_ = -1;
    }

    if (flag2.isSet(Mesh::CALL_EL_LEVEL)) {
      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)
    if (flag1.isSet(Mesh::CALL_LEAF_EL_LEVEL)) {
      flag1 &= ~Mesh::CALL_LEAF_EL_LEVEL;
      flag1 |= Mesh::CALL_MG_LEVEL;
      level1_ = level1;
      callLeafElLevel1_ = true;
    } else {
      level1_ = -1;
      callLeafElLevel1_ = false;
    }

    if (flag2.isSet(Mesh::CALL_LEAF_EL_LEVEL)) {
      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.traverseFirstOneMacro(mesh1, macroIndex, level1, flag1);
    while (*elInfo1 != NULL && skipEl1(*elInfo1)) {
      *elInfo1 = stack1.traverseNext(*elInfo1);
    }

    *elInfo2 = stack2.traverseFirstOneMacro(mesh2, macroIndex, level2, flag2);
    while (*elInfo2 != NULL && skipEl2(*elInfo2)) {
      *elInfo2 = stack2.traverseNext(*elInfo2);
    }

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

    rest = 1.0;

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

    // check for non domain covering level traverse
    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);
    }

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

    return true;
  }


201
202
203
204
205
  bool DualTraverse::traverseNext(ElInfo **elInfo1,
				  ElInfo **elInfo2,
				  ElInfo **elInfoSmall,
				  ElInfo **elInfoLarge)
  {
206
    FUNCNAME("DualTraverse::traverseNext()");
207
    // call standard traverse
208
209
210
    if (inc1) {
      do {
	*elInfo1 = stack1.traverseNext(*elInfo1);
211
      } while(*elInfo1 != NULL && skipEl1(*elInfo1));
212
213
214
215
    }
    if (inc2) {
      do {
	*elInfo2 = stack2.traverseNext(*elInfo2);
216
      } while (*elInfo2 != NULL && skipEl2(*elInfo2));
217
    }
Naumann, Andreas's avatar
Naumann, Andreas committed
218
219

    // finished ?
220
    if (*elInfo1 == NULL || *elInfo2 == NULL) {
221
222
223
      TEST_EXIT(*elInfo1 == *elInfo2)("invalid dual traverse\n");
      return false;
    }
224
225

    // finished ?
226
    if (*elInfo1 == NULL || *elInfo2 == NULL) {
227
      TEST_EXIT(*elInfo1 == *elInfo2)("invalid dual traverse\n");
228
229
230
      return false;
    }

231
    bool accepted = check(elInfo1, elInfo2, elInfoSmall, elInfoLarge);
232
233

    // check for non domain covering level traverse
234
235
236
237
238
239
240
    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);
    }
241
242
243

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

244
245
246
    return true;
  }

247

248
249
250
251
252
253
  void DualTraverse::prepareNextStep(ElInfo **elInfo1,
				     ElInfo **elInfo2,
				     ElInfo **elInfoSmall,
				     ElInfo **elInfoLarge)
  {
    // which is the smaller element ?
254
    *elInfoSmall =
255
256
257
      (*elInfo1)->getLevel() > (*elInfo2)->getLevel() ?
      *elInfo1 :
      *elInfo2;
258
    *elInfoLarge =
259
260
261
262
263
      (*elInfo1)->getLevel() <= (*elInfo2)->getLevel() ?
      *elInfo1 :
      *elInfo2;

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

266
    if (rest < 1e-32) {
267
      // large element finished -> increment both elements
268
      rest = 1.0;
269
270
      inc1 = true;
      inc2 = true;
271
272
    } else {
      // increment only small element
273
274
      inc1 = (*elInfo1 == *elInfoSmall) ? true : false;
      inc2 = (*elInfo2 == *elInfoSmall) ? true : false;
275
276
277
    }
  }

278

279
  void DualTraverse::fillSubElInfo(ElInfo *elInfo1,
280
281
282
283
				   ElInfo *elInfo2,
				   ElInfo *elInfoSmall,
				   ElInfo *elInfoLarge)
  {
284
    if (!fillSubElemMat)
285
286
      return;

287
    if (elInfo1 == elInfoSmall)
288
      stack1.fillRefinementPath(*elInfoSmall, *elInfo2);
289
    else
290
      stack2.fillRefinementPath(*elInfoSmall, *elInfo1);
291
292
293
294
295
  }


  int DualTraverse::getFace(DualElInfo *dualElInfo, int largeFace)
  {
296
    FUNCNAME_DBG("DualTraverse::getFace()");
297
298
299
300
301
302
303
304
305
306
307
308

    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;

309
    TEST_EXIT_DBG(smallElInfo->getRefinementPathLength() != 0)
310
311
312
      ("Refinement path is empty. This should not happen!\n");

    return -1;
313
  }
314
}