RefinementManager2d.cc 10.8 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
24
//     if (MPI::COMM_WORLD.Get_rank() == 0)
//       std::cout << "IN EL = " << el_info->getElement()->getIndex() << std::endl;

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

29
    if (el_info->getElement()->getMark() <= 0) {
30
31
32
//       if (MPI::COMM_WORLD.Get_rank() == 0)
// 	std::cout << "RETURN!" << std::endl;
      
33
      delete refineList;
34
35
36

      // Element may not be refined.
      return el_info;    
37
    }
38
39

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

42
    if (el_info->getProjection(0) && 
43
	el_info->getProjection(0)->getType() == VOLUME_PROJECTION)
44
      newCoords = true;
45
  
46
    // === Give the refinement edge the right orientation. ===
47

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

56
57
    // === Get the refinement patch. ===

58
    if (getRefinePatch(&el_info, edge, 0, refineList, &n_neigh)) {
59
      // Domain's boundary was reached while looping around the refinement edge.
60
61
62
63
      getRefinePatch(&el_info, edge, 1, refineList, &n_neigh);
      bound = true;
    }
    
64
    // === Check for periodic boundary ===
65
66
67
68
69
70
71
72
73
74
75

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

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

83
      TEST_EXIT_DBG(periodicList)("periodicList = NULL\n");
84

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

87
      delete periodicList;
88

89
      if (firstNewDOF == -1)
90
	firstNewDOF = newDOF;
91
92
      
      if (lastNewDOF != -1) {
93
	int i = 0;
94
	for (it = mesh->getPeriodicAssociations().begin(); it != end; ++it) {
95
96
97
// 	  if (MPI::COMM_WORLD.Get_rank() == 0)
// 	    std::cout << "i = " << i++ << std::endl;

98
	  if (it->second) {
99
100
101
// 	    if (MPI::COMM_WORLD.Get_rank() == 0)
// 	      std::cout << "*: " << newDOF << " " << lastNewDOF << std::endl;

102
103
104
	    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] &&
105
106
107
108
		 (*(it->second))[edge[1][0]] == last_edge[0][0])) {
	      (*(it->second))[lastNewDOF] = newDOF;
	      (*(it->second))[newDOF] = lastNewDOF;
	    } 
109
110
111
112
113
114
115
116
117
118
119
120
	  }
	}
      }
      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
    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] &&
127
128
129
130
131
132
133
	       (*(it->second))[first_edge[1][0]] == last_edge[0][0])) {
// 	    if (MPI::COMM_WORLD.Get_rank() == 0)
// 	      std::cout << "**: " << newDOF << " " << lastNewDOF << std::endl;

	    (*(it->second))[lastNewDOF] = firstNewDOF;
	    (*(it->second))[firstNewDOF] = lastNewDOF;
	  }
134
135
136
137
	}
      }
    }
  
138
    delete refineList;
139

140
    return el_info;
141
142
143
  }


144
  void RefinementManager2d::newCoordsFct(ElInfo *el_info)
145
146
147
148
149
150
  {
    Element *el = el_info->getElement();
    int dow = Global::getGeo(WORLD);

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

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

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

  void RefinementManager2d::setNewCoords()
  {
167
168
169
170
171
172
173
    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);
    }
174
175
176
177
178
179
180
181
182
  }


  DegreeOfFreedom RefinementManager2d::refinePatch(DegreeOfFreedom *edge[2], 
						   RCNeighbourList* refineList,
						   int n_neigh, bool bound)
  {
    FUNCNAME("RefinementManager2d::refinePatch()");
    DegreeOfFreedom *dof[3] = {NULL, NULL, NULL};
183
184
185
186
    Triangle *el = 
      dynamic_cast<Triangle*>(const_cast<Element*>(refineList->getElement(0)));
    Triangle *neigh = 
      dynamic_cast<Triangle*>(const_cast<Element*>(refineList->getElement(1)));
187

188
    // === There is one new vertex in the refinement edge. ===
189
190
191
192
193
194

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

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

195
    if (mesh->getNumberOfDOFs(EDGE)) {
196
      // There are two additional dofs in the refinement edge.
197
198
199
      dof[1] = mesh->getDOF(EDGE);
      dof[2] = mesh->getDOF(EDGE);
    }
200

201
    // === First refine the element. ===
202
203
204
205
206

    bisectTriangle(el, dof);

    if (neigh) {
      DegreeOfFreedom *tmp = dof[1];
207
208
209

      // === There is a neighbour; refine it also, but first exchange dof[1] and ===
      // === dof[2]; thus, dof[1] is always added on child[0]!                   ===
210
211
212
213
      dof[1] = dof[2];
      dof[2] = tmp;

      bisectTriangle(neigh, dof);
214
215
    } else {
      newCoords = true;
216
217
    }

218
    // === If there are functions to interpolate data to the finer grid, do so.
219
220
  
    int nrAdmin = mesh->getNumberOfDOFAdmin();
221
    for(int iadmin = 0; iadmin < nrAdmin; iadmin++) {
222
      DOFAdmin* admin = const_cast<DOFAdmin*>(&mesh->getDOFAdmin(iadmin));
223
      std::list<DOFIndexedBase*>::iterator end = admin->endDOFIndexed();
224
225
      for (std::list<DOFIndexedBase*>::iterator it = admin->beginDOFIndexed(); 
	   it != end; it++)
226
227
228
229
	(*it)->refineInterpol(*refineList, n_neigh);
    }


230
    if (!mesh->queryCoarseDOFs()) {
231
232
      // === If there should be no dof information on interior leaf elements ===
      // === remove dofs from edges and the centers of parents.              ===
233
234
235
      if (mesh->getNumberOfDOFs(EDGE)) {
	int node = mesh->getNode(EDGE);
	
236
237
	// === The only DOF that can be freed is that in the refinement edge; all ===
	// === other DOFs are handed on the children.                             ===
238
	
239
	mesh->freeDOF(const_cast<int*>(el->getDOF(node+2)), EDGE);
240
      }
241
      if (mesh->getNumberOfDOFs(EDGE) || mesh->getNumberOfDOFs(CENTER))
242
243
	refineList->removeDOFParents(n_neigh);
    }
244
245
246
247
248

    return dof[0][0];
  }


249
  void RefinementManager2d::bisectTriangle(Triangle *el, DegreeOfFreedom* newDOFs[3])
250
  {
251
    FUNCNAME("RefinementManager2d::bisectTriangle()");
252
    TEST_EXIT_DBG(mesh)("no mesh!\n");
253
254
255
256
257
  
    Triangle *child[2];

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

    signed char newMark = max(0, el->getMark() - 1);
260
261
262
263
    child[0]->setMark(newMark);
    child[1]->setMark(newMark);
    el->setMark(0);

264
    // === Transfer hidden data from parent to children. ===
265

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

271
272
    if (newMark > 0) 
      doMoreRecursiveRefine = true;
273

274
    // === Vertex 2 is the newest vertex. ===
275
276
277
278

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

279
    // === The other vertices are handed on from the parent. ===
280

281
    for (int i_child = 0; i_child < 2; i_child++) {
282
283
      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)));
284
    }
285

286
287
    // === There is one more leaf element, two hierachical elements and one ===
    // === more edge.                                                       ===
288
289
290
291
292

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

293
    if (mesh->getNumberOfDOFs(EDGE)) {
294
295
296
      DegreeOfFreedom* newEdgeDOFs = mesh->getDOF(EDGE);     

      // There are dofs in the midpoint of the edges.
297
298
299
      child[0]->setDOF(4, newEdgeDOFs);
      child[1]->setDOF(3, newEdgeDOFs);
      
300
301
302
      // 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)));
303
      
304
      // Dofs in the refinement edge.
305
306
307
308
309
      child[0]->setDOF(3, newDOFs[1]);
      child[1]->setDOF(4, newDOFs[2]);
    }
    
    if (mesh->getNumberOfDOFs(CENTER)) {
310
      int node = mesh->getNode(CENTER);
311
      
312
      // There are dofs at the barycenter of the triangles.
313
314
315
      child[0]->setDOF(node, mesh->getDOF(CENTER));
      child[1]->setDOF(node, mesh->getDOF(CENTER));
    }
316
317
318
319
320
321
322
  }

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

325
    if ((*el_info)->getNeighbour(2) && (*el_info)->getOppVertex(2) != 2) {
326
327
      // Neighbour is not compatible devisible; refine neighbour first, store the
      // opp_vertex to traverse back to el.
328
      int opp_vertex = (*el_info)->getOppVertex(2);
329
      
330
      ElInfo *neigh_info = stack->traverseNeighbour2d(*el_info, 2);
331
      neigh_info->getElement()->setMark(max(neigh_info->getElement()->getMark(), 1));
332
333
      neigh_info = refineFunction(neigh_info);
      
334
      // Now go back to the original element and refine the patch.
335
      *el_info = neigh_info = stack->traverseNeighbour2d(neigh_info, opp_vertex);
336
337
338
      TEST_EXIT_DBG(neigh_info->getElement() == 
		    dynamic_cast<Triangle*>(const_cast<Element*>((*el_info)->getElement())))
	("invalid traverse_neighbour1\n");
339
    }
340
341
  
    if (refineList->setElement(1, (*el_info)->getNeighbour(2))) {
342
      TEST_EXIT_DBG((*el_info)->getOppVertex(2) == 2)
343
344
345
346
	("no compatible ref. edge after recursive refinement of neighbour\n");
      *n_neigh = 2;
    }

347
    return 0;
348
349
350
  }

}