RefinementManager2d.cc 11.9 KB
Newer Older
1
2
3
4
5
6
7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9
10
11
12
13
14
15
16
17
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
18
 *
19
 ******************************************************************************/
20
21


22
23
24
25
26
27
28
29
30
31
32
33
34
35
#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"
36
#include "Debug.h"
37
38
39

namespace AMDiS {

40
  ElInfo* RefinementManager2d::refineFunction(ElInfo* elInfo)
41
  {
42
    if (elInfo->getElement()->getMark() <= 0)
43
      return elInfo;
44

45
46
    bool bound = false;
    DegreeOfFreedom *edge[2];
47
48
    RCNeighbourList refineList(2);
    refineList.setElement(0, elInfo->getElement());
49
    int n_neigh = 1;
50

51
    if (elInfo->getProjection(0) &&
52
	elInfo->getProjection(0)->getType() == VOLUME_PROJECTION)
53
      newCoords = true;
54

55
    // === Give the refinement edge the right orientation. ===
56

57
    if (elInfo->getElement()->getDof(0, 0) < elInfo->getElement()->getDof(1, 0)) {
58
59
      edge[0] = const_cast<DegreeOfFreedom*>(elInfo->getElement()->getDof(0));
      edge[1] = const_cast<DegreeOfFreedom*>(elInfo->getElement()->getDof(1));
60
    } else {
61
62
      edge[1] = const_cast<DegreeOfFreedom*>(elInfo->getElement()->getDof(0));
      edge[0] = const_cast<DegreeOfFreedom*>(elInfo->getElement()->getDof(1));
63
    }
64

65
66
    // === Get the refinement patch. ===

67
    getRefinePatch(&elInfo, edge, 0, refineList, &n_neigh);
68

69

70
    // === Check for periodic boundary ===
71

72
    DegreeOfFreedom *next_edge[2] = {NULL, NULL};
73
    DegreeOfFreedom *first_edge[2] = {edge[0], edge[1]};
74
    DegreeOfFreedom *last_edge[2] = {NULL, NULL};
75
76
    int n_neigh_periodic;

77
78
79
    DegreeOfFreedom *newDOF = NULL;
    DegreeOfFreedom *lastNewDOF = NULL;
    DegreeOfFreedom *firstNewDOF = NULL;
80

81
    RCNeighbourList periodicList;
82

83
    while (edge[0] != NULL) {
84
85
86
      refineList.periodicSplit(edge, next_edge,
			       &n_neigh, &n_neigh_periodic,
			       periodicList);
87

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

90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
      if (firstNewDOF == NULL)
        firstNewDOF = newDOF;

      if (lastNewDOF != NULL) {
        for (std::map<int, std::vector<VertexVector*> >::iterator it = mesh->getPeriodicAssociations().begin();
            it != mesh->getPeriodicAssociations().end(); ++it) {
          if (!it->second.empty()) {
            for (size_t iadmin = 0; iadmin < mesh->getNumberOfDOFAdmin(); iadmin++) {
              if (((*(it->second)[iadmin])[edge[0][iadmin]] == last_edge[0][iadmin] &&
                   (*(it->second)[iadmin])[edge[1][iadmin]] == last_edge[1][iadmin]) ||
                   ((*(it->second[iadmin]))[edge[0][iadmin]] == last_edge[1][iadmin] &&
                   (*(it->second)[iadmin])[edge[1][iadmin]] == last_edge[0][iadmin])) {
                (*(it->second)[iadmin])[lastNewDOF[iadmin]] = newDOF[iadmin];
                (*(it->second)[iadmin])[newDOF[iadmin]] = lastNewDOF[iadmin];
              }
            } // for
          } // if
        } // for
      } // if lastNewDOF
109
110
111
112
113
114
115
116
117
      lastNewDOF = newDOF;

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

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

118
    if (lastNewDOF != firstNewDOF) {
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
      for (std::map<int, std::vector<VertexVector*> >::iterator it = mesh->getPeriodicAssociations().begin();
           it != mesh->getPeriodicAssociations().end(); ++it) {
        if (!it->second.empty()) {
          for (size_t iadmin = 0; iadmin < mesh->getNumberOfDOFAdmin(); iadmin++) {
            if (((*(it->second)[iadmin])[first_edge[0][iadmin]] == last_edge[0][iadmin] &&
                 (*(it->second)[iadmin])[first_edge[1][iadmin]] == last_edge[1][iadmin]) ||
                 ((*(it->second)[iadmin])[first_edge[0][iadmin]] == last_edge[1][iadmin] &&
                 (*(it->second)[iadmin])[first_edge[1][iadmin]] == last_edge[0][iadmin])) {
              (*(it->second)[iadmin])[lastNewDOF[iadmin]] = firstNewDOF[iadmin];
              (*(it->second)[iadmin])[firstNewDOF[iadmin]] = lastNewDOF[iadmin];
            }
          } // for
        } // if
      } // for
    } // if !=
134

135
    return elInfo;
136
137
138
  }


139
  void RefinementManager2d::newCoordsFct(ElInfo *elInfo)
140
  {
141
    Element *el = elInfo->getElement();
142
143
    int dow = Global::getGeo(WORLD);

144
    Projection *projector = elInfo->getProjection(0);
145

146
    if (!projector || projector->getType() != VOLUME_PROJECTION)
147
      projector = elInfo->getProjection(2);
148

149
    if (el->getFirstChild() && projector && (!el->isNewCoordSet())) {
150
      WorldVector<double> *new_coord = new WorldVector<double>;
151

152
      for (int j = 0; j < dow; j++)
153
	(*new_coord)[j] = (elInfo->getCoord(0)[j] + elInfo->getCoord(1)[j]) * 0.5;
154

155
156
157
      projector->project(*new_coord);
      el->setNewCoord(new_coord);
    }
158
159
  }

160

161
  void RefinementManager2d::setNewCoords(int macroEl)
162
  {
163
    TraverseStack stack;
164
165
166
    ElInfo *elInfo;

    if (macroEl == -1)
167
168
      elInfo = stack.traverseFirst(mesh, -1,
				   Mesh::CALL_EVERY_EL_PREORDER |
169
170
171
				   Mesh::FILL_BOUND | Mesh::FILL_COORDS);
    else
      elInfo = stack.traverseFirstOneMacro(mesh, macroEl, -1,
172
					   Mesh::CALL_EVERY_EL_PREORDER |
173
					   Mesh::FILL_BOUND | Mesh::FILL_COORDS);
174

175
176
177
178
    while (elInfo) {
      newCoordsFct(elInfo);
      elInfo = stack.traverseNext(elInfo);
    }
179
180
181
  }


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

192
    // === There is one new vertex in the refinement edge. ===
193

194
    dof[0] = mesh->getDof(VERTEX);
195
196
197
198

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

199
    if (mesh->getNumberOfDofs(EDGE)) {
200
      // There are two additional dofs in the refinement edge.
201
202
      dof[1] = mesh->getDof(EDGE);
      dof[2] = mesh->getDof(EDGE);
203
    }
204

205
    // === First refine the element. ===
206
207
208
209
210

    bisectTriangle(el, dof);

    if (neigh) {
      DegreeOfFreedom *tmp = dof[1];
211
212
213

      // === There is a neighbour; refine it also, but first exchange dof[1] and ===
      // === dof[2]; thus, dof[1] is always added on child[0]!                   ===
214
215
216
217
      dof[1] = dof[2];
      dof[2] = tmp;

      bisectTriangle(neigh, dof);
218
219
    } else {
      newCoords = true;
220
221
    }

222
    // === If there are functions to interpolate data to the finer grid, do so.
223

224
    int nrAdmin = mesh->getNumberOfDOFAdmin();
225
    for(int iadmin = 0; iadmin < nrAdmin; iadmin++) {
226
      DOFAdmin* admin = const_cast<DOFAdmin*>(&mesh->getDofAdmin(iadmin));
227
      std::list<DOFIndexedBase*>::iterator end = admin->endDOFIndexed();
228
      for (std::list<DOFIndexedBase*>::iterator it = admin->beginDOFIndexed();
229
230
           it != end; it++)
        (*it)->refineInterpol(refineList, n_neigh);
231
232
233
    }


234
    if (!mesh->queryCoarseDOFs()) {
235
236
      // === If there should be no dof information on interior leaf elements ===
      // === remove dofs from edges and the centers of parents.              ===
237
      if (mesh->getNumberOfDofs(EDGE)) {
238
        int node = mesh->getNode(EDGE);
239

240
241
        // === The only DOF that can be freed is that in the refinement edge; all ===
        // === other DOFs are handed on the children.                             ===
242

243
        mesh->freeDof(const_cast<DegreeOfFreedom*>(el->getDof(node+2)), EDGE);
244
      }
245
      if (mesh->getNumberOfDofs(EDGE) || mesh->getNumberOfDofs(CENTER))
246
        refineList.removeDOFParents(n_neigh);
247
    }
248

249
    return dof[0];
250
251
252
  }


253
  void RefinementManager2d::bisectTriangle(Triangle *el, DegreeOfFreedom* newDOFs[3])
254
  {
255
    FUNCNAME_DBG("RefinementManager2d::bisectTriangle()");
256

257
258
    TEST_EXIT_DBG(mesh)("No mesh!\n");

259
260
261
    Triangle *child[2];
    child[0] = dynamic_cast<Triangle*>(mesh->createNewElement(el));
    child[1] = dynamic_cast<Triangle*>(mesh->createNewElement(el));
262

263
    int newMark = std::max(0, el->getMark() - 1);
264
265
266
267
    child[0]->setMark(newMark);
    child[1]->setMark(newMark);
    el->setMark(0);

268

269
    // === Transfer hidden data from parent to children. ===
270

271
272
    // call of subclass-method
    el->refineElementData(child[0], child[1]);
273
274
275
    el->setFirstChild(child[0]);
    el->setSecondChild(child[1]);

276
    if (newMark > 0)
277
      doMoreRecursiveRefine = true;
278

279

280
    // === Vertex 2 is the newest vertex. ===
281

282
283
    child[0]->setDof(2, newDOFs[0]);
    child[1]->setDof(2, newDOFs[0]);
284

285

286
    // === The other vertices are handed on from the parent. ===
287

288
    for (int i_child = 0; i_child < 2; i_child++) {
289
290
      child[i_child]->setDof(i_child, const_cast<DegreeOfFreedom*>(el->getDof(2)));
      child[i_child]->setDof(1 - i_child, const_cast<DegreeOfFreedom*>(el->getDof(i_child)));
291
    }
292

293

294
295
    // === There is one more leaf element, two hierachical elements and one ===
    // === more edge.                                                       ===
296
297
298
299
300

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

301
    if (mesh->getNumberOfDofs(EDGE)) {
302
      DegreeOfFreedom* newEdgeDOFs = mesh->getDof(EDGE);
303
304

      // There are dofs in the midpoint of the edges.
305
306
      child[0]->setDof(4, newEdgeDOFs);
      child[1]->setDof(3, newEdgeDOFs);
307

308
      // Dofs handed on by the parent.
309
310
      child[0]->setDof(5, const_cast<DegreeOfFreedom*>(el->getDof(4)));
      child[1]->setDof(5, const_cast<DegreeOfFreedom*>(el->getDof(3)));
311

312
      // Dofs in the refinement edge.
313
314
      child[0]->setDof(3, newDOFs[1]);
      child[1]->setDof(4, newDOFs[2]);
315
    }
316

317
    if (mesh->getNumberOfDofs(CENTER)) {
318
      int node = mesh->getNode(CENTER);
319

320
      // There are dofs at the barycenter of the triangles.
321
322
      child[0]->setDof(node, mesh->getDof(CENTER));
      child[1]->setDof(node, mesh->getDof(CENTER));
323
    }
324
325
  }

326

327
  void RefinementManager2d::getRefinePatch(ElInfo **elInfo,
328
					  DegreeOfFreedom *edge[2],
329
					  int dir, RCNeighbourList &refineList,
330
331
					  int *n_neigh)
  {
332
    FUNCNAME_DBG("RefinementManager2d::getRefinePatch()");
333

334
    if ((*elInfo)->getNeighbour(2) && (*elInfo)->getOppVertex(2) != 2) {
335
336
      // Neighbour is not compatible devisible; refine neighbour first, store the
      // opp_vertex to traverse back to el.
337
      int opp_vertex = (*elInfo)->getOppVertex(2);
338

339
      ElInfo *neigh_info = stack->traverseNeighbour2d(*elInfo, 2);
340
      neigh_info->getElement()->setMark(std::max(neigh_info->getElement()->getMark(), 1));
341
      neigh_info = refineFunction(neigh_info);
342
343
344
345
346
347
348
349
350
351
352

      // === 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())
353
	("Got wrong neighbour! Should be %d, but is %d!\n",
354
355
	 testIndex, (*elInfo)->getElement()->getIndex());

356
      TEST_EXIT_DBG(neigh_info->getElement() ==
357
		    dynamic_cast<Triangle*>(const_cast<Element*>((*elInfo)->getElement())))
358
	("invalid traverse_neighbour1\n");
359
    }
360

361
    if (refineList.setElement(1, (*elInfo)->getNeighbour(2))) {
362
      TEST_EXIT_DBG((*elInfo)->getOppVertex(2) == 2)
363
364
365
366
367
368
	("no compatible ref. edge after recursive refinement of neighbour\n");
      *n_neigh = 2;
    }
  }

}