Tetrahedron.cc 11 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
76
77
  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
  {
78
79
    FUNCNAME("Tetrahedron::hasSide()");
    TEST_EXIT_DBG(sideElem->isTriangle())("called for sideElem-type != Triangle\n");
80
81
82
83
84
    ERROR_EXIT("not yet\n");
    return false;
  }

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

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

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

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

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

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

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

187

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

195
    switch (bound.subObj) {
Thomas Witkowski's avatar
Thomas Witkowski committed
196
197
198
    case VERTEX:
      dofs.push_back(dof[bound.ithObj]);
      break;
199
    case EDGE:
200
      getVertexDofsAtEdge(feSpace, bound, dofs, parentVertices);
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
      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
  {

217
    if (parentVertices) {
218
      ERROR_EXIT("Einbau excludedSubstructures!\n");
219
220

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

238
239
240
241
    if (!child[0])
      return;

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

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

294
295
296
297
298
299

  void Tetrahedron::getVertexDofsAtEdge(FiniteElemSpace* feSpace, 
					BoundaryObject bound,
					DofContainer& dofs, 
					bool parentVertices) const
  {
300
301
302
    FUNCNAME("Tetrahedron::getVertexDofsAtEdge()");

    if (parentVertices) {
303
      ERROR_EXIT("Einbau excludedSubstructures!\n");
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
    }

    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);
    }
324
325
326
327
328
  }


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

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


353
  void Tetrahedron::getNonVertexDofs(FiniteElemSpace* feSpace,
354
				     BoundaryObject bound,
355
356
357
				     DofContainer& dofs) const
  {
    FUNCNAME("Tetrahedron::getNonVertexDofs()");
358

359
360
361
    BoundaryObject nextBound = bound;
    nextBound.elType = (bound.elType + 1) % 3;

362
    if (child[0]) {
363
364
      int childFace0 = sideOfChild[bound.elType][0][bound.ithObj];
      int childFace1 = sideOfChild[bound.elType][1][bound.ithObj];
365
366
367
368

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

369
370
371
372
      if (childFace0 != -1) {
	nextBound.ithObj = childFace0;
	child[0]->getNonVertexDofs(feSpace, nextBound, dofs);
      }
373

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

391
}