ElInfo1d.cc 10.5 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
#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 {
28
29
30
31
32
33
34
35
36
  
  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);

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

41
42
43
44
45
46
47
48
49
50
51
52
53
54
  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);
55
 
56

57
58
59
60
61
  void ElInfo1d::fillMacroInfo(const MacroElement * mel)
  {
    Element *nb;
    MacroElement *mnb;

Thomas Witkowski's avatar
Thomas Witkowski committed
62
63
    macroElement = const_cast<MacroElement*>(mel);
    element = const_cast<Element*>(mel->getElement());
Thomas Witkowski's avatar
Thomas Witkowski committed
64
    parent = NULL;
65
    level = 0;
66

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

Thomas Witkowski's avatar
Thomas Witkowski committed
69
70
    if (fillFlag.isSet(Mesh::FILL_COORDS) || fillFlag.isSet(Mesh::FILL_DET) ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
71
      
72
      for (int i = 0; i < vertices; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
73
	coord[i] = mel->coord[i];
74
    }
75

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

Thomas Witkowski's avatar
Thomas Witkowski committed
87
	  nb = const_cast<Element*>(mnb->getElement());
88
89

	  while (!(nb->isLeaf())) { // make nb nearest element
Thomas Witkowski's avatar
Thomas Witkowski committed
90
	    if (fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {
91
	      if (nb->isNewCoordSet()) {
92
93
94
		oppC = *(nb->getNewCoord());
	      } else {
		oppC = (mel->coord[i] + oppC) * 0.5;
95
	      }
96
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
97
	    nb = const_cast<Element*>(nb->getChild(1-i));
98
99
	  }

Thomas Witkowski's avatar
Thomas Witkowski committed
100
101
	  if (fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {
	    oppCoord[i] = oppC;
102
	  }
103
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
104
	neighbour[i] = nb;
105
	oppVertex[i] = nb ? i : -1;
106
      }
107
    }
108

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

Thomas Witkowski's avatar
Thomas Witkowski committed
113
114
      for (int i = 0; i < element->getGeo(PROJECTION); i++)
	projection[i] = mel->getProjection(i);
115
116
117
118
119
120
121
122
    }
  }

  /****************************************************************************/
  /*  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
123
  double ElInfo1d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
124
  {
125
    FUNCNAME("ElInfo1d::calcGrdLambda()");
126
127
128
129

    testFlag(Mesh::FILL_COORDS);

    WorldVector<double> e;
Thomas Witkowski's avatar
Thomas Witkowski committed
130
131
    e = coord[1]; 
    e -= coord[0];
132
    double adet2 = e * e;
133
134

    if (adet2 < 1.0E-15) {
135
      MSG("det*det = %g\n", adet2);
136
      grd_lam[0] = grd_lam[1] = 0.0;
137
    } else {
138
      grd_lam[1] = e * (1.0 / adet2);
139
      grd_lam[0] = grd_lam[1] * -1.0;
140
141
142
143
144
    }

    return sqrt(adet2);
  }

145
  int ElInfo1d::worldToCoord(const WorldVector<double>& x,
146
147
				   DimVec<double>* lambda) const
  {
148
149
    FUNCNAME("ElInfo1d::worldToCoord()");

150
    double lmin;
Thomas Witkowski's avatar
Thomas Witkowski committed
151
152
153
    double a = coord[0][0];
    double length = (coord[1][0] - a);
    int dim = mesh->getDim();
154
155
156

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

157
158
    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
159
    TEST_EXIT_DBG(dimOfWorld == dim)("not yet for DIM != DIM_OF_WORLD\n");
160
161
162
163
164
165

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

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

169
    int k = -1;
170
    lmin = 0.0;
171
    for (int i = 0; i <= dim; i++) {
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
      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                                 */
  /****************************************************************************/
188
  double ElInfo1d::getNormal(int side, WorldVector<double> &normal) const
189
  {
190
    FUNCNAME_DBG("ElInfo::getNormal()");
Thomas Witkowski's avatar
Thomas Witkowski committed
191
    normal = coord[side] - coord[(side + 1) % 2];
192
    double det = norm(&normal);
193
    TEST_EXIT_DBG(det > 1.e-30)("det = 0 on side %d\n", side);
194
195
    normal *= 1.0 / det;

196
    return det;
197
198
199
200
201
202
203
204
  }


  /****************************************************************************/
  /*  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                                 */
  /****************************************************************************/
Thomas Witkowski's avatar
Thomas Witkowski committed
205
  double ElInfo1d::getElementNormal(WorldVector<double> &elementNormal) const
206
  {
207
    FUNCNAME_DBG("ElInfo::getElementNormal()");
208

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

Thomas Witkowski's avatar
Thomas Witkowski committed
212
213
    elementNormal[0] = coord[1][1] - coord[0][1];
    elementNormal[1] = coord[0][0] - coord[1][0];
214

Thomas Witkowski's avatar
Thomas Witkowski committed
215
    double det = norm(&elementNormal);
216

217
    TEST_EXIT_DBG(det > 1.e-30)("det = 0");
218
219
220

    elementNormal *= 1.0 / det;
    
221
    return det;
222
223
224
  }


225
  void ElInfo1d::fillElInfo(int ichild, const ElInfo *elInfoOld)
226
  {
227
    FUNCNAME_DBG("ElInfo1d::fillElInfo()");
228

229
    Element *nb;
Thomas Witkowski's avatar
Thomas Witkowski committed
230
    Element *elem = elInfoOld->element;
231

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
237
238
239
    macroElement = elInfoOld->macroElement;
    fillFlag = elInfoOld->fillFlag;
    parent = elem;
240
    level = elInfoOld->level + 1;
241
    iChild = ichild;
242

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

Thomas Witkowski's avatar
Thomas Witkowski committed
245
246
    if (fillFlag.isSet(Mesh::FILL_COORDS) || fillFlag.isSet(Mesh::FILL_DET) ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
247

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

Thomas Witkowski's avatar
Thomas Witkowski committed
250
      coord[ichild] = (*old_coord)[ichild];
251
      if (elem->isNewCoordSet())
Thomas Witkowski's avatar
Thomas Witkowski committed
252
	coord[1 - ichild] = *(elem->getNewCoord());
253
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
254
	coord[1 - ichild] = ((*old_coord)[0] + (*old_coord)[1]) * 0.5;
255
    }
256

Thomas Witkowski's avatar
Thomas Witkowski committed
257
    if (fillFlag.isSet(Mesh::FILL_NEIGH) || fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {
258
259
      WorldVector<double> oppC;
      
Thomas Witkowski's avatar
Thomas Witkowski committed
260
      TEST_EXIT_DBG(fillFlag.isSet(Mesh::FILL_COORDS))
261
262
263
264
	("FILL_OPP_COORDS only with FILL_COORDS\n");
      
      for (int i = 0; i < neighbours; i++) {
	if (i != ichild) {
265
	  nb = const_cast<Element*>(elem->getChild(1-ichild));  
Thomas Witkowski's avatar
Thomas Witkowski committed
266
267
	  if (fillFlag.isSet(Mesh::FILL_OPP_COORDS))
	    oppC = elInfoOld->coord[i];
268
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
269
	  nb = const_cast<Element*>(elInfoOld->getNeighbour(i));
270

Thomas Witkowski's avatar
Thomas Witkowski committed
271
272
	  if (nb && fillFlag.isSet(Mesh::FILL_OPP_COORDS))
	    oppC = elInfoOld->oppCoord[i];
273
	}
274

275
276
	if (nb) {
	  while (nb->getChild(0)) {  // make nb nearest element
Thomas Witkowski's avatar
Thomas Witkowski committed
277
	    if (fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {
278
	      if (nb->isNewCoordSet())
279
		oppC = *(nb->getNewCoord());
280
	      else
Thomas Witkowski's avatar
Thomas Witkowski committed
281
		oppC = (coord[i] + oppC) * 0.5;
282
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
283
	    nb = const_cast<Element*>(nb->getChild(1-i));
284
	  }
285

Thomas Witkowski's avatar
Thomas Witkowski committed
286
287
	  if (fillFlag.isSet(Mesh::FILL_OPP_COORDS))
	    oppCoord[i] = oppC;
288
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
289
	neighbour[i] = nb;
290
	oppVertex[i] = nb ? i : -1;
291
      }
292
    }
293

Thomas Witkowski's avatar
Thomas Witkowski committed
294
295
296
    if (fillFlag.isSet(Mesh::FILL_BOUND)) {
      boundary[ichild] = elInfoOld->getBoundary(ichild);
      boundary[1 - ichild] = INTERIOR;
297

298
      if (elInfoOld->getProjection(0) && 
Thomas Witkowski's avatar
Thomas Witkowski committed
299
300
	  elInfoOld->getProjection(0)->getType() == VOLUME_PROJECTION)
	projection[0] = elInfoOld->getProjection(0);      
301
    }   
302
303
  }

304

305
  mtl::dense2D<double>& ElInfo1d::getSubElemCoordsMat(int degree) const
306
  {
307
308
    FUNCNAME("ElInfo1d::getSubElemCoordsMat()");

309
310
    using namespace mtl;

Thomas Witkowski's avatar
Thomas Witkowski committed
311
    if (subElemMatrices[degree].count(std::make_pair(refinementPathLength, refinementPath)) == 0) {
312
313
314
315
316
317
318
319
      switch (degree) {
      case 1:
	{
	  dense2D<double> mat(mat_d1);
	  dense2D<double> tmpMat(num_rows(mat), num_rows(mat));
	  
	  for (int i = 0; i < refinementPathLength; i++) {
	    if (refinementPath & (1 << i)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
320
	      tmpMat = mat * mat_d1_right;
321
322
	      mat = tmpMat;
	    } else  {
Thomas Witkowski's avatar
Thomas Witkowski committed
323
	      tmpMat = mat * mat_d1_left;
324
325
326
	      mat = tmpMat;
	    }
	  }
327

Thomas Witkowski's avatar
Thomas Witkowski committed
328
	  subElemMatrices[1][std::make_pair(refinementPathLength, refinementPath)] = mat;
329
330
331
332
333
334
335
336
337
338
	}
	
	break;	
      case 2:
	{
	  dense2D<double> mat(mat_d2);
	  dense2D<double> tmpMat(num_rows(mat), num_rows(mat));
	  
	  for (int i = 0; i < refinementPathLength; i++) {
	    if (refinementPath & (1 << i)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
339
	      tmpMat = mat * mat_d2_right;
340
341
	      mat = tmpMat;
	    } else  {
Thomas Witkowski's avatar
Thomas Witkowski committed
342
	      tmpMat = mat * mat_d2_left;
343
344
345
	      mat = tmpMat;
	    }
	  }
346

Thomas Witkowski's avatar
Thomas Witkowski committed
347
	  subElemMatrices[2][std::make_pair(refinementPathLength, refinementPath)] = mat;  
348
349
350
351
	}
	break;
      default:
	ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
352
      }
353
    }
354

Thomas Witkowski's avatar
Thomas Witkowski committed
355
    return subElemMatrices[degree][std::make_pair(refinementPathLength, refinementPath)];
356
357
  }

358

359
  mtl::dense2D<double>& ElInfo1d::getSubElemGradCoordsMat(int degree) const
360
  {
361
    FUNCNAME("ElInfo1d::getSubElemGradCoordsMat()");
362

363
    TEST_EXIT(degree == 1)("Not supported for basis functions with degree > 1!\n");
364

365
    using namespace mtl;
366

Thomas Witkowski's avatar
Thomas Witkowski committed
367
    if (subElemGradMatrices[degree].count(std::make_pair(refinementPathLength, refinementPath)) == 0) {
368
      dense2D<double> mat(mat_d1);
369

370
371
      for (int i = 0; i < refinementPathLength; i++)
	mat *= 0.5;
372

Thomas Witkowski's avatar
Thomas Witkowski committed
373
      subElemGradMatrices[1][std::make_pair(refinementPathLength, refinementPath)] = mat;
374
    }
375

Thomas Witkowski's avatar
Thomas Witkowski committed
376
    return subElemGradMatrices[degree][std::make_pair(refinementPathLength, refinementPath)];
377
378
  }

379

380
}