MeshManipulation.cc 16.8 KB
Newer Older
1
2
3
4
5
6
7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9
10
11
12
13
14
15
16
17
 * 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.
18
 *
19
 ******************************************************************************/
20
21


22
#include "parallel/MeshManipulation.h"
23
#include "DOFVector.h"
24
#include "Mesh.h"
25
#include "MeshStructure.h"
26
#include "BasisFunction.h"
27
#include "Traverse.h"
28
#include "Debug.h"
29
#include "FiniteElemSpace.h"
30

31
32
using namespace std;

33
namespace AMDiS { namespace Parallel {
34

35
36
  MeshManipulation::MeshManipulation(Mesh *m)
    : mesh(m)
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
  {
    switch (mesh->getDim()) {
    case 2:
      refineManager = new RefinementManager2d();
      break;
    case 3:
      refineManager = new RefinementManager3d();
      break;
    default:
      ERROR_EXIT("invalid dim!\n");
    }
  }


  MeshManipulation::~MeshManipulation()
  {
    delete refineManager;
  }

56
57
58
59
60
61
62
63
  void MeshManipulation::deleteDoubleDofs(std::vector<const FiniteElemSpace*>& feSpaces,
			  std::vector<MacroElement*>& newMacroEl,
			  ElementObjectDatabase &elObjDb)
  {
    std::set<MacroElement*> newMacroElSet(newMacroEl.begin(), newMacroEl.end());
    deleteDoubleDofs(feSpaces, newMacroElSet, elObjDb);
  }

64

65
  void MeshManipulation::deleteDoubleDofs(vector<const FiniteElemSpace*>& feSpaces,
66
					  std::set<MacroElement*>& newMacroEl,
67
					  ElementObjectDatabase &elObjDb)
68
  {
69
70
    FUNCNAME("MeshManipulation::deleteDoubleDofs()");

71
    // Search for the FE space with the highest degree of polynomials. Using this
72
    // FE space ensures that deleting DOFs defined on it, also DOFs of lower
73
74
75
    // order FE spaces will be deleted correctly.
    const FiniteElemSpace *feSpace = FiniteElemSpace::getHighest(feSpaces);

76
    // Create data structure that maps macro element indices to the
77
    // corresponding pointers.
78
    map<int, MacroElement*> macroIndexMap;
79
80
81
    for (std::set<MacroElement*>::iterator it = newMacroEl.begin();
	 it != newMacroEl.end(); ++it)
      macroIndexMap[(*it)->getIndex()] = *it;
82

83
    // === Traverse mesh and put all "old" macro element to macrosProcessed  ===
84
    // === that stores all macro elements which are really connected to the  ===
85
    // === overall mesh structure.                                           ===
86

87
    std::set<int> macrosProcessed;
88
    TraverseStack stack;
89
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
90
    while (elInfo) {
91
92
93
94
95
96
97
      if (newMacroEl.count(elInfo->getMacroElement()) == 0) {
	int index = elInfo->getMacroElement()->getIndex();

	macrosProcessed.insert(index);
	macroIndexMap[index] = elInfo->getMacroElement();
      }

98
99
100
      elInfo = stack.traverseNext(elInfo);
    }

101
#ifndef NDEBUG
102
    DOFVector<WorldVector<double> > coords(feSpace, "dofCorrds");
Thomas Witkowski's avatar
a    
Thomas Witkowski committed
103
    feSpace->getMesh()->getDofIndexCoords(coords);
104
105
#endif

106
107
108
109
    // === Traverse all new elements and connect them to the overall mesh     ===
    // === structure by deleting all double DOFs on macro element's vertices, ===
    // === edges and faces.                                                   ===

110
111
112
    map<const DegreeOfFreedom*, const DegreeOfFreedom*> mapDelDofs;
    map<const DegreeOfFreedom*, GeoIndex> dofPosIndex;

113
    for (std::set<MacroElement*>::iterator it = newMacroEl.begin();
114
	 it != newMacroEl.end(); ++it) {
115

116
117
      // === Traverse all vertices of the new element. ===

118
      for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
119
	vector<ElementObjectData> &vertexEl =
120
	  elObjDb.getElementsVertex((*it)->getIndex(), i);
121

122
	for (vector<ElementObjectData>::iterator elIt = vertexEl.begin();
123
124
125
	     elIt != vertexEl.end(); ++elIt) {
	  if (elIt->elIndex == (*it)->getIndex())
	    continue;
126

127
128
	  if (macrosProcessed.count(elIt->elIndex)) {
	    TEST_EXIT_DBG(macroIndexMap.count(elIt->elIndex))
129
130
131
132
133
134
135
136
137
	      ("Should not happen!\n");

	    Element *el0 = (*it)->getElement();
	    Element *el1 = macroIndexMap[elIt->elIndex]->getElement();

	    const DegreeOfFreedom *dof0 = el0->getDof(i);
	    const DegreeOfFreedom *dof1 = el1->getDof(elIt->ithObject);

	    mapDelDofs[dof0] = dof1;
138
139
	    dofPosIndex[dof0] = VERTEX;

140
	    break;
141
	  }
142
143
144
145
146
147
	}
      }

      for (int i = 0; i < mesh->getGeo(EDGE); i++) {
	ElementObjectData elObj((*it)->getIndex(), i);

148
	vector<ElementObjectData> &edgeEl =
149
	  elObjDb.getElementsEdge((*it)->getIndex(), i);
150

151
      	for (vector<ElementObjectData>::iterator elIt = edgeEl.begin();
152
153
154
155
	     elIt != edgeEl.end(); ++elIt) {
	  if (elIt->elIndex == (*it)->getIndex())
	    continue;

156
	  if (macrosProcessed.count(elIt->elIndex)) {
157
158
	    TEST_EXIT_DBG(macroIndexMap.count(elIt->elIndex))
	      ("Should not happen!\n");
159
160

	    Element *el0 = (*it)->getElement();
161
162
	    Element *el1 = macroIndexMap[elIt->elIndex]->getElement();

163
164
	    TEST_EXIT_DBG(el0->getMesh() == el1->getMesh())("Mesh is different.\n");

165
	    bool reverseMode = elObjDb.getEdgeReverseMode(elObj, *elIt);
166
167

	    BoundaryObject b0(el0, 0, EDGE, i, reverseMode);
168
169
170
	    BoundaryObject b1(el1, 0, EDGE, elIt->ithObject, false);

	    DofContainer dofs0, dofs1;
171
172
173
	    vector<GeoIndex> dofGeoIndex0, dofGeoIndex1;
	    el0->getAllDofs(feSpace, b0, dofs0, true, &dofGeoIndex0);
	    el1->getAllDofs(feSpace, b1, dofs1, true, &dofGeoIndex1);
174
175


176
#ifndef NDEBUG
177
	    if (feSpaces.size())
178
	      debug::testDofsByCoords(coords, dofs0, dofs1);
179
	    else
180
181
	      TEST_EXIT_DBG(dofs0.size() == dofs1.size())
		("Should not happen!\n");
182
#endif
183
 	    for (unsigned int i = 0; i < dofs0.size(); i++) {
184
 	      mapDelDofs[dofs0[i]] = dofs1[i];
185
186
187
	      dofPosIndex[dofs0[i]] = dofGeoIndex0[i];

	      TEST_EXIT_DBG(dofGeoIndex0[i] == dofGeoIndex1[i])
188
		("Should not happen: %d %d\n", dofGeoIndex0[i], dofGeoIndex1[i]);
189
	    }
190
191

	    break;
192
193
194
195
	  }
	}
      }

196

197
198
199
      for (int i = 0; i < mesh->getGeo(FACE); i++) {
	ElementObjectData elObj((*it)->getIndex(), i);

200
	vector<ElementObjectData> &faceEl =
201
	  elObjDb.getElementsFace((*it)->getIndex(), i);
202
203
204
205
206
207

	for (vector<ElementObjectData>::iterator elIt = faceEl.begin();
	     elIt != faceEl.end(); ++elIt) {
	  if (elIt->elIndex == (*it)->getIndex())
	    continue;

208
	  if (macrosProcessed.count(elIt->elIndex)) {
209
210
	    TEST_EXIT_DBG(macroIndexMap.count(elIt->elIndex))
	      ("Should not happen!\n");
211

212
	    Element *el0 = (*it)->getElement();
213
214
	    Element *el1 = macroIndexMap[elIt->elIndex]->getElement();

215
	    bool reverseMode = elObjDb.getFaceReverseMode(elObj, *elIt);
216
217
218
219

	    BoundaryObject b0(el0, 0, FACE, i, reverseMode);
	    BoundaryObject b1(el1, 0, FACE, elIt->ithObject, false);

220
221
222
223
	    DofContainer dofs0, dofs1;
	    vector<GeoIndex> dofGeoIndex0, dofGeoIndex1;
	    el0->getAllDofs(feSpace, b0, dofs0, true, &dofGeoIndex0);
	    el1->getAllDofs(feSpace, b1, dofs1, true, &dofGeoIndex1);
224

225
#ifndef NDEBUG
226
	    if (feSpaces.size())
227
	      debug::testDofsByCoords(coords, dofs0, dofs1);
228
	    else
229
230
	      TEST_EXIT_DBG(dofs0.size() == dofs1.size())
		("Should not happen!\n");
231
232
#endif

233
 	    for (unsigned int i = 0; i < dofs0.size(); i++) {
234
 	      mapDelDofs[dofs0[i]] = dofs1[i];
235
236
237
	      dofPosIndex[dofs0[i]] = dofGeoIndex1[i];

	      TEST_EXIT_DBG(dofGeoIndex0[i] == dofGeoIndex1[i])
238
		("Should not happen: %d %d\n",
239
		 dofGeoIndex0[i], dofGeoIndex1[i]);
240
	    }
241
242
243
244
245

	    break;
	  }
	}
      }
246
      macrosProcessed.insert((*it)->getIndex());
247
248
    }

249

250
251
    // === Remove all DOF replacments of the form A -> B, B -> C by A -> C. ===

252
253
254
    bool changed = false;
    do {
      changed = false;
255
      for (map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator it = mapDelDofs.begin();
256
	   it != mapDelDofs.end(); ++it) {
257
	map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator findIt = mapDelDofs.find(it->second);
258
	if (findIt != mapDelDofs.end()) {
259
260
261
	  TEST_EXIT_DBG(it->first != findIt->second)
	    ("Cycle %d -> %d and %d -> %d! Should not happen!\n",
	     *(it->first), *(it->second), *(findIt->first), *(findIt->second));
262
	  it->second = findIt->second;
263
264
	  changed = true;
	}
265
      }
266
267
268
    } while (changed);


269
270
    // === Set new DOF pointers in all elements of the mesh. ===

271
    elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
272
    while (elInfo) {
273
      for (int i = 0; i < mesh->getNumberOfNodes(); i++)
274
	if (mapDelDofs.count(elInfo->getElement()->getDof(i)) == 1)
275
	  elInfo->getElement()->setDof(i, const_cast<DegreeOfFreedom*>(mapDelDofs[elInfo->getElement()->getDof(i)]));
276
277
278
279

      elInfo = stack.traverseNext(elInfo);
    }

280

281
282
    // === And delete all double DOFs. ===

283
    for (map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator it = mapDelDofs.begin();
284
	 it != mapDelDofs.end(); ++it)
285
      mesh->freeDof(const_cast<DegreeOfFreedom*>(it->first),
286
		    dofPosIndex[it->first]);
287
288
  }

289

290
291
  bool MeshManipulation::fitElementToMeshCode(MeshStructure &code,
					      BoundaryObject &boundEl)
292
  {
293
    FUNCNAME("MeshManipulation::fitElementToMeshCode()");
294

295
    TEST_EXIT_DBG(boundEl.el)("No element given!\n");
296
297
298
299
300
301
302
303
304

    // If the code is empty, the element does not matter and the function can
    // return without chaning the mesh.
    if (code.empty())
      return false;

    // s0 and s1 are the number of the edge/face in both child of the element,
    // which contain the edge/face the function has to traverse through. If the
    // edge/face is not contained in one of the children, s0 or s1 is -1.
305
    int s0 = boundEl.el->getSubObjOfChild(0, boundEl.subObj,
306
					  boundEl.ithObj, boundEl.elType);
307
    int s1 = boundEl.el->getSubObjOfChild(1, boundEl.subObj,
308
					  boundEl.ithObj, boundEl.elType);
309
310
311
312

    TEST_EXIT_DBG(s0 != -1 || s1 != -1)("This should not happen!\n");

    bool meshChanged = false;
313
    Flag traverseFlag =
314
315
316
317
      Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND;

    // Test for reverse mode, in which the left and right children of elements
    // are flipped.
318
    if (boundEl.reverseMode)
319
      traverseFlag |= Mesh::CALL_REVERSE_MODE;
320
321
322
323


    // === If the edge/face is contained in both children. ===

324
325
    Mesh *mesh = boundEl.el->getMesh();

326
327
328
329
    if (s0 != -1 && s1 != -1) {
      // Create traverse stack and traverse within the mesh until the element,
      // which should be fitted to the mesh structure code, is reached.
      TraverseStack stack;
330
#ifndef NDEBUG
331
      ElInfo *elInfo =
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
332
	stack.traverseFirstOneMacro(mesh, boundEl.elIndex, -1, traverseFlag);
333

334
335
336
337
338
      TEST_EXIT(elInfo->getElement() == boundEl.el)
	("This should not happen!\n");
#else
      // stack must be initialized before passed to fitElementToMeshCode()
      stack.traverseFirstOneMacro(mesh, boundEl.elIndex, -1, traverseFlag);
339
#endif
340
341


342

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
343
344
345
346
      pcode = &code;
      pstack = &stack;
      rMode = boundEl.reverseMode;
      return fitElementToMeshCode(boundEl.subObj, boundEl.ithObj);
347
348
349
350
351
    }


    // === The edge/face is contained in only on of the both children. ===

352
    if (boundEl.el->isLeaf()) {
353
354
355

      // Create traverse stack and traverse the mesh to the element.
      TraverseStack stack;
356
      ElInfo *elInfo =
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
357
	stack.traverseFirstOneMacro(mesh, boundEl.elIndex, -1, traverseFlag);
358
359
360
361

      TEST_EXIT_DBG(elInfo)("This should not happen!\n");

      // Code is not leaf, therefore refine the element.
362
363
      boundEl.el->setMark(1);
      refineManager->setMesh(boundEl.el->getMesh());
364
365
366
      refineManager->setStack(&stack);
      refineManager->refineFunction(elInfo);
      meshChanged = true;
367
368
      // If element is leaf and code contains only one leaf element, we are finished.
      if (code.getNumElements() == 1 && code.isLeafElement())
369
	return true;
370
371


372
373
    }

374
375
376
    Element *child0 = boundEl.el->getFirstChild();
    Element *child1 = boundEl.el->getSecondChild();
    if (boundEl.reverseMode) {
377
      using std::swap;
378
      swap(s0, s1);
379
      swap(child0, child1);
380
381
382
383
384
385
386
    }

    // === We know that the edge/face is contained in only one of the children. ===
    // === Therefore, traverse the mesh to this children and fit this element   ===
    // === To the mesh structure code.                                          ===

    TraverseStack stack;
387
    ElInfo *elInfo =
388
      stack.traverseFirstOneMacro(mesh, boundEl.elIndex, -1, traverseFlag);
389
390
391

    if (s0 != -1) {
      while (elInfo && elInfo->getElement() != child0)
392
	elInfo = stack.traverseNext(elInfo);
393

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
394
395
396
397
      pcode = &code;
      pstack = &stack;
      rMode = boundEl.reverseMode;
      meshChanged |= fitElementToMeshCode(boundEl.subObj, s0);
398
    } else {
399
400
      while (elInfo && elInfo->getElement() != child1)
	elInfo = stack.traverseNext(elInfo);
401

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
402
403
404
405
      pcode = &code;
      pstack = &stack;
      rMode = boundEl.reverseMode;
      meshChanged |= fitElementToMeshCode(boundEl.subObj, s1);
406
407
408
409
410
411
412
    }


    return meshChanged;
  }


Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
413
  bool MeshManipulation::fitElementToMeshCode(GeoIndex subObj, int ithObj)
414
415
416
417
418
  {
    FUNCNAME("MeshManipulation::fitElementToMeshCode()");

    // === Test if there are more elements in stack to check with the code. ===

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
419
    ElInfo *elInfo = pstack->getElInfo();
420
421
422
423
424
425
426
427
    if (!elInfo)
      return false;


    // === Test if code contains a leaf element. If this is the case, the ===
    // === current element is finished. Traverse the mesh to the next     ===
    // === coarser element.                                               ===

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
428
    if (pcode->isLeafElement()) {
429
430
431
      int level = elInfo->getLevel();

      do {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
432
	elInfo = pstack->traverseNext(elInfo);
433
434
435
436
437
438
439
440
      } while (elInfo && elInfo->getLevel() > level);

      return false;
    }


    bool meshChanged = false;
    Element *el = elInfo->getElement();
441
442
    int s0 = el->getSubObjOfChild(0, subObj, ithObj, elInfo->getType());
    int s1 = el->getSubObjOfChild(1, subObj, ithObj, elInfo->getType());
443
444
445

    // === If element is leaf (and the code is not), refine the element. ===

446
447
    bool elementRefined = true;

448
449
450
    if (el->isLeaf()) {
      TEST_EXIT_DBG(elInfo->getLevel() < 255)("This should not happen!\n");

451
      // In some situations refinement of an element can be omitted, altough
452
453
454
455
456
457
458
459
460
461
462
463
      // it is defined in the mesh structure code. One case is when the
      // following holds:
      //   - the code is used to adapt a face of a tetrahedron
      //   - only one of the children of the current element contains this face
      //   - the face to be adapted of the current element is either the
      //     local face 0 or 1
      //   - the next structure code is 0, i.e., we would be finished after
      //     the refinement of this element
      // In this scenario, the element must be refined due to the structure
      // code, but the refinement does not introduce new DOFs on the face,
      // that should be adapted. Thus, we can ommit the refinement.
      if (subObj == FACE) {
464
  	if ((s0 != -1 && s1 == -1) || (s0 == -1 && s1 != -1)) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
465
  	  if (ithObj <= 1 && pcode->lookAhead() == 0) {
466
  	    elementRefined = false;
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
467
468
 	    pcode->nextElement();
 	    pstack->traverseNext(elInfo);
469
470
471
  	  }
  	}
      }
472

473
474
475
      if (elementRefined) {
	el->setMark(1);
	refineManager->setMesh(el->getMesh());
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
476
	refineManager->setStack(pstack);
477
478
	refineManager->refineFunction(elInfo);
	meshChanged = true;
479
      }
480
481
482
    }


483
484
485
    if (elementRefined) {
      // === Continue fitting the mesh structure code to the children of the ===
      // === current element.                                                ===
486

487
488
      Element *child0 = el->getFirstChild();
      Element *child1 = el->getSecondChild();
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
489
      if (rMode) {
490
	using std::swap;
491
492
493
	swap(s0, s1);
	swap(child0, child1);
      }
494

495

496
      // === Traverse left child. ===
497

498
499
500
      if (s0 != -1) {
	// The edge/face is contained in the left child, therefore fit this
	// child to the mesh structure code.
501

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
502
503
504
505
	pstack->traverseNext(elInfo);
	pcode->nextElement();
	meshChanged |= fitElementToMeshCode(subObj, s0);
	elInfo = pstack->getElInfo();
506
507
508
509
      } else {
	// The edge/face is not contained in the left child. Hence we need
	// to traverse through all subelements of the left child until we get
	// the second child of the current element.
510

511
	do {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
512
	  elInfo = pstack->traverseNext(elInfo);
513
514
	} while (elInfo && elInfo->getElement() != child1);

515
	TEST_EXIT_DBG(elInfo != NULL)("This should not happen!\n");
516
517
      }

518
519
520
      TEST_EXIT_DBG(elInfo->getElement() == child1)
	("Got wrong child with idx = %d! Searched for child idx = %d\n",
	 elInfo->getElement()->getIndex(), child1->getIndex());
521
522


523
      // === Traverse right child. ===
524

525
526
527
      if (s1 != -1) {
	// The edge/face is contained in the right child, therefore fit this
	// child to the mesh structure code.
528

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
529
530
	pcode->nextElement();
	meshChanged |= fitElementToMeshCode(subObj, s1);
531
532
533
534
      } else {
	// The edge/face is not contained in the right child. Hence we need
	// to traverse through all subelements of the right child until we have
	// finished traversing the current element with all its subelements.
535

536
	int level = elInfo->getLevel();
537

538
	do {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
539
	  elInfo = pstack->traverseNext(elInfo);
540
541
	} while (elInfo && elInfo->getLevel() > level);
      }
542

543
    }
544

545
546
    return meshChanged;
  }
547

548
} }