Tetrahedron.cc 11.3 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
  const unsigned char Tetrahedron::adjacentChild[2][2] = {{0,1}, {1,0}};

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

33
34
35
36
  const unsigned char Tetrahedron::edgeOfDofs[4][4] =  {{255, 0, 1, 2},
							{0, 255, 3, 4},
							{1, 3, 255, 5},
							{2, 4, 5, 255}};
37
38
39
40
41
42
43

  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
    FUNCNAME("Tetrahedron::getVertexDofsAtFace()");

    //    MSG("go at el %d with face %d\n", this->getIndex(), bound.ithObj);

221
222
223
224
    if (!child[0])
      return;

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

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

277
278
279

  void Tetrahedron::getVertexDofsAtEdge(FiniteElemSpace* feSpace, 
					BoundaryObject bound,
280
					DofContainer& dofs) const
281
  {
282
283
284
    FUNCNAME("Tetrahedron::getVertexDofsAtEdge()");

    if (!child[0])
285
      return;    
286
287
288
289
290

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

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

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

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

      if (nextBound0.ithObj != -1)
311
312
313
314
315
316
317
318
319
320
	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)
321
322
	child[1]->getVertexDofs(feSpace, nextBound1, dofs);
    }
323
324
325
326
327
  }


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

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


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

358
359
360
361
362
363
364
365
366
367
    switch (bound.subObj) {
    case VERTEX:
      return;
      break;
    case EDGE:
      break;
    case FACE:
      {      
	BoundaryObject nextBound = bound;
	nextBound.elType = (bound.elType + 1) % 3;
368

369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
	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;
	    child[0]->getNonVertexDofs(feSpace, nextBound, dofs);
	  }
	  
	  if (childFace1 != -1) {
	    nextBound.ithObj = childFace1;
	    child[1]->getNonVertexDofs(feSpace, nextBound, dofs);
	  }
	} 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());      
	}
396
      }
397
      break;
398
    default:
399
      ERROR_EXIT("Should not happen!\n");
400
    }
401
402
  }

403
}