ElInfo2d.cc 15.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#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 {

  void ElInfo2d::fillMacroInfo(const MacroElement * mel)
  {
19
20
    FUNCNAME("ElInfo::fillMacroInfo()");
 
21
    macroElement_ = const_cast<MacroElement*>(mel);
22
23
24
    element_  = const_cast<Element*>(mel->getElement());
    parent_  = NULL;
    level_ = 0;
25
26
27
28

    if (fillFlag_.isSet(Mesh::FILL_COORDS) || 
	fillFlag_.isSet(Mesh::FILL_DET)    ||
	fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
29
30
31

      int vertices = mesh_->getGeo(VERTEX);
      for (int i = 0; i < vertices; i++) {
32
33
34
35
36
37
38
39
40
41
	coord_[i] = mel->coord[i];
      }
    }

    //   if(fillFlag_.isSet(Mesh::FILL_DET) || fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
    //     det = calcGrdLambda(*Lambda);
    //   }

    int neighbours = mesh_->getGeo(NEIGH);

42
43
44
45
    if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS) || 
	fillFlag_.isSet(Mesh::FILL_NEIGH)) {

      bool fill_opp_coords = (fillFlag_.isSet(Mesh::FILL_OPP_COORDS));
46
    
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
      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];
80
		} else {
81
		  ERROR_EXIT("should not happen!\n");
82
		}
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
	
		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;
104
105
106
	      }
	    }
	  } else {
107
108
109
110
111
	    if (fill_opp_coords) {
	      oppCoord_[i] = macroNeighbour->coord[edgeNo];

	      neighbourCoord_[i] = macroNeighbour->coord;	      
	    }
112
	  }
113
114
115
	} else {
	  neighbour_[i] = NULL;
        }
116
      }
117
    }
118
    
119
120
    if (fillFlag_.isSet(Mesh::FILL_BOUND)) {   
      for (int i = 0; i < element_->getGeo(BOUNDARY); i++) {
121
122
123
	boundary_[i] = mel->getBoundary(i);
      }

124
      for (int i = 0; i < element_->getGeo(PROJECTION); i++) {
125
126
127
128
129
130
131
132
133
134
135
136
	projection_[i] = mel->getProjection(i);
      }
    }
  }


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

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

139
140
141
142
    Element *elem = elinfo_old->element_;
    Element *nb;

    Flag fill_flag = elinfo_old->fillFlag_;
143
144

    TEST_EXIT(elem->getFirstChild())("no children?\n");
145
    TEST_EXIT(element_ = const_cast<Element*>((ichild == 0) ? 
146
147
148
149
150
151
					      elem->getFirstChild() : 
					      elem->getSecondChild()))
      ("missing child %d?\n", ichild);

    macroElement_  = elinfo_old->macroElement_;
    fillFlag_ = fill_flag;
152
153
    parent_ = elem;
    level_ = elinfo_old->level_ + 1;
154
155
156
157
158

    int dow = Global::getGeo(WORLD);

    if (fillFlag_.isSet(Mesh::FILL_COORDS) || 
	fillFlag_.isSet(Mesh::FILL_DET)    ||
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
	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];
189
	}
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
	
	
	// 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]);
	    }
209

210
211
212
	    neighbourCoord_[1][0] = coord_[0];
	    neighbourCoord_[1][1] = coord_[2];
	    neighbourCoord_[1][2] = oppCoord_[1];  
213
214
	  }
	} else {
215
216
217
218
219
220
221
222
223
	  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];
224
225
226
227
	  }
	}


228
229
230
231
232
233
234
235
236
237
238
239
240
241
	// 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"); 
	  TEST_EXIT(nb->getFirstChild())("missing first child?\n");
	  TEST_EXIT(nb->getSecondChild())("missing second child?\n");
	 
	  nb = nb->getSecondChild();

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

242
	    if (fill_opp_coords) {
243
244
245
246
247
	      if (nb->getNewCoord(-1)) {
		oppCoord_[0] = *(nb->getNewCoord());
	      } else {
		oppCoord_[0].setMidpoint(elinfo_old->neighbourCoord_[2][1],
					 elinfo_old->neighbourCoord_[2][2]);
248
	      }
249
250
251
252
253
254
255
256
257
258
259

	      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;

260
	    if (fill_opp_coords) {
261
262
263
264
265
266
	      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]);
267
268
	    }
	  }
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
	}
	
	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());
295
	    } else {
296
297
	      oppCoord_[0].setMidpoint(elinfo_old->coord_[0], 
				       elinfo_old->coord_[2]);
298
	    }
299
300
301
302

	    neighbourCoord_[0][0] = coord_[2];
	    neighbourCoord_[0][1] = coord_[1];
	    neighbourCoord_[0][2] = oppCoord_[0];
303
	  }
304
305
306
307
308
309
310
311
312
313
	} 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];
314
	  }
315
316
317
318
319
320
321
322
323
324
325
326
327
	}

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

328
	    if (fill_opp_coords) {
329
330
	      if (nb->getNewCoord(-1)) {
		oppCoord_[1] = *(nb->getNewCoord());
331
	      } else {
332
333
		oppCoord_[1].setMidpoint(elinfo_old->neighbourCoord_[2][0],
					 elinfo_old->neighbourCoord_[2][2]);
334
	      }
335
336
337
338
339

	      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];
340
	    }
341
342
	    nb = nb->getSecondChild();

343
	  } else {
344
345
	    oppVertex_[1] = 0;

346
	    if (fill_opp_coords) {
347
348
349
350
351
352
	      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]);
353
354
355
	    }
	  }
	}
356
357
358
359
360
	neighbour_[1] = nb;
      } // if (ichild == 0) {} else
    } // if (fill_flag.isSet(Mesh::FILL_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS))
    

361
    if (fill_flag.isSet(Mesh::FILL_BOUND)) {
362
      if (elinfo_old->getBoundary(2)) {
363
	boundary_[5] = elinfo_old->getBoundary(2);
364
      } else {
365
	boundary_[5] = INTERIOR;
366
      }
367

368
      if (ichild == 0) {
369
370
371
372
373
374
375
376
377
378
379
380
381
	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);
      }

382
383
384
385
386
387
      if (elinfo_old->getProjection(0) && 
	  elinfo_old->getProjection(0)->getType() == VOLUME_PROJECTION) {
	
	projection_[0] = elinfo_old->getProjection(0);
      } else { // boundary projection
	if (ichild == 0) {
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
	  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);
	}
      }
    }
  }

  double ElInfo2d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam) const
  {
    FUNCNAME("ElInfo2d::calcGrdLambda");
    //   TEST_EXIT(Global::getGeo(WORLD) == 2)
    //     ("dim != dim_of_world ! use parametric elements!\n");
  
    WorldVector<double> e1, e2;
407
408
    double adet = 0.0;
    const WorldVector<double> v0 = coord_[0];
409
410
411
412
413

    testFlag(Mesh::FILL_COORDS);
    int dow = Global::getGeo(WORLD);
    int dim = mesh_->getDim();

414
    for (int i = 0; i < dow; i++) {
415
416
417
418
      e1[i] = coord_[1][i] - v0[i];
      e2[i] = coord_[2][i] - v0[i];
    }

419
420
    if (dow == 2) {
      double sdet = e1[0] * e2[1] - e1[1] * e2[0];
421
422
423
424
      adet = abs(sdet);

      if (adet < 1.0E-25) {
	MSG("abs(det) = %f\n", adet);
425
426
	for (int i = 0; i <= dim; i++)
	  for (int j = 0; j < dow; j++)
427
	    grd_lam[i][j] = 0.0;
428
429
430
431
432
433
434
435
436
437
438
439
440
      } else {
	double det1 = 1.0 / sdet;
	double a11 =  e2[1] * det1;       /* (a_ij) = A^{-T} */
	double a21 = -e2[0] * det1;
	double a12 = -e1[1] * det1;
	double a22 =  e1[0] * det1;
	
	grd_lam[1][0] = a11;
	grd_lam[1][1] = a21;
	grd_lam[2][0] = a12;
	grd_lam[2][1] = a22;
	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
444
445
446
447
448
449
450
      }
    } else {
      WorldVector<double> normal;
    
      vectorProduct(e1, e2, normal);

      adet = norm(&normal);

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

460
	for (int i = 0; i < dow; i++) {
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
	  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
  {
477
    FUNCNAME("ElInfo::worldToCoord()");
478
479
480

    TEST_EXIT(lambda)("lambda must not be NULL\n");

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

485
486
487
    int dim = mesh_->getDim();
    int dimOfWorld = Global::getGeo(WORLD);

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

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

    return k;
  }


  double ElInfo2d::getNormal(int side, WorldVector<double> &normal) const
  {
527
    FUNCNAME("ElInfo::getNormal()");
528

529
530
    int i0 = (side + 1) % 3;
    int i1 = (side + 2) % 3;
531

532
    if (Global::getGeo(WORLD) == 2){
533
534
535
536
537
      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;

538
539
540
541
542
543
      e0 = coord_[i1]; 
      e0 -= coord_[i0];
      e1 = coord_[i1]; 
      e1 -= coord_[side];
      e2 = coord_[i0]; 
      e2 -= coord_[side];
544
545
546
547
548

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

549
    double det = norm(&normal);
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565

    TEST_EXIT(det > 1.e-30)("det = 0 on face %d\n", side);

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

568
    TEST_EXIT(Global::getGeo(WORLD) == 3)(" element normal only well defined for  DIM_OF_WORLD = DIM + 1 !!");
569

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

    vectorProduct(e0, e1, elementNormal);

575
    double det = norm(&elementNormal);
576
577
578
579
580
581
582
583

    TEST_EXIT(det > 1.e-30)("det = 0");

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