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

namespace AMDiS {

36
  const int MeshStructure::structureSize = 64;
37

38

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

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

54
55
    pos++;
    nElements++;
56
57
  }

58

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

68

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

73
    TraverseStack stack;
74

Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
75
76
77
78
79
80
81
82
83
    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());
84
85
86
87
88
89
      elInfo = stack.traverseNext(elInfo);
    }
  
    commit();
  }

90
91
92
93
94
95
96

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

    Element *el = bound.el;

97
98
    TEST_EXIT_DBG(el)("No element!\n");

99
100
    clear();

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

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

    commit();    
  }

127
128
129

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

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

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

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

171

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

    if (code.size() > 0)
      currentCode = code[0];
    else 
      currentCode = 0;
182
183
  }

184

185
186
  bool MeshStructure::nextElement(MeshStructure *insert)
  {
187
    FUNCNAME_DBG("MeshStructure::nextElement()");
188
189

    if (insert)
190
191
      insert->insertElement(isLeafElement());

192
193
    pos++;
    currentElement++;
194

195
196
    if (currentElement >= nElements) 
      return false;
197

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

208
209
210
    return true;
  }

211

212
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
  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;
  }


240
241
  bool MeshStructure::skipBranch(MeshStructure *insert)
  {
242
    FUNCNAME_DBG("MeshStructure::skipBranch()");
243

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

255

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

262
263
264
265
266
    result->clear();
    structure1->reset();
    structure2->reset();

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

    result->commit();
  }

297

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

    bool cont = true;

    // decorate leaf data
    reset();
310
    TraverseStack stack;
311
312
313
314
315
    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);
316

317

318
    while (elInfo) {
319
320
      Element *element = elInfo->getElement();

321
322
      TEST_EXIT(cont)("unexpected structure code end!\n");

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

335
336
      TEST_EXIT_DBG(element)("Should not happen!\n");

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

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

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

    // refine mesh
354
    bool finished = true;
355

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

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

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

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

391

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

    return oss.str();
  }


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


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


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

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

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

469
470
471
472
473
      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 = 
474
	  elInfo->getElement()->getChild(0)->getDof(mesh->getDim(), nVertexPreDofs);
475
476
477
478
479
480
481
482
483
484
485
486
487
	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());
	}
      }
488
489
490
491
492
493

      elInfo = stack.traverseNext(elInfo);
    }
  }


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

500
    TEST_EXIT_DBG(vec)("No DOFVector defined!\n");
501

502
503
504
    const FiniteElemSpace *feSpace = vec->getFeSpace();
    Mesh *mesh = feSpace->getMesh();
    bool feSpaceHasNonVertexDofs = (feSpace->getBasisFcts()->getDegree() > 1);
505
    int nVertexPreDofs = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
506
    unsigned int counter = 0;
507
508
509
510
511
512
#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
513

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

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

528
      if (!elInfo->getElement()->isLeaf()) {
529
530
	// If no leaf element set the value of the "new" DOF that is created
	// by bisectioning of this element.
531
532
	TEST_EXIT_DBG(counter < values.size())("Should not happen!\n");
	
533
	(*vec)[elInfo->getElement()->getChild(0)->getDof(mesh->getDim(), nVertexPreDofs)] =
534
	  values[counter++];      
535
536
537
538
539
540
541
542
543
544
545
      } 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());
	}
546
547
548
549
550
      }

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

554
}