ElInfo2d.cc 17.4 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
    element_  = const_cast<Element*>(mel->getElement());
    parent_  = NULL;
39
    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
	projection_[i] = mel->getProjection(i);
      }
    }
  }


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

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

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

165
    Flag fill_flag = elInfoOld->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
    macroElement_  = elInfoOld->macroElement_;
174
    fillFlag_ = fill_flag;
175
    parent_ = elem;
176
    level = elInfoOld->level + 1;
177
    iChild = ichild;
178
179
180

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

    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.

205
206
	neighbour_[2] = elInfoOld->neighbour_[1];
	oppVertex_[2] = elInfoOld->oppVertex_[1];
207
208
	
	if (neighbour_[2] && fill_opp_coords) {
209
210
	  oppCoord_[2] = elInfoOld->oppCoord_[1];
	  neighbourCoord_[2] = elInfoOld->neighbourCoord_[1];
211
	}
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
	
	
	// 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 {      
228
229
	      oppCoord_[1].setMidpoint(elInfoOld->coord_[1], 
				       elInfoOld->coord_[2]);
230
	    }
231

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

	  if (fill_opp_coords) {
241
	    oppCoord_[1] = elInfoOld->coord_[1];
242

243
244
	    neighbourCoord_[1][0] = elInfoOld->coord_[1];
	    neighbourCoord_[1][1] = elInfoOld->coord_[2];
245
	    neighbourCoord_[1][2] = coord_[2];
246
247
248
249
	  }
	}


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

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

264
	    if (fill_opp_coords) {
265
266
267
	      if (nb->getNewCoord(-1)) {
		oppCoord_[0] = *(nb->getNewCoord());
	      } else {
268
269
		oppCoord_[0].setMidpoint(elInfoOld->neighbourCoord_[2][1],
					 elInfoOld->neighbourCoord_[2][2]);
270
	      }
271

272
273
274
	      neighbourCoord_[0][0].setMidpoint(elInfoOld->neighbourCoord_[2][0],
						elInfoOld->neighbourCoord_[2][1]);
	      neighbourCoord_[0][1] = elInfoOld->neighbourCoord_[2][1];
275
276
277
278
279
280
281
	      neighbourCoord_[0][2] = oppCoord_[0];
	    }	   
 
	    nb = nb->getFirstChild();
	  } else {
	    oppVertex_[0] = 1;

282
	    if (fill_opp_coords) {
283
	      oppCoord_[0] = elInfoOld->oppCoord_[2];    
284

285
286
287
288
	      neighbourCoord_[0][0] = elInfoOld->neighbourCoord_[2][0];
	      neighbourCoord_[0][1] = elInfoOld->neighbourCoord_[2][2];
	      neighbourCoord_[0][2].setMidpoint(elInfoOld->neighbourCoord_[2][0],
						elInfoOld->neighbourCoord_[2][1]);
289
290
	    }
	  }
291
292
293
294
295
296
297
	}
	
	neighbour_[0] = nb;
      } else {   /* ichild == 1 */
	// Calculation of the neighbour 2, its oppCoords and the
	// cooresponding oppVertex.

298
299
	neighbour_[2] = elInfoOld->neighbour_[0];
	oppVertex_[2] = elInfoOld->oppVertex_[0];
300
301

	if (neighbour_[2] && fill_opp_coords) {
302
303
	  oppCoord_[2] = elInfoOld->oppCoord_[0];
	  neighbourCoord_[2] = elInfoOld->neighbourCoord_[0];
304
305
306
307
308
309
310
311
312
313
314
315
316
	}
	

	// 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());
317
	    } else {
318
319
	      oppCoord_[0].setMidpoint(elInfoOld->coord_[0], 
				       elInfoOld->coord_[2]);
320
	    }
321
322
323
324

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

	  if (fill_opp_coords) {
331
	    oppCoord_[0] = elInfoOld->coord_[0];
332

333
334
	    neighbourCoord_[0][0] = elInfoOld->coord_[2];
	    neighbourCoord_[0][1] = elInfoOld->coord_[0];
335
	    neighbourCoord_[0][2] = coord_[2];
336
	  }
337
338
339
340
341
	}

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

342
	nb = elInfoOld->neighbour_[2];
343
	if (nb) {
344
	  TEST(elInfoOld->oppVertex_[2] == 2)("invalid neighbour\n"); 
345
346
347
348
349
	  TEST((nb = nb->getFirstChild()))("missing child?\n");

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

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

358
359
360
	      neighbourCoord_[1][0] = elInfoOld->neighbourCoord_[2][0];
	      neighbourCoord_[1][1].setMidpoint(elInfoOld->neighbourCoord_[2][0],
						elInfoOld->neighbourCoord_[2][1]);
361
	      neighbourCoord_[1][2] = oppCoord_[1];
362
	    }
363
364
	    nb = nb->getSecondChild();

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

368
	    if (fill_opp_coords) {
369
	      oppCoord_[1] = elInfoOld->oppCoord_[2];
370

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

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

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

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

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

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

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

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

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

	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
452
453
	grd_lam[0][0] = - grd_lam[1][0] - grd_lam[2][0];
	grd_lam[0][1] = - grd_lam[1][1] - grd_lam[2][1];
454
      }
455
456
    } else {  
      vectorProduct(*e1, *e2, *normal);
457

458
      adet = norm(normal);
459
460
461

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

Thomas Witkowski's avatar
Thomas Witkowski committed
471
	for (int i = 0; i < dimOfWorld; i++) {
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
	  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
  {
488
    FUNCNAME("ElInfo::worldToCoord()");
489

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

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

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

498
499
    for (int j = 0; j < dimOfWorld; j++) {
      double x0 = coord_[dim][j];
500
      x[j] = xy[j] - x0;
501
      for (int i = 0; i < dim; i++)
502
503
504
	edge[i][j] = coord_[i][j] - x0;
    }
  
505
506
507
    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]; 
508
509
510

    if (abs(det) < DBL_TOL) {
      ERROR("det = %le; abort\n", det);
511
512
      for (int i = 0; i <= dim; i++) 
	(*lambda)[i] = 1.0/dim;
513
514
515
516
517
518
519
520
521
      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;
522
    for (int i = 0; i <= dim; i++) {
523
524
525
526
527
528
529
530
531
532
533
534
      if ((*lambda)[i] < -1.E-5) {
	if ((*lambda)[i] < lmin) {
	  k = i;
	  lmin = (*lambda)[i];
	}
      }
    }

    return k;
  }


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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
542
    if (dimOfWorld == 2){
543
544
545
546
547
      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;

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

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

559
    double det = norm(&normal);
560

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

    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
  {
576
    FUNCNAME("ElInfo::getElementNormal()");
577

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

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

    vectorProduct(e0, e1, elementNormal);

586
    double det = norm(&elementNormal);
587

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

    elementNormal *= 1.0 / det;
    
    return(det);
  }
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639

  void ElInfo2d::getRefSimplexCoords(DimMat<double> *coords) const
  {
    (*coords)[0][0] = 1.0;
    (*coords)[1][0] = 0.0;
    (*coords)[2][0] = 0.0;

    (*coords)[0][1] = 0.0;
    (*coords)[1][1] = 1.0;
    (*coords)[2][1] = 0.0;

    (*coords)[0][2] = 0.0;
    (*coords)[1][2] = 0.0;
    (*coords)[2][2] = 1.0;
  }

  void ElInfo2d::getSubElementCoords(DimMat<double> *coords,
				     int iChild) const
  {
    double c0 = ((*coords)[0][0] + (*coords)[0][1]) * 0.5;
    double c1 = ((*coords)[1][0] + (*coords)[1][1]) * 0.5;
    double c2 = ((*coords)[2][0] + (*coords)[2][1]) * 0.5;

    if (iChild == 0) {
      (*coords)[0][1] = (*coords)[0][0];
      (*coords)[1][1] = (*coords)[1][0];
      (*coords)[2][1] = (*coords)[2][0];

      (*coords)[0][0] = (*coords)[0][2];
      (*coords)[1][0] = (*coords)[1][2];
      (*coords)[2][0] = (*coords)[2][2];
    } else {
      (*coords)[0][1] = (*coords)[0][2];
      (*coords)[1][1] = (*coords)[1][2];
      (*coords)[2][1] = (*coords)[2][2];

      (*coords)[0][0] = (*coords)[0][1];
      (*coords)[1][0] = (*coords)[1][1];
      (*coords)[2][0] = (*coords)[2][1];      
    }

    (*coords)[0][2] = c0;
    (*coords)[1][2] = c1;
    (*coords)[2][2] = c2;
  }

640
}