CoarseningManager2d.cc 6.83 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * 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.
 * 
 ******************************************************************************/
20
21


22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
#include "CoarseningManager2d.h"
#include "Mesh.h"
#include "AdaptStationary.h"
#include "AdaptInstationary.h"
#include "Traverse.h"
#include "MacroElement.h"
#include "RCNeighbourList.h"
#include "FixVec.h"
#include "DOFIndexed.h"

namespace AMDiS {

  /****************************************************************************/
  /*  coarseTriangle:  coarses a single element of the coarsening patch; dofs */
  /*  in the interior of the element are removed; dofs for higher order       */
  /*  at the boundary or the coarsening patch still belong to                 */
  /*  the parent. Do not remove them form the mesh!!!                         */
  /****************************************************************************/

  void CoarseningManager2d::coarsenTriangle(Triangle *el)
  {
43
    FUNCNAME_DBG("CoarseningManager2d::coarseTriangle()");
44

45
46
47
    Triangle *child[2];
    child[0] = dynamic_cast<Triangle*>(const_cast<Element*>(el->getChild(0)));
    child[1] = dynamic_cast<Triangle*>(const_cast<Element*>(el->getChild(1))); 
48

49
    TEST_EXIT_DBG(child[0]->getMark() < 0  &&  child[1]->getMark() < 0)
50
51
      ("element %d with children[%d,%d] must not be coarsend!\n",
       el->getIndex(), child[0]->getIndex(), child[1]->getIndex());
52
53

    // remove dof from common edge of child[0] and child[1]  
54
    if (mesh->getNumberOfDofs(EDGE))
55
      mesh->freeDof(const_cast<int*>(child[0]->getDof(4)), EDGE);
56

57
    // remove dof from the barycenters of child[0] and child[1] 
58
    if (mesh->getNumberOfDofs(CENTER)) {
59
      int node = mesh->getNode(CENTER);
60
      
61
62
      mesh->freeDof(const_cast<int*>(child[0]->getDof(node)), CENTER);
      mesh->freeDof(const_cast<int*>(child[1]->getDof(node)), CENTER);
63
64
    }
    
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
    el->coarsenElementData(child[0], child[1]);

    el->setFirstChild(NULL);
    el->setSecondChild(NULL);

    mesh->freeElement(child[0]);
    mesh->freeElement(child[1]);

    el->incrementMark();

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

  /****************************************************************************/
  /*  coarsenPatch: first rebuild the dofs on the parents then do restriction */
  /*  of data (if possible) and finally coarsen the patch elements            */
  /****************************************************************************/

85
  void CoarseningManager2d::coarsenPatch(RCNeighbourList &coarsenList, 
86
87
					 int n_neigh, 
					 int bound)
88
  {
89
    Triangle *el = 
90
      dynamic_cast<Triangle*>(const_cast<Element*>(coarsenList.getElement(0)));
91
    Triangle *neigh = 
92
      dynamic_cast<Triangle*>(const_cast<Element*>(coarsenList.getElement(1)));
93
94
    DegreeOfFreedom *dof[3];

95
    dof[0] = const_cast<int*>(el->getChild(0)->getDof(2));
96
    if (mesh->getNumberOfDofs(EDGE)) {
97
98
      dof[1] = const_cast<int*>(el->getChild(0)->getDof(3));
      dof[2] = const_cast<int*>(el->getChild(1)->getDof(4));
99
100
    } else {
      dof[1] = dof[2] = 0;
101
    }
102

103
    if (mesh->getNumberOfDofs(EDGE)) {
104
105
106
      int node = mesh->getNode(EDGE);
      // get new dof on el at the midpoint of the coarsening edge

107
108
      if (!el->getDof(node + 2)) {
	el->setDof(node + 2, mesh->getDof(EDGE));
109
	if (neigh)
110
	  neigh->setDof(node + 2, const_cast<int*>(el->getDof(node + 2)));
111
      }
112
    }
113

114
    if (mesh->getNumberOfDofs(EDGE) || mesh->getNumberOfDofs(CENTER))
115
      coarsenList.addDOFParents(n_neigh);    
116
117

    // restrict dof vectors to the parents on the patch 
118
119

    int nrAdmin = mesh->getNumberOfDOFAdmin();
120
    for (int iadmin = 0; iadmin < nrAdmin; iadmin++) {
121
      DOFAdmin* admin = const_cast<DOFAdmin*>(&mesh->getDofAdmin(iadmin));
122
123
      for (std::list<DOFIndexedBase*>::iterator it = admin->beginDOFIndexed(); 
	   it != admin->endDOFIndexed(); ++it)	
124
	(*it)->coarseRestrict(coarsenList, n_neigh);            
125
126
127
128
    }

    coarsenTriangle(el);

129
130
    if (neigh) 
      coarsenTriangle(neigh);
131

132
133
    // now, remove those dofs in the coarcening edge 

134
    mesh->freeDof(dof[0], VERTEX);
135
    if (mesh->getNumberOfDofs(EDGE)) {
136
137
      mesh->freeDof(dof[1], EDGE);
      mesh->freeDof(dof[2], EDGE);
138
    }
139
140
141
142
143

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

144

145
  void CoarseningManager2d::coarsenFunction(ElInfo *el_info)
146
  {
147
    Triangle *el = dynamic_cast<Triangle*>(const_cast<Element*>(el_info->getElement()));
148
    DegreeOfFreedom *edge[2];
149
    int n_neigh, bound = 0;
150
151
152
153
    RCNeighbourList coarse_list(2);

    coarse_list.setCoarseningManager(this);

154
    if (el->getMark() >= 0)
155
      return; // el must not be coarsend, return
156
    if (!(el->getChild(0))) 
157
      return;  // single leaves don't get coarsened
158
159

    if (el->getChild(0)->getMark() >= 0  || el->getChild(1)->getMark() >= 0) {
160
      // one of the children must not be coarsend; return :( 
161
      el->setMark(0);
162
      return;
163
    }
164

165
    if (!el->getChild(0)->isLeaf() || !el->getChild(1)->isLeaf()) {
166
      // one of the children is not a leaf element; try again later on 
167
      doMore = true;
168
      return;
169
    }
170

171
    // give the refinement edge the right orientation
172

173
174
175
    if (el->getDof(0,0) < el->getDof(1,0)) {
      edge[0] = const_cast<int*>(el->getDof(0));
      edge[1] = const_cast<int*>(el->getDof(1));
176
    } else {
177
178
      edge[1] = const_cast<int*>(el->getDof(0));
      edge[0] = const_cast<int*>(el->getDof(1));
179
    }
180
181
182
183

    coarse_list.setElement(0, el, true);

    n_neigh = 1;
184
185
186
187
    if (coarse_list.setElement(1, el_info->getNeighbour(2))) {
      n_neigh = 2;
      coarse_list.setCoarsePatch(1, el_info->getOppVertex(2) == 2);
    }
188
  
189
    // check wether we can coarsen the patch or not                            
190
191
192
193
194

    // ==========================================================================
    // === check for periodic boundary ==========================================
    // ==========================================================================

195
    if (coarse_list.doCoarsePatch(n_neigh)) {
196
197
      int n_neigh_periodic;
      DegreeOfFreedom *next_edge[2];
198
      RCNeighbourList periodicList;
199

200
      while (edge[0] != NULL) {
201
202
203
	coarse_list.periodicSplit(edge, next_edge,
				  &n_neigh, &n_neigh_periodic,
				  periodicList);
204
205
206
207
208
209
210
211
212
213

	coarsenPatch(periodicList, n_neigh_periodic, bound);

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

}