DualTraverse.cc 5.32 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
59
60
    while (*elInfo1 != NULL && skipEl1(*elInfo1)) {
	    *elInfo1 = stack1.traverseNext(*elInfo1);
    }
61
    *elInfo2 = stack2.traverseFirst(mesh2, level2, flag2);
Naumann, Andreas's avatar
Naumann, Andreas committed
62
63
64
    while (*elInfo2 != NULL && skipEl2(*elInfo2)) {
	    *elInfo2 = stack2.traverseNext(*elInfo2);
    }
65
 
66
    // finished ?
67
    if (*elInfo1 == NULL || *elInfo2 == NULL) {
68
      TEST_EXIT(*elInfo1 == *elInfo2)("invalid dual traverse\n");
69
70
71
      return false;
    }

72
    rest = 1.0;
73

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

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

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

87
88
89
90
91
92
93
94
95
    return true;
  }

  bool DualTraverse::traverseNext(ElInfo **elInfo1,
				  ElInfo **elInfo2,
				  ElInfo **elInfoSmall,
				  ElInfo **elInfoLarge)
  {
    // call standard traverse
Naumann, Andreas's avatar
Naumann, Andreas committed
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
	  if (inc1) {
		  do {
			  *elInfo1 = stack1.traverseNext(*elInfo1);
		  } while(*elInfo1 != NULL && skipEl1(*elInfo1));
	  }
	  if (inc2) {
		  do {
			  *elInfo2 = stack2.traverseNext(*elInfo2);
		  } while (*elInfo2 != NULL && skipEl2(*elInfo2));
	  }

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

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

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

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

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

132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
    return true;
  }

  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
151
    rest -= 1.0 / (1 << ((*elInfoSmall)->getLevel() - (*elInfoLarge)->getLevel()));
152

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

165
166
167
168
169
  void DualTraverse::fillSubElInfo(ElInfo *elInfo1, 
				   ElInfo *elInfo2,
				   ElInfo *elInfoSmall,
				   ElInfo *elInfoLarge)
  {
170
    if (!fillSubElemMat)
171
172
      return;

173
174
175
176
177
//     mtl::dense2D<double>& subCoordsMat = 
//       const_cast<mtl::dense2D<double>&>(elInfoSmall->getSubElemCoordsMat());
//     mtl::dense2D<double>& subCoordsMat_so = 
//       const_cast<mtl::dense2D<double>&>(elInfoSmall->getSubElemCoordsMat_so());

178
    if (elInfo1 == elInfoSmall) {
179
180
181
//       stack1.getCoordsInElem(elInfo2, basisFcts, subCoordsMat);
//       stack1.getCoordsInElem_so(elInfo2, basisFcts, subCoordsMat_so);
      stack1.fillRefinementPath(*elInfoSmall, *elInfo2);
182
    } else {
183
184
185
//       stack2.getCoordsInElem(elInfo1, basisFcts, subCoordsMat);
//       stack2.getCoordsInElem_so(elInfo1, basisFcts, subCoordsMat_so);
      stack2.fillRefinementPath(*elInfoSmall, *elInfo1);
186
187
    }
  }
188
}