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.58 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
    macroElement_ = const_cast<MacroElement*>( mel);
    element_ = const_cast<Element*>( mel->getElement());
    parent_ = NULL;
26
    level = 0;
27
28
29
30

    int vertices = mesh_->getGeo(VERTEX);

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

37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
    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)) {
52
	      if (nb->isNewCoordSet()) {
53
54
55
		oppC = *(nb->getNewCoord());
	      } else {
		oppC = (mel->coord[i] + oppC) * 0.5;
56
	      }
57
58
59
60
61
62
	    }
	    nb = const_cast<Element*>( nb->getChild(1-i));
	  }

	  if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
	    oppCoord_[i] = oppC;
63
	  }
64
65
66
	}
	neighbour_[i] = nb;
	oppVertex_[i] = nb ? i : -1;
67
      }
68
    }
69
70

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

74
      for (int i = 0; i < element_->getGeo(PROJECTION); i++)
75
76
77
78
79
80
81
82
83
	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
84
  double ElInfo1d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
85
  {
86
    FUNCNAME("ElInfo1d::calcGrdLambda()");
87
88
89
90

    testFlag(Mesh::FILL_COORDS);

    WorldVector<double> e;
91
92
93
    e = coord_[1]; 
    e -= coord_[0];
    double adet2 = e * e;
94
95
96

    if (adet2 < 1.0E-15) {
      MSG("det*det = %lf\n", adet2);
97
      for (int i = 0; i <= 1; i++)
98
99
	grd_lam[i] = 0.0;
    } else {
100
      grd_lam[1] = e * (1.0 / adet2);
101
102
103
104
105
106
107
108
109
      grd_lam[0] = grd_lam[1] * (-1.0);
    }

    return sqrt(adet2);
  }

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

112
113
114
115
116
117
118
    double lmin;
    double a = coord_[0][0];
    double length = (coord_[1][0] - a);
    int dim = mesh_->getDim();

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

119
120
    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
121
    TEST_EXIT_DBG(dimOfWorld == dim)("not yet for DIM != DIM_OF_WORLD\n");
122
123
124
125
126
127

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

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

131
    int k = -1;
132
    lmin = 0.0;
133
    for (int i = 0; i <= dim; i++) {
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
      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
150
  double ElInfo1d::getNormal(int side, WorldVector<double> &normal)
151
152
  {
    normal = coord_[side] - coord_[(side + 1) % 2];
153
    double det = norm(&normal);
154
    TEST_EXIT_DBG(det > 1.e-30)("det = 0 on side %d\n", side);
155
156
    normal *= 1.0 / det;

157
    return det;
158
159
160
161
162
163
164
165
166
167
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
170
171
    TEST_EXIT_DBG(dimOfWorld == 2)
      (" element normal only well defined for  DIM_OF_WORLD = DIM + 1 !!");
172
173
174
175

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

Thomas Witkowski's avatar
Thomas Witkowski committed
176
    double det = norm(&elementNormal);
177

178
    TEST_EXIT_DBG(det > 1.e-30)("det = 0");
179
180
181

    elementNormal *= 1.0 / det;
    
182
    return det;
183
184
185
  }


186
  void ElInfo1d::fillElInfo(int ichild, const ElInfo *elInfoOld)
187
  {
188
189
    FUNCNAME("ElInfo1d::fillElInfo()");

190
    Element *nb;
191
    Element *elem = elInfoOld->element_;
192

193
194
195
196
    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);
197

198
199
    macroElement_ = elInfoOld->macroElement_;
    fillFlag_ = elInfoOld->fillFlag_;
200
    parent_ = elem;
201
    level = elInfoOld->level + 1;
202
    iChild = ichild;
203
204
205
206

    int neighbours = mesh_->getGeo(NEIGH);

    if (fillFlag_.isSet(Mesh::FILL_COORDS) || fillFlag_.isSet(Mesh::FILL_DET) ||
207
	fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
208

209
      const FixVec<WorldVector<double>, VERTEX> *old_coord = &(elInfoOld->coord_);
210
211

      coord_[ichild] = (*old_coord)[ichild];
212
      if (elem->isNewCoordSet()) {
213
214
215
	coord_[1 - ichild] = *(elem->getNewCoord());
      } else {
	coord_[1 - ichild] = ((*old_coord)[0] + (*old_coord)[1]) * 0.5;
216
      }
217
    }
218

219
220
221
    if (fillFlag_.isSet(Mesh::FILL_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
      WorldVector<double> oppC;
      
222
      TEST_EXIT_DBG(fillFlag_.isSet(Mesh::FILL_COORDS))
223
224
225
226
227
228
	("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)) {
229
	    oppC = elInfoOld->coord_[i];
230
231
	  }
	} else {
232
	  nb = const_cast<Element*>( elInfoOld->getNeighbour(i));
233

234
	  if (nb && fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
235
	    oppC = elInfoOld->oppCoord_[i];
236
	  }
237
	}
238

239
240
241
	if (nb) {
	  while (nb->getChild(0)) {  // make nb nearest element
	    if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
242
	      if (nb->isNewCoordSet()) {
243
244
245
246
247
248
249
		oppC = *(nb->getNewCoord());
	      } else {
		oppC = (coord_[i] + oppC) * 0.5;
	      }
	    }
	    nb = const_cast<Element*>( nb->getChild(1-i));
	  }
250

251
252
	  if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
	    oppCoord_[i] = oppC;
253
	  }
254
255
256
	}
	neighbour_[i] = nb;
	oppVertex_[i] = nb ? i : -1;
257
      }
258
    }
259

260
    if (fillFlag_.isSet(Mesh::FILL_BOUND)) {
261
262
      boundary_[ichild] = elInfoOld->getBoundary(ichild);
      boundary_[1 - ichild] = INTERIOR;
263

264
265
266
      if (elInfoOld->getProjection(0) && 
	  elInfoOld->getProjection(0)->getType() == VOLUME_PROJECTION) {
	projection_[0] = elInfoOld->getProjection(0);
267
      }
268
    }   
269
270
  }

271
272
  void ElInfo1d::getRefSimplexCoords(const BasisFunction *basisFcts,
				     DimMat<double> *coords) const
273
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
274
    FUNCNAME("ElInfo1d::getRefSimplexCoords()");
275

Thomas Witkowski's avatar
Thomas Witkowski committed
276
277
278
279
280
281
282
283
284
285
286
287
288
289
    switch (basisFcts->getDegree()) {
    case 1:
      (*coords)[0][0] = 1.0;
      (*coords)[1][0] = 0.0;
      
      (*coords)[0][1] = 0.0;
      (*coords)[1][1] = 1.0;
      break;
    case 2:
      ERROR_EXIT("DDA\n");
      break;
    default:
      ERROR_EXIT("Not yet implemented!\n");
    }
290
291
  }

292
293
294
  void ElInfo1d::getSubElementCoords(const BasisFunction *basisFcts,
				     int iChild,
				     DimMat<double> *coords) const
295
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
    FUNCNAME("ElInfo1d::getSubElementCoords()");

    switch (basisFcts->getDegree()) {
    case 1:
      if (iChild == 0) {
	(*coords)[0][1] = ((*coords)[0][0] + (*coords)[0][1]) * 0.5;
	(*coords)[1][1] = ((*coords)[1][0] + (*coords)[1][1]) * 0.5;
      } else {
	(*coords)[0][0] = ((*coords)[0][0] + (*coords)[0][1]) * 0.5;
	(*coords)[1][0] = ((*coords)[1][0] + (*coords)[1][1]) * 0.5;
      }
      break;
    case 2:
      break;
    default:
      ERROR_EXIT("Not yet implemented!\n");
312
313
314
    }
  }

315
}