MeshStructure.cc 10.1 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
#include "DOFVector.h"
10
11
12

namespace AMDiS {

13
  const int MeshStructure::structureSize = 64;
14

15

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

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

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

35

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

45

46
  void MeshStructure::init(Mesh *mesh, int macroElIndex) 
47
48
  {
    clear();
49

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

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

63
64
65
66
67
68
69

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

    Element *el = bound.el;

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

72
73
    clear();

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

    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);
    }
85
86
87
    
    if (!el->isLeaf()) {
      if (s1 == -1)
88
	addAlongSide(el->getSecondChild(), bound.subObj, s2, 
89
90
		     el->getChildType(bound.elType), bound.reverseMode);
      else if (s2 == -1)
91
	addAlongSide(el->getFirstChild(), bound.subObj, s1, 
92
93
		     el->getChildType(bound.elType), bound.reverseMode);
      else
94
	addAlongSide(el, bound.subObj, bound.ithObj, bound.elType, bound.reverseMode);
95
    }
96
97
98
99

    commit();    
  }

100
101
102

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

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

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

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

140

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

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

153

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

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

161
162
    pos++;
    currentElement++;
163

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

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

177
178
179
    return true;
  }

180

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

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

196

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

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

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

    result->commit();
  }

229

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

    bool cont = true;

    // decorate leaf data
    reset();
242
    TraverseStack stack;
243
244
245
246
247
    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);
248

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

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

254
      if (isLeafElement()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
255
256
257
258
259
260
	if (ignoreFinerMesh && !element->isLeaf()) {
	  int level = elInfo->getLevel();
	  while (elInfo && level >= elInfo->getLevel())
	    elInfo = stack.traverseNext(elInfo);
	} else {
	  TEST_EXIT(element->isLeaf())
261
262
	    ("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
263
264
	}
      } 
265

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

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

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

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

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

      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
      }
317
    } while (!finished);
318
319
  }

320

321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
  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();
      }
    }
343
344
345
346
347
348
349
350
351
352
   
    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
    }
353
354
355
  }


356
357
358
359
  bool MeshStructure::compare(MeshStructure &other)
  {
    return (other.getCode() == code);
  }
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419


  void MeshStructure::getMeshStructureValues(Mesh *mesh,
					     int macroElIndex,
					     const DOFVector<double>* vec,
					     std::vector<double>& values)
  {
    FUNCNAME("MeshStructure::getMeshStructureValues()");

    values.clear();

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1, 
						 Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
      if (elInfo->getLevel() == 0) {
	for (int i = 0; i < mesh->getGeo(VERTEX); i++)
	  values.push_back((*vec)[elInfo->getElement()->getDof(i, 0)]);
      } else {
	if (!elInfo->getElement()->isLeaf())
	  values.push_back((*vec)[elInfo->getElement()->getChild(0)->getDof(mesh->getDim(), 0)]);
      }

      elInfo = stack.traverseNext(elInfo);
    }
  }


  void MeshStructure::setMeshStructureValues(Mesh *mesh,
					     int macroElIndex,
					     DOFVector<double>* vec,
					     const std::vector<double>& values)
  {
    FUNCNAME("MeshStructure::setMeshStructureValues()");

    TEST_EXIT_DBG(values.size() >= mesh->getGeo(VERTEX))("Should not happen!\n");

    unsigned int counter = 0;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1, 
						 Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
      if (elInfo->getLevel() == 0) {
	for (int i = 0; i < mesh->getGeo(VERTEX); i++)
	  (*vec)[elInfo->getElement()->getDof(i, 0)] = values[counter++];
      } else {
	if (!elInfo->getElement()->isLeaf()) {
	  TEST_EXIT_DBG(counter < values.size())("Should not happen!\n");

	  (*vec)[elInfo->getElement()->getChild(0)->getDof(mesh->getDim(), 0)] =
	    values[counter++];
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }
      
  }

420
}