RefinementManager3d.cc 29.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
20
21


22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#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"
37
#include "Debug.h"
38
39
40

namespace AMDiS {

41
42
  FixRefinementPatch::ConnectedEdges FixRefinementPatch::connectedEdges(0);

43
  void RefinementManager3d::bisectTetrahedron(RCNeighbourList& refineList, 
44
45
46
47
					      int index,
					      DegreeOfFreedom* dof[3],
					      DegreeOfFreedom *edge[2])
  {
48
    Tetrahedron *el = 
49
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(refineList.getElement(index)));
50
    Tetrahedron *child[2];
51
    int el_type = refineList.getType(index);
52
53
54
55

    child[0] = dynamic_cast<Tetrahedron*>(mesh->createNewElement(el));
    child[1] = dynamic_cast<Tetrahedron*>(mesh->createNewElement(el));
  
56
    int mark = std::max(0, el->getMark() - 1);
57
58
59
60
61
62
63
64
65
66
67
68
69
    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]);

70
71
    if (child[0]->getMark() > 0)
      doMoreRecursiveRefine = true;
72
73

    int n_vertices = mesh->getGeo(VERTEX);
74
75
    child[0]->setDof(n_vertices - 1, dof[0]);
    child[1]->setDof(n_vertices - 1, dof[0]);
76

77
    for (int i = 0; i < n_vertices - 1; i++) {
78
      child[0]->
79
	setDof(i, const_cast<int*>(el->getDof(Tetrahedron::childVertex[el_type][0][i])));
80
      child[1]->
81
	setDof(i, const_cast<int*>(el->getDof(Tetrahedron::childVertex[el_type][1][i])));
82
    }
83
84
85
86
87
88
89
90
91
92
93
94
    /****************************************************************************/
    /*  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                                                              */
    /****************************************************************************/

95
    if (mesh->getNumberOfDofs(EDGE)) {
96
      int node = mesh->getNode(EDGE);
97

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

159
    if (mesh->getNumberOfDofs(EDGE) || mesh->getNumberOfDofs(FACE))
160
      fillPatchConnectivity(refineList, index);
161
  }
162

163

164
  void RefinementManager3d::fillPatchConnectivity(RCNeighbourList &refineList,
165
166
						  int index)
  {
167
    FUNCNAME_DBG("RefinementManager3d::fillPatchConnectivity");
Thomas Witkowski's avatar
Thomas Witkowski committed
168

169
170
    Element *el = refineList.getElement(index);
    int el_type = refineList.getType(index);
171
    int n_type = 0;
172
    int dir, adjc, i_neigh, j_neigh;
173
    int node0, node1, oppVertex = 0;
174

175
    for (dir = 0; dir < 2; dir++) {
176
      Element *neigh = refineList.getNeighbourElement(index, dir);
177
      if (neigh) {
178
179
	n_type = refineList.getType(refineList.getNeighbourNr(index, dir));
	oppVertex = refineList.getOppVertex(index, dir);
180
181
      }

182
      if (!neigh || neigh->isLeaf()) {
183
	/****************************************************************************/
184
185
186
	/*  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)                      */
187
188
	/****************************************************************************/

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

217
218
	for (int i = 0; i < 2; i++) {
	  int j = Tetrahedron::adjacentChild[adjc][i];
219

220
	  i_neigh = Tetrahedron::nChildFace[el_type][i][dir];
221
	  j_neigh = Tetrahedron::nChildFace[n_type][j][oppVertex - 2];
222

223
224
225
226
	  /****************************************************************************/
	  /*  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!)   */
	  /****************************************************************************/
227

228
	  if (mesh->getNumberOfDofs(EDGE)) {
229
	    node0 = mesh->getNode(EDGE) + Tetrahedron::nChildEdge[el_type][i][dir];
230
	    node1 = mesh->getNode(EDGE) + Tetrahedron::nChildEdge[n_type][j][oppVertex - 2];
231

232
	    TEST_EXIT_DBG(neigh->getChild(j)->getDof(node1))
233
234
	      ("no dof on neighbour %d at node %d\n", 
	       neigh->getChild(j)->getIndex(), node1);
235

236
	    (const_cast<Element*>(el->getChild(i)))->
237
	      setDof(node0, const_cast<int*>(neigh->getChild(j)->getDof(node1)));
238
	  }
239
	  if (mesh->getNumberOfDofs(FACE)) {
240
241
	    node0 = mesh->getNode(FACE) + i_neigh;
	    node1 = mesh->getNode(FACE) + j_neigh;
242

243
	    TEST_EXIT_DBG(neigh->getChild(j)->getDof(node1))
Thomas Witkowski's avatar
Thomas Witkowski committed
244
	      ("No DOF on neighbour %d at node %d!\n",
245
	       neigh->getChild(j)->getIndex(), node1);
246

247
	    (const_cast<Element*>(el->getChild(i)))->
248
	      setDof(node0, const_cast<int*>(neigh->getChild(j)->getDof(node1)));
249
250
	  }

251
252
253
	}  /*   for (i = 0; i < 2; i++)                                       */
      }    /*   else of   if (!neigh  ||  !neigh->child[0])                   */
    }      /*   for (dir = 0; dir < 2; dir++)                                 */
254
255
256
  }


257
  void RefinementManager3d::newCoordsFct(ElInfo *elInfo, RCNeighbourList &refineList)
258
  {
259
260
    FUNCNAME("RefinementManager3d::newCoordsFct()");

261
    Element *el = elInfo->getElement();
262
    DegreeOfFreedom *edge[2];
263
    ElInfo *elinfo = elInfo;
264
    int dow = Global::getGeo(WORLD);
265
    Projection *projector = elInfo->getProjection(0);
266

267
    if (!projector || projector->getType() != VOLUME_PROJECTION)
268
      projector = elInfo->getProjection(4);    
269
270

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

273
      for (int j = 0; j < dow; j++)
274
	(*new_coord)[j] = (elInfo->getCoord(0)[j] + elInfo->getCoord(1)[j]) * 0.5;
275
276
277
278
279
280
281
282

      projector->project(*new_coord);

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

287
      for (int i = 0; i < 2; i++)
288
	edge[i] = const_cast<int*>(elInfo->getElement()->getDof(i));
289

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

293
	getRefinePatch(&elinfo, edge, 1, refineList, &n_neigh);
294
295
      }

296
      for (int i = 1; i < n_neigh; i++) {            /* start with 1, as list[0]=el */
297
298
299
	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);
300
	
301
	refineList.getElement(i)->setNewCoord(el->getNewCoord());
302
      }
303
304
305
    }
  }

306

307
  void RefinementManager3d::setNewCoords(int macroEl)
308
  {
309
    RCNeighbourList refineList(mesh->getMaxEdgeNeigh());
310
311
312
313
314
315
316
317
318
319
320
321
322
    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);
    
323

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


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

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

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

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

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

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

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


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

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


406
407
408
409
410
  bool RefinementManager3d::getRefinePatch(ElInfo **elInfo, 
					   DegreeOfFreedom *edge[2], 
					   int direction,
					   RCNeighbourList &refineList, 
					   int *n_neigh)
411
  {
412
    FUNCNAME_DBG("RefinementManager3d::getRefinePatch()");
413

414
    int localNeighbour = 3 - direction;
415
    Tetrahedron *el = 
416
      dynamic_cast<Tetrahedron*>(const_cast<Element*>((*elInfo)->getElement()));
417

418
419
    if ((*elInfo)->getNeighbour(localNeighbour) == NULL)
      return true;    
420
  
421
    int oppVertex = (*elInfo)->getOppVertex(localNeighbour);
422
#if DEBUG
423
    int testIndex = (*elInfo)->getNeighbour(localNeighbour)->getIndex();
424
#endif
425
    ElInfo *neighInfo = stack->traverseNeighbour3d((*elInfo), localNeighbour);
426
427
428
429
    int neighElType = neighInfo->getType();

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

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

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

      int edgeDof0, edgeDof1;
439
440
441

      // get local/elementwise DOF indices of Start and End Vertices of EDGE
      // on Neighbour Element
442
      for (edgeDof0 = 0; edgeDof0 < vertices; edgeDof0++)
443
	if (neigh->getDof(edgeDof0) == edge[0])
444
	  break;
445
      for (edgeDof1 = 0; edgeDof1 < vertices; edgeDof1++)
446
	if (neigh->getDof(edgeDof1) == edge[1])
447
	  break;
448

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

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

477
      TEST_EXIT_DBG(edgeDof0 < vertices && edgeDof1 < vertices)
478
	("DOF %d or DOF %d not found on element %d with nodes (%d %d %d %d)\n", 
479
480
	 edge[0][0], edge[1][0], neigh->getIndex(), neigh->getDof(0, 0),
	 neigh->getDof(1, 0), neigh->getDof(2, 0), neigh->getDof(3, 0));
481
482
483
484
485
486

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

488
	neigh->setMark(std::max(neigh->getMark(), 1));
489
490
491
492
493
494

	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.                                                         ===
495
	
496
497
498
	neighInfo = stack->traverseNext(neighInfo);
	neighElType = neighInfo->getType();
	bool reverseMode = stack->getTraverseFlag().isSet(Mesh::CALL_REVERSE_MODE);
499
	
500
	switch (edgeNo) {
501
	case 1: 
502
503
504
505
506
507
	  if (reverseMode) {
	    neighInfo = stack->traverseNext(neighInfo);
	    neighElType = neighInfo->getType();
	  }
	    
	  oppVertex = (oppVertex == 1 ? 3 : 2);
508
509
	  break;
	case 2: 
510
511
512
513
514
515
	  if (reverseMode) {
	    neighInfo = stack->traverseNext(neighInfo);
	    neighElType = neighInfo->getType();
	  }

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

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

	  if (neighElType != 1)
	    oppVertex = (oppVertex == 0 ? 3 : 1);
536
	  else
537
	    oppVertex = (oppVertex == 0 ? 3 : 2);
538
539
	  break;
	case 5:
540
541
542
543
544
	  if (neighElType != 1) {
	    if (!reverseMode) {
	      neighInfo = stack->traverseNext(neighInfo);
	      neighElType = neighInfo->getType();
	    }
545
	  }
546
	  oppVertex = 3;
547
	  break;
548
	}
549

550
	neigh = 
551
	  dynamic_cast<Tetrahedron*>(const_cast<Element*>(neighInfo->getElement()));
552
      } else {
553
554
555
	// Neigh is compatible devisible. Put neigh to the list of patch elements 
	// and go to next neighbour.

556
557
	TEST_EXIT_DBG(*n_neigh < mesh->getMaxEdgeNeigh())
	  ("too many neighbours %d in refine patch\n", *n_neigh);
558
		
559
560
	refineList.setElement(*n_neigh, neigh);
	refineList.setElType(*n_neigh, neighElType);
561
	
562
	// We have to go back to the starting element via oppVertex values.
563

564
	refineList.setOppVertex(*n_neigh, 0, oppVertex); 
565
	
566
567
568
569
570
571
	(*n_neigh)++;

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

	if (neighInfo->getNeighbour(i)) {
	  oppVertex = neighInfo->getOppVertex(i);
572
    
573
#if DEBUG
574
	  int testIndex = neighInfo->getNeighbour(i)->getIndex();
575
#endif
576
577
578
579
580
581
582

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

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

	  neighElType = neighInfo->getType();
583
	  neigh = 
584
585
	    dynamic_cast<Tetrahedron*>(const_cast<Element*>(neighInfo->getElement()));
	} else {
586
	  break;
587
	}
588
      }
589
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
590
   
591
    if (neigh == el) {
592
      (*elInfo) = neighInfo;
593

594
      return false;
595
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
596
597

   
598
    // The domain's boundary is reached. Loop back to the starting el.
599
    
600
    int i = *n_neigh - 1;
601
    oppVertex = refineList.getOppVertex(i, 0);
602
    do {
603
      TEST_EXIT_DBG(neighInfo->getNeighbour(oppVertex) && i > 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
604
	("While looping back domains boundary was reached or i == 0\n");
605
      oppVertex = refineList.getOppVertex(i--, 0);
606

607
#if DEBUG
608
      int testIndex = neighInfo->getNeighbour(oppVertex)->getIndex();
609
#endif
610
611
612
613
      neighInfo = stack->traverseNeighbour3d(neighInfo, oppVertex);

      int edgeDof0, edgeDof1;
      for (edgeDof0 = 0; edgeDof0 < vertices; edgeDof0++)
614
	if (neigh->getDof(edgeDof0) == edge[0])
615
616
	  break;
      for (edgeDof1 = 0; edgeDof1 < vertices; edgeDof1++)
617
	if (neigh->getDof(edgeDof1) == edge[1])
618
619
620
621
622
	  break;

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

624
    (*elInfo) = neighInfo;
625
    
626
    return true;
627
628
629
  }


630
  ElInfo* RefinementManager3d::refineFunction(ElInfo* elInfo)
631
  {
632
    FUNCNAME_DBG("RefinementManager3d::refineFunction()");
Thomas Witkowski's avatar
Thomas Witkowski committed
633

634
    Element *el = elInfo->getElement();
635
636

    if (el->getMark() <= 0)  
637
638
      return elInfo;

639
    int bound = false;
640
    DegreeOfFreedom *edge[2];
641

642
643
    // === Get memory for a list of all elements at the refinement edge. ===

644
    RCNeighbourList refineList(mesh->getMaxEdgeNeigh());
645
    refineList.setElType(0, elInfo->getType());
646
    refineList.setElement(0, el);
647
648
    int n_neigh = 1;

649
650
    if (elInfo->getProjection(0) && 
	elInfo->getProjection(0)->getType() == VOLUME_PROJECTION)
651
      newCoords = true;
652
653


654
655
    // === Give the refinement edge the right orientation. ===

656
657
658
    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));
659
    } else {
660
661
      edge[1] = const_cast<DegreeOfFreedom*>(el->getDof(0));
      edge[0] = const_cast<DegreeOfFreedom*>(el->getDof(1));
662
    }
663

664

665
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
666
667
    vector<FixRefinementPatch::EdgeInEl> refineEdges;
    FixRefinementPatch::getOtherEl(stack, refineEdges);
668
669
670
671
672
673
//     if (refineEdges.size()) {
//       MSG("FIX REFINEMENT PATH ON ELEMENT %d %d: %d additional edges\n", 
// 	  elInfo->getElement()->getIndex(),
// 	  elInfo->getMacroElement()->getIndex(),
// 	  refineEdges.size());
//     }
674
#endif
675
676
677

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

678
    if (getRefinePatch(&elInfo, edge, 0, refineList, &n_neigh)) {      
679
      // Domain's boundary was reached while looping around the refinement edge
680
      getRefinePatch(&elInfo, edge, 1, refineList, &n_neigh);
681
682
      bound = true;
    }
683

684
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
685
686
687
688
    // === If the refinement edge must be fixed, add also the other part  ===
    // === of this edge to the refinement patch.                          ===

    for (int edgeIndex = 0; 
689
 	 edgeIndex < static_cast<unsigned int>(refineEdges.size()); edgeIndex++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
690
      Element *otherEl = refineEdges[edgeIndex].first;   
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
      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])) {

706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
	    // 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;
729
	      foundEdge = true;
730
731
732
	      break;
	    } else {
	      // Edge i not refinement edge, so refine the element further.
733

734
	      Element *el2 = elInfo2->getElement();
735
736
737
738
739
740
741
742
743
744
	      el2->setMark(std::max(el2->getMark(), 1));

	      TraverseStack *tmpStack = stack;
	      stack = &stack2;

	      elInfo2 = refineFunction(elInfo2);

	      stack = tmpStack;
	      break;
	    }
745
746
747
	  }
	}

748
749
750
	if (foundEdge)
	  break;

751
752
753
754
755
756
757
	elInfo2 = stack2.traverseNext(elInfo2);
      }

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

758
    // fill neighbour information inside the patch in the refinement list
759
    refineList.fillNeighbourRelations(n_neigh, bound);
760

761
    // ============ Check for periodic boundary ============
762

763
    DegreeOfFreedom *next_edge[2] = {NULL, NULL};
764
765
    DegreeOfFreedom *first_edge[2] = {edge[0], edge[1]};
    DegreeOfFreedom *last_edge[2] = {NULL, NULL};
766
    int n_neigh_periodic = 0;
767

768
769
    DegreeOfFreedom lastNewDof = -1;
    DegreeOfFreedom firstNewDof = -1;
770

771
    RCNeighbourList periodicList;
772

Thomas Witkowski's avatar
Thomas Witkowski committed
773
    while (edge[0] != NULL) {
774
775
776
      refineList.periodicSplit(edge, next_edge, 
			       &n_neigh, &n_neigh_periodic, 
			       periodicList);
777

778
779
      DegreeOfFreedom newDof = 
	refinePatch(edge, periodicList, n_neigh_periodic, bound);
780

781
782
      if (firstNewDof == -1)
	firstNewDof = newDof;
783

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

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

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

807
808
809
    if (lastNewDof != firstNewDof) {
      for (std::map<int, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
	   it != mesh->getPeriodicAssociations().end(); ++it) {
810
811
	if (it->second) {
	  if (((*(it->second))[first_edge[0][0]] == last_edge[0][0] &&
Thomas Witkowski's avatar
Thomas Witkowski committed
812
813
814
	       (*(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])) {	    
815
816
	    (*(it->second))[lastNewDof] = firstNewDof;
	    (*(it->second))[firstNewDof] = lastNewDof;
Thomas Witkowski's avatar
Thomas Witkowski committed
817
	  }
818
819
820
	}
      }
    }
821

822
823
    stack->update();

824
825
826
827
828
    return elInfo;
  }


  void FixRefinementPatch::getOtherEl(TraverseStack *stack, 
Thomas Witkowski's avatar
Thomas Witkowski committed
829
				      vector<EdgeInEl>& refineEdges)
830
  {
831
    FUNCNAME_DBG("FixRefinementPatch::getOtherEl()");
832
833
834

    if (!FixRefinementPatch::connectedEdges.empty()) {
      // === Get stack of current traverse. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
835
836
      vector<ElInfo*> elInfos;
      vector<int> infos;      
837
838
839
840
841
842
843
844
845
846
      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");
	
847

848
849
850
851
852
	int parentType = elInfos[checkIndex - 1]->getType();
	int ithParentChild = 0;
	if (elInfos[checkIndex - 1]->getElement()->getChild(1) ==
	    elInfos[checkIndex]->getElement())
	  ithParentChild = 1;
853
854
855
856