RefinementManager2d.cc 14.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#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 "LeafData.h"
#include "VertexVector.h"

namespace AMDiS {

  ElInfo* RefinementManager2d::refineFunction(ElInfo* el_info)
  {
20
21
22
    FUNCNAME("RefinementManager::refineFunction()");

    int n_neigh;
23
24
25
26
27
28
    bool bound = false;

    DegreeOfFreedom *edge[2];

    RCNeighbourList* refineList = NEW RCNeighbourList(2);

29
30
31
32
    if (el_info->getElement()->getMark() <= 0) {
      DELETE refineList;
      return el_info;    /*   element may not be refined   */
    }
33
34
35
36

    refineList->setElement(0, el_info->getElement());
    n_neigh = 1;

37
38
39
40
    if (el_info->getProjection(0) && 
	el_info->getProjection(0)->getType() == VOLUME_PROJECTION) {
      newCoords = true;
    }
41
42
43
44
45
  
    /****************************************************************************/
    /*  give the refinement edge the right orientation                          */
    /****************************************************************************/

46
47
48
49
50
51
52
    if (el_info->getElement()->getDOF(0,0) < el_info->getElement()->getDOF(1, 0)) {
      edge[0] = const_cast<int*>( el_info->getElement()->getDOF(0));
      edge[1] = const_cast<int*>( el_info->getElement()->getDOF(1));
    } else {
      edge[1] = const_cast<int*>( el_info->getElement()->getDOF(0));
      edge[0] = const_cast<int*>( el_info->getElement()->getDOF(1));
    }
53
54
55
56

    /****************************************************************************/
    /*  get the refinement patch                                                */
    /****************************************************************************/
57
58
59
60
61
62
63
64
    if (getRefinePatch(&el_info, edge, 0, refineList, &n_neigh)) {
      /****************************************************************************/
      /*  domain's boundary was reached while looping around the refinement edge  */
      /****************************************************************************/
      getRefinePatch(&el_info, edge, 1, refineList, &n_neigh);
      bound = true;
    }
    
65
66
67
68
69
70
71
72
73
74
75
76
77
78
    // ==========================================================================
    // === 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;
79
80
    std::map<int, VertexVector*>::iterator it;
    std::map<int, VertexVector*>::iterator end = mesh->getPeriodicAssociations().end();
81

82
    while (edge[0] != NULL) {
83
84
85
86
87
      periodicList = refineList->periodicSplit(edge, 
					       next_edge,
					       &n_neigh,
					       &n_neigh_periodic);

88
      TEST_EXIT_DBG(periodicList)("periodicList = NULL\n");
89

90
      newDOF = refinePatch(edge, periodicList, n_neigh_periodic, bound);
91

92
93
      DELETE periodicList;

94
      if (firstNewDOF == -1) {
95
96
	firstNewDOF = newDOF;
      }
97
98
99
100
101
102
103
104
      
      if (lastNewDOF != -1) {
	for (it = mesh->getPeriodicAssociations().begin(); it != end; ++it) {
	  if (it->second) {
	    if (((*(it->second))[edge[0][0]] == last_edge[0][0] &&
		 (*(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]))
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
	      {
		(*(it->second))[lastNewDOF] = newDOF;
		(*(it->second))[newDOF] = lastNewDOF;
	      } 
	  }
	}
      }
      lastNewDOF = newDOF;

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

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

121
122
123
124
125
126
127
    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] &&
	       (*(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]))
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
	    {
	      (*(it->second))[lastNewDOF] = firstNewDOF;
	      (*(it->second))[firstNewDOF] = lastNewDOF;
	    }
	}
      }
    }
  

    /****************************************************************************/
    /*  and now refine the patch                                                */
    /****************************************************************************/

    DELETE refineList;

143
    return el_info;
144
145
146
147
148
149
150
151
152
153
  }


  int RefinementManager2d::newCoordsFct(ElInfo *el_info)
  {
    Element *el = el_info->getElement();
    int dow = Global::getGeo(WORLD);

    Projection *projector = el_info->getProjection(0);

154
    if (!projector || projector->getType() != VOLUME_PROJECTION) {
155
156
157
      projector = el_info->getProjection(2);
    }

158
159
160
161
162
163
164
    if (el->getFirstChild() && projector && (!el->isNewCoordSet())) {
      WorldVector<double> *new_coord = NEW WorldVector<double>;
      
      for (int j = 0; j < dow; j++)
	(*new_coord)[j] = (el_info->getCoord(0)[j] + el_info->getCoord(1)[j])*0.5;
      
      projector->project(*new_coord);
165

166
167
168
      el->setNewCoord(new_coord);
    }
    
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
    return 0;
  }

  void RefinementManager2d::setNewCoords()
  {
    Flag fillFlag = Mesh::CALL_EVERY_EL_PREORDER|
      Mesh::FILL_BOUND|
      Mesh::FILL_COORDS;
  
    mesh->traverse(-1, fillFlag, newCoordsFct);
  }


  DegreeOfFreedom RefinementManager2d::refinePatch(DegreeOfFreedom *edge[2], 
						   RCNeighbourList* refineList,
						   int n_neigh, bool bound)
  {
    FUNCNAME("RefinementManager2d::refinePatch()");
    DegreeOfFreedom *dof[3] = {NULL, NULL, NULL};
188
189
    Triangle *el = dynamic_cast<Triangle*>(const_cast<Element*>(refineList->getElement(0)));
    Triangle *neigh = dynamic_cast<Triangle*>(const_cast<Element*>(refineList->getElement(1)));
190
191
192
193
194
195
196
197
198
199

    /****************************************************************************/
    /*  there is one new vertex in the refinement edge                          */
    /****************************************************************************/

    dof[0] = mesh->getDOF(VERTEX);

    mesh->incrementNumberOfVertices(1);
    mesh->incrementNumberOfEdges(1);

200
201
202
203
204
205
206
    if (mesh->getNumberOfDOFs(EDGE)) {
      /****************************************************************************/
      /*  there are two additional dofs in the refinement edge                    */
      /****************************************************************************/
      dof[1] = mesh->getDOF(EDGE);
      dof[2] = mesh->getDOF(EDGE);
    }
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223

    /****************************************************************************/
    /*  first refine the element                                                */
    /****************************************************************************/

    bisectTriangle(el, dof);

    if (neigh) {
      DegreeOfFreedom *tmp = dof[1];
      /****************************************************************************/
      /*  there is a neighbour; refine it also, but first exchange dof[1] and     */
      /*  dof[2]; thus, dof[1] is always added on child[0]!                       */
      /****************************************************************************/
      dof[1] = dof[2];
      dof[2] = tmp;

      bisectTriangle(neigh, dof);
224
225
    } else {
      newCoords = true;
226
227
228
229
230
231
232
233
234
    }

    /****************************************************************************/
    /*  if there are functions to interpolate data to the finer grid, do so     */
    /****************************************************************************/
  
    int iadmin;
    int nrAdmin = mesh->getNumberOfDOFAdmin();
    for(iadmin = 0; iadmin < nrAdmin; iadmin++) {
235
      std::list<DOFIndexedBase*>::iterator it;
236
      DOFAdmin* admin = const_cast<DOFAdmin*>(&mesh->getDOFAdmin(iadmin));
237
      std::list<DOFIndexedBase*>::iterator end = admin->endDOFIndexed();
238
239
240
241
242
243
      for(it = admin->beginDOFIndexed(); it != end; it++) {
	(*it)->refineInterpol(*refineList, n_neigh);
      }
    }


244
245
246
247
248
249
250
251
252
    if (!mesh->queryCoarseDOFs()) {
      /****************************************************************************/
      /*  if there should be no dof information on interior leaf elements remove  */
      /*  dofs from edges and the centers of parents                              */
      /****************************************************************************/
      if (mesh->getNumberOfDOFs(EDGE)) {
	
	int node = mesh->getNode(EDGE);
	
253
	/****************************************************************************/
254
255
	/*  the only DOF that can be freed is that in the refinement edge; all other*/
	/*  DOFs are handed on the children                                         */
256
	/****************************************************************************/
257
258
	
	mesh->freeDOF(const_cast<int*>( el->getDOF(node+2)), EDGE);
259
      }
260
261
262
263
      if (mesh->getNumberOfDOFs(EDGE)  ||  mesh->getNumberOfDOFs(CENTER)) {
	refineList->removeDOFParents(n_neigh);
      }
    }
264
265
266
267
268
269
270
271

    return dof[0][0];
  }


  void RefinementManager2d::bisectTriangle(Triangle *el, 
					   DegreeOfFreedom* newDOFs[3])
  {
272
    FUNCNAME("RefinementManager2d::bisectTriangle()");
273
    TEST_EXIT_DBG(mesh)("no mesh!\n");
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
  
    Triangle *child[2];

    child[0] = dynamic_cast<Triangle*>(mesh->createNewElement(el));
    child[1] = dynamic_cast<Triangle*>(mesh->createNewElement(el));
  
    signed char newMark = max(0, el->getMark()-1);
    child[0]->setMark(newMark);
    child[1]->setMark(newMark);
    el->setMark(0);

    /****************************************************************************/
    /*  transfer hidden data from parent to children                            */
    /****************************************************************************/

    el->refineElementData(child[0], child[1]); // call of subclass-method
    el->setFirstChild(child[0]);
    el->setSecondChild(child[1]);

293
294
    if (newMark > 0) 
      doMoreRecursiveRefine = true;
295
296
297
298
299
300
301
302
303
304
305
306

    /****************************************************************************/
    /*  vertex 2 is the newest vertex                                           */
    /****************************************************************************/

    child[0]->setDOF(2, newDOFs[0]);
    child[1]->setDOF(2, newDOFs[0]);

    /****************************************************************************/
    /*  the other vertices are handed on from the parent                        */
    /****************************************************************************/

307
308
309
310
    for (int i_child = 0; i_child < 2; i_child++) {
      child[i_child]->setDOF(i_child, const_cast<int*>( el->getDOF(2)));
      child[i_child]->setDOF(1-i_child, const_cast<int*>( el->getDOF(i_child)));
    }
311
312
313
314
315
316
317
318
319
320

    /****************************************************************************/
    /*  there is one more leaf element, two hierachical elements and one more   */
    /*  edge                                        			    */
    /****************************************************************************/

    mesh->incrementNumberOfEdges(1);
    mesh->incrementNumberOfLeaves(1);
    mesh->incrementNumberOfElements(2);

321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
    if (mesh->getNumberOfDOFs(EDGE)) {
      /****************************************************************************/
      /*  there are dof's in the midpoint of the edges			    */
      /****************************************************************************/
      DegreeOfFreedom* newEdgeDOFs = mesh->getDOF(EDGE);
      
      child[0]->setDOF(4, newEdgeDOFs);
      child[1]->setDOF(3, newEdgeDOFs);
      
      /****************************************************************************/
      /*  dofs handed on by the parent                                            */
      /****************************************************************************/
      child[0]->setDOF(5, const_cast<int*>( el->getDOF(4)));
      child[1]->setDOF(5, const_cast<int*>( el->getDOF(3)));
      
      /****************************************************************************/
      /*  dofs in the refinement edge                                             */
      /****************************************************************************/
      child[0]->setDOF(3, newDOFs[1]);
      child[1]->setDOF(4, newDOFs[2]);
    }
    
    if (mesh->getNumberOfDOFs(CENTER)) {
344
      int node = mesh->getNode(CENTER);
345
346
347
348
349
350
351
      
      /****************************************************************************/
      /* there are dofs at the barycenter of the triangles                        */
      /****************************************************************************/
      child[0]->setDOF(node, mesh->getDOF(CENTER));
      child[1]->setDOF(node, mesh->getDOF(CENTER));
    }
352
353
354
355
356
357
358
  }

  int RefinementManager2d::getRefinePatch(ElInfo **el_info, 
					  DegreeOfFreedom *edge[2],
					  int dir, RCNeighbourList* refineList, 
					  int *n_neigh)
  {
359
    FUNCNAME("RefinementManager2d::getRefinePatch()");
360
361

    Triangle *el = dynamic_cast<Triangle*>(const_cast<Element*>( (*el_info)->getElement()));
362
    int opp_vertex = 0;
363

364
365
366
367
368
369
370
    if ((*el_info)->getNeighbour(2) && (*el_info)->getOppVertex(2) != 2) {
      /****************************************************************************/
      /*  neighbour is not compatible devisible; refine neighbour first, store the*/
      /*  opp_vertex to traverse back to el                                       */
      /****************************************************************************/
      opp_vertex = (*el_info)->getOppVertex(2);
      
371
      ElInfo *neigh_info = stack->traverseNeighbour2d(*el_info, 2);
372
      neigh_info->getElement()->setMark(max(neigh_info->getElement()->getMark(), 1));
373
374
375
376
377
378
      neigh_info = refineFunction(neigh_info);
      
      /****************************************************************************/
      /*  now go back to the original element and refine the patch                */
      /****************************************************************************/
      *el_info = neigh_info = stack->traverseNeighbour2d(neigh_info, opp_vertex);
379
      TEST_EXIT_DBG(neigh_info->getElement() == el)("invalid traverse_neighbour1\n");
380
    }
381
382
  
    if (refineList->setElement(1, (*el_info)->getNeighbour(2))) {
383
      TEST_EXIT_DBG((*el_info)->getOppVertex(2) == 2)
384
385
386
387
388
389
390
391
	("no compatible ref. edge after recursive refinement of neighbour\n");
      *n_neigh = 2;
    }

    return(0);
  }

}