ElInfo2d.cc 16.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#include "ElInfo2d.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 {

17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
  ElInfo2d::ElInfo2d(Mesh *aMesh) 
    : ElInfo(aMesh) 
  {
    e1 = NEW WorldVector<double>;
    e2 = NEW WorldVector<double>;
    normal = NEW WorldVector<double>;
  }

  ElInfo2d::~ElInfo2d()
  {
    DELETE e1;
    DELETE e2;
    DELETE normal;
  }

32
33
  void ElInfo2d::fillMacroInfo(const MacroElement * mel)
  {
34
35
    FUNCNAME("ElInfo::fillMacroInfo()");
 
36
    macroElement_ = const_cast<MacroElement*>(mel);
37
38
39
    element_  = const_cast<Element*>(mel->getElement());
    parent_  = NULL;
    level_ = 0;
40
41
42
43

    if (fillFlag_.isSet(Mesh::FILL_COORDS) || 
	fillFlag_.isSet(Mesh::FILL_DET)    ||
	fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
44
45
46

      int vertices = mesh_->getGeo(VERTEX);
      for (int i = 0; i < vertices; i++) {
47
48
49
50
51
52
	coord_[i] = mel->coord[i];
      }
    }

    int neighbours = mesh_->getGeo(NEIGH);

53
54
55
56
    if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS) || 
	fillFlag_.isSet(Mesh::FILL_NEIGH)) {

      bool fill_opp_coords = (fillFlag_.isSet(Mesh::FILL_OPP_COORDS));
57
    
58
59
60
61
      for (int i = 0; i < neighbours; i++) {
	MacroElement *macroNeighbour = mel->getNeighbour(i);

	if (macroNeighbour) {
Thomas Witkowski's avatar
Thomas Witkowski committed
62
	  neighbour_[i] = macroNeighbour->getElement();	  
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
	  Element *nb = const_cast<Element*>(neighbour_[i]);

	  int edgeNo = oppVertex_[i] = mel->getOppVertex(i);
	  if (nb->getFirstChild() && (edgeNo != 2)) {   // make nb nearest el.
	    if (edgeNo == 0) {
	      nb = neighbour_[i] = nb->getSecondChild();
	    } else {
	      nb = neighbour_[i] = nb->getFirstChild();
	    }

	    oppVertex_[i] = 2;

	    if (fill_opp_coords) {
	      if (nb->getNewCoord(-1)) {
		oppCoord_[i] = *(nb->getNewCoord());
	      } else {
		oppCoord_[i] = (macroNeighbour->coord[0] + macroNeighbour->coord[1]) * 0.5;
	      }	    	      
	      
	      switch (i) {
	      case 0:
		if (*(macroNeighbour->getElement()->getDOF(2)) == *(element_->getDOF(2))) {
		  neighbourCoord_[i][0] = macroNeighbour->coord[2];
		  neighbourCoord_[i][1] = macroNeighbour->coord[0];
		} else if (*(macroNeighbour->getElement()->getDOF(2)) == *(element_->getDOF(1))) {
		  neighbourCoord_[i][0] = macroNeighbour->coord[1];
		  neighbourCoord_[i][1] = macroNeighbour->coord[2];
90
		} else {
91
		  ERROR_EXIT("should not happen!\n");
92
		}
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
	
		neighbourCoord_[i][2] = oppCoord_[i];
		break;
		
	      case 1:
		if (*(macroNeighbour->getElement()->getDOF(2)) == *(element_->getDOF(2))) {
		  neighbourCoord_[i][0] = macroNeighbour->coord[1];
		  neighbourCoord_[i][1] = macroNeighbour->coord[2];
		} else if (*(macroNeighbour->getElement()->getDOF(2)) == *(element_->getDOF(0))) {
		  neighbourCoord_[i][0] = macroNeighbour->coord[2];
		  neighbourCoord_[i][1] = macroNeighbour->coord[0];
		} else {
		  ERROR_EXIT("should not happen!\n");
		}
		
		neighbourCoord_[i][2] = oppCoord_[i];
		break;
		
	      default:
Thomas Witkowski's avatar
Thomas Witkowski committed
112
113
114
115
116
117
118
119
120
121
122
123
124
		std::cout << "------------- Error --------------" << std::endl;
		std::cout << "  Element index = " << element_->getIndex() << "\n\n";
		for (int j = 0; j < neighbours; j++) {
		  if (mel->getNeighbour(j)) {
		    std::cout << "  Neighbour " << j << ": " 
			      << mel->getNeighbour(j)->getElement()->getIndex() 
			      << std::endl;
		  } else {
		    std::cout << "  Neighbour " << j << ": not existing" << std::endl;
		  }
		  std::cout << "  OppVertex " << j << ": " << static_cast<int>(mel->getOppVertex(j)) << std::endl;
		  std::cout << std::endl;
		}
125
126
		ERROR_EXIT("should not happen!\n");
		break;
127
128
129
	      }
	    }
	  } else {
130
131
132
133
134
	    if (fill_opp_coords) {
	      oppCoord_[i] = macroNeighbour->coord[edgeNo];

	      neighbourCoord_[i] = macroNeighbour->coord;	      
	    }
135
	  }
136
137
138
	} else {
	  neighbour_[i] = NULL;
        }
139
      }
140
    }
141
    
142
143
    if (fillFlag_.isSet(Mesh::FILL_BOUND)) {   
      for (int i = 0; i < element_->getGeo(BOUNDARY); i++) {
144
145
146
	boundary_[i] = mel->getBoundary(i);
      }

147
      for (int i = 0; i < element_->getGeo(PROJECTION); i++) {
148
149
150
151
152
153
154
155
156
157
158
159
	projection_[i] = mel->getProjection(i);
      }
    }
  }


  /****************************************************************************/
  /*   fill ElInfo structure for one child of an element   		    */
  /****************************************************************************/

  void ElInfo2d::fillElInfo(int ichild, const ElInfo *elinfo_old)
  {
160
    FUNCNAME("ElInfo::fillElInfo()");
161

162
163
164
165
    Element *elem = elinfo_old->element_;
    Element *nb;

    Flag fill_flag = elinfo_old->fillFlag_;
166

167
168
169
170
171
    TEST_EXIT_DBG(elem->getFirstChild())("no children?\n");
    element_ = const_cast<Element*>((ichild == 0) ? 
				    elem->getFirstChild() : 
				    elem->getSecondChild());
    TEST_EXIT_DBG(element_)("missing child %d?\n", ichild);
172
173
174

    macroElement_  = elinfo_old->macroElement_;
    fillFlag_ = fill_flag;
175
176
    parent_ = elem;
    level_ = elinfo_old->level_ + 1;
177
178
179

    if (fillFlag_.isSet(Mesh::FILL_COORDS) || 
	fillFlag_.isSet(Mesh::FILL_DET)    ||
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
	fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
      
      if (elem->getNewCoord(-1)) {
	coord_[2] = *(elem->getNewCoord());
      } else {
	coord_[2].setMidpoint(elinfo_old->coord_[0], elinfo_old->coord_[1]);
      }
      
      if (ichild == 0) {
	coord_[0] = elinfo_old->coord_[2];
	coord_[1] = elinfo_old->coord_[0];
      } else {
	coord_[0] = elinfo_old->coord_[1];
	coord_[1] = elinfo_old->coord_[2];
      }
    }

    bool fill_opp_coords = (fill_flag.isSet(Mesh::FILL_OPP_COORDS));

    if (fill_flag.isSet(Mesh::FILL_NEIGH) || fill_opp_coords) {     
      if (ichild == 0) {
	// Calculation of the neighbour 2, its oppCoords and the
	// cooresponding oppVertex.

	neighbour_[2] = elinfo_old->neighbour_[1];
	oppVertex_[2] = elinfo_old->oppVertex_[1];
	
	if (neighbour_[2] && fill_opp_coords) {
	  oppCoord_[2] = elinfo_old->oppCoord_[1];
	  neighbourCoord_[2] = elinfo_old->neighbourCoord_[1];
210
	}
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
	
	
	// Calculation of the neighbour 1, its oppCoords and the
	// cooresponding oppVertex.
	
	if (elem->getFirstChild()  &&  
	    elem->getSecondChild()->getFirstChild()  &&  
	    elem->getSecondChild()->getFirstChild()) {
	  
	  neighbour_[1] = elem->getSecondChild()->getSecondChild();
	  oppVertex_[1] = 2;
	  
	  if (fill_opp_coords) {
	    if (elem->getSecondChild()->getNewCoord(-1)) {
	      oppCoord_[1] = *(elem->getSecondChild()->getNewCoord());
	    } else {      
	      oppCoord_[1].setMidpoint(elinfo_old->coord_[1], 
				       elinfo_old->coord_[2]);
	    }
230

231
232
233
	    neighbourCoord_[1][0] = coord_[0];
	    neighbourCoord_[1][1] = coord_[2];
	    neighbourCoord_[1][2] = oppCoord_[1];  
234
235
	  }
	} else {
236
237
238
239
240
241
242
243
244
	  neighbour_[1] = elem->getSecondChild();
	  oppVertex_[1] = 0;

	  if (fill_opp_coords) {
	    oppCoord_[1] = elinfo_old->coord_[1];

	    neighbourCoord_[1][0] = elinfo_old->coord_[1];
	    neighbourCoord_[1][1] = elinfo_old->coord_[2];
	    neighbourCoord_[1][2] = coord_[2];
245
246
247
248
	  }
	}


249
250
251
252
253
254
	// Calculation of the neighbour 0, its oppCoords and the
	// cooresponding oppVertex.
	
	nb = elinfo_old->neighbour_[2];
	if (nb) {
	  TEST(elinfo_old->oppVertex_[2] == 2)("invalid neighbour\n"); 
255
256
	  TEST_EXIT_DBG(nb->getFirstChild())("missing first child?\n");
	  TEST_EXIT_DBG(nb->getSecondChild())("missing second child?\n");
257
258
259
260
261
262
	 
	  nb = nb->getSecondChild();

	  if (nb->getFirstChild()) {
	    oppVertex_[0] = 2;

263
	    if (fill_opp_coords) {
264
265
266
267
268
	      if (nb->getNewCoord(-1)) {
		oppCoord_[0] = *(nb->getNewCoord());
	      } else {
		oppCoord_[0].setMidpoint(elinfo_old->neighbourCoord_[2][1],
					 elinfo_old->neighbourCoord_[2][2]);
269
	      }
270
271
272
273
274
275
276
277
278
279
280

	      neighbourCoord_[0][0].setMidpoint(elinfo_old->neighbourCoord_[2][0],
						elinfo_old->neighbourCoord_[2][1]);
	      neighbourCoord_[0][1] = elinfo_old->neighbourCoord_[2][1];
	      neighbourCoord_[0][2] = oppCoord_[0];
	    }	   
 
	    nb = nb->getFirstChild();
	  } else {
	    oppVertex_[0] = 1;

281
	    if (fill_opp_coords) {
282
283
284
285
286
287
	      oppCoord_[0] = elinfo_old->oppCoord_[2];    

	      neighbourCoord_[0][0] = elinfo_old->neighbourCoord_[2][0];
	      neighbourCoord_[0][1] = elinfo_old->neighbourCoord_[2][2];
	      neighbourCoord_[0][2].setMidpoint(elinfo_old->neighbourCoord_[2][0],
						elinfo_old->neighbourCoord_[2][1]);
288
289
	    }
	  }
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
	}
	
	neighbour_[0] = nb;
      } else {   /* ichild == 1 */
	// Calculation of the neighbour 2, its oppCoords and the
	// cooresponding oppVertex.

	neighbour_[2] = elinfo_old->neighbour_[0];
	oppVertex_[2] = elinfo_old->oppVertex_[0];

	if (neighbour_[2] && fill_opp_coords) {
	  oppCoord_[2] = elinfo_old->oppCoord_[0];
	  neighbourCoord_[2] = elinfo_old->neighbourCoord_[0];
	}
	

	// Calculation of the neighbour 0, its oppCoords and the
	// cooresponding oppVertex.

	if (elem->getFirstChild()->getFirstChild()) {
	  neighbour_[0] = elem->getFirstChild()->getFirstChild();
	  oppVertex_[0] = 2;

	  if (fill_opp_coords) {
	    if (elem->getFirstChild()->getNewCoord(-1)) {
	      oppCoord_[0] = *(elem->getFirstChild()->getNewCoord());
316
	    } else {
317
318
	      oppCoord_[0].setMidpoint(elinfo_old->coord_[0], 
				       elinfo_old->coord_[2]);
319
	    }
320
321
322
323

	    neighbourCoord_[0][0] = coord_[2];
	    neighbourCoord_[0][1] = coord_[1];
	    neighbourCoord_[0][2] = oppCoord_[0];
324
	  }
325
326
327
328
329
330
331
332
333
334
	} else {
	  neighbour_[0] = elem->getFirstChild();
	  oppVertex_[0] = 1;

	  if (fill_opp_coords) {
	    oppCoord_[0] = elinfo_old->coord_[0];

	    neighbourCoord_[0][0] = elinfo_old->coord_[2];
	    neighbourCoord_[0][1] = elinfo_old->coord_[0];
	    neighbourCoord_[0][2] = coord_[2];
335
	  }
336
337
338
339
340
341
342
343
344
345
346
347
348
	}

	// Calculation of the neighbour 1, its oppCoords and the
	// cooresponding oppVertex.

	nb = elinfo_old->neighbour_[2];
	if (nb) {
	  TEST(elinfo_old->oppVertex_[2] == 2)("invalid neighbour\n"); 
	  TEST((nb = nb->getFirstChild()))("missing child?\n");

	  if (nb->getFirstChild()) {
	    oppVertex_[1] = 2;

349
	    if (fill_opp_coords) {
350
351
	      if (nb->getNewCoord(-1)) {
		oppCoord_[1] = *(nb->getNewCoord());
352
	      } else {
353
354
		oppCoord_[1].setMidpoint(elinfo_old->neighbourCoord_[2][0],
					 elinfo_old->neighbourCoord_[2][2]);
355
	      }
356
357
358
359
360

	      neighbourCoord_[1][0] = elinfo_old->neighbourCoord_[2][0];
	      neighbourCoord_[1][1].setMidpoint(elinfo_old->neighbourCoord_[2][0],
						elinfo_old->neighbourCoord_[2][1]);
	      neighbourCoord_[1][2] = oppCoord_[1];
361
	    }
362
363
	    nb = nb->getSecondChild();

364
	  } else {
365
366
	    oppVertex_[1] = 0;

367
	    if (fill_opp_coords) {
368
369
370
371
372
373
	      oppCoord_[1] = elinfo_old->oppCoord_[2];

	      neighbourCoord_[1][0] = elinfo_old->neighbourCoord_[2][2];	      
	      neighbourCoord_[1][1] = elinfo_old->neighbourCoord_[2][0];
	      neighbourCoord_[1][2].setMidpoint(elinfo_old->neighbourCoord_[2][0],
						elinfo_old->neighbourCoord_[2][1]);
374
375
376
	    }
	  }
	}
377
378
379
380
381
	neighbour_[1] = nb;
      } // if (ichild == 0) {} else
    } // if (fill_flag.isSet(Mesh::FILL_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS))
    

382
    if (fill_flag.isSet(Mesh::FILL_BOUND)) {
383
      if (elinfo_old->getBoundary(2)) {
384
	boundary_[5] = elinfo_old->getBoundary(2);
385
      } else {
386
	boundary_[5] = INTERIOR;
387
      }
388

389
      if (ichild == 0) {
390
391
392
393
394
395
396
397
398
399
400
401
402
	boundary_[3] = elinfo_old->getBoundary(5);
	boundary_[4] = elinfo_old->getBoundary(3);
	boundary_[0] = elinfo_old->getBoundary(2);
	boundary_[1] = INTERIOR;
	boundary_[2] = elinfo_old->getBoundary(1);
      } else {
	boundary_[3] = elinfo_old->getBoundary(4);
	boundary_[4] = elinfo_old->getBoundary(5);
	boundary_[0] = INTERIOR;
	boundary_[1] = elinfo_old->getBoundary(2);
	boundary_[2] = elinfo_old->getBoundary(0);
      }

403
404
405
406
407
408
      if (elinfo_old->getProjection(0) && 
	  elinfo_old->getProjection(0)->getType() == VOLUME_PROJECTION) {
	
	projection_[0] = elinfo_old->getProjection(0);
      } else { // boundary projection
	if (ichild == 0) {
409
410
411
412
413
414
415
416
417
418
419
420
	  projection_[0] = elinfo_old->getProjection(2);
	  projection_[1] = NULL;
	  projection_[2] = elinfo_old->getProjection(1);
	} else {
	  projection_[0] = NULL;
	  projection_[1] = elinfo_old->getProjection(2);
	  projection_[2] = elinfo_old->getProjection(0);
	}
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
421
  double ElInfo2d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
422
  {
423
    FUNCNAME("ElInfo2d::calcGrdLambda()");
424
425

    testFlag(Mesh::FILL_COORDS);
426
427
  
    double adet = 0.0;
428
429
    int dim = mesh_->getDim();

Thomas Witkowski's avatar
Thomas Witkowski committed
430
    for (int i = 0; i < dimOfWorld; i++) {
431
432
      (*e1)[i] = coord_[1][i] - coord_[0][i];
      (*e2)[i] = coord_[2][i] - coord_[0][i];
433
434
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
435
    if (dimOfWorld == 2) {
436
      double sdet = (*e1)[0] * (*e2)[1] - (*e1)[1] * (*e2)[0];
437
438
439
440
      adet = abs(sdet);

      if (adet < 1.0E-25) {
	MSG("abs(det) = %f\n", adet);
441
	for (int i = 0; i <= dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
442
	  for (int j = 0; j < dimOfWorld; j++)
443
	    grd_lam[i][j] = 0.0;
444
445
      } else {
	double det1 = 1.0 / sdet;
446
447
448
449
450

	grd_lam[1][0] = (*e2)[1] * det1;  // a11: (a_ij) = A^{-T}
	grd_lam[1][1] = -(*e2)[0] * det1; // a21
	grd_lam[2][0] = -(*e1)[1] * det1; // a12
	grd_lam[2][1] = (*e1)[0] * det1;  // a22
451
452
	grd_lam[0][0] = - grd_lam[1][0] - grd_lam[2][0];
	grd_lam[0][1] = - grd_lam[1][1] - grd_lam[2][1];
453
      }
454
455
    } else {  
      vectorProduct(*e1, *e2, *normal);
456

457
      adet = norm(normal);
458
459
460

      if (adet < 1.0E-15) {
	MSG("abs(det) = %lf\n", adet);
461
	for (int i = 0; i <= dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
462
	  for (int j = 0; j < dimOfWorld; j++)
463
464
	    grd_lam[i][j] = 0.0;
      } else {
465
466
	vectorProduct(*e2, *normal, grd_lam[1]);
	vectorProduct(*normal, *e1, grd_lam[2]);
467
      
468
	double adet2 = 1.0 / (adet * adet);
469

Thomas Witkowski's avatar
Thomas Witkowski committed
470
	for (int i = 0; i < dimOfWorld; i++) {
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
	  grd_lam[1][i] *= adet2;
	  grd_lam[2][i] *= adet2;
	}

	grd_lam[0][0] = - grd_lam[1][0] - grd_lam[2][0];
	grd_lam[0][1] = - grd_lam[1][1] - grd_lam[2][1];
	grd_lam[0][2] = - grd_lam[1][2] - grd_lam[2][2];
      }
    }

    return adet;
  }

  const int ElInfo2d::worldToCoord(const WorldVector<double>& xy,
				   DimVec<double>* lambda) const
  {
487
    FUNCNAME("ElInfo::worldToCoord()");
488

489
    TEST_EXIT_DBG(lambda)("lambda must not be NULL\n");
490

491
492
493
494
    DimVec<WorldVector<double> > edge(mesh_->getDim(), NO_INIT);
    WorldVector<double> x; 
    static DimVec<double> vec(mesh_->getDim(), NO_INIT);

495
496
    int dim = mesh_->getDim();

497
498
    for (int j = 0; j < dimOfWorld; j++) {
      double x0 = coord_[dim][j];
499
      x[j] = xy[j] - x0;
500
      for (int i = 0; i < dim; i++)
501
502
503
	edge[i][j] = coord_[i][j] - x0;
    }
  
504
505
506
    double det  = edge[0][0] * edge[1][1] - edge[0][1] * edge[1][0]; 
    double det0 =       x[0] * edge[1][1] -       x[1] * edge[1][0]; 
    double det1 = edge[0][0] * x[1]       - edge[0][1] * x[0]; 
507
508
509

    if (abs(det) < DBL_TOL) {
      ERROR("det = %le; abort\n", det);
510
511
      for (int i = 0; i <= dim; i++) 
	(*lambda)[i] = 1.0/dim;
512
513
514
515
516
517
518
519
520
      return 0;
    }

    (*lambda)[0] = det0 / det;
    (*lambda)[1] = det1 / det;
    (*lambda)[2] = 1.0 - (*lambda)[0] - (*lambda)[1];

    int k = -1;
    double lmin = 0.0;
521
    for (int i = 0; i <= dim; i++) {
522
523
524
525
526
527
528
529
530
531
532
533
      if ((*lambda)[i] < -1.E-5) {
	if ((*lambda)[i] < lmin) {
	  k = i;
	  lmin = (*lambda)[i];
	}
      }
    }

    return k;
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
534
  double ElInfo2d::getNormal(int side, WorldVector<double> &normal)
535
  {
536
    FUNCNAME("ElInfo::getNormal()");
537

538
539
    int i0 = (side + 1) % 3;
    int i1 = (side + 2) % 3;
540

Thomas Witkowski's avatar
Thomas Witkowski committed
541
    if (dimOfWorld == 2){
542
543
544
545
546
      normal[0] = coord_[i1][1] - coord_[i0][1];
      normal[1] = coord_[i0][0] - coord_[i1][0];
    } else { // dow == 3
      WorldVector<double> e0, e1,e2, elementNormal;

547
548
549
550
551
552
      e0 = coord_[i1]; 
      e0 -= coord_[i0];
      e1 = coord_[i1]; 
      e1 -= coord_[side];
      e2 = coord_[i0]; 
      e2 -= coord_[side];
553
554
555
556
557

      vectorProduct(e1, e2, elementNormal);
      vectorProduct(elementNormal, e0, normal);
    }

558
    double det = norm(&normal);
559

560
    TEST_EXIT_DBG(det > 1.e-30)("det = 0 on face %d\n", side);
561
562
563
564
565
566
567
568
569
570
571
572
573
574

    normal *= 1.0 / det;
    
    return(det);
  }


  /****************************************************************************/
  /*  calculate the normal of the element for dim of world = dim + 1          */
  /*  return the absulute value of the determinant from the                   */
  /*  transformation to the reference element                                 */
  /****************************************************************************/
  double ElInfo2d::getElementNormal( WorldVector<double> &elementNormal) const
  {
575
    FUNCNAME("ElInfo::getElementNormal()");
576

Thomas Witkowski's avatar
Thomas Witkowski committed
577
578
    TEST_EXIT_DBG(dimOfWorld == 3)
      (" element normal only well defined for  DIM_OF_WORLD = DIM + 1 !!");
579

580
581
    WorldVector<double> e0 = coord_[1] - coord_[0];
    WorldVector<double> e1 = coord_[2] - coord_[0];
582
583
584

    vectorProduct(e0, e1, elementNormal);

585
    double det = norm(&elementNormal);
586

587
    TEST_EXIT_DBG(det > 1.e-30)("det = 0");
588
589
590
591
592
593

    elementNormal *= 1.0 / det;
    
    return(det);
  }
}