RefinementManager3d.cc 29.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
#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

namespace AMDiS {

32
33
  FixRefinementPatch::ConnectedEdges FixRefinementPatch::connectedEdges(0);

34
  void RefinementManager3d::bisectTetrahedron(RCNeighbourList& refineList, 
35
36
37
38
					      int index,
					      DegreeOfFreedom* dof[3],
					      DegreeOfFreedom *edge[2])
  {
39
40
    FUNCNAME("RefinementManager3d::bisectTetrahedron()");

41
    Tetrahedron *el = 
42
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(refineList.getElement(index)));
43
    Tetrahedron *child[2];
44
    int el_type = refineList.getType(index);
45
46
47
48

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

63
64
    if (child[0]->getMark() > 0)
      doMoreRecursiveRefine = true;
65
66

    int n_vertices = mesh->getGeo(VERTEX);
67
68
    child[0]->setDof(n_vertices - 1, dof[0]);
    child[1]->setDof(n_vertices - 1, dof[0]);
69

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

88
    if (mesh->getNumberOfDofs(EDGE)) {
89
      int node = mesh->getNode(EDGE);
90

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

152
    if (mesh->getNumberOfDofs(EDGE) || mesh->getNumberOfDofs(FACE))
153
      fillPatchConnectivity(refineList, index);
154
  }
155

156

157
  void RefinementManager3d::fillPatchConnectivity(RCNeighbourList &refineList,
158
159
160
						  int index)
  {
    FUNCNAME("RefinementManager3d::fillPatchConnectivity");
Thomas Witkowski's avatar
Thomas Witkowski committed
161

162
163
    Element *el = refineList.getElement(index);
    int el_type = refineList.getType(index);
164
    int n_type = 0;
165
    int dir, adjc, i_neigh, j_neigh;
166
    int node0, node1, oppVertex = 0;
167

168
    for (dir = 0; dir < 2; dir++) {
169
      Element *neigh = refineList.getNeighbourElement(index, dir);
170
      if (neigh) {
171
172
	n_type = refineList.getType(refineList.getNeighbourNr(index, dir));
	oppVertex = refineList.getOppVertex(index, dir);
173
174
      }

175
      if (!neigh || neigh->isLeaf()) {
176
	/****************************************************************************/
177
178
179
	/*  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)                      */
180
181
	/****************************************************************************/

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

210
211
	for (int i = 0; i < 2; i++) {
	  int 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
  void RefinementManager3d::newCoordsFct(ElInfo *elInfo, RCNeighbourList &refineList)
251
  {
252
253
    FUNCNAME("RefinementManager3d::newCoordsFct()");

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

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

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

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

      projector->project(*new_coord);

      el->setNewCoord(new_coord);
      /****************************************************************************/
      /*  now, information should be passed on to patch neighbours...             */
      /*  get the refinement patch                                                */
      /****************************************************************************/
276
      refineList.setElement(0, el);
277
      refineList.setElType(0, elInfo->getType());
278
      int n_neigh = 1;
279

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

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

286
	getRefinePatch(&elinfo, edge, 1, refineList, &n_neigh);
287
288
      }

289
      for (int i = 1; i < n_neigh; i++) {            /* start with 1, as list[0]=el */
290
291
292
	TEST(!refineList.getElement(i)->isNewCoordSet())
	  ("non-nil new_coord in el %d refineList[%d] el %d (n_neigh=%d)\n",
	   el->getIndex(), i, refineList.getElement(i)->getIndex(), n_neigh);
293
	
294
	refineList.getElement(i)->setNewCoord(el->getNewCoord());
295
      }
296
297
298
    }
  }

299

300
  void RefinementManager3d::setNewCoords(int macroEl)
301
  {
302
    RCNeighbourList refineList(mesh->getMaxEdgeNeigh());
303
304
305
306
307
308
309
310
311
312
313
314
315
    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);
    
316

317
    while (elInfo) {
318
      newCoordsFct(elInfo, refineList);
319
      elInfo = stack->traverseNext(elInfo);
320
321
322
323
324
    }
  }


  DegreeOfFreedom RefinementManager3d::refinePatch(DegreeOfFreedom *edge[2], 
325
						   RCNeighbourList &refineList,
326
327
						   int n_neigh, bool bound)
  {
328
329
    FUNCNAME("RefinementManager3d::refinePatch()");

330
    Tetrahedron *el = 
331
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(refineList.getElement(0)));
332
333
334
335
336
337
338
    /* first element in the list */
    DegreeOfFreedom *dof[3] = {NULL, NULL, NULL};

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

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

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

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

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

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


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

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


401
402
403
404
405
  bool RefinementManager3d::getRefinePatch(ElInfo **elInfo, 
					   DegreeOfFreedom *edge[2], 
					   int direction,
					   RCNeighbourList &refineList, 
					   int *n_neigh)
406
  {
407
    FUNCNAME("RefinementManager3d::getRefinePatch()");
408

409
    int localNeighbour = 3 - direction;
410
    Tetrahedron *el = 
411
      dynamic_cast<Tetrahedron*>(const_cast<Element*>((*elInfo)->getElement()));
412

413
414
    if ((*elInfo)->getNeighbour(localNeighbour) == NULL)
      return true;    
415
  
416
    int oppVertex = (*elInfo)->getOppVertex(localNeighbour);
417
#if DEBUG
418
    int testIndex = (*elInfo)->getNeighbour(localNeighbour)->getIndex();
419
#endif
420
    ElInfo *neighInfo = stack->traverseNeighbour3d((*elInfo), localNeighbour);
421
422
423
424
    int neighElType = neighInfo->getType();

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

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

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

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

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

462
	  TEST_EXIT_DBG(edgeDof1 < vertices)
463
 	    ("DOF %d not found on element %d with nodes (%d %d %d %d)\n", 
464
465
 	     edge[1][0], neigh->getIndex(), neigh->getDof(0, 0),
 	     neigh->getDof(1, 0), neigh->getDof(2, 0), neigh->getDof(3, 0));
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
552
	refineList.setElement(*n_neigh, neigh);
	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
	(*n_neigh)++;

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

	if (neighInfo->getNeighbour(i)) {
	  oppVertex = neighInfo->getOppVertex(i);
564
    
565
#if DEBUG
566
	  int testIndex = neighInfo->getNeighbour(i)->getIndex();
567
#endif
568
569
570
571
572
573
574

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

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

	  neighElType = neighInfo->getType();
575
	  neigh = 
576
577
	    dynamic_cast<Tetrahedron*>(const_cast<Element*>(neighInfo->getElement()));
	} else {
578
	  break;
579
	}
580
      }
581
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
582
   
583
    if (neigh == el) {
584
      (*elInfo) = neighInfo;
585

586
      return false;
587
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
588
589

   
590
    // The domain's boundary is reached. Loop back to the starting el.
591
    
592
    int i = *n_neigh - 1;
593
    oppVertex = refineList.getOppVertex(i, 0);
594
    do {
595
      TEST_EXIT_DBG(neighInfo->getNeighbour(oppVertex) && i > 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
596
	("While looping back domains boundary was reached or i == 0\n");
597
      oppVertex = refineList.getOppVertex(i--, 0);
598

599
#if DEBUG
600
      int testIndex = neighInfo->getNeighbour(oppVertex)->getIndex();
601
#endif
602
603
604
605
      neighInfo = stack->traverseNeighbour3d(neighInfo, oppVertex);

      int edgeDof0, edgeDof1;
      for (edgeDof0 = 0; edgeDof0 < vertices; edgeDof0++)
606
	if (neigh->getDof(edgeDof0) == edge[0])
607
608
	  break;
      for (edgeDof1 = 0; edgeDof1 < vertices; edgeDof1++)
609
	if (neigh->getDof(edgeDof1) == edge[1])
610
611
612
613
614
	  break;

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

616
    (*elInfo) = neighInfo;
617
    
618
    return true;
619
620
621
  }


622
  ElInfo* RefinementManager3d::refineFunction(ElInfo* elInfo)
623
624
  {
    FUNCNAME("RefinementManager3d::refineFunction()");
Thomas Witkowski's avatar
Thomas Witkowski committed
625

626
    Element *el = elInfo->getElement();
627
628

    if (el->getMark() <= 0)  
629
630
      return elInfo;

631
    int bound = false;
632
    DegreeOfFreedom *edge[2];
633

634
635
    // === Get memory for a list of all elements at the refinement edge. ===

636
    RCNeighbourList refineList(mesh->getMaxEdgeNeigh());
637
    refineList.setElType(0, elInfo->getType());
638
    refineList.setElement(0, el);
639
640
    int n_neigh = 1;

641
642
    if (elInfo->getProjection(0) && 
	elInfo->getProjection(0)->getType() == VOLUME_PROJECTION)
643
      newCoords = true;
644
645


646
647
    // === Give the refinement edge the right orientation. ===

648
649
650
    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));
651
    } else {
652
653
      edge[1] = const_cast<DegreeOfFreedom*>(el->getDof(0));
      edge[0] = const_cast<DegreeOfFreedom*>(el->getDof(1));
654
    }
655

656
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
657
658
659
660
    Element *otherEl = NULL;
    int otherEdge = -1;
    FixRefinementPatch::getOtherEl(stack, &otherEl, otherEdge);
#endif
661
662
663

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

664
    if (getRefinePatch(&elInfo, edge, 0, refineList, &n_neigh)) {      
665
      // Domain's boundary was reached while looping around the refinement edge
666
      getRefinePatch(&elInfo, edge, 1, refineList, &n_neigh);
667
668
      bound = true;
    }
669

670
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
    // === If the refinement edge must be fixed, add also the other part of this ===
    // === edge to the refinement patch.                                         ===

    if (otherEl) {
      TraverseStack stack2;
      ElInfo *elInfo2 = 
	stack2.traverseFirstOneMacro(mesh, otherEl->getIndex(), -1, 
				     Mesh::CALL_LEAF_EL | 
				     Mesh::FILL_NEIGH | 
				     Mesh::FILL_BOUND);

      bool foundEdge = false;

      while (elInfo2) {
	for (int i = 0; i < 6; i++) {
	  DofEdge edge2 = elInfo2->getElement()->getEdge(i);
	  if (edge2.first == *(edge[0]) &&
	      edge2.second == *(edge[1])) {

690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
	    // We found the edge on the other element side on leaf level. Two
	    // cases may occure: either it is already a refinement edge, i.e,
	    // it has local edge number 0, or it is not a refinement edge. In
	    // the first case, we are fine and this leaf element may be set
	    // to all elements on the refinement edge. Otherwise, the element
	    // must be refined at least once to get a refinement edge.

	    if (i == 0) {
	      // Edge is refinement edge, so add it to refine list.

	      refineList.setElType(n_neigh, elInfo2->getType());
	      refineList.setElement(n_neigh, elInfo2->getElement());
	      n_neigh++;
	      
	      TraverseStack *tmpStack = stack;
	      stack = &stack2;
	      
	      if (getRefinePatch(&elInfo2, edge, 0, refineList, &n_neigh)) {
		getRefinePatch(&elInfo2, edge, 1, refineList, &n_neigh);
		bound = true;
	      }
	      
	      stack = tmpStack;
713
	      foundEdge = true;
714
715
716
	      break;
	    } else {
	      // Edge i not refinement edge, so refine the element further.
717

718
	      Element *el2 = elInfo2->getElement();
719
720
721
722
723
724
725
726
727
728
	      el2->setMark(std::max(el2->getMark(), 1));

	      TraverseStack *tmpStack = stack;
	      stack = &stack2;

	      elInfo2 = refineFunction(elInfo2);

	      stack = tmpStack;
	      break;
	    }
729
730
731
	  }
	}

732
733
734
	if (foundEdge)
	  break;

735
736
737
738
739
740
741
	elInfo2 = stack2.traverseNext(elInfo2);
      }

      TEST_EXIT_DBG(foundEdge)("Should not happen!\n");
    }
#endif

742
    // fill neighbour information inside the patch in the refinement list
743
    refineList.fillNeighbourRelations(n_neigh, bound);
744

745
    // ============ Check for periodic boundary ============
746
747
748
749
750
751

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

752
753
    DegreeOfFreedom lastNewDof = -1;
    DegreeOfFreedom firstNewDof = -1;
754

755
    RCNeighbourList periodicList;
756

Thomas Witkowski's avatar
Thomas Witkowski committed
757
    while (edge[0] != NULL) {
758
759
760
      refineList.periodicSplit(edge, next_edge, 
			       &n_neigh, &n_neigh_periodic, 
			       periodicList);
761

762
      DegreeOfFreedom newDof = refinePatch(edge, periodicList, n_neigh_periodic, bound);
763

764
765
      if (firstNewDof == -1)
	firstNewDof = newDof;
766

767
768
769
      if (lastNewDof != -1) {
	for (std::map<int, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
	     it != mesh->getPeriodicAssociations().end(); ++it) {
770
771
	  if (it->second) {
	    if (((*(it->second))[edge[0][0]] == last_edge[0][0] &&
Thomas Witkowski's avatar
Thomas Witkowski committed
772
773
774
		 (*(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])) {
775
776
	      (*(it->second))[lastNewDof] = newDof;
	      (*(it->second))[newDof] = lastNewDof;
Thomas Witkowski's avatar
Thomas Witkowski committed
777
	    } 
778
779
780
	  }
	}
      }
781
      lastNewDof = newDof;
782
783
784
785
786
787
788
789

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

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

790
791
792
    if (lastNewDof != firstNewDof) {
      for (std::map<int, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
	   it != mesh->getPeriodicAssociations().end(); ++it) {
793
794
	if (it->second) {
	  if (((*(it->second))[first_edge[0][0]] == last_edge[0][0] &&
Thomas Witkowski's avatar
Thomas Witkowski committed
795
796
797
	       (*(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])) {	    
798
799
	    (*(it->second))[lastNewDof] = firstNewDof;
	    (*(it->second))[firstNewDof] = lastNewDof;
Thomas Witkowski's avatar
Thomas Witkowski committed
800
	  }
801
802
803
	}
      }
    }
804

805
806
    stack->update();

807
808
809
810
811
    return elInfo;
  }


  void FixRefinementPatch::getOtherEl(TraverseStack *stack, 
812
813
				      Element **otherEl, 
				      int &otherEdge)
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
  {
    FUNCNAME("FixRefinementPatch::getOtherEl()");

    if (!FixRefinementPatch::connectedEdges.empty()) {
      // === Get stack of current traverse. ===
      std::vector<ElInfo*> elInfos;
      std::vector<int> infos;      
      int stackUsed = stack->getStackData(elInfos, infos);         
      int checkIndex = stackUsed;
      int localEdgeNo = 0;
      ElInfo *elInfo = elInfos[stackUsed];

      // === Calculate the refinement edge number on the macro element level. ===
      
      for (int i = 0; i < elInfo->getLevel(); i++) {
	TEST_EXIT_DBG(checkIndex >= 1)("Should not happen!\n");
	
	int parentType = elInfos[checkIndex - 1]->getType();
	int ithParentChild = 0;
	if (elInfos[checkIndex - 1]->getElement()->getChild(1) ==
	    elInfos[checkIndex]->getElement())
	  ithParentChild = 1;
	localEdgeNo = Tetrahedron::childEdge[parentType][ithParentChild][localEdgeNo];
	
	if (localEdgeNo >= 4) {
	  localEdgeNo = -1;
	  break;
	}
	
	checkIndex--;
      }
      
      // If the refinement edge is part of an edge on the macro level, we must check
      // if the refinement edge is part of an edge that must be fixed.
      if (localEdgeNo >= 0) {
	TEST_EXIT_DBG(elInfos[checkIndex]->getLevel() == 0)("Should not happen!\n");
	TEST_EXIT_DBG(elInfos[checkIndex]->getElement()->getIndex() ==
		      elInfo->getMacroElement()->getIndex())("Should not happen!\n");
	TEST_EXIT_DBG(localEdgeNo <= 5)("Should not happen!\n");
	
	for (unsigned int i = 0; i < FixRefinementPatch::connectedEdges.size(); i++) {
	  if (FixRefinementPatch::connectedEdges[i].first.first->getIndex() ==
	      elInfos[checkIndex]->getElement()->getIndex() &&
	      FixRefinementPatch::connectedEdges[i].first.second ==
	      localEdgeNo) {
	    // Okay, we have found that this edge must be fixed.
	    *otherEl = FixRefinementPatch::connectedEdges[i].second.first;
	    otherEdge = FixRefinementPatch::connectedEdges[i].second.second;
	    
	    break;
	  }
	}
      }
    }
868
869
  }
}