RCNeighbourList.cc 11.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#include "RCNeighbourList.h"
#include "Element.h"
#include "Line.h"
#include "Triangle.h"
#include "Tetrahedron.h"
#include "MacroElement.h"
#include "Mesh.h"
#include "Traverse.h"
#include "FixVec.h"
#include "CoarseningManager.h"
#include "PeriodicBC.h"
#include "DOFVector.h"
#include "LeafData.h"

namespace AMDiS {

17
18
  RCNeighbourList::RCNeighbourList(int maxEdgeNeigh) 
  { 
19
    rclist.resize(maxEdgeNeigh); 
20
    for (int i = 0; i < maxEdgeNeigh; i++)
21
22
23
      rclist[i] = new RCListElement;
  }

24
25
26
  RCNeighbourList::~RCNeighbourList() 
  {
    for (int i = 0; i < static_cast<int>(rclist.size()); i++)
27
28
29
30
31
32
33
34
35
36
      delete rclist[i];
  }

  /****************************************************************************/
  /*  do_coarse_patch:  if patch can be coarsend return true, else false and  */
  /*  reset the element marks                                                 */
  /****************************************************************************/

  bool RCNeighbourList::doCoarsePatch(int n_neigh)
  {
37
    FUNCNAME("RCNeighbourList::doCoarsePatch()");
38
  
39
40
    for (int i = 0; i < n_neigh; i++) {
      Element *lel = rclist[i]->el;
41

42
      if (lel->getMark() >= 0 || lel->isLeaf()) {
43
44
45
46
47
48
49
50
51
	/****************************************************************************/
	/*  element must not be coarsend or element is a leaf element, reset the    */
	/*  the coarsening flag on all those elements that have to be coarsend with */
	/*  this element                                                            */
	/****************************************************************************/
	lel->setMark(0);
	for (int j = 0; j < n_neigh; j++)
	  if (rclist[j]->flag)  
	    rclist[j]->el->setMark(0);
52

53
	return false;
54
      } else if (lel->getFirstChild()->getMark() >= 0 ||  
55
56
57
58
59
60
61
62
63
64
65
66
		 lel->getSecondChild()->getMark() >= 0) {

	/****************************************************************************/
	/*  one of the element's children must not be coarsend; reset the coarsening*/
	/*  flag on all those elements that have to be coarsend with this element   */
	/****************************************************************************/
	lel->setMark(0);
	for (int j = 0; j < n_neigh; j++)
	  if (rclist[j]->flag)  
	    rclist[j]->el->setMark(0);

	return false;
67
68
      } else if (!lel->getFirstChild()->isLeaf() ||  
		 !lel->getSecondChild()->isLeaf()) {
69
70
71
72
73
74
75
76
77
78
79
80
81
	/****************************************************************************/
	/*  one of the element's children is not a leaf element;                    */
	/*  element may be coarsend after coarsening one of the children; try again */
	/****************************************************************************/
	coarseningManager->doMore = true;

	return false;
      } else {
	/****************************************************************************/ 
	/*  either one element is a macro element or we can coarsen the patch       */
	/****************************************************************************/
	if (rclist[i]->flag == 0) {
	  Mesh* coarse_mesh = coarseningManager->getMesh(); 
82
	  std::deque<MacroElement*>::const_iterator mel;
83
84
	  // set in Mesh::coarsen()
	  for (mel = coarse_mesh->firstMacroElement(); 
85
	       mel!=coarse_mesh->endOfMacroElements(); ++mel)
86
87
88
	    if ((*mel)->getElement() == lel)  
	      break;

89
	  TEST_EXIT(mel != coarse_mesh->endOfMacroElements())
90
91
	    ("incompatible coarsening patch found\n");
	}
92
      }
93
    }
94

95
    return true;
96
97
98
99
100
  }


  void  RCNeighbourList::getNeighOnPatch(int n_neigh, int bound)
  {
101
102
103
    FUNCNAME("RCNeighbourList::getNeighOnPatch()");

    for (int i = 0; i < n_neigh; i++) {
104
105
      Element *el = rclist[i]->el;
      Element *neigh = NULL;
106
107
      rclist[i]->no = i;
      
108
109
110
      for (int dir = 0; dir < 2; dir++) {
	int j = 0;
	for (; j < n_neigh; j++) {
111
112
113
	  if ((neigh = rclist[j]->el) == el)  
	    continue;
	  
114
115
	  int k = 0;
	  for (; k < 2; k++) {
116
117
118
119
120
	    if (neigh->getDOF(2+k) == el->getDOF(3-dir)) {
	      rclist[i]->neigh[dir] = rclist[j];
	      rclist[i]->oppVertex[dir] = 3-k;
	      break;
	    }
121
	  }
122
123
124
125
126
127
128
129
130
	  
	  if (k < 2) 
	    break;
	}

	if (j >= n_neigh) {
	  rclist[i]->neigh[dir] = NULL;
	  rclist[i]->oppVertex[dir] = -1;
	}
131
      }
132
    }
133
134
135
136
137
  }

  void RCNeighbourList::addDOFParent(int elIndex, DegreeOfFreedom* dof)  // 3d
  {
    Element *el = rclist[elIndex]->el;
138
    RCListElement *neighbour = NULL; 
139
140
141
    Mesh *coarse_mesh = coarseningManager->getMesh();
    RCListElement *coarse_list = rclist[elIndex];
  
142
143
    if (coarse_mesh->getNumberOfDOFs(EDGE)) {
      int node = coarse_mesh->getNode(EDGE);
144

145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
      /****************************************************************************/
      /*  set the dof in the coarsening edge                                      */
      /****************************************************************************/
      
      el->setDOF(node, dof);
      
      /****************************************************************************/
      /*  and now those handed on by the children                                 */
      /****************************************************************************/
      
      el->setDOF(node + 1, const_cast<int*>( el->getFirstChild()->getDOF(node)));
      el->setDOF(node + 2, const_cast<int*>( el->getFirstChild()->getDOF(node + 1)));
      el->setDOF(node + 5, const_cast<int*>( el->getFirstChild()->getDOF(node + 3)));
      
      if (coarse_list->elType) {
	el->setDOF(node + 3, const_cast<int*>( el->getSecondChild()->getDOF(node)));
	el->setDOF(node + 4, const_cast<int*>( el->getSecondChild()->getDOF(node + 1)));
      } else {
	el->setDOF(node + 3, const_cast<int*>( el->getSecondChild()->getDOF(node + 1)));
	el->setDOF(node + 4, const_cast<int*>( el->getSecondChild()->getDOF(node)));
165
      }
166
    }
167

168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
    if (coarse_mesh->getNumberOfDOFs(FACE)) {
      int node = coarse_mesh->getNode(FACE);
      /****************************************************************************/
      /*  dof's at the faces within the patch: add new dof's if it is a boundary  */
      /*  face or neighbour is an element behind this element in the corsen list  */
      /****************************************************************************/
      
      neighbour = coarse_list->neigh[0];
      if (!neighbour || neighbour > coarse_list) {
	if (!el->getDOF(node + 2)) {
	  // face 2
	  el->setDOF(node + 2, const_cast<int*>(coarse_mesh->getDOF(FACE)));
	  if (neighbour)
	    neighbour->el->setDOF(node + coarse_list->oppVertex[0], 
				  const_cast<int*>(el->getDOF(node + 2)));
	}
184
      }
185
186
187
188
189
190
191
192
193
194
      
      neighbour = coarse_list->neigh[1];
      if (!neighbour || neighbour > coarse_list) {
	if (!el->getDOF(node + 3)) {
	  // face 3
	  el->setDOF(node + 3, const_cast<int*>(coarse_mesh->getDOF(FACE)));
	  if (neighbour)
	    neighbour->el->setDOF(node + coarse_list->oppVertex[1], 
				  const_cast<int*>(el->getDOF(node + 3)));
	}
195
      }
196
197
198
199
200
201
202
203
204
205
206
207
208
      /****************************************************************************/
      /*  and now those handed on by the children                                 */
      /****************************************************************************/

      el->setDOF(node, const_cast<int*>(el->getSecondChild()->getDOF(node + 3)));
      el->setDOF(node + 1, const_cast<int*>(el->getFirstChild()->getDOF(node + 3)));
    }
    
    if (coarse_mesh->getNumberOfDOFs(CENTER)) {
      int node = coarse_mesh->getNode(CENTER);
      if (!el->getDOF(node))
	el->setDOF(node, const_cast<int*>(coarse_mesh->getDOF(CENTER)));
    }
209
210
211
212
213
214
  }

  void RCNeighbourList::addDOFParents(int n_neigh) // 2d 
  {
    Mesh *coarse_mesh = coarseningManager->getMesh();

215
216
217
218
219
220
221
222
223
    if (coarse_mesh->getNumberOfDOFs(EDGE)) {
      int node = coarse_mesh->getNode(EDGE);
      
      /****************************************************************************/
      /*  get dofs on the boundary of the coarsening patch from the children      */
      /****************************************************************************/
      for (int i = 0; i < n_neigh; i++) {
	rclist[i]->el->setDOF(node, const_cast<int*>(rclist[i]->el->getSecondChild()->getDOF(node + 2)));
	rclist[i]->el->setDOF(node + 1, const_cast<int*>(rclist[i]->el->getFirstChild()->getDOF(node + 2)));
224
      }
225
    }
226

227
228
229
230
231
232
233
234
235
236
    if (coarse_mesh->getNumberOfDOFs(CENTER)) {
      int node = coarse_mesh->getNode(CENTER);
      
      /****************************************************************************/
      /*  get new dof on parents at the barycenter                                */
      /****************************************************************************/
      for (int i = 0; i < n_neigh; i++)
	if (!rclist[i]->el->getDOF(node))
	  rclist[i]->el->setDOF(node, const_cast<int*>(coarse_mesh->getDOF(CENTER)));
    }
237
238
239
240
241
242
243
244
  }


  void RCNeighbourList::removeDOFParents(int n_neigh)
  {
    Mesh *mesh = rclist[0]->el->getMesh();
    int edges = mesh->getGeo(EDGE);

245
246
    if (mesh->getNumberOfDOFs(EDGE)) {      
      int node = mesh->getNode(EDGE);
247

248
249
250
251
	for (int i = 0; i < n_neigh; i++)
	  for (int j = 0; j < edges; j++)
	    rclist[i]->el->setDOF(node + j, NULL);
    }
252

253
254
255
    if (mesh->getNumberOfDOFs(CENTER)) {
      int node = mesh->getNode(CENTER);
      for (int i = 0; i < n_neigh; i++) {
256
	mesh->freeDof(const_cast<int*>(rclist[i]->el->getDOF(node)), CENTER);
257
	rclist[i]->el->setDOF(node, NULL);
258
      }
259
    }
260
261
  }

262

263
264
265
266
267
268
269
270
  void RCNeighbourList::removeDOFParent(int index)
  {
    Tetrahedron *el = dynamic_cast<Tetrahedron*>(rclist[index]->el);
    Mesh *mesh = el->getMesh();
    int edges = mesh->getGeo(EDGE);
    int faces = mesh->getGeo(FACE);


271
272
273
274
275
    if (mesh->getNumberOfDOFs(EDGE)) {
      int node = mesh->getNode(EDGE);
      for (int j = 0; j < edges; j++)
	el->setDOF(node + j, NULL);
    }
276

277
278
279
    if (mesh->getNumberOfDOFs(FACE)) {
      int node = mesh->getNode(FACE);
      RCListElement *neigh = rclist[index]->neigh[0];
280

281
282
      // face 2
      if (!neigh || neigh > rclist[index])
283
	mesh->freeDof(const_cast<int*>(el->getDOF(node + 2)), FACE);
284
285
286
287
      
      neigh = rclist[index]->neigh[1];
      // face 3
      if (!neigh || neigh > rclist[index])
288
	mesh->freeDof(const_cast<int*>(el->getDOF(node + 3)), FACE);
289
290
291
292
      
      for (int j = 0; j < faces; j++)
	el->setDOF(node + j, NULL);
    }
293
  
294
295
    if (mesh->getNumberOfDOFs(CENTER)) {
      int node = mesh->getNode(CENTER);
296
      mesh->freeDof(const_cast<int*>(el->getDOF(node)), CENTER);
297
298
      el->setDOF(node, NULL);
    }
299
300
  }

301

302
303
304
305
306
  RCNeighbourList *RCNeighbourList::periodicSplit(DegreeOfFreedom *edge[2],
						  DegreeOfFreedom *nextEdge[2],
						  int *n_neigh,
						  int *n_neigh_periodic)
  {
307
    RCNeighbourList *periodicList = NULL;
308
    *n_neigh_periodic = 0;
309
310
    int count = 0;
    int n_neigh_old = *n_neigh;
311
312
313
314
315
    bool secondPart = false;

    nextEdge[0] = NULL;
    nextEdge[1] = NULL;

316
317
    std::vector<RCListElement*>::iterator it = rclist.begin();
    std::vector<RCListElement*>::iterator insertIt;
318

319
    while (count < n_neigh_old) {
320
321
322
      DegreeOfFreedom *dof0 = const_cast<DegreeOfFreedom*>((*it)->el->getDOF(0));
      DegreeOfFreedom *dof1 = const_cast<DegreeOfFreedom*>((*it)->el->getDOF(1));

323
      if (dof0 != edge[0] && dof0 != edge[1]) {
324
	secondPart = true;
325
	if (!nextEdge[0]) {
326
327
328
329
330
331
332
333
	  nextEdge[0] = dof0;
	  nextEdge[1] = dof1;
	}
	++it;
      } else {
	(*n_neigh)--;
	(*n_neigh_periodic)++;

334
	if (!periodicList) {
Thomas Witkowski's avatar
Thomas Witkowski committed
335
	  periodicList = new RCNeighbourList;
336
337
338
	  insertIt = periodicList->rclist.end();
	  periodicList->coarseningManager = coarseningManager;
	} else {
339
340
341
	  if (edge[0]) {
	    TEST_EXIT_DBG((dof0 == edge[0] && dof1 == edge[1]) ||
			  (dof1 == edge[0] && dof0 == edge[1]))
342
343
344
345
	      ("invalid macro file?\n");
	  }
	}

346
	if (secondPart) {
347
348
349
350
351
352
353
354
355
	  insertIt = periodicList->rclist.begin();
	  secondPart = false;
	} 

	insertIt = periodicList->rclist.insert(insertIt, *it);
	++insertIt;

	it = rclist.erase(it);
      } 
356
      count++;
357
358
    }

359
    for (int i = 0; i < *n_neigh_periodic; i++)
360
361
362
363
364
365
      periodicList->rclist[i]->no = i;

    return periodicList;
  }

}