RefinementManager3d.cc 30.2 KB
Newer Older
</
1
2
3
4
5
6
7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9
10
11
12
13
14
15
16
17
 * 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.
18
 *
19
 ******************************************************************************/
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
using namespace std;

41
42
namespace AMDiS {

43
  map<Mesh*, FixRefinementPatch::ConnectedEdges> FixRefinementPatch::connectedEdges;
44

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

    child[0] = dynamic_cast<Tetrahedron*>(mesh->createNewElement(el));
    child[1] = dynamic_cast<Tetrahedron*>(mesh->createNewElement(el));
57

58
    int mark = std::max(0, el->getMark() - 1);
59
60
61
62
63
64
65
66
67
68
69
70
71
    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]);

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

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

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

97
    if (mesh->getNumberOfDofs(EDGE)) {
98
      int node = mesh->getNode(EDGE);
99

100
101
102
      /****************************************************************************/
      /*  set pointers to those dof's that are handed on from the parant          */
      /****************************************************************************/
103

104
      child[0]->
105
	setDof(node,
106
	       const_cast<DegreeOfFreedom*>(el->getDof(node + Tetrahedron::childEdge[el_type][0][0])));
107
      child[1]->
108
	setDof(node,
109
	       const_cast<DegreeOfFreedom*>(el->getDof(node + Tetrahedron::childEdge[el_type][1][0])));
110
      child[0]->
111
	setDof(node + 1,
112
	       const_cast<DegreeOfFreedom*>(el->getDof(node + Tetrahedron::childEdge[el_type][0][1])));
113
      child[1]->
114
	setDof(node + 1,
115
	       const_cast<DegreeOfFreedom*>(el->getDof(node + Tetrahedron::childEdge[el_type][1][1])));
116
      child[0]->
117
	setDof(node + 3,
118
	       const_cast<DegreeOfFreedom*>(el->getDof(node + Tetrahedron::childEdge[el_type][0][3])));
119
      child[1]->
120
	setDof(node + 3,
121
	       const_cast<DegreeOfFreedom*>(el->getDof(node + Tetrahedron::childEdge[el_type][1][3])));
122

123
124
125
      /****************************************************************************/
      /*  adjust pointers to the dof's in the refinement edge                     */
      /****************************************************************************/
126

127
128
129
      if (el->getDof(0) == edge[0]) {
	child[0]->setDof(node + 2, dof[1]);
	child[1]->setDof(node + 2, dof[2]);
130
      } else {
131
132
	child[0]->setDof(node + 2, dof[2]);
	child[1]->setDof(node + 2, dof[1]);
133
134
      }
    }
135

136
    if (mesh->getNumberOfDofs(FACE)) {
137
138
      int node = mesh->getNode(FACE);

139
140
141
      /****************************************************************************/
      /*  set pointers to those dof's that are handed on from the parant          */
      /****************************************************************************/
142

143
144
      child[0]->setDof(node + 3, const_cast<DegreeOfFreedom*>(el->getDof(node + 1)));
      child[1]->setDof(node + 3, const_cast<DegreeOfFreedom*>(el->getDof(node + 0)));
145

146
147
148
      /****************************************************************************/
      /*  get new dof for the common face of child0 and child1                    */
      /****************************************************************************/
149

150
      DegreeOfFreedom *newDOF = mesh->getDof(FACE);
151
152
      child[0]->setDof(node, static_cast<DegreeOfFreedom*>(newDOF));
      child[1]->setDof(node, static_cast<DegreeOfFreedom*>(newDOF));
153
    }
154

155
    if (mesh->getNumberOfDofs(CENTER)) {
156
      int node = mesh->getNode(CENTER);
157
158
      child[0]->setDof(node, const_cast<DegreeOfFreedom*>(mesh->getDof(CENTER)));
      child[1]->setDof(node, const_cast<DegreeOfFreedom*>(mesh->getDof(CENTER)));
159
    }
160

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

165

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

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

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

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

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

211
	if (el->getDof(0) == neigh->getDof(0)) {
212
	  // Same position at refinement edge.
213
214
	  adjc = 0;
	} else {
215
	  // Different position at refinement edge.
216
217
	  adjc = 1;
	}
218

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

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

225
226
227
228
	  /****************************************************************************/
	  /*  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!)   */
	  /****************************************************************************/
229

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

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

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

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

249
	    (const_cast<Element*>(el->getChild(i)))->
250
	      setDof(node0, const_cast<DegreeOfFreedom*>(neigh->getChild(j)->getDof(node1)));
251
252
	  }

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


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

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

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

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

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

      projector->project(*new_coord);

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

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

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

295
	getRefinePatch(&elinfo, edge, 1, refineList, &n_neigh);
296
297
      }

298
      for (int i = 1; i < n_neigh; i++) {            /* start with 1, as list[0]=el */
299
300
301
	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);
302

303
	refineList.getElement(i)->setNewCoord(el->getNewCoord());
304
      }
305
306
307
    }
  }

308

309
  void RefinementManager3d::setNewCoords(int macroEl)
310
  {
311
    RCNeighbourList refineList(mesh->getMaxEdgeNeigh());
312
313
314
    ElInfo *elInfo;

    if (macroEl == -1)
315
316
317
      elInfo = stack->traverseFirst(mesh, -1,
				    Mesh::CALL_EVERY_EL_PREORDER |
				    Mesh::FILL_BOUND | Mesh::FILL_COORDS |
318
319
320
				    Mesh::FILL_NEIGH);
    else
      elInfo = stack->traverseFirstOneMacro(mesh, macroEl, -1,
321
322
					    Mesh::CALL_EVERY_EL_PREORDER |
					    Mesh::FILL_BOUND | Mesh::FILL_COORDS |
323
					    Mesh::FILL_NEIGH);
324

325

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


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

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

346
    dof[0] = mesh->getDof(VERTEX);
347
    mesh->incrementNumberOfVertices(1);
348

349
    if (mesh->getNumberOfDofs(EDGE)) {
350
351
      dof[1] = mesh->getDof(EDGE);
      dof[2] = mesh->getDof(EDGE);
352
    }
353

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

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

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

378
	el = dynamic_cast<Tetrahedron*>(const_cast<Element*>(refineList.getElement(0)));
379
	mesh->freeDof(const_cast<DegreeOfFreedom*>(el->getDof(mesh->getNode(EDGE))), EDGE);
380
      }
381
382

      if (mesh->getNumberOfDofs(EDGE) ||
383
384
	  mesh->getNumberOfDofs(FACE) ||
	  mesh->getNumberOfDofs(CENTER)) {
385
	for (int i = 0; i < n_neigh; i++)
386
	  refineList.removeDOFParent(i);
387
388
      }
    }
389
390
391
392
393
394
395


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

396
397
    if (bound) {
      mesh->incrementNumberOfEdges(n_neigh + 2);
398
      mesh->incrementNumberOfFaces(2 * n_neigh + 1);
399
      newCoords = true; // added to allow BOUNDARY_PROJECTION
400
401
    } else {
      mesh->incrementNumberOfEdges(n_neigh + 1);
402
      mesh->incrementNumberOfFaces(2 * n_neigh);
403
    }
404

405
406
407
408
    return dof[0][0];
  }


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

417
    int localNeighbour = 3 - direction;
418
    Tetrahedron *el =
419
      dynamic_cast<Tetrahedron*>(const_cast<Element*>((*elInfo)->getElement()));
420

421
    if ((*elInfo)->getNeighbour(localNeighbour) == NULL)
422
423
      return true;

424
    int oppVertex = (*elInfo)->getOppVertex(localNeighbour);
425
#ifndef NDEBUG
426
    int testIndex = (*elInfo)->getNeighbour(localNeighbour)->getIndex();
427
#endif
428
    ElInfo *neighInfo = stack->traverseNeighbour3d((*elInfo), localNeighbour);
429
430
431
432
    int neighElType = neighInfo->getType();

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

434
435
    Tetrahedron *neigh =
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(neighInfo->getElement()));
436
437
438
    int vertices = mesh->getGeo(VERTEX);

    while (neigh != el) {
439
440
441
      // === Determine the common edge of the element and its neighbour. ===

      int edgeDof0, edgeDof1;
442
443
444

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

452
453
      if (edgeDof0 > 3 || edgeDof1 > 3) {
	for (edgeDof0 = 0; edgeDof0 < vertices; edgeDof0++)
454
	  if (mesh->associated(neigh->getDof(edgeDof0, 0), edge[0][0]))
455
	    break;
456
	for (edgeDof1 = 0; edgeDof1 < vertices; edgeDof1++)
457
	  if (mesh->associated(neigh->getDof(edgeDof1, 0), edge[1][0]))
458
	    break;
459

460
461
	if (edgeDof0 > 3 || edgeDof1 > 3) {
	  for (edgeDof0 = 0; edgeDof0 < vertices; edgeDof0++)
462
	    if (mesh->indirectlyAssociated(neigh->getDof(edgeDof0, 0), edge[0][0]))
463
	      break;
464
	  for (edgeDof1 = 0; edgeDof1 < vertices; edgeDof1++)
465
	    if (mesh->indirectlyAssociated(neigh->getDof(edgeDof1, 0), edge[1][0]))
466
	      break;
467

468
	  TEST_EXIT_DBG(edgeDof0 < vertices)
469
 	    ("DOF %d not found on element %d with nodes (%d %d %d %d)\n",
470
 	     edge[0][0], neigh->getIndex(), neigh->getDof(0, 0),
471
 	     neigh->getDof(1, 0), neigh->getDof(2, 0), neigh->getDof(3, 0));
472

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

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

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

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

491
	neigh->setMark(std::max(neigh->getMark(), 1));
492
493
494
495
496
497

	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.                                                         ===
498

499
500
501
	neighInfo = stack->traverseNext(neighInfo);
	neighElType = neighInfo->getType();
	bool reverseMode = stack->getTraverseFlag().isSet(Mesh::CALL_REVERSE_MODE);
502

503
	switch (edgeNo) {
504
	case 1:
505
506
507
508
	  if (reverseMode) {
	    neighInfo = stack->traverseNext(neighInfo);
	    neighElType = neighInfo->getType();
	  }
509

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

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

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

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

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

559
560
	TEST_EXIT_DBG(*n_neigh < mesh->getMaxEdgeNeigh())
	  ("too many neighbours %d in refine patch\n", *n_neigh);
561

562
563
	refineList.setElement(*n_neigh, neigh);
	refineList.setElType(*n_neigh, neighElType);
564

565
	// We have to go back to the starting element via oppVertex values.
566

567
568
	refineList.setOppVertex(*n_neigh, 0, oppVertex);

569
570
571
572
573
574
	(*n_neigh)++;

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

	if (neighInfo->getNeighbour(i)) {
	  oppVertex = neighInfo->getOppVertex(i);
575

576
#ifndef NDEBUG
577
	  int testIndex = neighInfo->getNeighbour(i)->getIndex();
578
#endif
579
580
581
582
583
584
585

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

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

	  neighElType = neighInfo->getType();
586
	  neigh =
587
588
	    dynamic_cast<Tetrahedron*>(const_cast<Element*>(neighInfo->getElement()));
	} else {
589
	  break;
590
	}
591
      }
592
    }
593

594
    if (neigh == el) {
595
      (*elInfo) = neighInfo;
596

597
      return false;
598
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
599

600

601
    // The domain's boundary is reached. Loop back to the starting el.
602

603
    int i = *n_neigh - 1;
604
    oppVertex = refineList.getOppVertex(i, 0);
605
    do {
606
      TEST_EXIT_DBG(neighInfo->getNeighbour(oppVertex) && i > 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
607
	("While looping back domains boundary was reached or i == 0\n");
608
      oppVertex = refineList.getOppVertex(i--, 0);
609

610
#ifndef NDEBUG
611
      int testIndex = neighInfo->getNeighbour(oppVertex)->getIndex();
612
#endif
613
614
615
616
      neighInfo = stack->traverseNeighbour3d(neighInfo, oppVertex);

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

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

627
    (*elInfo) = neighInfo;
628

629
    return true;
630
631
632
  }


633
  ElInfo* RefinementManager3d::refineFunction(ElInfo* elInfo)
634
  {
635
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
636
    FUNCNAME_DBG("RefinementManager3d::refineFunction()");
637
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
638

639
    Element *el = elInfo->getElement();
640

641
    if (el->getMark() <= 0)
642
643
      return elInfo;

644
    int bound = false;
645
    DegreeOfFreedom *edge[2];
646

647
648
    // === Get memory for a list of all elements at the refinement edge. ===

649
    RCNeighbourList refineList(mesh->getMaxEdgeNeigh());
650
    refineList.setElType(0, elInfo->getType());
651
    refineList.setElement(0, el);
652
653
    int n_neigh = 1;

654
    if (elInfo->getProjection(0) &&
655
	elInfo->getProjection(0)->getType() == VOLUME_PROJECTION)
656
      newCoords = true;
657
658


659
660
    // === Give the refinement edge the right orientation. ===

661
662
663
    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));
664
    } else {
665
666
      edge[1] = const_cast<DegreeOfFreedom*>(el->getDof(0));
      edge[0] = const_cast<DegreeOfFreedom*>(el->getDof(1));
667
    }
668

669

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

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

683
    if (getRefinePatch(&elInfo, edge, 0, refineList, &n_neigh)) {
684
      // Domain's boundary was reached while looping around the refinement edge
685
      getRefinePatch(&elInfo, edge, 1, refineList, &n_neigh);
686
687
      bound = true;
    }
688

689
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
690
691
692
    // === If the refinement edge must be fixed, add also the other part  ===
    // === of this edge to the refinement patch.                          ===

693
    for (int edgeIndex = 0;
694
 	 edgeIndex < static_cast<int>(refineEdges.size()); edgeIndex++) {
695
      Element *otherEl = refineEdges[edgeIndex].first;