Tetrahedron.cc 11.7 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
205
				BoundaryObject bound,
				DofContainer& dofs) const
206
  {
207
    FUNCNAME("Tetrahedron::getNodeDofs()");
208

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


228
  void Tetrahedron::getNodeDofsAtFace(const FiniteElemSpace* feSpace, 
229
230
				      BoundaryObject bound,
				      DofContainer& dofs) const
231
  {
232
    FUNCNAME("Tetrahedron::getNodeDofsAtFace()");
233

234
235
236
237
    if (!child[0])
      return;

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

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

291

292
  void Tetrahedron::getNodeDofsAtEdge(const FiniteElemSpace* feSpace, 
293
294
				      BoundaryObject bound,
				      DofContainer& dofs) const
295
  {
296
    FUNCNAME("Tetrahedron::getNodeDofsAtEdge()");
297
298

    if (!child[0])
299
      return;    
300

301
    int n0 = 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
	child[1]->getNodeDofs(feSpace, nextBound0, dofs);
313
	dofs.push_back(&(child[0]->getDof()[3][n0]));
314
	child[0]->getNodeDofs(feSpace, nextBound1, dofs);
315
      } else {
316
	child[0]->getNodeDofs(feSpace, nextBound0, dofs);
317
	dofs.push_back(&(child[0]->getDof()[3][n0]));
318
	child[1]->getNodeDofs(feSpace, nextBound1, dofs);
319
320
321
322
      }

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

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

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


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

346
347
    switch (bound.subObj) {
    case FACE:      
348
349
      for (ExcludeList::iterator it = bound.excludedSubstructures.begin(); 
   	   it != bound.excludedSubstructures.end(); ++it) {
350
351
	if (it->first == EDGE && it->second != -1)
	  it->second = edgeOfChild[bound.elType][ithChild][it->second];	
352
      }
353
354
355
356
357
358
359
360
361
362
363
364
      
      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");
    }
365
366
367
  }


368
  void Tetrahedron::getHigherOrderDofs(const FiniteElemSpace* feSpace,
369
370
				       BoundaryObject bound,
				       DofContainer& dofs) const
371
  {
372
    FUNCNAME("Tetrahedron::getHigherOrderDofs()");
373

374
375
376
377
378
379
380
381
382
383
    switch (bound.subObj) {
    case VERTEX:
      return;
      break;
    case EDGE:
      break;
    case FACE:
      {      
	BoundaryObject nextBound = bound;
	nextBound.elType = (bound.elType + 1) % 3;
384

385
386
387
388
389
390
391
392
393
	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;
394
	    child[0]->getHigherOrderDofs(feSpace, nextBound, dofs);
395
396
397
398
	  }
	  
	  if (childFace1 != -1) {
	    nextBound.ithObj = childFace1;
399
	    child[1]->getHigherOrderDofs(feSpace, nextBound, dofs);
400
401
402
403
404
405
406
407
408
409
410
411
	  }
	} else { 
	  ElementDofIterator elDofIter(feSpace, true);
	  elDofIter.reset(this);
	  do {
	    if (elDofIter.getCurrentPos() == 2 && 
		elDofIter.getCurrentElementPos() == bound.ithObj) {
	      ERROR_EXIT("Check this, if it will really work!\n");
	      dofs.push_back(elDofIter.getDofPtr());	
	    }	
	  } while(elDofIter.next());      
	}
412
      }
413
      break;
414
    default:
415
      ERROR_EXIT("Should not happen!\n");
416
    }
417
418
  }

419
}