Am Montag, 13. Mai 2022, finden Wartungsarbeiten am Gitlab-Server (Update auf neue Version statt). Der Dienst wird daher am Montag für einige Zeit nicht verfügbar sein.
On Monday, May 13th 2022, the Gitlab server will be updated. The service will therefore not be accessible for some time on Monday.

ElInfo1d.cc 8.55 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

  void ElInfo1d::fillMacroInfo(const MacroElement * mel)
  {
31
    FUNCNAME("ElInfo1d::fillMacroInfo()");
32
33
34
    Element *nb;
    MacroElement *mnb;

35
36
37
    macroElement_ = const_cast<MacroElement*>( mel);
    element_ = const_cast<Element*>( mel->getElement());
    parent_ = NULL;
38
    level = 0;
39
40
41
42

    int vertices = mesh_->getGeo(VERTEX);

    if (fillFlag_.isSet(Mesh::FILL_COORDS) || fillFlag_.isSet(Mesh::FILL_DET) ||
43
44
	fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
      
45
      for (int i = 0; i < vertices; i++)
46
	coord_[i] = mel->coord[i];
47
    }
48

49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
    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)) {
64
	      if (nb->isNewCoordSet()) {
65
66
67
		oppC = *(nb->getNewCoord());
	      } else {
		oppC = (mel->coord[i] + oppC) * 0.5;
68
	      }
69
70
71
72
73
74
	    }
	    nb = const_cast<Element*>( nb->getChild(1-i));
	  }

	  if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
	    oppCoord_[i] = oppC;
75
	  }
76
77
78
	}
	neighbour_[i] = nb;
	oppVertex_[i] = nb ? i : -1;
79
      }
80
    }
81
82

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

86
      for (int i = 0; i < element_->getGeo(PROJECTION); i++)
87
88
89
90
91
92
93
94
95
	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                                                                 */
  /****************************************************************************/
Thomas Witkowski's avatar
Thomas Witkowski committed
96
  double ElInfo1d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
97
  {
98
    FUNCNAME("ElInfo1d::calcGrdLambda()");
99
100
101
102

    testFlag(Mesh::FILL_COORDS);

    WorldVector<double> e;
103
104
105
    e = coord_[1]; 
    e -= coord_[0];
    double adet2 = e * e;
106
107
108

    if (adet2 < 1.0E-15) {
      MSG("det*det = %lf\n", adet2);
109
      for (int i = 0; i <= 1; i++)
110
111
	grd_lam[i] = 0.0;
    } else {
112
      grd_lam[1] = e * (1.0 / adet2);
113
114
115
116
117
118
119
120
121
      grd_lam[0] = grd_lam[1] * (-1.0);
    }

    return sqrt(adet2);
  }

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

124
125
126
127
128
129
130
    double lmin;
    double a = coord_[0][0];
    double length = (coord_[1][0] - a);
    int dim = mesh_->getDim();

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

131
132
    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
133
    TEST_EXIT_DBG(dimOfWorld == dim)("not yet for DIM != DIM_OF_WORLD\n");
134
135
136
137
138
139

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

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

143
    int k = -1;
144
    lmin = 0.0;
145
    for (int i = 0; i <= dim; i++) {
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
      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
162
  double ElInfo1d::getNormal(int side, WorldVector<double> &normal)
163
164
  {
    normal = coord_[side] - coord_[(side + 1) % 2];
165
    double det = norm(&normal);
166
    TEST_EXIT_DBG(det > 1.e-30)("det = 0 on side %d\n", side);
167
168
    normal *= 1.0 / det;

169
    return det;
170
171
172
173
174
175
176
177
178
179
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
182
183
    TEST_EXIT_DBG(dimOfWorld == 2)
      (" element normal only well defined for  DIM_OF_WORLD = DIM + 1 !!");
184
185
186
187

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

Thomas Witkowski's avatar
Thomas Witkowski committed
188
    double det = norm(&elementNormal);
189

190
    TEST_EXIT_DBG(det > 1.e-30)("det = 0");
191
192
193

    elementNormal *= 1.0 / det;
    
194
    return det;
195
196
197
  }


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

202
    Element *nb;
203
    Element *elem = elInfoOld->element_;
204

205
206
207
208
    TEST_EXIT_DBG(elem->getChild(0))("no children?\n");
    element_ = const_cast<Element*>(elem->getChild(ichild));

    TEST_EXIT_DBG(element_)("missing child %d?\n", ichild);
209

210
211
    macroElement_ = elInfoOld->macroElement_;
    fillFlag_ = elInfoOld->fillFlag_;
212
    parent_ = elem;
213
    level = elInfoOld->level + 1;
214
    iChild = ichild;
215
216
217
218

    int neighbours = mesh_->getGeo(NEIGH);

    if (fillFlag_.isSet(Mesh::FILL_COORDS) || fillFlag_.isSet(Mesh::FILL_DET) ||
219
	fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
220

221
      const FixVec<WorldVector<double>, VERTEX> *old_coord = &(elInfoOld->coord_);
222
223

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

231
232
233
    if (fillFlag_.isSet(Mesh::FILL_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
      WorldVector<double> oppC;
      
234
      TEST_EXIT_DBG(fillFlag_.isSet(Mesh::FILL_COORDS))
235
236
237
238
239
240
	("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)) {
241
	    oppC = elInfoOld->coord_[i];
242
243
	  }
	} else {
244
	  nb = const_cast<Element*>( elInfoOld->getNeighbour(i));
245

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

251
252
253
	if (nb) {
	  while (nb->getChild(0)) {  // make nb nearest element
	    if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
254
	      if (nb->isNewCoordSet()) {
255
256
257
258
259
260
261
		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
    if (fillFlag_.isSet(Mesh::FILL_BOUND)) {
273
274
      boundary_[ichild] = elInfoOld->getBoundary(ichild);
      boundary_[1 - ichild] = INTERIOR;
275

276
277
278
      if (elInfoOld->getProjection(0) && 
	  elInfoOld->getProjection(0)->getType() == VOLUME_PROJECTION) {
	projection_[0] = elInfoOld->getProjection(0);
279
      }
280
    }   
281
282
  }

283
284
285
286
287
  void ElInfo1d::getRefSimplexCoords(const BasisFunction *basisFcts,
				     mtl::dense2D<double>& coords) const
  {
    TEST_EXIT(basisFcts->getDegree() == 1)("Wrong basis function degree!\n");

288
    coords = mat_d1;
289
290
  }

291
292
293
294
295
296
297
298
299
  void ElInfo1d::getSubElementCoords(const BasisFunction *basisFcts,
				     int iChild,
				     mtl::dense2D<double>& coords) const
  {
    FUNCNAME("ElInfo1d::getSubElementCoords()");

    switch (basisFcts->getDegree()) {
    case 1:
      if (iChild == 0)
300
	coords *= mat_d1_left;
301
      else
302
	coords *= mat_d1_right;
303
304
305
306
307
308
309
310
311

      break;
    case 2:
      break;
    default:
      ERROR_EXIT("Not yet implemented!\n");
    }
  }

312
}