RefinementManager2d.cc 14.3 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
29
30
31
32
33
34
    bool bound = false;

    DegreeOfFreedom *edge[2];

    RCNeighbourList* refineList = NEW RCNeighbourList(2);

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

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

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

44
45
46
47
48
49
50
    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));
    }
51
52
53
54

    /****************************************************************************/
    /*  get the refinement patch                                                */
    /****************************************************************************/
55
56
57
58
59
60
61
62
    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;
    }
    
63
64
65
66
67
68
69
70
71
72
73
74
75
76
    // ==========================================================================
    // === 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;
77
78
    std::map<int, VertexVector*>::iterator it;
    std::map<int, VertexVector*>::iterator end = mesh->getPeriodicAssociations().end();
79

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

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

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

91
      if (firstNewDOF == -1) {
92
93
	firstNewDOF = newDOF;
      }
94
95
96
97
98
99
100
101
      
      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]))
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
	      {
		(*(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];
    }

118
119
120
121
122
123
124
    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]))
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
	    {
	      (*(it->second))[lastNewDOF] = firstNewDOF;
	      (*(it->second))[firstNewDOF] = lastNewDOF;
	    }
	}
      }
    }
  

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

    DELETE refineList;

    return(el_info);
  }


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

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

151
    if (!projector || projector->getType() != VOLUME_PROJECTION) {
152
153
154
      projector = el_info->getProjection(2);
    }

155
156
157
158
159
160
161
    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);
162

163
164
165
      el->setNewCoord(new_coord);
    }
    
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
    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};
    Triangle        *el = dynamic_cast<Triangle*>(const_cast<Element*>( refineList->getElement(0)));
    Triangle        *neigh = dynamic_cast<Triangle*>(const_cast<Element*>( refineList->getElement(1)));

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

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

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

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

    /****************************************************************************/
    /*  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);
221
222
    } else {
      newCoords = true;
223
224
225
226
227
228
229
230
231
    }

    /****************************************************************************/
    /*  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++) {
232
      std::list<DOFIndexedBase*>::iterator it;
233
      DOFAdmin* admin = const_cast<DOFAdmin*>(&mesh->getDOFAdmin(iadmin));
234
      std::list<DOFIndexedBase*>::iterator end = admin->endDOFIndexed();
235
236
237
238
239
240
      for(it = admin->beginDOFIndexed(); it != end; it++) {
	(*it)->refineInterpol(*refineList, n_neigh);
      }
    }


241
242
243
244
245
246
247
248
249
    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);
	
250
	/****************************************************************************/
251
252
	/*  the only DOF that can be freed is that in the refinement edge; all other*/
	/*  DOFs are handed on the children                                         */
253
	/****************************************************************************/
254
255
	
	mesh->freeDOF(const_cast<int*>( el->getDOF(node+2)), EDGE);
256
      }
257
258
259
260
      if (mesh->getNumberOfDOFs(EDGE)  ||  mesh->getNumberOfDOFs(CENTER)) {
	refineList->removeDOFParents(n_neigh);
      }
    }
261
262
263
264
265
266
267
268

    return dof[0][0];
  }


  void RefinementManager2d::bisectTriangle(Triangle *el, 
					   DegreeOfFreedom* newDOFs[3])
  {
269
    FUNCNAME("RefinementManager2d::bisectTriangle()");
270
    TEST_EXIT_DBG(mesh)("no mesh!\n");
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
  
    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]);

290
291
    if (newMark > 0) 
      doMoreRecursiveRefine = true;
292
293
294
295
296
297
298
299
300
301
302
303

    /****************************************************************************/
    /*  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                        */
    /****************************************************************************/

304
305
306
307
    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)));
    }
308
309
310
311
312
313
314
315
316
317

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

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

318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
    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)) {
341
      int node = mesh->getNode(CENTER);
342
343
344
345
346
347
348
      
      /****************************************************************************/
      /* there are dofs at the barycenter of the triangles                        */
      /****************************************************************************/
      child[0]->setDOF(node, mesh->getDOF(CENTER));
      child[1]->setDOF(node, mesh->getDOF(CENTER));
    }
349
350
351
352
353
354
355
  }

  int RefinementManager2d::getRefinePatch(ElInfo **el_info, 
					  DegreeOfFreedom *edge[2],
					  int dir, RCNeighbourList* refineList, 
					  int *n_neigh)
  {
356
    FUNCNAME("RefinementManager2d::getRefinePatch()");
357
  
358
    int opp_vertex = 0;
359

360
361
362
363
364
365
366
    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);
      
367
      ElInfo *neigh_info = stack->traverseNeighbour2d(*el_info, 2);
368
369
370
371
372
373
374
375
      neigh_info->getElement()->setMark(max(neigh_info->getElement()->getMark(), 
					    1));
      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);
376
      TEST_EXIT_DBG(neigh_info->getElement() == el)("invalid traverse_neighbour1\n");
377
    }
378
379
  
    if (refineList->setElement(1, (*el_info)->getNeighbour(2))) {
380
      TEST_EXIT_DBG((*el_info)->getOppVertex(2) == 2)
381
382
383
384
385
386
387
388
	("no compatible ref. edge after recursive refinement of neighbour\n");
      *n_neigh = 2;
    }

    return(0);
  }

}