CoarseningManager3d.cc 18.5 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
27
28
29
30
#include "CoarseningManager3d.h"
#include "Mesh.h"
#include "AdaptStationary.h"
#include "AdaptInstationary.h"
#include "Traverse.h"
#include "MacroElement.h"
#include "RCNeighbourList.h"
#include "FixVec.h"
#include "DOFIndexed.h"
31
#include "Debug.h"
32

33
34
using namespace std;

35
36
namespace AMDiS {

37
  void CoarseningManager3d::coarsenFunction(ElInfo *elInfo)
38
  {
39
40
41
#if HAVE_PARALLEL_DOMAIN_AMDIS
    FUNCNAME_DBG("CoarseningManager3d::coarsenFunction()");
#endif
42

43
    Tetrahedron *el = 
44
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(elInfo->getElement()));
45

46
    // If element must not be coarsend, return.
47
    if (el->getMark() >= 0) 
48
      return; 
49

50
    // Single leaves don't get coarsened.
51
    if (el->isLeaf()) 
52
      return;
53

54
    if (el->getChild(0)->getMark() >= 0 || el->getChild(1)->getMark() >= 0) {
55
      // One of the children must not be coarsend, so return.
56

57
      el->setMark(0);
58
      return;
59
60
61
    }

    if (!(el->getChild(0)->isLeaf()) || !(el->getChild(1)->isLeaf())) {
62
      // One of the children is not a leaf element. Try again later on.
63
      doMore = true;
64
      return;
65
    }
66

67
68
69
    DegreeOfFreedom *edge[2];
    int n_neigh, bound = 0;

70
71
72
    /****************************************************************************/
    /*  get a list for storing all elements at the coarsening edge and fill it  */
    /****************************************************************************/
73
74
    RCNeighbourList coarsenList(mesh->getMaxEdgeNeigh());
    coarsenList.setCoarseningManager(this);
75
76
77
78
79

    /****************************************************************************/
    /*  give the refinement edge the right orientation                          */
    /****************************************************************************/

80
81
82
    if (el->getDof(0, 0) < el->getDof(1, 0)) {
      edge[0] = const_cast<DegreeOfFreedom*>(el->getDof(0));
      edge[1] = const_cast<DegreeOfFreedom*>(el->getDof(1));
83
    } else {
84
85
      edge[1] = const_cast<DegreeOfFreedom*>(el->getDof(0));
      edge[0] = const_cast<DegreeOfFreedom*>(el->getDof(1));
86
87
    }

88
    coarsenList.setElement(0, el, true);
89
90
    n_neigh = 1;

91
    coarsenList.setOppVertex(0, 0, 0);
92
    coarsenList.setElType(0, elInfo->getType());
93
    bound = false;
94
95
    if (getCoarsenPatch(elInfo, edge, 0, coarsenList, &n_neigh)) {
      getCoarsenPatch(elInfo, edge, 1, coarsenList, &n_neigh);
96
97
      bound = true;
    }
98
99
100


#if HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
101
    vector<FixRefinementPatch::EdgeInEl> refineEdges;
102
    FixRefinementPatch::getOtherEl(mesh, stack, refineEdges);
Thomas Witkowski's avatar
Thomas Witkowski committed
103

104
105
106
    // === If the refinement edge must be fixed, add also the other part of this ===
    // === edge to the refinement patch.                                         ===
    
Thomas Witkowski's avatar
Thomas Witkowski committed
107
    if (refineEdges.size()) {
108
109
110
      // TODO: Remove these two lines and make something more meaningful!!
      el->setMark(0);
      return;
Thomas Witkowski's avatar
Thomas Witkowski committed
111
112

      Element *otherEl = refineEdges[0].first;
113
      TEST_EXIT_DBG(otherEl->getMesh() == mesh)("Something is wrong.\n");
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
      
      TraverseStack stack2;
      ElInfo *elInfo2 = 
	stack2.traverseFirstOneMacro(mesh, otherEl->getIndex(), -1, 
				     Mesh::CALL_EVERY_EL_PREORDER | 
				     Mesh::FILL_NEIGH |
				     Mesh::FILL_BOUND);
      
      bool foundEdge = false;
      
      while (elInfo2) {
	Element *el2 = elInfo2->getElement();
	for (int i = 0; i < 6; i++) {
	  DofEdge edge2 = elInfo2->getElement()->getEdge(i);
	  if (edge2.first == *(edge[0]) && edge2.second == *(edge[1]) && !el2->isLeaf()) {
	    
	    if (!el2->getChild(0)->isLeaf() || !el2->getChild(1)->isLeaf()) {
	      int edgeNo0 = el->getEdgeOfChild(0, i, elInfo2->getType());
	      int edgeNo1 = el->getEdgeOfChild(1, i, elInfo2->getType());
	      
	      bool refineChildFirst = 
135
		!((i > 0 && (edgeNo0 >= 0 && !el2->getChild(0)->isLeaf())) || 
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
		  (edgeNo1 >= 0 && !el2->getChild(1)->isLeaf()));
	      
	      if (refineChildFirst) {
		// TODO: WAS SOLL ICH NUR MACHEN???
		el->setMark(0);
		return;
	      }
	    } else {
	      coarsenList.setElType(n_neigh, elInfo2->getType());
	      coarsenList.setElement(n_neigh, elInfo2->getElement(), true);
	      n_neigh++;
	      
	      foundEdge = true;
	      
	      TraverseStack *tmpStack = stack;
	      stack = &stack2;
	      
	      if (getCoarsenPatch(elInfo2, edge, 0, coarsenList, &n_neigh)) {
		getCoarsenPatch(elInfo2, edge, 1, coarsenList, &n_neigh);
		bound = true;
	      }
	      stack = tmpStack;
	      break;
	    }
	  }
	}
	
	elInfo2 = stack2.traverseNext(elInfo2);
      }
      
      TEST_EXIT_DBG(foundEdge)("Should not happen!\n");
    }  
#endif
    
170
    coarsenList.fillNeighbourRelations(n_neigh, bound);
171
172
173
174
175
176
177
178
179

    /****************************************************************************/
    /*  check wether we can coarsen the patch or not                            */
    /****************************************************************************/

    // ==========================================================================
    // === check for periodic boundary ==========================================
    // ==========================================================================

180
    if (coarsenList.doCoarsePatch(n_neigh)) {
181
182
      int n_neigh_periodic;
      DegreeOfFreedom *next_edge[2];
183
      RCNeighbourList periodicList;
184

185
      while (edge[0] != NULL) {
186
	coarsenList.periodicSplit(edge, next_edge,
187
188
				  &n_neigh, &n_neigh_periodic,
				  periodicList);
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204

	coarsenPatch(periodicList, n_neigh_periodic, bound);

	edge[0] = next_edge[0];
	edge[1] = next_edge[1];
      }
    }
  }

  /*****************************************************************************/
  /*  coarsenTetrahedron:  coarses a single element of the coarsening patch;   */
  /*  dofs in the interior of the element are removed and dofs in the faces of */
  /*  the element are removed if neighbour has been coarsend or if the face    */
  /*  is part of the domains boundary                                          */
  /*****************************************************************************/

205
  void CoarseningManager3d::coarsenTetrahedron(RCNeighbourList &coarsenList, 
206
					       int index)
207
  {
208
    Tetrahedron *el = 
209
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(coarsenList.getElement(index)));
210
    Tetrahedron *child[2];
211

212
213
    child[0] = dynamic_cast<Tetrahedron*>(const_cast<Element*>(el->getChild(0)));
    child[1] = dynamic_cast<Tetrahedron*>(const_cast<Element*>(el->getChild(1)));
214
    int el_type = coarsenList.getType(index);
215
216
217
218
219
220

    /****************************************************************************/
    /*  Information about patch neighbours is still valid! But edge and face    */
    /*  dof's in a common face of patch neighbours have to be removed           */
    /****************************************************************************/

221
222
    for (int dir = 0; dir < 2; dir++) {
      Tetrahedron *neigh = 
223
	dynamic_cast<Tetrahedron*>(const_cast<Element*>(coarsenList.getNeighbourElement(index, dir)));
224

225
      if (!neigh || neigh->isLeaf()) {
226
227
228
229
	/****************************************************************************/
	/*  boundary face or  neigh has been coarsend: free the dof's in the face   */
	/****************************************************************************/

230
	if (mesh->getNumberOfDofs(EDGE)) {
231
	  int node = mesh->getNode(EDGE) + Tetrahedron::nChildEdge[el_type][0][dir];
232
	  mesh->freeDof(const_cast<DegreeOfFreedom*>(child[0]->getDof(node)), EDGE);
233
	}
234
	if (mesh->getNumberOfDofs(FACE)) {
235
	  int node = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][0][dir];
236
	  mesh->freeDof(const_cast<DegreeOfFreedom*>(child[0]->getDof(node)), FACE);
237
	  node = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][1][dir];
238
	  mesh->freeDof(const_cast<DegreeOfFreedom*>(child[1]->getDof(node)), FACE);
239
240
241
242
243
244
245
246
247
	}
      }
    }

    /****************************************************************************/
    /*  finally remove the interior dof's: in the common face of child[0] and   */
    /*  child[1] and at the two barycenter                                      */
    /****************************************************************************/

248
    if (mesh->getNumberOfDofs(FACE)) {
249
      int node = mesh->getNode(FACE);
250
      mesh->freeDof(const_cast<DegreeOfFreedom*>(child[0]->getDof(node)), FACE);
251
252
253
    }


254
    if (mesh->getNumberOfDofs(CENTER)) {
255
256
      int node = mesh->getNode(CENTER);
      for (int i = 0; i < 2; i++)
257
	mesh->freeDof(const_cast<DegreeOfFreedom*>(child[i]->getDof(node)), CENTER);
258
259
260
261
262
263
264
265
    }

    /****************************************************************************/
    /*  get new data on parent and transfer data from children to parent        */
    /****************************************************************************/

    el->coarsenElementData(child[0], child[1], el_type);

266
267
    el->setFirstChild(NULL);
    el->setSecondChild(NULL);
268
269
270
271
272
273
274
275
276
277
278
279

    mesh->freeElement(child[0]);
    mesh->freeElement(child[1]);

    el->incrementMark();
  
    mesh->incrementNumberOfLeaves(-1);
    mesh->incrementNumberOfElements(-2);
  }

  /****************************************************************************/
  /*  get_coarse_patch:  gets the patch for coarsening starting on element    */
280
  /*  elInfo->el in direction of neighbour [3-dir]; returns 1 if a boundary  */
281
282
283
284
285
286
287
288
289
290
291
292
293
  /*  reached and 0 if we come back to the starting element.                  */
  /*                                                                          */
  /*  if NEIGH_IN_EL we only can find the complete coarsening patch if the    */
  /*  can be coarsend; otherwise neighbour information is not valid for       */
  /*  parents; in such situation we stop looping around the edge and return 0 */
  /*                                                                          */
  /*  if !NEIGH_IN_EL we complete the loop also in the case of a incompatible */
  /*  coarsening patch since then all marks of patch elements are reset by    */
  /*  do_coarse_patch() and this minimizes calls of traverse_neighbour();     */
  /*  if we reach a boundary while looping around the edge we loop back to    */
  /*  the starting element before we return                                   */
  /****************************************************************************/
 
294
  bool CoarseningManager3d::getCoarsenPatch(ElInfo *elInfo, 
295
					    DegreeOfFreedom *edge[2],
296
					    int dir, 
297
					    RCNeighbourList &coarsenList, 
298
					    int *n_neigh)
299
  {
300
    FUNCNAME_DBG("CoarseningManager3d::getCoarsenPatch()");
301

302
    static const unsigned char next_el[6][2] = {{3,2},
303
304
305
306
307
308
					  {1,3},
					  {1,2},
					  {0,3},
					  {0,2},
					  {0,1}};

309
310
311
312
    Tetrahedron *el = 
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(elInfo->getElement()));
    Tetrahedron *neigh = 
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(elInfo->getNeighbour(3 - dir)));
313
    if (neigh == NULL)
314
315
      return true;

316
317
    int opp_v = elInfo->getOppVertex(3 - dir);
    ElInfo *neigh_info = stack->traverseNeighbour3d(elInfo, 3 - dir);
318

319
    TEST_EXIT_DBG(neigh == neigh_info->getElement())
320
321
322
323
324
325
      ("neigh %d and neigh_info->el %d are not identical\n", 
       neigh->getIndex(), neigh_info->getElement()->getIndex());
    /****************************************************************************/
    /*  we have to go back to the starting element via opp_v values             */
    /*  correct information is produce by get_neigh_on_patch()                  */
    /****************************************************************************/
326
327
328
    coarsenList.setOppVertex(*n_neigh, 0, opp_v);  
    coarsenList.setElement(*n_neigh, neigh);
    coarsenList.setElType(*n_neigh, neigh_info->getType()); 
329
  
330
331
    int nVertices = mesh->getGeo(VERTEX);			 

332
    while (neigh != el) {
333
334
335
336
337
338
339
      // === Find the coarsening edge on the current neighbour element. ===

      int i, j, k;

      // === First, try to identify the edge DOFs directly. ===

      for (j = 0; j < nVertices; j++)
340
341
	if (neigh->getDof(j) == edge[0])
	  break;
342
343

      for (k = 0; k < nVertices; k++)
344
345
	if (neigh->getDof(k) == edge[1])
	  break;
346

347
348
349
350
351
352
353
      // === If one of the edge DOFs was not found, try to make use of periodic ===
      // === DOF associations. First, use the direct associations between DOFs. ===
      // === If this will not work, continue with testing on indirect           ===
      // === associations.                                                      ===

      if (j >= nVertices) {
	for (j = 0; j < nVertices; j++)
354
355
	  if (mesh->associated(neigh->getDof(j, 0), edge[0][0]))
	    break;
356
357
358
359
360
361
362
363
364
365
366
367
368
369
	
	if (j >= nVertices)
	  for (j = 0; j < nVertices; j++)
	    if (mesh->indirectlyAssociated(neigh->getDof(j, 0), edge[0][0]))
	      break;

	TEST_EXIT_DBG(j < nVertices)
	  ("Process element %d: DOF %d not found on element %d with nodes (%d %d %d %d)\n", 
	   el->getIndex(), edge[0][0], neigh->getIndex(), neigh->getDof(0, 0),
	   neigh->getDof(1, 0), neigh->getDof(2, 0), neigh->getDof(3, 0));
      }
      
      if (k >= nVertices) {
	for (k = 0; k < nVertices; k++)
370
371
	  if (mesh->associated(neigh->getDof(k, 0), edge[1][0]))
	    break;
372
373
374
375
376

	if (k >= nVertices)
	  for (k = 0; k < nVertices; k++)
	    if (mesh->indirectlyAssociated(neigh->getDof(k, 0), edge[1][0]))
	      break;
377
     
378
379
380
	TEST_EXIT_DBG(k < nVertices)
	  ("Process element %d: DOF %d not found on element %d with nodes (%d %d %d %d)\n", 
	   el->getIndex(), edge[1][0], neigh->getIndex(), neigh->getDof(0, 0),
381
	   neigh->getDof(1, 0), neigh->getDof(2, 0), neigh->getDof(3, 0));
382
383
      }

384
385
386
387
388
      int edgeNo = Tetrahedron::edgeOfDofs[j][k];
      coarsenList.setCoarsePatch(*n_neigh, edgeNo == 0);


      // === Get the direction of the next neighbour. ===
389

390
391
      if (next_el[edgeNo][0] != opp_v)
	i = next_el[edgeNo][0];
392
      else
393
	i = next_el[edgeNo][1];
394
395
396
397

      ++*n_neigh;

      opp_v = neigh_info->getOppVertex(i);
398
399
      neigh = 
	dynamic_cast<Tetrahedron*>(const_cast<Element*>(neigh_info->getNeighbour(i)));
400
      if (neigh) {
401
	neigh_info = stack->traverseNeighbour3d(neigh_info, i);
402
	TEST_EXIT_DBG(neigh == neigh_info->getElement())
403
404
405
406
407
408
	  ("neigh %d and neigh_info->el %d are not identical\n", 
	   neigh->getIndex(), neigh_info->getElement()->getIndex());
	/****************************************************************************/
	/*  we have to go back to the starting element via opp_v values             */
	/*  correct information is produce by get_neigh_on_patch()                  */
	/****************************************************************************/
409
410
411
	coarsenList.setOppVertex(*n_neigh, 0, opp_v);  
	coarsenList.setElement(*n_neigh, neigh);
	coarsenList.setElType(*n_neigh, neigh_info->getType());
412
      } else {
413
	break;
414
      }
415
416
417
418
419
420
421
422
423
    }

    if (neigh == el)  
      return false;

    /****************************************************************************/
    /*  the domain's boundary is reached; loop back to the starting el          */
    /****************************************************************************/

424
    int i = *n_neigh - 1;
425
    opp_v = coarsenList.getOppVertex(i, 0);
426
    do {
427
      TEST_EXIT_DBG(neigh_info->getNeighbour(opp_v)  &&  i > 0)
428
	("while looping back domains boundary was reached or i == 0\n");
429
      opp_v = coarsenList.getOppVertex(i--, 0);
430
      neigh_info = stack->traverseNeighbour3d(neigh_info, opp_v);
431
    } while (neigh_info->getElement() != el);
432
433
434
435
436
437
438
439
440

    return true;
  }

  /****************************************************************************/
  /*  coarse_patch: first rebuild the dofs on the parents then do restriction */
  /*  of data (if possible) and finally coarsen the patch elements            */
  /****************************************************************************/

441
  void CoarseningManager3d::coarsenPatch(RCNeighbourList &coarsenList, 
442
443
					 int n_neigh, 
					 int bound)
444
  {
445
    FUNCNAME_DBG("CoarseningManager3d::coarsenPatch()");
446
447

    Tetrahedron *el = 
448
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(coarsenList.getElement(0)));
449
    DegreeOfFreedom *dof = NULL;
450
451
452

    TEST_EXIT_DBG(el)("No element!\n");
    TEST_EXIT_DBG(el->getChild(0))("No child in element!\n");
453

454
    if (mesh->getNumberOfDofs(EDGE)) {
455
456
457
      // === Get dof for coarsening edge. ===

      int node = mesh->getNode(EDGE);
458
      if (!(dof = const_cast<DegreeOfFreedom*>(el->getDof(node))))
459
	dof = mesh->getDof(EDGE);
460
    }
461

462
463
464
    if (mesh->getNumberOfDofs(EDGE) || 
	mesh->getNumberOfDofs(FACE) || 
	mesh->getNumberOfDofs(CENTER)) {     
465
      for (int i = 0; i < n_neigh; i++)
466
	coarsenList.addDOFParent(i, dof);
467
468
    }
      
469

470
    // === Restrict dof vectors to the parents on the patch. ===
471
472

    int nrAdmin = mesh->getNumberOfDOFAdmin();
473
    for (int iadmin = 0; iadmin < nrAdmin; iadmin++) {
474
      std::list<DOFIndexedBase*>::iterator it;
475
      DOFAdmin* admin = const_cast<DOFAdmin*>(&mesh->getDofAdmin(iadmin));
476
      std::list<DOFIndexedBase*>::iterator end = admin->endDOFIndexed();
477
      for (it = admin->beginDOFIndexed(); it != end; ++it)
478
	(*it)->coarseRestrict(coarsenList, n_neigh);
479
480
    } 

481
    // === And now start to coarsen the patch: remove dof's of the coarsening edge. ===
482

483
    mesh->freeDof(const_cast<DegreeOfFreedom*>(el->getChild(0)->getDof(3)), VERTEX);
484
485
    mesh->incrementNumberOfVertices(-1);

486
    if (mesh->getNumberOfDofs(EDGE)) {
487
      int node = mesh->getNode(EDGE) + 2;
488
489
      mesh->freeDof(const_cast<DegreeOfFreedom*>(el->getChild(0)->getDof(node)), EDGE);
      mesh->freeDof(const_cast<DegreeOfFreedom*>(el->getChild(1)->getDof(node)), EDGE);
490
491
    }

492
493
    if (coarsenList.getElement(0)->isNewCoordSet())
      coarsenList.getElement(0)->eraseNewCoord();   
494

495
    for (int i = 0; i < n_neigh; i++) {
496
      coarsenList.getElement(i)->setNewCoord(NULL);
497
498
      coarsenTetrahedron(coarsenList, i);
    }
499

500
501
502
    // === If the coarsening edge is an interior edge there are  n_neigh + 1 edges ===
    // === and 2 * n_neigh + 1 faces removed; if it is a boundary edge it is one   ===
    // === more edge and one more face.                                            ===
503
504
505

    if (bound) {
      mesh->incrementNumberOfEdges(-(n_neigh + 2));
506
      mesh->incrementNumberOfFaces(-(2 * n_neigh + 1));
507
508
    } else {
      mesh->incrementNumberOfEdges(-(n_neigh + 1));
509
      mesh->incrementNumberOfFaces(-(2 * n_neigh));
510
511
512
513
    }
  }

}