CoarseningManager3d.cc 17.1 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
    FUNCNAME("CoarseningManager3d::coarsenFunction()");

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

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

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

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

44
      el->setMark(0);
45
      return;
46
47
48
    }

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

54
55
56
    DegreeOfFreedom *edge[2];
    int n_neigh, bound = 0;

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

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

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

75
    coarsenList.setElement(0, el, true);
76
77
    n_neigh = 1;

78
    coarsenList.setOppVertex(0, 0, 0);
79
    coarsenList.setElType(0, elInfo->getType());
80
    bound = false;
81
82
    if (getCoarsenPatch(elInfo, edge, 0, coarsenList, &n_neigh)) {
      getCoarsenPatch(elInfo, edge, 1, coarsenList, &n_neigh);
83
84
      bound = true;
    }
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
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


#if HAVE_PARALLEL_DOMAIN_AMDIS
    Element *otherEl = NULL;
    int otherEdge = -1;      
    FixRefinementPatch::getOtherEl(stack, &otherEl, otherEdge);
    
    // === If the refinement edge must be fixed, add also the other part of this ===
    // === edge to the refinement patch.                                         ===
    
    if (otherEl) {
      // TODO: Remove these two lines and make something more meaningful!!
      el->setMark(0);
      return;
      
      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
    
156
    coarsenList.fillNeighbourRelations(n_neigh, bound);
157
158
159
160
161
162
163
164
165

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

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

166
    if (coarsenList.doCoarsePatch(n_neigh)) {
167
168
      int n_neigh_periodic;
      DegreeOfFreedom *next_edge[2];
169
      RCNeighbourList periodicList;
170

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

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

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

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

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

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

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

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

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

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


240
    if (mesh->getNumberOfDofs(CENTER)) {
241
242
      int node = mesh->getNode(CENTER);
      for (int i = 0; i < 2; i++)
243
	mesh->freeDof(const_cast<int*>(child[i]->getDof(node)), CENTER);
244
245
246
247
248
249
250
251
252
253
254
255
256
257
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);

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

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

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

302
303
    int opp_v = elInfo->getOppVertex(3 - dir);
    ElInfo *neigh_info = stack->traverseNeighbour3d(elInfo, 3 - dir);
304

305
    TEST_EXIT_DBG(neigh == neigh_info->getElement())
306
307
308
309
310
311
      ("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()                  */
    /****************************************************************************/
312
313
314
    coarsenList.setOppVertex(*n_neigh, 0, opp_v);  
    coarsenList.setElement(*n_neigh, neigh);
    coarsenList.setElType(*n_neigh, neigh_info->getType()); 
315
316
  
    int n_vertices = mesh->getGeo(VERTEX);
317
    int i, j, k, edge_no;
318
319
320
			 
    while (neigh != el) {
      for (j = 0; j < n_vertices; j++)
321
322
	if (neigh->getDof(j) == edge[0])
	  break;
323
      for (k = 0; k < n_vertices; k++)
324
325
	if (neigh->getDof(k) == edge[1])
	  break;
326

327
      if (j > 3 || k > 3) {
328
	for (j = 0; j < n_vertices; j++)
329
330
	  if (mesh->associated(neigh->getDof(j, 0), edge[0][0]))
	    break;
331
	for (k = 0; k < n_vertices; k++)
332
333
	  if (mesh->associated(neigh->getDof(k, 0), edge[1][0]))
	    break;
334
     
335
336
337
338
339
340
341
342
343
	TEST_EXIT_DBG(j < n_vertices)
	  ("dof %d not found on element %d with nodes (%d %d %d %d)\n", 
	   edge[0][0], neigh->getIndex(), neigh->getDof(0, 0),
	   neigh->getDof(1, 0), neigh->getDof(2, 0), neigh->getDof(3, 0));
	
	TEST_EXIT_DBG(k < n_vertices)
	  ("dof %d not found on element %d with nodes (%d %d %d %d)\n", 
	   edge[1][0], neigh->getIndex(), neigh->getDof(0, 0),
	   neigh->getDof(1, 0), neigh->getDof(2, 0), neigh->getDof(3, 0));
344
      }
345
      edge_no = Tetrahedron::edgeOfDofs[j][k];
346
      coarsenList.setCoarsePatch(*n_neigh, edge_no == 0);
347
348
349
350
351
352
353
354
355
356
357
358
359

      /****************************************************************************/
      /*  get the direction of the next neighbour				    */
      /****************************************************************************/

      if (next_el[edge_no][0] != opp_v)
	i = next_el[edge_no][0];
      else
	i = next_el[edge_no][1];

      ++*n_neigh;

      opp_v = neigh_info->getOppVertex(i);
360
361
      neigh = 
	dynamic_cast<Tetrahedron*>(const_cast<Element*>(neigh_info->getNeighbour(i)));
362
      if (neigh) {
363
	neigh_info = stack->traverseNeighbour3d(neigh_info, i);
364
	TEST_EXIT_DBG(neigh == neigh_info->getElement())
365
366
367
368
369
370
	  ("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()                  */
	/****************************************************************************/
371
372
373
	coarsenList.setOppVertex(*n_neigh, 0, opp_v);  
	coarsenList.setElement(*n_neigh, neigh);
	coarsenList.setElType(*n_neigh, neigh_info->getType());
374
      } else {
375
	break;
376
      }
377
378
379
380
381
382
383
384
385
    }

    if (neigh == el)  
      return false;

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

386
    i = *n_neigh - 1;
387
    opp_v = coarsenList.getOppVertex(i, 0);
388
    do {
389
      TEST_EXIT_DBG(neigh_info->getNeighbour(opp_v)  &&  i > 0)
390
	("while looping back domains boundary was reached or i == 0\n");
391
      opp_v = coarsenList.getOppVertex(i--, 0);
392
      neigh_info = stack->traverseNeighbour3d(neigh_info, opp_v);
393
    } while (neigh_info->getElement() != el);
394
395
396
397
398
399
400
401
402

    return true;
  }

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

403
  void CoarseningManager3d::coarsenPatch(RCNeighbourList &coarsenList, 
404
405
					 int n_neigh, 
					 int bound)
406
  {
407
408
409
    FUNCNAME("CoarseningManager3d::coarsenPatch()");

    Tetrahedron *el = 
410
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(coarsenList.getElement(0)));
411
412
413
414
    DegreeOfFreedom *dof = NULL;

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

416
    if (mesh->getNumberOfDofs(EDGE)) {
417
418
419
420
      // === Get dof for coarsening edge. ===

      int node = mesh->getNode(EDGE);
      if (!(dof = const_cast<int*>(el->getDof(node))))
421
	dof = mesh->getDof(EDGE);
422
    }
423

424
425
426
    if (mesh->getNumberOfDofs(EDGE) || 
	mesh->getNumberOfDofs(FACE) || 
	mesh->getNumberOfDofs(CENTER)) {     
427
      for (int i = 0; i < n_neigh; i++)
428
	coarsenList.addDOFParent(i, dof);
429
430
    }
      
431

432
    // === Restrict dof vectors to the parents on the patch. ===
433
434

    int nrAdmin = mesh->getNumberOfDOFAdmin();
435
    for (int iadmin = 0; iadmin < nrAdmin; iadmin++) {
436
      std::list<DOFIndexedBase*>::iterator it;
437
      DOFAdmin* admin = const_cast<DOFAdmin*>(&mesh->getDofAdmin(iadmin));
438
      std::list<DOFIndexedBase*>::iterator end = admin->endDOFIndexed();
439
      for (it = admin->beginDOFIndexed(); it != end; ++it)
440
	(*it)->coarseRestrict(coarsenList, n_neigh);
441
442
    } 

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

445
    mesh->freeDof(const_cast<int*>(el->getChild(0)->getDof(3)), VERTEX);
446
447
    mesh->incrementNumberOfVertices(-1);

448
    if (mesh->getNumberOfDofs(EDGE)) {
449
450
451
      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);
452
453
    }

454
455
    if (coarsenList.getElement(0)->isNewCoordSet())
      coarsenList.getElement(0)->eraseNewCoord();   
456

457
    for (int i = 0; i < n_neigh; i++) {
458
      coarsenList.getElement(i)->setNewCoord(NULL);
459
460
      coarsenTetrahedron(coarsenList, i);
    }
461

462
463
464
    // === 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.                                            ===
465
466
467

    if (bound) {
      mesh->incrementNumberOfEdges(-(n_neigh + 2));
468
      mesh->incrementNumberOfFaces(-(2 * n_neigh + 1));
469
470
    } else {
      mesh->incrementNumberOfEdges(-(n_neigh + 1));
471
      mesh->incrementNumberOfFaces(-(2 * n_neigh));
472
473
474
475
    }
  }

}