ElInfo2d.cc 15.6 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

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

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


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

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

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

	      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;

254
	    if (fill_opp_coords) {
255
256
257
258
259
260
	      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]);
261
262
	    }
	  }
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
	}
	
	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());
289
	    } else {
290
291
	      oppCoord_[0].setMidpoint(elinfo_old->coord_[0], 
				       elinfo_old->coord_[2]);
292
	    }
293
294
295
296

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

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

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

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

337
	  } else {
338
339
	    oppVertex_[1] = 0;

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

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

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

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

    testFlag(Mesh::FILL_COORDS);
    int dim = mesh_->getDim();

Thomas Witkowski's avatar
Thomas Witkowski committed
405
    for (int i = 0; i < dimOfWorld; i++) {
406
407
408
409
      e1[i] = coord_[1][i] - v0[i];
      e2[i] = coord_[2][i] - v0[i];
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
410
    if (dimOfWorld == 2) {
411
      double sdet = e1[0] * e2[1] - e1[1] * e2[0];
412
413
414
415
      adet = abs(sdet);

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

      adet = norm(&normal);

      if (adet < 1.0E-15) {
	MSG("abs(det) = %lf\n", adet);
442
	for (int i = 0; i <= dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
443
	  for (int j = 0; j < dimOfWorld; j++)
444
445
446
447
448
	    grd_lam[i][j] = 0.0;
      } else {
	vectorProduct(e2, normal, grd_lam[1]);
	vectorProduct(normal, e1, grd_lam[2]);
      
449
	double adet2 = 1.0 / (adet * adet);
450

Thomas Witkowski's avatar
Thomas Witkowski committed
451
	for (int i = 0; i < dimOfWorld; i++) {
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
	  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
  {
468
    FUNCNAME("ElInfo::worldToCoord()");
469

470
    TEST_EXIT_DBG(lambda)("lambda must not be NULL\n");
471

472
473
474
475
    DimVec<WorldVector<double> > edge(mesh_->getDim(), NO_INIT);
    WorldVector<double> x; 
    static DimVec<double> vec(mesh_->getDim(), NO_INIT);

476
477
    int dim = mesh_->getDim();

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

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

519
520
    int i0 = (side + 1) % 3;
    int i1 = (side + 2) % 3;
521

Thomas Witkowski's avatar
Thomas Witkowski committed
522
    if (dimOfWorld == 2){
523
524
525
526
527
      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;

528
529
530
531
532
533
      e0 = coord_[i1]; 
      e0 -= coord_[i0];
      e1 = coord_[i1]; 
      e1 -= coord_[side];
      e2 = coord_[i0]; 
      e2 -= coord_[side];
534
535
536
537
538

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

539
    double det = norm(&normal);
540

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

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

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

561
562
    WorldVector<double> e0 = coord_[1] - coord_[0];
    WorldVector<double> e1 = coord_[2] - coord_[0];
563
564
565

    vectorProduct(e0, e1, elementNormal);

566
    double det = norm(&elementNormal);
567

568
    TEST_EXIT_DBG(det > 1.e-30)("det = 0");
569
570
571
572
573
574

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