RefinementManager3d.cc 25.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#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"

namespace AMDiS {

  void RefinementManager3d::bisectTetrahedron(RCNeighbourList* ref_list, 
					      int index,
					      DegreeOfFreedom* dof[3],
					      DegreeOfFreedom *edge[2])
  {
24
25
26
    Tetrahedron *el = 
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(ref_list->getElement(index)));
    Tetrahedron *child[2];
27
    int el_type = ref_list->getType(index);
28
29
30
31

    child[0] = dynamic_cast<Tetrahedron*>(mesh->createNewElement(el));
    child[1] = dynamic_cast<Tetrahedron*>(mesh->createNewElement(el));
  
32
    int mark = max(0, el->getMark() - 1);
33
34
35
36
37
38
39
40
41
42
43
44
45
    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]);

46
47
    if (child[0]->getMark() > 0)
      doMoreRecursiveRefine = true;
48
49

    int n_vertices = mesh->getGeo(VERTEX);
50
51
    child[0]->setDOF(n_vertices - 1, dof[0]);
    child[1]->setDOF(n_vertices - 1, dof[0]);
52

53
    for (int i = 0; i < n_vertices - 1; i++) {
54
      child[0]->
55
	setDOF(i, const_cast<int*>(el->getDOF(Tetrahedron::childVertex[el_type][0][i])));
56
      child[1]->
57
	setDOF(i, const_cast<int*>(el->getDOF(Tetrahedron::childVertex[el_type][1][i])));
58
    }
59
60
61
62
63
64
65
66
67
68
69
70
    /****************************************************************************/
    /*  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                                                              */
    /****************************************************************************/

71
    if (mesh->getNumberOfDOFs(EDGE)) {
72
      int node = mesh->getNode(EDGE);
73

74
75
76
77
78
79
      /****************************************************************************/
      /*  set pointers to those dof's that are handed on from the parant          */
      /****************************************************************************/
      
      child[0]->
	setDOF(node, 
80
	       const_cast<int*>(el->getDOF(node + Tetrahedron::childEdge[el_type][0][0])));
81
82
      child[1]->
	setDOF(node, 
83
	       const_cast<int*>(el->getDOF(node + Tetrahedron::childEdge[el_type][1][0])));
84
      child[0]->
85
	setDOF(node + 1, 
86
	       const_cast<int*>(el->getDOF(node + Tetrahedron::childEdge[el_type][0][1])));
87
      child[1]->
88
	setDOF(node + 1, 
89
	       const_cast<int*>(el->getDOF(node + Tetrahedron::childEdge[el_type][1][1])));
90
      child[0]->
91
	setDOF(node + 3, 
92
	       const_cast<int*>(el->getDOF(node + Tetrahedron::childEdge[el_type][0][3])));
93
      child[1]->
94
	setDOF(node + 3, 
95
	       const_cast<int*>(el->getDOF(node + Tetrahedron::childEdge[el_type][1][3])));
96
97
98
99
100
101
      
      /****************************************************************************/
      /*  adjust pointers to the dof's in the refinement edge                     */
      /****************************************************************************/
      
      if (el->getDOF(0) == edge[0]) {
102
103
	child[0]->setDOF(node + 2, dof[1]);
	child[1]->setDOF(node + 2, dof[2]);
104
      } else {
105
106
	child[0]->setDOF(node + 2, dof[2]);
	child[1]->setDOF(node + 2, dof[1]);
107
108
109
110
      }
    }
    
    if (mesh->getNumberOfDOFs(FACE)) {
111
      int node = mesh->getNode(FACE);    
112
113
114
115
116
      
      /****************************************************************************/
      /*  set pointers to those dof's that are handed on from the parant          */
      /****************************************************************************/
      
117
118
      child[0]->setDOF(node + 3, const_cast<int*>(el->getDOF(node + 1)));
      child[1]->setDOF(node + 3, const_cast<int*>(el->getDOF(node + 0)));
119
120
121
122
123
124
      
      /****************************************************************************/
      /*  get new dof for the common face of child0 and child1                    */
      /****************************************************************************/
      
      DegreeOfFreedom *newDOF = mesh->getDOF(FACE);
125
126
      child[0]->setDOF(node, static_cast<int*>(newDOF));
      child[1]->setDOF(node, static_cast<int*>(newDOF));
127
128
129
    }
    
    if (mesh->getNumberOfDOFs(CENTER)) {
130
131
132
      int node = mesh->getNode(CENTER);
      child[0]->setDOF(node, const_cast<int*>(mesh->getDOF(CENTER)));
      child[1]->setDOF(node, const_cast<int*>(mesh->getDOF(CENTER)));
133
    }
134

135
    if (mesh->getNumberOfDOFs(EDGE) || mesh->getNumberOfDOFs(FACE))
136
137
      fillPatchConnectivity(ref_list, index);
  }
138

139

140
141
142
143
  void RefinementManager3d::fillPatchConnectivity(RCNeighbourList* ref_list, 
						  int index)
  {
    FUNCNAME("RefinementManager3d::fillPatchConnectivity");
Thomas Witkowski's avatar
Thomas Witkowski committed
144

145
    Element *el = ref_list->getElement(index), *neigh;
146
147
148
149
    int el_type = ref_list->getType(index);
    int n_type = 0;
    int dir, adjc, i, j, i_neigh, j_neigh;
    int node0, node1, opp_v = 0;
150

151
152
153
    for (dir = 0; dir < 2; dir++) {
      neigh = ref_list->getNeighbourElement(index, dir);
      if (neigh) {
154
	n_type = ref_list->getType(ref_list->getNeighbourNr(index, dir));
155
	opp_v = ref_list->getOppVertex(index, dir);
156
157
      }

158
      if (!neigh || neigh->isLeaf()) {
159
160

	/****************************************************************************/
161
162
163
	/*  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)                      */
164
165
	/****************************************************************************/

166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
	if (mesh->getNumberOfDOFs(EDGE)) {
	  node0 = node1 = mesh->getNode(EDGE);
	  node0 += Tetrahedron::nChildEdge[el_type][0][dir];
	  node1 += Tetrahedron::nChildEdge[el_type][1][dir];
	  DegreeOfFreedom *newDOF = mesh->getDOF(EDGE);
	  (const_cast<Element*>(el->getFirstChild()))->setDOF(node0, newDOF);
	  (const_cast<Element*>(el->getSecondChild()))->setDOF(node1, newDOF);
	}
	if (mesh->getNumberOfDOFs(FACE)) {
	  node0 = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][0][dir];
	  (const_cast<Element*>(el->getFirstChild()))->setDOF(node0, mesh->getDOF(FACE));
	  node1 = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][1][dir];
	  (const_cast<Element*>(el->getSecondChild()))->setDOF(node1, mesh->getDOF(FACE));
	}
      } else {     /*   if (!neigh  ||  !neigh->child[0])                         */
181
	/****************************************************************************/
182
183
	/*  interior face and neighbour has been refined, look for position at the  */
	/*  refinement edge                                                         */
184
	/****************************************************************************/
185
186
187
188
189
190
191
192
193
194
195
196
      
	if (el->getDOF(0) == neigh->getDOF(0)) {
	  /****************************************************************************/
	  /* same position at refinement edge                                         */
	  /****************************************************************************/
	  adjc = 0;
	} else {
	  /****************************************************************************/
	  /* different position at refinement edge                                    */
	  /****************************************************************************/
	  adjc = 1;
	}
197

198
199
	for (i = 0; i < 2; i++) {
	  j = Tetrahedron::adjacentChild[adjc][i];
200

201
	  i_neigh = Tetrahedron::nChildFace[el_type][i][dir];
202
	  j_neigh = Tetrahedron::nChildFace[n_type][j][opp_v - 2];
203

204
205
206
207
	  /****************************************************************************/
	  /*  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!)   */
	  /****************************************************************************/
208

209
	  if (mesh->getNumberOfDOFs(EDGE)) {
210
211
	    node0 = mesh->getNode(EDGE) + Tetrahedron::nChildEdge[el_type][i][dir];
	    node1 = mesh->getNode(EDGE) + Tetrahedron::nChildEdge[n_type][j][opp_v - 2];
212

213
214
215
	    TEST_EXIT_DBG(neigh->getChild(j)->getDOF(node1))
	      ("no dof on neighbour %d at node %d\n", 
	       neigh->getChild(j)->getIndex(), node1);
216

217
	    (const_cast<Element*>(el->getChild(i)))->
218
	      setDOF(node0, const_cast<int*>(neigh->getChild(j)->getDOF(node1)));
219
220
221
222
	  }
	  if (mesh->getNumberOfDOFs(FACE)) {
	    node0 = mesh->getNode(FACE) + i_neigh;
	    node1 = mesh->getNode(FACE) + j_neigh;
223

224
	    TEST_EXIT_DBG(neigh->getChild(j)->getDOF(node1))
Thomas Witkowski's avatar
Thomas Witkowski committed
225
	      ("No DOF on neighbour %d at node %d!\n",
226
	       neigh->getChild(j)->getIndex(), node1);
227

228
229
	    (const_cast<Element*>(el->getChild(i)))->
	      setDOF(node0, const_cast<int*>(neigh->getChild(j)->getDOF(node1)));
230
231
	  }

232
233
234
	}  /*   for (i = 0; i < 2; i++)                                       */
      }    /*   else of   if (!neigh  ||  !neigh->child[0])                   */
    }      /*   for (dir = 0; dir < 2; dir++)                                 */
235
236
237
238
239
240
241
  }


  int RefinementManager3d::newCoordsFct(ElInfo *el_info)
  {
    Element *el = el_info->getElement();
    DegreeOfFreedom *edge[2];
242
    ElInfo *elinfo = el_info;
243
244
245
    int dow = Global::getGeo(WORLD);
    Projection *projector = el_info->getProjection(0);

246
247
    if (!projector || projector->getType() != VOLUME_PROJECTION)
      projector = el_info->getProjection(4);    
248
249

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

252
253
      for (int j = 0; j < dow; j++)
	(*new_coord)[j] = (el_info->getCoord(0)[j] + el_info->getCoord(1)[j]) * 0.5;
254
255
256
257
258
259
260
261
262

      projector->project(*new_coord);

      el->setNewCoord(new_coord);
      /****************************************************************************/
      /*  now, information should be passed on to patch neighbours...             */
      /*  get the refinement patch                                                */
      /****************************************************************************/
      refList->setElement(0, el);
263
264
      refList->setElType(0, el_info->getType());
      int n_neigh = 1;
265

266
267
      for (int i = 0; i < 2; i++)
	edge[i] = const_cast<int*>(el_info->getElement()->getDOF(i));
268

269
270
271
272
273
274
275
      if (getRefinePatch(&elinfo, edge, 0, refList, &n_neigh)) {
	/****************************************************************************/
	/*  domain's boundary was reached while looping around the refinement edge  */
	/****************************************************************************/
	getRefinePatch(&elinfo, edge, 1, refList, &n_neigh);
      }

276
      for (int i = 1; i < n_neigh; i++) {            /* start with 1, as list[0]=el */
277
278
279
280
281
282
	TEST(!refList->getElement(i)->isNewCoordSet())
	  ("non-nil new_coord in el %d ref_list[%d] el %d (n_neigh=%d)\n",
	   el->getIndex(), i, refList->getElement(i)->getIndex(), n_neigh);
	
	refList->getElement(i)->setNewCoord(el->getNewCoord());
      }
283
284
285
286
287
    }

    return 0;
  }

288

289
290
291
  void RefinementManager3d::setNewCoords()
  {
    if (refList)
Thomas Witkowski's avatar
Thomas Witkowski committed
292
      delete refList;
293
 
Thomas Witkowski's avatar
Thomas Witkowski committed
294
    refList = new RCNeighbourList(mesh->getMaxEdgeNeigh());
295

296
297
298
299
300
301
302
    ElInfo *elInfo = stack->traverseFirst(mesh, -1, 
					  Mesh::CALL_EVERY_EL_PREORDER | 
					  Mesh::FILL_BOUND | Mesh::FILL_COORDS | 
					  Mesh::FILL_NEIGH);
    while (elInfo) {
      newCoordsFct(elInfo);
      elInfo = stack->traverseNext(elInfo);
303
304
305
306
307
308
309
310
    }
  }


  DegreeOfFreedom RefinementManager3d::refinePatch(DegreeOfFreedom *edge[2], 
						   RCNeighbourList* refineList,
						   int n_neigh, bool bound)
  {
311
312
    Tetrahedron *el = 
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(refineList->getElement(0)));
313
314
315
316
317
318
319
320
321
322
    /* first element in the list */
    DegreeOfFreedom *dof[3] = {NULL, NULL, NULL};

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

    dof[0] = mesh->getDOF(VERTEX);
    mesh->incrementNumberOfVertices(1);
  
323
324
325
326
    if (mesh->getNumberOfDOFs(EDGE)) {
      dof[1] = mesh->getDOF(EDGE);
      dof[2] = mesh->getDOF(EDGE);
    }
327

328
    for (int i = 0; i < n_neigh; i++)
329
330
331
332
333
334
335
      bisectTetrahedron(refineList, i, dof, edge);

    /****************************************************************************/
    /*  if there are functions to interpolate data to the finer grid, do so     */
    /****************************************************************************/
    int iadmin;
    int nrAdmin = mesh->getNumberOfDOFAdmin();
336
    for (iadmin = 0; iadmin < nrAdmin; iadmin++) {
337
      std::list<DOFIndexedBase*>::iterator it;
338
      DOFAdmin* admin = const_cast<DOFAdmin*>(&mesh->getDOFAdmin(iadmin));
339
      std::list<DOFIndexedBase*>::iterator end = admin->endDOFIndexed();
340
      for (it = admin->beginDOFIndexed(); it != end; it++)
341
342
343
	(*it)->refineInterpol(*refineList, n_neigh);
    }

344
345
346
347
348
349
    if (!mesh->queryCoarseDOFs()) {
      /****************************************************************************/
      /*  if there should be no dof information on interior leaf elements remove  */
      /*  dofs from edges, faces and the centers of parents                       */
      /****************************************************************************/
      if (mesh->getNumberOfDOFs(EDGE)) {
350
	/****************************************************************************/
351
	/*  remove dof of the midpoint of the common refinement edge                */
352
353
	/****************************************************************************/

354
355
	el = dynamic_cast<Tetrahedron*>(const_cast<Element*>(refineList->getElement(0)));
	mesh->freeDOF(const_cast<int*>(el->getDOF(mesh->getNode(EDGE))), EDGE);
356
      }
357
      
358
359
      if (mesh->getNumberOfDOFs(EDGE) || 
	  mesh->getNumberOfDOFs(FACE) ||
360
	  mesh->getNumberOfDOFs(CENTER)) {
361
	for (int i = 0; i < n_neigh; i++)
362
363
364
	  refineList->removeDOFParent(i);
      }
    }
365
366
367
368
369
370
371


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

372
373
    if (bound) {
      mesh->incrementNumberOfEdges(n_neigh + 2);
374
      mesh->incrementNumberOfFaces(2 * n_neigh + 1);
375
376
    } else {
      mesh->incrementNumberOfEdges(n_neigh + 1);
377
      mesh->incrementNumberOfFaces(2 * n_neigh);
378
379
    }
    
380
381
382
383
384
385
386
387
388
389
    return dof[0][0];
  }


  int RefinementManager3d::getRefinePatch(ElInfo **el_info, 
					  DegreeOfFreedom *edge[2], 
					  int dir,
					  RCNeighbourList* refineList, 
					  int *n_neigh)
  {
390
    FUNCNAME("RefinementManager3d::getRefinePatch()");
391

392
    int i, j, k;
393

394
395
396
397
    Tetrahedron *el = 
      dynamic_cast<Tetrahedron*>(const_cast<Element*>((*el_info)->getElement()));
    
    if ((*el_info)->getNeighbour(3 - dir) == NULL)
Thomas Witkowski's avatar
Thomas Witkowski committed
398
      return 1;    
399
  
400
401
402
    int opp_v = (*el_info)->getOppVertex(3-dir);
    ElInfo *neigh_info = stack->traverseNeighbour3d((*el_info), 3 - dir);
    int neigh_el_type = neigh_info->getType();
403

404
    Tetrahedron *neigh = 
405
      dynamic_cast<Tetrahedron*>(const_cast<Element*>(neigh_info->getElement()));
406
407
408
409
410
  
    int vertices = mesh->getGeo(VERTEX);

    while (neigh != el) {
      for (j = 0; j < vertices; j++)
411
412
	if (neigh->getDOF(j) == edge[0])
	  break;
413
      for (k = 0; k < vertices; k++)
414
415
	if (neigh->getDOF(k) == edge[1])
	  break;
416
417


418
      if (j > 3 || k > 3) {
419
	for (j = 0; j < vertices; j++)
420
421
	  if (mesh->associated(neigh->getDOF(j, 0), edge[0][0]))  
	    break;
422
	for (k = 0; k < vertices; k++)
423
424
	  if (mesh->associated(neigh->getDOF(k, 0), edge[1][0]))  
	    break;
425
	    
426
	if (j > 3 || k > 3) {
427
	  for (j = 0; j < vertices; j++)
428
429
	    if (mesh->indirectlyAssociated(neigh->getDOF(j, 0), edge[0][0]))  
	      break;
430
	  for (k = 0; k < vertices; k++)
431
432
	    if (mesh->indirectlyAssociated(neigh->getDOF(k, 0), edge[1][0]))  
	      break;
433
	    
434
	  TEST_EXIT_DBG(j < vertices && k < vertices)
435
436
437
438
439
440
441
	    ("dof %d or dof %d not found on element %d with nodes (%d %d %d %d)\n", 
	     edge[0][0], edge[1][0], neigh->getIndex(), neigh->getDOF(0,0),
	     neigh->getDOF(1,0), neigh->getDOF(2,0), neigh->getDOF(3,0));

	}
      }

442
      TEST_EXIT_DBG(j <  mesh->getGeo(VERTEX) &&
443
444
445
446
447
		k < mesh->getGeo(VERTEX))
	("dof %d or dof %d not found on element %d with nodes (%d %d %d %d)\n", 
	 edge[0][0], edge[1][0], neigh->getIndex(), neigh->getDOF(0,0),
	 neigh->getDOF(1,0), neigh->getDOF(2,0), neigh->getDOF(3,0));

448
449
      int edge_no = Tetrahedron::edgeOfDOFs[j][k];
      if (edge_no) {
450
451
452
453
454
455
456
457
458
459
460
461
	/****************************************************************************/
	/*  neigh has not a compatible refinement edge                              */
	/****************************************************************************/
	neigh->setMark(max(neigh->getMark(), 1));
	neigh_info = refineFunction(neigh_info);
	
	/****************************************************************************/
	/*  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      */
	/****************************************************************************/
	
	neigh_info = stack->traverseNext(neigh_info);
462
	neigh_el_type = neigh_info->getType();
463
464
465
466
467
468
469
470
471
	
	switch (edge_no) {
	case 1: 
	  opp_v = opp_v == 1 ? 3 : 2;
	  break;
	case 2: 
	  opp_v = opp_v == 2 ? 1 : 3;
	  break;
	case 3: 
472
	  neigh_info = stack->traverseNext(neigh_info);
473
	  neigh_el_type = neigh_info->getType();
474
475
	  if (neigh_el_type != 1)
	    opp_v = opp_v == 0 ? 3 : 2;
476
	  else
477
478
479
480
	    opp_v = opp_v == 0 ? 3 : 1;
	  break;
	case 4:
	  neigh_info = stack->traverseNext(neigh_info);
481
	  neigh_el_type = neigh_info->getType();
482
483
	  if (neigh_el_type != 1)
	    opp_v = opp_v == 0 ? 3 : 1;
484
	  else
485
486
487
488
489
	    opp_v = opp_v == 0 ? 3 : 2;
	  break;
	case 5:
	  if (neigh_el_type != 1) {
	    neigh_info = stack->traverseNext(neigh_info);
490
	    neigh_el_type = neigh_info->getType();
491
492
493
	  }
	  opp_v = 3;
	  break;
494
	}
495
496
	neigh = 
	  dynamic_cast<Tetrahedron*>(const_cast<Element*>(neigh_info->getElement()));
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
      } else {
	/****************************************************************************/
	/*  neigh is compatible devisible; put neigh to the list of patch elements; */
	/*  go to next neighbour                                                    */
	/****************************************************************************/
	TEST_EXIT_DBG(*n_neigh < mesh->getMaxEdgeNeigh())
	  ("too many neighbours %d in refine patch\n", *n_neigh);
	
	
	refineList->setElement(*n_neigh, neigh);
	refineList->setElType(*n_neigh, neigh_el_type);
	
	/****************************************************************************/
	/*  we have to go back to the starting element via opp_v values             */
	/*  correct information is produce by get_neigh_on_patch()                  */
	/****************************************************************************/
	refineList->setOppVertex(*n_neigh, 0, opp_v); 
	
	++*n_neigh;
516

517
518
519
520
521
522
523
524
	if (opp_v != 3)
	  i = 3;
	else
	  i = 2;
	
	if (neigh_info->getNeighbour(i)) {
	  opp_v = neigh_info->getOppVertex(i);
	  neigh_info = stack->traverseNeighbour3d(neigh_info, i);
525
526
527
	  neigh_el_type = neigh_info->getType();
	  neigh = 
	    dynamic_cast<Tetrahedron*>(const_cast<Element*>(neigh_info->getElement()));
528
529
	} else
	  break;
530
      }
531
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
532
   
533
534
    if (neigh == el) {
      (*el_info) = neigh_info;
535
      return 0;
536
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
537
538

   
539
540
541
    /****************************************************************************/
    /*  the domain's boundary is reached; loop back to the starting el          */
    /****************************************************************************/
542
    
Thomas Witkowski's avatar
Thomas Witkowski committed
543
    i = *n_neigh - 1;
544
    opp_v = refineList->getOppVertex(i, 0);
545
    do {
Thomas Witkowski's avatar
Thomas Witkowski committed
546
547
      TEST_EXIT_DBG(neigh_info->getNeighbour(opp_v) && i > 0)
	("While looping back domains boundary was reached or i == 0\n");
548
549
550
      opp_v = refineList->getOppVertex(i--, 0);
      neigh_info = stack->traverseNeighbour3d(neigh_info, opp_v);
    } while (neigh_info->getElement() != el);
Thomas Witkowski's avatar
Thomas Witkowski committed
551

552
    (*el_info) = neigh_info;
553
554
    
    return 1;
555
556
557
558
559
560
  }


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

562
    int bound = false;
563
    DegreeOfFreedom *edge[2];
564
565

    if (el_info->getElement()->getMark() <= 0)  
566
      return el_info;    /*   element may not be refined   */
567
568
569
570

    /****************************************************************************/
    /*  get memory for a list of all elements at the refinement edge            */
    /****************************************************************************/
571
    RCNeighbourList *ref_list = new RCNeighbourList(mesh->getMaxEdgeNeigh());
572
    ref_list->setElType(0, el_info->getType());
573
    ref_list->setElement(0, el_info->getElement());
574
575
576
577
578
    int n_neigh = 1;

    if (el_info->getProjection(0) && 
	el_info->getProjection(0)->getType() == VOLUME_PROJECTION)
      newCoords = true;
579
580
581
582
583

    /****************************************************************************/
    /*  give the refinement edge the right orientation                          */
    /****************************************************************************/

584
    if (el_info->getElement()->getDOF(0,0) < el_info->getElement()->getDOF(1,0)) {
585
586
      edge[0] = const_cast<int*>(el_info->getElement()->getDOF(0));
      edge[1] = const_cast<int*>(el_info->getElement()->getDOF(1));
587
    } else {
588
589
      edge[1] = const_cast<int*>(el_info->getElement()->getDOF(0));
      edge[0] = const_cast<int*>(el_info->getElement()->getDOF(1));
590
    }
591
592
593
594

    /****************************************************************************/
    /*  get the refinement patch                                                */
    /****************************************************************************/
595
596
597
598
599
600
601
    if (getRefinePatch(&el_info, edge, 0, ref_list, &n_neigh)) {      
      /****************************************************************************/
      /*  domain's boundary was reached while looping around the refinement edge  */
      /****************************************************************************/
      getRefinePatch(&el_info, edge, 1, ref_list, &n_neigh);
      bound = true;
    }
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621

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

    // ==========================================================================
    // === check for periodic boundary ==========================================
    // ==========================================================================

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

    DegreeOfFreedom newDOF = -1;
    DegreeOfFreedom lastNewDOF = -1;
    DegreeOfFreedom firstNewDOF = -1;

    RCNeighbourList *periodicList;
622
623
    std::map<int, VertexVector*>::iterator it;
    std::map<int, VertexVector*>::iterator end = mesh->getPeriodicAssociations().end();
624

Thomas Witkowski's avatar
Thomas Witkowski committed
625
    while (edge[0] != NULL) {
626
627
628
629
630
      periodicList = ref_list->periodicSplit(edge, 
					     next_edge,
					     &n_neigh,
					     &n_neigh_periodic);

631
      TEST_EXIT_DBG(periodicList)("periodicList = NULL\n");
632

633
      newDOF = refinePatch(edge, periodicList, n_neigh_periodic, bound);
634

635
      if (firstNewDOF == -1)
636
637
	firstNewDOF = newDOF;

638
639
640
641
      if (lastNewDOF != -1) {
	for (it = mesh->getPeriodicAssociations().begin(); it != end; ++it) {
	  if (it->second) {
	    if (((*(it->second))[edge[0][0]] == last_edge[0][0] &&
Thomas Witkowski's avatar
Thomas Witkowski committed
642
643
644
645
646
647
		 (*(it->second))[edge[1][0]] == last_edge[1][0]) || 
		((*(it->second))[edge[0][0]] == last_edge[1][0] &&
		 (*(it->second))[edge[1][0]] == last_edge[0][0])) {
	      (*(it->second))[lastNewDOF] = newDOF;
	      (*(it->second))[newDOF] = lastNewDOF;
	    } 
648
649
650
651
652
653
654
655
656
657
658
659
	  }
	}
      }
      lastNewDOF = newDOF;

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

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

660
661
662
663
    if (lastNewDOF != firstNewDOF) {
      for (it = mesh->getPeriodicAssociations().begin(); it != end; ++it) {
	if (it->second) {
	  if (((*(it->second))[first_edge[0][0]] == last_edge[0][0] &&
Thomas Witkowski's avatar
Thomas Witkowski committed
664
665
666
667
668
669
	       (*(it->second))[first_edge[1][0]] == last_edge[1][0]) || 
	      ((*(it->second))[first_edge[0][0]] == last_edge[1][0] &&
	       (*(it->second))[first_edge[1][0]] == last_edge[0][0])) {	    
	    (*(it->second))[lastNewDOF] = firstNewDOF;
	    (*(it->second))[firstNewDOF] = lastNewDOF;
	  }
670
671
672
673
674
675
676
677
678
679
680
681
	}
      }
    }
  

    /****************************************************************************/
    /*  and now refine the patch                                                */
    /****************************************************************************/
    //refinePatch(edge, ref_list, n_neigh, bound);

    stack->update();

Thomas Witkowski's avatar
Thomas Witkowski committed
682
    delete ref_list;
683

Thomas Witkowski's avatar
Thomas Witkowski committed
684
    return el_info;
685
686
687
  }

}