ElInfo1d.cc 9.43 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#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 {
16
17
18
19
20
21
22
23
24
25
26
27
  
  double ElInfo1d::mat_d1_val[2][2] = {{1.0, 0.0}, 
				       {0.0, 1.0}};
  mtl::dense2D<double> ElInfo1d::mat_d1(mat_d1_val);

  double ElInfo1d::mat_d1_left_val[2][2] = {{1.0, 0.5}, 
					    {0.0, 0.5}};
  mtl::dense2D<double> ElInfo1d::mat_d1_left(mat_d1_left_val);

  double ElInfo1d::mat_d1_right_val[2][2] = {{0.5, 1.0}, 
					     {0.5, 0.0}};
  mtl::dense2D<double> ElInfo1d::mat_d1_right(mat_d1_right_val);
28

29
30
31
32
33
34
35
36
37
38
39
40
41
42
  double ElInfo1d::mat_d2_val[3][3] = {{1.0, 0.0, 0.0},
				       {0.0, 1.0, 0.0},
				       {0.0, 0.0, 1.0}};
  mtl::dense2D<double> ElInfo1d::mat_d2(mat_d2_val);

  double ElInfo1d::mat_d2_left_val[3][3] = {{1.0, 0.375, 0.0},
					    {0.0, 0.75, 1.0},
					    {0.0, -0.125, 0.0}};
  mtl::dense2D<double> ElInfo1d::mat_d2_left(mat_d2_left_val);
  
  double ElInfo1d::mat_d2_right_val[3][3] = {{0.0, -0.125, 0.0},
					     {1.0, 0.75, 0.0},
					     {0.0, 0.375, 1.0}};
  mtl::dense2D<double> ElInfo1d::mat_d2_right(mat_d2_right_val);
43
44
45
46

  double ElInfo1d::test1_val[2][2] = {{1.0, 0.0},
				      {0.0, 1.0}};
  mtl::dense2D<double> ElInfo1d::test2_val(test1_val);
47
48
  

49
50
  void ElInfo1d::fillMacroInfo(const MacroElement * mel)
  {
51
    FUNCNAME("ElInfo1d::fillMacroInfo()");
52
53
54
    Element *nb;
    MacroElement *mnb;

Thomas Witkowski's avatar
Thomas Witkowski committed
55
56
57
    macroElement = const_cast<MacroElement*>( mel);
    element = const_cast<Element*>( mel->getElement());
    parent = NULL;
58
    level = 0;
59

Thomas Witkowski's avatar
Thomas Witkowski committed
60
    int vertices = mesh->getGeo(VERTEX);
61

Thomas Witkowski's avatar
Thomas Witkowski committed
62
63
    if (fillFlag.isSet(Mesh::FILL_COORDS) || fillFlag.isSet(Mesh::FILL_DET) ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
64
      
65
      for (int i = 0; i < vertices; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
66
	coord[i] = mel->coord[i];
67
    }
68

Thomas Witkowski's avatar
Thomas Witkowski committed
69
    if (fillFlag.isSet(Mesh::FILL_NEIGH) || fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {
70
71
      WorldVector<double> oppC;
      
Thomas Witkowski's avatar
Thomas Witkowski committed
72
      int neighbours =  mesh->getGeo(NEIGH);
73
74
75
      for (int i = 0; i < neighbours; i++) {
	nb = NULL;
	if ((mnb = const_cast<MacroElement*>( mel->getNeighbour(i)))) {
Thomas Witkowski's avatar
Thomas Witkowski committed
76
	  if (fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {
77
78
79
80
81
82
	    oppC = mnb->coord[i];
	  }

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

	  while (!(nb->isLeaf())) { // make nb nearest element
Thomas Witkowski's avatar
Thomas Witkowski committed
83
	    if (fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {
84
	      if (nb->isNewCoordSet()) {
85
86
87
		oppC = *(nb->getNewCoord());
	      } else {
		oppC = (mel->coord[i] + oppC) * 0.5;
88
	      }
89
90
91
92
	    }
	    nb = const_cast<Element*>( nb->getChild(1-i));
	  }

Thomas Witkowski's avatar
Thomas Witkowski committed
93
94
	  if (fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {
	    oppCoord[i] = oppC;
95
	  }
96
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
97
	neighbour[i] = nb;
98
	oppVertex[i] = nb ? i : -1;
99
      }
100
    }
101

Thomas Witkowski's avatar
Thomas Witkowski committed
102
    if (fillFlag.isSet(Mesh::FILL_BOUND) ) {
103
      for (int i = 0; i < vertices; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
104
	boundary[i] = mel->getBoundary(i);
105

Thomas Witkowski's avatar
Thomas Witkowski committed
106
107
      for (int i = 0; i < element->getGeo(PROJECTION); i++)
	projection[i] = mel->getProjection(i);
108
109
110
111
112
113
114
115
    }
  }

  /****************************************************************************/
  /*  compute gradients of basis functions on element; return the absulute    */
  /*  value of the determinante from the transformation to the reference      */
  /*  element                                                                 */
  /****************************************************************************/
Thomas Witkowski's avatar
Thomas Witkowski committed
116
  double ElInfo1d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
117
  {
118
    FUNCNAME("ElInfo1d::calcGrdLambda()");
119
120
121
122

    testFlag(Mesh::FILL_COORDS);

    WorldVector<double> e;
Thomas Witkowski's avatar
Thomas Witkowski committed
123
124
    e = coord[1]; 
    e -= coord[0];
125
    double adet2 = e * e;
126
127
128

    if (adet2 < 1.0E-15) {
      MSG("det*det = %lf\n", adet2);
129
      grd_lam[0] = grd_lam[1] = 0.0;
130
    } else {
131
      grd_lam[1] = e * (1.0 / adet2);
132
      grd_lam[0] = grd_lam[1] * -1.0;
133
134
135
136
137
138
139
140
    }

    return sqrt(adet2);
  }

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

143
    double lmin;
Thomas Witkowski's avatar
Thomas Witkowski committed
144
145
146
    double a = coord[0][0];
    double length = (coord[1][0] - a);
    int dim = mesh->getDim();
147
148
149

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

150
151
    TEST_EXIT_DBG(lambda)("lambda must not be NULL\n");
    TEST_EXIT_DBG(dim == 1)("dim!=1\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
152
    TEST_EXIT_DBG(dimOfWorld == dim)("not yet for DIM != DIM_OF_WORLD\n");
153
154
155
156
157
158

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

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

162
    int k = -1;
163
    lmin = 0.0;
164
    for (int i = 0; i <= dim; i++) {
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
      if ((*lambda)[i] < -1.E-5) {
	if ((*lambda)[i] < lmin) {
	  k = i;
	  lmin = (*lambda)[i];
	}
      }
    }

    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                                 */
  /****************************************************************************/
Thomas Witkowski's avatar
Thomas Witkowski committed
181
  double ElInfo1d::getNormal(int side, WorldVector<double> &normal)
182
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
183
    normal = coord[side] - coord[(side + 1) % 2];
184
    double det = norm(&normal);
185
    TEST_EXIT_DBG(det > 1.e-30)("det = 0 on side %d\n", side);
186
187
    normal *= 1.0 / det;

188
    return det;
189
190
191
192
193
194
195
196
197
198
  }


  /****************************************************************************/
  /*  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
  {
199
200
    FUNCNAME("ElInfo::getElementNormal()");

Thomas Witkowski's avatar
Thomas Witkowski committed
201
202
    TEST_EXIT_DBG(dimOfWorld == 2)
      (" element normal only well defined for  DIM_OF_WORLD = DIM + 1 !!");
203

Thomas Witkowski's avatar
Thomas Witkowski committed
204
205
    elementNormal[0] = coord[1][1] - coord[0][1];
    elementNormal[1] = coord[0][0] - coord[1][0];
206

Thomas Witkowski's avatar
Thomas Witkowski committed
207
    double det = norm(&elementNormal);
208

209
    TEST_EXIT_DBG(det > 1.e-30)("det = 0");
210
211
212

    elementNormal *= 1.0 / det;
    
213
    return det;
214
215
216
  }


217
  void ElInfo1d::fillElInfo(int ichild, const ElInfo *elInfoOld)
218
  {
219
220
    FUNCNAME("ElInfo1d::fillElInfo()");

221
    Element *nb;
Thomas Witkowski's avatar
Thomas Witkowski committed
222
    Element *elem = elInfoOld->element;
223

224
    TEST_EXIT_DBG(elem->getChild(0))("no children?\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
225
    element = const_cast<Element*>(elem->getChild(ichild));
226

Thomas Witkowski's avatar
Thomas Witkowski committed
227
    TEST_EXIT_DBG(element)("missing child %d?\n", ichild);
228

Thomas Witkowski's avatar
Thomas Witkowski committed
229
230
231
    macroElement = elInfoOld->macroElement;
    fillFlag = elInfoOld->fillFlag;
    parent = elem;
232
    level = elInfoOld->level + 1;
233
    iChild = ichild;
234

Thomas Witkowski's avatar
Thomas Witkowski committed
235
    int neighbours = mesh->getGeo(NEIGH);
236

Thomas Witkowski's avatar
Thomas Witkowski committed
237
238
    if (fillFlag.isSet(Mesh::FILL_COORDS) || fillFlag.isSet(Mesh::FILL_DET) ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
239

Thomas Witkowski's avatar
Thomas Witkowski committed
240
      const FixVec<WorldVector<double>, VERTEX> *old_coord = &(elInfoOld->coord);
241

Thomas Witkowski's avatar
Thomas Witkowski committed
242
      coord[ichild] = (*old_coord)[ichild];
243
      if (elem->isNewCoordSet())
Thomas Witkowski's avatar
Thomas Witkowski committed
244
	coord[1 - ichild] = *(elem->getNewCoord());
245
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
246
	coord[1 - ichild] = ((*old_coord)[0] + (*old_coord)[1]) * 0.5;
247
    }
248

Thomas Witkowski's avatar
Thomas Witkowski committed
249
    if (fillFlag.isSet(Mesh::FILL_NEIGH) || fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {
250
251
      WorldVector<double> oppC;
      
Thomas Witkowski's avatar
Thomas Witkowski committed
252
      TEST_EXIT_DBG(fillFlag.isSet(Mesh::FILL_COORDS))
253
254
255
256
	("FILL_OPP_COORDS only with FILL_COORDS\n");
      
      for (int i = 0; i < neighbours; i++) {
	if (i != ichild) {
257
	  nb = const_cast<Element*>(elem->getChild(1-ichild));  
Thomas Witkowski's avatar
Thomas Witkowski committed
258
259
	  if (fillFlag.isSet(Mesh::FILL_OPP_COORDS))
	    oppC = elInfoOld->coord[i];
260
	} else {
261
	  nb = const_cast<Element*>( elInfoOld->getNeighbour(i));
262

Thomas Witkowski's avatar
Thomas Witkowski committed
263
264
	  if (nb && fillFlag.isSet(Mesh::FILL_OPP_COORDS))
	    oppC = elInfoOld->oppCoord[i];
265
	}
266

267
268
	if (nb) {
	  while (nb->getChild(0)) {  // make nb nearest element
Thomas Witkowski's avatar
Thomas Witkowski committed
269
	    if (fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {
270
	      if (nb->isNewCoordSet())
271
		oppC = *(nb->getNewCoord());
272
	      else
Thomas Witkowski's avatar
Thomas Witkowski committed
273
		oppC = (coord[i] + oppC) * 0.5;
274
275
276
	    }
	    nb = const_cast<Element*>( nb->getChild(1-i));
	  }
277

Thomas Witkowski's avatar
Thomas Witkowski committed
278
279
	  if (fillFlag.isSet(Mesh::FILL_OPP_COORDS))
	    oppCoord[i] = oppC;
280
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
281
	neighbour[i] = nb;
282
	oppVertex[i] = nb ? i : -1;
283
      }
284
    }
285

Thomas Witkowski's avatar
Thomas Witkowski committed
286
287
288
    if (fillFlag.isSet(Mesh::FILL_BOUND)) {
      boundary[ichild] = elInfoOld->getBoundary(ichild);
      boundary[1 - ichild] = INTERIOR;
289

290
291
      if (elInfoOld->getProjection(0) && 
	  elInfoOld->getProjection(0)->getType() == VOLUME_PROJECTION) {
Thomas Witkowski's avatar
Thomas Witkowski committed
292
	projection[0] = elInfoOld->getProjection(0);
293
      }
294
    }   
295
296
  }

297
298
  void ElInfo1d::getSubElemCoordsMat(mtl::dense2D<double>& mat, 
				     int basisFctsDegree) const
299
  {
300
    switch (basisFctsDegree) {
301
    case 1:
302
303
304
305
306
307
308
309
310
      mat = mat_d1;

      for (int i = 0; i < refinementPathLength; i++) {
	if (refinementPath & (1 << i)) 
	  mat *= mat_d1_left;
	else 
	  mat *= mat_d1_right;
      }

311
      break;
312
      
313
    case 2:
314
315
316
317
318
319
320
321
322
      mat = mat_d2;

      for (int i = 0; i < refinementPathLength; i++) {
	if (refinementPath & (1 << i)) 
	  mat *= mat_d2_left;
	else 
	  mat *= mat_d2_right;
      }

323
      break;
324

325
    default:
326
      ERROR_EXIT("Not supported for basis function degree: %d\n", basisFctsDegree);
327
    }
328
329
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
330
331
  void ElInfo1d::getSubElemCoordsMat_so(mtl::dense2D<double>& mat, 
					int basisFctsDegree) const
332
  {
333
    switch (basisFctsDegree) {
334
    case 1:
335
      mat = test2_val;
336

337
338
      for (int i = 0; i < refinementPathLength; i++)
	mat *= 0.5;
339

340
341
      break;
    default:
342
      ERROR_EXIT("Not supported for basis function degree: %d\n", basisFctsDegree);
343
344
345
    }
  }

346

347
}