ElInfo2d.cc 15.7 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
	coord_[i] = mel->coord[i];
      }
    }

    int neighbours = mesh_->getGeo(NEIGH);

38
39
40
41
    if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS) || 
	fillFlag_.isSet(Mesh::FILL_NEIGH)) {

      bool fill_opp_coords = (fillFlag_.isSet(Mesh::FILL_OPP_COORDS));
42
    
43
44
45
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
      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];
76
		} else {
77
		  ERROR_EXIT("should not happen!\n");
78
		}
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
	
		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;
100
101
102
	      }
	    }
	  } else {
103
104
105
106
107
	    if (fill_opp_coords) {
	      oppCoord_[i] = macroNeighbour->coord[edgeNo];

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

120
      for (int i = 0; i < element_->getGeo(PROJECTION); i++) {
121
122
123
124
125
126
127
128
129
130
131
132
	projection_[i] = mel->getProjection(i);
      }
    }
  }


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

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

135
136
137
138
    Element *elem = elinfo_old->element_;
    Element *nb;

    Flag fill_flag = elinfo_old->fillFlag_;
139

140
141
142
143
144
    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);
145
146
147

    macroElement_  = elinfo_old->macroElement_;
    fillFlag_ = fill_flag;
148
149
    parent_ = elem;
    level_ = elinfo_old->level_ + 1;
150
151
152
153
154

    int dow = Global::getGeo(WORLD);

    if (fillFlag_.isSet(Mesh::FILL_COORDS) || 
	fillFlag_.isSet(Mesh::FILL_DET)    ||
155
156
157
158
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
	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];
185
	}
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
	
	
	// 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]);
	    }
205

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


224
225
226
227
228
229
	// 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"); 
230
231
	  TEST_EXIT_DBG(nb->getFirstChild())("missing first child?\n");
	  TEST_EXIT_DBG(nb->getSecondChild())("missing second child?\n");
232
233
234
235
236
237
	 
	  nb = nb->getSecondChild();

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

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

	      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;

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

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

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

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

	      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];
336
	    }
337
338
	    nb = nb->getSecondChild();

339
	  } else {
340
341
	    oppVertex_[1] = 0;

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

357
    if (fill_flag.isSet(Mesh::FILL_BOUND)) {
358
      if (elinfo_old->getBoundary(2)) {
359
	boundary_[5] = elinfo_old->getBoundary(2);
360
      } else {
361
	boundary_[5] = INTERIOR;
362
      }
363

364
      if (ichild == 0) {
365
366
367
368
369
370
371
372
373
374
375
376
377
	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);
      }

378
379
380
381
382
383
      if (elinfo_old->getProjection(0) && 
	  elinfo_old->getProjection(0)->getType() == VOLUME_PROJECTION) {
	
	projection_[0] = elinfo_old->getProjection(0);
      } else { // boundary projection
	if (ichild == 0) {
384
385
386
387
388
389
390
391
392
393
394
395
396
397
	  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
  {
398
    FUNCNAME("ElInfo2d::calcGrdLambda()");
399
400
  
    WorldVector<double> e1, e2;
401
402
    double adet = 0.0;
    const WorldVector<double> v0 = coord_[0];
403
404
405
406
407

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

408
    for (int i = 0; i < dow; i++) {
409
410
411
412
      e1[i] = coord_[1][i] - v0[i];
      e2[i] = coord_[2][i] - v0[i];
    }

413
414
    if (dow == 2) {
      double sdet = e1[0] * e2[1] - e1[1] * e2[0];
415
416
417
418
      adet = abs(sdet);

      if (adet < 1.0E-25) {
	MSG("abs(det) = %f\n", adet);
419
420
	for (int i = 0; i <= dim; i++)
	  for (int j = 0; j < dow; j++)
421
	    grd_lam[i][j] = 0.0;
422
423
424
425
426
427
428
429
430
431
432
433
434
      } 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];
435
436
437
438
439
440
441
442
443
444
      }
    } else {
      WorldVector<double> normal;
    
      vectorProduct(e1, e2, normal);

      adet = norm(&normal);

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

454
	for (int i = 0; i < dow; i++) {
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
	  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
  {
471
    FUNCNAME("ElInfo::worldToCoord()");
472

473
    TEST_EXIT_DBG(lambda)("lambda must not be NULL\n");
474

475
476
477
478
    DimVec<WorldVector<double> > edge(mesh_->getDim(), NO_INIT);
    WorldVector<double> x; 
    static DimVec<double> vec(mesh_->getDim(), NO_INIT);

479
480
481
    int dim = mesh_->getDim();
    int dimOfWorld = Global::getGeo(WORLD);

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

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

523
524
    int i0 = (side + 1) % 3;
    int i1 = (side + 2) % 3;
525

526
    if (Global::getGeo(WORLD) == 2){
527
528
529
530
531
      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;

532
533
534
535
536
537
      e0 = coord_[i1]; 
      e0 -= coord_[i0];
      e1 = coord_[i1]; 
      e1 -= coord_[side];
      e2 = coord_[i0]; 
      e2 -= coord_[side];
538
539
540
541
542

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

543
    double det = norm(&normal);
544

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

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

562
    TEST_EXIT_DBG(Global::getGeo(WORLD) == 3)(" element normal only well defined for  DIM_OF_WORLD = DIM + 1 !!");
563

564
565
    WorldVector<double> e0 = coord_[1] - coord_[0];
    WorldVector<double> e1 = coord_[2] - coord_[0];
566
567
568

    vectorProduct(e0, e1, elementNormal);

569
    double det = norm(&elementNormal);
570

571
    TEST_EXIT_DBG(det > 1.e-30)("det = 0");
572
573
574
575
576
577

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