Tetrahedron.cc 15.6 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
  void Tetrahedron::sortFaceIndices(int face, FixVec<int,WORLD> &vec) const
124
125
126
  {
    static MatrixOfFixVecs<FixVec<int,WORLD> > *sorted_3d = NULL;

127
    if (sorted_3d == NULL) {
128
      sorted_3d = new MatrixOfFixVecs<FixVec<int,WORLD> >(3, 4, 7, NO_INIT);
129
 
130
131
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
      (*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;
175
176
    }

177
    const int *vof = vertexOfFace[face];
178
    int no = 0;
179
180
181
182
183
184
185
    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;

186
187
    TEST_EXIT_DBG(no >= 1 && no <= 6)
      ("Cannot sort face indices of element %d at face %d\n", getIndex(), face);
188

189
    vec = (*sorted_3d)[face][no];
190
191
  }

192

193
  void Tetrahedron::getNodeDofs(const FiniteElemSpace* feSpace, 
194
				BoundaryObject bound,
195
196
				DofContainer& dofs,
				bool baseDofPtr) const
197
  {
198
    FUNCNAME("Tetrahedron::getNodeDofs()");
199

200
    switch (bound.subObj) {
Thomas Witkowski's avatar
Thomas Witkowski committed
201
    case VERTEX:
202
      {
203
204
205
206
207
208
	if (baseDofPtr) {
	  dofs.push_back(dof[bound.ithObj]);
	} else {
	  int n0 = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
	  dofs.push_back(&(dof[bound.ithObj][n0]));
	}
209
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
210
      break;
211
    case EDGE:
212
      getNodeDofsAtEdge(feSpace, bound, dofs, baseDofPtr);
213
214
      break;
    case FACE:
215
      getNodeDofsAtFace(feSpace, bound, dofs, baseDofPtr);
216
217
218
219
220
221
222
      break;      
    default:
      ERROR_EXIT("Should not happen!\n");
    }
  }


223
  void Tetrahedron::getNodeDofsAtFace(const FiniteElemSpace* feSpace, 
224
				      BoundaryObject bound,
225
226
				      DofContainer& dofs,
				      bool baseDofPtr) const
227
  {
228
    FUNCNAME("Tetrahedron::getNodeDofsAtFace()");
229

230
    if (!child[0]){
231
      return;
232
	}
233
234

    switch (bound.ithObj) {
235
236
    case 0:
      {
237
238
	BoundaryObject nextBound = bound;
	prepareNextBound(nextBound, 1);
239
	child[1]->getNodeDofs(feSpace, nextBound, dofs, baseDofPtr);
240
241
242
243
      }
      break;
    case 1:
      {
244
245
	BoundaryObject nextBound = bound;
	prepareNextBound(nextBound, 0);
246
	child[0]->getNodeDofs(feSpace, nextBound, dofs, baseDofPtr);
247
248
249
250
251
252
      }
      break;

    case 2:
    case 3:
      {
Thomas Witkowski's avatar
Thomas Witkowski committed
253
	int n0 = (baseDofPtr ? 0 : feSpace->getAdmin()->getNumberOfPreDofs(VERTEX));
254
255
256
257
258
259
260
	BoundaryObject nextBound0 = bound, nextBound1 = bound;
	prepareNextBound(nextBound0, 0);
	prepareNextBound(nextBound1, 1);
	nextBound0.reverseMode = false;
	nextBound1.reverseMode = false;
	
	bool addDof = true;
261
262
	for (ExcludeList::iterator it = bound.excludedSubstructures.begin(); 
	     it != bound.excludedSubstructures.end(); ++it)
263
	  if (it->first == EDGE && it->second == 0)
264
	    addDof = false;	
265
266
	
	if (bound.reverseMode) {
267
	  child[1]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr);
268
269
	  
	  if (addDof)
270
	    dofs.push_back(&(child[0]->getDof()[3][n0]));
271
	  
272
	  child[0]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr);
273
	} else {
274
	  child[0]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr);
275
276
	  
	  if (addDof)
277
	    dofs.push_back(&(child[0]->getDof()[3][n0]));
278
	  
279
	  child[1]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr);
280
281
282
283
	}
      }
      break;
    default:
284
      ERROR_EXIT("Should never happen: %d\n", bound.ithObj);
285
    }
286
287
  }

288

289
  void Tetrahedron::getNodeDofsAtEdge(const FiniteElemSpace* feSpace, 
290
				      BoundaryObject bound,
291
292
				      DofContainer& dofs,
				      bool baseDofPtr) const
293
  {
294
    FUNCNAME("Tetrahedron::getNodeDofsAtEdge()");
295

296
297
298
299

    if (!child[0]){
      return;
	}
300

Thomas Witkowski's avatar
Thomas Witkowski committed
301
    int n0 = (baseDofPtr ? 0 : feSpace->getAdmin()->getNumberOfPreDofs(VERTEX));
302
303
304
305
    BoundaryObject nextBound0 = bound, nextBound1 = bound;
    prepareNextBound(nextBound0, 0);
    prepareNextBound(nextBound1, 1);

306
307
308
309
310
311
    switch (bound.ithObj) {
    case 0:      
      nextBound0.reverseMode = false;
      nextBound1.reverseMode = false;

      if (bound.reverseMode) {
312
313
314
315
			child[1]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr);
			dofs.push_back(&(child[0]->getDof()[3][n0]));
			child[0]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr);

316
      } else {
317
318
319
320
			child[0]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr);
			dofs.push_back(&(child[0]->getDof()[3][n0]));
			child[1]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr);

321
322
323
324
      }

      break;
    case 5:
325
326
      TEST_EXIT_DBG(nextBound0.ithObj == nextBound1.ithObj)
	("Should not happen!\n");
327
328

      if (nextBound0.ithObj != -1)
329
	child[0]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr);
330
331
332
333
334
335
      break;
    default:
      TEST_EXIT_DBG(nextBound0.ithObj == -1 || nextBound1.ithObj == -1)
	("This should not happen!\n");
      
      if (nextBound0.ithObj != -1)
336
	child[0]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr);
337
338

      if (nextBound1.ithObj != -1)
339
	child[1]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr);
340
    }
341

342
343
344
  }


345
  void Tetrahedron::getHigherOrderDofs(const FiniteElemSpace* feSpace,
346
				       BoundaryObject bound,
347
				       DofContainer& dofs,
348
349
				       bool baseDofPtr,
				       vector<GeoIndex>* dofGeoIndex) const
350
  {
351
    FUNCNAME("Tetrahedron::getHigherOrderDofs()");
352

353
354
355
356
357
    switch (bound.subObj) {
    case VERTEX:
      return;
      break;
    case EDGE:
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
      {
	// === 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.

374
375
376
377
378
379
	 	// Additonal Check if Boundary Edge is NOT Refinement Edge of BoundaryObject
	 	// AND the complete Edge is part of both childs
		// therefore only recursion on ONE child is necessary
	  if ( (bound.ithObj > 0) && (nextBound0.ithObj >= 0 && nextBound1.ithObj >= 0) ){
			child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs,
						   baseDofPtr, dofGeoIndex);
380
	  } else {
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
	  	if (bound.reverseMode) {
	    	if (nextBound1.ithObj >= 0)
	    	  child[1]->getHigherOrderDofs(feSpace, nextBound1, dofs, 
						   baseDofPtr, dofGeoIndex);
	    	if (nextBound0.ithObj >= 0)
	    	  child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs,
						   baseDofPtr, dofGeoIndex);
		  } else {
		    if (nextBound0.ithObj >= 0)
		      child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs,
						   baseDofPtr, dofGeoIndex);
		    if (nextBound1.ithObj >= 0)
		      child[1]->getHigherOrderDofs(feSpace, nextBound1, dofs,
						   baseDofPtr, dofGeoIndex);
		  }
396
	  }
397

398
399
400
401
402
403
404
	} 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);
405
406
407
408

	  if (baseDofPtr) {
	    do {
	      if (elDofIter.getCurrentPos() == 1 && 
Thomas Witkowski's avatar
Thomas Witkowski committed
409
		  elDofIter.getCurrentElementPos() == bound.ithObj) {
410
		dofs.push_back(elDofIter.getBaseDof());	
Thomas Witkowski's avatar
Thomas Witkowski committed
411
412
413
		if (dofGeoIndex != NULL)
		  dofGeoIndex->push_back(EDGE);
	      }
414
415
416
417
	    } while (elDofIter.nextStrict());
	  } else {
	    do {
	      if (elDofIter.getCurrentPos() == 1 && 
Thomas Witkowski's avatar
Thomas Witkowski committed
418
		  elDofIter.getCurrentElementPos() == bound.ithObj) {
419
420
421
422
				dofs.push_back(elDofIter.getDofPtr());	

				if (dofGeoIndex != NULL)
				  dofGeoIndex->push_back(EDGE);
Thomas Witkowski's avatar
Thomas Witkowski committed
423
	      }
424
425
	    } while (elDofIter.next());
	  }
426
427
428
	}
      }

429
430
431
432
      break;
    case FACE:
      {      
	if (child[0]) {
433
434
435
436
437
	  BoundaryObject nextBound0 = bound, nextBound1 = bound;
	  prepareNextBound(nextBound0, 0);
	  prepareNextBound(nextBound1, 1);

	  TEST_EXIT_DBG(nextBound0.ithObj != -1 || nextBound1.ithObj != 1)
438
	    ("No new face for child elements!\n");
439
440

	  TEST_EXIT_DBG(!bound.reverseMode)("Not yet implemented!\n");
441
	  
442
443
	  if (nextBound0.ithObj != -1)
	    child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs,
444
					 baseDofPtr, dofGeoIndex);
445
	  
446
447
	  if (nextBound1.ithObj != -1)
	    child[1]->getHigherOrderDofs(feSpace, nextBound1, dofs,
448
					 baseDofPtr, dofGeoIndex);
449
	} else { 
Thomas Witkowski's avatar
Thomas Witkowski committed
450

451
452
453
	  // === On leaf elements iterate over all DOFs and collect the  === 
	  // === right ones.                                             ===

454
455
	  ElementDofIterator elDofIter(feSpace, true);
	  elDofIter.reset(this);
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493

	  bool next = true;
	  do {
	    bool addDof = false;

	    switch (elDofIter.getCurrentPos()) {
	    case 0:
	      // VERTEX is never a higher order DOF
	      addDof = false;
	      break;
	    case 1:
	      // On EDGE, first check whether this edge is part of the FACE and
	      // then check if this EDGE is not inside of excludedSubstructures.

	      {
		int localEdge = elDofIter.getCurrentElementPos();
		if (edgeOfFace[bound.ithObj][0] == localEdge ||
		    edgeOfFace[bound.ithObj][1] == localEdge ||
		    edgeOfFace[bound.ithObj][2] == localEdge) {
		  addDof = true;

		  for (ExcludeList::iterator it = bound.excludedSubstructures.begin(); 
		       it != bound.excludedSubstructures.end(); ++it)
		    if (it->first == EDGE && it->second == localEdge)
		      addDof = false;	
		} else {
		  addDof = false;
		}
	      }

	      break;
	    case 2:
	      // On FACE, check for the right one.
	      if (elDofIter.getCurrentElementPos() == bound.ithObj)
		addDof = true;
	      else
		addDof = false;
	      break;
494
495
496
497
		case 3:
		  // On CENTER, dof can not be part on internal boundary
		addDof = false;
		  break;
498
499
500
501
502
503
504
	    default:
	      ERROR_EXIT("Should not happen!\n");
	    }


	    if (addDof) {
	      if (baseDofPtr)
505
			dofs.push_back(elDofIter.getBaseDof());	
506
	      else
507
			dofs.push_back(elDofIter.getDofPtr());	
508
509
	      
	      if (dofGeoIndex != NULL)
510
511
			dofGeoIndex->push_back(elDofIter.getPosIndex());

512
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
513

514
515
	    next = (baseDofPtr ? elDofIter.nextStrict() : elDofIter.next());	    
	  } while (next);
516
	}
517
      }
518
      break;
519
    default:
520
      ERROR_EXIT("Should not happen!\n");
521
    }
522

523
524
  }

525
526
527

  void Tetrahedron::getSubBoundary(BoundaryObject bound, 
				   vector<BoundaryObject> &subBound) const
528
  {}
529
530
531
532
533
534
535
536
537


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

    switch (bound.subObj) {
    case FACE:      
      for (ExcludeList::iterator it = bound.excludedSubstructures.begin(); 
538
   	   it != bound.excludedSubstructures.end(); ++it)
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
	if (it->first == EDGE && it->second != -1)
	  it->second = edgeOfChild[bound.elType][ithChild][it->second];	
      
      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");
    }
  }

555
}