ElInfo.cc 4.56 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
#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 {

  int ElInfo::traverseId = -1;
  int ElInfo::traverseIdCounter = 0;

  ElInfo::ElInfo(Mesh *aMesh) 
    : mesh_(aMesh),
      element_(NULL),
      parent_(NULL),
      macroElement_(NULL),
      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
32
33
34
35
36
      oppVertex_(mesh_->getDim(), NO_INIT),
      grdLambda_(mesh_->getDim(), NO_INIT)
		//parametric_(false)
  
  {
    projection_.set(NULL);
37
38
39
40

    for (int i = 0; i < neighbourCoord_.getSize(); i++) {
      neighbourCoord_[i].init(mesh_->getDim());
    }
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
  }

  ElInfo::~ElInfo()
  {
  }


  /****************************************************************************/
  /*  transform local coordintes l to world coordinates; if w is non NULL      */
  /*  store them at w otherwise return a pointer to some local static         */
  /*  area containing world coordintes                                        */
  /****************************************************************************/
  const WorldVector<double> *ElInfo::coordToWorld(const DimVec<double>& l, 
						  WorldVector<double>* w) const

  {
    return coordToWorld(l, &coord_, w);
  }

  const WorldVector<double> *ElInfo::coordToWorld(const DimVec<double>& l,
						  const FixVec<WorldVector<double>, VERTEX> *coords,
						  WorldVector<double> *w)

  {
    int dim = l.getSize() - 1;
    int dimOfWorld = Global::getGeo(WORLD);

68
69
    static WorldVector<double> world[2];
    WorldVector<double> *ret = w ? w : &world[omp_get_thread_num()];
70
71
    WorldVector<double> v = (*coords)[0];
    double c = l[0];
72

73
    for (int j = 0; j < dimOfWorld; j++)
74
      (*ret)[j] = c * v[j];
75
76
77
  
    int vertices =  Global::getGeo(VERTEX, dim);

78
    for (int i = 1; i < vertices; i++) {
79
80
      v = (*coords)[i];
      c = l[i];
81
      for (int j = 0; j < dimOfWorld; j++)
82
	(*ret)[j] += c * v[j];
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
    }
    return ret;
  }

  /****************************************************************************/
  /*  compute volume of an element                                            */
  /****************************************************************************/

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

  double ElInfo::calcDet(const FixVec<WorldVector<double>, VERTEX> &coords)
  {
99
    FUNCNAME("ElInfo::calcDet()");
100
101
102
103

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

104
    TEST_EXIT_DBG(dim <= dow)("dim > dow\n");
105
106
107
108
109
110
111
112
113
114
115

    double det = 0.0;

    if (dim == 0)
      return 1.0;

    // dim = dim of world
    WorldVector<double> e1, e2, e3, v0;

    v0 = coords[0];

116
    for (int i = 0; i < dow; i++) {
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
      e1[i] = coords[1][i] - v0[i];
      if(dim > 1)
	e2[i] = coords[2][i] - v0[i];
      if(dim > 2)
	e3[i] = coords[3][i] - v0[i];
    }
  
    switch(dow) {
    case 1:
      det = coords[1][0] - coords[0][0];
      break;
    case 2:
      if(dim == 1) {
	det = norm(&e1);
      } else {
	det = e1[0]*e2[1] - e1[1]*e2[0];
      }
      break;
    case 3: 
      {
	WorldVector<double> n;

139
	if (dim > 1) {
140
141
142
143
144
	  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];
	}

145
	if (dim == 1) {
146
147
148
149
	  det = norm(&e1);
	} else if (dim == 2) {
	  det = norm(&n);
	} else if (dim == 3) {
150
	  det = n[0] * e3[0] + n[1] * e3[1] + n[2] * e3[2];
151
152
153
154
155
156
157
158
159
160
161
162
163
	} else
	  ERROR_EXIT("not yet for problem dimension = %d",dim);
	break; 
      }
    default:
      ERROR_EXIT("not yet for Global::getGeo(WORLD) = %d",dow);
    }

    return abs(det);
  }

  void ElInfo::testFlag(const Flag& flags) const
  {
164
    TEST_EXIT_DBG(fillFlag_.isSet(flags))("flag not set\n");
165
166
167
168
169
  }


  void ElInfo::fillDetGrdLambda() 
  { 
170
    if (fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
171
172
      det_ = calcGrdLambda(grdLambda_);
    } else {
173
      if (fillFlag_.isSet(Mesh::FILL_DET)) {
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
	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);
189
    int offset = indexOffset[dim - 1][posIndex];
190
191
192
193
194
  
    return boundary_[offset + i];
  }

}