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
102
103
    vector<FixRefinementPatch::EdgeInEl> refineEdges;
    FixRefinementPatch::getOtherEl(stack, refineEdges);

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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
      
      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 = 
134
		!((i > 0 && (edgeNo0 >= 0 && !el2->getChild(0)->isLeaf())) || 
135
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
		  (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
    
169
    coarsenList.fillNeighbourRelations(n_neigh, bound);
170
171
172
173
174
175
176
177
178

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

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

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

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

	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                                          */
  /*****************************************************************************/

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

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

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

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

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

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

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

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


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

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

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

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

    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    */
279
  /*  elInfo->el in direction of neighbour [3-dir]; returns 1 if a boundary  */
280
281
282
283
284
285
286
287
288
289
290
291
292
  /*  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                                   */
  /****************************************************************************/
 
293
  bool CoarseningManager3d::getCoarsenPatch(ElInfo *elInfo, 
294
					    DegreeOfFreedom *edge[2],
295
					    int dir, 
296
					    RCNeighbourList &coarsenList, 
297
					    int *n_neigh)
298
  {
299
    FUNCNAME_DBG("CoarseningManager3d::getCoarsenPatch()");
300

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

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

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

318
    TEST_EXIT_DBG(neigh == neigh_info->getElement())
319
320
321
322
323
324
      ("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()                  */
    /****************************************************************************/
325
326
327
    coarsenList.setOppVertex(*n_neigh, 0, opp_v);  
    coarsenList.setElement(*n_neigh, neigh);
    coarsenList.setElType(*n_neigh, neigh_info->getType()); 
328
  
329
330
    int nVertices = mesh->getGeo(VERTEX);			 

331
    while (neigh != el) {
332
333
334
335
336
337
338
      // === 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++)
339
340
	if (neigh->getDof(j) == edge[0])
	  break;
341
342

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

346
347
348
349
350
351
352
      // === 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++)
353
354
	  if (mesh->associated(neigh->getDof(j, 0), edge[0][0]))
	    break;
355
356
357
358
359
360
361
362
363
364
365
366
367
368
	
	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++)
369
370
	  if (mesh->associated(neigh->getDof(k, 0), edge[1][0]))
	    break;
371
372
373
374
375

	if (k >= nVertices)
	  for (k = 0; k < nVertices; k++)
	    if (mesh->indirectlyAssociated(neigh->getDof(k, 0), edge[1][0]))
	      break;
376
     
377
378
379
	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),
380
	   neigh->getDof(1, 0), neigh->getDof(2, 0), neigh->getDof(3, 0));
381
382
      }

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


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

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

      ++*n_neigh;

      opp_v = neigh_info->getOppVertex(i);
397
398
      neigh = 
	dynamic_cast<Tetrahedron*>(const_cast<Element*>(neigh_info->getNeighbour(i)));
399
      if (neigh) {
400
	neigh_info = stack->traverseNeighbour3d(neigh_info, i);
401
	TEST_EXIT_DBG(neigh == neigh_info->getElement())
402
403
404
405
406
407
	  ("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()                  */
	/****************************************************************************/
408
409
410
	coarsenList.setOppVertex(*n_neigh, 0, opp_v);  
	coarsenList.setElement(*n_neigh, neigh);
	coarsenList.setElType(*n_neigh, neigh_info->getType());
411
      } else {
412
	break;
413
      }
414
415
416
417
418
419
420
421
422
    }

    if (neigh == el)  
      return false;

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

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

    return true;
  }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

499
500
501
    // === 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.                                            ===
502
503
504

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

}