RCNeighbourList.cc 12.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
#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 {

29
30
  RCNeighbourList::RCNeighbourList(int maxEdgeNeigh) 
  { 
31
    rclist.resize(maxEdgeNeigh); 
32
    for (int i = 0; i < maxEdgeNeigh; i++)
33
34
35
      rclist[i] = new RCListElement;
  }

36

37
38
  RCNeighbourList::~RCNeighbourList() 
  {
39
    clearList();
40
41
  }

42

43
44
45
46
47
48
49
  /****************************************************************************/
  /*  do_coarse_patch:  if patch can be coarsend return true, else false and  */
  /*  reset the element marks                                                 */
  /****************************************************************************/

  bool RCNeighbourList::doCoarsePatch(int n_neigh)
  {
50
    FUNCNAME("RCNeighbourList::doCoarsePatch()");
51
  
52
53
    for (int i = 0; i < n_neigh; i++) {
      Element *lel = rclist[i]->el;
54

55
      if (lel->getMark() >= 0 || lel->isLeaf()) {
56
57
58
59
60
61
62
63
64
	/****************************************************************************/
	/*  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);
65

66
	return false;
67
      } else if (lel->getFirstChild()->getMark() >= 0 ||  
68
69
70
71
72
73
74
75
76
77
78
79
		 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;
80
81
      } else if (!lel->getFirstChild()->isLeaf() ||  
		 !lel->getSecondChild()->isLeaf()) {
82
83
84
85
86
87
88
89
90
91
92
93
94
	/****************************************************************************/
	/*  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(); 
95
	  std::deque<MacroElement*>::const_iterator mel;
96
97
	  // set in Mesh::coarsen()
	  for (mel = coarse_mesh->firstMacroElement(); 
98
	       mel!=coarse_mesh->endOfMacroElements(); ++mel)
99
100
101
	    if ((*mel)->getElement() == lel)  
	      break;

102
	  TEST_EXIT(mel != coarse_mesh->endOfMacroElements())
103
104
	    ("incompatible coarsening patch found\n");
	}
105
      }
106
    }
107

108
    return true;
109
110
111
  }


112
  void  RCNeighbourList::fillNeighbourRelations(int n_neigh, int bound)
113
  {
114
115
    FUNCNAME("RCNeighbourList::getNeighOnPatch()");

116
	// for all RC Elements in list
117
    for (int i = 0; i < n_neigh; i++) {
118
      Element *el = rclist[i]->el;
119
      rclist[i]->ith = i;
120
      
121
	  // for all directions of refinementEdges =? Sides of Refinement Edge
122
123
      for (int dir = 0; dir < 2; dir++) {
	int j = 0;
124
125

	//Test if other elements in List are Neighbours of active Element
126
	for (; j < n_neigh; j++) {
127
	  Element *neigh = rclist[j]->el;
128
129

	  //test-neighbour Element is active element -> skip tests for this elements
130
	  if (neigh == el)  
131
132
	    continue;
	  
133
	  //Test if test-neighbour Element is neighbour of active Element
134
135
	  int k = 0;
	  for (; k < 2; k++) {
136
	    if (neigh->getDof(2 + k) == el->getDof(3 - dir)) {
137
138
139
140
				//Test FACE neighbourhood via oppVertex
				// Test only Neighbourhood of Faces 2 & 3(local Indexing)
				//			Faces next to refinement Edge
				// check works by comparing DOFs (rank Indexing)
141
	      rclist[i]->neigh[dir] = rclist[j];
142
	      rclist[i]->oppVertex[dir] = 3 - k;
143
	      break;
144
145
146
	    } else {
	      rclist[i]->neigh[dir] = NULL;
	      rclist[i]->oppVertex[dir] = -1;				
147
	    }
148
	  }
149
150
151
152
153
	  
	  if (k < 2) 
	    break;
	}

154
	//none of the Elements in rcList is FACE neighbour of active Element
155
156
157
158
	if (j >= n_neigh) {
	  rclist[i]->neigh[dir] = NULL;
	  rclist[i]->oppVertex[dir] = -1;
	}
159
      }
160
    }
161
162
  }

163

164
165
166
  void RCNeighbourList::addDOFParent(int elIndex, DegreeOfFreedom* dof)  // 3d
  {
    Element *el = rclist[elIndex]->el;
167
    RCListElement *neighbour = NULL; 
168
169
170
    Mesh *coarse_mesh = coarseningManager->getMesh();
    RCListElement *coarse_list = rclist[elIndex];
  
171
    if (coarse_mesh->getNumberOfDofs(EDGE)) {
172
      int node = coarse_mesh->getNode(EDGE);
173

174
175
176
177
      /****************************************************************************/
      /*  set the dof in the coarsening edge                                      */
      /****************************************************************************/
      
178
      el->setDof(node, dof);
179
180
181
182
183
      
      /****************************************************************************/
      /*  and now those handed on by the children                                 */
      /****************************************************************************/
      
184
185
186
      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)));
187
188
      
      if (coarse_list->elType) {
189
190
	el->setDof(node + 3, const_cast<int*>(el->getSecondChild()->getDof(node)));
	el->setDof(node + 4, const_cast<int*>(el->getSecondChild()->getDof(node + 1)));
191
      } else {
192
193
	el->setDof(node + 3, const_cast<int*>(el->getSecondChild()->getDof(node + 1)));
	el->setDof(node + 4, const_cast<int*>(el->getSecondChild()->getDof(node)));
194
      }
195
    }
196

197
    if (coarse_mesh->getNumberOfDofs(FACE)) {
198
199
200
201
202
203
204
205
      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) {
206
	if (!el->getDof(node + 2)) {
207
	  // face 2
208
	  el->setDof(node + 2, const_cast<int*>(coarse_mesh->getDof(FACE)));
209
	  if (neighbour)
210
211
	    neighbour->el->setDof(node + coarse_list->oppVertex[0], 
				  const_cast<int*>(el->getDof(node + 2)));
212
	}
213
      }
214
215
216
      
      neighbour = coarse_list->neigh[1];
      if (!neighbour || neighbour > coarse_list) {
217
	if (!el->getDof(node + 3)) {
218
	  // face 3
219
	  el->setDof(node + 3, const_cast<int*>(coarse_mesh->getDof(FACE)));
220
	  if (neighbour)
221
222
	    neighbour->el->setDof(node + coarse_list->oppVertex[1], 
				  const_cast<int*>(el->getDof(node + 3)));
223
	}
224
      }
225
226
227
228
      /****************************************************************************/
      /*  and now those handed on by the children                                 */
      /****************************************************************************/

229
230
      el->setDof(node, const_cast<int*>(el->getSecondChild()->getDof(node + 3)));
      el->setDof(node + 1, const_cast<int*>(el->getFirstChild()->getDof(node + 3)));
231
232
    }
    
233
    if (coarse_mesh->getNumberOfDofs(CENTER)) {
234
      int node = coarse_mesh->getNode(CENTER);
235
236
      if (!el->getDof(node))
	el->setDof(node, const_cast<int*>(coarse_mesh->getDof(CENTER)));
237
    }
238
239
  }

240

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

245
    if (coarse_mesh->getNumberOfDofs(EDGE)) {
246
247
248
249
250
251
      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++) {
252
253
	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)));
254
      }
255
    }
256

257
    if (coarse_mesh->getNumberOfDofs(CENTER)) {
258
259
260
261
262
263
      int node = coarse_mesh->getNode(CENTER);
      
      /****************************************************************************/
      /*  get new dof on parents at the barycenter                                */
      /****************************************************************************/
      for (int i = 0; i < n_neigh; i++)
264
265
	if (!rclist[i]->el->getDof(node))
	  rclist[i]->el->setDof(node, const_cast<int*>(coarse_mesh->getDof(CENTER)));
266
    }
267
268
269
270
271
272
273
274
  }


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

275
    if (mesh->getNumberOfDofs(EDGE)) {      
276
      int node = mesh->getNode(EDGE);
277

278
279
	for (int i = 0; i < n_neigh; i++)
	  for (int j = 0; j < edges; j++)
280
	    rclist[i]->el->setDof(node + j, NULL);
281
    }
282

283
    if (mesh->getNumberOfDofs(CENTER)) {
284
285
      int node = mesh->getNode(CENTER);
      for (int i = 0; i < n_neigh; i++) {
286
287
	mesh->freeDof(const_cast<int*>(rclist[i]->el->getDof(node)), CENTER);
	rclist[i]->el->setDof(node, NULL);
288
      }
289
    }
290
291
  }

292

293
294
295
296
297
298
299
300
  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);


301
    if (mesh->getNumberOfDofs(EDGE)) {
302
303
      int node = mesh->getNode(EDGE);
      for (int j = 0; j < edges; j++)
304
	el->setDof(node + j, NULL);
305
    }
306

307
    if (mesh->getNumberOfDofs(FACE)) {
308
309
      int node = mesh->getNode(FACE);
      RCListElement *neigh = rclist[index]->neigh[0];
310

311
312
      // face 2
      if (!neigh || neigh > rclist[index])
313
	mesh->freeDof(const_cast<int*>(el->getDof(node + 2)), FACE);
314
315
316
317
      
      neigh = rclist[index]->neigh[1];
      // face 3
      if (!neigh || neigh > rclist[index])
318
	mesh->freeDof(const_cast<int*>(el->getDof(node + 3)), FACE);
319
320
      
      for (int j = 0; j < faces; j++)
321
	el->setDof(node + j, NULL);
322
    }
323
  
324
    if (mesh->getNumberOfDofs(CENTER)) {
325
      int node = mesh->getNode(CENTER);
326
327
      mesh->freeDof(const_cast<int*>(el->getDof(node)), CENTER);
      el->setDof(node, NULL);
328
    }
329
330
  }

331

332
333
334
335
336
  void RCNeighbourList::periodicSplit(DegreeOfFreedom *edge[2],
				      DegreeOfFreedom *nextEdge[2],
				      int *n_neigh,
				      int *n_neigh_periodic,
				      RCNeighbourList &periodicList)
337
  {
338
339
    FUNCNAME("RCNeighbourList::periodicSplit()");

340
    *n_neigh_periodic = 0;
341
342
    int count = 0;
    int n_neigh_old = *n_neigh;
343
    bool secondPart = false;
344
    bool firstSplit = true;
345
346
347
348

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

349
350
    std::vector<RCListElement*>::iterator it = rclist.begin();
    std::vector<RCListElement*>::iterator insertIt;
351

352
    while (count < n_neigh_old) {
353
354
      DegreeOfFreedom *dof0 = const_cast<DegreeOfFreedom*>((*it)->el->getDof(0));
      DegreeOfFreedom *dof1 = const_cast<DegreeOfFreedom*>((*it)->el->getDof(1));
355

356
      if (dof0 != edge[0] && dof0 != edge[1]) {
357
	secondPart = true;
358
	if (!nextEdge[0]) {
359
360
361
362
363
364
365
366
	  nextEdge[0] = dof0;
	  nextEdge[1] = dof1;
	}
	++it;
      } else {
	(*n_neigh)--;
	(*n_neigh_periodic)++;

367
368
369
370
371
	if (firstSplit) {
	  periodicList.clearList();
	  insertIt = periodicList.rclist.end();
	  periodicList.coarseningManager = coarseningManager;
	  firstSplit = false;
372
	} else {
373
374
375
	  if (edge[0]) {
	    TEST_EXIT_DBG((dof0 == edge[0] && dof1 == edge[1]) ||
			  (dof1 == edge[0] && dof0 == edge[1]))
376
	      ("Invalid macro file?\n");
377
378
379
	  }
	}

380
	if (secondPart) {
381
	  insertIt = periodicList.rclist.begin();
382
383
384
	  secondPart = false;
	} 

385
	insertIt = periodicList.rclist.insert(insertIt, *it);
386
387
388
389
	++insertIt;

	it = rclist.erase(it);
      } 
390
      count++;
391
392
    }

393
    for (int i = 0; i < *n_neigh_periodic; i++)
394
395
396
397
398
399
400
401
      periodicList.rclist[i]->ith = i;
  }


  void RCNeighbourList::clearList()
  {    
    for (unsigned int i = 0; i < rclist.size(); i++)
      delete rclist[i];
402

403
    rclist.clear();
404
405
406
  }

}