ElInfo1d.cc 7.85 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#include "ElInfo1d.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 {

  void ElInfo1d::fillMacroInfo(const MacroElement * mel)
  {
19
    FUNCNAME("ElInfo1d::fillMacroInfo()");
20
21
22
    Element *nb;
    MacroElement *mnb;

23
24
25
26
    macroElement_ = const_cast<MacroElement*>( mel);
    element_ = const_cast<Element*>( mel->getElement());
    parent_ = NULL;
    level_  = 0;
27
28
29
30

    int vertices = mesh_->getGeo(VERTEX);

    if (fillFlag_.isSet(Mesh::FILL_COORDS) || fillFlag_.isSet(Mesh::FILL_DET) ||
31
32
33
34
35
	fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
      
      for (int i = 0; i < vertices; i++) {
	coord_[i] = mel->coord[i];
      }
36
37
38
39
40
41
      }

    //   if(fillFlag_.isSet(Mesh::FILL_DET) || fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
    //     det = elGrdLambda(*Lambda);
    //   }

42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
    if (fillFlag_.isSet(Mesh::FILL_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
      WorldVector<double> oppC;
      
      int neighbours =  mesh_->getGeo(NEIGH);
      for (int i = 0; i < neighbours; i++) {
	nb = NULL;
	if ((mnb = const_cast<MacroElement*>( mel->getNeighbour(i)))) {
	  if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
	    oppC = mnb->coord[i];
	  }

	  nb = const_cast<Element*>( mnb->getElement());

	  while (!(nb->isLeaf())) { // make nb nearest element
	    if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
	      if (nb->getNewCoord(-1)) {
		oppC = *(nb->getNewCoord());
	      } else {
		oppC = (mel->coord[i] + oppC) * 0.5;
61
	      }
62
63
64
65
66
67
	    }
	    nb = const_cast<Element*>( nb->getChild(1-i));
	  }

	  if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
	    oppCoord_[i] = oppC;
68
	  }
69
70
71
	}
	neighbour_[i] = nb;
	oppVertex_[i] = nb ? i : -1;
72
      }
73
    }
74
75

    if (fillFlag_.isSet(Mesh::FILL_BOUND) ) {
76
      for (int i = 0; i < vertices; i++)
77
78
	boundary_[i] = mel->getBoundary(i);

79
      for (int i = 0; i < element_->getGeo(PROJECTION); i++) {
80
81
82
83
84
85
86
87
88
89
90
91
	projection_[i] = mel->getProjection(i);
      }
    }
  }

  /****************************************************************************/
  /*  compute gradients of basis functions on element; return the absulute    */
  /*  value of the determinante from the transformation to the reference      */
  /*  element                                                                 */
  /****************************************************************************/
  double ElInfo1d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam) const
  {
92
    FUNCNAME("ElInfo1d::calcGrdLambda()");
93
94
95
96

    testFlag(Mesh::FILL_COORDS);

    WorldVector<double> e;
97
98
99
    e = coord_[1]; 
    e -= coord_[0];
    double adet2 = e * e;
100
101
102

    if (adet2 < 1.0E-15) {
      MSG("det*det = %lf\n", adet2);
103
      for (int i = 0; i <= 1; i++)
104
105
	grd_lam[i] = 0.0;
    } else {
106
      grd_lam[1] = e * (1.0 / adet2);
107
108
109
110
111
112
113
114
115
      grd_lam[0] = grd_lam[1] * (-1.0);
    }

    return sqrt(adet2);
  }

  const int ElInfo1d::worldToCoord(const WorldVector<double>& x,
				   DimVec<double>* lambda) const
  {
116
117
    FUNCNAME("ElInfo1d::worldToCoord()");

118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
    double lmin;
    double a = coord_[0][0];
    double length = (coord_[1][0] - a);
    int dim = mesh_->getDim();

    static DimVec<double> vec(dim, NO_INIT);

    TEST_EXIT(lambda)("lambda must not be NULL\n");
    TEST_EXIT(dim==1)("dim!=1\n");
    TEST_EXIT(Global::getGeo(WORLD)==dim)("not yet for DIM != DIM_OF_WORLD\n");

    if (abs(length) < DBL_TOL) {
      ERROR_EXIT("length = %le; abort\n", length);
      return 0;
    }

    (*lambda)[1] = (x[0]-a)/length;
    (*lambda)[0] = 1.0 - (*lambda)[1];

137
    int k = -1;
138
    lmin = 0.0;
139
140
    int j = 0;
    for (int i = 0; i <= dim; i++) {
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
      if ((*lambda)[i] < -1.E-5) {
	if ((*lambda)[i] < lmin) {
	  k = i;
	  lmin = (*lambda)[i];
	}
	j++;
      }
    }

    return k;
  }

  /****************************************************************************/
  /*  calculate a facenormal of edge side of a triangle with coordinates      */
  /*  coord; return the absulute value of the determinant from the           */
  /*  transformation to the reference element                                 */
  /****************************************************************************/
  double ElInfo1d::getNormal(int side, WorldVector<double> &normal) const
  {
    normal = coord_[side] - coord_[(side + 1) % 2];
161
    double det = norm(&normal);
162
    TEST_EXIT(det > 1.e-30)("det = 0 on side %d\n", side);
163
164
    normal *= 1.0 / det;

165
166
167
168
169
170
171
172
173
174
175
    return(det);
  }


  /****************************************************************************/
  /*  calculate the normal of the element for dim of world = 2                */
  /*  return the absulute value of the determinant from the                   */
  /*  transformation to the reference element                                 */
  /****************************************************************************/
  double ElInfo1d::getElementNormal( WorldVector<double> &elementNormal) const
  {
176
177
178
    FUNCNAME("ElInfo::getElementNormal()");

    double det = 0.0;
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
    int dow = Global::getGeo(WORLD);

    TEST_EXIT(dow =  2)(" element normal only well defined for  DIM_OF_WORLD = DIM + 1 !!");

    elementNormal[0] = coord_[1][1] - coord_[0][1];
    elementNormal[1] = coord_[0][0] - coord_[1][0];

    det = norm(&elementNormal);

    TEST_EXIT(det > 1.e-30)("det = 0");

    elementNormal *= 1.0 / det;
    
    return(det);
  }


  void ElInfo1d::fillElInfo(int ichild, const ElInfo *elinfo_old)
  {
198
199
    FUNCNAME("ElInfo1d::fillElInfo()");

200
201
202
203
204
205
206
    Element *nb;
    Element *elem = elinfo_old->element_;

    TEST_EXIT(elem->getChild(0))("no children?\n");
    TEST_EXIT((element_ = const_cast<Element*>( elem->getChild(ichild))))
      ("missing child %d?\n", ichild);

207
208
209
210
    macroElement_ = elinfo_old->macroElement_;
    fillFlag_ = elinfo_old->fillFlag_;
    parent_ = elem;
    level_ = elinfo_old->level_ + 1;
211
212
213
214

    int neighbours = mesh_->getGeo(NEIGH);

    if (fillFlag_.isSet(Mesh::FILL_COORDS) || fillFlag_.isSet(Mesh::FILL_DET) ||
215
	fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
216

217
218
219
220
221
222
223
      const FixVec<WorldVector<double>, VERTEX> *old_coord = &(elinfo_old->coord_);

      coord_[ichild] = (*old_coord)[ichild];
      if (elem->getNewCoord(-1)) {
	coord_[1 - ichild] = *(elem->getNewCoord());
      } else {
	coord_[1 - ichild] = ((*old_coord)[0] + (*old_coord)[1]) * 0.5;
224
      }
225
    }
226
227
228
229
230

    //   if(fillFlag_.isSet(Mesh::FILL_DET) || fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
    //     det = calcGrdLambda(*Lambda);
    //   }

231
232
233
234
235
236
237
238
239
240
241
242
243
244
    if (fillFlag_.isSet(Mesh::FILL_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
      WorldVector<double> oppC;
      
      TEST_EXIT(fillFlag_.isSet(Mesh::FILL_COORDS))
	("FILL_OPP_COORDS only with FILL_COORDS\n");
      
      for (int i = 0; i < neighbours; i++) {
	if (i != ichild) {
	  nb = const_cast<Element*>( elem->getChild(1-ichild));  
	  if ( fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
	    oppC = elinfo_old->coord_[i];
	  }
	} else {
	  nb = const_cast<Element*>( elinfo_old->getNeighbour(i));
245

246
247
	  if (nb && fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
	    oppC = elinfo_old->oppCoord_[i];
248
	  }
249
	}
250

251
252
253
254
255
256
257
258
259
260
261
	if (nb) {
	  while (nb->getChild(0)) {  // make nb nearest element
	    if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
	      if (nb->getNewCoord(-1)) {
		oppC = *(nb->getNewCoord());
	      } else {
		oppC = (coord_[i] + oppC) * 0.5;
	      }
	    }
	    nb = const_cast<Element*>( nb->getChild(1-i));
	  }
262

263
264
	  if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
	    oppCoord_[i] = oppC;
265
	  }
266
267
268
	}
	neighbour_[i] = nb;
	oppVertex_[i] = nb ? i : -1;
269
      }
270
    }
271

272
273
274
275
276
277
278
279
280
281
    if (fillFlag_.isSet(Mesh::FILL_BOUND)) {
      boundary_[ichild] = elinfo_old->getBoundary(ichild);
      boundary_[1-ichild] = INTERIOR;

      if (elinfo_old->getProjection(0) && 
	  elinfo_old->getProjection(0)->getType() == VOLUME_PROJECTION) {
	projection_[0] = elinfo_old->getProjection(0);
      }
    }
    
282
283
284
285
    return;
  }

}