DualTraverse.cc 5.61 KB
Newer Older
1
2
3
4
5
6
#include "DualTraverse.h"
#include "Mesh.h"
#include "ElInfo.h"

namespace AMDiS {

7
8
9
10
11
12
  bool DualTraverse::traverseFirst(Mesh *mesh1, 
				   Mesh *mesh2, 
				   int level1, 
				   int level2, 
				   Flag flag1,
				   Flag flag2,
13
14
15
16
17
18
				   ElInfo **elInfo1,
				   ElInfo **elInfo2,
				   ElInfo **elInfoSmall,
				   ElInfo **elInfoLarge)
  {
    // replace CALL_EL_LEVEL by CALL_MG_LEVEL (covers whole domain)
19
    if (flag1.isSet(Mesh::CALL_EL_LEVEL)) {
20
21
22
23
24
25
26
      flag1 &= ~Mesh::CALL_EL_LEVEL;
      flag1 |= Mesh::CALL_MG_LEVEL;
      level1_ = level1;
    } else {
      level1_ = -1;
    }

27
    if (flag2.isSet(Mesh::CALL_EL_LEVEL)) {
28
29
30
31
32
33
34
35
      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)
36
    if (flag1.isSet(Mesh::CALL_LEAF_EL_LEVEL)) {
37
38
39
40
41
42
43
44
45
      flag1 &= ~Mesh::CALL_LEAF_EL_LEVEL;
      flag1 |= Mesh::CALL_MG_LEVEL;
      level1_ = level1;
      callLeafElLevel1_ = true;
    } else {
      level1_ = -1;
      callLeafElLevel1_ = false;
    }

46
    if (flag2.isSet(Mesh::CALL_LEAF_EL_LEVEL)) {
47
48
49
50
51
52
53
54
55
56
57
      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
58
    while (*elInfo1 != NULL && skipEl1(*elInfo1)) {
59
      *elInfo1 = stack1.traverseNext(*elInfo1);
Naumann, Andreas's avatar
Naumann, Andreas committed
60
    }
61

62
    *elInfo2 = stack2.traverseFirst(mesh2, level2, flag2);
Naumann, Andreas's avatar
Naumann, Andreas committed
63
    while (*elInfo2 != NULL && skipEl2(*elInfo2)) {
64
      *elInfo2 = stack2.traverseNext(*elInfo2);
Naumann, Andreas's avatar
Naumann, Andreas committed
65
    }
66
 
67
    // finished ?
68
    if (*elInfo1 == NULL || *elInfo2 == NULL) {
69
      TEST_EXIT(*elInfo1 == *elInfo2)("invalid dual traverse\n");
70
71
72
      return false;
    }

73
    rest = 1.0;
74

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

    // check for non domain covering level traverse
78
79
80
81
82
83
84
    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);
    }
85

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

88
89
90
    return true;
  }

91

92
93
94
95
96
97
  bool DualTraverse::traverseNext(ElInfo **elInfo1,
				  ElInfo **elInfo2,
				  ElInfo **elInfoSmall,
				  ElInfo **elInfoLarge)
  {
    // call standard traverse
98
99
100
101
102
103
104
105
106
107
    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
108
109

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

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

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

    // check for non domain covering level traverse
124
125
126
127
128
129
130
    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);
    }
131
132
133

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

134
135
136
    return true;
  }

137

138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
  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
154
    rest -= 1.0 / (1 << ((*elInfoSmall)->getLevel() - (*elInfoLarge)->getLevel()));
155

156
    if (rest < 1e-32) {
157
      // large element finished -> increment both elements
158
159
160
      rest = 1.0; 
      inc1 = true;
      inc2 = true;
161
162
    } else {
      // increment only small element
163
164
      inc1 = (*elInfo1 == *elInfoSmall) ? true : false;
      inc2 = (*elInfo2 == *elInfoSmall) ? true : false;
165
166
167
    }
  }

168

169
170
171
172
173
  void DualTraverse::fillSubElInfo(ElInfo *elInfo1, 
				   ElInfo *elInfo2,
				   ElInfo *elInfoSmall,
				   ElInfo *elInfoLarge)
  {
174
    if (!fillSubElemMat)
175
176
      return;

177
    if (elInfo1 == elInfoSmall)
178
      stack1.fillRefinementPath(*elInfoSmall, *elInfo2);
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
    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;
208
  }
209
}