Tetrahedron.cc 10.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
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
				  DofContainer& dofs) const
191
192
  {
    FUNCNAME("Tetrahedron::getVertexDofs()");
193

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


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

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

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

270
271
272

  void Tetrahedron::getVertexDofsAtEdge(FiniteElemSpace* feSpace, 
					BoundaryObject bound,
273
					DofContainer& dofs) const
274
  {
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
    FUNCNAME("Tetrahedron::getVertexDofsAtEdge()");

    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);
    }
295
296
297
298
299
  }


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

302
303
    switch (bound.subObj) {
    case FACE:      
304
305
      for (ExcludeList::iterator it = bound.excludedSubstructures.begin(); 
   	   it != bound.excludedSubstructures.end(); ++it) {
306
307
308
	if (it->first == EDGE)
	  it->second = edgeOfChild[bound.elType][ithChild][it->second];
      }
309
310
311
312
313
314
315
316
317
318
319
320
      
      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");
    }
321
322
323
  }


324
  void Tetrahedron::getNonVertexDofs(FiniteElemSpace* feSpace,
325
				     BoundaryObject bound,
326
327
328
				     DofContainer& dofs) const
  {
    FUNCNAME("Tetrahedron::getNonVertexDofs()");
329

330
331
332
    BoundaryObject nextBound = bound;
    nextBound.elType = (bound.elType + 1) % 3;

333
    if (child[0]) {
334
335
      int childFace0 = sideOfChild[bound.elType][0][bound.ithObj];
      int childFace1 = sideOfChild[bound.elType][1][bound.ithObj];
336
337
338
339

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

340
341
342
343
      if (childFace0 != -1) {
	nextBound.ithObj = childFace0;
	child[0]->getNonVertexDofs(feSpace, nextBound, dofs);
      }
344

345
346
347
348
      if (childFace1 != -1) {
	nextBound.ithObj = childFace1;
	child[1]->getNonVertexDofs(feSpace, nextBound, dofs);
      }
349
350
351
352
353
    } else { 
      ElementDofIterator elDofIter(feSpace, true);
      elDofIter.reset(this);
      do {
	if (elDofIter.getCurrentPos() == 2 && 
354
	    elDofIter.getCurrentElementPos() == bound.ithObj) {
355
356
357
	  ERROR_EXIT("Check this, if it will really work!\n");
	  dofs.push_back(elDofIter.getDofPtr());	
	}	
358
359
      } while(elDofIter.next());      
    }
360
361
  }

362
}