ElInfo2d.cc 20.1 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
  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);

22
  /*
23
  double ElInfo2d::mat_d1_left_val[3][3] = {{0.0, 1.0, 0.5}, 
24
25
26
27
28
29
30
31
  					    {0.0, 0.0, 0.5},
  					    {1.0, 0.0, 0.0}};
  */

  double ElInfo2d::mat_d1_left_val[3][3] = {{0.0, 0.0, 1.0}, 
  					    {1.0, 0.0, 0.0},
  					    {0.5, 0.5, 0.0}};

32
33
  mtl::dense2D<double> ElInfo2d::mat_d1_left(mat_d1_left_val);

34
  /*
35
  double ElInfo2d::mat_d1_right_val[3][3] = {{0.0, 0.0, 0.5}, 
36
37
38
39
40
41
42
43
  					     {1.0, 0.0, 0.5},
  					     {0.0, 1.0, 0.0}};
  */

  double ElInfo2d::mat_d1_right_val[3][3] = {{0.0, 1.0, 0.0}, 
  					     {0.0, 0.0, 1.0},
  					     {0.5, 0.5, 0.0}};

44
45
  mtl::dense2D<double> ElInfo2d::mat_d1_right(mat_d1_right_val);

46

47
48
49
  ElInfo2d::ElInfo2d(Mesh *aMesh) 
    : ElInfo(aMesh) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
50
51
52
    e1 = new WorldVector<double>;
    e2 = new WorldVector<double>;
    normal = new WorldVector<double>;
53
54
  }

55

56
57
  ElInfo2d::~ElInfo2d()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
58
59
60
    delete e1;
    delete e2;
    delete normal;
61
62
  }

63

64
65
  void ElInfo2d::fillMacroInfo(const MacroElement * mel)
  {
66
67
    FUNCNAME("ElInfo::fillMacroInfo()");
 
Thomas Witkowski's avatar
Thomas Witkowski committed
68
69
70
    macroElement = const_cast<MacroElement*>(mel);
    element = const_cast<Element*>(mel->getElement());
    parent = NULL;
71
    level = 0;
72

Thomas Witkowski's avatar
Thomas Witkowski committed
73
74
75
    if (fillFlag.isSet(Mesh::FILL_COORDS) || 
	fillFlag.isSet(Mesh::FILL_DET)    ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
76

Thomas Witkowski's avatar
Thomas Witkowski committed
77
      int vertices = mesh->getGeo(VERTEX);
78
      for (int i = 0; i < vertices; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
79
	coord[i] = mel->coord[i];
80
81
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
82
    int neighbours = mesh->getGeo(NEIGH);
83

Thomas Witkowski's avatar
Thomas Witkowski committed
84
85
    if (fillFlag.isSet(Mesh::FILL_OPP_COORDS) || 
	fillFlag.isSet(Mesh::FILL_NEIGH)) {
86

Thomas Witkowski's avatar
Thomas Witkowski committed
87
      bool fill_opp_coords = (fillFlag.isSet(Mesh::FILL_OPP_COORDS));
88
    
89
90
91
92
      for (int i = 0; i < neighbours; i++) {
	MacroElement *macroNeighbour = mel->getNeighbour(i);

	if (macroNeighbour) {
Thomas Witkowski's avatar
Thomas Witkowski committed
93
94
	  neighbour[i] = macroNeighbour->getElement();	  
	  Element *nb = const_cast<Element*>(neighbour[i]);
95

96
	  int edgeNo = oppVertex[i] = mel->getOppVertex(i);
97
98


99
	  if (nb->getFirstChild() && edgeNo != 2) {  
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119

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


120
	    if (edgeNo == 0) {
121
122
123
124
125
126
127
128
129
130
131
132
	      // The situation is as follows:
	      //
	      //          -------
	      //          \    /|\
	      //           \  / | \
	      //            \/  |  \
	      //             \  |   \
	      //              \ |    \ 
	      //               \|     \
	      //                -------
	      //
	      //            nb     el
133
              //
134
135
136
137
	      // 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.

Thomas Witkowski's avatar
Thomas Witkowski committed
138
	      nb = neighbour[i] = nb->getSecondChild();
139
	    } else {
140
141
	      // The situation is as shown in the picture above. So the next
	      // neighbour is the first child of the same level neighbour element.
Thomas Witkowski's avatar
Thomas Witkowski committed
142
	      nb = neighbour[i] = nb->getFirstChild();
143
144
	    }

145
146
	    // In both cases the opp vertex number is 2, as one can see in the 
	    // pictures above.
147
	    oppVertex[i] = 2;
148
149

	    if (fill_opp_coords) {
150
	      if (nb->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
151
		oppCoord[i] = *(nb->getNewCoord());
152
	      } else {
153
154
155
		// 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.
Thomas Witkowski's avatar
Thomas Witkowski committed
156
		oppCoord[i] = (macroNeighbour->coord[0] + 
157
158
				macroNeighbour->coord[1]) * 0.5;
	      }
159
160
161
	      
	      switch (i) {
	      case 0:
162
163
164
165
		// The common edge is the edge 0 of this element.

		switch (edgeNo) {
		case 1:
166
167
		  neighbourCoord[i][0] = macroNeighbour->coord[2];
		  neighbourCoord[i][1] = macroNeighbour->coord[0];
168
169
		  break;
		case 0:		  
170
171
		  neighbourCoord[i][0] = macroNeighbour->coord[1];
		  neighbourCoord[i][1] = macroNeighbour->coord[2];
172
173
		  break;
		default:
Thomas Witkowski's avatar
Thomas Witkowski committed
174
		  ERROR_EXIT("Should not happen!\n");
175
		}
176
	
Thomas Witkowski's avatar
Thomas Witkowski committed
177
		neighbourCoord[i][2] = oppCoord[i];
178
179
180
		break;
		
	      case 1:
181
182
183
		// The commonedge is the edge 1 of this element.
		switch (edgeNo) {
		case 0:
184
185
		  neighbourCoord[i][0] = macroNeighbour->coord[1];
		  neighbourCoord[i][1] = macroNeighbour->coord[2];
186
187
		  break;
		case 1:
188
189
		  neighbourCoord[i][0] = macroNeighbour->coord[2];
		  neighbourCoord[i][1] = macroNeighbour->coord[0];
190
191
		  break;
		default:
Thomas Witkowski's avatar
Thomas Witkowski committed
192
		  ERROR_EXIT("Should not happen!\n");
193
194
		}
		
Thomas Witkowski's avatar
Thomas Witkowski committed
195
		neighbourCoord[i][2] = oppCoord[i];
196
197
		break;
		
Thomas Witkowski's avatar
Thomas Witkowski committed
198
	      case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
199
		if (*(macroNeighbour->getElement()->getDOF(2)) == *(element->getDOF(0))) {
200
201
		  neighbourCoord[i][0] = macroNeighbour->coord[2];
		  neighbourCoord[i][1] = macroNeighbour->coord[1];
Thomas Witkowski's avatar
Thomas Witkowski committed
202
		} else if (*(macroNeighbour->getElement()->getDOF(2)) == *(element->getDOF(1))) {
203
204
		  neighbourCoord[i][0] = macroNeighbour->coord[0];
		  neighbourCoord[i][1] = macroNeighbour->coord[2];		 
205
206
207
208
		} else {
		  ERROR_EXIT("Should not happen!\n");
		}

209
210
211
		// 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.
212
213
		//		ERROR_EXIT("Should not happen!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
214
215
		break;

216
	      default:
Thomas Witkowski's avatar
Thomas Witkowski committed
217
		std::cout << "------------- Error --------------" << std::endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
218
		std::cout << "  Neighbour counter = " << i << "\n";
Thomas Witkowski's avatar
Thomas Witkowski committed
219
		std::cout << "  Element index     = " << element->getIndex() << "\n\n";
Thomas Witkowski's avatar
Thomas Witkowski committed
220
221
222
223
224
225
226
227
		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;
		  }
228
229
230
		  std::cout << "  OppVertex " << j << ": " 
			    << static_cast<int>(mel->getOppVertex(j)) 
			    << std::endl << std::endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
231
		}
232
233
		ERROR_EXIT("should not happen!\n");
		break;
234
235
236
	      }
	    }
	  } else {
237
238
239
240
241
242
243

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

244
	    if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
245
	      oppCoord[i] = macroNeighbour->coord[edgeNo];
246
	      neighbourCoord[i] = macroNeighbour->coord;	      
247
	    }
248
	  }
249
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
250
	  neighbour[i] = NULL;
251
        }
252
      }
253
    }
254
    
Thomas Witkowski's avatar
Thomas Witkowski committed
255
256
257
258
259
    if (fillFlag.isSet(Mesh::FILL_BOUND)) {   
      for (int i = 0; i < element->getGeo(BOUNDARY); i++)
	boundary[i] = mel->getBoundary(i);
      for (int i = 0; i < element->getGeo(PROJECTION); i++)
	projection[i] = mel->getProjection(i);      
260
261
262
263
264
265
266
267
    }
  }


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

268
  void ElInfo2d::fillElInfo(int ichild, const ElInfo *elInfoOld)
269
  {
270
    FUNCNAME("ElInfo::fillElInfo()");
271

Thomas Witkowski's avatar
Thomas Witkowski committed
272
    Element *elem = elInfoOld->element;
273
274
    Element *nb;

Thomas Witkowski's avatar
Thomas Witkowski committed
275
    Flag fill_flag = elInfoOld->fillFlag;
276

277
    TEST_EXIT_DBG(elem->getFirstChild())("no children?\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
278
    element = const_cast<Element*>((ichild == 0) ? 
279
280
				    elem->getFirstChild() : 
				    elem->getSecondChild());
Thomas Witkowski's avatar
Thomas Witkowski committed
281
    TEST_EXIT_DBG(element)("missing child %d?\n", ichild);
282

Thomas Witkowski's avatar
Thomas Witkowski committed
283
284
285
    macroElement  = elInfoOld->macroElement;
    fillFlag = fill_flag;
    parent = elem;
286
    level = elInfoOld->level + 1;
287
    iChild = ichild;
288

Thomas Witkowski's avatar
Thomas Witkowski committed
289
290
291
    if (fillFlag.isSet(Mesh::FILL_COORDS) || 
	fillFlag.isSet(Mesh::FILL_DET)    ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
292
      
293
      if (elem->isNewCoordSet())
Thomas Witkowski's avatar
Thomas Witkowski committed
294
	coord[2] = *(elem->getNewCoord());
295
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
296
	coord[2].setMidpoint(elInfoOld->coord[0], elInfoOld->coord[1]);      
297
298
      
      if (ichild == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
299
300
	coord[0] = elInfoOld->coord[2];
	coord[1] = elInfoOld->coord[0];
301
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
302
303
	coord[0] = elInfoOld->coord[1];
	coord[1] = elInfoOld->coord[2];
304
305
306
307
308
309
310
311
312
313
      }
    }

    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.

Thomas Witkowski's avatar
Thomas Witkowski committed
314
	neighbour[2] = elInfoOld->neighbour[1];
315
	oppVertex[2] = elInfoOld->oppVertex[1];
316
	
Thomas Witkowski's avatar
Thomas Witkowski committed
317
318
	if (neighbour[2] && fill_opp_coords) {
	  oppCoord[2] = elInfoOld->oppCoord[1];
319
	  neighbourCoord[2] = elInfoOld->neighbourCoord[1];
320
	}
321
322
323
324
325
326
327
328
329
	
	
	// Calculation of the neighbour 1, its oppCoords and the
	// cooresponding oppVertex.
	
	if (elem->getFirstChild()  &&  
	    elem->getSecondChild()->getFirstChild()  &&  
	    elem->getSecondChild()->getFirstChild()) {
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
330
	  neighbour[1] = elem->getSecondChild()->getSecondChild();
331
	  oppVertex[1] = 2;
332
333
	  
	  if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
334
335
336
337
338
339
340
341
            if (elem->getSecondChild()->isNewCoordSet())
	      oppCoord[1] = *(elem->getSecondChild()->getNewCoord());
	    else
	      oppCoord[1].setMidpoint(elInfoOld->coord[1], elInfoOld->coord[2]);

	    neighbourCoord[1][0] = coord[0];
	    neighbourCoord[1][1] = coord[2];
	    neighbourCoord[1][2] = oppCoord[1];  
342
343
	  }
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
344
	  neighbour[1] = elem->getSecondChild();
345
	  oppVertex[1] = 0;
346
347

	  if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
348
	    oppCoord[1] = elInfoOld->coord[1];
349

Thomas Witkowski's avatar
Thomas Witkowski committed
350
351
352
	    neighbourCoord[1][0] = elInfoOld->coord[1];
	    neighbourCoord[1][1] = elInfoOld->coord[2];
	    neighbourCoord[1][2] = coord[2];
353
354
355
356
	  }
	}


357
358
359
	// Calculation of the neighbour 0, its oppCoords and the
	// cooresponding oppVertex.
	
Thomas Witkowski's avatar
Thomas Witkowski committed
360
	nb = elInfoOld->neighbour[2];
361
	if (nb) {
362
	  TEST(elInfoOld->oppVertex[2] == 2)("invalid neighbour\n"); 
363
364
	  TEST_EXIT_DBG(nb->getFirstChild())("missing first child?\n");
	  TEST_EXIT_DBG(nb->getSecondChild())("missing second child?\n");
365
366
367
368
	 
	  nb = nb->getSecondChild();

	  if (nb->getFirstChild()) {
369
	    oppVertex[0] = 2;
370

371
	    if (fill_opp_coords) {
372
	      if (nb->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
373
		oppCoord[0] = *(nb->getNewCoord());
374
	      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
375
		oppCoord[0].setMidpoint(elInfoOld->neighbourCoord[2][1],
376
					 elInfoOld->neighbourCoord[2][2]);
377
	      }
378

379
380
381
	      neighbourCoord[0][0].setMidpoint(elInfoOld->neighbourCoord[2][0],
						elInfoOld->neighbourCoord[2][1]);
	      neighbourCoord[0][1] = elInfoOld->neighbourCoord[2][1];
Thomas Witkowski's avatar
Thomas Witkowski committed
382
	      neighbourCoord[0][2] = oppCoord[0];
383
384
385
386
	    }	   
 
	    nb = nb->getFirstChild();
	  } else {
387
	    oppVertex[0] = 1;
388

389
	    if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
390
	      oppCoord[0] = elInfoOld->oppCoord[2];    
391

392
393
394
395
	      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]);
396
397
	    }
	  }
398
399
	}
	
Thomas Witkowski's avatar
Thomas Witkowski committed
400
	neighbour[0] = nb;
401
402
403
404
      } else {   /* ichild == 1 */
	// Calculation of the neighbour 2, its oppCoords and the
	// cooresponding oppVertex.

Thomas Witkowski's avatar
Thomas Witkowski committed
405
	neighbour[2] = elInfoOld->neighbour[0];
406
	oppVertex[2] = elInfoOld->oppVertex[0];
407

Thomas Witkowski's avatar
Thomas Witkowski committed
408
409
	if (neighbour[2] && fill_opp_coords) {
	  oppCoord[2] = elInfoOld->oppCoord[0];
410
	  neighbourCoord[2] = elInfoOld->neighbourCoord[0];
411
412
413
414
415
416
417
	}
	

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

	if (elem->getFirstChild()->getFirstChild()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
418
	  neighbour[0] = elem->getFirstChild()->getFirstChild();
419
	  oppVertex[0] = 2;
420
421

	  if (fill_opp_coords) {
422
            if (elem->getFirstChild()->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
423
	      oppCoord[0] = *(elem->getFirstChild()->getNewCoord());
424
	    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
425
426
	      oppCoord[0].setMidpoint(elInfoOld->coord[0], 
				       elInfoOld->coord[2]);
427
	    }
428

Thomas Witkowski's avatar
Thomas Witkowski committed
429
430
431
	    neighbourCoord[0][0] = coord[2];
	    neighbourCoord[0][1] = coord[1];
	    neighbourCoord[0][2] = oppCoord[0];
432
	  }
433
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
434
	  neighbour[0] = elem->getFirstChild();
435
	  oppVertex[0] = 1;
436
437

	  if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
438
	    oppCoord[0] = elInfoOld->coord[0];
439

Thomas Witkowski's avatar
Thomas Witkowski committed
440
441
442
	    neighbourCoord[0][0] = elInfoOld->coord[2];
	    neighbourCoord[0][1] = elInfoOld->coord[0];
	    neighbourCoord[0][2] = coord[2];
443
	  }
444
445
446
447
448
	}

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

Thomas Witkowski's avatar
Thomas Witkowski committed
449
	nb = elInfoOld->neighbour[2];
450
	if (nb) {
451
	  TEST(elInfoOld->oppVertex[2] == 2)("invalid neighbour\n"); 
452
453
454
	  TEST((nb = nb->getFirstChild()))("missing child?\n");

	  if (nb->getFirstChild()) {
455
	    oppVertex[1] = 2;
456

457
	    if (fill_opp_coords) {
458
	      if (nb->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
459
		oppCoord[1] = *(nb->getNewCoord());
460
	      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
461
		oppCoord[1].setMidpoint(elInfoOld->neighbourCoord[2][0],
462
					 elInfoOld->neighbourCoord[2][2]);
463
	      }
464

465
466
467
	      neighbourCoord[1][0] = elInfoOld->neighbourCoord[2][0];
	      neighbourCoord[1][1].setMidpoint(elInfoOld->neighbourCoord[2][0],
					       elInfoOld->neighbourCoord[2][1]);
Thomas Witkowski's avatar
Thomas Witkowski committed
468
	      neighbourCoord[1][2] = oppCoord[1];
469
	    }
470
471
	    nb = nb->getSecondChild();

472
	  } else {
473
	    oppVertex[1] = 0;
474

475
	    if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
476
	      oppCoord[1] = elInfoOld->oppCoord[2];
477

478
479
480
481
	      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]);
482
483
484
	    }
	  }
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
485
	neighbour[1] = nb;
486
      } // if (ichild == 0) {} else
Thomas Witkowski's avatar
Thomas Witkowski committed
487
    } // if (fill_flag.isSet(Mesh::FILL_NEIGH) || fillFlag.isSet(Mesh::FILL_OPP_COORDS))
488
489
    

490
    if (fill_flag.isSet(Mesh::FILL_BOUND)) {
491
      if (elInfoOld->getBoundary(2))
Thomas Witkowski's avatar
Thomas Witkowski committed
492
	boundary[5] = elInfoOld->getBoundary(2);
493
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
494
	boundary[5] = INTERIOR;
495

496
      if (ichild == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
497
498
499
500
501
	boundary[3] = elInfoOld->getBoundary(5);
	boundary[4] = elInfoOld->getBoundary(3);
	boundary[0] = elInfoOld->getBoundary(2);
	boundary[1] = INTERIOR;
	boundary[2] = elInfoOld->getBoundary(1);
502
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
503
504
505
506
507
	boundary[3] = elInfoOld->getBoundary(4);
	boundary[4] = elInfoOld->getBoundary(5);
	boundary[0] = INTERIOR;
	boundary[1] = elInfoOld->getBoundary(2);
	boundary[2] = elInfoOld->getBoundary(0);
508
509
      }

510
511
      if (elInfoOld->getProjection(0) && 
	  elInfoOld->getProjection(0)->getType() == VOLUME_PROJECTION) {
512
	
Thomas Witkowski's avatar
Thomas Witkowski committed
513
	projection[0] = elInfoOld->getProjection(0);
514
515
      } else { // boundary projection
	if (ichild == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
516
517
518
	  projection[0] = elInfoOld->getProjection(2);
	  projection[1] = NULL;
	  projection[2] = elInfoOld->getProjection(1);
519
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
520
521
522
	  projection[0] = NULL;
	  projection[1] = elInfoOld->getProjection(2);
	  projection[2] = elInfoOld->getProjection(0);
523
524
525
526
527
	}
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
528
  double ElInfo2d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
529
  {
530
    FUNCNAME("ElInfo2d::calcGrdLambda()");
531
532

    testFlag(Mesh::FILL_COORDS);
533
534
  
    double adet = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
535
    int dim = mesh->getDim();
536

Thomas Witkowski's avatar
Thomas Witkowski committed
537
    for (int i = 0; i < dimOfWorld; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
538
539
      (*e1)[i] = coord[1][i] - coord[0][i];
      (*e2)[i] = coord[2][i] - coord[0][i];
540
541
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
542
    if (dimOfWorld == 2) {
543
      double sdet = (*e1)[0] * (*e2)[1] - (*e1)[1] * (*e2)[0];
544
545
546
547
      adet = abs(sdet);

      if (adet < 1.0E-25) {
	MSG("abs(det) = %f\n", adet);
548
	for (int i = 0; i <= dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
549
	  for (int j = 0; j < dimOfWorld; j++)
550
	    grd_lam[i][j] = 0.0;
551
552
      } else {
	double det1 = 1.0 / sdet;
553
554
555
556
557

	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
558
559
	grd_lam[0][0] = - grd_lam[1][0] - grd_lam[2][0];
	grd_lam[0][1] = - grd_lam[1][1] - grd_lam[2][1];
560
      }
561
562
    } else {  
      vectorProduct(*e1, *e2, *normal);
563

564
      adet = norm(normal);
565
566
567

      if (adet < 1.0E-15) {
	MSG("abs(det) = %lf\n", adet);
568
	for (int i = 0; i <= dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
569
	  for (int j = 0; j < dimOfWorld; j++)
570
571
	    grd_lam[i][j] = 0.0;
      } else {
572
573
	vectorProduct(*e2, *normal, grd_lam[1]);
	vectorProduct(*normal, *e1, grd_lam[2]);
574
      
575
	double adet2 = 1.0 / (adet * adet);
576

Thomas Witkowski's avatar
Thomas Witkowski committed
577
	for (int i = 0; i < dimOfWorld; i++) {
578
579
580
581
582
583
584
585
586
587
588
589
590
	  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;
  }

591
592

  const int ElInfo2d::worldToCoord(const WorldVector<double>& xy, 
593
594
				   DimVec<double>* lambda) const
  {
595
    FUNCNAME("ElInfo::worldToCoord()");
596

597
    TEST_EXIT_DBG(lambda)("lambda must not be NULL\n");
598

Thomas Witkowski's avatar
Thomas Witkowski committed
599
    DimVec<WorldVector<double> > edge(mesh->getDim(), NO_INIT);
600
    WorldVector<double> x; 
Thomas Witkowski's avatar
Thomas Witkowski committed
601
    static DimVec<double> vec(mesh->getDim(), NO_INIT);
602

Thomas Witkowski's avatar
Thomas Witkowski committed
603
    int dim = mesh->getDim();
604

605
    for (int j = 0; j < dimOfWorld; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
606
      double x0 = coord[dim][j];
607
      x[j] = xy[j] - x0;
608
      for (int i = 0; i < dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
609
	edge[i][j] = coord[i][j] - x0;
610
611
    }
  
612
613
614
    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]; 
615
616
617

    if (abs(det) < DBL_TOL) {
      ERROR("det = %le; abort\n", det);
618
619
      for (int i = 0; i <= dim; i++) 
	(*lambda)[i] = 1.0/dim;
620
621
622
623
624
625
626
627
628
      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;
629
    for (int i = 0; i <= dim; i++) {
630
631
632
633
634
635
636
637
638
639
640
641
      if ((*lambda)[i] < -1.E-5) {
	if ((*lambda)[i] < lmin) {
	  k = i;
	  lmin = (*lambda)[i];
	}
      }
    }

    return k;
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
642
  double ElInfo2d::getNormal(int side, WorldVector<double> &normal)
643
  {
644
    FUNCNAME("ElInfo::getNormal()");
645

646
647
    int i0 = (side + 1) % 3;
    int i1 = (side + 2) % 3;
648

Thomas Witkowski's avatar
Thomas Witkowski committed
649
    if (dimOfWorld == 2){
Thomas Witkowski's avatar
Thomas Witkowski committed
650
651
      normal[0] = coord[i1][1] - coord[i0][1];
      normal[1] = coord[i0][0] - coord[i1][0];
652
653
654
    } else { // dow == 3
      WorldVector<double> e0, e1,e2, elementNormal;

Thomas Witkowski's avatar
Thomas Witkowski committed
655
656
657
658
659
660
      e0 = coord[i1]; 
      e0 -= coord[i0];
      e1 = coord[i1]; 
      e1 -= coord[side];
      e2 = coord[i0]; 
      e2 -= coord[side];
661
662
663
664
665

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

666
    double det = norm(&normal);
667

668
    TEST_EXIT_DBG(det > 1.e-30)("det = 0 on face %d\n", side);
669
670
671

    normal *= 1.0 / det;
    
Thomas Witkowski's avatar
Thomas Witkowski committed
672
    return det;
673
674
  }

675

676
677
678
679
680
  /****************************************************************************/
  /*  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
681
  double ElInfo2d::getElementNormal(WorldVector<double> &elementNormal) const
682
  {
683
    FUNCNAME("ElInfo::getElementNormal()");
684

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

Thomas Witkowski's avatar
Thomas Witkowski committed
688
689
    WorldVector<double> e0 = coord[1] - coord[0];
    WorldVector<double> e1 = coord[2] - coord[0];
690
691
692

    vectorProduct(e0, e1, elementNormal);

693
    double det = norm(&elementNormal);
694

695
    TEST_EXIT_DBG(det > 1.e-30)("det = 0");
696
697
698

    elementNormal *= 1.0 / det;
    
Thomas Witkowski's avatar
Thomas Witkowski committed
699
    return det;
700
  }
701

702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737

  void ElInfo2d::getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const
  {
    FUNCNAME("ElInfo2d::getSubElemCoordsMat()");

    using namespace mtl;

    switch (degree) {
    case 1:
      {
      mat = mat_d1;
      dense2D<double> tmpMat(num_rows(mat), num_rows(mat));

      for (int i = 0; i < refinementPathLength; i++) {
	if (refinementPath & (1 << i)) {
	  tmpMat = mat_d1_right * mat;
	  mat = tmpMat;
	} else  {
	  tmpMat = mat_d1_left * mat;
	  mat = tmpMat;
	}
      }
      }
      break;
    default:
      ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
    }
  }


  void ElInfo2d::getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const
  {
    FUNCNAME("ElInfo2d::getSubElemCoordsMat_so()");

    ERROR_EXIT("Not yet implemented!\n");
  }
738
}