ElInfo.cc 4.59 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
18
19
20
21
22
23
24
25
26
27
28
#include "ElInfo.h"
#include "BasisFunction.h"
#include "Element.h"
#include "Line.h"
#include "Triangle.h"
#include "Tetrahedron.h"
#include "FiniteElemSpace.h"
#include "Flag.h"
#include "MacroElement.h"
#include "Mesh.h"
#include "Global.h"
#include "FixVec.h"
#include "DOFVector.h"

namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
29
  std::vector<std::map<std::pair<int, unsigned long>, mtl::dense2D<double> > > ElInfo::subElemMatrices(5);
30

Thomas Witkowski's avatar
Thomas Witkowski committed
31
  std::vector<std::map<std::pair<int, unsigned long>, mtl::dense2D<double> > > ElInfo::subElemGradMatrices(5);
32

33
  ElInfo::ElInfo(Mesh *aMesh) 
Thomas Witkowski's avatar
Thomas Witkowski committed
34
35
36
37
    : mesh(aMesh),
      element(NULL),
      parent(NULL),
      macroElement(NULL),
38
      level(0),
39
      elType(0),
40
      iChild(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
41
42
43
44
45
46
47
48
      coord(mesh->getDim(), NO_INIT),
      boundary(mesh->getDim(), DEFAULT_VALUE, INTERIOR),
      projection(mesh->getDim(), NO_INIT),
      oppCoord(mesh->getDim(), NO_INIT),
      neighbour(mesh->getDim(), NO_INIT),
      neighbourCoord(mesh->getDim(), NO_INIT),
      oppVertex(mesh->getDim(), NO_INIT),
      grdLambda(mesh->getDim(), NO_INIT),
49
50
      refinementPath(0),
      refinementPathLength(0)
51
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
52
    projection.set(NULL);
53

54
    for (int i = 0; i < neighbourCoord.getSize(); i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
55
      neighbourCoord[i].init(mesh->getDim());
Thomas Witkowski's avatar
Thomas Witkowski committed
56
57

    dimOfWorld = Global::getGeo(WORLD);
58
  }
59
 
60
61

  ElInfo::~ElInfo()
62
  {}
63

64

65
66
  void ElInfo::coordToWorld(const DimVec<double>& l, 
			    WorldVector<double>& w) const
67
68

  {
69
70
    testFlag(Mesh::FILL_COORDS);

71
    double c = l[0];
72

73
    for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
74
      w[j] = c * coord[0][j];
75
  
76
    int vertices = Global::getGeo(VERTEX, l.getSize() - 1);
77

78
    for (int i = 1; i < vertices; i++) {
79
      c = l[i];
80
      for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
81
	w[j] += c * coord[i][j];
82
83
84
    }
  }

85

86
87
88
  double ElInfo::calcDet() const
  {
    testFlag(Mesh::FILL_COORDS);
Thomas Witkowski's avatar
Thomas Witkowski committed
89
    return calcDet(coord);
90
91
  }

92

Thomas Witkowski's avatar
Thomas Witkowski committed
93
  double ElInfo::calcDet(const FixVec<WorldVector<double>, VERTEX> &coords) const
94
  {
95
    FUNCNAME("ElInfo::calcDet()");
96
97
98
99

    int dim = coords.getSize() - 1;
    int dow = Global::getGeo(WORLD);

100
    TEST_EXIT_DBG(dim <= dow)("dim > dow\n");
101
102
103
104
105
106

    double det = 0.0;

    if (dim == 0)
      return 1.0;

Thomas Witkowski's avatar
Thomas Witkowski committed
107
    switch (dow) {
108
109
110
111
    case 1:
      det = coords[1][0] - coords[0][0];
      break;
    case 2:
112
      if (dim == 1) {
Thomas Witkowski's avatar
Thomas Witkowski committed
113
114
115
116
117

	WorldVector<double> e1;

	e1[0] = coords[1][0] - coords[0][0];
	e1[1] = coords[1][1] - coords[0][1];
118
	det = norm(&e1);
Thomas Witkowski's avatar
Thomas Witkowski committed
119

120
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
121
122
	det = (coords[1][0] - coords[0][0]) * (coords[2][1] - coords[0][1]) - 
	  (coords[1][1] - coords[0][1]) * (coords[2][0] - coords[0][0]);
123
124
125
126
      }
      break;
    case 3: 
      {
Thomas Witkowski's avatar
Thomas Witkowski committed
127
	WorldVector<double> e1, e2, e3, n;
128
	
Thomas Witkowski's avatar
Thomas Witkowski committed
129
130
131
132
133
134
135
136
	for (int i = 0; i < dow; i++) {
	  e1[i] = coords[1][i] - coords[0][i];
	  if (dim > 1)
	    e2[i] = coords[2][i] - coords[0][i];
	  if (dim > 2)
	    e3[i] = coords[3][i] - coords[0][i];
	}
  
137
	if (dim > 1) {
Thomas Witkowski's avatar
Thomas Witkowski committed
138
139
140
	  n[0] = e1[1] * e2[2] - e1[2] * e2[1];
	  n[1] = e1[2] * e2[0] - e1[0] * e2[2];
	  n[2] = e1[0] * e2[1] - e1[1] * e2[0];
141
	}
142
	
143
	if (dim == 1) {
144
145
146
147
	  det = norm(&e1);
	} else if (dim == 2) {
	  det = norm(&n);
	} else if (dim == 3) {
148
	  det = n[0] * e3[0] + n[1] * e3[1] + n[2] * e3[2];
149
	} else
Thomas Witkowski's avatar
Thomas Witkowski committed
150
	  ERROR_EXIT("not yet for problem dimension = %d", dim);
151
152
	break;
      } 
153
    default:
Thomas Witkowski's avatar
Thomas Witkowski committed
154
      ERROR_EXIT("not yet for Global::getGeo(WORLD) = %d", dow);
155
    }
156
    
157
158
159
    return abs(det);
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
160

161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
  double ElInfo::calcSurfaceDet(VectorOfFixVecs<DimVec<double> > &surfVert) const
  {
    double surfDet;
    int dim = surfVert[0].getSize() - 1;
    FixVec<WorldVector<double>, VERTEX> worldCoords(dim - 1, NO_INIT);
  
    // transform barycentric coords to world coords
    for (int i = 0; i < dim; i++) {
      coordToWorld(surfVert[i], worldCoords[i]);
    }
  
    // calculate determinant for surface
    surfDet = calcDet(worldCoords);

    return surfDet;
  }


179
  void ElInfo::testFlag(const Flag& flags) const
Thomas Witkowski's avatar
Thomas Witkowski committed
180
  {    
181
    TEST_EXIT_DBG(fillFlag.isSet(flags))("flag not set\n");
182
183
184
185
186
  }


  void ElInfo::fillDetGrdLambda() 
  { 
Thomas Witkowski's avatar
Thomas Witkowski committed
187
188
    if (fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
      det = calcGrdLambda(grdLambda);
189
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
190
191
      if (fillFlag.isSet(Mesh::FILL_DET))
	det = calcDet();
192
193
194
    }
  }

195

196
197
198
199
200
201
202
203
  BoundaryType ElInfo::getBoundary(GeoIndex pos, int i)
  {
    static int indexOffset[3][3] = {
      { 0, 0, 0},
      { 3, 0, 0},
      {10, 4, 0}
    };

Thomas Witkowski's avatar
Thomas Witkowski committed
204
    int dim = mesh->getDim();
205
    int posIndex = DIM_OF_INDEX(pos, dim);
206
    int offset = indexOffset[dim - 1][posIndex];
207
    
Thomas Witkowski's avatar
Thomas Witkowski committed
208
    return boundary[offset + i];
209
210
211
  }

}