ElInfo.cc 3.68 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
#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 {

  ElInfo::ElInfo(Mesh *aMesh) 
    : mesh_(aMesh),
      element_(NULL),
      parent_(NULL),
      macroElement_(NULL),
22
      level(0),
23
      elType(0),
24
      iChild(0),
25
26
27
28
29
      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),
30
      neighbourCoord(mesh_->getDim(), NO_INIT),
31
      oppVertex(mesh_->getDim(), NO_INIT),
Thomas Witkowski's avatar
Thomas Witkowski committed
32
      grdLambda(mesh_->getDim(), NO_INIT),
33
34
      refinementPath(0),
      refinementPathLength(0)
35
36
  {
    projection_.set(NULL);
37

38
39
    for (int i = 0; i < neighbourCoord.getSize(); i++)
      neighbourCoord[i].init(mesh_->getDim());
Thomas Witkowski's avatar
Thomas Witkowski committed
40
41

    dimOfWorld = Global::getGeo(WORLD);
42
43
44
  }

  ElInfo::~ElInfo()
45
  {}
46

47
48
  void ElInfo::coordToWorld(const DimVec<double>& l, 
			    WorldVector<double>& w) const
49
50

  {
51
52
    testFlag(Mesh::FILL_COORDS);

53
    double c = l[0];
54

55
    for (int j = 0; j < dimOfWorld; j++)
56
      w[j] = c * coord_[0][j];
57
  
58
    int vertices = Global::getGeo(VERTEX, l.getSize() - 1);
59

60
    for (int i = 1; i < vertices; i++) {
61
      c = l[i];
62
      for (int j = 0; j < dimOfWorld; j++)
63
	w[j] += c * coord_[i][j];
64
65
66
67
68
69
70
71
72
    }
  }

  double ElInfo::calcDet() const
  {
    testFlag(Mesh::FILL_COORDS);
    return calcDet(coord_);
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
73
  double ElInfo::calcDet(const FixVec<WorldVector<double>, VERTEX> &coords) const
74
  {
75
    FUNCNAME("ElInfo::calcDet()");
76
77
78
79

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

80
    TEST_EXIT_DBG(dim <= dow)("dim > dow\n");
81
82
83
84
85
86

    double det = 0.0;

    if (dim == 0)
      return 1.0;

Thomas Witkowski's avatar
Thomas Witkowski committed
87
    switch (dow) {
88
89
90
91
    case 1:
      det = coords[1][0] - coords[0][0];
      break;
    case 2:
92
      if (dim == 1) {
Thomas Witkowski's avatar
Thomas Witkowski committed
93
94
95
96
97

	WorldVector<double> e1;

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

100
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
101
102
103
104

	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]);

105
106
107
108
      }
      break;
    case 3: 
      {
Thomas Witkowski's avatar
Thomas Witkowski committed
109
	WorldVector<double> e1, e2, e3, n;
110
	
Thomas Witkowski's avatar
Thomas Witkowski committed
111
112
113
114
115
116
117
118
	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];
	}
  
119
	if (dim > 1) {
Thomas Witkowski's avatar
Thomas Witkowski committed
120
121
122
	  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];
123
	}
124
	
125
	if (dim == 1) {
126
127
128
129
	  det = norm(&e1);
	} else if (dim == 2) {
	  det = norm(&n);
	} else if (dim == 3) {
130
	  det = n[0] * e3[0] + n[1] * e3[1] + n[2] * e3[2];
131
	} else
Thomas Witkowski's avatar
Thomas Witkowski committed
132
	  ERROR_EXIT("not yet for problem dimension = %d", dim);
133
134
	break;
      } 
135
    default:
Thomas Witkowski's avatar
Thomas Witkowski committed
136
      ERROR_EXIT("not yet for Global::getGeo(WORLD) = %d", dow);
137
    }
138
    
139
140
141
142
    return abs(det);
  }

  void ElInfo::testFlag(const Flag& flags) const
Thomas Witkowski's avatar
Thomas Witkowski committed
143
  {    
144
    TEST_EXIT_DBG(fillFlag_.isSet(flags))("flag not set\n");
145
146
147
148
149
  }


  void ElInfo::fillDetGrdLambda() 
  { 
150
    if (fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
151
      det_ = calcGrdLambda(grdLambda);
152
    } else {
153
      if (fillFlag_.isSet(Mesh::FILL_DET))
154
155
156
157
158
159
160
161
162
163
164
165
166
167
	det_ = calcDet();
    }
  }

  BoundaryType ElInfo::getBoundary(GeoIndex pos, int i)
  {
    static int indexOffset[3][3] = {
      { 0, 0, 0},
      { 3, 0, 0},
      {10, 4, 0}
    };

    int dim = mesh_->getDim();
    int posIndex = DIM_OF_INDEX(pos, dim);
168
    int offset = indexOffset[dim - 1][posIndex];
169
    
170
171
172
173
    return boundary_[offset + i];
  }

}