Am Montag, 13. Mai 2022, finden Wartungsarbeiten am Gitlab-Server (Update auf neue Version statt). Der Dienst wird daher am Montag für einige Zeit nicht verfügbar sein.
On Monday, May 13th 2022, the Gitlab server will be updated. The service will therefore not be accessible for some time on Monday.

RefinementManager3d.cc 24.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#include "RefinementManager.h"
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "DOFAdmin.h"
#include "AdaptStationary.h"
#include "AdaptInstationary.h"
#include "FixVec.h"
#include "RCNeighbourList.h"
#include "ProblemStatBase.h"
#include "DOFIndexed.h"
#include "Projection.h"
#include "DOFVector.h"
#include "PeriodicBC.h"
#include "VertexVector.h"

namespace AMDiS {

  void RefinementManager3d::bisectTetrahedron(RCNeighbourList* ref_list, 
					      int index,
					      DegreeOfFreedom* dof[3],
					      DegreeOfFreedom *edge[2])
  {
24
25
26
    Tetrahedron *el = 
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(ref_list->getElement(index)));
    Tetrahedron *child[2];
27
    int el_type = ref_list->getType(index);
28
29
30
31

    child[0] = dynamic_cast<Tetrahedron*>(mesh->createNewElement(el));
    child[1] = dynamic_cast<Tetrahedron*>(mesh->createNewElement(el));
  
32
    int mark = max(0, el->getMark() - 1);
33
34
35
36
37
38
39
40
41
42
43
44
45
    child[0]->setMark(mark);
    child[1]->setMark(mark);
    el->setMark(0);

    /****************************************************************************/
    /*  transfer hidden data from parent to children                            */
    /****************************************************************************/

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

    el->setFirstChild(child[0]);
    el->setSecondChild(child[1]);

46
47
    if (child[0]->getMark() > 0)
      doMoreRecursiveRefine = true;
48
49

    int n_vertices = mesh->getGeo(VERTEX);
50
51
    child[0]->setDof(n_vertices - 1, dof[0]);
    child[1]->setDof(n_vertices - 1, dof[0]);
52

53
    for (int i = 0; i < n_vertices - 1; i++) {
54
      child[0]->
55
	setDof(i, const_cast<int*>(el->getDof(Tetrahedron::childVertex[el_type][0][i])));
56
      child[1]->
57
	setDof(i, const_cast<int*>(el->getDof(Tetrahedron::childVertex[el_type][1][i])));
58
    }
59
60
61
62
63
64
65
66
67
68
69
70
    /****************************************************************************/
    /*  there is one more leaf element and two more hierachical elements        */
    /****************************************************************************/

    mesh->incrementNumberOfLeaves(1);
    mesh->incrementNumberOfElements(2);

    /****************************************************************************/
    /* first set those dof pointers for higher order without neighbour          */
    /* information                                                              */
    /****************************************************************************/

71
    if (mesh->getNumberOfDofs(EDGE)) {
72
      int node = mesh->getNode(EDGE);
73

74
75
76
77
78
      /****************************************************************************/
      /*  set pointers to those dof's that are handed on from the parant          */
      /****************************************************************************/
      
      child[0]->
79
80
	setDof(node, 
	       const_cast<int*>(el->getDof(node + Tetrahedron::childEdge[el_type][0][0])));
81
      child[1]->
82
83
	setDof(node, 
	       const_cast<int*>(el->getDof(node + Tetrahedron::childEdge[el_type][1][0])));
84
      child[0]->
85
86
	setDof(node + 1, 
	       const_cast<int*>(el->getDof(node + Tetrahedron::childEdge[el_type][0][1])));
87
      child[1]->
88
89
	setDof(node + 1, 
	       const_cast<int*>(el->getDof(node + Tetrahedron::childEdge[el_type][1][1])));
90
      child[0]->
91
92
	setDof(node + 3, 
	       const_cast<int*>(el->getDof(node + Tetrahedron::childEdge[el_type][0][3])));
93
      child[1]->
94
95
	setDof(node + 3, 
	       const_cast<int*>(el->getDof(node + Tetrahedron::childEdge[el_type][1][3])));
96
97
98
99
100
      
      /****************************************************************************/
      /*  adjust pointers to the dof's in the refinement edge                     */
      /****************************************************************************/
      
101
102
103
      if (el->getDof(0) == edge[0]) {
	child[0]->setDof(node + 2, dof[1]);
	child[1]->setDof(node + 2, dof[2]);
104
      } else {
105
106
	child[0]->setDof(node + 2, dof[2]);
	child[1]->setDof(node + 2, dof[1]);
107
108
109
      }
    }
    
110
    if (mesh->getNumberOfDofs(FACE)) {
111
      int node = mesh->getNode(FACE);    
112
113
114
115
116
      
      /****************************************************************************/
      /*  set pointers to those dof's that are handed on from the parant          */
      /****************************************************************************/
      
117
118
      child[0]->setDof(node + 3, const_cast<int*>(el->getDof(node + 1)));
      child[1]->setDof(node + 3, const_cast<int*>(el->getDof(node + 0)));
119
120
121
122
123
      
      /****************************************************************************/
      /*  get new dof for the common face of child0 and child1                    */
      /****************************************************************************/
      
124
      DegreeOfFreedom *newDOF = mesh->getDof(FACE);
125
126
      child[0]->setDof(node, static_cast<int*>(newDOF));
      child[1]->setDof(node, static_cast<int*>(newDOF));
127
128
    }
    
129
    if (mesh->getNumberOfDofs(CENTER)) {
130
      int node = mesh->getNode(CENTER);
131
132
      child[0]->setDof(node, const_cast<int*>(mesh->getDof(CENTER)));
      child[1]->setDof(node, const_cast<int*>(mesh->getDof(CENTER)));
133
    }
134

135
    if (mesh->getNumberOfDofs(EDGE) || mesh->getNumberOfDofs(FACE))
136
137
      fillPatchConnectivity(ref_list, index);
  }
138

139

140
141
142
143
  void RefinementManager3d::fillPatchConnectivity(RCNeighbourList* ref_list, 
						  int index)
  {
    FUNCNAME("RefinementManager3d::fillPatchConnectivity");
Thomas Witkowski's avatar
Thomas Witkowski committed
144

145
    Element *el = ref_list->getElement(index), *neigh;
146
147
148
    int el_type = ref_list->getType(index);
    int n_type = 0;
    int dir, adjc, i, j, i_neigh, j_neigh;
149
    int node0, node1, oppVertex = 0;
150

151
152
153
    for (dir = 0; dir < 2; dir++) {
      neigh = ref_list->getNeighbourElement(index, dir);
      if (neigh) {
154
	n_type = ref_list->getType(ref_list->getNeighbourNr(index, dir));
155
	oppVertex = ref_list->getOppVertex(index, dir);
156
157
      }

158
      if (!neigh || neigh->isLeaf()) {
159
160

	/****************************************************************************/
161
162
163
	/*  get new dof's in the midedge of the face of el and for the two midpoints*/
	/*  of the sub-faces. If face is an interior face those pointers have to be */
	/*  adjusted by the neighbour element also (see below)                      */
164
165
	/****************************************************************************/

166
	if (mesh->getNumberOfDofs(EDGE)) {
167
168
169
	  node0 = node1 = mesh->getNode(EDGE);
	  node0 += Tetrahedron::nChildEdge[el_type][0][dir];
	  node1 += Tetrahedron::nChildEdge[el_type][1][dir];
170
	  DegreeOfFreedom *newDOF = mesh->getDof(EDGE);
171
172
	  (const_cast<Element*>(el->getFirstChild()))->setDof(node0, newDOF);
	  (const_cast<Element*>(el->getSecondChild()))->setDof(node1, newDOF);
173
	}
174
	if (mesh->getNumberOfDofs(FACE)) {
175
	  node0 = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][0][dir];
176
	  (const_cast<Element*>(el->getFirstChild()))->setDof(node0, mesh->getDof(FACE));
177
	  node1 = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][1][dir];
178
	  (const_cast<Element*>(el->getSecondChild()))->setDof(node1, mesh->getDof(FACE));
179
180
	}
      } else {     /*   if (!neigh  ||  !neigh->child[0])                         */
181
	/****************************************************************************/
182
183
	/*  interior face and neighbour has been refined, look for position at the  */
	/*  refinement edge                                                         */
184
	/****************************************************************************/
185
      
186
	if (el->getDof(0) == neigh->getDof(0)) {
187
188
189
190
191
192
193
194
195
196
	  /****************************************************************************/
	  /* same position at refinement edge                                         */
	  /****************************************************************************/
	  adjc = 0;
	} else {
	  /****************************************************************************/
	  /* different position at refinement edge                                    */
	  /****************************************************************************/
	  adjc = 1;
	}
197

198
199
	for (i = 0; i < 2; i++) {
	  j = Tetrahedron::adjacentChild[adjc][i];
200

201
	  i_neigh = Tetrahedron::nChildFace[el_type][i][dir];
202
	  j_neigh = Tetrahedron::nChildFace[n_type][j][oppVertex - 2];
203

204
205
206
207
	  /****************************************************************************/
	  /*  adjust dof pointer in the edge in the common face of el and neigh and   */
	  /*  the dof pointer in the sub-face child_i-child_j (allocated by neigh!)   */
	  /****************************************************************************/
208

209
	  if (mesh->getNumberOfDofs(EDGE)) {
210
	    node0 = mesh->getNode(EDGE) + Tetrahedron::nChildEdge[el_type][i][dir];
211
	    node1 = mesh->getNode(EDGE) + Tetrahedron::nChildEdge[n_type][j][oppVertex - 2];
212

213
	    TEST_EXIT_DBG(neigh->getChild(j)->getDof(node1))
214
215
	      ("no dof on neighbour %d at node %d\n", 
	       neigh->getChild(j)->getIndex(), node1);
216

217
	    (const_cast<Element*>(el->getChild(i)))->
218
	      setDof(node0, const_cast<int*>(neigh->getChild(j)->getDof(node1)));
219
	  }
220
	  if (mesh->getNumberOfDofs(FACE)) {
221
222
	    node0 = mesh->getNode(FACE) + i_neigh;
	    node1 = mesh->getNode(FACE) + j_neigh;
223

224
	    TEST_EXIT_DBG(neigh->getChild(j)->getDof(node1))
Thomas Witkowski's avatar
Thomas Witkowski committed
225
	      ("No DOF on neighbour %d at node %d!\n",
226
	       neigh->getChild(j)->getIndex(), node1);
227

228
	    (const_cast<Element*>(el->getChild(i)))->
229
	      setDof(node0, const_cast<int*>(neigh->getChild(j)->getDof(node1)));
230
231
	  }

232
233
234
	}  /*   for (i = 0; i < 2; i++)                                       */
      }    /*   else of   if (!neigh  ||  !neigh->child[0])                   */
    }      /*   for (dir = 0; dir < 2; dir++)                                 */
235
236
237
238
239
240
241
  }


  int RefinementManager3d::newCoordsFct(ElInfo *el_info)
  {
    Element *el = el_info->getElement();
    DegreeOfFreedom *edge[2];
242
    ElInfo *elinfo = el_info;
243
244
245
    int dow = Global::getGeo(WORLD);
    Projection *projector = el_info->getProjection(0);

246
247
    if (!projector || projector->getType() != VOLUME_PROJECTION)
      projector = el_info->getProjection(4);    
248
249

    if (el->getFirstChild() && projector && (!el->isNewCoordSet())) {
Thomas Witkowski's avatar
Thomas Witkowski committed
250
      WorldVector<double> *new_coord = new WorldVector<double>;
251

252
253
      for (int j = 0; j < dow; j++)
	(*new_coord)[j] = (el_info->getCoord(0)[j] + el_info->getCoord(1)[j]) * 0.5;
254
255
256
257
258
259
260
261
262

      projector->project(*new_coord);

      el->setNewCoord(new_coord);
      /****************************************************************************/
      /*  now, information should be passed on to patch neighbours...             */
      /*  get the refinement patch                                                */
      /****************************************************************************/
      refList->setElement(0, el);
263
264
      refList->setElType(0, el_info->getType());
      int n_neigh = 1;
265

266
      for (int i = 0; i < 2; i++)
267
	edge[i] = const_cast<int*>(el_info->getElement()->getDof(i));
268

269
      if (getRefinePatch(&elinfo, edge, 0, refList, &n_neigh)) {
270
271
	// Domain's boundary was reached while looping around the refinement edge.

272
273
274
	getRefinePatch(&elinfo, edge, 1, refList, &n_neigh);
      }

275
      for (int i = 1; i < n_neigh; i++) {            /* start with 1, as list[0]=el */
276
277
278
279
280
281
	TEST(!refList->getElement(i)->isNewCoordSet())
	  ("non-nil new_coord in el %d ref_list[%d] el %d (n_neigh=%d)\n",
	   el->getIndex(), i, refList->getElement(i)->getIndex(), n_neigh);
	
	refList->getElement(i)->setNewCoord(el->getNewCoord());
      }
282
283
284
285
286
    }

    return 0;
  }

287

288
  void RefinementManager3d::setNewCoords(int macroEl)
289
290
  {
    if (refList)
Thomas Witkowski's avatar
Thomas Witkowski committed
291
      delete refList;
292
 
Thomas Witkowski's avatar
Thomas Witkowski committed
293
    refList = new RCNeighbourList(mesh->getMaxEdgeNeigh());
294
295
296
297
298
299
300
301
302
303
304
305
306
    ElInfo *elInfo;

    if (macroEl == -1)
      elInfo = stack->traverseFirst(mesh, -1, 
				    Mesh::CALL_EVERY_EL_PREORDER | 
				    Mesh::FILL_BOUND | Mesh::FILL_COORDS | 
				    Mesh::FILL_NEIGH);
    else
      elInfo = stack->traverseFirstOneMacro(mesh, macroEl, -1,
					    Mesh::CALL_EVERY_EL_PREORDER | 
					    Mesh::FILL_BOUND | Mesh::FILL_COORDS | 
					    Mesh::FILL_NEIGH);
    
307

308
309
310
    while (elInfo) {
      newCoordsFct(elInfo);
      elInfo = stack->traverseNext(elInfo);
311
312
313
314
315
316
317
318
    }
  }


  DegreeOfFreedom RefinementManager3d::refinePatch(DegreeOfFreedom *edge[2], 
						   RCNeighbourList* refineList,
						   int n_neigh, bool bound)
  {
319
320
    Tetrahedron *el = 
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(refineList->getElement(0)));
321
322
323
324
325
326
327
    /* first element in the list */
    DegreeOfFreedom *dof[3] = {NULL, NULL, NULL};

    /****************************************************************************/
    /*  get new dof's in the refinement edge                                    */
    /****************************************************************************/

328
    dof[0] = mesh->getDof(VERTEX);
329
330
    mesh->incrementNumberOfVertices(1);
  
331
    if (mesh->getNumberOfDofs(EDGE)) {
332
333
      dof[1] = mesh->getDof(EDGE);
      dof[2] = mesh->getDof(EDGE);
334
    }
335

336
    for (int i = 0; i < n_neigh; i++)
337
338
339
340
341
342
343
      bisectTetrahedron(refineList, i, dof, edge);

    /****************************************************************************/
    /*  if there are functions to interpolate data to the finer grid, do so     */
    /****************************************************************************/
    int iadmin;
    int nrAdmin = mesh->getNumberOfDOFAdmin();
344
    for (iadmin = 0; iadmin < nrAdmin; iadmin++) {
345
      std::list<DOFIndexedBase*>::iterator it;
346
      DOFAdmin* admin = const_cast<DOFAdmin*>(&mesh->getDofAdmin(iadmin));
347
      std::list<DOFIndexedBase*>::iterator end = admin->endDOFIndexed();
348
      for (it = admin->beginDOFIndexed(); it != end; it++)
349
350
351
	(*it)->refineInterpol(*refineList, n_neigh);
    }

352
353
354
355
356
    if (!mesh->queryCoarseDOFs()) {
      /****************************************************************************/
      /*  if there should be no dof information on interior leaf elements remove  */
      /*  dofs from edges, faces and the centers of parents                       */
      /****************************************************************************/
357
      if (mesh->getNumberOfDofs(EDGE)) {
358
	/****************************************************************************/
359
	/*  remove dof of the midpoint of the common refinement edge                */
360
361
	/****************************************************************************/

362
	el = dynamic_cast<Tetrahedron*>(const_cast<Element*>(refineList->getElement(0)));
363
	mesh->freeDof(const_cast<int*>(el->getDof(mesh->getNode(EDGE))), EDGE);
364
      }
365
      
366
367
368
      if (mesh->getNumberOfDofs(EDGE) || 
	  mesh->getNumberOfDofs(FACE) ||
	  mesh->getNumberOfDofs(CENTER)) {
369
	for (int i = 0; i < n_neigh; i++)
370
371
372
	  refineList->removeDOFParent(i);
      }
    }
373
374
375
376
377
378
379


    /****************************************************************************/
    /*  update the number of edges and faces; depends whether refinement edge   */
    /*  is a boundary edge or not                                               */
    /****************************************************************************/

380
381
    if (bound) {
      mesh->incrementNumberOfEdges(n_neigh + 2);
382
      mesh->incrementNumberOfFaces(2 * n_neigh + 1);
383
384
    } else {
      mesh->incrementNumberOfEdges(n_neigh + 1);
385
      mesh->incrementNumberOfFaces(2 * n_neigh);
386
387
    }
    
388
389
390
391
392
393
    return dof[0][0];
  }


  int RefinementManager3d::getRefinePatch(ElInfo **el_info, 
					  DegreeOfFreedom *edge[2], 
394
					  int direction,
395
396
397
					  RCNeighbourList* refineList, 
					  int *n_neigh)
  {
398
    FUNCNAME("RefinementManager3d::getRefinePatch()");
399

400
    int localNeighbour = 3 - direction;
401
402
403
    Tetrahedron *el = 
      dynamic_cast<Tetrahedron*>(const_cast<Element*>((*el_info)->getElement()));
    
404
    if ((*el_info)->getNeighbour(localNeighbour) == NULL)
Thomas Witkowski's avatar
Thomas Witkowski committed
405
      return 1;    
406
  
407
408
409
410
411
412
413
    int oppVertex = (*el_info)->getOppVertex(localNeighbour);
    int testIndex = (*el_info)->getNeighbour(localNeighbour)->getIndex();
    ElInfo *neighInfo = stack->traverseNeighbour3d((*el_info), localNeighbour);
    int neighElType = neighInfo->getType();

    TEST_EXIT_DBG(neighInfo->getElement()->getIndex() == testIndex)
      ("Should not happen!\n");
414

415
    Tetrahedron *neigh = 
416
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(neighInfo->getElement())); 
417
418
419
    int vertices = mesh->getGeo(VERTEX);

    while (neigh != el) {
420
421
422
423
      // === Determine the common edge of the element and its neighbour. ===

      int edgeDof0, edgeDof1;
      for (edgeDof0 = 0; edgeDof0 < vertices; edgeDof0++)
424
	if (neigh->getDof(edgeDof0) == edge[0])
425
	  break;
426
      for (edgeDof1 = 0; edgeDof1 < vertices; edgeDof1++)
427
	if (neigh->getDof(edgeDof1) == edge[1])
428
	  break;
429

430
431
      if (edgeDof0 > 3 || edgeDof1 > 3) {
	for (edgeDof0 = 0; edgeDof0 < vertices; edgeDof0++)
432
	  if (mesh->associated(neigh->getDof(edgeDof0, 0), edge[0][0]))  
433
	    break;
434
	for (edgeDof1 = 0; edgeDof1 < vertices; edgeDof1++)
435
	  if (mesh->associated(neigh->getDof(edgeDof1, 0), edge[1][0]))  
436
	    break;
437
	    
438
439
	if (edgeDof0 > 3 || edgeDof1 > 3) {
	  for (edgeDof0 = 0; edgeDof0 < vertices; edgeDof0++)
440
	    if (mesh->indirectlyAssociated(neigh->getDof(edgeDof0, 0), edge[0][0]))  
441
	      break;
442
	  for (edgeDof1 = 0; edgeDof1 < vertices; edgeDof1++)
443
	    if (mesh->indirectlyAssociated(neigh->getDof(edgeDof1, 0), edge[1][0]))  
444
	      break;
445
	    
446
447
	  TEST_EXIT_DBG(edgeDof0 < vertices && edgeDof1 < vertices)
 	    ("dof %d or dof %d not found on element %d with nodes (%d %d %d %d)\n", 
448
449
 	     edge[0][0], edge[1][0], neigh->getIndex(), neigh->getDof(0,0),
 	     neigh->getDof(1, 0), neigh->getDof(2, 0), neigh->getDof(3, 0));
450
451
452
453

	}
      }

454
      TEST_EXIT_DBG(edgeDof0 < vertices && edgeDof1 < vertices)
455
	("dof %d or dof %d not found on element %d with nodes (%d %d %d %d)\n", 
456
457
	 edge[0][0], edge[1][0], neigh->getIndex(), neigh->getDof(0, 0),
	 neigh->getDof(1, 0), neigh->getDof(2, 0), neigh->getDof(3, 0));
458
459
460
461
462
463

      int edgeNo = Tetrahedron::edgeOfDofs[edgeDof0][edgeDof1];

      if (edgeNo) {
	// Only 0 can be a compatible commen refinement edge. Thus, neigh has not 
	// a compatible refinement edge. Refine it first.
464

465
	neigh->setMark(max(neigh->getMark(), 1));
466
467
468
469
470
471

	neighInfo = refineFunction(neighInfo);

	// === Now, go to a child at the edge and determine the opposite vertex ===
	// === for  this child; continue the looping around the edge with this  ===
	// === element.                                                         ===
472
	
473
474
475
	neighInfo = stack->traverseNext(neighInfo);
	neighElType = neighInfo->getType();
	bool reverseMode = stack->getTraverseFlag().isSet(Mesh::CALL_REVERSE_MODE);
476
	
477
	switch (edgeNo) {
478
	case 1: 
479
480
481
482
483
484
	  if (reverseMode) {
	    neighInfo = stack->traverseNext(neighInfo);
	    neighElType = neighInfo->getType();
	  }
	    
	  oppVertex = (oppVertex == 1 ? 3 : 2);
485
486
	  break;
	case 2: 
487
488
489
490
491
492
	  if (reverseMode) {
	    neighInfo = stack->traverseNext(neighInfo);
	    neighElType = neighInfo->getType();
	  }

	  oppVertex = (oppVertex == 2 ? 1 : 3);
493
494
	  break;
	case 3: 
495
496
497
498
499
500
501
	  if (!reverseMode) {
	    neighInfo = stack->traverseNext(neighInfo);
	    neighElType = neighInfo->getType();
	  }

	  if (neighElType != 1)
	    oppVertex = (oppVertex == 0 ? 3 : 2);
502
	  else
503
	    oppVertex = (oppVertex == 0 ? 3 : 1);
504
505
	  break;
	case 4:
506
507
508
509
510
511
512
	  if (!reverseMode) {
	    neighInfo = stack->traverseNext(neighInfo);
	    neighElType = neighInfo->getType();
	  }

	  if (neighElType != 1)
	    oppVertex = (oppVertex == 0 ? 3 : 1);
513
	  else
514
	    oppVertex = (oppVertex == 0 ? 3 : 2);
515
516
	  break;
	case 5:
517
518
519
520
521
	  if (neighElType != 1) {
	    if (!reverseMode) {
	      neighInfo = stack->traverseNext(neighInfo);
	      neighElType = neighInfo->getType();
	    }
522
	  }
523
	  oppVertex = 3;
524
	  break;
525
	}
526

527
	neigh = 
528
	  dynamic_cast<Tetrahedron*>(const_cast<Element*>(neighInfo->getElement()));
529
      } else {
530
531
532
	// Neigh is compatible devisible. Put neigh to the list of patch elements 
	// and go to next neighbour.

533
534
	TEST_EXIT_DBG(*n_neigh < mesh->getMaxEdgeNeigh())
	  ("too many neighbours %d in refine patch\n", *n_neigh);
535
		
536
	refineList->setElement(*n_neigh, neigh);
537
	refineList->setElType(*n_neigh, neighElType);
538
	
539
	// We have to go back to the starting element via oppVertex values.
540

541
	refineList->setOppVertex(*n_neigh, 0, oppVertex); 
542
	
543
544
545
546
547
548
549
550
551
552
553
554
555
556
	(*n_neigh)++;

	int i = (oppVertex != 3 ? 3 : 2);

	if (neighInfo->getNeighbour(i)) {
	  oppVertex = neighInfo->getOppVertex(i);
	  int testIndex = neighInfo->getNeighbour(i)->getIndex();

	  neighInfo = stack->traverseNeighbour3d(neighInfo, i);

	  TEST_EXIT_DBG(neighInfo->getElement()->getIndex() == testIndex)
	    ("Should not happen!\n");

	  neighElType = neighInfo->getType();
557
	  neigh = 
558
559
	    dynamic_cast<Tetrahedron*>(const_cast<Element*>(neighInfo->getElement()));
	} else {
560
	  break;
561
	}
562
      }
563
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
564
   
565
    if (neigh == el) {
566
567
      (*el_info) = neighInfo;

568
      return 0;
569
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
570
571

   
572
    // The domain's boundary is reached. Loop back to the starting el.
573
    
574
575
    int i = *n_neigh - 1;
    oppVertex = refineList->getOppVertex(i, 0);
576
    do {
577
      TEST_EXIT_DBG(neighInfo->getNeighbour(oppVertex) && i > 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
578
	("While looping back domains boundary was reached or i == 0\n");
579
580
581
582
583
584
585
      oppVertex = refineList->getOppVertex(i--, 0);

      int testIndex = neighInfo->getNeighbour(oppVertex)->getIndex();
      neighInfo = stack->traverseNeighbour3d(neighInfo, oppVertex);

      int edgeDof0, edgeDof1;
      for (edgeDof0 = 0; edgeDof0 < vertices; edgeDof0++)
586
	if (neigh->getDof(edgeDof0) == edge[0])
587
588
	  break;
      for (edgeDof1 = 0; edgeDof1 < vertices; edgeDof1++)
589
	if (neigh->getDof(edgeDof1) == edge[1])
590
591
592
593
594
	  break;

      TEST_EXIT_DBG(neighInfo->getElement()->getIndex() == testIndex)
	("Should not happen!\n");
    } while (neighInfo->getElement() != el);
Thomas Witkowski's avatar
Thomas Witkowski committed
595

596
    (*el_info) = neighInfo;
597
598
    
    return 1;
599
600
601
602
603
604
  }


  ElInfo* RefinementManager3d::refineFunction(ElInfo* el_info)
  {
    FUNCNAME("RefinementManager3d::refineFunction()");
Thomas Witkowski's avatar
Thomas Witkowski committed
605

606
607
608
    if (el_info->getElement()->getMark() <= 0)  
      return el_info;

609
    int bound = false;
610
    DegreeOfFreedom *edge[2];
611

612
613
    // === Get memory for a list of all elements at the refinement edge. ===

614
    RCNeighbourList *ref_list = new RCNeighbourList(mesh->getMaxEdgeNeigh());
615
    ref_list->setElType(0, el_info->getType());
616
    ref_list->setElement(0, el_info->getElement());
617
618
619
620
621
    int n_neigh = 1;

    if (el_info->getProjection(0) && 
	el_info->getProjection(0)->getType() == VOLUME_PROJECTION)
      newCoords = true;
622
623


624
625
    // === Give the refinement edge the right orientation. ===

626
627
628
    if (el_info->getElement()->getDof(0, 0) < el_info->getElement()->getDof(1, 0)) {
      edge[0] = const_cast<DegreeOfFreedom*>(el_info->getElement()->getDof(0));
      edge[1] = const_cast<DegreeOfFreedom*>(el_info->getElement()->getDof(1));
629
    } else {
630
631
      edge[1] = const_cast<DegreeOfFreedom*>(el_info->getElement()->getDof(0));
      edge[0] = const_cast<DegreeOfFreedom*>(el_info->getElement()->getDof(1));
632
    }
633

634
635
636

    // === Traverse and refine the refinement patch. ====

637
    if (getRefinePatch(&el_info, edge, 0, ref_list, &n_neigh)) {      
638
      // Domain's boundary was reached while looping around the refinement edge
639
640
641
      getRefinePatch(&el_info, edge, 1, ref_list, &n_neigh);
      bound = true;
    }
642

643
    // fill neighbour information inside the patch in the refinement list
644
645
    ref_list->getNeighOnPatch(n_neigh, bound);

646
647

    // ============ Check for periodic boundary ============
648
649
650
651
652
653
654
655
656
657
658

    DegreeOfFreedom *next_edge[2];
    DegreeOfFreedom *first_edge[2] = {edge[0], edge[1]};
    DegreeOfFreedom *last_edge[2] = {NULL, NULL};
    int n_neigh_periodic;

    DegreeOfFreedom newDOF = -1;
    DegreeOfFreedom lastNewDOF = -1;
    DegreeOfFreedom firstNewDOF = -1;

    RCNeighbourList *periodicList;
659
660
    std::map<int, VertexVector*>::iterator it;
    std::map<int, VertexVector*>::iterator end = mesh->getPeriodicAssociations().end();
661

Thomas Witkowski's avatar
Thomas Witkowski committed
662
    while (edge[0] != NULL) {
663
664
665
666
667
      periodicList = ref_list->periodicSplit(edge, 
					     next_edge,
					     &n_neigh,
					     &n_neigh_periodic);

668
      TEST_EXIT_DBG(periodicList)("periodicList = NULL\n");
669

670
      newDOF = refinePatch(edge, periodicList, n_neigh_periodic, bound);
671

672
      if (firstNewDOF == -1)
673
674
	firstNewDOF = newDOF;

675
676
677
678
      if (lastNewDOF != -1) {
	for (it = mesh->getPeriodicAssociations().begin(); it != end; ++it) {
	  if (it->second) {
	    if (((*(it->second))[edge[0][0]] == last_edge[0][0] &&
Thomas Witkowski's avatar
Thomas Witkowski committed
679
680
681
682
683
684
		 (*(it->second))[edge[1][0]] == last_edge[1][0]) || 
		((*(it->second))[edge[0][0]] == last_edge[1][0] &&
		 (*(it->second))[edge[1][0]] == last_edge[0][0])) {
	      (*(it->second))[lastNewDOF] = newDOF;
	      (*(it->second))[newDOF] = lastNewDOF;
	    } 
685
686
687
688
689
690
691
692
693
694
695
696
	  }
	}
      }
      lastNewDOF = newDOF;

      last_edge[0] = edge[0];
      last_edge[1] = edge[1];

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

697
698
699
700
    if (lastNewDOF != firstNewDOF) {
      for (it = mesh->getPeriodicAssociations().begin(); it != end; ++it) {
	if (it->second) {
	  if (((*(it->second))[first_edge[0][0]] == last_edge[0][0] &&
Thomas Witkowski's avatar
Thomas Witkowski committed
701
702
703
704
705
706
	       (*(it->second))[first_edge[1][0]] == last_edge[1][0]) || 
	      ((*(it->second))[first_edge[0][0]] == last_edge[1][0] &&
	       (*(it->second))[first_edge[1][0]] == last_edge[0][0])) {	    
	    (*(it->second))[lastNewDOF] = firstNewDOF;
	    (*(it->second))[firstNewDOF] = lastNewDOF;
	  }
707
708
709
710
711
712
	}
      }
    }
  
    stack->update();

Thomas Witkowski's avatar
Thomas Witkowski committed
713
    delete ref_list;
714

Thomas Witkowski's avatar
Thomas Witkowski committed
715
    return el_info;
716
717
718
  }

}