MeshStructure.cc 11.4 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
141
	  addAlongSide(el->getFirstChild(), subObj, s1, 
		       el->getChildType(elType), reverseOrder);
142
	if (s2 != -1)
143
144
	  addAlongSide(el->getSecondChild(), subObj, s2, 
		       el->getChildType(elType), reverseOrder);
145
146
      } else {
	if (s2 != -1)
147
148
	  addAlongSide(el->getSecondChild(), subObj, s2, 
		       el->getChildType(elType), reverseOrder);
149
	if (s1 != -1) 
150
151
	  addAlongSide(el->getFirstChild(), subObj, s1, 
		       el->getChildType(elType), reverseOrder);
152
      }
153
154
155
    }
  }

156

157
158
159
160
161
  void MeshStructure::reset() 
  {
    currentIndex = 0;
    pos = 0;
    currentElement = 0;
162
163
164
165
166

    if (code.size() > 0)
      currentCode = code[0];
    else 
      currentCode = 0;
167
168
  }

169

170
171
  bool MeshStructure::nextElement(MeshStructure *insert)
  {
172
173
174
    FUNCNAME("MeshStructure::nextElement()");

    if (insert)
175
176
      insert->insertElement(isLeafElement());

177
178
    pos++;
    currentElement++;
179

180
181
    if (currentElement >= nElements) 
      return false;
182

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

193
194
195
    return true;
  }

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
223
224
225
226
  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;
  }


227
228
  bool MeshStructure::skipBranch(MeshStructure *insert)
  {
229
230
    FUNCNAME("MeshStructure::skipBranch()");

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

242

243
244
245
246
  void MeshStructure::merge(MeshStructure *structure1,
			    MeshStructure *structure2,
			    MeshStructure *result)
  {
247
248
    FUNCNAME("MeshStructure::merge()");

249
250
251
252
253
    result->clear();
    structure1->reset();
    structure2->reset();

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

    result->commit();
  }

275

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

    bool cont = true;

    // decorate leaf data
    reset();
288
    TraverseStack stack;
289
290
291
292
293
    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);
294

295

296
    while (elInfo) {
297
298
      Element *element = elInfo->getElement();

299
300
      TEST_EXIT(cont)("unexpected structure code end!\n");

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

313
314
      TEST_EXIT_DBG(element)("Should not happen!\n");

315
      if (element->isLeaf() && !isLeafElement()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
316
	MeshStructure *structure = new MeshStructure();
317
318
319
	cont = skipBranch(structure);
	structure->commit();

320
321
322
	MeshStructure_ED *elData = new MeshStructure_ED(element->getElementData());
	elData->setStructure(structure);
	element->setElementData(elData);
323
324
325
326
      } else {
	cont = nextElement();
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
327
328
      if (elInfo)
	elInfo = stack.traverseNext(elInfo);
329
330
331
    }

    // refine mesh
332
    bool finished = true;
333

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

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

356
357
358
359
	if (macroElIndex == -1)
	  manager->refineMesh(mesh);
	else
	  manager->refineMacroElement(mesh, macroElIndex);
360

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

369

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

    return oss.str();
  }


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


413
414
415
416
  bool MeshStructure::compare(MeshStructure &other)
  {
    return (other.getCode() == code);
  }
417
418
419
420
421
422
423
424
425


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

426
427
428
    TEST_EXIT_DBG(mesh)("No mesh defined!\n");
    TEST_EXIT_DBG(vec)("No DOFVector defined!\n");

429
430
431
432
433
434
    values.clear();

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

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

      elInfo = stack.traverseNext(elInfo);
    }
  }


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

454
455
456
457
    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");
458
459
460
461
462
463
464

    unsigned int counter = 0;

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

469
470
471
472
473
      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++];      
474
475
476
477
478
479
480
      }

      elInfo = stack.traverseNext(elInfo);
    }
      
  }

481
}