CoarseningManager2d.cc 7.94 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
#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)
  {
22
    FUNCNAME("CoarseningManager2d::coarseTriangle()");
23

24
25
26
    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))); 
27

28
    TEST_EXIT_DBG(child[0]->getMark() < 0  &&  child[1]->getMark() < 0)
29
30
31
      ("element %d with children[%d,%d] must not be coarsend!\n",
       el->getIndex(), child[0]->getIndex(), child[1]->getIndex());
  
32
33
34
35
36
37
    if (mesh->getNumberOfDOFs(EDGE)) {
      /****************************************************************************/
      /*  remove dof from common edge of child[0] and child[1]                    */
      /****************************************************************************/
      mesh->freeDOF(const_cast<int*>(child[0]->getDOF(4)), EDGE);
    }
38

39
40
41
42
43
44
45
46
47
48
    if (mesh->getNumberOfDOFs(CENTER)) {
      /****************************************************************************/
      /*  remove dof from the barycenters of child[0] and child[1]                */
      /****************************************************************************/
      int  node = mesh->getNode(CENTER);
      
      mesh->freeDOF(const_cast<int*>(child[0]->getDOF(node)), CENTER);
      mesh->freeDOF(const_cast<int*>(child[1]->getDOF(node)), CENTER);
    }
    
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
    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);

    return;
  }

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

  void CoarseningManager2d::coarsenPatch(RCNeighbourList *coarsenList, 
72
73
					 int n_neigh, 
					 int bound)
74
75
76
77
78
79
  {
    Triangle *el = dynamic_cast<Triangle*>(const_cast<Element*>( coarsenList->getElement(0)));
    Triangle *neigh = dynamic_cast<Triangle*>(const_cast<Element*>( coarsenList->getElement(1)));
    DegreeOfFreedom *dof[3];

    dof[0] = const_cast<int*>( el->getChild(0)->getDOF(2));
80
81
82
83
    if (mesh->getNumberOfDOFs(EDGE)) {
      dof[1] = const_cast<int*>( el->getChild(0)->getDOF(3));
      dof[2] = const_cast<int*>( el->getChild(1)->getDOF(4));
    }
84
85

    int node = mesh->getNode(EDGE);
86
87
88
89
90
91
92
93
94
    if (mesh->getNumberOfDOFs(EDGE)) {
      /****************************************************************************/
      /*  get new dof on el at the midpoint of the coarsening edge                */
      /****************************************************************************/
      if (!el->getDOF(node+2)) {
	el->setDOF(node+2, mesh->getDOF(EDGE));
	if (neigh) {
	  neigh->setDOF(node+2, const_cast<int*>( el->getDOF(node+2)));
	}
95
      }
96
    }
97
98
99
100
101
102
103
104
105

    if (mesh->getNumberOfDOFs(EDGE)  ||  mesh->getNumberOfDOFs(CENTER)) {
      coarsenList->addDOFParents(n_neigh);
    }

    /****************************************************************************/
    /*  restrict dof vectors to the parents on the patch                        */
    /****************************************************************************/
    int nrAdmin = mesh->getNumberOfDOFAdmin();
106
    for (int iadmin = 0; iadmin < nrAdmin; iadmin++) {
107
      std::list<DOFIndexedBase*>::iterator it;
108
      DOFAdmin* admin = const_cast<DOFAdmin*>(&mesh->getDOFAdmin(iadmin));
109
      std::list<DOFIndexedBase*>::iterator end = admin->endDOFIndexed();    
110
      for (it = admin->beginDOFIndexed(); it != end; ++it) {
111
112
113
114
115
116
	(*it)->coarseRestrict(*coarsenList, n_neigh);
      }
    }

    coarsenTriangle(el);

117
118
    if (neigh) 
      coarsenTriangle(neigh);
119
120
121
122
123

    /****************************************************************************/
    /*  now, remove those dofs in the coarcening edge                           */
    /****************************************************************************/
    mesh->freeDOF(dof[0], VERTEX);
124
125
126
127
    if (mesh->getNumberOfDOFs(EDGE)) {
      mesh->freeDOF(dof[1], EDGE);
      mesh->freeDOF(dof[2], EDGE);
    }
128
129
130
131
132
133
134
135
136

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

    return;
  }

  int CoarseningManager2d::coarsenFunction(ElInfo *el_info)
  {
137
    Triangle *el = dynamic_cast<Triangle*>(const_cast<Element*>( el_info->getElement()));
138
    DegreeOfFreedom *edge[2];
139
    int n_neigh, bound = 0;
140
141
142
143
    RCNeighbourList coarse_list(2);

    coarse_list.setCoarseningManager(this);

144
145
146
147
148
149
150
151
152
153
154
155
    if (el->getMark() >= 0)
      return 0; // el must not be coarsend, return
    if (!(el->getChild(0))) 
      return 0;  // single leaves don't get coarsened

    if (el->getChild(0)->getMark() >= 0  || el->getChild(1)->getMark() >= 0) {
      /****************************************************************************/
      /*  one of the children must not be coarsend; return :-(                    */
      /****************************************************************************/
      el->setMark(0);
      return 0;
    }
156

157
158
159
160
161
162
163
    if (!el->getChild(0)->isLeaf() || !el->getChild(1)->isLeaf()) {
      /****************************************************************************/
      /*  one of the children is not a leaf element; try again later on           */
      /****************************************************************************/
      doMore = true;
      return 0;
    }
164
165
166
167
168

    /****************************************************************************/
    /*  give the refinement edge the right orientation                          */
    /****************************************************************************/

169
170
171
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));
    } else {
      edge[1] = const_cast<int*>( el->getDOF(0));
      edge[0] = const_cast<int*>( el->getDOF(1));
    }
176
177
178
179

    coarse_list.setElement(0, el, true);

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

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

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

198
      while (edge[0] != NULL) {
199
200
201
202
203
	periodicList = coarse_list.periodicSplit(edge, 
						 next_edge,
						 &n_neigh,
						 &n_neigh_periodic);

204
	TEST_EXIT_DBG(periodicList)("periodicList = NULL\n");
205
206
207
208
209
210
211
212
213
214
215
216

	coarsenPatch(periodicList, n_neigh_periodic, bound);

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

    return 0;
  }

}