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

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


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

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

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

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

287

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

    if (!child[0])
295
      return;    
296
297
298
299
300

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

301
302
303
304
305
306
    switch (bound.ithObj) {
    case 0:      
      nextBound0.reverseMode = false;
      nextBound1.reverseMode = false;

      if (bound.reverseMode) {
307
	child[1]->getNodeDofs(feSpace, nextBound0, dofs);
308
	dofs.push_back(child[0]->getDof(3));
309
	child[0]->getNodeDofs(feSpace, nextBound1, dofs);
310
      } else {
311
	child[0]->getNodeDofs(feSpace, nextBound0, dofs);
312
	dofs.push_back(child[0]->getDof(3));
313
	child[1]->getNodeDofs(feSpace, nextBound1, dofs);
314
315
316
317
318
      }

      break;
    case 5:
      TEST_EXIT_DBG(nextBound0.ithObj == nextBound1.ithObj)("Should not happen!\n");
319
320

      if (nextBound0.ithObj != -1)
321
	child[0]->getNodeDofs(feSpace, nextBound0, dofs);      
322
323
324
325
326
327
      break;
    default:
      TEST_EXIT_DBG(nextBound0.ithObj == -1 || nextBound1.ithObj == -1)
	("This should not happen!\n");
      
      if (nextBound0.ithObj != -1)
328
	child[0]->getNodeDofs(feSpace, nextBound0, dofs);      
329
330

      if (nextBound1.ithObj != -1)
331
	child[1]->getNodeDofs(feSpace, nextBound1, dofs);
332
    }
333
334
335
336
337
  }


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

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


362
363
364
  void Tetrahedron::getHigherOrderDofs(FiniteElemSpace* feSpace,
				       BoundaryObject bound,
				       DofContainer& dofs) const
365
  {
366
    FUNCNAME("Tetrahedron::getHigherOrderDofs()");
367

368
369
370
371
372
373
374
375
376
377
    switch (bound.subObj) {
    case VERTEX:
      return;
      break;
    case EDGE:
      break;
    case FACE:
      {      
	BoundaryObject nextBound = bound;
	nextBound.elType = (bound.elType + 1) % 3;
378

379
380
381
382
383
384
385
386
387
	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;
388
	    child[0]->getHigherOrderDofs(feSpace, nextBound, dofs);
389
390
391
392
	  }
	  
	  if (childFace1 != -1) {
	    nextBound.ithObj = childFace1;
393
	    child[1]->getHigherOrderDofs(feSpace, nextBound, dofs);
394
395
396
397
398
399
400
401
402
403
404
405
	  }
	} 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());      
	}
406
      }
407
      break;
408
    default:
409
      ERROR_EXIT("Should not happen!\n");
410
    }
411
412
  }

413
}