Triangle.cc 6.49 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 "Triangle.h"
#include "DOFAdmin.h"
#include "Mesh.h"
#include "CoarseningManager.h"
#include "FixVec.h"
18
#include "ElementDofIterator.h"
19 20 21

namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
22 23 24
  const int Triangle::vertexOfEdge[3][2] = {{1, 2}, {2, 0}, {0, 1}};
  const int Triangle::sideOfChild[2][3] = {{-1, 2, 0}, {2, -1, 1}};
  const int Triangle::vertexOfParent[2][3] = {{2, 0, -1}, {1, 2, -1}};
25 26 27 28

  bool Triangle::hasSide(Element* sideElem) const
  {
    FUNCNAME("Triangle::hasSide");
29
    TEST_EXIT_DBG(sideElem->isLine())("called for sideElem-type != Line\n");
30 31 32 33
    ERROR_EXIT("not yet\n");
    return false;
  }

34 35
  int Triangle::getVertexOfPosition(GeoIndex position, int positionIndex,
				    int vertexIndex) const 
36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
  {
    FUNCNAME("Triangle::getVertexOfPosition");
    switch(position) {
    case VERTEX:
      return positionIndex;
      break;
    case EDGE:
      return vertexOfEdge[positionIndex][vertexIndex];
      break;
    case CENTER:
      return vertexIndex;
      break;
    default:
      ERROR_EXIT("invalid position\n");
      return 0;
    }
  }

54
  const FixVec<int, WORLD>& Triangle::sortFaceIndices(int face, 
55
						      FixVec<int, WORLD> *vec) const
56
  {
57
    static MatrixOfFixVecs<FixVec<int, WORLD> > *sorted_2d = NULL;
58 59

    int no = 0;
60
    FixVec<int, WORLD> *val = NULL;
61
    const int *vof = vertexOfEdge[face];
62

63
    if (NULL == sorted_2d) {
64
      sorted_2d = new MatrixOfFixVecs<FixVec<int, WORLD> >(2, 3, 2, NO_INIT);
65

66 67 68 69 70 71
      (*sorted_2d)[1][0][1] = (*sorted_2d)[1][1][0] =
	(*sorted_2d)[2][0][0] = (*sorted_2d)[2][1][1] = 0;
      (*sorted_2d)[0][0][0] = (*sorted_2d)[0][1][1] =
	(*sorted_2d)[2][0][1] = (*sorted_2d)[2][1][0] = 1;
      (*sorted_2d)[0][0][1] = (*sorted_2d)[0][1][0] =
	(*sorted_2d)[1][0][0] = (*sorted_2d)[1][1][1] = 2;
72 73 74 75 76 77 78 79 80 81
    }

    if (dof[vof[0]][0] < dof[vof[1]][0])
      no = 0;
    else
      no = 1;

    if (vec) {
      val = vec;
      *val = (*sorted_2d)[face][no];
82
    } else {
83
      val = &((*sorted_2d)[face][no]);
84
    }
85

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

89

90
  void Triangle::getNodeDofs(const FiniteElemSpace* feSpace, 
91 92
			     BoundaryObject bound,
			     DofContainer& dofs) const
93
  {
94
    FUNCNAME("Triangle::getNodeDofs()");
95

96 97 98
    // Get displacement if more than one FE space is defined on mesh.
    int n0 = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);

99
    if (bound.subObj == VERTEX) {
100
      dofs.push_back(&(dof[bound.ithObj][n0]));
101 102 103 104 105
      return;
    }

    TEST_EXIT_DBG(bound.subObj == EDGE)("This should not happen!\n");

106 107 108
    BoundaryObject nextBound = bound;

    switch (bound.ithObj) {
109 110
    case 0: 
      {
111
	if (child[1] && child[1]->getFirstChild()) {
112 113
	  const DegreeOfFreedom** elDofs = child[1]->getFirstChild()->getDof();

114 115
 	  if (bound.reverseMode) {
 	    nextBound.ithObj = 1;
116
 	    child[1]->getSecondChild()->getNodeDofs(feSpace, nextBound, dofs);
117
 	    dofs.push_back(&(elDofs[2][n0]));
118
 	    nextBound.ithObj = 0;
119
 	    child[1]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs);
120 121
 	  } else {
	    nextBound.ithObj = 0;
122
	    child[1]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs);
123
	    dofs.push_back(&(elDofs[2][n0]));
124
	    nextBound.ithObj = 1;
125
	    child[1]->getSecondChild()->getNodeDofs(feSpace, nextBound, dofs);
126
	  }
127 128 129 130 131
	}
      }
      break;
    case 1:
      {
132
	if (child[0] && child[0]->getFirstChild()) {
133 134
	  const DegreeOfFreedom** elDofs = child[0]->getFirstChild()->getDof();

135 136
	  if (bound.reverseMode) {
	    nextBound.ithObj = 1;
137
	    child[0]->getSecondChild()->getNodeDofs(feSpace, nextBound, dofs);
138
	    dofs.push_back(&(elDofs[2][n0]));
139
	    nextBound.ithObj = 0;
140
	    child[0]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs);
141 142
	  } else {
	    nextBound.ithObj = 0;
143
	    child[0]->getFirstChild()->getNodeDofs(feSpace, nextBound, dofs);
144
	    dofs.push_back(&(elDofs[2][n0]));
145
	    nextBound.ithObj = 1;
146
	    child[0]->getSecondChild()->getNodeDofs(feSpace, nextBound, dofs);
147
	  }
148 149 150 151 152
	}
      }
      break;
    case 2:      
      if (child[0]) {
153 154
	  const DegreeOfFreedom** elDofs = child[0]->getDof();

155 156
	if (bound.reverseMode) {
	  nextBound.ithObj = 1;
157
	  child[1]->getNodeDofs(feSpace, nextBound, dofs);
158
	  dofs.push_back(&(elDofs[2][n0]));
159
	  nextBound.ithObj = 0;
160
	  child[0]->getNodeDofs(feSpace, nextBound, dofs);
161 162
	} else {
	  nextBound.ithObj = 0;
163
	  child[0]->getNodeDofs(feSpace, nextBound, dofs);	  
164
	  dofs.push_back(&(elDofs[2][n0]));
165
	  nextBound.ithObj = 1;
166
	  child[1]->getNodeDofs(feSpace, nextBound, dofs);
167
	}
168 169 170 171 172 173 174 175
      }
      break;      
    default:
      ERROR_EXIT("Should never happen!\n");
    }
  }


176
  void Triangle::getHigherOrderDofs(const FiniteElemSpace* feSpace,
177 178
				    BoundaryObject bound,
				    DofContainer& dofs) const
179
  {
180
    FUNCNAME("Triange::getHigherOrderDofs()");
181

182 183 184 185 186
    if (bound.subObj == VERTEX)
      return;

    TEST_EXIT_DBG(bound.subObj == EDGE)("This should not happen!\n");

187
    bool addThisEdge = false;
188
    BoundaryObject nextBound = bound;
189

190
    switch (bound.ithObj) {
191
    case 0:
192 193
      if (child[1]) {
	nextBound.ithObj = 2;
194
	child[1]->getHigherOrderDofs(feSpace, nextBound, dofs);
195
      } else {
196
	addThisEdge = true;
197
      }
198 199 200

      break;
    case 1:
201 202
      if (child[0]) {
	nextBound.ithObj = 2;
203
	child[0]->getHigherOrderDofs(feSpace, nextBound, dofs);
204
      } else {
205
	addThisEdge = true;
206
      }
207 208 209 210

      break;
    case 2:
      if (child[0]) {
211 212
	if (bound.reverseMode) {
	  nextBound.ithObj = 1;
213
	  child[1]->getHigherOrderDofs(feSpace, nextBound, dofs);
214
	  nextBound.ithObj = 0;
215
	  child[0]->getHigherOrderDofs(feSpace, nextBound, dofs);
216 217
	} else {
	  nextBound.ithObj = 0;
218
	  child[0]->getHigherOrderDofs(feSpace, nextBound, dofs);
219
	  nextBound.ithObj = 1;
220
	  child[1]->getHigherOrderDofs(feSpace, nextBound, dofs);
221
	}
222 223 224 225 226 227 228 229 230 231
      } else {
	addThisEdge = true;
      } 

      break;
    default:
      ERROR_EXIT("Should never happen!\n");
    }

    if (addThisEdge) {
232
      DofContainer addDofs;
233 234 235 236
      ElementDofIterator elDofIter(feSpace, true);
      elDofIter.reset(this);
      do {
	if (elDofIter.getCurrentPos() == 1 && 
237
	    elDofIter.getCurrentElementPos() == bound.ithObj)
238
	  addDofs.push_back(elDofIter.getDofPtr());	
239
      } while(elDofIter.next());      
240 241 242 243 244 245 246 247


      if (bound.reverseMode)
 	for (int i = addDofs.size() - 1; i >= 0; i--)
 	  dofs.push_back(addDofs[i]);
      else
 	for (unsigned int i = 0; i < addDofs.size(); i++) 
 	  dofs.push_back(addDofs[i]);      
248 249 250 251
    }
  }


252
}