Tetrahedron.cc 14.1 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
#include "Tetrahedron.h"
#include "DOFAdmin.h"
#include "Mesh.h"
#include "CoarseningManager.h"
#include "FixVec.h"
18
#include "ElementDofIterator.h"
19
#include "Global.h"
20
21
22
23
24
25

namespace AMDiS {

  const unsigned char Tetrahedron::nChildEdge[3][2][2] = {{{5,4},{4,5}},
							  {{5,4},{5,4}},
							  {{5,4},{5,4}}};
26

27
28
29
  const unsigned char Tetrahedron::nChildFace[3][2][2] = {{{1,2},{2,1}},
							  {{1,2},{1,2}},
							  {{1,2},{1,2}}};
30

31
32
33
  const int Tetrahedron::childVertex[3][2][4] = {{{0,2,3,4},{1,3,2,4}},
						 {{0,2,3,4},{1,2,3,4}},
						 {{0,2,3,4},{1,2,3,4}}};
34

35
36
37
  const int Tetrahedron::childEdge[3][2][6] = {{{1,2,0,5,5,4}, {4,3,0,5,5,4}},
					       {{1,2,0,5,4,5}, {3,4,0,5,4,5}},
					       {{1,2,0,5,4,5}, {3,4,0,5,4,5}}};
38

39
40
41
42
43
44
  const unsigned char Tetrahedron::adjacentChild[2][2] = {{0,1}, {1,0}};

  const signed char Tetrahedron::childOrientation[3][2] = {{1,1}, 
							   {1,-1}, 
							   {1,-1}};

45
46
47
48
  const unsigned char Tetrahedron::edgeOfDofs[4][4] =  {{255, 0, 1, 2},
							{0, 255, 3, 4},
							{1, 3, 255, 5},
							{2, 4, 5, 255}};
49
50
51
52
53
54
55

  const int Tetrahedron::vertexOfEdge[6][2] = {{0,1},
					       {0,2},
					       {0,3},
					       {1,2},
					       {1,3},
					       {2,3}};
56

57
58
59
60
61
62
63
64
65
66
67
68
  const int Tetrahedron::vertexOfFace[4][3] = {{1,2,3},
					       {0,2,3},
					       {0,1,3},
					       {0,1,2}};

  const int Tetrahedron::sideOfChild[3][2][4] = {{{-1, 3, 1, 2},  // type 0
						  {3, -1, 2, 1}},
						 {{-1, 3, 1, 2},  // type 1
						  {3, -1, 1, 2}},
						 {{-1, 3, 1, 2},  // type 2
						  {3, -1, 1, 2}}};

69
70
71
72
73
74
75
  const int Tetrahedron::edgeOfChild[3][2][6] = {{{2, 0, 1, -1, -1, 3},   // type 0
						  {2, -1, -1, 1, 0, 3}},
						 {{2, 0, 1, -1, -1, 3},   // type 1
						  {2, -1, -1, 0, 1, 3}},
						 {{2, 0, 1, -1, -1, 3},   // type 2
						  {2, -1, -1, 0, 1, 3}}};

76
77
78
79
80
81
82
83
84
85
86
87
  const int Tetrahedron::vertexOfParent[3][2][4] = {{{0, 2, 3, -1},  // type 0
						     {1, 3, 2, -1}},
						    {{0, 2, 3, -1},  // type 1
						     {1, 2, 3, -1}},
						    {{0, 2, 3, -1},  // type 2
						     {1, 2, 3, -1}}};

  const int Tetrahedron::edgeOfFace[4][3] = {{5, 4, 3},  // face 0
					     {5, 2, 1},  // face 1
					     {4, 2, 0},  // face 2
					     {3, 1, 0}}; // face 3

88

89
90
  bool Tetrahedron::hasSide(Element* sideElem) const
  {
91
92
    FUNCNAME("Tetrahedron::hasSide()");
    TEST_EXIT_DBG(sideElem->isTriangle())("called for sideElem-type != Triangle\n");
93
94
95
96
    ERROR_EXIT("not yet\n");
    return false;
  }

97

98
  int Tetrahedron::getVertexOfPosition(GeoIndex position,
99
100
				       int positionIndex,
				       int vertexIndex) const 
101
  {
102
    FUNCNAME("Triangle::getVertexOfPosition()");
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
    switch(position) {
    case VERTEX:
      return positionIndex;
      break;
    case EDGE:
      return vertexOfEdge[positionIndex][vertexIndex];
      break;
    case FACE:
      return vertexOfFace[positionIndex][vertexIndex];
      break;
    case CENTER:
      return vertexIndex;
      break;
    default:
      ERROR_EXIT("invalid position\n");
      return 0;
    }
  }

122

123
124
  const FixVec<int, WORLD>& Tetrahedron::sortFaceIndices(int face, 
							 FixVec<int,WORLD> *vec) const
125
126
127
  {
    static MatrixOfFixVecs<FixVec<int,WORLD> > *sorted_3d = NULL;
    FixVec<int,WORLD> *val = NULL;
128
    const int *vof = vertexOfFace[face];
129

130
131
    if (NULL == sorted_3d) {
      sorted_3d = new MatrixOfFixVecs<FixVec<int,WORLD> >(3, 4, 7, NO_INIT);
132
 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
      (*sorted_3d)[0][0][0] = (*sorted_3d)[0][0][1] =
	(*sorted_3d)[0][0][2] = (*sorted_3d)[1][0][0] =
	(*sorted_3d)[1][0][1] = (*sorted_3d)[1][0][2] =
	(*sorted_3d)[1][1][0] = (*sorted_3d)[1][2][1] =
	(*sorted_3d)[1][3][0] = (*sorted_3d)[1][4][2] =
	(*sorted_3d)[1][5][1] = (*sorted_3d)[1][6][2] =
	(*sorted_3d)[2][0][0] = (*sorted_3d)[2][0][1] =
	(*sorted_3d)[2][0][2] = (*sorted_3d)[2][1][0] =
	(*sorted_3d)[2][2][1] = (*sorted_3d)[2][3][0] =
	(*sorted_3d)[2][4][2] = (*sorted_3d)[2][5][1] =
	(*sorted_3d)[2][6][2] = (*sorted_3d)[3][0][0] =
	(*sorted_3d)[3][0][1] = (*sorted_3d)[3][0][2] =
	(*sorted_3d)[3][1][0] = (*sorted_3d)[3][2][1] =
	(*sorted_3d)[3][3][0] = (*sorted_3d)[3][4][2] =
	(*sorted_3d)[3][5][1] = (*sorted_3d)[3][6][2] = 0;

      (*sorted_3d)[0][1][0] = (*sorted_3d)[0][2][1] =
	(*sorted_3d)[0][3][0] = (*sorted_3d)[0][4][2] =
	(*sorted_3d)[0][5][1] = (*sorted_3d)[0][6][2] =
	(*sorted_3d)[2][1][2] = (*sorted_3d)[2][2][0] =
	(*sorted_3d)[2][3][1] = (*sorted_3d)[2][4][1] =
	(*sorted_3d)[2][5][2] = (*sorted_3d)[2][6][0] =
	(*sorted_3d)[3][1][2] = (*sorted_3d)[3][2][0] =
	(*sorted_3d)[3][3][1] = (*sorted_3d)[3][4][1] =
	(*sorted_3d)[3][5][2] = (*sorted_3d)[3][6][0] = 1;

      (*sorted_3d)[0][1][2] = (*sorted_3d)[0][2][0] =
	(*sorted_3d)[0][3][1] = (*sorted_3d)[0][4][1] =
	(*sorted_3d)[0][5][2] = (*sorted_3d)[0][6][0] =
	(*sorted_3d)[1][1][2] = (*sorted_3d)[1][2][0] =
	(*sorted_3d)[1][3][1] = (*sorted_3d)[1][4][1] =
	(*sorted_3d)[1][5][2] = (*sorted_3d)[1][6][0] =
	(*sorted_3d)[3][1][1] = (*sorted_3d)[3][2][2] =
	(*sorted_3d)[3][3][2] = (*sorted_3d)[3][4][0] =
	(*sorted_3d)[3][5][0] = (*sorted_3d)[3][6][1] = 2;

      (*sorted_3d)[0][1][1] = (*sorted_3d)[0][2][2] =
	(*sorted_3d)[0][3][2] = (*sorted_3d)[0][4][0] =
	(*sorted_3d)[0][5][0] = (*sorted_3d)[0][6][1] =
	(*sorted_3d)[1][1][1] = (*sorted_3d)[1][2][2] =
	(*sorted_3d)[1][3][2] = (*sorted_3d)[1][4][0] =
	(*sorted_3d)[1][5][0] = (*sorted_3d)[1][6][1] =
	(*sorted_3d)[2][1][1] = (*sorted_3d)[2][2][2] =
	(*sorted_3d)[2][3][2] = (*sorted_3d)[2][4][0] =
	(*sorted_3d)[2][5][0] = (*sorted_3d)[2][6][1] = 3;
178
179
    }

180
    int no = 0;
181
182
183
184
185
186
187
    if (dof[vof[0]][0] < dof[vof[1]][0])
      no++;
    if (dof[vof[1]][0] < dof[vof[2]][0])
      no += 2;
    if (dof[vof[2]][0] < dof[vof[0]][0])
      no += 4;

188
    if (!(no >= 1  &&  no <= 6))
189
      ERROR_EXIT("Cannot sort face indices of element %d at face %d\n", 
190
191
192
193
194
		 getIndex(), face);

    if (vec) {
      val = vec;
      *val = (*sorted_3d)[face][no];
195
    } else {
196
      val = &((*sorted_3d)[face][no]);
197
    }
198
199
200
201

    return(*(const_cast<const FixVec<int,WORLD>* >(val)));
  }

202

203
  void Tetrahedron::getNodeDofs(const FiniteElemSpace* feSpace, 
204
				BoundaryObject bound,
205
206
				DofContainer& dofs,
				bool baseDofPtr) const
207
  {
208
    FUNCNAME("Tetrahedron::getNodeDofs()");
209

210
    switch (bound.subObj) {
Thomas Witkowski's avatar
Thomas Witkowski committed
211
    case VERTEX:
212
      {
213
214
215
216
217
218
	if (baseDofPtr) {
	  dofs.push_back(dof[bound.ithObj]);
	} else {
	  int n0 = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
	  dofs.push_back(&(dof[bound.ithObj][n0]));
	}
219
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
220
      break;
221
    case EDGE:
222
      getNodeDofsAtEdge(feSpace, bound, dofs, baseDofPtr);
223
224
      break;
    case FACE:
225
      getNodeDofsAtFace(feSpace, bound, dofs, baseDofPtr);
226
227
228
229
230
231
232
      break;      
    default:
      ERROR_EXIT("Should not happen!\n");
    }
  }


233
  void Tetrahedron::getNodeDofsAtFace(const FiniteElemSpace* feSpace, 
234
				      BoundaryObject bound,
235
236
				      DofContainer& dofs,
				      bool baseDofPtr) const
237
  {
238
    FUNCNAME("Tetrahedron::getNodeDofsAtFace()");
239

240
241
242
243
    if (!child[0])
      return;

    switch (bound.ithObj) {
244
245
    case 0:
      {
246
247
	BoundaryObject nextBound = bound;
	prepareNextBound(nextBound, 1);
248
	child[1]->getNodeDofs(feSpace, nextBound, dofs, baseDofPtr);
249
250
251
252
      }
      break;
    case 1:
      {
253
254
	BoundaryObject nextBound = bound;
	prepareNextBound(nextBound, 0);
255
	child[0]->getNodeDofs(feSpace, nextBound, dofs, baseDofPtr);
256
257
258
259
260
261
      }
      break;

    case 2:
    case 3:
      {
262
263
	int n0 = 
	  (baseDofPtr ? 0 : feSpace->getAdmin()->getNumberOfPreDofs(FACE));
264
265
266
267
268
269
270
	BoundaryObject nextBound0 = bound, nextBound1 = bound;
	prepareNextBound(nextBound0, 0);
	prepareNextBound(nextBound1, 1);
	nextBound0.reverseMode = false;
	nextBound1.reverseMode = false;
	
	bool addDof = true;
271
272
	for (ExcludeList::iterator it = bound.excludedSubstructures.begin(); 
	     it != bound.excludedSubstructures.end(); ++it)
273
	  if (it->first == EDGE && it->second == 0)
274
	    addDof = false;	
275
276
	
	if (bound.reverseMode) {
277
	  child[1]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr);
278
279
	  
	  if (addDof)
280
	    dofs.push_back(&(child[0]->getDof()[3][n0]));
281
	  
282
	  child[0]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr);
283
	} else {
284
	  child[0]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr);
285
286
	  
	  if (addDof)
287
	    dofs.push_back(&(child[0]->getDof()[3][n0]));
288
	  
289
	  child[1]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr);
290
291
292
293
	}
      }
      break;
    default:
294
      ERROR_EXIT("Should never happen: %d\n", bound.ithObj);
295
    }
296
297
  }

298

299
  void Tetrahedron::getNodeDofsAtEdge(const FiniteElemSpace* feSpace, 
300
				      BoundaryObject bound,
301
302
				      DofContainer& dofs,
				      bool baseDofPtr) const
303
  {
304
    FUNCNAME("Tetrahedron::getNodeDofsAtEdge()");
305
306

    if (!child[0])
307
      return;    
308

309
    int n0 = (baseDofPtr ? 0 : feSpace->getAdmin()->getNumberOfPreDofs(EDGE));
310
311
312
313
    BoundaryObject nextBound0 = bound, nextBound1 = bound;
    prepareNextBound(nextBound0, 0);
    prepareNextBound(nextBound1, 1);

314
315
316
317
318
319
    switch (bound.ithObj) {
    case 0:      
      nextBound0.reverseMode = false;
      nextBound1.reverseMode = false;

      if (bound.reverseMode) {
320
	child[1]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr);
321
	dofs.push_back(&(child[0]->getDof()[3][n0]));
322
	child[0]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr);
323
      } else {
324
	child[0]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr);
325
	dofs.push_back(&(child[0]->getDof()[3][n0]));
326
	child[1]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr);
327
328
329
330
      }

      break;
    case 5:
331
332
      TEST_EXIT_DBG(nextBound0.ithObj == nextBound1.ithObj)
	("Should not happen!\n");
333
334

      if (nextBound0.ithObj != -1)
335
	child[0]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr);
336
337
338
339
340
341
      break;
    default:
      TEST_EXIT_DBG(nextBound0.ithObj == -1 || nextBound1.ithObj == -1)
	("This should not happen!\n");
      
      if (nextBound0.ithObj != -1)
342
	child[0]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr);
343
344

      if (nextBound1.ithObj != -1)
345
	child[1]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr);
346
    }
347
348
349
350
351
  }


  void Tetrahedron::prepareNextBound(BoundaryObject &bound, int ithChild) const
  {
352
    FUNCNAME("Tetrahedron::prepareNextBound()");
353

354
355
    switch (bound.subObj) {
    case FACE:      
356
357
      for (ExcludeList::iterator it = bound.excludedSubstructures.begin(); 
   	   it != bound.excludedSubstructures.end(); ++it) {
358
359
	if (it->first == EDGE && it->second != -1)
	  it->second = edgeOfChild[bound.elType][ithChild][it->second];	
360
      }
361
362
363
364
365
366
367
368
369
370
371
372
      
      bound.ithObj = sideOfChild[bound.elType][ithChild][bound.ithObj];
      bound.elType = (bound.elType + 1) % 3;

      break;
    case EDGE:
      bound.ithObj = edgeOfChild[bound.elType][ithChild][bound.ithObj];
      bound.elType = (bound.elType + 1) % 3;
      break;
    default:
      ERROR_EXIT("Should not happen!\n");
    }
373
374
375
  }


376
  void Tetrahedron::getHigherOrderDofs(const FiniteElemSpace* feSpace,
377
				       BoundaryObject bound,
378
				       DofContainer& dofs,
379
380
				       bool baseDofPtr,
				       vector<GeoIndex>* dofGeoIndex) const
381
  {
382
    FUNCNAME("Tetrahedron::getHigherOrderDofs()");
383

384
385
    TEST_EXIT(dofGeoIndex != NULL)("Not yet implemented!\n");

386
387
388
389
390
    switch (bound.subObj) {
    case VERTEX:
      return;
      break;
    case EDGE:
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
      {
	// === Create boundary information objects for children elements. ===

	BoundaryObject nextBound0 = bound;
	prepareNextBound(nextBound0, 0);

	BoundaryObject nextBound1 = bound;
	prepareNextBound(nextBound1, 1);

	// === Check for boundary on children elements. ===

	if ((nextBound0.ithObj >= 0 || nextBound1.ithObj >= 0) &&  child[0]) {
	  // So, the edge is contained in at least on of the children and the
	  // element is also refined. Then we have go down further in refinement
	  // hierarchie.

	  if (bound.reverseMode) {
	    if (nextBound1.ithObj >= 0)
409
410
	      child[1]->getHigherOrderDofs(feSpace, nextBound1, dofs, 
					   baseDofPtr, dofGeoIndex);
411
	    if (nextBound0.ithObj >= 0)
412
413
	      child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs,
					   baseDofPtr, dofGeoIndex);
414
415
	  } else {
	    if (nextBound0.ithObj >= 0)
416
417
	      child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs,
					   baseDofPtr, dofGeoIndex);
418
	    if (nextBound1.ithObj >= 0)
419
420
	      child[1]->getHigherOrderDofs(feSpace, nextBound1, dofs,
					   baseDofPtr, dofGeoIndex);
421
422
423
424
425
426
427
428
	  }
	} else {
	  // Either the edge is not contained in further refined children, or
	  // the element is not refined further on this edge. Then we can get
	  // all the DOFs on this edge.

	  ElementDofIterator elDofIter(feSpace, true);
	  elDofIter.reset(this);
429
430
431
432
433
434
435
436
437
438
439
440
441
442

	  if (baseDofPtr) {
	    do {
	      if (elDofIter.getCurrentPos() == 1 && 
		  elDofIter.getCurrentElementPos() == bound.ithObj)
		dofs.push_back(elDofIter.getBaseDof());	
	    } while (elDofIter.nextStrict());
	  } else {
	    do {
	      if (elDofIter.getCurrentPos() == 1 && 
		  elDofIter.getCurrentElementPos() == bound.ithObj)
		dofs.push_back(elDofIter.getDofPtr());	
	    } while (elDofIter.next());
	  }
443
444
445
	}
      }

446
447
448
449
450
      break;
    case FACE:
      {      
	BoundaryObject nextBound = bound;
	nextBound.elType = (bound.elType + 1) % 3;
451

452
453
454
455
456
457
458
459
460
	if (child[0]) {
	  int childFace0 = sideOfChild[bound.elType][0][bound.ithObj];
	  int childFace1 = sideOfChild[bound.elType][1][bound.ithObj];
	  
	  TEST_EXIT(childFace0 != -1 || childFace1 != -1)
	    ("No new face for child elements!\n");
	  
	  if (childFace0 != -1) {
	    nextBound.ithObj = childFace0;
461
462
	    child[0]->getHigherOrderDofs(feSpace, nextBound, dofs,
					 baseDofPtr, dofGeoIndex);
463
464
465
466
	  }
	  
	  if (childFace1 != -1) {
	    nextBound.ithObj = childFace1;
467
468
	    child[1]->getHigherOrderDofs(feSpace, nextBound, dofs,
					 baseDofPtr, dofGeoIndex);
469
470
471
472
	  }
	} else { 
	  ElementDofIterator elDofIter(feSpace, true);
	  elDofIter.reset(this);
473
474
475
476
477
478
479
480
481
482
483
484
485
	  if (baseDofPtr) {
	    do {
	      if (elDofIter.getCurrentPos() == 2 && 
		  elDofIter.getCurrentElementPos() == bound.ithObj)
		dofs.push_back(elDofIter.getBaseDof());	
	    } while (elDofIter.nextStrict());
	  } else {
	    do {
	      if (elDofIter.getCurrentPos() == 2 && 
		  elDofIter.getCurrentElementPos() == bound.ithObj)
		dofs.push_back(elDofIter.getDofPtr());	
	    } while (elDofIter.next());      
	  }
486
	}
487
      }
488
      break;
489
    default:
490
      ERROR_EXIT("Should not happen!\n");
491
    }
492
493
  }

494
}