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

namespace AMDiS {

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

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

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

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

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

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

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

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

53
    getRefinePatch(&elInfo, edge, 0, refineList, &n_neigh);
54

55

56
    // === Check for periodic boundary ===
57

58
    DegreeOfFreedom *next_edge[2] = {NULL, NULL};
59
60
61
62
63
64
65
66
67
    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;
68
69
    std::map<int, VertexVector*>::iterator it;
    std::map<int, VertexVector*>::iterator end = mesh->getPeriodicAssociations().end();
70

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

75
      TEST_EXIT_DBG(periodicList)("periodicList = NULL\n");
76

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

79
      delete periodicList;
80

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

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

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

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

122
    return elInfo;
123
124
125
  }


126
  void RefinementManager2d::newCoordsFct(ElInfo *elInfo)
127
  {
128
    Element *el = elInfo->getElement();
129
130
    int dow = Global::getGeo(WORLD);

131
    Projection *projector = elInfo->getProjection(0);
132

133
    if (!projector || projector->getType() != VOLUME_PROJECTION)
134
      projector = elInfo->getProjection(2);   
135

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

147

148
149
  void RefinementManager2d::setNewCoords()
  {
150
151
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, 
152
153
					 Mesh::CALL_EVERY_EL_PREORDER | 
					 Mesh::FILL_BOUND | Mesh::FILL_COORDS);
154
155
156
157
    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
 
237
    TEST_EXIT_DBG(mesh)("no mesh!\n");
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
  void RefinementManager2d::getRefinePatch(ElInfo **elInfo, 
304
305
306
307
					  DegreeOfFreedom *edge[2],
					  int dir, RCNeighbourList* refineList, 
					  int *n_neigh)
  {
308
    FUNCNAME("RefinementManager2d::getRefinePatch()");
309

310
    if ((*elInfo)->getNeighbour(2) && (*elInfo)->getOppVertex(2) != 2) {
311
312
      // Neighbour is not compatible devisible; refine neighbour first, store the
      // opp_vertex to traverse back to el.
313
      int opp_vertex = (*elInfo)->getOppVertex(2);
314
      
315
      ElInfo *neigh_info = stack->traverseNeighbour2d(*elInfo, 2);
316
      neigh_info->getElement()->setMark(max(neigh_info->getElement()->getMark(), 1));
317
      neigh_info = refineFunction(neigh_info);
318
319
320
321
322
323
324
325
326
327
328
329
330
331

      // === Now go back to the original element and refine the patch. ===


      // In Debug mode we test if traverNeighbour2d returns the expected element.
      int testIndex = neigh_info->getNeighbour(opp_vertex)->getIndex();

      *elInfo = neigh_info = stack->traverseNeighbour2d(neigh_info, opp_vertex);


      TEST_EXIT_DBG(testIndex == (*elInfo)->getElement()->getIndex())
	("Got wrong neighbour! Should be %d, but is %d!\n", 
	 testIndex, (*elInfo)->getElement()->getIndex());

332
      TEST_EXIT_DBG(neigh_info->getElement() == 
333
		    dynamic_cast<Triangle*>(const_cast<Element*>((*elInfo)->getElement())))
334
	("invalid traverse_neighbour1\n");
335
    }
336
337
338
 
    if (refineList->setElement(1, (*elInfo)->getNeighbour(2))) {
      TEST_EXIT_DBG((*elInfo)->getOppVertex(2) == 2)
339
340
341
342
343
344
	("no compatible ref. edge after recursive refinement of neighbour\n");
      *n_neigh = 2;
    }
  }

}