ElInfo1d.cc 10.3 KB
Newer Older
1 2 3 4 5 6 7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9 10 11 12 13 14 15 16 17
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
18
 *
19
 ******************************************************************************/
20 21


22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
#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 {
37 38

  double ElInfo1d::mat_d1_val[2][2] = {{1.0, 0.0},
39 40 41
				       {0.0, 1.0}};
  mtl::dense2D<double> ElInfo1d::mat_d1(mat_d1_val);

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

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

50 51 52 53 54 55 56 57 58
  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);
59

60 61 62 63
  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);
64

65

66 67 68 69 70
  void ElInfo1d::fillMacroInfo(const MacroElement * mel)
  {
    Element *nb;
    MacroElement *mnb;

Thomas Witkowski's avatar
Thomas Witkowski committed
71 72
    macroElement = const_cast<MacroElement*>(mel);
    element = const_cast<Element*>(mel->getElement());
73
    parent = NULL;
74
    level = 0;
75

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

Thomas Witkowski's avatar
Thomas Witkowski committed
78 79
    if (fillFlag.isSet(Mesh::FILL_COORDS) || fillFlag.isSet(Mesh::FILL_DET) ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
80

81
      for (int i = 0; i < vertices; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
82
	coord[i] = mel->coord[i];
83
    }
84

Thomas Witkowski's avatar
Thomas Witkowski committed
85
    if (fillFlag.isSet(Mesh::FILL_NEIGH) || fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {
86
      WorldVector<double> oppC;
87

Thomas Witkowski's avatar
Thomas Witkowski committed
88
      int neighbours =  mesh->getGeo(NEIGH);
89
      for (int i = 0; i < neighbours; i++) {
90
	nb = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
91
	if ((mnb = const_cast<MacroElement*>(mel->getNeighbour(i)))) {
Thomas Witkowski's avatar
Thomas Witkowski committed
92
	  if (fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {
93 94 95
	    oppC = mnb->coord[i];
	  }

Thomas Witkowski's avatar
Thomas Witkowski committed
96
	  nb = const_cast<Element*>(mnb->getElement());
97 98

	  while (!(nb->isLeaf())) { // make nb nearest element
Thomas Witkowski's avatar
Thomas Witkowski committed
99
	    if (fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {
100
	      if (nb->isNewCoordSet()) {
101 102 103
		oppC = *(nb->getNewCoord());
	      } else {
		oppC = (mel->coord[i] + oppC) * 0.5;
104
	      }
105
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
106
	    nb = const_cast<Element*>(nb->getChild(1-i));
107 108
	  }

Thomas Witkowski's avatar
Thomas Witkowski committed
109 110
	  if (fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {
	    oppCoord[i] = oppC;
111
	  }
112
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
113
	neighbour[i] = nb;
114
	oppVertex[i] = nb ? i : -1;
115
      }
116
    }
117

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

Thomas Witkowski's avatar
Thomas Witkowski committed
122 123
      for (int i = 0; i < element->getGeo(PROJECTION); i++)
	projection[i] = mel->getProjection(i);
124 125 126 127 128 129 130 131
    }
  }

  /****************************************************************************/
  /*  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
132
  double ElInfo1d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
133
  {
134
    FUNCNAME("ElInfo1d::calcGrdLambda()");
135 136 137 138

    testFlag(Mesh::FILL_COORDS);

    WorldVector<double> e;
139
    e = coord[1];
Thomas Witkowski's avatar
Thomas Witkowski committed
140
    e -= coord[0];
141
    double adet2 = e * e;
142 143

    if (adet2 < 1.0E-15) {
144
      MSG("det*det = %g\n", adet2);
145
      grd_lam[0] = grd_lam[1] = 0.0;
146
    } else {
147
      grd_lam[1] = e * (1.0 / adet2);
148
      grd_lam[0] = grd_lam[1] * -1.0;
149 150 151 152 153
    }

    return sqrt(adet2);
  }

154
  int ElInfo1d::worldToCoord(const WorldVector<double>& x,
155 156
				   DimVec<double>* lambda) const
  {
157 158
    FUNCNAME("ElInfo1d::worldToCoord()");

159
    double lmin;
Thomas Witkowski's avatar
Thomas Witkowski committed
160 161 162
    double a = coord[0][0];
    double length = (coord[1][0] - a);
    int dim = mesh->getDim();
163

164
    TEST_EXIT_DBG(lambda)("lambda must not be NULL\n");
165
    TEST_EXIT_DBG(dim == 1)("dim!=1\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
166
    TEST_EXIT_DBG(dimOfWorld == dim)("not yet for DIM != DIM_OF_WORLD\n");
167 168 169 170 171 172

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

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

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

203
    return det;
204 205 206 207 208 209 210 211
  }


  /****************************************************************************/
  /*  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
212
  double ElInfo1d::getElementNormal(WorldVector<double> &elementNormal) const
213
  {
214
    FUNCNAME_DBG("ElInfo::getElementNormal()");
215

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

Thomas Witkowski's avatar
Thomas Witkowski committed
219 220
    elementNormal[0] = coord[1][1] - coord[0][1];
    elementNormal[1] = coord[0][0] - coord[1][0];
221

Thomas Witkowski's avatar
Thomas Witkowski committed
222
    double det = norm(&elementNormal);
223

224
    TEST_EXIT_DBG(det > 1.e-30)("det = 0");
225 226

    elementNormal *= 1.0 / det;
227

228
    return det;
229 230 231
  }


232
  void ElInfo1d::fillElInfo(int ichild, const ElInfo *elInfoOld)
233
  {
234
    FUNCNAME_DBG("ElInfo1d::fillElInfo()");
235

236
    Element *nb;
Thomas Witkowski's avatar
Thomas Witkowski committed
237
    Element *elem = elInfoOld->element;
238

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
244 245 246
    macroElement = elInfoOld->macroElement;
    fillFlag = elInfoOld->fillFlag;
    parent = elem;
247
    level = elInfoOld->level + 1;
248
    iChild = ichild;
249

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

Thomas Witkowski's avatar
Thomas Witkowski committed
252 253
    if (fillFlag.isSet(Mesh::FILL_COORDS) || fillFlag.isSet(Mesh::FILL_DET) ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
254

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

Thomas Witkowski's avatar
Thomas Witkowski committed
257
      coord[ichild] = (*old_coord)[ichild];
258
      if (elem->isNewCoordSet())
Thomas Witkowski's avatar
Thomas Witkowski committed
259
	coord[1 - ichild] = *(elem->getNewCoord());
260
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
261
	coord[1 - ichild] = ((*old_coord)[0] + (*old_coord)[1]) * 0.5;
262
    }
263

Thomas Witkowski's avatar
Thomas Witkowski committed
264
    if (fillFlag.isSet(Mesh::FILL_NEIGH) || fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {
265
      WorldVector<double> oppC;
266

Thomas Witkowski's avatar
Thomas Witkowski committed
267
      TEST_EXIT_DBG(fillFlag.isSet(Mesh::FILL_COORDS))
268
	("FILL_OPP_COORDS only with FILL_COORDS\n");
269

270 271
      for (int i = 0; i < neighbours; i++) {
	if (i != ichild) {
272
	  nb = const_cast<Element*>(elem->getChild(1-ichild));
Thomas Witkowski's avatar
Thomas Witkowski committed
273 274
	  if (fillFlag.isSet(Mesh::FILL_OPP_COORDS))
	    oppC = elInfoOld->coord[i];
275
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
276
	  nb = const_cast<Element*>(elInfoOld->getNeighbour(i));
277

Thomas Witkowski's avatar
Thomas Witkowski committed
278 279
	  if (nb && fillFlag.isSet(Mesh::FILL_OPP_COORDS))
	    oppC = elInfoOld->oppCoord[i];
280
	}
281

282 283
	if (nb) {
	  while (nb->getChild(0)) {  // make nb nearest element
Thomas Witkowski's avatar
Thomas Witkowski committed
284
	    if (fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {
285
	      if (nb->isNewCoordSet())
286
		oppC = *(nb->getNewCoord());
287
	      else
Thomas Witkowski's avatar
Thomas Witkowski committed
288
		oppC = (coord[i] + oppC) * 0.5;
289
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
290
	    nb = const_cast<Element*>(nb->getChild(1-i));
291
	  }
292

Thomas Witkowski's avatar
Thomas Witkowski committed
293 294
	  if (fillFlag.isSet(Mesh::FILL_OPP_COORDS))
	    oppCoord[i] = oppC;
295
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
296
	neighbour[i] = nb;
297
	oppVertex[i] = nb ? i : -1;
298
      }
299
    }
300

Thomas Witkowski's avatar
Thomas Witkowski committed
301 302 303
    if (fillFlag.isSet(Mesh::FILL_BOUND)) {
      boundary[ichild] = elInfoOld->getBoundary(ichild);
      boundary[1 - ichild] = INTERIOR;
304

305
      if (elInfoOld->getProjection(0) &&
Thomas Witkowski's avatar
Thomas Witkowski committed
306
	  elInfoOld->getProjection(0)->getType() == VOLUME_PROJECTION)
307 308
	projection[0] = elInfoOld->getProjection(0);
    }
309 310
  }

311

312
  mtl::dense2D<double>& ElInfo1d::getSubElemCoordsMat(int degree) const
313
  {
314 315
    FUNCNAME("ElInfo1d::getSubElemCoordsMat()");

316 317
    using namespace mtl;

Thomas Witkowski's avatar
Thomas Witkowski committed
318
    if (subElemMatrices[degree].count(std::make_pair(refinementPathLength, refinementPath)) == 0) {
319 320 321 322 323
      switch (degree) {
      case 1:
	{
	  dense2D<double> mat(mat_d1);
	  dense2D<double> tmpMat(num_rows(mat), num_rows(mat));
324

325 326
	  for (int i = 0; i < refinementPathLength; i++) {
	    if (refinementPath & (1 << i)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
327
	      tmpMat = mat * mat_d1_right;
328 329
	      mat = tmpMat;
	    } else  {
Thomas Witkowski's avatar
Thomas Witkowski committed
330
	      tmpMat = mat * mat_d1_left;
331 332 333
	      mat = tmpMat;
	    }
	  }
334

Thomas Witkowski's avatar
Thomas Witkowski committed
335
	  subElemMatrices[1][std::make_pair(refinementPathLength, refinementPath)] = mat;
336
	}
337 338

	break;
339 340 341 342
      case 2:
	{
	  dense2D<double> mat(mat_d2);
	  dense2D<double> tmpMat(num_rows(mat), num_rows(mat));
343

344 345
	  for (int i = 0; i < refinementPathLength; i++) {
	    if (refinementPath & (1 << i)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
346
	      tmpMat = mat * mat_d2_right;
347 348
	      mat = tmpMat;
	    } else  {
Thomas Witkowski's avatar
Thomas Witkowski committed
349
	      tmpMat = mat * mat_d2_left;
350 351 352
	      mat = tmpMat;
	    }
	  }
353

354
	  subElemMatrices[2][std::make_pair(refinementPathLength, refinementPath)] = mat;
355 356 357 358
	}
	break;
      default:
	ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
359
      }
360
    }
361

Thomas Witkowski's avatar
Thomas Witkowski committed
362
    return subElemMatrices[degree][std::make_pair(refinementPathLength, refinementPath)];
363 364
  }

365

366
  mtl::dense2D<double>& ElInfo1d::getSubElemGradCoordsMat(int degree) const
367
  {
368
    return getSubElemCoordsMat(degree);
369 370
  }

371

372
}