RefinementManager2d.cc 10.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
    FUNCNAME("RefinementManager::refineFunction()");

22
23
    bool bound = false;
    DegreeOfFreedom *edge[2];
24
    RCNeighbourList* refineList = new RCNeighbourList(2);
25

26
    if (el_info->getElement()->getMark() <= 0) {
27
      delete refineList;
28
29
30

      // Element may not be refined.
      return el_info;    
31
    }
32
33

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

36
    if (el_info->getProjection(0) && 
37
	el_info->getProjection(0)->getType() == VOLUME_PROJECTION)
38
      newCoords = true;
39
  
40
    // === Give the refinement edge the right orientation. ===
41

42
    if (el_info->getElement()->getDOF(0,0) < el_info->getElement()->getDOF(1, 0)) {
43
44
      edge[0] = const_cast<int*>(el_info->getElement()->getDOF(0));
      edge[1] = const_cast<int*>(el_info->getElement()->getDOF(1));
45
    } else {
46
47
      edge[1] = const_cast<int*>(el_info->getElement()->getDOF(0));
      edge[0] = const_cast<int*>(el_info->getElement()->getDOF(1));
48
    }
49

50
51
    // === Get the refinement patch. ===

52
    if (getRefinePatch(&el_info, edge, 0, refineList, &n_neigh)) {
53
      // Domain's boundary was reached while looping around the refinement edge.
54
55
56
57
      getRefinePatch(&el_info, edge, 1, refineList, &n_neigh);
      bound = true;
    }
    
58
    // === Check for periodic boundary ===
59
60
61
62
63
64
65
66
67
68
69

    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;
70
71
    std::map<int, VertexVector*>::iterator it;
    std::map<int, VertexVector*>::iterator end = mesh->getPeriodicAssociations().end();
72

73
    while (edge[0] != NULL) {
74
75
      periodicList = refineList->periodicSplit(edge, next_edge,
					       &n_neigh, &n_neigh_periodic);
76

77
      TEST_EXIT_DBG(periodicList)("periodicList = NULL\n");
78

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

81
      delete periodicList;
82

83
      if (firstNewDOF == -1)
84
	firstNewDOF = newDOF;
85
86
87
88
89
90
91
      
      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] &&
92
93
94
95
		 (*(it->second))[edge[1][0]] == last_edge[0][0])) {
	      (*(it->second))[lastNewDOF] = newDOF;
	      (*(it->second))[newDOF] = lastNewDOF;
	    } 
96
97
98
99
100
101
102
103
104
105
106
107
	  }
	}
      }
      lastNewDOF = newDOF;

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

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

108
109
110
111
112
113
    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] &&
114
115
116
117
	       (*(it->second))[first_edge[1][0]] == last_edge[0][0])) {
	    (*(it->second))[lastNewDOF] = firstNewDOF;
	    (*(it->second))[firstNewDOF] = lastNewDOF;
	  }
118
119
120
121
	}
      }
    }
  
122
    delete refineList;
123

124
    return el_info;
125
126
127
  }


128
  void RefinementManager2d::newCoordsFct(ElInfo *el_info)
129
130
131
132
133
134
  {
    Element *el = el_info->getElement();
    int dow = Global::getGeo(WORLD);

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

135
136
    if (!projector || projector->getType() != VOLUME_PROJECTION)
      projector = el_info->getProjection(2);   
137

138
    if (el->getFirstChild() && projector && (!el->isNewCoordSet())) {
139
      WorldVector<double> *new_coord = new WorldVector<double>;
140
141
      
      for (int j = 0; j < dow; j++)
142
	(*new_coord)[j] = (el_info->getCoord(0)[j] + el_info->getCoord(1)[j]) * 0.5;
143
144
145
146
      
      projector->project(*new_coord);
      el->setNewCoord(new_coord);
    }
147
148
149
150
  }

  void RefinementManager2d::setNewCoords()
  {
151
152
153
154
155
156
157
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, 
					 Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_BOUND | Mesh::FILL_COORDS);
    while (elInfo) {
      newCoordsFct(elInfo);
      elInfo = stack.traverseNext(elInfo);
    }
158
159
160
161
162
163
164
165
166
  }


  DegreeOfFreedom RefinementManager2d::refinePatch(DegreeOfFreedom *edge[2], 
						   RCNeighbourList* refineList,
						   int n_neigh, bool bound)
  {
    FUNCNAME("RefinementManager2d::refinePatch()");
    DegreeOfFreedom *dof[3] = {NULL, NULL, NULL};
167
168
169
170
    Triangle *el = 
      dynamic_cast<Triangle*>(const_cast<Element*>(refineList->getElement(0)));
    Triangle *neigh = 
      dynamic_cast<Triangle*>(const_cast<Element*>(refineList->getElement(1)));
171

172
    // === There is one new vertex in the refinement edge. ===
173
174
175
176
177
178

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

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

179
    if (mesh->getNumberOfDOFs(EDGE)) {
180
      // There are two additional dofs in the refinement edge.
181
182
183
      dof[1] = mesh->getDOF(EDGE);
      dof[2] = mesh->getDOF(EDGE);
    }
184

185
    // === First refine the element. ===
186
187
188
189
190

    bisectTriangle(el, dof);

    if (neigh) {
      DegreeOfFreedom *tmp = dof[1];
191
192
193

      // === There is a neighbour; refine it also, but first exchange dof[1] and ===
      // === dof[2]; thus, dof[1] is always added on child[0]!                   ===
194
195
196
197
      dof[1] = dof[2];
      dof[2] = tmp;

      bisectTriangle(neigh, dof);
198
199
    } else {
      newCoords = true;
200
201
    }

202
    // === If there are functions to interpolate data to the finer grid, do so.
203
204
  
    int nrAdmin = mesh->getNumberOfDOFAdmin();
205
    for(int iadmin = 0; iadmin < nrAdmin; iadmin++) {
206
      DOFAdmin* admin = const_cast<DOFAdmin*>(&mesh->getDOFAdmin(iadmin));
207
      std::list<DOFIndexedBase*>::iterator end = admin->endDOFIndexed();
208
209
      for (std::list<DOFIndexedBase*>::iterator it = admin->beginDOFIndexed(); 
	   it != end; it++)
210
211
212
213
	(*it)->refineInterpol(*refineList, n_neigh);
    }


214
    if (!mesh->queryCoarseDOFs()) {
215
216
      // === If there should be no dof information on interior leaf elements ===
      // === remove dofs from edges and the centers of parents.              ===
217
218
219
      if (mesh->getNumberOfDOFs(EDGE)) {
	int node = mesh->getNode(EDGE);
	
220
221
	// === The only DOF that can be freed is that in the refinement edge; all ===
	// === other DOFs are handed on the children.                             ===
222
	
223
	mesh->freeDOF(const_cast<int*>(el->getDOF(node+2)), EDGE);
224
      }
225
      if (mesh->getNumberOfDOFs(EDGE) || mesh->getNumberOfDOFs(CENTER))
226
227
	refineList->removeDOFParents(n_neigh);
    }
228
229
230
231
232

    return dof[0][0];
  }


233
  void RefinementManager2d::bisectTriangle(Triangle *el, DegreeOfFreedom* newDOFs[3])
234
  {
235
    FUNCNAME("RefinementManager2d::bisectTriangle()");
236
    TEST_EXIT_DBG(mesh)("no mesh!\n");
237
238
239
240
241
  
    Triangle *child[2];

    child[0] = dynamic_cast<Triangle*>(mesh->createNewElement(el));
    child[1] = dynamic_cast<Triangle*>(mesh->createNewElement(el));
242
243

    signed char newMark = max(0, el->getMark() - 1);
244
245
246
247
    child[0]->setMark(newMark);
    child[1]->setMark(newMark);
    el->setMark(0);

248
    // === Transfer hidden data from parent to children. ===
249

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

255
256
    if (newMark > 0) 
      doMoreRecursiveRefine = true;
257

258
    // === Vertex 2 is the newest vertex. ===
259
260
261
262

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

263
    // === The other vertices are handed on from the parent. ===
264

265
    for (int i_child = 0; i_child < 2; i_child++) {
266
267
      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)));
268
    }
269

270
271
    // === There is one more leaf element, two hierachical elements and one ===
    // === more edge.                                                       ===
272
273
274
275
276

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

277
    if (mesh->getNumberOfDOFs(EDGE)) {
278
279
280
      DegreeOfFreedom* newEdgeDOFs = mesh->getDOF(EDGE);     

      // There are dofs in the midpoint of the edges.
281
282
283
      child[0]->setDOF(4, newEdgeDOFs);
      child[1]->setDOF(3, newEdgeDOFs);
      
284
285
286
      // 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)));
287
      
288
      // Dofs in the refinement edge.
289
290
291
292
293
      child[0]->setDOF(3, newDOFs[1]);
      child[1]->setDOF(4, newDOFs[2]);
    }
    
    if (mesh->getNumberOfDOFs(CENTER)) {
294
      int node = mesh->getNode(CENTER);
295
      
296
      // There are dofs at the barycenter of the triangles.
297
298
299
      child[0]->setDOF(node, mesh->getDOF(CENTER));
      child[1]->setDOF(node, mesh->getDOF(CENTER));
    }
300
301
302
303
304
305
306
  }

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

309
    if ((*el_info)->getNeighbour(2) && (*el_info)->getOppVertex(2) != 2) {
310
311
      // Neighbour is not compatible devisible; refine neighbour first, store the
      // opp_vertex to traverse back to el.
312
      int opp_vertex = (*el_info)->getOppVertex(2);
313
      
314
      ElInfo *neigh_info = stack->traverseNeighbour2d(*el_info, 2);
315
      neigh_info->getElement()->setMark(max(neigh_info->getElement()->getMark(), 1));
316
317
      neigh_info = refineFunction(neigh_info);
      
318
      // Now go back to the original element and refine the patch.
319
      *el_info = neigh_info = stack->traverseNeighbour2d(neigh_info, opp_vertex);
320
321
322
      TEST_EXIT_DBG(neigh_info->getElement() == 
		    dynamic_cast<Triangle*>(const_cast<Element*>((*el_info)->getElement())))
	("invalid traverse_neighbour1\n");
323
    }
324
325
  
    if (refineList->setElement(1, (*el_info)->getNeighbour(2))) {
326
      TEST_EXIT_DBG((*el_info)->getOppVertex(2) == 2)
327
328
329
330
	("no compatible ref. edge after recursive refinement of neighbour\n");
      *n_neigh = 2;
    }

331
    return 0;
332
333
334
  }

}