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 25.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
22
23
24
25
26
27
28
29
30
31
32
33
34
35
#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])
  {
36
37
38
    Tetrahedron *el = 
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(ref_list->getElement(index)));
    Tetrahedron *child[2];
39
    int el_type = ref_list->getType(index);
40
41
42
43

    child[0] = dynamic_cast<Tetrahedron*>(mesh->createNewElement(el));
    child[1] = dynamic_cast<Tetrahedron*>(mesh->createNewElement(el));
  
44
    int mark = max(0, el->getMark() - 1);
45
46
47
48
49
50
51
52
53
54
55
56
57
    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]);

58
59
    if (child[0]->getMark() > 0)
      doMoreRecursiveRefine = true;
60
61

    int n_vertices = mesh->getGeo(VERTEX);
62
63
    child[0]->setDof(n_vertices - 1, dof[0]);
    child[1]->setDof(n_vertices - 1, dof[0]);
64

65
    for (int i = 0; i < n_vertices - 1; i++) {
66
      child[0]->
67
	setDof(i, const_cast<int*>(el->getDof(Tetrahedron::childVertex[el_type][0][i])));
68
      child[1]->
69
	setDof(i, const_cast<int*>(el->getDof(Tetrahedron::childVertex[el_type][1][i])));
70
    }
71
72
73
74
75
76
77
78
79
80
81
82
    /****************************************************************************/
    /*  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                                                              */
    /****************************************************************************/

83
    if (mesh->getNumberOfDofs(EDGE)) {
84
      int node = mesh->getNode(EDGE);
85

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

147
    if (mesh->getNumberOfDofs(EDGE) || mesh->getNumberOfDofs(FACE))
148
149
      fillPatchConnectivity(ref_list, index);
  }
150

151

152
153
154
155
  void RefinementManager3d::fillPatchConnectivity(RCNeighbourList* ref_list, 
						  int index)
  {
    FUNCNAME("RefinementManager3d::fillPatchConnectivity");
Thomas Witkowski's avatar
Thomas Witkowski committed
156

157
    Element *el = ref_list->getElement(index), *neigh;
158
159
160
    int el_type = ref_list->getType(index);
    int n_type = 0;
    int dir, adjc, i, j, i_neigh, j_neigh;
161
    int node0, node1, oppVertex = 0;
162

163
164
165
    for (dir = 0; dir < 2; dir++) {
      neigh = ref_list->getNeighbourElement(index, dir);
      if (neigh) {
166
	n_type = ref_list->getType(ref_list->getNeighbourNr(index, dir));
167
	oppVertex = ref_list->getOppVertex(index, dir);
168
169
      }

170
      if (!neigh || neigh->isLeaf()) {
171
172

	/****************************************************************************/
173
174
175
	/*  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)                      */
176
177
	/****************************************************************************/

178
	if (mesh->getNumberOfDofs(EDGE)) {
179
180
181
	  node0 = node1 = mesh->getNode(EDGE);
	  node0 += Tetrahedron::nChildEdge[el_type][0][dir];
	  node1 += Tetrahedron::nChildEdge[el_type][1][dir];
182
	  DegreeOfFreedom *newDOF = mesh->getDof(EDGE);
183
184
	  (const_cast<Element*>(el->getFirstChild()))->setDof(node0, newDOF);
	  (const_cast<Element*>(el->getSecondChild()))->setDof(node1, newDOF);
185
	}
186
	if (mesh->getNumberOfDofs(FACE)) {
187
	  node0 = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][0][dir];
188
	  (const_cast<Element*>(el->getFirstChild()))->setDof(node0, mesh->getDof(FACE));
189
	  node1 = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][1][dir];
190
	  (const_cast<Element*>(el->getSecondChild()))->setDof(node1, mesh->getDof(FACE));
191
192
	}
      } else {     /*   if (!neigh  ||  !neigh->child[0])                         */
193
	/****************************************************************************/
194
195
	/*  interior face and neighbour has been refined, look for position at the  */
	/*  refinement edge                                                         */
196
	/****************************************************************************/
197
      
198
	if (el->getDof(0) == neigh->getDof(0)) {
199
200
201
202
203
204
205
206
207
208
	  /****************************************************************************/
	  /* same position at refinement edge                                         */
	  /****************************************************************************/
	  adjc = 0;
	} else {
	  /****************************************************************************/
	  /* different position at refinement edge                                    */
	  /****************************************************************************/
	  adjc = 1;
	}
209

210
211
	for (i = 0; i < 2; i++) {
	  j = Tetrahedron::adjacentChild[adjc][i];
212

213
	  i_neigh = Tetrahedron::nChildFace[el_type][i][dir];
214
	  j_neigh = Tetrahedron::nChildFace[n_type][j][oppVertex - 2];
215

216
217
218
219
	  /****************************************************************************/
	  /*  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!)   */
	  /****************************************************************************/
220

221
	  if (mesh->getNumberOfDofs(EDGE)) {
222
	    node0 = mesh->getNode(EDGE) + Tetrahedron::nChildEdge[el_type][i][dir];
223
	    node1 = mesh->getNode(EDGE) + Tetrahedron::nChildEdge[n_type][j][oppVertex - 2];
224

225
	    TEST_EXIT_DBG(neigh->getChild(j)->getDof(node1))
226
227
	      ("no dof on neighbour %d at node %d\n", 
	       neigh->getChild(j)->getIndex(), node1);
228

229
	    (const_cast<Element*>(el->getChild(i)))->
230
	      setDof(node0, const_cast<int*>(neigh->getChild(j)->getDof(node1)));
231
	  }
232
	  if (mesh->getNumberOfDofs(FACE)) {
233
234
	    node0 = mesh->getNode(FACE) + i_neigh;
	    node1 = mesh->getNode(FACE) + j_neigh;
235

236
	    TEST_EXIT_DBG(neigh->getChild(j)->getDof(node1))
Thomas Witkowski's avatar
Thomas Witkowski committed
237
	      ("No DOF on neighbour %d at node %d!\n",
238
	       neigh->getChild(j)->getIndex(), node1);
239

240
	    (const_cast<Element*>(el->getChild(i)))->
241
	      setDof(node0, const_cast<int*>(neigh->getChild(j)->getDof(node1)));
242
243
	  }

244
245
246
	}  /*   for (i = 0; i < 2; i++)                                       */
      }    /*   else of   if (!neigh  ||  !neigh->child[0])                   */
    }      /*   for (dir = 0; dir < 2; dir++)                                 */
247
248
249
250
251
252
253
  }


  int RefinementManager3d::newCoordsFct(ElInfo *el_info)
  {
    Element *el = el_info->getElement();
    DegreeOfFreedom *edge[2];
254
    ElInfo *elinfo = el_info;
255
256
257
    int dow = Global::getGeo(WORLD);
    Projection *projector = el_info->getProjection(0);

258
259
    if (!projector || projector->getType() != VOLUME_PROJECTION)
      projector = el_info->getProjection(4);    
260
261

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

264
265
      for (int j = 0; j < dow; j++)
	(*new_coord)[j] = (el_info->getCoord(0)[j] + el_info->getCoord(1)[j]) * 0.5;
266
267
268
269
270
271
272
273
274

      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);
275
276
      refList->setElType(0, el_info->getType());
      int n_neigh = 1;
277

278
      for (int i = 0; i < 2; i++)
279
	edge[i] = const_cast<int*>(el_info->getElement()->getDof(i));
280

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

284
285
286
	getRefinePatch(&elinfo, edge, 1, refList, &n_neigh);
      }

287
      for (int i = 1; i < n_neigh; i++) {            /* start with 1, as list[0]=el */
288
289
290
291
292
293
	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());
      }
294
295
296
297
298
    }

    return 0;
  }

299

300
  void RefinementManager3d::setNewCoords(int macroEl)
301
302
  {
    if (refList)
Thomas Witkowski's avatar
Thomas Witkowski committed
303
      delete refList;
304
 
Thomas Witkowski's avatar
Thomas Witkowski committed
305
    refList = new RCNeighbourList(mesh->getMaxEdgeNeigh());
306
307
308
309
310
311
312
313
314
315
316
317
318
    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);
    
319

320
321
322
    while (elInfo) {
      newCoordsFct(elInfo);
      elInfo = stack->traverseNext(elInfo);
323
324
325
326
327
328
329
330
    }
  }


  DegreeOfFreedom RefinementManager3d::refinePatch(DegreeOfFreedom *edge[2], 
						   RCNeighbourList* refineList,
						   int n_neigh, bool bound)
  {
331
332
    Tetrahedron *el = 
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(refineList->getElement(0)));
333
334
335
336
337
338
339
    /* first element in the list */
    DegreeOfFreedom *dof[3] = {NULL, NULL, NULL};

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

340
    dof[0] = mesh->getDof(VERTEX);
341
342
    mesh->incrementNumberOfVertices(1);
  
343
    if (mesh->getNumberOfDofs(EDGE)) {
344
345
      dof[1] = mesh->getDof(EDGE);
      dof[2] = mesh->getDof(EDGE);
346
    }
347

348
    for (int i = 0; i < n_neigh; i++)
349
350
351
352
353
354
355
      bisectTetrahedron(refineList, i, dof, edge);

    /****************************************************************************/
    /*  if there are functions to interpolate data to the finer grid, do so     */
    /****************************************************************************/
    int iadmin;
    int nrAdmin = mesh->getNumberOfDOFAdmin();
356
    for (iadmin = 0; iadmin < nrAdmin; iadmin++) {
357
      std::list<DOFIndexedBase*>::iterator it;
358
      DOFAdmin* admin = const_cast<DOFAdmin*>(&mesh->getDofAdmin(iadmin));
359
      std::list<DOFIndexedBase*>::iterator end = admin->endDOFIndexed();
360
      for (it = admin->beginDOFIndexed(); it != end; it++)
361
362
363
	(*it)->refineInterpol(*refineList, n_neigh);
    }

364
365
366
367
368
    if (!mesh->queryCoarseDOFs()) {
      /****************************************************************************/
      /*  if there should be no dof information on interior leaf elements remove  */
      /*  dofs from edges, faces and the centers of parents                       */
      /****************************************************************************/
369
      if (mesh->getNumberOfDofs(EDGE)) {
370
	/****************************************************************************/
371
	/*  remove dof of the midpoint of the common refinement edge                */
372
373
	/****************************************************************************/

374
	el = dynamic_cast<Tetrahedron*>(const_cast<Element*>(refineList->getElement(0)));
375
	mesh->freeDof(const_cast<int*>(el->getDof(mesh->getNode(EDGE))), EDGE);
376
      }
377
      
378
379
380
      if (mesh->getNumberOfDofs(EDGE) || 
	  mesh->getNumberOfDofs(FACE) ||
	  mesh->getNumberOfDofs(CENTER)) {
381
	for (int i = 0; i < n_neigh; i++)
382
383
384
	  refineList->removeDOFParent(i);
      }
    }
385
386
387
388
389
390
391


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

392
393
    if (bound) {
      mesh->incrementNumberOfEdges(n_neigh + 2);
394
      mesh->incrementNumberOfFaces(2 * n_neigh + 1);
395
396
    } else {
      mesh->incrementNumberOfEdges(n_neigh + 1);
397
      mesh->incrementNumberOfFaces(2 * n_neigh);
398
399
    }
    
400
401
402
403
404
405
    return dof[0][0];
  }


  int RefinementManager3d::getRefinePatch(ElInfo **el_info, 
					  DegreeOfFreedom *edge[2], 
406
					  int direction,
407
408
409
					  RCNeighbourList* refineList, 
					  int *n_neigh)
  {
410
    FUNCNAME("RefinementManager3d::getRefinePatch()");
411

412
    int localNeighbour = 3 - direction;
413
414
415
    Tetrahedron *el = 
      dynamic_cast<Tetrahedron*>(const_cast<Element*>((*el_info)->getElement()));
    
416
    if ((*el_info)->getNeighbour(localNeighbour) == NULL)
Thomas Witkowski's avatar
Thomas Witkowski committed
417
      return 1;    
418
  
419
420
421
422
423
424
425
    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");
426

427
    Tetrahedron *neigh = 
428
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(neighInfo->getElement())); 
429
430
431
    int vertices = mesh->getGeo(VERTEX);

    while (neigh != el) {
432
433
434
435
      // === Determine the common edge of the element and its neighbour. ===

      int edgeDof0, edgeDof1;
      for (edgeDof0 = 0; edgeDof0 < vertices; edgeDof0++)
436
	if (neigh->getDof(edgeDof0) == edge[0])
437
	  break;
438
      for (edgeDof1 = 0; edgeDof1 < vertices; edgeDof1++)
439
	if (neigh->getDof(edgeDof1) == edge[1])
440
	  break;
441

442
443
      if (edgeDof0 > 3 || edgeDof1 > 3) {
	for (edgeDof0 = 0; edgeDof0 < vertices; edgeDof0++)
444
	  if (mesh->associated(neigh->getDof(edgeDof0, 0), edge[0][0]))  
445
	    break;
446
	for (edgeDof1 = 0; edgeDof1 < vertices; edgeDof1++)
447
	  if (mesh->associated(neigh->getDof(edgeDof1, 0), edge[1][0]))  
448
	    break;
449
	    
450
451
	if (edgeDof0 > 3 || edgeDof1 > 3) {
	  for (edgeDof0 = 0; edgeDof0 < vertices; edgeDof0++)
452
	    if (mesh->indirectlyAssociated(neigh->getDof(edgeDof0, 0), edge[0][0]))  
453
	      break;
454
	  for (edgeDof1 = 0; edgeDof1 < vertices; edgeDof1++)
455
	    if (mesh->indirectlyAssociated(neigh->getDof(edgeDof1, 0), edge[1][0]))  
456
	      break;
457
	    
458
459
	  TEST_EXIT_DBG(edgeDof0 < vertices && edgeDof1 < vertices)
 	    ("dof %d or dof %d not found on element %d with nodes (%d %d %d %d)\n", 
460
461
 	     edge[0][0], edge[1][0], neigh->getIndex(), neigh->getDof(0,0),
 	     neigh->getDof(1, 0), neigh->getDof(2, 0), neigh->getDof(3, 0));
462
463
464
465

	}
      }

466
      TEST_EXIT_DBG(edgeDof0 < vertices && edgeDof1 < vertices)
467
	("dof %d or dof %d not found on element %d with nodes (%d %d %d %d)\n", 
468
469
	 edge[0][0], edge[1][0], neigh->getIndex(), neigh->getDof(0, 0),
	 neigh->getDof(1, 0), neigh->getDof(2, 0), neigh->getDof(3, 0));
470
471
472
473
474
475

      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.
476

477
	neigh->setMark(max(neigh->getMark(), 1));
478
479
480
481
482
483

	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.                                                         ===
484
	
485
486
487
	neighInfo = stack->traverseNext(neighInfo);
	neighElType = neighInfo->getType();
	bool reverseMode = stack->getTraverseFlag().isSet(Mesh::CALL_REVERSE_MODE);
488
	
489
	switch (edgeNo) {
490
	case 1: 
491
492
493
494
495
496
	  if (reverseMode) {
	    neighInfo = stack->traverseNext(neighInfo);
	    neighElType = neighInfo->getType();
	  }
	    
	  oppVertex = (oppVertex == 1 ? 3 : 2);
497
498
	  break;
	case 2: 
499
500
501
502
503
504
	  if (reverseMode) {
	    neighInfo = stack->traverseNext(neighInfo);
	    neighElType = neighInfo->getType();
	  }

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

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

	  if (neighElType != 1)
	    oppVertex = (oppVertex == 0 ? 3 : 1);
525
	  else
526
	    oppVertex = (oppVertex == 0 ? 3 : 2);
527
528
	  break;
	case 5:
529
530
531
532
533
	  if (neighElType != 1) {
	    if (!reverseMode) {
	      neighInfo = stack->traverseNext(neighInfo);
	      neighElType = neighInfo->getType();
	    }
534
	  }
535
	  oppVertex = 3;
536
	  break;
537
	}
538

539
	neigh = 
540
	  dynamic_cast<Tetrahedron*>(const_cast<Element*>(neighInfo->getElement()));
541
      } else {
542
543
544
	// Neigh is compatible devisible. Put neigh to the list of patch elements 
	// and go to next neighbour.

545
546
	TEST_EXIT_DBG(*n_neigh < mesh->getMaxEdgeNeigh())
	  ("too many neighbours %d in refine patch\n", *n_neigh);
547
		
548
	refineList->setElement(*n_neigh, neigh);
549
	refineList->setElType(*n_neigh, neighElType);
550
	
551
	// We have to go back to the starting element via oppVertex values.
552

553
	refineList->setOppVertex(*n_neigh, 0, oppVertex); 
554
	
555
556
557
558
559
560
561
562
563
564
565
566
567
568
	(*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();
569
	  neigh = 
570
571
	    dynamic_cast<Tetrahedron*>(const_cast<Element*>(neighInfo->getElement()));
	} else {
572
	  break;
573
	}
574
      }
575
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
576
   
577
    if (neigh == el) {
578
579
      (*el_info) = neighInfo;

580
      return 0;
581
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
582
583

   
584
    // The domain's boundary is reached. Loop back to the starting el.
585
    
586
587
    int i = *n_neigh - 1;
    oppVertex = refineList->getOppVertex(i, 0);
588
    do {
589
      TEST_EXIT_DBG(neighInfo->getNeighbour(oppVertex) && i > 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
590
	("While looping back domains boundary was reached or i == 0\n");
591
592
593
594
595
596
597
      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++)
598
	if (neigh->getDof(edgeDof0) == edge[0])
599
600
	  break;
      for (edgeDof1 = 0; edgeDof1 < vertices; edgeDof1++)
601
	if (neigh->getDof(edgeDof1) == edge[1])
602
603
604
605
606
	  break;

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

608
    (*el_info) = neighInfo;
609
610
    
    return 1;
611
612
613
614
615
616
  }


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

618
619
620
    if (el_info->getElement()->getMark() <= 0)  
      return el_info;

621
    int bound = false;
622
    DegreeOfFreedom *edge[2];
623

624
625
    // === Get memory for a list of all elements at the refinement edge. ===

626
    RCNeighbourList *ref_list = new RCNeighbourList(mesh->getMaxEdgeNeigh());
627
    ref_list->setElType(0, el_info->getType());
628
    ref_list->setElement(0, el_info->getElement());
629
630
631
632
633
    int n_neigh = 1;

    if (el_info->getProjection(0) && 
	el_info->getProjection(0)->getType() == VOLUME_PROJECTION)
      newCoords = true;
634
635


636
637
    // === Give the refinement edge the right orientation. ===

638
639
640
    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));
641
    } else {
642
643
      edge[1] = const_cast<DegreeOfFreedom*>(el_info->getElement()->getDof(0));
      edge[0] = const_cast<DegreeOfFreedom*>(el_info->getElement()->getDof(1));
644
    }
645

646
647
648

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

649
    if (getRefinePatch(&el_info, edge, 0, ref_list, &n_neigh)) {      
650
      // Domain's boundary was reached while looping around the refinement edge
651
652
653
      getRefinePatch(&el_info, edge, 1, ref_list, &n_neigh);
      bound = true;
    }
654

655
    // fill neighbour information inside the patch in the refinement list
656
657
    ref_list->getNeighOnPatch(n_neigh, bound);

658
659

    // ============ Check for periodic boundary ============
660
661
662
663
664
665
666
667
668
669
670

    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;
671
672
    std::map<int, VertexVector*>::iterator it;
    std::map<int, VertexVector*>::iterator end = mesh->getPeriodicAssociations().end();
673

Thomas Witkowski's avatar
Thomas Witkowski committed
674
    while (edge[0] != NULL) {
675
676
677
678
679
      periodicList = ref_list->periodicSplit(edge, 
					     next_edge,
					     &n_neigh,
					     &n_neigh_periodic);

680
      TEST_EXIT_DBG(periodicList)("periodicList = NULL\n");
681

682
      newDOF = refinePatch(edge, periodicList, n_neigh_periodic, bound);
683

684
      if (firstNewDOF == -1)
685
686
	firstNewDOF = newDOF;

687
688
689
690
      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
691
692
693
694
695
696
		 (*(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;
	    } 
697
698
699
700
701
702
703
704
705
706
707
708
	  }
	}
      }
      lastNewDOF = newDOF;

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

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

709
710
711
712
    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
713
714
715
716
717
718
	       (*(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;
	  }
719
720
721
722
723
724
	}
      }
    }
  
    stack->update();

Thomas Witkowski's avatar
Thomas Witkowski committed
725
    delete ref_list;
726

Thomas Witkowski's avatar
Thomas Witkowski committed
727
    return el_info;
728
729
730
  }

}