MeshManipulation.cc 15.9 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
#include "parallel/MeshManipulation.h"
#include "Mesh.h"
15
#include "MeshStructure.h"
16
#include "BasisFunction.h"
17
#include "Traverse.h"
18
#include "Debug.h"
19
20
21

namespace AMDiS {

22
23
  MeshManipulation::MeshManipulation(Mesh *m)
    : mesh(m)
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
  {
    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;
  }


44
  void MeshManipulation::deleteDoubleDofs(vector<const FiniteElemSpace*>& feSpaces,
45
					  std::set<MacroElement*>& newMacroEl, 
46
					  ElementObjectDatabase &elObjDb)
47
  {
48
49
    FUNCNAME("MeshManipulation::deleteDoubleDofs()");

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

55
56
    // Create data structure that maps macro element indices to the 
    // corresponding pointers.
57
    map<int, MacroElement*> macroIndexMap;
58
59
60
61
    for (std::set<MacroElement*>::iterator it = newMacroEl.begin();
	 it != newMacroEl.end(); ++it)
      macroIndexMap[(*it)->getIndex()] = *it;

62
    // === Traverse mesh and put all "old" macro element to macrosProcessed  ===
63
    // === that stores all macro elements which are really connected to the  ===
64
    // === overall mesh structure.                                           ===
65

66
    std::set<int> macrosProcessed;   
67
    TraverseStack stack;
68
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
69
    while (elInfo) {
70
71
72
73
74
75
76
      if (newMacroEl.count(elInfo->getMacroElement()) == 0) {
	int index = elInfo->getMacroElement()->getIndex();

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

77
78
79
80
      elInfo = stack.traverseNext(elInfo);
    }


81
82
83
84
85
86
#if (DEBUG != 0)
    DOFVector<WorldVector<double> > coords(feSpace, "dofCorrds");
    feSpace->getMesh()->getDofIndexCoords(feSpace, coords);
#endif


87
88
89
90
    // === 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.                                                   ===

91
92
93
    map<const DegreeOfFreedom*, const DegreeOfFreedom*> mapDelDofs;
    map<const DegreeOfFreedom*, GeoIndex> dofPosIndex;

94
95
    for (std::set<MacroElement*>::iterator it = newMacroEl.begin(); 
	 it != newMacroEl.end(); ++it) {
96

97
98
      // === Traverse all vertices of the new element. ===

99
      for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
100
	vector<ElementObjectData> &vertexEl = 
101
	  elObjDb.getElementsVertex((*it)->getIndex(), i);
102

103
	for (vector<ElementObjectData>::iterator elIt = vertexEl.begin();
104
105
106
	     elIt != vertexEl.end(); ++elIt) {
	  if (elIt->elIndex == (*it)->getIndex())
	    continue;
107

108
109
110
111
112
113
114
115
116
117
118
	  if (macrosProcessed.count(elIt->elIndex) == 1) {
	    TEST_EXIT_DBG(macroIndexMap.count(elIt->elIndex) == 1)
	      ("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;
119
120
	    dofPosIndex[dof0] = VERTEX;

121
122
123
124
125
126
127
128
	    break;
	  } 
	}
      }

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

129
	vector<ElementObjectData> &edgeEl = 
130
	  elObjDb.getElementsEdge((*it)->getIndex(), i);
131
132

	for (vector<ElementObjectData>::iterator elIt = edgeEl.begin();
133
134
135
136
137
	     elIt != edgeEl.end(); ++elIt) {
	  if (elIt->elIndex == (*it)->getIndex())
	    continue;

	  if (macrosProcessed.count(elIt->elIndex) == 1) {
138
139
	    TEST_EXIT_DBG(macroIndexMap.count(elIt->elIndex))
	      ("Should not happen!\n");
140
	    
141
142
143
	    Element *el0 = (*it)->getElement();	    
	    Element *el1 = macroIndexMap[elIt->elIndex]->getElement();

144
	    bool reverseMode = elObjDb.getEdgeReverseMode(elObj, *elIt);
145
146

	    BoundaryObject b0(el0, 0, EDGE, i, reverseMode);
147
148
149
	    BoundaryObject b1(el1, 0, EDGE, elIt->ithObject, false);

	    DofContainer dofs0, dofs1;
150
151
152
	    vector<GeoIndex> dofGeoIndex0, dofGeoIndex1;
	    el0->getAllDofs(feSpace, b0, dofs0, true, &dofGeoIndex0);
	    el1->getAllDofs(feSpace, b1, dofs1, true, &dofGeoIndex1);
153
154

#if (DEBUG != 0)
155
	    if (feSpaces.size() == 1)
156
	      debug::testDofsByCoords(coords, dofs0, dofs1);
157
	    else 
158
159
	      TEST_EXIT_DBG(dofs0.size() == dofs1.size())
		("Should not happen!\n");
160
161
#endif

162
 	    for (unsigned int i = 0; i < dofs0.size(); i++) {
163
 	      mapDelDofs[dofs0[i]] = dofs1[i];
164
165
166
	      dofPosIndex[dofs0[i]] = dofGeoIndex0[i];

	      TEST_EXIT_DBG(dofGeoIndex0[i] == dofGeoIndex1[i])
167
		("Should not happen: %d %d\n", dofGeoIndex0[i], dofGeoIndex1[i]);
168
	    }
169
170

	    break;
171
172
173
174
	  }
	}
      }

175

176
177
178
179
      for (int i = 0; i < mesh->getGeo(FACE); i++) {
	ElementObjectData elObj((*it)->getIndex(), i);

	vector<ElementObjectData> &faceEl = 
180
	  elObjDb.getElementsFace((*it)->getIndex(), i);
181
182
183
184
185
186
187

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

	  if (macrosProcessed.count(elIt->elIndex) == 1) {
188
189
	    TEST_EXIT_DBG(macroIndexMap.count(elIt->elIndex))
	      ("Should not happen!\n");
190
191
192
193

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

194
	    bool reverseMode = elObjDb.getFaceReverseMode(elObj, *elIt);
195
196
197
198

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

199
200
201
202
	    DofContainer dofs0, dofs1;
	    vector<GeoIndex> dofGeoIndex0, dofGeoIndex1;
	    el0->getAllDofs(feSpace, b0, dofs0, true, &dofGeoIndex0);
	    el1->getAllDofs(feSpace, b1, dofs1, true, &dofGeoIndex1);
203
204

#if (DEBUG != 0)
205
	    if (feSpaces.size() == 1)
206
	      debug::testDofsByCoords(coords, dofs0, dofs1);
207
	    else 
208
209
	      TEST_EXIT_DBG(dofs0.size() == dofs1.size())
		("Should not happen!\n");
210
211
#endif

212
 	    for (unsigned int i = 0; i < dofs0.size(); i++) {
213
 	      mapDelDofs[dofs0[i]] = dofs1[i];
214
215
216
	      dofPosIndex[dofs0[i]] = dofGeoIndex1[i];

	      TEST_EXIT_DBG(dofGeoIndex0[i] == dofGeoIndex1[i])
217
218
		("Should not happen: %d %d\n", 
		 dofGeoIndex0[i], dofGeoIndex1[i]);
219
	    }
220
221
222
223
224

	    break;
	  }
	}
      }
225
226

      macrosProcessed.insert((*it)->getIndex());
227
228
    }

229

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

232
233
234
    bool changed = false;
    do {
      changed = false;
235
      for (map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator it = mapDelDofs.begin();
236
	   it != mapDelDofs.end(); ++it) {
237
	map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator findIt = mapDelDofs.find(it->second);
238
	if (findIt != mapDelDofs.end()) {
239
240
241
	  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));
242
243
244
245
246
247
248
	  it->second = findIt->second;   
	  changed = true;
	}
      } 
    } while (changed);


249
250
    // === Set new DOF pointers in all elements of the mesh. ===

251
    elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
252
    while (elInfo) {
253
      for (int i = 0; i < mesh->getNumberOfNodes(); i++)
254
	if (mapDelDofs.count(elInfo->getElement()->getDof(i)) == 1)
255
	  elInfo->getElement()->setDof(i, const_cast<DegreeOfFreedom*>(mapDelDofs[elInfo->getElement()->getDof(i)]));
256
257
258
259

      elInfo = stack.traverseNext(elInfo);
    }

260

261
262
    // === And delete all double DOFs. ===

263
    for (map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator it = mapDelDofs.begin();
264
	 it != mapDelDofs.end(); ++it)
265
266
      mesh->freeDof(const_cast<DegreeOfFreedom*>(it->first), 
		    dofPosIndex[it->first]);
267
268
  }

269

270
271
  bool MeshManipulation::fitElementToMeshCode(MeshStructure &code,
					      BoundaryObject &boundEl)
272
  {
273
    FUNCNAME("MeshManipulation::fitElementToMeshCode()");
274

275
    TEST_EXIT_DBG(boundEl.el)("No element given!\n");
276
277
278
279
280
281
282
283
284

    // 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.
285
286
287
288
    int s0 = boundEl.el->getSubObjOfChild(0, boundEl.subObj, 
					  boundEl.ithObj, boundEl.elType);
    int s1 = boundEl.el->getSubObjOfChild(1, boundEl.subObj, 
					  boundEl.ithObj, boundEl.elType);
289
290
291
292
293
294
295
296
297

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

    bool meshChanged = false;
    Flag traverseFlag = 
      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.
298
    if (boundEl.reverseMode)
299
300
301
302
303
      traverseFlag |= Mesh::CALL_REVERSE_MODE;    


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

304
305
    Mesh *mesh = boundEl.el->getMesh();

306
307
308
309
    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;
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
310
311
      ElInfo *elInfo = 
	stack.traverseFirstOneMacro(mesh, boundEl.elIndex, -1, traverseFlag);
312

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
313
314
      TEST_EXIT_DBG(elInfo->getElement() == boundEl.el)
	("This should not happen!\n");
315

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
316
317
318
319
      pcode = &code;
      pstack = &stack;
      rMode = boundEl.reverseMode;
      return fitElementToMeshCode(boundEl.subObj, boundEl.ithObj);
320
321
322
323
324
    }


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

325
    if (boundEl.el->isLeaf()) {
326
327
328
329
330
331
332

      // If element is leaf and code contains only one leaf element, we are finished.
      if (code.getNumElements() == 1 && code.isLeafElement())
	return false;     

      // Create traverse stack and traverse the mesh to the element.
      TraverseStack stack;
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
333
334
      ElInfo *elInfo = 
	stack.traverseFirstOneMacro(mesh, boundEl.elIndex, -1, traverseFlag);
335
336
337
338

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

      // Code is not leaf, therefore refine the element.
339
340
      boundEl.el->setMark(1);
      refineManager->setMesh(boundEl.el->getMesh());
341
342
343
344
345
      refineManager->setStack(&stack);
      refineManager->refineFunction(elInfo);
      meshChanged = true;
    }

346
347
348
    Element *child0 = boundEl.el->getFirstChild();
    Element *child1 = boundEl.el->getSecondChild();
    if (boundEl.reverseMode) {
349
350
351
352
353
354
355
356
357
      swap(s0, s1);
      swap(child0, child1);    
    }

    // === 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;
358
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
359
360
361
362
363

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

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
364
365
366
367
      pcode = &code;
      pstack = &stack;
      rMode = boundEl.reverseMode;
      meshChanged |= fitElementToMeshCode(boundEl.subObj, s0);
368
369
370
371
    } else {
      while (elInfo && elInfo->getElement() != child1) 
	elInfo = stack.traverseNext(elInfo);      

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
372
373
374
375
      pcode = &code;
      pstack = &stack;
      rMode = boundEl.reverseMode;
      meshChanged |= fitElementToMeshCode(boundEl.subObj, s1);
376
377
378
379
380
381
382
    }


    return meshChanged;
  }


Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
383
  bool MeshManipulation::fitElementToMeshCode(GeoIndex subObj, int ithObj)
384
385
386
387
388
  {
    FUNCNAME("MeshManipulation::fitElementToMeshCode()");

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

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
389
    ElInfo *elInfo = pstack->getElInfo();
390
391
392
393
394
395
396
397
    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
398
    if (pcode->isLeafElement()) {
399
400
401
      int level = elInfo->getLevel();

      do {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
402
	elInfo = pstack->traverseNext(elInfo);
403
404
405
406
407
408
409
410
      } while (elInfo && elInfo->getLevel() > level);

      return false;
    }


    bool meshChanged = false;
    Element *el = elInfo->getElement();
411
412
    int s0 = el->getSubObjOfChild(0, subObj, ithObj, elInfo->getType());
    int s1 = el->getSubObjOfChild(1, subObj, ithObj, elInfo->getType());
413
414
415

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

416
417
    bool elementRefined = true;

418
419
420
    if (el->isLeaf()) {
      TEST_EXIT_DBG(elInfo->getLevel() < 255)("This should not happen!\n");

421
422
423
424
425
426
427
428
429
430
431
432
433
434
      // In some situations refinement of an element can be omitted, altough 
      // 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) {
  	if (s0 != -1 && s1 == -1 || s0 == -1 && s1 != -1) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
435
  	  if (ithObj <= 1 && pcode->lookAhead() == 0) {
436
  	    elementRefined = false;
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
437
438
 	    pcode->nextElement();
 	    pstack->traverseNext(elInfo);
439
440
441
  	  }
  	}
      }
442

443
444
445
      if (elementRefined) {
	el->setMark(1);
	refineManager->setMesh(el->getMesh());
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
446
	refineManager->setStack(pstack);
447
448
449
	refineManager->refineFunction(elInfo);
	meshChanged = true;
      } 
450
451
452
    }


453
454
455
456
457
458
    if (elementRefined) {
      // === Continue fitting the mesh structure code to the children of the ===
      // === current element.                                                ===
      
      Element *child0 = el->getFirstChild();
      Element *child1 = el->getSecondChild();
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
459
      if (rMode) {
460
461
462
	swap(s0, s1);
	swap(child0, child1);
      }
463

464
465
466
467
468
469
470
    
      // === Traverse left child. ===
      
      if (s0 != -1) {
	// The edge/face is contained in the left child, therefore fit this
	// child to the mesh structure code.
	
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
471
472
473
474
	pstack->traverseNext(elInfo);
	pcode->nextElement();
	meshChanged |= fitElementToMeshCode(subObj, s0);
	elInfo = pstack->getElInfo();
475
476
477
478
479
480
      } 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.
	
	do {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
481
	  elInfo = pstack->traverseNext(elInfo);
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
	} while (elInfo && elInfo->getElement() != child1); 
	
	TEST_EXIT_DBG(elInfo != NULL)("This should not happen!\n");
      }  
      
      TEST_EXIT_DBG(elInfo->getElement() == child1)
	("Got wrong child with idx = %d! Searched for child idx = %d\n",
	 elInfo->getElement()->getIndex(), child1->getIndex());
      
      
      // === Traverse right child. ===
      
      if (s1 != -1) {
	// The edge/face is contained in the right child, therefore fit this
	// child to the mesh structure code.
	
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
498
499
	pcode->nextElement();
	meshChanged |= fitElementToMeshCode(subObj, s1);
500
501
502
503
504
505
506
507
      } 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.
	
	int level = elInfo->getLevel();
	
	do {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
508
	  elInfo = pstack->traverseNext(elInfo);
509
510
511
	} while (elInfo && elInfo->getLevel() > level);
      }
      
512
    }
513
      
514
515
    return meshChanged;
  }
516
  
517
}