MeshStructure.cc 14.2 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
117
118
119
120
121
122
123
124
125
126
127
128
129
130
    /*    switch (bound.subObj)
    {
      case EDGE:
    */
	if (!el->isLeaf()) {
	  if (s1 == -1)
	    addAlongSide(el->getSecondChild(), bound.subObj, s2, 
			el->getChildType(bound.elType), bound.reverseMode);
	  else if (s2 == -1)
	    addAlongSide(el->getFirstChild(), bound.subObj, s1, 
			el->getChildType(bound.elType), bound.reverseMode);
	  else
	    addAlongSide(el, bound.subObj, bound.ithObj, bound.elType, bound.reverseMode);
	}
	  /*	break;
      case FACE:
131
	addAlongSide(el, bound.subObj, bound.ithObj, bound.elType, bound.reverseMode);
132
133
134
	break;
      default:
	ERROR_EXIT("What is this?\n");
135
    }
136
137
	  */
    
138
139
140
    commit();    
  }

141
142
143

  void MeshStructure::addAlongSide(Element *el, GeoIndex subObj, int ithObj, 
				   int elType, bool reverseOrder)
144
  {
145
146
147
148
    FUNCNAME("MeshStructure::addAlongSide()");

    if (debugMode) {
      MSG("addAlondSide(%d, %d, %d, %d)\n",
149
	  el->getIndex(), ithObj, elType, reverseOrder);
150
151
152
      MSG("Element is leaf: %d\n", el->isLeaf());
    }
    
153
154
155
    insertElement(el->isLeaf());
    
    if (!el->isLeaf()) {
156
157
      int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
      int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
158

159
160
161
162
163
164
165
166
      if (debugMode) {
	MSG("Child index %d  %d\n", 
	    el->getFirstChild()->getIndex(),
	    el->getSecondChild()->getIndex());
	MSG("s1 = %d    s2 = %d\n", s1, s2);
	MSG("   \n");
      }

167
168
      if (!reverseOrder) {
	if (s1 != -1) 
169
170
	  addAlongSide(el->getFirstChild(), subObj, s1, 
		       el->getChildType(elType), reverseOrder);
171
	if (s2 != -1)
172
173
	  addAlongSide(el->getSecondChild(), subObj, s2, 
		       el->getChildType(elType), reverseOrder);
174
175
      } else {
	if (s2 != -1)
176
177
	  addAlongSide(el->getSecondChild(), subObj, s2, 
		       el->getChildType(elType), reverseOrder);
178
	if (s1 != -1) 
179
180
	  addAlongSide(el->getFirstChild(), subObj, s1, 
		       el->getChildType(elType), reverseOrder);
181
      }
182
183
184
    }
  }

185

186
187
188
189
190
  void MeshStructure::reset() 
  {
    currentIndex = 0;
    pos = 0;
    currentElement = 0;
191
192
193
194
195

    if (code.size() > 0)
      currentCode = code[0];
    else 
      currentCode = 0;
196
197
  }

198

199
200
  bool MeshStructure::nextElement(MeshStructure *insert)
  {
201
    FUNCNAME_DBG("MeshStructure::nextElement()");
202
203

    if (insert)
204
205
      insert->insertElement(isLeafElement());

206
207
    pos++;
    currentElement++;
208

209
210
    if (currentElement >= nElements) 
      return false;
211

212
    if (pos >= structureSize) {
213
214
215
216
217
      currentIndex++;
      TEST_EXIT_DBG(currentIndex < static_cast<int>(code.size()))
	("End of structure reached!\n");
      pos = 0;
      currentCode = code[currentIndex];
218
    } else {
219
      currentCode >>= 1;
220
    }
221

222
223
224
    return true;
  }

225

226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
  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;
  }


254
255
  bool MeshStructure::skipBranch(MeshStructure *insert)
  {
256
    FUNCNAME_DBG("MeshStructure::skipBranch()");
257

258
    if (isLeafElement()) {
259
260
261
262
      return nextElement(insert);
    } else {
      bool cont = nextElement(insert);
      cont = skipBranch(insert); // left branch
263
      TEST_EXIT_DBG(cont)("Invalid structure!\n");
264
265
266
267
268
      cont = skipBranch(insert); // righ branch
      return cont;
    }
  }

269

270
271
272
273
  void MeshStructure::merge(MeshStructure *structure1,
			    MeshStructure *structure2,
			    MeshStructure *result)
  {
274
    FUNCNAME_DBG("MeshStructure::merge()");
275

276
277
278
279
280
    result->clear();
    structure1->reset();
    structure2->reset();

    bool cont = true;
281
    while (cont) {
282
283
284
285
      bool cont1;
#if DEBUG != 0
      bool cont2;
#endif
286
      if (structure1->isLeafElement() == structure2->isLeafElement()) {
287
	cont1 = structure1->nextElement(result);
288
#if DEBUG != 0
289
	cont2 = structure2->nextElement();
290
#endif
291
      } else {
292
	if (structure1->isLeafElement()) {
293
	  cont1 = structure1->nextElement();
294
#if DEBUG != 0
295
	  cont2 = structure2->skipBranch(result);
296
#endif
297
298
	} else {
	  cont1 = structure1->skipBranch(result);
299
#if DEBUG != 0
300
	  cont2 = structure2->nextElement();
301
#endif
302
303
	}
      }
304
      TEST_EXIT_DBG(cont1 == cont2)("Structures don't match!\n");
305
306
307
308
309
310
      cont = cont1;
    }

    result->commit();
  }

311

312
313
  void MeshStructure::fitMeshToStructure(Mesh *mesh,
					 RefinementManager *manager,
314
					 bool debugMode,
Thomas Witkowski's avatar
Thomas Witkowski committed
315
316
					 int macroElIndex,
					 bool ignoreFinerMesh) 
317
318
319
320
321
322
323
  {
    FUNCNAME("MeshStructure::fitMeshToStructure()");

    bool cont = true;

    // decorate leaf data
    reset();
324
    TraverseStack stack;
325
    ElInfo *elInfo = nullptr;
326
327
328
329
    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);
330

331

332
    while (elInfo) {
333
334
      Element *element = elInfo->getElement();

335
336
      TEST_EXIT(cont)("unexpected structure code end!\n");

337
      if (isLeafElement()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
338
339
340
341
342
343
	if (ignoreFinerMesh && !element->isLeaf()) {
	  int level = elInfo->getLevel();
	  while (elInfo && level >= elInfo->getLevel())
	    elInfo = stack.traverseNext(elInfo);
	} else {
	  TEST_EXIT(element->isLeaf())
344
345
	    ("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
346
347
	}
      } 
348

349
350
      TEST_EXIT_DBG(element)("Should not happen!\n");

351
      if (element->isLeaf() && !isLeafElement()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
352
	MeshStructure *structure = new MeshStructure();
353
354
355
	cont = skipBranch(structure);
	structure->commit();

356
357
358
	MeshStructure_ED *elData = new MeshStructure_ED(element->getElementData());
	elData->setStructure(structure);
	element->setElementData(elData);
359
360
361
362
      } else {
	cont = nextElement();
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
363
364
      if (elInfo)
	elInfo = stack.traverseNext(elInfo);
365
366
367
    }

    // refine mesh
368
    bool finished = true;
369

370
    do {
371
      finished = true;
372
373
374
375
      if (macroElIndex == -1)
	elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
      else
	elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1, Mesh::CALL_LEAF_EL);      
376
      while (elInfo) {
377
	Element *element = elInfo->getElement();
378
	if (element->getElementData(MESH_STRUCTURE) != nullptr) {
379
380
381
382
383
384
385
	  element->setMark(1);
	  finished = false;
	} else {
	  element->setMark(0);
	}
	elInfo = stack.traverseNext(elInfo);
      }
386
387
388
389
390

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

392
393
394
395
	if (macroElIndex == -1)
	  manager->refineMesh(mesh);
	else
	  manager->refineMacroElement(mesh, macroElIndex);
396

397
398
399
400
401
#if (DEBUG != 0)
	TEST_EXIT(oldMeshIndex != mesh->getChangeIndex())
	  ("Mesh has not been changed by refinement procedure!\n");
#endif
      }
402
    } while (!finished);
403
404
  }

405

Thomas Witkowski's avatar
Thomas Witkowski committed
406
  string MeshStructure::toStr(bool resetCode)
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
  {
    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
426
427
428
429
430
431
432
433
434
435

    return oss.str();
  }


  void MeshStructure::print(bool resetCode)
  {
    FUNCNAME("MeshStructure::print()");
    
    string str = toStr(resetCode);
436
   
Thomas Witkowski's avatar
Thomas Witkowski committed
437
438
    if (str.length() < 255) {
      MSG("Mesh structure code: %s\n", str.c_str());
439
440
    } else {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
441
      std::cout << "[" << MPI::COMM_WORLD.Get_rank() << "]                Mesh structure code: " << str << "\n";
442
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
443
      std::cout << "                Mesh structure code: " << str << "\n";
444
445
#endif
    }
446
447
448
  }


449
450
451
452
  bool MeshStructure::compare(MeshStructure &other)
  {
    return (other.getCode() == code);
  }
453
454


455
  void MeshStructure::getMeshStructureValues(int macroElIndex,
456
457
458
					     const DOFVector<double>* vec,
					     std::vector<double>& values)
  {
459
    FUNCNAME_DBG("MeshStructure::getMeshStructureValues()");
460

461
    TEST_EXIT_DBG(vec)("No DOFVector defined!\n");
462
463
464
  
    const FiniteElemSpace *feSpace = vec->getFeSpace();
    Mesh *mesh = feSpace->getMesh();
465
    int nVertexPreDofs = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
466
    bool feSpaceHasNonVertexDofs = (feSpace->getBasisFcts()->getDegree() > 1);
467
    values.clear();
468
469
470
471
#if (DEBUG != 0)    
    // In debug mode we add the macro element index to the value code.
    values.push_back(static_cast<double>(macroElIndex));
#endif
472

473
    ElementDofIterator elDofIter(feSpace);
474
475
476
477
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1, 
						 Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
478
      // For the macro element the mesh structure code stores all vertex values.
479
      if (elInfo->getLevel() == 0)
480
	for (int i = 0; i < mesh->getGeo(VERTEX); i++)
481
	  values.push_back((*vec)[elInfo->getElement()->getDof(i, nVertexPreDofs)]);
482

483
484
485
486
487
      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 = 
488
	  elInfo->getElement()->getChild(0)->getDof(mesh->getDim(), nVertexPreDofs);
489
490
491
492
493
494
495
496
497
498
499
500
501
	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());
	}
      }
502
503
504
505
506
507

      elInfo = stack.traverseNext(elInfo);
    }
  }


508
  void MeshStructure::setMeshStructureValues(int macroElIndex,
509
510
511
					     DOFVector<double>* vec,
					     const std::vector<double>& values)
  {
512
    FUNCNAME_DBG("MeshStructure::setMeshStructureValues()");
513

514
    TEST_EXIT_DBG(vec)("No DOFVector defined!\n");
515

516
517
518
    const FiniteElemSpace *feSpace = vec->getFeSpace();
    Mesh *mesh = feSpace->getMesh();
    bool feSpaceHasNonVertexDofs = (feSpace->getBasisFcts()->getDegree() > 1);
519
    int nVertexPreDofs = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
520
    unsigned int counter = 0;
521
522
523
524
525
526
#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
527

528
529
530
531
    TEST_EXIT_DBG(static_cast<int>(values.size()) >= mesh->getGeo(VERTEX))
      ("Should not happen!\n");

    ElementDofIterator elDofIter(feSpace);
532
533
534
535
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1, 
						 Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
536
      // For the macro element all vertex nodes are set first.
537
      if (elInfo->getLevel() == 0)
538
	for (int i = 0; i < mesh->getGeo(VERTEX); i++)
539
540
	  (*vec)[elInfo->getElement()->getDof(i, nVertexPreDofs)] = 
	    values[counter++];
541

542
      if (!elInfo->getElement()->isLeaf()) {
543
544
	// If no leaf element set the value of the "new" DOF that is created
	// by bisectioning of this element.
545
546
	TEST_EXIT_DBG(counter < values.size())("Should not happen!\n");
	
547
	(*vec)[elInfo->getElement()->getChild(0)->getDof(mesh->getDim(), nVertexPreDofs)] =
548
	  values[counter++];      
549
550
551
552
553
554
555
556
557
558
559
      } 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());
	}
560
561
562
563
564
      }

      elInfo = stack.traverseNext(elInfo);
    }
      
565
    TEST_EXIT_DBG(values.size() == counter)("Should not happen!\n");
566
567
  }

568
}