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
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
640

  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][0] = (*coords)[0][1];
      (*coords)[1][0] = (*coords)[1][1];
Thomas Witkowski's avatar
Thomas Witkowski committed
641
642
643
644
645
      (*coords)[2][0] = (*coords)[2][1];  

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

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

653
}