Tetrahedron.cc 9.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
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
    switch (bound.subObj) {
    case EDGE:
      //      ERROR_EXIT("Not yet implemented: %d\n", bound.elIndex);
      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
296
297
298
299
300
301
302
303
304
305
306
307
308
309

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


  void Tetrahedron::prepareNextBound(BoundaryObject &bound, int ithChild) const
  {
    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;
  }


310
  void Tetrahedron::getNonVertexDofs(FiniteElemSpace* feSpace,
311
				     BoundaryObject bound,
312
313
314
				     DofContainer& dofs) const
  {
    FUNCNAME("Tetrahedron::getNonVertexDofs()");
315

316
317
318
    BoundaryObject nextBound = bound;
    nextBound.elType = (bound.elType + 1) % 3;

319
    if (child[0]) {
320
321
      int childFace0 = sideOfChild[bound.elType][0][bound.ithObj];
      int childFace1 = sideOfChild[bound.elType][1][bound.ithObj];
322
323
324
325

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

326
327
328
329
      if (childFace0 != -1) {
	nextBound.ithObj = childFace0;
	child[0]->getNonVertexDofs(feSpace, nextBound, dofs);
      }
330

331
332
333
334
      if (childFace1 != -1) {
	nextBound.ithObj = childFace1;
	child[1]->getNonVertexDofs(feSpace, nextBound, dofs);
      }
335
336
337
338
339
    } else { 
      ElementDofIterator elDofIter(feSpace, true);
      elDofIter.reset(this);
      do {
	if (elDofIter.getCurrentPos() == 2 && 
340
	    elDofIter.getCurrentElementPos() == bound.ithObj) {
341
342
343
	  ERROR_EXIT("Check this, if it will really work!\n");
	  dofs.push_back(elDofIter.getDofPtr());	
	}	
344
345
      } while(elDofIter.next());      
    }
346
347
  }

348
}