CoarseningManager3d.cc 18.4 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

namespace AMDiS {

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

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

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

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

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

55
      el->setMark(0);
56
      return;
57
58
59
    }

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

65
66
67
    DegreeOfFreedom *edge[2];
    int n_neigh, bound = 0;

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

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

78
79
80
    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));
81
    } else {
82
83
      edge[1] = const_cast<DegreeOfFreedom*>(el->getDof(0));
      edge[0] = const_cast<DegreeOfFreedom*>(el->getDof(1));
84
85
    }

86
    coarsenList.setElement(0, el, true);
87
88
    n_neigh = 1;

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


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

102
103
104
    // === 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
105
    if (refineEdges.size()) {
106
107
108
      // TODO: Remove these two lines and make something more meaningful!!
      el->setMark(0);
      return;
Thomas Witkowski's avatar
Thomas Witkowski committed
109
110

      Element *otherEl = refineEdges[0].first;
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
      
      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 = 
132
		!((i > 0 && (edgeNo0 >= 0 && !el2->getChild(0)->isLeaf())) || 
133
134
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
		  (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
    
167
    coarsenList.fillNeighbourRelations(n_neigh, bound);
168
169
170
171
172
173
174
175
176

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

      ++*n_neigh;

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

    if (neigh == el)  
      return false;

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

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

    return true;
  }

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

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

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

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

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

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

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

467
    // === Restrict dof vectors to the parents on the patch. ===
468
469

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

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

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

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

489
490
    if (coarsenList.getElement(0)->isNewCoordSet())
      coarsenList.getElement(0)->eraseNewCoord();   
491

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

497
498
499
    // === 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.                                            ===
500
501
502

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

}