RefinementManager3d.cc 25.2 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
#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"
28
#include "Debug.h"
29
30
31
32
33
34
35
36

namespace AMDiS {

  void RefinementManager3d::bisectTetrahedron(RCNeighbourList* ref_list, 
					      int index,
					      DegreeOfFreedom* dof[3],
					      DegreeOfFreedom *edge[2])
  {
37
38
39
    Tetrahedron *el = 
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(ref_list->getElement(index)));
    Tetrahedron *child[2];
40
    int el_type = ref_list->getType(index);
41
42
43
44

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

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

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

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

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

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

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

152

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


  int RefinementManager3d::newCoordsFct(ElInfo *el_info)
  {
253
254
    FUNCNAME("RefinementManager3d::newCoordsFct()");

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

261
262
    if (!projector || projector->getType() != VOLUME_PROJECTION)
      projector = el_info->getProjection(4);    
263
264

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

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

      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);
278
279
      refList->setElType(0, el_info->getType());
      int n_neigh = 1;
280

281
      for (int i = 0; i < 2; i++)
282
	edge[i] = const_cast<int*>(el_info->getElement()->getDof(i));
283

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

287
288
289
	getRefinePatch(&elinfo, edge, 1, refList, &n_neigh);
      }

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

    return 0;
  }

302

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

323
324
325
    while (elInfo) {
      newCoordsFct(elInfo);
      elInfo = stack->traverseNext(elInfo);
326
327
328
329
330
331
332
333
    }
  }


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

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

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

351
    for (int i = 0; i < n_neigh; i++)
352
353
354
355
356
357
358
      bisectTetrahedron(refineList, i, dof, edge);

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

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

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


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

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


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

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

430
    Tetrahedron *neigh = 
431
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(neighInfo->getElement())); 
432
433
434
    int vertices = mesh->getGeo(VERTEX);

    while (neigh != el) {
435
436
437
438
      // === Determine the common edge of the element and its neighbour. ===

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

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

	}
      }

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

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

480
	neigh->setMark(std::max(neigh->getMark(), 1));
481
482
483
484
485
486

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

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

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

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

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

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

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

583
      return 0;
584
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
585
586

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

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

611
    (*el_info) = neighInfo;
612
613
    
    return 1;
614
615
616
617
618
619
  }


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

621
622
623
    if (el_info->getElement()->getMark() <= 0)  
      return el_info;

624
    int bound = false;
625
    DegreeOfFreedom *edge[2];
626

627
628
    // === Get memory for a list of all elements at the refinement edge. ===

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

    if (el_info->getProjection(0) && 
	el_info->getProjection(0)->getType() == VOLUME_PROJECTION)
      newCoords = true;
637
638


639
640
    // === Give the refinement edge the right orientation. ===

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

649
650
651

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

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

658
    // fill neighbour information inside the patch in the refinement list
659
660
    ref_list->getNeighOnPatch(n_neigh, bound);

661
662

    // ============ Check for periodic boundary ============
663
664
665
666
667
668
669
670
671
672
673

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

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

683
      TEST_EXIT_DBG(periodicList)("periodicList = NULL\n");
684

685
      newDOF = refinePatch(edge, periodicList, n_neigh_periodic, bound);
686

687
      if (firstNewDOF == -1)
688
689
	firstNewDOF = newDOF;

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
728
    delete ref_list;
729

Thomas Witkowski's avatar
Thomas Witkowski committed
730
    return el_info;
731
732
733
  }

}