MeshStructure.cc 11.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
14
15
16
17
18
19
#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
20
#include "Debug.h"
21
#include "DOFVector.h"
22
23
24

namespace AMDiS {

25
  const int MeshStructure::structureSize = 64;
26

27

28
29
  void MeshStructure::insertElement(bool isLeaf) 
  {
30
    // overflow? -> next index
31
    if (pos >= structureSize) {
32
33
34
      code.push_back(currentCode);
      pos = 0;
      currentCode = 0;
35
36
37
    }

    // insert element in binary code
38
    if (!isLeaf) {
39
      uint64_t one = 1;
40
      currentCode += (one << pos);
41
42
    } 

43
44
    pos++;
    nElements++;
45
46
  }

47

48
49
  void MeshStructure::clear()
  {
50
51
52
53
54
    currentCode = 0;
    code.resize(0);
    pos = 0;
    nElements = 0;
    currentElement = 0;
55
56
  }

57

58
  void MeshStructure::init(Mesh *mesh, int macroElIndex) 
59
60
  {
    clear();
61

62
63
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
64
    while (elInfo) {
65
66
67
68
      if (macroElIndex == -1 ||
	  elInfo->getMacroElement()->getIndex() == macroElIndex)
	insertElement(elInfo->getElement()->isLeaf());

69
70
71
72
73
74
      elInfo = stack.traverseNext(elInfo);
    }
  
    commit();
  }

75
76
77
78
79
80
81

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

    Element *el = bound.el;

82
83
    TEST_EXIT_DBG(el)("No element!\n");

84
85
    clear();

86
87
    int s1 = el->getSubObjOfChild(0, bound.subObj, bound.ithObj, bound.elType);
    int s2 = el->getSubObjOfChild(1, bound.subObj, bound.ithObj, bound.elType);
88
89
    
    TEST_EXIT(s1 != -1 || s2 != -1)("This should not happen!\n");
90
91
92
93
94
95
96

    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);
    }
97
98
99
    
    if (!el->isLeaf()) {
      if (s1 == -1)
100
	addAlongSide(el->getSecondChild(), bound.subObj, s2, 
101
102
		     el->getChildType(bound.elType), bound.reverseMode);
      else if (s2 == -1)
103
	addAlongSide(el->getFirstChild(), bound.subObj, s1, 
104
105
		     el->getChildType(bound.elType), bound.reverseMode);
      else
106
	addAlongSide(el, bound.subObj, bound.ithObj, bound.elType, bound.reverseMode);
107
    }
108
109
110
111

    commit();    
  }

112
113
114

  void MeshStructure::addAlongSide(Element *el, GeoIndex subObj, int ithObj, 
				   int elType, bool reverseOrder)
115
  {
116
117
118
119
    FUNCNAME("MeshStructure::addAlongSide()");

    if (debugMode) {
      MSG("addAlondSide(%d, %d, %d, %d)\n",
120
	  el->getIndex(), ithObj, elType, reverseOrder);
121
122
123
      MSG("Element is leaf: %d\n", el->isLeaf());
    }
    
124
125
126
    insertElement(el->isLeaf());
    
    if (!el->isLeaf()) {
127
128
      int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
      int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
129

130
131
132
133
134
135
136
137
      if (debugMode) {
	MSG("Child index %d  %d\n", 
	    el->getFirstChild()->getIndex(),
	    el->getSecondChild()->getIndex());
	MSG("s1 = %d    s2 = %d\n", s1, s2);
	MSG("   \n");
      }

138
139
      if (!reverseOrder) {
	if (s1 != -1) 
140
	  addAlongSide(el->getFirstChild(), subObj, s1, el->getChildType(elType), reverseOrder);
141
	if (s2 != -1)
142
	  addAlongSide(el->getSecondChild(), subObj, s2, el->getChildType(elType), reverseOrder);
143
144
      } else {
	if (s2 != -1)
145
	  addAlongSide(el->getSecondChild(), subObj, s2, el->getChildType(elType), reverseOrder);
146
	if (s1 != -1) 
147
	  addAlongSide(el->getFirstChild(), subObj, s1, el->getChildType(elType), reverseOrder);
148
      }
149
150
151
    }
  }

152

153
154
155
156
157
  void MeshStructure::reset() 
  {
    currentIndex = 0;
    pos = 0;
    currentElement = 0;
158
159
160
161
162

    if (code.size() > 0)
      currentCode = code[0];
    else 
      currentCode = 0;
163
164
  }

165

166
167
  bool MeshStructure::nextElement(MeshStructure *insert)
  {
168
169
170
    FUNCNAME("MeshStructure::nextElement()");

    if (insert)
171
172
      insert->insertElement(isLeafElement());

173
174
    pos++;
    currentElement++;
175

176
177
    if (currentElement >= nElements) 
      return false;
178

179
    if (pos >= structureSize) {
180
181
182
183
184
      currentIndex++;
      TEST_EXIT_DBG(currentIndex < static_cast<int>(code.size()))
	("End of structure reached!\n");
      pos = 0;
      currentCode = code[currentIndex];
185
    } else {
186
      currentCode >>= 1;
187
    }
188

189
190
191
    return true;
  }

192

193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
  int MeshStructure::lookAhead(unsigned int n)
  {
    FUNCNAME("MeshStructure::lookAhead()");

    int returnValue = 0;

    int tmp_pos = pos;
    int tmp_currentElement = currentElement;
    int tmp_currentIndex = currentIndex;
    uint64_t tmp_currentCode = currentCode;

    for (unsigned int i = 0; i < n; i++) {
      if (nextElement() == false) {
	returnValue = -1;      
	break;
      }
    }

    if (returnValue != -1)
      returnValue = static_cast<int>(!isLeafElement());

    pos = tmp_pos;
    currentElement = tmp_currentElement;
    currentIndex = tmp_currentIndex;
    currentCode = tmp_currentCode;

    return returnValue;
  }


223
224
  bool MeshStructure::skipBranch(MeshStructure *insert)
  {
225
226
    FUNCNAME("MeshStructure::skipBranch()");

227
    if (isLeafElement()) {
228
229
230
231
      return nextElement(insert);
    } else {
      bool cont = nextElement(insert);
      cont = skipBranch(insert); // left branch
232
      TEST_EXIT_DBG(cont)("Invalid structure!\n");
233
234
235
236
237
      cont = skipBranch(insert); // righ branch
      return cont;
    }
  }

238

239
240
241
242
  void MeshStructure::merge(MeshStructure *structure1,
			    MeshStructure *structure2,
			    MeshStructure *result)
  {
243
244
    FUNCNAME("MeshStructure::merge()");

245
246
247
248
249
    result->clear();
    structure1->reset();
    structure2->reset();

    bool cont = true;
250
    while (cont) {
251
      bool cont1, cont2;
252
      if (structure1->isLeafElement() == structure2->isLeafElement()) {
253
254
255
	cont1 = structure1->nextElement(result);
	cont2 = structure2->nextElement();
      } else {
256
	if (structure1->isLeafElement()) {
257
258
259
260
261
262
263
	  cont1 = structure1->nextElement();
	  cont2 = structure2->skipBranch(result);
	} else {
	  cont1 = structure1->skipBranch(result);
	  cont2 = structure2->nextElement();
	}
      }
264
      TEST_EXIT_DBG(cont1 == cont2)("Structures don't match!\n");
265
266
267
268
269
270
      cont = cont1;
    }

    result->commit();
  }

271

272
273
  void MeshStructure::fitMeshToStructure(Mesh *mesh,
					 RefinementManager *manager,
274
					 bool debugMode,
Thomas Witkowski's avatar
Thomas Witkowski committed
275
276
					 int macroElIndex,
					 bool ignoreFinerMesh) 
277
278
279
280
281
282
283
  {
    FUNCNAME("MeshStructure::fitMeshToStructure()");

    bool cont = true;

    // decorate leaf data
    reset();
284
    TraverseStack stack;
285
286
287
288
289
    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);
290

291

292
    while (elInfo) {
293
294
      Element *element = elInfo->getElement();

295
296
      TEST_EXIT(cont)("unexpected structure code end!\n");

297
      if (isLeafElement()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
298
299
300
301
302
303
	if (ignoreFinerMesh && !element->isLeaf()) {
	  int level = elInfo->getLevel();
	  while (elInfo && level >= elInfo->getLevel())
	    elInfo = stack.traverseNext(elInfo);
	} else {
	  TEST_EXIT(element->isLeaf())
304
305
	    ("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
306
307
	}
      } 
308

309
310
      TEST_EXIT_DBG(element)("Should not happen!\n");

311
      if (element->isLeaf() && !isLeafElement()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
312
	MeshStructure *structure = new MeshStructure();
313
314
315
	cont = skipBranch(structure);
	structure->commit();

316
317
318
	MeshStructure_ED *elData = new MeshStructure_ED(element->getElementData());
	elData->setStructure(structure);
	element->setElementData(elData);
319
320
321
322
      } else {
	cont = nextElement();
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
323
324
      if (elInfo)
	elInfo = stack.traverseNext(elInfo);
325
326
327
    }

    // refine mesh
328
    bool finished = true;
329

330
    do {
331
      finished = true;
332
333
334
335
      if (macroElIndex == -1)
	elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
      else
	elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1, Mesh::CALL_LEAF_EL);      
336
      while (elInfo) {
337
	Element *element = elInfo->getElement();
338
	if (element->getElementData(MESH_STRUCTURE) != NULL) {
339
340
341
342
343
344
345
	  element->setMark(1);
	  finished = false;
	} else {
	  element->setMark(0);
	}
	elInfo = stack.traverseNext(elInfo);
      }
346
347
348
349
350

      if (!finished) {
#if (DEBUG != 0)
	int oldMeshIndex = mesh->getChangeIndex();
#endif
351

352
353
354
355
	if (macroElIndex == -1)
	  manager->refineMesh(mesh);
	else
	  manager->refineMacroElement(mesh, macroElIndex);
356

357
358
359
360
361
#if (DEBUG != 0)
	TEST_EXIT(oldMeshIndex != mesh->getChangeIndex())
	  ("Mesh has not been changed by refinement procedure!\n");
#endif
      }
362
    } while (!finished);
363
364
  }

365

Thomas Witkowski's avatar
Thomas Witkowski committed
366
  string MeshStructure::toStr(bool resetCode)
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
  {
    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();
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
386
387
388
389
390
391
392
393
394
395

    return oss.str();
  }


  void MeshStructure::print(bool resetCode)
  {
    FUNCNAME("MeshStructure::print()");
    
    string str = toStr(resetCode);
396
   
Thomas Witkowski's avatar
Thomas Witkowski committed
397
398
    if (str.length() < 255) {
      MSG("Mesh structure code: %s\n", str.c_str());
399
400
    } else {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
401
      std::cout << "[" << MPI::COMM_WORLD.Get_rank() << "]                Mesh structure code: " << str << "\n";
402
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
403
      std::cout << "                Mesh structure code: " << str << "\n";
404
405
#endif
    }
406
407
408
  }


409
410
411
412
  bool MeshStructure::compare(MeshStructure &other)
  {
    return (other.getCode() == code);
  }
413
414
415
416
417
418
419
420
421


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

422
423
424
    TEST_EXIT_DBG(mesh)("No mesh defined!\n");
    TEST_EXIT_DBG(vec)("No DOFVector defined!\n");

425
426
427
428
429
430
    values.clear();

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1, 
						 Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
431
      if (elInfo->getLevel() == 0)
432
433
	for (int i = 0; i < mesh->getGeo(VERTEX); i++)
	  values.push_back((*vec)[elInfo->getElement()->getDof(i, 0)]);
434
435
436

      if (!elInfo->getElement()->isLeaf())
	values.push_back((*vec)[elInfo->getElement()->getChild(0)->getDof(mesh->getDim(), 0)]);      
437
438
439
440
441
442
443
444
445
446
447
448
449

      elInfo = stack.traverseNext(elInfo);
    }
  }


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

450
451
452
453
    TEST_EXIT_DBG(mesh)("No mesh defined!\n");
    TEST_EXIT_DBG(vec)("No DOFVector defined!\n");
    TEST_EXIT_DBG(static_cast<int>(values.size()) >= mesh->getGeo(VERTEX))
      ("Should not happen!\n");
454
455
456
457
458
459
460

    unsigned int counter = 0;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1, 
						 Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
461
      if (elInfo->getLevel() == 0)
462
	for (int i = 0; i < mesh->getGeo(VERTEX); i++)
463
	  (*vec)[elInfo->getElement()->getDof(i, 0)] = values[counter++];      
464

465
466
467
468
469
      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++];      
470
471
472
473
474
475
476
      }

      elInfo = stack.traverseNext(elInfo);
    }
      
  }

477
}