MeshStructure.cc 8.41 KB
Newer Older
1
2
3
4
5
6
7
#include "MeshStructure.h"
#include "MeshStructure_ED.h"
#include "Mesh.h"
#include "Element.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "RefinementManager.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
8
#include "Debug.h"
9
10
11

namespace AMDiS {

12
  const int MeshStructure::structureSize = 64;
13

14

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

    // insert element in binary code
25
    if (!isLeaf) {
26
      uint64_t 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
  void MeshStructure::init(Mesh *mesh, int macroElIndex) 
46
47
  {
    clear();
48

49
50
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
51
    while (elInfo) {
52
53
54
55
      if (macroElIndex == -1 ||
	  elInfo->getMacroElement()->getIndex() == macroElIndex)
	insertElement(elInfo->getElement()->isLeaf());

56
57
58
59
60
61
      elInfo = stack.traverseNext(elInfo);
    }
  
    commit();
  }

62
63
64
65
66
67
68

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

    Element *el = bound.el;

69
70
    TEST_EXIT_DBG(el)("No element!\n");

71
72
    clear();

73
74
    int s1 = el->getSubObjOfChild(0, bound.subObj, bound.ithObj, bound.elType);
    int s2 = el->getSubObjOfChild(1, bound.subObj, bound.ithObj, bound.elType);
75
76
    
    TEST_EXIT(s1 != -1 || s2 != -1)("This should not happen!\n");
77
78
79
80
81
82
83

    if (debugMode) {
      MSG("addAlondSide(%d, %d, %d, %d)\n",
	  bound.elIndex, bound.ithObj, bound.elType, bound.reverseMode);
      MSG("Element is leaf: %d\n", el->isLeaf());
      MSG("s1 = %d    s2 = %d\n", s1, s2);
    }
84
85
86
    
    if (!el->isLeaf()) {
      if (s1 == -1)
87
	addAlongSide(el->getSecondChild(), bound.subObj, s2, 
88
89
		     el->getChildType(bound.elType), bound.reverseMode);
      else if (s2 == -1)
90
	addAlongSide(el->getFirstChild(), bound.subObj, s1, 
91
92
		     el->getChildType(bound.elType), bound.reverseMode);
      else
93
	addAlongSide(el, bound.subObj, bound.ithObj, bound.elType, bound.reverseMode);
94
    }
95
96
97
98

    commit();    
  }

99
100
101

  void MeshStructure::addAlongSide(Element *el, GeoIndex subObj, int ithObj, 
				   int elType, bool reverseOrder)
102
  {
103
104
105
106
    FUNCNAME("MeshStructure::addAlongSide()");

    if (debugMode) {
      MSG("addAlondSide(%d, %d, %d, %d)\n",
107
	  el->getIndex(), ithObj, elType, reverseOrder);
108
109
110
      MSG("Element is leaf: %d\n", el->isLeaf());
    }
    
111
112
113
    insertElement(el->isLeaf());
    
    if (!el->isLeaf()) {
114
115
      int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
      int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
116

117
118
119
120
121
122
123
124
      if (debugMode) {
	MSG("Child index %d  %d\n", 
	    el->getFirstChild()->getIndex(),
	    el->getSecondChild()->getIndex());
	MSG("s1 = %d    s2 = %d\n", s1, s2);
	MSG("   \n");
      }

125
126
      if (!reverseOrder) {
	if (s1 != -1) 
127
	  addAlongSide(el->getFirstChild(), subObj, s1, el->getChildType(elType), reverseOrder);
128
	if (s2 != -1)
129
	  addAlongSide(el->getSecondChild(), subObj, s2, el->getChildType(elType), reverseOrder);
130
131
      } else {
	if (s2 != -1)
132
	  addAlongSide(el->getSecondChild(), subObj, s2, el->getChildType(elType), reverseOrder);
133
	if (s1 != -1) 
134
	  addAlongSide(el->getFirstChild(), subObj, s1, el->getChildType(elType), reverseOrder);
135
      }
136
137
138
    }
  }

139

140
141
142
143
144
  void MeshStructure::reset() 
  {
    currentIndex = 0;
    pos = 0;
    currentElement = 0;
145
146
147
148
149

    if (code.size() > 0)
      currentCode = code[0];
    else 
      currentCode = 0;
150
151
  }

152

153
154
  bool MeshStructure::nextElement(MeshStructure *insert)
  {
155
156
157
    FUNCNAME("MeshStructure::nextElement()");

    if (insert)
158
159
      insert->insertElement(isLeafElement());

160
161
    pos++;
    currentElement++;
162

163
164
    if (currentElement >= nElements) 
      return false;
165

166
    if (pos >= structureSize) {
167
168
169
170
171
      currentIndex++;
      TEST_EXIT_DBG(currentIndex < static_cast<int>(code.size()))
	("End of structure reached!\n");
      pos = 0;
      currentCode = code[currentIndex];
172
    } else {
173
      currentCode >>= 1;
174
    }
175

176
177
178
    return true;
  }

179

180
181
  bool MeshStructure::skipBranch(MeshStructure *insert)
  {
182
183
    FUNCNAME("MeshStructure::skipBranch()");

184
    if (isLeafElement()) {
185
186
187
188
      return nextElement(insert);
    } else {
      bool cont = nextElement(insert);
      cont = skipBranch(insert); // left branch
189
      TEST_EXIT_DBG(cont)("Invalid structure!\n");
190
191
192
193
194
      cont = skipBranch(insert); // righ branch
      return cont;
    }
  }

195

196
197
198
199
  void MeshStructure::merge(MeshStructure *structure1,
			    MeshStructure *structure2,
			    MeshStructure *result)
  {
200
201
    FUNCNAME("MeshStructure::merge()");

202
203
204
205
206
    result->clear();
    structure1->reset();
    structure2->reset();

    bool cont = true;
207
    while (cont) {
208
      bool cont1, cont2;
209
      if (structure1->isLeafElement() == structure2->isLeafElement()) {
210
211
212
	cont1 = structure1->nextElement(result);
	cont2 = structure2->nextElement();
      } else {
213
	if (structure1->isLeafElement()) {
214
215
216
217
218
219
220
	  cont1 = structure1->nextElement();
	  cont2 = structure2->skipBranch(result);
	} else {
	  cont1 = structure1->skipBranch(result);
	  cont2 = structure2->nextElement();
	}
      }
221
      TEST_EXIT_DBG(cont1 == cont2)("Structures don't match!\n");
222
223
224
225
226
227
      cont = cont1;
    }

    result->commit();
  }

228

229
230
  void MeshStructure::fitMeshToStructure(Mesh *mesh,
					 RefinementManager *manager,
231
					 bool debugMode,
Thomas Witkowski's avatar
Thomas Witkowski committed
232
233
					 int macroElIndex,
					 bool ignoreFinerMesh) 
234
235
236
237
238
239
240
  {
    FUNCNAME("MeshStructure::fitMeshToStructure()");

    bool cont = true;

    // decorate leaf data
    reset();
241
    TraverseStack stack;
242
243
244
245
246
    ElInfo *elInfo = NULL;
    if (macroElIndex == -1)
      elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
    else
      elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1, Mesh::CALL_EVERY_EL_PREORDER);
247

248
    while (elInfo) {
249
250
      Element *element = elInfo->getElement();

251
252
      TEST_EXIT(cont)("unexpected structure code end!\n");

253
      if (isLeafElement()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
254
255
256
257
258
259
	if (ignoreFinerMesh && !element->isLeaf()) {
	  int level = elInfo->getLevel();
	  while (elInfo && level >= elInfo->getLevel())
	    elInfo = stack.traverseNext(elInfo);
	} else {
	  TEST_EXIT(element->isLeaf())
260
261
	    ("Mesh is finer than strucutre code! (Element index: %d Macro element index: %d)\n", 
	     element->getIndex(), elInfo->getMacroElement()->getIndex());
Thomas Witkowski's avatar
Thomas Witkowski committed
262
263
	}
      } 
264

265
      if (element->isLeaf() && !isLeafElement()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
266
	MeshStructure *structure = new MeshStructure();
267
268
269
	cont = skipBranch(structure);
	structure->commit();

270
271
272
	MeshStructure_ED *elData = new MeshStructure_ED(element->getElementData());
	elData->setStructure(structure);
	element->setElementData(elData);
273
274
275
276
      } else {
	cont = nextElement();
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
277
278
      if (elInfo)
	elInfo = stack.traverseNext(elInfo);
279
280
281
    }

    // refine mesh
282
    bool finished = true;
283
284
285

    do {
      finished = true;
286
287
288
289
      if (macroElIndex == -1)
	elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
      else
	elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1, Mesh::CALL_LEAF_EL);      
290
      while (elInfo) {
291
	Element *element = elInfo->getElement();
292
	if (element->getElementData(MESH_STRUCTURE) != NULL) {
293
294
295
296
297
298
299
	  element->setMark(1);
	  finished = false;
	} else {
	  element->setMark(0);
	}
	elInfo = stack.traverseNext(elInfo);
      }
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315

      if (!finished) {
#if (DEBUG != 0)
	int oldMeshIndex = mesh->getChangeIndex();
#endif
	
	if (macroElIndex == -1)
	  manager->refineMesh(mesh);
	else
	  manager->refineMacroElement(mesh, macroElIndex);
	
#if (DEBUG != 0)
	TEST_EXIT(oldMeshIndex != mesh->getChangeIndex())
	  ("Mesh has not been changed by refinement procedure!\n");
#endif
      }
316
    } while (!finished);
317
318
  }

319

320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
  void MeshStructure::print(bool resetCode)
  {
    FUNCNAME("MeshStructure::print()");
    
    std::stringstream oss;
    
    if (empty()) {
      oss << "-" << std::endl;
    }	else {	
      if (resetCode)
	reset();

      bool cont = true;
      while (cont) {
	if (isLeafElement())
	  oss << "0";
	else
	  oss << "1";
	
	cont = nextElement();
      }
    }
342
343
344
345
346
347
348
349
350
351
   
    if (oss.str().length() < 255) {
      MSG("Mesh structure code: %s\n", oss.str().c_str());
    } else {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      std::cout << "[" << MPI::COMM_WORLD.Get_rank() << "]                Mesh structure code: " << oss.str() << "\n";
#else
      std::cout << "                Mesh structure code: " << oss.str() << "\n";
#endif
    }
352
353
354
  }


355
356
357
358
  bool MeshStructure::compare(MeshStructure &other)
  {
    return (other.getCode() == code);
  }
359
}