DualTraverse.cc 6.19 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * 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.
 * 
 ******************************************************************************/
20
21


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

namespace AMDiS {

28
29
30
31
32
33
  bool DualTraverse::traverseFirst(Mesh *mesh1, 
				   Mesh *mesh2, 
				   int level1, 
				   int level2, 
				   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

114
115
116
117
118
  bool DualTraverse::traverseNext(ElInfo **elInfo1,
				  ElInfo **elInfo2,
				  ElInfo **elInfoSmall,
				  ElInfo **elInfoLarge)
  {
119
    FUNCNAME("DualTraverse::traverseNext()");
120
    // call standard traverse
121
122
123
    if (inc1) {
      do {
	*elInfo1 = stack1.traverseNext(*elInfo1);
124
      } while(*elInfo1 != NULL && skipEl1(*elInfo1));
125
126
127
128
    }
    if (inc2) {
      do {
	*elInfo2 = stack2.traverseNext(*elInfo2);
129
      } while (*elInfo2 != NULL && skipEl2(*elInfo2));
130
    }
Naumann, Andreas's avatar
Naumann, Andreas committed
131
132

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

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

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

    // check for non domain covering level traverse
147
148
149
150
151
152
153
    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);
    }
154
155
156

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

157
158
159
    return true;
  }

160

161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
  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
177
    rest -= 1.0 / (1 << ((*elInfoSmall)->getLevel() - (*elInfoLarge)->getLevel()));
178

179
    if (rest < 1e-32) {
180
      // large element finished -> increment both elements
181
182
183
      rest = 1.0; 
      inc1 = true;
      inc2 = true;
184
185
    } else {
      // increment only small element
186
187
      inc1 = (*elInfo1 == *elInfoSmall) ? true : false;
      inc2 = (*elInfo2 == *elInfoSmall) ? true : false;
188
189
190
    }
  }

191

192
193
194
195
196
  void DualTraverse::fillSubElInfo(ElInfo *elInfo1, 
				   ElInfo *elInfo2,
				   ElInfo *elInfoSmall,
				   ElInfo *elInfoLarge)
  {
197
    if (!fillSubElemMat)
198
199
      return;

200
    if (elInfo1 == elInfoSmall)
201
      stack1.fillRefinementPath(*elInfoSmall, *elInfo2);
202
203
204
205
206
207
208
    else
      stack2.fillRefinementPath(*elInfoSmall, *elInfo1);    
  }


  int DualTraverse::getFace(DualElInfo *dualElInfo, int largeFace)
  {
209
    FUNCNAME_DBG("DualTraverse::getFace()");
210
211
212
213
214
215
216
217
218
219
220
221

    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;

222
    TEST_EXIT_DBG(smallElInfo->getRefinementPathLength() != 0)
223
224
225
      ("Refinement path is empty. This should not happen!\n");

    return -1;
226
  }
227
}