MeshStructure.cc 14 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
20
21


22
23
24
25
26
#include "Debug.h"
#include "DOFVector.h"
#include "Element.h"
#include "ElementDofIterator.h"
#include "ElInfo.h"
27
28
29
30
#include "MeshStructure.h"
#include "MeshStructure_ED.h"
#include "Mesh.h"
#include "RefinementManager.h"
31
32
#include "Traverse.h"

33

34
35
using namespace std;

36
37
namespace AMDiS {

38
  const int MeshStructure::structureSize = 64;
39

40

41
42
  void MeshStructure::insertElement(bool isLeaf) 
  {
43
    // overflow? -> next index
44
    if (pos >= structureSize) {
45
46
47
      code.push_back(currentCode);
      pos = 0;
      currentCode = 0;
48
49
50
    }

    // insert element in binary code
51
    if (!isLeaf) {
52
      uint64_t one = 1;
53
      currentCode += (one << pos);
54
55
    } 

56
57
    pos++;
    nElements++;
58
59
  }

60

61
62
  void MeshStructure::clear()
  {
63
64
65
66
67
    currentCode = 0;
    code.resize(0);
    pos = 0;
    nElements = 0;
    currentElement = 0;
68
69
  }

70

71
  void MeshStructure::init(Mesh *mesh, int macroElIndex) 
72
73
  {
    clear();
74

75
    TraverseStack stack;
76

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
77
78
79
80
81
82
83
84
85
    ElInfo *elInfo;
    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);

    while (elInfo) {
      insertElement(elInfo->getElement()->isLeaf());
86
87
88
89
90
91
      elInfo = stack.traverseNext(elInfo);
    }
  
    commit();
  }

92
93
94
95
96
97
98

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

    Element *el = bound.el;

99
100
    TEST_EXIT_DBG(el)("No element!\n");

101
102
    clear();

103
104
    int s1 = el->getSubObjOfChild(0, bound.subObj, bound.ithObj, bound.elType);
    int s2 = el->getSubObjOfChild(1, bound.subObj, bound.ithObj, bound.elType);
105
106
    
    TEST_EXIT(s1 != -1 || s2 != -1)("This should not happen!\n");
107
108
109
110
111
112
113

    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);
    }
114
115
116
    
    if (!el->isLeaf()) {
      if (s1 == -1)
117
	addAlongSide(el->getSecondChild(), bound.subObj, s2, 
118
119
		     el->getChildType(bound.elType), bound.reverseMode);
      else if (s2 == -1)
120
	addAlongSide(el->getFirstChild(), bound.subObj, s1, 
121
122
		     el->getChildType(bound.elType), bound.reverseMode);
      else
123
	addAlongSide(el, bound.subObj, bound.ithObj, bound.elType, bound.reverseMode);
124
    }
125
126
127
128

    commit();    
  }

129
130
131

  void MeshStructure::addAlongSide(Element *el, GeoIndex subObj, int ithObj, 
				   int elType, bool reverseOrder)
132
  {
133
134
135
136
    FUNCNAME("MeshStructure::addAlongSide()");

    if (debugMode) {
      MSG("addAlondSide(%d, %d, %d, %d)\n",
137
	  el->getIndex(), ithObj, elType, reverseOrder);
138
139
140
      MSG("Element is leaf: %d\n", el->isLeaf());
    }
    
141
142
143
    insertElement(el->isLeaf());
    
    if (!el->isLeaf()) {
144
145
      int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
      int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
146

147
148
149
150
151
152
153
154
      if (debugMode) {
	MSG("Child index %d  %d\n", 
	    el->getFirstChild()->getIndex(),
	    el->getSecondChild()->getIndex());
	MSG("s1 = %d    s2 = %d\n", s1, s2);
	MSG("   \n");
      }

155
156
      if (!reverseOrder) {
	if (s1 != -1) 
157
158
	  addAlongSide(el->getFirstChild(), subObj, s1, 
		       el->getChildType(elType), reverseOrder);
159
	if (s2 != -1)
160
161
	  addAlongSide(el->getSecondChild(), subObj, s2, 
		       el->getChildType(elType), reverseOrder);
162
163
      } else {
	if (s2 != -1)
164
165
	  addAlongSide(el->getSecondChild(), subObj, s2, 
		       el->getChildType(elType), reverseOrder);
166
	if (s1 != -1) 
167
168
	  addAlongSide(el->getFirstChild(), subObj, s1, 
		       el->getChildType(elType), reverseOrder);
169
      }
170
171
172
    }
  }

173

174
175
176
177
178
  void MeshStructure::reset() 
  {
    currentIndex = 0;
    pos = 0;
    currentElement = 0;
179
180
181
182
183

    if (code.size() > 0)
      currentCode = code[0];
    else 
      currentCode = 0;
184
185
  }

186

187
188
  bool MeshStructure::nextElement(MeshStructure *insert)
  {
189
    FUNCNAME_DBG("MeshStructure::nextElement()");
190
191

    if (insert)
192
193
      insert->insertElement(isLeafElement());

194
195
    pos++;
    currentElement++;
196

197
198
    if (currentElement >= nElements) 
      return false;
199

200
    if (pos >= structureSize) {
201
202
203
204
205
      currentIndex++;
      TEST_EXIT_DBG(currentIndex < static_cast<int>(code.size()))
	("End of structure reached!\n");
      pos = 0;
      currentCode = code[currentIndex];
206
    } else {
207
      currentCode >>= 1;
208
    }
209

210
211
212
    return true;
  }

213

214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
  int MeshStructure::lookAhead(unsigned int n)
  {
    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;
  }


242
243
  bool MeshStructure::skipBranch(MeshStructure *insert)
  {
244
    FUNCNAME_DBG("MeshStructure::skipBranch()");
245

246
    if (isLeafElement()) {
247
248
249
250
      return nextElement(insert);
    } else {
      bool cont = nextElement(insert);
      cont = skipBranch(insert); // left branch
251
      TEST_EXIT_DBG(cont)("Invalid structure!\n");
252
253
254
255
256
      cont = skipBranch(insert); // righ branch
      return cont;
    }
  }

257

258
259
260
261
  void MeshStructure::merge(MeshStructure *structure1,
			    MeshStructure *structure2,
			    MeshStructure *result)
  {
262
    FUNCNAME_DBG("MeshStructure::merge()");
263

264
265
266
267
268
    result->clear();
    structure1->reset();
    structure2->reset();

    bool cont = true;
269
    while (cont) {
270
271
272
273
      bool cont1;
#if DEBUG != 0
      bool cont2;
#endif
274
      if (structure1->isLeafElement() == structure2->isLeafElement()) {
275
	cont1 = structure1->nextElement(result);
276
#if DEBUG != 0
277
	cont2 = structure2->nextElement();
278
#endif
279
      } else {
280
	if (structure1->isLeafElement()) {
281
	  cont1 = structure1->nextElement();
282
#if DEBUG != 0
283
	  cont2 = structure2->skipBranch(result);
284
#endif
285
286
	} else {
	  cont1 = structure1->skipBranch(result);
287
#if DEBUG != 0
288
	  cont2 = structure2->nextElement();
289
#endif
290
291
	}
      }
292
      TEST_EXIT_DBG(cont1 == cont2)("Structures don't match!\n");
293
294
295
296
297
298
      cont = cont1;
    }

    result->commit();
  }

299

300
301
  void MeshStructure::fitMeshToStructure(Mesh *mesh,
					 RefinementManager *manager,
302
					 bool debugMode,
Thomas Witkowski's avatar
Thomas Witkowski committed
303
304
					 int macroElIndex,
					 bool ignoreFinerMesh) 
305
306
307
308
309
310
311
  {
    FUNCNAME("MeshStructure::fitMeshToStructure()");

    bool cont = true;

    // decorate leaf data
    reset();
312
    TraverseStack stack;
313
    ElInfo *elInfo = nullptr;
314
315
316
317
    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);
318

319

320
    while (elInfo) {
321
322
      Element *element = elInfo->getElement();

323
324
      TEST_EXIT(cont)("unexpected structure code end!\n");

325
      if (isLeafElement()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
326
327
328
329
330
331
	if (ignoreFinerMesh && !element->isLeaf()) {
	  int level = elInfo->getLevel();
	  while (elInfo && level >= elInfo->getLevel())
	    elInfo = stack.traverseNext(elInfo);
	} else {
	  TEST_EXIT(element->isLeaf())
332
333
	    ("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
334
335
	}
      } 
336

337
338
      TEST_EXIT_DBG(element)("Should not happen!\n");

339
      if (element->isLeaf() && !isLeafElement()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
340
	MeshStructure *structure = new MeshStructure();
341
342
343
	cont = skipBranch(structure);
	structure->commit();

344
345
346
	MeshStructure_ED *elData = new MeshStructure_ED(element->getElementData());
	elData->setStructure(structure);
	element->setElementData(elData);
347
348
349
350
      } else {
	cont = nextElement();
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
351
352
      if (elInfo)
	elInfo = stack.traverseNext(elInfo);
353
354
355
    }

    // refine mesh
356
    bool finished = true;
357

358
    do {
359
      finished = true;
360
361
362
363
      if (macroElIndex == -1)
	elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
      else
	elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1, Mesh::CALL_LEAF_EL);      
364
      while (elInfo) {
365
	Element *element = elInfo->getElement();
366
	if (element->getElementData(MESH_STRUCTURE) != nullptr) {
367
368
369
370
371
372
373
	  element->setMark(1);
	  finished = false;
	} else {
	  element->setMark(0);
	}
	elInfo = stack.traverseNext(elInfo);
      }
374
375
376
377
378

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

380
381
382
383
	if (macroElIndex == -1)
	  manager->refineMesh(mesh);
	else
	  manager->refineMacroElement(mesh, macroElIndex);
384

385
386
387
388
389
#if (DEBUG != 0)
	TEST_EXIT(oldMeshIndex != mesh->getChangeIndex())
	  ("Mesh has not been changed by refinement procedure!\n");
#endif
      }
390
    } while (!finished);
391
392
  }

393

Thomas Witkowski's avatar
Thomas Witkowski committed
394
  string MeshStructure::toStr(bool resetCode)
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
  {
    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
414
415
416
417
418
419
420
421
422
423

    return oss.str();
  }


  void MeshStructure::print(bool resetCode)
  {
    FUNCNAME("MeshStructure::print()");
    
    string str = toStr(resetCode);
424
   
Thomas Witkowski's avatar
Thomas Witkowski committed
425
426
    if (str.length() < 255) {
      MSG("Mesh structure code: %s\n", str.c_str());
427
428
    } else {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
429
      std::cout << "[" << MPI::COMM_WORLD.Get_rank() << "]                Mesh structure code: " << str << "\n";
430
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
431
      std::cout << "                Mesh structure code: " << str << "\n";
432
433
#endif
    }
434
435
436
  }


437
438
439
440
  bool MeshStructure::compare(MeshStructure &other)
  {
    return (other.getCode() == code);
  }
441
442


443
  void MeshStructure::getMeshStructureValues(int macroElIndex,
444
445
446
					     const DOFVector<double>* vec,
					     std::vector<double>& values)
  {
447
    FUNCNAME_DBG("MeshStructure::getMeshStructureValues()");
448

449
    TEST_EXIT_DBG(vec)("No DOFVector defined!\n");
450
451
452
  
    const FiniteElemSpace *feSpace = vec->getFeSpace();
    Mesh *mesh = feSpace->getMesh();
453
    int nVertexPreDofs = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
454
    bool feSpaceHasNonVertexDofs = (feSpace->getBasisFcts()->getDegree() > 1);
455
    values.clear();
456
457
458
459
#if (DEBUG != 0)    
    // In debug mode we add the macro element index to the value code.
    values.push_back(static_cast<double>(macroElIndex));
#endif
460

461
    ElementDofIterator elDofIter(feSpace);
462
463
464
465
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1, 
						 Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
466
      // For the macro element the mesh structure code stores all vertex values.
467
      if (elInfo->getLevel() == 0)
468
	for (int i = 0; i < mesh->getGeo(VERTEX); i++)
469
	  values.push_back((*vec)[elInfo->getElement()->getDof(i, nVertexPreDofs)]);
470

471
472
473
474
475
      if (!elInfo->getElement()->isLeaf()) {
	// If no leaf element store the value of the "new" DOF that is created
	// by bisectioning of this element.

	DegreeOfFreedom dof0 = 
476
	  elInfo->getElement()->getChild(0)->getDof(mesh->getDim(), nVertexPreDofs);
477
478
479
480
481
482
483
484
485
486
487
488
489
	values.push_back((*vec)[dof0]);
      } else {
	// If leaf element store all non vertex values of this element, thus
	// only relevant for higher order basis functions.

	if (feSpaceHasNonVertexDofs) {
	  elDofIter.reset(elInfo->getElement());
	  do {
	    if (elDofIter.getPosIndex() != VERTEX) 
	      values.push_back((*vec)[elDofIter.getDof()]);
	  } while (elDofIter.next());
	}
      }
490
491
492
493
494
495

      elInfo = stack.traverseNext(elInfo);
    }
  }


496
  void MeshStructure::setMeshStructureValues(int macroElIndex,
497
498
499
					     DOFVector<double>* vec,
					     const std::vector<double>& values)
  {
500
    FUNCNAME_DBG("MeshStructure::setMeshStructureValues()");
501

502
    TEST_EXIT_DBG(vec)("No DOFVector defined!\n");
503

504
505
506
    const FiniteElemSpace *feSpace = vec->getFeSpace();
    Mesh *mesh = feSpace->getMesh();
    bool feSpaceHasNonVertexDofs = (feSpace->getBasisFcts()->getDegree() > 1);
507
    int nVertexPreDofs = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
508
    unsigned int counter = 0;
509
510
511
512
513
514
#if (DEBUG != 0)    
    TEST_EXIT(static_cast<int>(values[0]) == macroElIndex)
      ("Value structure code was created for macro element index %d, but should be set to macro element index %d\n",
       static_cast<int>(values[0]), macroElIndex);
    counter++;
#endif
515

516
517
518
519
    TEST_EXIT_DBG(static_cast<int>(values.size()) >= mesh->getGeo(VERTEX))
      ("Should not happen!\n");

    ElementDofIterator elDofIter(feSpace);
520
521
522
523
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1, 
						 Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
524
      // For the macro element all vertex nodes are set first.
525
      if (elInfo->getLevel() == 0)
526
	for (int i = 0; i < mesh->getGeo(VERTEX); i++)
527
528
	  (*vec)[elInfo->getElement()->getDof(i, nVertexPreDofs)] = 
	    values[counter++];
529

530
      if (!elInfo->getElement()->isLeaf()) {
531
532
	// If no leaf element set the value of the "new" DOF that is created
	// by bisectioning of this element.
533
534
	TEST_EXIT_DBG(counter < values.size())("Should not happen!\n");
	
535
	(*vec)[elInfo->getElement()->getChild(0)->getDof(mesh->getDim(), nVertexPreDofs)] =
536
	  values[counter++];      
537
538
539
540
541
542
543
544
545
546
547
      } else {
	// On leaf elements set all non vertex values (thus DOFs of higher order
	// basis functions).

	if (feSpaceHasNonVertexDofs) {
	  elDofIter.reset(elInfo->getElement());
	  do {
	    if (elDofIter.getPosIndex() != VERTEX) 
	      (*vec)[elDofIter.getDof()] = values[counter++];
	  } while (elDofIter.next());
	}
548
549
550
551
552
      }

      elInfo = stack.traverseNext(elInfo);
    }
      
553
    TEST_EXIT_DBG(values.size() == counter)("Should not happen!\n");
554
555
  }

556
}