ElInfo2d.cc 18 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
	  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) {
76
	      if (nb->isNewCoordSet()) {
77
78
79
80
81
82
83
84
85
86
87
88
89
		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 {
Thomas Witkowski's avatar
Thomas Witkowski committed
91
		  ERROR_EXIT("Should not happen!\n");
92
		}
93
94
95
96
97
98
99
100
101
102
103
104
	
		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 {
Thomas Witkowski's avatar
Thomas Witkowski committed
105
		  ERROR_EXIT("Should not happen!\n");
106
107
108
109
110
		}
		
		neighbourCoord_[i][2] = oppCoord_[i];
		break;
		
Thomas Witkowski's avatar
Thomas Witkowski committed
111
112
113
114
115
116
117
118
119
120
121
122
	      case 2:
		if (*(macroNeighbour->getElement()->getDOF(2)) == *(element_->getDOF(0))) {
		  neighbourCoord_[i][0] = macroNeighbour->coord[2];
		  neighbourCoord_[i][1] = macroNeighbour->coord[1];
		} else if (*(macroNeighbour->getElement()->getDOF(2)) == *(element_->getDOF(1))) {
		  neighbourCoord_[i][0] = macroNeighbour->coord[0];
		  neighbourCoord_[i][1] = macroNeighbour->coord[2];		 
		} else {
		  ERROR_EXIT("Should not happen!\n");
		}
		break;

123
	      default:
Thomas Witkowski's avatar
Thomas Witkowski committed
124
		std::cout << "------------- Error --------------" << std::endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
125
126
		std::cout << "  Neighbour counter = " << i << "\n";
		std::cout << "  Element index     = " << element_->getIndex() << "\n\n";
Thomas Witkowski's avatar
Thomas Witkowski committed
127
128
129
130
131
132
133
134
135
136
137
		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;
		}
138
139
		ERROR_EXIT("should not happen!\n");
		break;
140
141
142
	      }
	    }
	  } else {
143
144
145
146
147
	    if (fill_opp_coords) {
	      oppCoord_[i] = macroNeighbour->coord[edgeNo];

	      neighbourCoord_[i] = macroNeighbour->coord;	      
	    }
148
	  }
149
150
151
	} else {
	  neighbour_[i] = NULL;
        }
152
      }
153
    }
154
    
155
156
    if (fillFlag_.isSet(Mesh::FILL_BOUND)) {   
      for (int i = 0; i < element_->getGeo(BOUNDARY); i++) {
157
158
159
	boundary_[i] = mel->getBoundary(i);
      }

160
      for (int i = 0; i < element_->getGeo(PROJECTION); i++) {
161
162
163
164
165
166
167
168
169
170
	projection_[i] = mel->getProjection(i);
      }
    }
  }


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

171
  void ElInfo2d::fillElInfo(int ichild, const ElInfo *elInfoOld)
172
  {
173
    FUNCNAME("ElInfo::fillElInfo()");
174

175
    Element *elem = elInfoOld->element_;
176
177
    Element *nb;

178
    Flag fill_flag = elInfoOld->fillFlag_;
179

180
181
182
183
184
    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);
185

186
    macroElement_  = elInfoOld->macroElement_;
187
    fillFlag_ = fill_flag;
188
    parent_ = elem;
189
    level = elInfoOld->level + 1;
190
    iChild = ichild;
191
192
193

    if (fillFlag_.isSet(Mesh::FILL_COORDS) || 
	fillFlag_.isSet(Mesh::FILL_DET)    ||
194
195
	fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
      
196
      if (elem->isNewCoordSet()) {
197
198
	coord_[2] = *(elem->getNewCoord());
      } else {
199
	coord_[2].setMidpoint(elInfoOld->coord_[0], elInfoOld->coord_[1]);
200
201
202
      }
      
      if (ichild == 0) {
203
204
	coord_[0] = elInfoOld->coord_[2];
	coord_[1] = elInfoOld->coord_[0];
205
      } else {
206
207
	coord_[0] = elInfoOld->coord_[1];
	coord_[1] = elInfoOld->coord_[2];
208
209
210
211
212
213
214
215
216
217
      }
    }

    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.

218
219
	neighbour_[2] = elInfoOld->neighbour_[1];
	oppVertex_[2] = elInfoOld->oppVertex_[1];
220
221
	
	if (neighbour_[2] && fill_opp_coords) {
222
223
	  oppCoord_[2] = elInfoOld->oppCoord_[1];
	  neighbourCoord_[2] = elInfoOld->neighbourCoord_[1];
224
	}
225
226
227
228
229
230
231
232
233
234
235
236
237
	
	
	// 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) {
238
            if (elem->getSecondChild()->isNewCoordSet()) {
239
240
	      oppCoord_[1] = *(elem->getSecondChild()->getNewCoord());
	    } else {      
241
242
	      oppCoord_[1].setMidpoint(elInfoOld->coord_[1], 
				       elInfoOld->coord_[2]);
243
	    }
244

245
246
247
	    neighbourCoord_[1][0] = coord_[0];
	    neighbourCoord_[1][1] = coord_[2];
	    neighbourCoord_[1][2] = oppCoord_[1];  
248
249
	  }
	} else {
250
251
252
253
	  neighbour_[1] = elem->getSecondChild();
	  oppVertex_[1] = 0;

	  if (fill_opp_coords) {
254
	    oppCoord_[1] = elInfoOld->coord_[1];
255

256
257
	    neighbourCoord_[1][0] = elInfoOld->coord_[1];
	    neighbourCoord_[1][1] = elInfoOld->coord_[2];
258
	    neighbourCoord_[1][2] = coord_[2];
259
260
261
262
	  }
	}


263
264
265
	// Calculation of the neighbour 0, its oppCoords and the
	// cooresponding oppVertex.
	
266
	nb = elInfoOld->neighbour_[2];
267
	if (nb) {
268
	  TEST(elInfoOld->oppVertex_[2] == 2)("invalid neighbour\n"); 
269
270
	  TEST_EXIT_DBG(nb->getFirstChild())("missing first child?\n");
	  TEST_EXIT_DBG(nb->getSecondChild())("missing second child?\n");
271
272
273
274
275
276
	 
	  nb = nb->getSecondChild();

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

277
	    if (fill_opp_coords) {
278
	      if (nb->isNewCoordSet()) {
279
280
		oppCoord_[0] = *(nb->getNewCoord());
	      } else {
281
282
		oppCoord_[0].setMidpoint(elInfoOld->neighbourCoord_[2][1],
					 elInfoOld->neighbourCoord_[2][2]);
283
	      }
284

285
286
287
	      neighbourCoord_[0][0].setMidpoint(elInfoOld->neighbourCoord_[2][0],
						elInfoOld->neighbourCoord_[2][1]);
	      neighbourCoord_[0][1] = elInfoOld->neighbourCoord_[2][1];
288
289
290
291
292
293
294
	      neighbourCoord_[0][2] = oppCoord_[0];
	    }	   
 
	    nb = nb->getFirstChild();
	  } else {
	    oppVertex_[0] = 1;

295
	    if (fill_opp_coords) {
296
	      oppCoord_[0] = elInfoOld->oppCoord_[2];    
297

298
299
300
301
	      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]);
302
303
	    }
	  }
304
305
306
307
308
309
310
	}
	
	neighbour_[0] = nb;
      } else {   /* ichild == 1 */
	// Calculation of the neighbour 2, its oppCoords and the
	// cooresponding oppVertex.

311
312
	neighbour_[2] = elInfoOld->neighbour_[0];
	oppVertex_[2] = elInfoOld->oppVertex_[0];
313
314

	if (neighbour_[2] && fill_opp_coords) {
315
316
	  oppCoord_[2] = elInfoOld->oppCoord_[0];
	  neighbourCoord_[2] = elInfoOld->neighbourCoord_[0];
317
318
319
320
321
322
323
324
325
326
327
	}
	

	// 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) {
328
            if (elem->getFirstChild()->isNewCoordSet()) {
329
	      oppCoord_[0] = *(elem->getFirstChild()->getNewCoord());
330
	    } else {
331
332
	      oppCoord_[0].setMidpoint(elInfoOld->coord_[0], 
				       elInfoOld->coord_[2]);
333
	    }
334
335
336
337

	    neighbourCoord_[0][0] = coord_[2];
	    neighbourCoord_[0][1] = coord_[1];
	    neighbourCoord_[0][2] = oppCoord_[0];
338
	  }
339
340
341
342
343
	} else {
	  neighbour_[0] = elem->getFirstChild();
	  oppVertex_[0] = 1;

	  if (fill_opp_coords) {
344
	    oppCoord_[0] = elInfoOld->coord_[0];
345

346
347
	    neighbourCoord_[0][0] = elInfoOld->coord_[2];
	    neighbourCoord_[0][1] = elInfoOld->coord_[0];
348
	    neighbourCoord_[0][2] = coord_[2];
349
	  }
350
351
352
353
354
	}

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

355
	nb = elInfoOld->neighbour_[2];
356
	if (nb) {
357
	  TEST(elInfoOld->oppVertex_[2] == 2)("invalid neighbour\n"); 
358
359
360
361
362
	  TEST((nb = nb->getFirstChild()))("missing child?\n");

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

363
	    if (fill_opp_coords) {
364
	      if (nb->isNewCoordSet()) {
365
		oppCoord_[1] = *(nb->getNewCoord());
366
	      } else {
367
368
		oppCoord_[1].setMidpoint(elInfoOld->neighbourCoord_[2][0],
					 elInfoOld->neighbourCoord_[2][2]);
369
	      }
370

371
372
373
	      neighbourCoord_[1][0] = elInfoOld->neighbourCoord_[2][0];
	      neighbourCoord_[1][1].setMidpoint(elInfoOld->neighbourCoord_[2][0],
						elInfoOld->neighbourCoord_[2][1]);
374
	      neighbourCoord_[1][2] = oppCoord_[1];
375
	    }
376
377
	    nb = nb->getSecondChild();

378
	  } else {
379
380
	    oppVertex_[1] = 0;

381
	    if (fill_opp_coords) {
382
	      oppCoord_[1] = elInfoOld->oppCoord_[2];
383

384
385
386
387
	      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]);
388
389
390
	    }
	  }
	}
391
392
393
394
395
	neighbour_[1] = nb;
      } // if (ichild == 0) {} else
    } // if (fill_flag.isSet(Mesh::FILL_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS))
    

396
    if (fill_flag.isSet(Mesh::FILL_BOUND)) {
397
398
      if (elInfoOld->getBoundary(2)) {
	boundary_[5] = elInfoOld->getBoundary(2);
399
      } else {
400
	boundary_[5] = INTERIOR;
401
      }
402

403
      if (ichild == 0) {
404
405
406
	boundary_[3] = elInfoOld->getBoundary(5);
	boundary_[4] = elInfoOld->getBoundary(3);
	boundary_[0] = elInfoOld->getBoundary(2);
407
	boundary_[1] = INTERIOR;
408
	boundary_[2] = elInfoOld->getBoundary(1);
409
      } else {
410
411
	boundary_[3] = elInfoOld->getBoundary(4);
	boundary_[4] = elInfoOld->getBoundary(5);
412
	boundary_[0] = INTERIOR;
413
414
	boundary_[1] = elInfoOld->getBoundary(2);
	boundary_[2] = elInfoOld->getBoundary(0);
415
416
      }

417
418
      if (elInfoOld->getProjection(0) && 
	  elInfoOld->getProjection(0)->getType() == VOLUME_PROJECTION) {
419
	
420
	projection_[0] = elInfoOld->getProjection(0);
421
422
      } else { // boundary projection
	if (ichild == 0) {
423
	  projection_[0] = elInfoOld->getProjection(2);
424
	  projection_[1] = NULL;
425
	  projection_[2] = elInfoOld->getProjection(1);
426
427
	} else {
	  projection_[0] = NULL;
428
429
	  projection_[1] = elInfoOld->getProjection(2);
	  projection_[2] = elInfoOld->getProjection(0);
430
431
432
433
434
	}
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
435
  double ElInfo2d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
436
  {
437
    FUNCNAME("ElInfo2d::calcGrdLambda()");
438
439

    testFlag(Mesh::FILL_COORDS);
440
441
  
    double adet = 0.0;
442
443
    int dim = mesh_->getDim();

Thomas Witkowski's avatar
Thomas Witkowski committed
444
    for (int i = 0; i < dimOfWorld; i++) {
445
446
      (*e1)[i] = coord_[1][i] - coord_[0][i];
      (*e2)[i] = coord_[2][i] - coord_[0][i];
447
448
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
449
    if (dimOfWorld == 2) {
450
      double sdet = (*e1)[0] * (*e2)[1] - (*e1)[1] * (*e2)[0];
451
452
453
454
      adet = abs(sdet);

      if (adet < 1.0E-25) {
	MSG("abs(det) = %f\n", adet);
455
	for (int i = 0; i <= dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
456
	  for (int j = 0; j < dimOfWorld; j++)
457
	    grd_lam[i][j] = 0.0;
458
459
      } else {
	double det1 = 1.0 / sdet;
460
461
462
463
464

	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
465
466
	grd_lam[0][0] = - grd_lam[1][0] - grd_lam[2][0];
	grd_lam[0][1] = - grd_lam[1][1] - grd_lam[2][1];
467
      }
468
469
    } else {  
      vectorProduct(*e1, *e2, *normal);
470

471
      adet = norm(normal);
472
473
474

      if (adet < 1.0E-15) {
	MSG("abs(det) = %lf\n", adet);
475
	for (int i = 0; i <= dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
476
	  for (int j = 0; j < dimOfWorld; j++)
477
478
	    grd_lam[i][j] = 0.0;
      } else {
479
480
	vectorProduct(*e2, *normal, grd_lam[1]);
	vectorProduct(*normal, *e1, grd_lam[2]);
481
      
482
	double adet2 = 1.0 / (adet * adet);
483

Thomas Witkowski's avatar
Thomas Witkowski committed
484
	for (int i = 0; i < dimOfWorld; i++) {
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
	  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
  {
501
    FUNCNAME("ElInfo::worldToCoord()");
502

503
    TEST_EXIT_DBG(lambda)("lambda must not be NULL\n");
504

505
506
507
508
    DimVec<WorldVector<double> > edge(mesh_->getDim(), NO_INIT);
    WorldVector<double> x; 
    static DimVec<double> vec(mesh_->getDim(), NO_INIT);

509
510
    int dim = mesh_->getDim();

511
512
    for (int j = 0; j < dimOfWorld; j++) {
      double x0 = coord_[dim][j];
513
      x[j] = xy[j] - x0;
514
      for (int i = 0; i < dim; i++)
515
516
517
	edge[i][j] = coord_[i][j] - x0;
    }
  
518
519
520
    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]; 
521
522
523

    if (abs(det) < DBL_TOL) {
      ERROR("det = %le; abort\n", det);
524
525
      for (int i = 0; i <= dim; i++) 
	(*lambda)[i] = 1.0/dim;
526
527
528
529
530
531
532
533
534
      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;
535
    for (int i = 0; i <= dim; i++) {
536
537
538
539
540
541
542
543
544
545
546
547
      if ((*lambda)[i] < -1.E-5) {
	if ((*lambda)[i] < lmin) {
	  k = i;
	  lmin = (*lambda)[i];
	}
      }
    }

    return k;
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
548
  double ElInfo2d::getNormal(int side, WorldVector<double> &normal)
549
  {
550
    FUNCNAME("ElInfo::getNormal()");
551

552
553
    int i0 = (side + 1) % 3;
    int i1 = (side + 2) % 3;
554

Thomas Witkowski's avatar
Thomas Witkowski committed
555
    if (dimOfWorld == 2){
556
557
558
559
560
      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;

561
562
563
564
565
566
      e0 = coord_[i1]; 
      e0 -= coord_[i0];
      e1 = coord_[i1]; 
      e1 -= coord_[side];
      e2 = coord_[i0]; 
      e2 -= coord_[side];
567
568
569
570
571

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

572
    double det = norm(&normal);
573

574
    TEST_EXIT_DBG(det > 1.e-30)("det = 0 on face %d\n", side);
575
576
577

    normal *= 1.0 / det;
    
Thomas Witkowski's avatar
Thomas Witkowski committed
578
    return det;
579
580
581
582
583
584
585
586
  }


  /****************************************************************************/
  /*  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                                 */
  /****************************************************************************/
Thomas Witkowski's avatar
Thomas Witkowski committed
587
  double ElInfo2d::getElementNormal(WorldVector<double> &elementNormal) const
588
  {
589
    FUNCNAME("ElInfo::getElementNormal()");
590

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

594
595
    WorldVector<double> e0 = coord_[1] - coord_[0];
    WorldVector<double> e1 = coord_[2] - coord_[0];
596
597
598

    vectorProduct(e0, e1, elementNormal);

599
    double det = norm(&elementNormal);
600

601
    TEST_EXIT_DBG(det > 1.e-30)("det = 0");
602
603
604

    elementNormal *= 1.0 / det;
    
Thomas Witkowski's avatar
Thomas Witkowski committed
605
    return det;
606
  }
607

608
609
  void ElInfo2d::getRefSimplexCoords(const BasisFunction *basisFcts,
				     DimMat<double> *coords) const
610
611
612
613
614
615
616
617
618
619
620
621
622
623
  {
    (*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;
  }

624
625
626
  void ElInfo2d::getSubElementCoords(const BasisFunction *basisFcts,
				     int iChild,
				     DimMat<double> *coords) const
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
  {
    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][0] = (*coords)[0][1];
      (*coords)[1][0] = (*coords)[1][1];
Thomas Witkowski's avatar
Thomas Witkowski committed
643
644
645
646
647
      (*coords)[2][0] = (*coords)[2][1];  

      (*coords)[0][1] = (*coords)[0][2];
      (*coords)[1][1] = (*coords)[1][2];
      (*coords)[2][1] = (*coords)[2][2];   
648
649
650
651
652
653
654
    }

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

655
}