RefinementManager2d.cc 14.4 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
      newDOF = refinePatch(edge, periodicList, n_neigh_periodic, bound);
89

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

117
118
119
120
121
122
123
    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]))
124
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
	    {
	      (*(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);

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

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

162
163
164
      el->setNewCoord(new_coord);
    }
    
165
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
    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);

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

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

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


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

    return dof[0][0];
  }


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

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

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

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

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

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

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

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

    Triangle *el = dynamic_cast<Triangle*>(const_cast<Element*>( (*el_info)->getElement()));
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
      neigh_info->getElement()->setMark(max(neigh_info->getElement()->getMark(), 1));
369
370
371
372
373
374
      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);
375
      TEST_EXIT_DBG(neigh_info->getElement() == el)("invalid traverse_neighbour1\n");
376
    }
377
378
  
    if (refineList->setElement(1, (*el_info)->getNeighbour(2))) {
379
      TEST_EXIT_DBG((*el_info)->getOppVertex(2) == 2)
380
381
382
383
384
385
386
387
	("no compatible ref. edge after recursive refinement of neighbour\n");
      *n_neigh = 2;
    }

    return(0);
  }

}