Tetrahedron.cc 11.7 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13 14 15 16 17
#include "Tetrahedron.h"
#include "DOFAdmin.h"
#include "Mesh.h"
#include "CoarseningManager.h"
#include "FixVec.h"
18
#include "ElementDofIterator.h"
19
#include "Global.h"
20 21 22 23 24 25

namespace AMDiS {

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

27 28 29
  const unsigned char Tetrahedron::nChildFace[3][2][2] = {{{1,2},{2,1}},
							  {{1,2},{1,2}},
							  {{1,2},{1,2}}};
30

31 32 33
  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}}};
34

35 36 37
  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}}};
38

39 40 41 42 43 44
  const unsigned char Tetrahedron::adjacentChild[2][2] = {{0,1}, {1,0}};

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

45 46 47 48
  const unsigned char Tetrahedron::edgeOfDofs[4][4] =  {{255, 0, 1, 2},
							{0, 255, 3, 4},
							{1, 3, 255, 5},
							{2, 4, 5, 255}};
49 50 51 52 53 54 55

  const int Tetrahedron::vertexOfEdge[6][2] = {{0,1},
					       {0,2},
					       {0,3},
					       {1,2},
					       {1,3},
					       {2,3}};
56

57 58 59 60 61 62 63 64 65 66 67 68
  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}}};

69 70 71 72 73 74 75
  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}}};

76 77 78 79 80 81 82 83 84 85 86 87
  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

88

89 90
  bool Tetrahedron::hasSide(Element* sideElem) const
  {
91 92
    FUNCNAME("Tetrahedron::hasSide()");
    TEST_EXIT_DBG(sideElem->isTriangle())("called for sideElem-type != Triangle\n");
93 94 95 96
    ERROR_EXIT("not yet\n");
    return false;
  }

97

98
  int Tetrahedron::getVertexOfPosition(GeoIndex position,
99 100
				       int positionIndex,
				       int vertexIndex) const 
101
  {
102
    FUNCNAME("Triangle::getVertexOfPosition()");
103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121
    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;
    }
  }

122

123 124
  const FixVec<int, WORLD>& Tetrahedron::sortFaceIndices(int face, 
							 FixVec<int,WORLD> *vec) const
125 126 127
  {
    static MatrixOfFixVecs<FixVec<int,WORLD> > *sorted_3d = NULL;
    FixVec<int,WORLD> *val = NULL;
128
    const int *vof = vertexOfFace[face];
129

130 131
    if (NULL == sorted_3d) {
      sorted_3d = new MatrixOfFixVecs<FixVec<int,WORLD> >(3, 4, 7, NO_INIT);
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 166 167 168 169 170 171 172 173 174 175 176 177
      (*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;
178 179
    }

180
    int no = 0;
181 182 183 184 185 186 187
    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;

188
    if (!(no >= 1  &&  no <= 6))
189
      ERROR_EXIT("Cannot sort face indices of element %d at face %d\n", 
190 191 192 193 194
		 getIndex(), face);

    if (vec) {
      val = vec;
      *val = (*sorted_3d)[face][no];
195
    } else {
196
      val = &((*sorted_3d)[face][no]);
197
    }
198 199 200 201

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

202

203
  void Tetrahedron::getNodeDofs(const FiniteElemSpace* feSpace, 
204 205
				BoundaryObject bound,
				DofContainer& dofs) const
206
  {
207
    FUNCNAME("Tetrahedron::getNodeDofs()");
208

209
    switch (bound.subObj) {
Thomas Witkowski's avatar
Thomas Witkowski committed
210
    case VERTEX:
211 212 213 214
      {
	int n0 = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
	dofs.push_back(&(dof[bound.ithObj][n0]));
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
215
      break;
216
    case EDGE:
217
      getNodeDofsAtEdge(feSpace, bound, dofs);
218 219
      break;
    case FACE:
220
      getNodeDofsAtFace(feSpace, bound, dofs);
221 222 223 224 225 226 227
      break;      
    default:
      ERROR_EXIT("Should not happen!\n");
    }
  }


228
  void Tetrahedron::getNodeDofsAtFace(const FiniteElemSpace* feSpace, 
229 230
				      BoundaryObject bound,
				      DofContainer& dofs) const
231
  {
232
    FUNCNAME("Tetrahedron::getNodeDofsAtFace()");
233

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

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

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

291

292
  void Tetrahedron::getNodeDofsAtEdge(const FiniteElemSpace* feSpace, 
293 294
				      BoundaryObject bound,
				      DofContainer& dofs) const
295
  {
296
    FUNCNAME("Tetrahedron::getNodeDofsAtEdge()");
297 298

    if (!child[0])
299
      return;    
300

301
    int n0 = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
302 303 304 305
    BoundaryObject nextBound0 = bound, nextBound1 = bound;
    prepareNextBound(nextBound0, 0);
    prepareNextBound(nextBound1, 1);

306 307 308 309 310 311
    switch (bound.ithObj) {
    case 0:      
      nextBound0.reverseMode = false;
      nextBound1.reverseMode = false;

      if (bound.reverseMode) {
312
	child[1]->getNodeDofs(feSpace, nextBound0, dofs);
313
	dofs.push_back(&(child[0]->getDof()[3][n0]));
314
	child[0]->getNodeDofs(feSpace, nextBound1, dofs);
315
      } else {
316
	child[0]->getNodeDofs(feSpace, nextBound0, dofs);
317
	dofs.push_back(&(child[0]->getDof()[3][n0]));
318
	child[1]->getNodeDofs(feSpace, nextBound1, dofs);
319 320 321 322
      }

      break;
    case 5:
323 324
      TEST_EXIT_DBG(nextBound0.ithObj == nextBound1.ithObj)
	("Should not happen!\n");
325 326

      if (nextBound0.ithObj != -1)
327
	child[0]->getNodeDofs(feSpace, nextBound0, dofs);      
328 329 330 331 332 333
      break;
    default:
      TEST_EXIT_DBG(nextBound0.ithObj == -1 || nextBound1.ithObj == -1)
	("This should not happen!\n");
      
      if (nextBound0.ithObj != -1)
334
	child[0]->getNodeDofs(feSpace, nextBound0, dofs);      
335 336

      if (nextBound1.ithObj != -1)
337
	child[1]->getNodeDofs(feSpace, nextBound1, dofs);
338
    }
339 340 341 342 343
  }


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

346 347
    switch (bound.subObj) {
    case FACE:      
348 349
      for (ExcludeList::iterator it = bound.excludedSubstructures.begin(); 
   	   it != bound.excludedSubstructures.end(); ++it) {
350 351
	if (it->first == EDGE && it->second != -1)
	  it->second = edgeOfChild[bound.elType][ithChild][it->second];	
352
      }
353 354 355 356 357 358 359 360 361 362 363 364
      
      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");
    }
365 366 367
  }


368
  void Tetrahedron::getHigherOrderDofs(const FiniteElemSpace* feSpace,
369 370
				       BoundaryObject bound,
				       DofContainer& dofs) const
371
  {
372
    FUNCNAME("Tetrahedron::getHigherOrderDofs()");
373

374 375 376 377 378 379 380 381 382 383
    switch (bound.subObj) {
    case VERTEX:
      return;
      break;
    case EDGE:
      break;
    case FACE:
      {      
	BoundaryObject nextBound = bound;
	nextBound.elType = (bound.elType + 1) % 3;
384

385 386 387 388 389 390 391 392 393
	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;
394
	    child[0]->getHigherOrderDofs(feSpace, nextBound, dofs);
395 396 397 398
	  }
	  
	  if (childFace1 != -1) {
	    nextBound.ithObj = childFace1;
399
	    child[1]->getHigherOrderDofs(feSpace, nextBound, dofs);
400 401 402 403 404 405 406 407 408 409 410 411
	  }
	} 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());      
	}
412
      }
413
      break;
414
    default:
415
      ERROR_EXIT("Should not happen!\n");
416
    }
417 418
  }

419
}