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
8
9
10
11
12

namespace AMDiS {

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

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

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

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

26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
  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}};
43

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

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

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

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

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

107
108
  const FixVec<int, WORLD>& Tetrahedron::sortFaceIndices(int face, 
							 FixVec<int,WORLD> *vec) const
109
110
111
  {
    static MatrixOfFixVecs<FixVec<int,WORLD> > *sorted_3d = NULL;
    FixVec<int,WORLD> *val = NULL;
112
    const int *vof = vertexOfFace[face];
113

114
115
    if (NULL == sorted_3d) {
      sorted_3d = new MatrixOfFixVecs<FixVec<int,WORLD> >(3, 4, 7, NO_INIT);
116
 
117
118
119
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
      (*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;
162
163
    }

164
    int no = 0;
165
166
167
168
169
170
171
    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;

172
    if (!(no >= 1  &&  no <= 6))
173
      ERROR_EXIT("Cannot sort face indices of element %d at face %d\n", 
174
175
176
177
178
		 getIndex(), face);

    if (vec) {
      val = vec;
      *val = (*sorted_3d)[face][no];
179
    } else {
180
      val = &((*sorted_3d)[face][no]);
181
    }
182
183
184
185

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

186

187
  void Tetrahedron::getVertexDofs(FiniteElemSpace* feSpace, 
188
				  BoundaryObject bound,
189
190
				  DofContainer& dofs, 
				  bool parentVertices) const
191
192
  {
    FUNCNAME("Tetrahedron::getVertexDofs()");
193

194
195
    switch (bound.subObj) {
    case EDGE:
196
      getVertexDofsAtEdge(feSpace, bound, dofs, parentVertices);
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
      break;
    case FACE:
      getVertexDofsAtFace(feSpace, bound, dofs, parentVertices);
      break;      
    default:
      ERROR_EXIT("Should not happen!\n");
    }
  }


  void Tetrahedron::getVertexDofsAtFace(FiniteElemSpace* feSpace, 
					BoundaryObject bound,
					DofContainer& dofs, 
					bool parentVertices) const
  {

213
    if (parentVertices) {
214
215
216
      ERROR_EXIT("Einbau notIncludedSubStructures!\n");

      switch (bound.ithObj) {
217
      case 0:
218
	dofs.push_back(dof[1]); dofs.push_back(dof[2]); dofs.push_back(dof[3]);
219
220
	break;
      case 1:
221
	dofs.push_back(dof[0]); dofs.push_back(dof[2]); dofs.push_back(dof[3]);
222
223
	break;
      case 2:
224
	dofs.push_back(dof[0]); dofs.push_back(dof[1]); dofs.push_back(dof[3]);
225
226
	break;
      case 3:
227
	dofs.push_back(dof[0]); dofs.push_back(dof[1]); dofs.push_back(dof[2]);
228
229
230
231
232
233
	break;
      default:
	ERROR_EXIT("Should never happen!\n");
      }
    }

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

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

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

290
291
292
293
294
295

  void Tetrahedron::getVertexDofsAtEdge(FiniteElemSpace* feSpace, 
					BoundaryObject bound,
					DofContainer& dofs, 
					bool parentVertices) const
  {
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
    FUNCNAME("Tetrahedron::getVertexDofsAtEdge()");

    if (parentVertices) {
      ERROR_EXIT("Einbau notIncludedSubStructures!\n");
    }

    if (!child[0])
      return;

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

    if (bound.ithObj == 0) {
      child[0]->getVertexDofs(feSpace, nextBound0, dofs);
      dofs.push_back(child[0]->getDOF(3));
      child[1]->getVertexDofs(feSpace, nextBound1, dofs);
    } else {
      if (nextBound0.ithObj != -1)
	child[0]->getVertexDofs(feSpace, nextBound0, dofs);

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


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

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


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

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

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

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

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

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

385
}