MeshStructure.cc 6.67 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
#include "MeshStructure.h"
#include "MeshStructure_ED.h"
#include "PartitionElementData.h"
#include "Mesh.h"
#include "Element.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "RefinementManager.h"

namespace AMDiS {

12
  const int MeshStructure::unsignedLongSize = sizeof(unsigned long int) * 8;
13

14

15
16
  void MeshStructure::insertElement(bool isLeaf) 
  {
17
    // overflow? -> next index
18
19
20
21
    if (pos >= unsignedLongSize) {
      code.push_back(currentCode);
      pos = 0;
      currentCode = 0;
22
23
24
    }

    // insert element in binary code
25
    if (!isLeaf) {
26
      unsigned long int one = 1;
27
      currentCode += (one << pos);
28
29
    } 

30
31
    pos++;
    nElements++;
32
33
  }

34

35
36
  void MeshStructure::clear()
  {
37
38
39
40
41
    currentCode = 0;
    code.resize(0);
    pos = 0;
    nElements = 0;
    currentElement = 0;
42
43
  }

44

45
46
47
  void MeshStructure::init(Mesh *mesh) 
  {
    clear();
48

49
50
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
51
    while (elInfo) {
52
53
54
55
56
57
58
      insertElement(elInfo->getElement()->isLeaf());
      elInfo = stack.traverseNext(elInfo);
    }
  
    commit();
  }

59
60
61
62
63
64
65
66
67

  void MeshStructure::init(BoundaryObject &bound)
  {
    FUNCNAME("MeshStructure::init()");

    Element *el = bound.el;

    clear();

68
69
    int s1 = el->getSubObjOfChild(0, bound.subObj, bound.ithObj, bound.elType);
    int s2 = el->getSubObjOfChild(1, bound.subObj, bound.ithObj, bound.elType);
70
71
72
73
74
    
    TEST_EXIT(s1 != -1 || s2 != -1)("This should not happen!\n");
    
    if (!el->isLeaf()) {
      if (s1 == -1)
75
	addAlongSide(el->getSecondChild(), bound.subObj, s2, 
76
77
		     el->getChildType(bound.elType), bound.reverseMode);
      else if (s2 == -1)
78
	addAlongSide(el->getFirstChild(), bound.subObj, s1, 
79
80
		     el->getChildType(bound.elType), bound.reverseMode);
      else
81
	addAlongSide(el, bound.subObj, bound.ithObj, bound.elType, bound.reverseMode);
82
    }
83
84
85
86

    commit();    
  }

87
88
89

  void MeshStructure::addAlongSide(Element *el, GeoIndex subObj, int ithObj, 
				   int elType, bool reverseOrder)
90
  {
91
92
93
94
    FUNCNAME("MeshStructure::addAlongSide()");

    if (debugMode) {
      MSG("addAlondSide(%d, %d, %d, %d)\n",
95
	  el->getIndex(), ithObj, elType, reverseOrder);
96
97
98
      MSG("Element is leaf: %d\n", el->isLeaf());
    }
    
99
100
101
    insertElement(el->isLeaf());
    
    if (!el->isLeaf()) {
102
103
      int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
      int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
104

105
106
      if (!reverseOrder) {
	if (s1 != -1) 
107
	  addAlongSide(el->getFirstChild(), subObj, s1, el->getChildType(elType), false);
108
	if (s2 != -1)
109
	  addAlongSide(el->getSecondChild(), subObj, s2, el->getChildType(elType), false);
110
111
      } else {
	if (s2 != -1)
112
	  addAlongSide(el->getSecondChild(), subObj, s2, el->getChildType(elType), false);
113
	if (s1 != -1) 
114
	  addAlongSide(el->getFirstChild(), subObj, s1, el->getChildType(elType), false);
115
      }
116
117
118
    }
  }

119

120
121
122
123
124
  void MeshStructure::reset() 
  {
    currentIndex = 0;
    pos = 0;
    currentElement = 0;
125
126
127
128
129

    if (code.size() > 0)
      currentCode = code[0];
    else 
      currentCode = 0;
130
131
  }

132

133
134
  bool MeshStructure::nextElement(MeshStructure *insert)
  {
135
136
137
    FUNCNAME("MeshStructure::nextElement()");

    if (insert)
138
139
      insert->insertElement(isLeafElement());

140
141
    pos++;
    currentElement++;
142

143
144
    if (currentElement >= nElements) 
      return false;
145

146
147
148
149
150
151
    if (pos >= unsignedLongSize) {
      currentIndex++;
      TEST_EXIT_DBG(currentIndex < static_cast<int>(code.size()))
	("End of structure reached!\n");
      pos = 0;
      currentCode = code[currentIndex];
152
    } else {
153
      currentCode >>= 1;
154
    }
155

156
157
158
    return true;
  }

159

160
161
  bool MeshStructure::skipBranch(MeshStructure *insert)
  {
162
163
    FUNCNAME("MeshStructure::skipBranch()");

164
    if (isLeafElement()) {
165
166
167
168
      return nextElement(insert);
    } else {
      bool cont = nextElement(insert);
      cont = skipBranch(insert); // left branch
169
      TEST_EXIT_DBG(cont)("Invalid structure!\n");
170
171
172
173
174
      cont = skipBranch(insert); // righ branch
      return cont;
    }
  }

175

176
177
178
179
  void MeshStructure::merge(MeshStructure *structure1,
			    MeshStructure *structure2,
			    MeshStructure *result)
  {
180
181
    FUNCNAME("MeshStructure::merge()");

182
183
184
185
186
    result->clear();
    structure1->reset();
    structure2->reset();

    bool cont = true;
187
    while (cont) {
188
      bool cont1, cont2;
189
      if (structure1->isLeafElement() == structure2->isLeafElement()) {
190
191
192
	cont1 = structure1->nextElement(result);
	cont2 = structure2->nextElement();
      } else {
193
	if (structure1->isLeafElement()) {
194
195
196
197
198
199
200
	  cont1 = structure1->nextElement();
	  cont2 = structure2->skipBranch(result);
	} else {
	  cont1 = structure1->skipBranch(result);
	  cont2 = structure2->nextElement();
	}
      }
201
      TEST_EXIT_DBG(cont1 == cont2)("Structures don't match!\n");
202
203
204
205
206
207
      cont = cont1;
    }

    result->commit();
  }

208

209
210
  void MeshStructure::fitMeshToStructure(Mesh *mesh,
					 RefinementManager *manager,
211
212
					 bool checkPartition,
					 bool debugMode) 
213
214
215
216
217
218
219
  {
    FUNCNAME("MeshStructure::fitMeshToStructure()");

    bool cont = true;

    // decorate leaf data
    reset();
220
221
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
222
    while (elInfo) {
223
      TEST_EXIT(cont)("unexpected structure code end!\n");
224
225
226

      Element *element = elInfo->getElement();

227
      if (isLeafElement()) {
228
	TEST_EXIT_DBG(element->isLeaf())("mesh finer than code\n");
229
      }
230

231
      if (element->isLeaf() && !isLeafElement()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
232
	MeshStructure *structure = new MeshStructure();
233
234
235
236
	cont = skipBranch(structure);
	structure->commit();

	bool decorate = true;
237
	if (checkPartition) {
238
239
	  PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
	    (element->getElementData(PARTITION_ED));
240
	  TEST_EXIT_DBG(partitionData)("no partition element data\n");
241
	  PartitionStatus status = partitionData->getPartitionStatus();
242
	  if (debugMode == false && (status == OUT || status == UNDEFINED))
243
244
245
	    decorate = false;
	}

246
	if (decorate) {
Thomas Witkowski's avatar
Thomas Witkowski committed
247
	  MeshStructure_ED *elData = new MeshStructure_ED(element->getElementData());
248
249
250
	  elData->setStructure(structure);
	  element->setElementData(elData);
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
251
	  delete structure;
252
253
254
255
256
257
258
259
260
	}
      } else {
	cont = nextElement();
      }

      elInfo = stack.traverseNext(elInfo);
    }

    // refine mesh
261
    bool finished = true;
262
263
264
265

    do {
      finished = true;
      elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
266
      while (elInfo) {
267
	Element *element = elInfo->getElement();
268
	if (element->getElementData(MESH_STRUCTURE) != NULL) {
269
270
271
272
273
274
275
276
	  element->setMark(1);
	  finished = false;
	} else {
	  element->setMark(0);
	}
	elInfo = stack.traverseNext(elInfo);
      }
      manager->refineMesh(mesh);
277
    } while (!finished);
278
279
  }

280

281
282
283
284
  bool MeshStructure::compare(MeshStructure &other)
  {
    return (other.getCode() == code);
  }
285
}