Tetrahedron.cc 10.9 KB
Newer Older
1
2
3
4
5
#include "Tetrahedron.h"
#include "DOFAdmin.h"
#include "Mesh.h"
#include "CoarseningManager.h"
#include "FixVec.h"
6
#include "ElementDofIterator.h"
7
#include "Global.h"
8
9
10
11
12
13

namespace AMDiS {

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

15
16
17
  const unsigned char Tetrahedron::nChildFace[3][2][2] = {{{1,2},{2,1}},
							  {{1,2},{1,2}},
							  {{1,2},{1,2}}};
18

19
20
21
  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}}};
22

23
24
25
  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}}};
26

27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
  const unsigned char Tetrahedron::adjacentChild[2][2] = {{0,1}, {1,0}};

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

  const unsigned char Tetrahedron::edgeOfDOFs[4][4] =  {{255,0,1,2},
							{0,255,3,4},
							{1,3,255,5},
							{2,4,5,255}};

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

45
46
47
48
49
50
51
52
53
54
55
56
  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}}};

57
58
59
60
61
62
63
  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}}};

64
65
66
67
68
69
70
71
72
73
74
75
  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

76

77
78
  bool Tetrahedron::hasSide(Element* sideElem) const
  {
79
80
    FUNCNAME("Tetrahedron::hasSide()");
    TEST_EXIT_DBG(sideElem->isTriangle())("called for sideElem-type != Triangle\n");
81
82
83
84
    ERROR_EXIT("not yet\n");
    return false;
  }

85

86
  int Tetrahedron::getVertexOfPosition(GeoIndex position,
87
88
				       int positionIndex,
				       int vertexIndex) const 
89
  {
90
    FUNCNAME("Triangle::getVertexOfPosition()");
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
    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;
    }
  }

110

111
112
  const FixVec<int, WORLD>& Tetrahedron::sortFaceIndices(int face, 
							 FixVec<int,WORLD> *vec) const
113
114
115
  {
    static MatrixOfFixVecs<FixVec<int,WORLD> > *sorted_3d = NULL;
    FixVec<int,WORLD> *val = NULL;
116
    const int *vof = vertexOfFace[face];
117

118
119
    if (NULL == sorted_3d) {
      sorted_3d = new MatrixOfFixVecs<FixVec<int,WORLD> >(3, 4, 7, NO_INIT);
120
 
121
122
123
124
125
126
127
128
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
      (*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;
166
167
    }

168
    int no = 0;
169
170
171
172
173
174
175
    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;

176
    if (!(no >= 1  &&  no <= 6))
177
      ERROR_EXIT("Cannot sort face indices of element %d at face %d\n", 
178
179
180
181
182
		 getIndex(), face);

    if (vec) {
      val = vec;
      *val = (*sorted_3d)[face][no];
183
    } else {
184
      val = &((*sorted_3d)[face][no]);
185
    }
186
187
188
189

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

190

191
  void Tetrahedron::getVertexDofs(FiniteElemSpace* feSpace, 
192
				  BoundaryObject bound,
193
				  DofContainer& dofs) const
194
195
  {
    FUNCNAME("Tetrahedron::getVertexDofs()");
196

197
    switch (bound.subObj) {
Thomas Witkowski's avatar
Thomas Witkowski committed
198
199
200
    case VERTEX:
      dofs.push_back(dof[bound.ithObj]);
      break;
201
    case EDGE:
202
      getVertexDofsAtEdge(feSpace, bound, dofs);
203
204
      break;
    case FACE:
205
      getVertexDofsAtFace(feSpace, bound, dofs);
206
207
208
209
210
211
212
213
214
      break;      
    default:
      ERROR_EXIT("Should not happen!\n");
    }
  }


  void Tetrahedron::getVertexDofsAtFace(FiniteElemSpace* feSpace, 
					BoundaryObject bound,
215
					DofContainer& dofs) const
216
217
218
219
220
  {
    if (!child[0])
      return;

    switch (bound.ithObj) {
221
222
    case 0:
      {
223
224
225
	BoundaryObject nextBound = bound;
	prepareNextBound(nextBound, 1);
	child[1]->getVertexDofs(feSpace, nextBound, dofs);	
226
227
228
229
      }
      break;
    case 1:
      {
230
231
232
	BoundaryObject nextBound = bound;
	prepareNextBound(nextBound, 0);
	child[0]->getVertexDofs(feSpace, nextBound, dofs);	
233
234
235
236
237
238
      }
      break;

    case 2:
    case 3:
      {
239
240
241
242
243
244
245
	BoundaryObject nextBound0 = bound, nextBound1 = bound;
	prepareNextBound(nextBound0, 0);
	prepareNextBound(nextBound1, 1);
	nextBound0.reverseMode = false;
	nextBound1.reverseMode = false;
	
	bool addDof = true;
246
247
	for (ExcludeList::iterator it = bound.excludedSubstructures.begin(); 
	     it != bound.excludedSubstructures.end(); ++it)
248
	  if (it->first == EDGE && it->second == 0)
249
250
251
252
253
254
	    addDof = false;
	
	if (bound.reverseMode) {
	  child[1]->getVertexDofs(feSpace, nextBound1, dofs);
	  
	  if (addDof)
255
	    dofs.push_back(child[0]->getDOF(3));
256
257
258
259
260
261
	  
	  child[0]->getVertexDofs(feSpace, nextBound0, dofs);
	} else {
	  child[0]->getVertexDofs(feSpace, nextBound0, dofs);
	  
	  if (addDof)
262
	    dofs.push_back(child[0]->getDOF(3));
263
264
	  
	  child[1]->getVertexDofs(feSpace, nextBound1, dofs);	  
265
266
267
268
	}
      }
      break;
    default:
269
      ERROR_EXIT("Should never happen: %d\n", bound.ithObj);
270
    }
271
272
  }

273
274
275

  void Tetrahedron::getVertexDofsAtEdge(FiniteElemSpace* feSpace, 
					BoundaryObject bound,
276
					DofContainer& dofs) const
277
  {
278
279
280
    FUNCNAME("Tetrahedron::getVertexDofsAtEdge()");

    if (!child[0])
281
      return;    
282
283
284
285
286

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

287
288
289
290
291
292
293
294
295
296
    switch (bound.ithObj) {
    case 0:      
      nextBound0.reverseMode = false;
      nextBound1.reverseMode = false;

      if (bound.reverseMode) {
	child[1]->getVertexDofs(feSpace, nextBound0, dofs);
	dofs.push_back(child[0]->getDOF(3));
	child[0]->getVertexDofs(feSpace, nextBound1, dofs);
      } else {
297
	child[0]->getVertexDofs(feSpace, nextBound0, dofs);
298
299
300
301
302
303
304
	dofs.push_back(child[0]->getDOF(3));
	child[1]->getVertexDofs(feSpace, nextBound1, dofs);
      }

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

      if (nextBound0.ithObj != -1)
307
308
309
310
311
312
313
314
315
316
	child[0]->getVertexDofs(feSpace, nextBound0, dofs);      
      break;
    default:
      TEST_EXIT_DBG(nextBound0.ithObj == -1 || nextBound1.ithObj == -1)
	("This should not happen!\n");
      
      if (nextBound0.ithObj != -1)
	child[0]->getVertexDofs(feSpace, nextBound0, dofs);      

      if (nextBound1.ithObj != -1)
317
318
	child[1]->getVertexDofs(feSpace, nextBound1, dofs);
    }
319
320
321
322
323
  }


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

326
327
    switch (bound.subObj) {
    case FACE:      
328
329
      for (ExcludeList::iterator it = bound.excludedSubstructures.begin(); 
   	   it != bound.excludedSubstructures.end(); ++it) {
330
331
332
	if (it->first == EDGE)
	  it->second = edgeOfChild[bound.elType][ithChild][it->second];
      }
333
334
335
336
337
338
339
340
341
342
343
344
      
      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");
    }
345
346
347
  }


348
  void Tetrahedron::getNonVertexDofs(FiniteElemSpace* feSpace,
349
				     BoundaryObject bound,
350
351
352
				     DofContainer& dofs) const
  {
    FUNCNAME("Tetrahedron::getNonVertexDofs()");
353

354
355
356
    BoundaryObject nextBound = bound;
    nextBound.elType = (bound.elType + 1) % 3;

357
    if (child[0]) {
358
359
      int childFace0 = sideOfChild[bound.elType][0][bound.ithObj];
      int childFace1 = sideOfChild[bound.elType][1][bound.ithObj];
360
361
362
363

      TEST_EXIT(childFace0 != -1 || childFace1 != -1)
	("No new face for child elements!\n");

364
365
366
367
      if (childFace0 != -1) {
	nextBound.ithObj = childFace0;
	child[0]->getNonVertexDofs(feSpace, nextBound, dofs);
      }
368

369
370
371
372
      if (childFace1 != -1) {
	nextBound.ithObj = childFace1;
	child[1]->getNonVertexDofs(feSpace, nextBound, dofs);
      }
373
374
375
376
377
    } else { 
      ElementDofIterator elDofIter(feSpace, true);
      elDofIter.reset(this);
      do {
	if (elDofIter.getCurrentPos() == 2 && 
378
	    elDofIter.getCurrentElementPos() == bound.ithObj) {
379
380
381
	  ERROR_EXIT("Check this, if it will really work!\n");
	  dofs.push_back(elDofIter.getDofPtr());	
	}	
382
383
      } while(elDofIter.next());      
    }
384
385
  }

386
}