CoarseningManager3d.cc 17.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
15
16
17
18
19
20
21
#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"
22
#include "Debug.h"
23
24
25

namespace AMDiS {

26
  void CoarseningManager3d::coarsenFunction(ElInfo *elInfo)
27
  {
28
29
30
#if HAVE_PARALLEL_DOMAIN_AMDIS
    FUNCNAME_DBG("CoarseningManager3d::coarsenFunction()");
#endif
31

32
    Tetrahedron *el = 
33
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(elInfo->getElement()));
34

35
    // If element must not be coarsend, return.
36
    if (el->getMark() >= 0) 
37
      return; 
38

39
    // Single leaves don't get coarsened.
40
    if (el->isLeaf()) 
41
      return;
42

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

46
      el->setMark(0);
47
      return;
48
49
50
    }

    if (!(el->getChild(0)->isLeaf()) || !(el->getChild(1)->isLeaf())) {
51
      // One of the children is not a leaf element. Try again later on.
52
      doMore = true;
53
      return;
54
    }
55

56
57
58
    DegreeOfFreedom *edge[2];
    int n_neigh, bound = 0;

59
60
61
    /****************************************************************************/
    /*  get a list for storing all elements at the coarsening edge and fill it  */
    /****************************************************************************/
62
63
    RCNeighbourList coarsenList(mesh->getMaxEdgeNeigh());
    coarsenList.setCoarseningManager(this);
64
65
66
67
68

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

69
70
71
    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));
72
    } else {
73
74
      edge[1] = const_cast<DegreeOfFreedom*>(el->getDof(0));
      edge[0] = const_cast<DegreeOfFreedom*>(el->getDof(1));
75
76
    }

77
    coarsenList.setElement(0, el, true);
78
79
    n_neigh = 1;

80
    coarsenList.setOppVertex(0, 0, 0);
81
    coarsenList.setElType(0, elInfo->getType());
82
    bound = false;
83
84
    if (getCoarsenPatch(elInfo, edge, 0, coarsenList, &n_neigh)) {
      getCoarsenPatch(elInfo, edge, 1, coarsenList, &n_neigh);
85
86
      bound = true;
    }
87
88
89


#if HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
90
91
92
    vector<FixRefinementPatch::EdgeInEl> refineEdges;
    FixRefinementPatch::getOtherEl(stack, refineEdges);

93
94
95
    // === 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
96
    if (refineEdges.size()) {
97
98
99
      // TODO: Remove these two lines and make something more meaningful!!
      el->setMark(0);
      return;
Thomas Witkowski's avatar
Thomas Witkowski committed
100
101

      Element *otherEl = refineEdges[0].first;
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
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
      
      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 = 
		!(i > 0 &&
		  (edgeNo0 >= 0 && !el2->getChild(0)->isLeaf()) || 
		  (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
    
159
    coarsenList.fillNeighbourRelations(n_neigh, bound);
160
161
162
163
164
165
166
167
168

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

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

169
    if (coarsenList.doCoarsePatch(n_neigh)) {
170
171
      int n_neigh_periodic;
      DegreeOfFreedom *next_edge[2];
172
      RCNeighbourList periodicList;
173

174
      while (edge[0] != NULL) {
175
	coarsenList.periodicSplit(edge, next_edge,
176
177
				  &n_neigh, &n_neigh_periodic,
				  periodicList);
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193

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

194
  void CoarseningManager3d::coarsenTetrahedron(RCNeighbourList &coarsenList, 
195
					       int index)
196
  {
197
    Tetrahedron *el = 
198
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(coarsenList.getElement(index)));
199
    Tetrahedron *child[2];
200

201
202
    child[0] = dynamic_cast<Tetrahedron*>(const_cast<Element*>(el->getChild(0)));
    child[1] = dynamic_cast<Tetrahedron*>(const_cast<Element*>(el->getChild(1)));
203
    int el_type = coarsenList.getType(index);
204
205
206
207
208
209

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

210
211
    for (int dir = 0; dir < 2; dir++) {
      Tetrahedron *neigh = 
212
	dynamic_cast<Tetrahedron*>(const_cast<Element*>(coarsenList.getNeighbourElement(index, dir)));
213

214
      if (!neigh || neigh->isLeaf()) {
215
216
217
218
	/****************************************************************************/
	/*  boundary face or  neigh has been coarsend: free the dof's in the face   */
	/****************************************************************************/

219
	if (mesh->getNumberOfDofs(EDGE)) {
220
	  int node = mesh->getNode(EDGE) + Tetrahedron::nChildEdge[el_type][0][dir];
221
	  mesh->freeDof(const_cast<int*>(child[0]->getDof(node)), EDGE);
222
	}
223
	if (mesh->getNumberOfDofs(FACE)) {
224
	  int node = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][0][dir];
225
	  mesh->freeDof(const_cast<int*>(child[0]->getDof(node)), FACE);
226
	  node = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][1][dir];
227
	  mesh->freeDof(const_cast<int*>(child[1]->getDof(node)), FACE);
228
229
230
231
232
233
234
235
236
	}
      }
    }

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

237
    if (mesh->getNumberOfDofs(FACE)) {
238
      int node = mesh->getNode(FACE);
239
      mesh->freeDof(const_cast<int*>(child[0]->getDof(node)), FACE);
240
241
242
    }


243
    if (mesh->getNumberOfDofs(CENTER)) {
244
245
      int node = mesh->getNode(CENTER);
      for (int i = 0; i < 2; i++)
246
	mesh->freeDof(const_cast<int*>(child[i]->getDof(node)), CENTER);
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
    }

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

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

    el->setFirstChild(NULL);
    el->setSecondChild(NULL);

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

    static unsigned char next_el[6][2] = {{3,2},
					  {1,3},
					  {1,2},
					  {0,3},
					  {0,2},
					  {0,1}};

298
299
300
301
    Tetrahedron *el = 
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(elInfo->getElement()));
    Tetrahedron *neigh = 
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(elInfo->getNeighbour(3 - dir)));
302
    if (neigh == NULL)
303
304
      return true;

305
306
    int opp_v = elInfo->getOppVertex(3 - dir);
    ElInfo *neigh_info = stack->traverseNeighbour3d(elInfo, 3 - dir);
307

308
    TEST_EXIT_DBG(neigh == neigh_info->getElement())
309
310
311
312
313
314
      ("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()                  */
    /****************************************************************************/
315
316
317
    coarsenList.setOppVertex(*n_neigh, 0, opp_v);  
    coarsenList.setElement(*n_neigh, neigh);
    coarsenList.setElType(*n_neigh, neigh_info->getType()); 
318
  
319
320
    int nVertices = mesh->getGeo(VERTEX);			 

321
    while (neigh != el) {
322
323
324
325
326
327
328
      // === 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++)
329
330
	if (neigh->getDof(j) == edge[0])
	  break;
331
332

      for (k = 0; k < nVertices; k++)
333
334
	if (neigh->getDof(k) == edge[1])
	  break;
335

336
337
338
339
340
341
342
      // === 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++)
343
344
	  if (mesh->associated(neigh->getDof(j, 0), edge[0][0]))
	    break;
345
346
347
348
349
350
351
352
353
354
355
356
357
358
	
	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++)
359
360
	  if (mesh->associated(neigh->getDof(k, 0), edge[1][0]))
	    break;
361
362
363
364
365

	if (k >= nVertices)
	  for (k = 0; k < nVertices; k++)
	    if (mesh->indirectlyAssociated(neigh->getDof(k, 0), edge[1][0]))
	      break;
366
     
367
368
369
	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),
370
	   neigh->getDof(1, 0), neigh->getDof(2, 0), neigh->getDof(3, 0));
371
372
      }

373
374
375
376
377
      int edgeNo = Tetrahedron::edgeOfDofs[j][k];
      coarsenList.setCoarsePatch(*n_neigh, edgeNo == 0);


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

379
380
      if (next_el[edgeNo][0] != opp_v)
	i = next_el[edgeNo][0];
381
      else
382
	i = next_el[edgeNo][1];
383
384
385
386

      ++*n_neigh;

      opp_v = neigh_info->getOppVertex(i);
387
388
      neigh = 
	dynamic_cast<Tetrahedron*>(const_cast<Element*>(neigh_info->getNeighbour(i)));
389
      if (neigh) {
390
	neigh_info = stack->traverseNeighbour3d(neigh_info, i);
391
	TEST_EXIT_DBG(neigh == neigh_info->getElement())
392
393
394
395
396
397
	  ("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()                  */
	/****************************************************************************/
398
399
400
	coarsenList.setOppVertex(*n_neigh, 0, opp_v);  
	coarsenList.setElement(*n_neigh, neigh);
	coarsenList.setElType(*n_neigh, neigh_info->getType());
401
      } else {
402
	break;
403
      }
404
405
406
407
408
409
410
411
412
    }

    if (neigh == el)  
      return false;

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

413
    int i = *n_neigh - 1;
414
    opp_v = coarsenList.getOppVertex(i, 0);
415
    do {
416
      TEST_EXIT_DBG(neigh_info->getNeighbour(opp_v)  &&  i > 0)
417
	("while looping back domains boundary was reached or i == 0\n");
418
      opp_v = coarsenList.getOppVertex(i--, 0);
419
      neigh_info = stack->traverseNeighbour3d(neigh_info, opp_v);
420
    } while (neigh_info->getElement() != el);
421
422
423
424
425
426
427
428
429

    return true;
  }

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

430
  void CoarseningManager3d::coarsenPatch(RCNeighbourList &coarsenList, 
431
432
					 int n_neigh, 
					 int bound)
433
  {
434
    FUNCNAME_DBG("CoarseningManager3d::coarsenPatch()");
435
436

    Tetrahedron *el = 
437
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(coarsenList.getElement(0)));
438
439
440
441
    DegreeOfFreedom *dof = NULL;

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

443
    if (mesh->getNumberOfDofs(EDGE)) {
444
445
446
447
      // === Get dof for coarsening edge. ===

      int node = mesh->getNode(EDGE);
      if (!(dof = const_cast<int*>(el->getDof(node))))
448
	dof = mesh->getDof(EDGE);
449
    }
450

451
452
453
    if (mesh->getNumberOfDofs(EDGE) || 
	mesh->getNumberOfDofs(FACE) || 
	mesh->getNumberOfDofs(CENTER)) {     
454
      for (int i = 0; i < n_neigh; i++)
455
	coarsenList.addDOFParent(i, dof);
456
457
    }
      
458

459
    // === Restrict dof vectors to the parents on the patch. ===
460
461

    int nrAdmin = mesh->getNumberOfDOFAdmin();
462
    for (int iadmin = 0; iadmin < nrAdmin; iadmin++) {
463
      std::list<DOFIndexedBase*>::iterator it;
464
      DOFAdmin* admin = const_cast<DOFAdmin*>(&mesh->getDofAdmin(iadmin));
465
      std::list<DOFIndexedBase*>::iterator end = admin->endDOFIndexed();
466
      for (it = admin->beginDOFIndexed(); it != end; ++it)
467
	(*it)->coarseRestrict(coarsenList, n_neigh);
468
469
    } 

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

472
    mesh->freeDof(const_cast<int*>(el->getChild(0)->getDof(3)), VERTEX);
473
474
    mesh->incrementNumberOfVertices(-1);

475
    if (mesh->getNumberOfDofs(EDGE)) {
476
477
478
      int node = mesh->getNode(EDGE) + 2;
      mesh->freeDof(const_cast<int*>(el->getChild(0)->getDof(node)), EDGE);
      mesh->freeDof(const_cast<int*>(el->getChild(1)->getDof(node)), EDGE);
479
480
    }

481
482
    if (coarsenList.getElement(0)->isNewCoordSet())
      coarsenList.getElement(0)->eraseNewCoord();   
483

484
    for (int i = 0; i < n_neigh; i++) {
485
      coarsenList.getElement(i)->setNewCoord(NULL);
486
487
      coarsenTetrahedron(coarsenList, i);
    }
488

489
490
491
    // === 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.                                            ===
492
493
494

    if (bound) {
      mesh->incrementNumberOfEdges(-(n_neigh + 2));
495
      mesh->incrementNumberOfFaces(-(2 * n_neigh + 1));
496
497
    } else {
      mesh->incrementNumberOfEdges(-(n_neigh + 1));
498
      mesh->incrementNumberOfFaces(-(2 * n_neigh));
499
500
501
502
    }
  }

}