RefinementManager3d.cc 28.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
#include "RefinementManager.h"
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "DOFAdmin.h"
#include "AdaptStationary.h"
#include "AdaptInstationary.h"
#include "FixVec.h"
#include "RCNeighbourList.h"
#include "ProblemStatBase.h"
#include "DOFIndexed.h"
#include "Projection.h"
#include "DOFVector.h"
#include "PeriodicBC.h"
#include "VertexVector.h"
28
#include "Debug.h"
29
30
31

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
417
418
    int oppVertex = (*elInfo)->getOppVertex(localNeighbour);
    int testIndex = (*elInfo)->getNeighbour(localNeighbour)->getIndex();
    ElInfo *neighInfo = stack->traverseNeighbour3d((*elInfo), localNeighbour);
419
420
421
422
    int neighElType = neighInfo->getType();

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

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

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

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

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

460
461
462
463
	  TEST_EXIT_DBG(edgeDof1 < vertices)
 	    ("dof %d not found on element %d with nodes (%d %d %d %d)\n", 
 	     edge[1][0], neigh->getIndex(), neigh->getDof(0, 0),
 	     neigh->getDof(1, 0), neigh->getDof(2, 0), neigh->getDof(3, 0));
464
465
466
	}
      }

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

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

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

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

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

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

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

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

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

554
	refineList.setOppVertex(*n_neigh, 0, oppVertex); 
555
	
556
557
558
559
560
561
562
563
564
565
566
567
568
569
	(*n_neigh)++;

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

	if (neighInfo->getNeighbour(i)) {
	  oppVertex = neighInfo->getOppVertex(i);
	  int testIndex = neighInfo->getNeighbour(i)->getIndex();

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

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

	  neighElType = neighInfo->getType();
570
	  neigh = 
571
572
	    dynamic_cast<Tetrahedron*>(const_cast<Element*>(neighInfo->getElement()));
	} else {
573
	  break;
574
	}
575
      }
576
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
577
   
578
    if (neigh == el) {
579
      (*elInfo) = neighInfo;
580

581
      return false;
582
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
583
584

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

      int testIndex = neighInfo->getNeighbour(oppVertex)->getIndex();
      neighInfo = stack->traverseNeighbour3d(neighInfo, oppVertex);

      int edgeDof0, edgeDof1;
      for (edgeDof0 = 0; edgeDof0 < vertices; edgeDof0++)
599
	if (neigh->getDof(edgeDof0) == edge[0])
600
601
	  break;
      for (edgeDof1 = 0; edgeDof1 < vertices; edgeDof1++)
602
	if (neigh->getDof(edgeDof1) == edge[1])
603
604
605
606
607
	  break;

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

609
    (*elInfo) = neighInfo;
610
    
611
    return true;
612
613
614
  }


615
  ElInfo* RefinementManager3d::refineFunction(ElInfo* elInfo)
616
617
  {
    FUNCNAME("RefinementManager3d::refineFunction()");
Thomas Witkowski's avatar
Thomas Witkowski committed
618

619
    Element *el = elInfo->getElement();
620
621

    if (el->getMark() <= 0)  
622
623
624
      return elInfo;

    //    MSG("Refine: %d\n", elInfo->getElement()->getIndex());
625

626
    int bound = false;
627
    DegreeOfFreedom *edge[2];
628

629
630
    // === Get memory for a list of all elements at the refinement edge. ===

631
    RCNeighbourList refineList(mesh->getMaxEdgeNeigh());
632
    refineList.setElType(0, elInfo->getType());
633
    refineList.setElement(0, el);
634
635
    int n_neigh = 1;

636
637
    if (elInfo->getProjection(0) && 
	elInfo->getProjection(0)->getType() == VOLUME_PROJECTION)
638
      newCoords = true;
639
640


641
642
    // === Give the refinement edge the right orientation. ===

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

651
652
653
654
655
#if HAVE_PARALLEL_DOMAIN_AMDIS
    Element *otherEl = NULL;
    int otherEdge = -1;
    FixRefinementPatch::getOtherEl(stack, &otherEl, otherEdge);
#endif
656
657
658

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

659
    if (getRefinePatch(&elInfo, edge, 0, refineList, &n_neigh)) {      
660
      // Domain's boundary was reached while looping around the refinement edge
661
      getRefinePatch(&elInfo, edge, 1, refineList, &n_neigh);
662
663
      bound = true;
    }
664

665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
#if HAVE_PARALLEL_DOMAIN_AMDIS
    // === 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])) {
	    refineList.setElType(n_neigh, elInfo2->getType());
	    refineList.setElement(n_neigh, elInfo2->getElement());
	    n_neigh++;

	    foundEdge = true;
	    TraverseStack *tmpStack = stack;
	    stack = &stack2;
	    
	    if (getRefinePatch(&elInfo2, edge, 0, refineList, &n_neigh)) {
	      getRefinePatch(&elInfo2, edge, 1, refineList, &n_neigh);
	      bound = true;
	    }

	    stack = tmpStack;
	    break;
	  }
	}

	elInfo2 = stack2.traverseNext(elInfo2);
      }

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


710
    // fill neighbour information inside the patch in the refinement list
711
    refineList.fillNeighbourRelations(n_neigh, bound);
712

713
714

    // ============ Check for periodic boundary ============
715
716
717
718
719
720

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

721
722
    DegreeOfFreedom lastNewDof = -1;
    DegreeOfFreedom firstNewDof = -1;
723

724
    RCNeighbourList periodicList;
725

Thomas Witkowski's avatar
Thomas Witkowski committed
726
    while (edge[0] != NULL) {
727
728
729
      refineList.periodicSplit(edge, next_edge, 
			       &n_neigh, &n_neigh_periodic, 
			       periodicList);
730

731
      DegreeOfFreedom newDof = refinePatch(edge, periodicList, n_neigh_periodic, bound);
732

733
734
      if (firstNewDof == -1)
	firstNewDof = newDof;
735

736
737
738
      if (lastNewDof != -1) {
	for (std::map<int, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
	     it != mesh->getPeriodicAssociations().end(); ++it) {
739
740
	  if (it->second) {
	    if (((*(it->second))[edge[0][0]] == last_edge[0][0] &&
Thomas Witkowski's avatar
Thomas Witkowski committed
741
742
743
		 (*(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])) {
744
745
	      (*(it->second))[lastNewDof] = newDof;
	      (*(it->second))[newDof] = lastNewDof;
Thomas Witkowski's avatar
Thomas Witkowski committed
746
	    } 
747
748
749
	  }
	}
      }
750
      lastNewDof = newDof;
751
752
753
754
755
756
757
758

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

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

759
760
761
    if (lastNewDof != firstNewDof) {
      for (std::map<int, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
	   it != mesh->getPeriodicAssociations().end(); ++it) {
762
763
	if (it->second) {
	  if (((*(it->second))[first_edge[0][0]] == last_edge[0][0] &&
Thomas Witkowski's avatar
Thomas Witkowski committed
764
765
766
	       (*(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])) {	    
767
768
	    (*(it->second))[lastNewDof] = firstNewDof;
	    (*(it->second))[firstNewDof] = lastNewDof;
Thomas Witkowski's avatar
Thomas Witkowski committed
769
	  }
770
771
772
773
774
775
	}
      }
    }
  
    stack->update();

776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
    return elInfo;
  }


  void FixRefinementPatch::getOtherEl(TraverseStack *stack, 
				     Element **otherEl, 
				     int &otherEdge)
  {
    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;
	  }
	}
      }
    }
837
838
  }
}