ElInfo2d.cc 15.8 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
62
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
90
      for (int i = 0; i < neighbours; i++) {
	MacroElement *macroNeighbour = mel->getNeighbour(i);

	if (macroNeighbour) {
	  neighbour_[i] = macroNeighbour->getElement();
	  
	  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];
91
		} else {
92
		  ERROR_EXIT("should not happen!\n");
93
		}
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
	
		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:
		ERROR_EXIT("should not happen!\n");
		break;
115
116
117
	      }
	    }
	  } else {
118
119
120
121
122
	    if (fill_opp_coords) {
	      oppCoord_[i] = macroNeighbour->coord[edgeNo];

	      neighbourCoord_[i] = macroNeighbour->coord;	      
	    }
123
	  }
124
125
126
	} else {
	  neighbour_[i] = NULL;
        }
127
      }
128
    }
129
    
130
131
    if (fillFlag_.isSet(Mesh::FILL_BOUND)) {   
      for (int i = 0; i < element_->getGeo(BOUNDARY); i++) {
132
133
134
	boundary_[i] = mel->getBoundary(i);
      }

135
      for (int i = 0; i < element_->getGeo(PROJECTION); i++) {
136
137
138
139
140
141
142
143
144
145
146
147
	projection_[i] = mel->getProjection(i);
      }
    }
  }


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

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

150
151
152
153
    Element *elem = elinfo_old->element_;
    Element *nb;

    Flag fill_flag = elinfo_old->fillFlag_;
154

155
156
157
158
159
    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);
160
161
162

    macroElement_  = elinfo_old->macroElement_;
    fillFlag_ = fill_flag;
163
164
    parent_ = elem;
    level_ = elinfo_old->level_ + 1;
165
166
167

    if (fillFlag_.isSet(Mesh::FILL_COORDS) || 
	fillFlag_.isSet(Mesh::FILL_DET)    ||
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
	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];
198
	}
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
	
	
	// 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]);
	    }
218

219
220
221
	    neighbourCoord_[1][0] = coord_[0];
	    neighbourCoord_[1][1] = coord_[2];
	    neighbourCoord_[1][2] = oppCoord_[1];  
222
223
	  }
	} else {
224
225
226
227
228
229
230
231
232
	  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];
233
234
235
236
	  }
	}


237
238
239
240
241
242
	// 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"); 
243
244
	  TEST_EXIT_DBG(nb->getFirstChild())("missing first child?\n");
	  TEST_EXIT_DBG(nb->getSecondChild())("missing second child?\n");
245
246
247
248
249
250
	 
	  nb = nb->getSecondChild();

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

251
	    if (fill_opp_coords) {
252
253
254
255
256
	      if (nb->getNewCoord(-1)) {
		oppCoord_[0] = *(nb->getNewCoord());
	      } else {
		oppCoord_[0].setMidpoint(elinfo_old->neighbourCoord_[2][1],
					 elinfo_old->neighbourCoord_[2][2]);
257
	      }
258
259
260
261
262
263
264
265
266
267
268

	      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;

269
	    if (fill_opp_coords) {
270
271
272
273
274
275
	      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]);
276
277
	    }
	  }
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
	}
	
	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());
304
	    } else {
305
306
	      oppCoord_[0].setMidpoint(elinfo_old->coord_[0], 
				       elinfo_old->coord_[2]);
307
	    }
308
309
310
311

	    neighbourCoord_[0][0] = coord_[2];
	    neighbourCoord_[0][1] = coord_[1];
	    neighbourCoord_[0][2] = oppCoord_[0];
312
	  }
313
314
315
316
317
318
319
320
321
322
	} 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];
323
	  }
324
325
326
327
328
329
330
331
332
333
334
335
336
	}

	// 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;

337
	    if (fill_opp_coords) {
338
339
	      if (nb->getNewCoord(-1)) {
		oppCoord_[1] = *(nb->getNewCoord());
340
	      } else {
341
342
		oppCoord_[1].setMidpoint(elinfo_old->neighbourCoord_[2][0],
					 elinfo_old->neighbourCoord_[2][2]);
343
	      }
344
345
346
347
348

	      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];
349
	    }
350
351
	    nb = nb->getSecondChild();

352
	  } else {
353
354
	    oppVertex_[1] = 0;

355
	    if (fill_opp_coords) {
356
357
358
359
360
361
	      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]);
362
363
364
	    }
	  }
	}
365
366
367
368
369
	neighbour_[1] = nb;
      } // if (ichild == 0) {} else
    } // if (fill_flag.isSet(Mesh::FILL_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS))
    

370
    if (fill_flag.isSet(Mesh::FILL_BOUND)) {
371
      if (elinfo_old->getBoundary(2)) {
372
	boundary_[5] = elinfo_old->getBoundary(2);
373
      } else {
374
	boundary_[5] = INTERIOR;
375
      }
376

377
      if (ichild == 0) {
378
379
380
381
382
383
384
385
386
387
388
389
390
	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);
      }

391
392
393
394
395
396
      if (elinfo_old->getProjection(0) && 
	  elinfo_old->getProjection(0)->getType() == VOLUME_PROJECTION) {
	
	projection_[0] = elinfo_old->getProjection(0);
      } else { // boundary projection
	if (ichild == 0) {
397
398
399
400
401
402
403
404
405
406
407
408
	  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
409
  double ElInfo2d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
410
  {
411
    FUNCNAME("ElInfo2d::calcGrdLambda()");
412
413

    testFlag(Mesh::FILL_COORDS);
414
415
  
    double adet = 0.0;
416
417
    int dim = mesh_->getDim();

Thomas Witkowski's avatar
Thomas Witkowski committed
418
    for (int i = 0; i < dimOfWorld; i++) {
419
420
      (*e1)[i] = coord_[1][i] - coord_[0][i];
      (*e2)[i] = coord_[2][i] - coord_[0][i];
421
422
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
423
    if (dimOfWorld == 2) {
424
      double sdet = (*e1)[0] * (*e2)[1] - (*e1)[1] * (*e2)[0];
425
426
427
428
      adet = abs(sdet);

      if (adet < 1.0E-25) {
	MSG("abs(det) = %f\n", adet);
429
	for (int i = 0; i <= dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
430
	  for (int j = 0; j < dimOfWorld; j++)
431
	    grd_lam[i][j] = 0.0;
432
433
      } else {
	double det1 = 1.0 / sdet;
434
435
436
437
438

	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
439
440
	grd_lam[0][0] = - grd_lam[1][0] - grd_lam[2][0];
	grd_lam[0][1] = - grd_lam[1][1] - grd_lam[2][1];
441
      }
442
443
    } else {  
      vectorProduct(*e1, *e2, *normal);
444

445
      adet = norm(normal);
446
447
448

      if (adet < 1.0E-15) {
	MSG("abs(det) = %lf\n", adet);
449
	for (int i = 0; i <= dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
450
	  for (int j = 0; j < dimOfWorld; j++)
451
452
	    grd_lam[i][j] = 0.0;
      } else {
453
454
	vectorProduct(*e2, *normal, grd_lam[1]);
	vectorProduct(*normal, *e1, grd_lam[2]);
455
      
456
	double adet2 = 1.0 / (adet * adet);
457

Thomas Witkowski's avatar
Thomas Witkowski committed
458
	for (int i = 0; i < dimOfWorld; i++) {
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
	  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
  {
475
    FUNCNAME("ElInfo::worldToCoord()");
476

477
    TEST_EXIT_DBG(lambda)("lambda must not be NULL\n");
478

479
480
481
482
    DimVec<WorldVector<double> > edge(mesh_->getDim(), NO_INIT);
    WorldVector<double> x; 
    static DimVec<double> vec(mesh_->getDim(), NO_INIT);

483
484
    int dim = mesh_->getDim();

485
486
    for (int j = 0; j < dimOfWorld; j++) {
      double x0 = coord_[dim][j];
487
      x[j] = xy[j] - x0;
488
      for (int i = 0; i < dim; i++)
489
490
491
	edge[i][j] = coord_[i][j] - x0;
    }
  
492
493
494
    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]; 
495
496
497

    if (abs(det) < DBL_TOL) {
      ERROR("det = %le; abort\n", det);
498
499
      for (int i = 0; i <= dim; i++) 
	(*lambda)[i] = 1.0/dim;
500
501
502
503
504
505
506
507
508
      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;
509
    for (int i = 0; i <= dim; i++) {
510
511
512
513
514
515
516
517
518
519
520
521
      if ((*lambda)[i] < -1.E-5) {
	if ((*lambda)[i] < lmin) {
	  k = i;
	  lmin = (*lambda)[i];
	}
      }
    }

    return k;
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
522
  double ElInfo2d::getNormal(int side, WorldVector<double> &normal)
523
  {
524
    FUNCNAME("ElInfo::getNormal()");
525

526
527
    int i0 = (side + 1) % 3;
    int i1 = (side + 2) % 3;
528

Thomas Witkowski's avatar
Thomas Witkowski committed
529
    if (dimOfWorld == 2){
530
531
532
533
534
      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;

535
536
537
538
539
540
      e0 = coord_[i1]; 
      e0 -= coord_[i0];
      e1 = coord_[i1]; 
      e1 -= coord_[side];
      e2 = coord_[i0]; 
      e2 -= coord_[side];
541
542
543
544
545

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

546
    double det = norm(&normal);
547

548
    TEST_EXIT_DBG(det > 1.e-30)("det = 0 on face %d\n", side);
549
550
551
552
553
554
555
556
557
558
559
560
561
562

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

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

568
569
    WorldVector<double> e0 = coord_[1] - coord_[0];
    WorldVector<double> e1 = coord_[2] - coord_[0];
570
571
572

    vectorProduct(e0, e1, elementNormal);

573
    double det = norm(&elementNormal);
574

575
    TEST_EXIT_DBG(det > 1.e-30)("det = 0");
576
577
578
579
580
581

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