ElInfo2d.cc 19.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
  double ElInfo2d::mat_d1_val[3][3] = {{1.0, 0.0, 0.0}, 
				       {0.0, 1.0, 0.0}, 
				       {0.0, 0.0, 1.0}};
  mtl::dense2D<double> ElInfo2d::mat_d1(mat_d1_val);

  double ElInfo2d::mat_d1_left_val[3][3] = {{0.0, 1.0, 0.5}, 
					    {0.0, 0.0, 0.5},
					    {1.0, 0.0, 0.0}};
  mtl::dense2D<double> ElInfo2d::mat_d1_left(mat_d1_left_val);

  double ElInfo2d::mat_d1_right_val[3][3] = {{0.0, 0.0, 0.5}, 
					     {1.0, 0.0, 0.5},
					     {0.0, 1.0, 0.0}};
  mtl::dense2D<double> ElInfo2d::mat_d1_right(mat_d1_right_val);

32
33
34
  ElInfo2d::ElInfo2d(Mesh *aMesh) 
    : ElInfo(aMesh) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
35
36
37
    e1 = new WorldVector<double>;
    e2 = new WorldVector<double>;
    normal = new WorldVector<double>;
38
39
40
41
  }

  ElInfo2d::~ElInfo2d()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
42
43
44
    delete e1;
    delete e2;
    delete normal;
45
46
  }

47
48
  void ElInfo2d::fillMacroInfo(const MacroElement * mel)
  {
49
50
    FUNCNAME("ElInfo::fillMacroInfo()");
 
51
    macroElement_ = const_cast<MacroElement*>(mel);
52
53
    element_ = const_cast<Element*>(mel->getElement());
    parent_ = NULL;
54
    level = 0;
55
56
57
58

    if (fillFlag_.isSet(Mesh::FILL_COORDS) || 
	fillFlag_.isSet(Mesh::FILL_DET)    ||
	fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
59
60

      int vertices = mesh_->getGeo(VERTEX);
61
      for (int i = 0; i < vertices; i++)
62
63
64
65
66
	coord_[i] = mel->coord[i];
    }

    int neighbours = mesh_->getGeo(NEIGH);

67
68
69
70
    if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS) || 
	fillFlag_.isSet(Mesh::FILL_NEIGH)) {

      bool fill_opp_coords = (fillFlag_.isSet(Mesh::FILL_OPP_COORDS));
71
    
72
73
74
75
      for (int i = 0; i < neighbours; i++) {
	MacroElement *macroNeighbour = mel->getNeighbour(i);

	if (macroNeighbour) {
Thomas Witkowski's avatar
Thomas Witkowski committed
76
	  neighbour_[i] = macroNeighbour->getElement();	  
77
78
	  Element *nb = const_cast<Element*>(neighbour_[i]);

79
	  int edgeNo = oppVertex[i] = mel->getOppVertex(i);
80
81


82
	  if (nb->getFirstChild() && edgeNo != 2) {  
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102

	    // Search for the next neighbour. In many cases, the neighbour element 
	    // may be refinemed in a way, such that there is no new vertex on the 
	    // common edge. This situation is shown in the following picture: 
	    //
	    //               /|\
	    //              / | \
	    //             /  |  \
	    //            /\  |   \
	    //           /  \ |    \ 
	    //          /    \|     \
	    //          -------------
	    //
	    //            nb     el
	    //
	    // Note that we know (because of the last if statement), that the 
	    // neighbour element has children and the common edge is not the 
	    // refinement edge, which has always the number 2, of our element.


103
	    if (edgeNo == 0) {
104
105
106
107
108
109
110
111
112
113
114
115
	      // The situation is as follows:
	      //
	      //          -------
	      //          \    /|\
	      //           \  / | \
	      //            \/  |  \
	      //             \  |   \
	      //              \ |    \ 
	      //               \|     \
	      //                -------
	      //
	      //            nb     el
116
              //
117
118
119
120
	      // That means, the edge 0 of the same level neighbour is the common
	      // edge, i.e., the direct neighbour is the second child of the same
	      // level neighbour.

121
122
	      nb = neighbour_[i] = nb->getSecondChild();
	    } else {
123
124
	      // The situation is as shown in the picture above. So the next
	      // neighbour is the first child of the same level neighbour element.
125
126
127
	      nb = neighbour_[i] = nb->getFirstChild();
	    }

128
129
	    // In both cases the opp vertex number is 2, as one can see in the 
	    // pictures above.
130
	    oppVertex[i] = 2;
131
132

	    if (fill_opp_coords) {
133
	      if (nb->isNewCoordSet()) {
134
135
		oppCoord_[i] = *(nb->getNewCoord());
	      } else {
136
137
138
139
140
141
		// In both cases, that are shown in the pictures above, the opp
		// vertex of the neighbour edge is the midpoint of the vertex 0
		// and vertex 1 of the same level neighbour element.
		oppCoord_[i] = (macroNeighbour->coord[0] + 
				macroNeighbour->coord[1]) * 0.5;
	      }
142
143
144
	      
	      switch (i) {
	      case 0:
145
146
147
148
		// The common edge is the edge 0 of this element.

		switch (edgeNo) {
		case 1:
149
150
		  neighbourCoord_[i][0] = macroNeighbour->coord[2];
		  neighbourCoord_[i][1] = macroNeighbour->coord[0];
151
152
		  break;
		case 0:		  
153
154
		  neighbourCoord_[i][0] = macroNeighbour->coord[1];
		  neighbourCoord_[i][1] = macroNeighbour->coord[2];
155
156
		  break;
		default:
Thomas Witkowski's avatar
Thomas Witkowski committed
157
		  ERROR_EXIT("Should not happen!\n");
158
		}
159
160
161
162
163
	
		neighbourCoord_[i][2] = oppCoord_[i];
		break;
		
	      case 1:
164
165
166
		// The commonedge is the edge 1 of this element.
		switch (edgeNo) {
		case 0:
167
168
		  neighbourCoord_[i][0] = macroNeighbour->coord[1];
		  neighbourCoord_[i][1] = macroNeighbour->coord[2];
169
170
		  break;
		case 1:
171
172
		  neighbourCoord_[i][0] = macroNeighbour->coord[2];
		  neighbourCoord_[i][1] = macroNeighbour->coord[0];
173
174
		  break;
		default:
Thomas Witkowski's avatar
Thomas Witkowski committed
175
		  ERROR_EXIT("Should not happen!\n");
176
177
178
179
180
		}
		
		neighbourCoord_[i][2] = oppCoord_[i];
		break;
		
Thomas Witkowski's avatar
Thomas Witkowski committed
181
	      case 2:
182
183
184
185
186
187
188
189
190
191
		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");
		}

192
193
194
		// I've deleted here some code, be I think that this case is not
		// possible. If an error occurs in this line, please check AMDiS
		// revision <= 476 at the same position.
195
196
		//		ERROR_EXIT("Should not happen!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
197
198
		break;

199
	      default:
Thomas Witkowski's avatar
Thomas Witkowski committed
200
		std::cout << "------------- Error --------------" << std::endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
201
202
		std::cout << "  Neighbour counter = " << i << "\n";
		std::cout << "  Element index     = " << element_->getIndex() << "\n\n";
Thomas Witkowski's avatar
Thomas Witkowski committed
203
204
205
206
207
208
209
210
		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;
		  }
211
212
213
		  std::cout << "  OppVertex " << j << ": " 
			    << static_cast<int>(mel->getOppVertex(j)) 
			    << std::endl << std::endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
214
		}
215
216
		ERROR_EXIT("should not happen!\n");
		break;
217
218
219
	      }
	    }
	  } else {
220
221
222
223
224
225
226

	    // In this case, we know that the common edge is the refinement edge.
	    // This makes everything much more simpler, because we know that the
	    // next neighbour is equal to the samel level neighbour. If the same
	    // level neighbour would be refinement, also this element must to be 
	    // refinement, because they share the refinement edge.

227
228
229
230
	    if (fill_opp_coords) {
	      oppCoord_[i] = macroNeighbour->coord[edgeNo];
	      neighbourCoord_[i] = macroNeighbour->coord;	      
	    }
231
	  }
232
233
234
	} else {
	  neighbour_[i] = NULL;
        }
235
      }
236
    }
237
    
238
239
    if (fillFlag_.isSet(Mesh::FILL_BOUND)) {   
      for (int i = 0; i < element_->getGeo(BOUNDARY); i++) {
240
241
242
	boundary_[i] = mel->getBoundary(i);
      }

243
      for (int i = 0; i < element_->getGeo(PROJECTION); i++) {
244
245
246
247
248
249
250
251
252
253
	projection_[i] = mel->getProjection(i);
      }
    }
  }


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

254
  void ElInfo2d::fillElInfo(int ichild, const ElInfo *elInfoOld)
255
  {
256
    FUNCNAME("ElInfo::fillElInfo()");
257

258
    Element *elem = elInfoOld->element_;
259
260
    Element *nb;

261
    Flag fill_flag = elInfoOld->fillFlag_;
262

263
264
265
266
267
    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);
268

269
    macroElement_  = elInfoOld->macroElement_;
270
    fillFlag_ = fill_flag;
271
    parent_ = elem;
272
    level = elInfoOld->level + 1;
273
    iChild = ichild;
274
275
276

    if (fillFlag_.isSet(Mesh::FILL_COORDS) || 
	fillFlag_.isSet(Mesh::FILL_DET)    ||
277
278
	fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
      
279
      if (elem->isNewCoordSet())
280
	coord_[2] = *(elem->getNewCoord());
281
282
      else
	coord_[2].setMidpoint(elInfoOld->coord_[0], elInfoOld->coord_[1]);      
283
284
      
      if (ichild == 0) {
285
286
	coord_[0] = elInfoOld->coord_[2];
	coord_[1] = elInfoOld->coord_[0];
287
      } else {
288
289
	coord_[0] = elInfoOld->coord_[1];
	coord_[1] = elInfoOld->coord_[2];
290
291
292
293
294
295
296
297
298
299
      }
    }

    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.

300
	neighbour_[2] = elInfoOld->neighbour_[1];
301
	oppVertex[2] = elInfoOld->oppVertex[1];
302
303
	
	if (neighbour_[2] && fill_opp_coords) {
304
305
	  oppCoord_[2] = elInfoOld->oppCoord_[1];
	  neighbourCoord_[2] = elInfoOld->neighbourCoord_[1];
306
	}
307
308
309
310
311
312
313
314
315
316
	
	
	// 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();
317
	  oppVertex[1] = 2;
318
319
	  
	  if (fill_opp_coords) {
320
            if (elem->getSecondChild()->isNewCoordSet()) {
321
322
	      oppCoord_[1] = *(elem->getSecondChild()->getNewCoord());
	    } else {      
323
324
	      oppCoord_[1].setMidpoint(elInfoOld->coord_[1], 
				       elInfoOld->coord_[2]);
325
	    }
326

327
328
329
	    neighbourCoord_[1][0] = coord_[0];
	    neighbourCoord_[1][1] = coord_[2];
	    neighbourCoord_[1][2] = oppCoord_[1];  
330
331
	  }
	} else {
332
	  neighbour_[1] = elem->getSecondChild();
333
	  oppVertex[1] = 0;
334
335

	  if (fill_opp_coords) {
336
	    oppCoord_[1] = elInfoOld->coord_[1];
337

338
339
	    neighbourCoord_[1][0] = elInfoOld->coord_[1];
	    neighbourCoord_[1][1] = elInfoOld->coord_[2];
340
	    neighbourCoord_[1][2] = coord_[2];
341
342
343
344
	  }
	}


345
346
347
	// Calculation of the neighbour 0, its oppCoords and the
	// cooresponding oppVertex.
	
348
	nb = elInfoOld->neighbour_[2];
349
	if (nb) {
350
	  TEST(elInfoOld->oppVertex[2] == 2)("invalid neighbour\n"); 
351
352
	  TEST_EXIT_DBG(nb->getFirstChild())("missing first child?\n");
	  TEST_EXIT_DBG(nb->getSecondChild())("missing second child?\n");
353
354
355
356
	 
	  nb = nb->getSecondChild();

	  if (nb->getFirstChild()) {
357
	    oppVertex[0] = 2;
358

359
	    if (fill_opp_coords) {
360
	      if (nb->isNewCoordSet()) {
361
362
		oppCoord_[0] = *(nb->getNewCoord());
	      } else {
363
364
		oppCoord_[0].setMidpoint(elInfoOld->neighbourCoord_[2][1],
					 elInfoOld->neighbourCoord_[2][2]);
365
	      }
366

367
368
369
	      neighbourCoord_[0][0].setMidpoint(elInfoOld->neighbourCoord_[2][0],
						elInfoOld->neighbourCoord_[2][1]);
	      neighbourCoord_[0][1] = elInfoOld->neighbourCoord_[2][1];
370
371
372
373
374
	      neighbourCoord_[0][2] = oppCoord_[0];
	    }	   
 
	    nb = nb->getFirstChild();
	  } else {
375
	    oppVertex[0] = 1;
376

377
	    if (fill_opp_coords) {
378
	      oppCoord_[0] = elInfoOld->oppCoord_[2];    
379

380
381
382
383
	      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]);
384
385
	    }
	  }
386
387
388
389
390
391
392
	}
	
	neighbour_[0] = nb;
      } else {   /* ichild == 1 */
	// Calculation of the neighbour 2, its oppCoords and the
	// cooresponding oppVertex.

393
	neighbour_[2] = elInfoOld->neighbour_[0];
394
	oppVertex[2] = elInfoOld->oppVertex[0];
395
396

	if (neighbour_[2] && fill_opp_coords) {
397
398
	  oppCoord_[2] = elInfoOld->oppCoord_[0];
	  neighbourCoord_[2] = elInfoOld->neighbourCoord_[0];
399
400
401
402
403
404
405
406
	}
	

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

	if (elem->getFirstChild()->getFirstChild()) {
	  neighbour_[0] = elem->getFirstChild()->getFirstChild();
407
	  oppVertex[0] = 2;
408
409

	  if (fill_opp_coords) {
410
            if (elem->getFirstChild()->isNewCoordSet()) {
411
	      oppCoord_[0] = *(elem->getFirstChild()->getNewCoord());
412
	    } else {
413
414
	      oppCoord_[0].setMidpoint(elInfoOld->coord_[0], 
				       elInfoOld->coord_[2]);
415
	    }
416
417
418
419

	    neighbourCoord_[0][0] = coord_[2];
	    neighbourCoord_[0][1] = coord_[1];
	    neighbourCoord_[0][2] = oppCoord_[0];
420
	  }
421
422
	} else {
	  neighbour_[0] = elem->getFirstChild();
423
	  oppVertex[0] = 1;
424
425

	  if (fill_opp_coords) {
426
	    oppCoord_[0] = elInfoOld->coord_[0];
427

428
429
	    neighbourCoord_[0][0] = elInfoOld->coord_[2];
	    neighbourCoord_[0][1] = elInfoOld->coord_[0];
430
	    neighbourCoord_[0][2] = coord_[2];
431
	  }
432
433
434
435
436
	}

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

437
	nb = elInfoOld->neighbour_[2];
438
	if (nb) {
439
	  TEST(elInfoOld->oppVertex[2] == 2)("invalid neighbour\n"); 
440
441
442
	  TEST((nb = nb->getFirstChild()))("missing child?\n");

	  if (nb->getFirstChild()) {
443
	    oppVertex[1] = 2;
444

445
	    if (fill_opp_coords) {
446
	      if (nb->isNewCoordSet()) {
447
		oppCoord_[1] = *(nb->getNewCoord());
448
	      } else {
449
450
		oppCoord_[1].setMidpoint(elInfoOld->neighbourCoord_[2][0],
					 elInfoOld->neighbourCoord_[2][2]);
451
	      }
452

453
454
455
	      neighbourCoord_[1][0] = elInfoOld->neighbourCoord_[2][0];
	      neighbourCoord_[1][1].setMidpoint(elInfoOld->neighbourCoord_[2][0],
						elInfoOld->neighbourCoord_[2][1]);
456
	      neighbourCoord_[1][2] = oppCoord_[1];
457
	    }
458
459
	    nb = nb->getSecondChild();

460
	  } else {
461
	    oppVertex[1] = 0;
462

463
	    if (fill_opp_coords) {
464
	      oppCoord_[1] = elInfoOld->oppCoord_[2];
465

466
467
468
469
	      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]);
470
471
472
	    }
	  }
	}
473
474
475
476
477
	neighbour_[1] = nb;
      } // if (ichild == 0) {} else
    } // if (fill_flag.isSet(Mesh::FILL_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS))
    

478
    if (fill_flag.isSet(Mesh::FILL_BOUND)) {
479
      if (elInfoOld->getBoundary(2))
480
	boundary_[5] = elInfoOld->getBoundary(2);
481
      else
482
483
	boundary_[5] = INTERIOR;

484
      if (ichild == 0) {
485
486
487
	boundary_[3] = elInfoOld->getBoundary(5);
	boundary_[4] = elInfoOld->getBoundary(3);
	boundary_[0] = elInfoOld->getBoundary(2);
488
	boundary_[1] = INTERIOR;
489
	boundary_[2] = elInfoOld->getBoundary(1);
490
      } else {
491
492
	boundary_[3] = elInfoOld->getBoundary(4);
	boundary_[4] = elInfoOld->getBoundary(5);
493
	boundary_[0] = INTERIOR;
494
495
	boundary_[1] = elInfoOld->getBoundary(2);
	boundary_[2] = elInfoOld->getBoundary(0);
496
497
      }

498
499
      if (elInfoOld->getProjection(0) && 
	  elInfoOld->getProjection(0)->getType() == VOLUME_PROJECTION) {
500
	
501
	projection_[0] = elInfoOld->getProjection(0);
502
503
      } else { // boundary projection
	if (ichild == 0) {
504
	  projection_[0] = elInfoOld->getProjection(2);
505
	  projection_[1] = NULL;
506
	  projection_[2] = elInfoOld->getProjection(1);
507
508
	} else {
	  projection_[0] = NULL;
509
510
	  projection_[1] = elInfoOld->getProjection(2);
	  projection_[2] = elInfoOld->getProjection(0);
511
512
513
514
515
	}
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
516
  double ElInfo2d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
517
  {
518
    FUNCNAME("ElInfo2d::calcGrdLambda()");
519
520

    testFlag(Mesh::FILL_COORDS);
521
522
  
    double adet = 0.0;
523
524
    int dim = mesh_->getDim();

Thomas Witkowski's avatar
Thomas Witkowski committed
525
    for (int i = 0; i < dimOfWorld; i++) {
526
527
      (*e1)[i] = coord_[1][i] - coord_[0][i];
      (*e2)[i] = coord_[2][i] - coord_[0][i];
528
529
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
530
    if (dimOfWorld == 2) {
531
      double sdet = (*e1)[0] * (*e2)[1] - (*e1)[1] * (*e2)[0];
532
533
534
535
      adet = abs(sdet);

      if (adet < 1.0E-25) {
	MSG("abs(det) = %f\n", adet);
536
	for (int i = 0; i <= dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
537
	  for (int j = 0; j < dimOfWorld; j++)
538
	    grd_lam[i][j] = 0.0;
539
540
      } else {
	double det1 = 1.0 / sdet;
541
542
543
544
545

	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
546
547
	grd_lam[0][0] = - grd_lam[1][0] - grd_lam[2][0];
	grd_lam[0][1] = - grd_lam[1][1] - grd_lam[2][1];
548
      }
549
550
    } else {  
      vectorProduct(*e1, *e2, *normal);
551

552
      adet = norm(normal);
553
554
555

      if (adet < 1.0E-15) {
	MSG("abs(det) = %lf\n", adet);
556
	for (int i = 0; i <= dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
557
	  for (int j = 0; j < dimOfWorld; j++)
558
559
	    grd_lam[i][j] = 0.0;
      } else {
560
561
	vectorProduct(*e2, *normal, grd_lam[1]);
	vectorProduct(*normal, *e1, grd_lam[2]);
562
      
563
	double adet2 = 1.0 / (adet * adet);
564

Thomas Witkowski's avatar
Thomas Witkowski committed
565
	for (int i = 0; i < dimOfWorld; i++) {
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
	  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
  {
582
    FUNCNAME("ElInfo::worldToCoord()");
583

584
    TEST_EXIT_DBG(lambda)("lambda must not be NULL\n");
585

586
587
588
589
    DimVec<WorldVector<double> > edge(mesh_->getDim(), NO_INIT);
    WorldVector<double> x; 
    static DimVec<double> vec(mesh_->getDim(), NO_INIT);

590
591
    int dim = mesh_->getDim();

592
593
    for (int j = 0; j < dimOfWorld; j++) {
      double x0 = coord_[dim][j];
594
      x[j] = xy[j] - x0;
595
      for (int i = 0; i < dim; i++)
596
597
598
	edge[i][j] = coord_[i][j] - x0;
    }
  
599
600
601
    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]; 
602
603
604

    if (abs(det) < DBL_TOL) {
      ERROR("det = %le; abort\n", det);
605
606
      for (int i = 0; i <= dim; i++) 
	(*lambda)[i] = 1.0/dim;
607
608
609
610
611
612
613
614
615
      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;
616
    for (int i = 0; i <= dim; i++) {
617
618
619
620
621
622
623
624
625
626
627
628
      if ((*lambda)[i] < -1.E-5) {
	if ((*lambda)[i] < lmin) {
	  k = i;
	  lmin = (*lambda)[i];
	}
      }
    }

    return k;
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
629
  double ElInfo2d::getNormal(int side, WorldVector<double> &normal)
630
  {
631
    FUNCNAME("ElInfo::getNormal()");
632

633
634
    int i0 = (side + 1) % 3;
    int i1 = (side + 2) % 3;
635

Thomas Witkowski's avatar
Thomas Witkowski committed
636
    if (dimOfWorld == 2){
637
638
639
640
641
      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;

642
643
644
645
646
647
      e0 = coord_[i1]; 
      e0 -= coord_[i0];
      e1 = coord_[i1]; 
      e1 -= coord_[side];
      e2 = coord_[i0]; 
      e2 -= coord_[side];
648
649
650
651
652

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

653
    double det = norm(&normal);
654

655
    TEST_EXIT_DBG(det > 1.e-30)("det = 0 on face %d\n", side);
656
657
658

    normal *= 1.0 / det;
    
Thomas Witkowski's avatar
Thomas Witkowski committed
659
    return det;
660
661
662
663
664
665
666
  }

  /****************************************************************************/
  /*  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
667
  double ElInfo2d::getElementNormal(WorldVector<double> &elementNormal) const
668
  {
669
    FUNCNAME("ElInfo::getElementNormal()");
670

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

674
675
    WorldVector<double> e0 = coord_[1] - coord_[0];
    WorldVector<double> e1 = coord_[2] - coord_[0];
676
677
678

    vectorProduct(e0, e1, elementNormal);

679
    double det = norm(&elementNormal);
680

681
    TEST_EXIT_DBG(det > 1.e-30)("det = 0");
682
683
684

    elementNormal *= 1.0 / det;
    
Thomas Witkowski's avatar
Thomas Witkowski committed
685
    return det;
686
  }
687

688
}