ElInfo2d.cc 19.9 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
  					    {0.0, 0.0, 0.5},
25
  					    {1.0, 0.0, 0.0}}; 
26
27
  mtl::dense2D<double> ElInfo2d::mat_d1_left(mat_d1_left_val);

28
  
29
  double ElInfo2d::mat_d1_right_val[3][3] = {{0.0, 0.0, 0.5}, 
30
  					     {1.0, 0.0, 0.5},
31
  					     {0.0, 1.0, 0.0}}; 
32
33
  mtl::dense2D<double> ElInfo2d::mat_d1_right(mat_d1_right_val);

34

35
36
37
  ElInfo2d::ElInfo2d(Mesh *aMesh) 
    : ElInfo(aMesh) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
38
39
40
    e1 = new WorldVector<double>;
    e2 = new WorldVector<double>;
    normal = new WorldVector<double>;
41
42
  }

43

44
45
  ElInfo2d::~ElInfo2d()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
46
47
48
    delete e1;
    delete e2;
    delete normal;
49
50
  }

51

52
53
  void ElInfo2d::fillMacroInfo(const MacroElement * mel)
  {
54
55
    FUNCNAME("ElInfo::fillMacroInfo()");
 
Thomas Witkowski's avatar
Thomas Witkowski committed
56
57
58
    macroElement = const_cast<MacroElement*>(mel);
    element = const_cast<Element*>(mel->getElement());
    parent = NULL;
59
    level = 0;
60

Thomas Witkowski's avatar
Thomas Witkowski committed
61
62
63
    if (fillFlag.isSet(Mesh::FILL_COORDS) || 
	fillFlag.isSet(Mesh::FILL_DET)    ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
64

Thomas Witkowski's avatar
Thomas Witkowski committed
65
      int vertices = mesh->getGeo(VERTEX);
66
      for (int i = 0; i < vertices; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
67
	coord[i] = mel->coord[i];
68
69
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
72
73
    if (fillFlag.isSet(Mesh::FILL_OPP_COORDS) || 
	fillFlag.isSet(Mesh::FILL_NEIGH)) {
74

Thomas Witkowski's avatar
Thomas Witkowski committed
75
      bool fill_opp_coords = (fillFlag.isSet(Mesh::FILL_OPP_COORDS));
76
    
77
78
79
80
      for (int i = 0; i < neighbours; i++) {
	MacroElement *macroNeighbour = mel->getNeighbour(i);

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

84
	  int edgeNo = oppVertex[i] = mel->getOppVertex(i);
85
86


87
	  if (nb->getFirstChild() && edgeNo != 2) {  
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107

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


108
	    if (edgeNo == 0) {
109
110
111
112
113
114
115
116
117
118
119
120
	      // The situation is as follows:
	      //
	      //          -------
	      //          \    /|\
	      //           \  / | \
	      //            \/  |  \
	      //             \  |   \
	      //              \ |    \ 
	      //               \|     \
	      //                -------
	      //
	      //            nb     el
121
              //
122
123
124
125
	      // 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
126
	      nb = neighbour[i] = nb->getSecondChild();
127
	    } else {
128
129
	      // 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
130
	      nb = neighbour[i] = nb->getFirstChild();
131
132
	    }

133
134
	    // In both cases the opp vertex number is 2, as one can see in the 
	    // pictures above.
135
	    oppVertex[i] = 2;
136
137

	    if (fill_opp_coords) {
138
	      if (nb->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
139
		oppCoord[i] = *(nb->getNewCoord());
140
	      } else {
141
142
143
		// 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
144
		oppCoord[i] = (macroNeighbour->coord[0] + 
145
146
				macroNeighbour->coord[1]) * 0.5;
	      }
147
148
149
	      
	      switch (i) {
	      case 0:
150
151
152
153
		// The common edge is the edge 0 of this element.

		switch (edgeNo) {
		case 1:
154
155
		  neighbourCoord[i][0] = macroNeighbour->coord[2];
		  neighbourCoord[i][1] = macroNeighbour->coord[0];
156
157
		  break;
		case 0:		  
158
159
		  neighbourCoord[i][0] = macroNeighbour->coord[1];
		  neighbourCoord[i][1] = macroNeighbour->coord[2];
160
161
		  break;
		default:
Thomas Witkowski's avatar
Thomas Witkowski committed
162
		  ERROR_EXIT("Should not happen!\n");
163
		}
164
	
Thomas Witkowski's avatar
Thomas Witkowski committed
165
		neighbourCoord[i][2] = oppCoord[i];
166
167
168
		break;
		
	      case 1:
169
170
171
		// The commonedge is the edge 1 of this element.
		switch (edgeNo) {
		case 0:
172
173
		  neighbourCoord[i][0] = macroNeighbour->coord[1];
		  neighbourCoord[i][1] = macroNeighbour->coord[2];
174
175
		  break;
		case 1:
176
177
		  neighbourCoord[i][0] = macroNeighbour->coord[2];
		  neighbourCoord[i][1] = macroNeighbour->coord[0];
178
179
		  break;
		default:
Thomas Witkowski's avatar
Thomas Witkowski committed
180
		  ERROR_EXIT("Should not happen!\n");
181
182
		}
		
Thomas Witkowski's avatar
Thomas Witkowski committed
183
		neighbourCoord[i][2] = oppCoord[i];
184
185
		break;
		
Thomas Witkowski's avatar
Thomas Witkowski committed
186
	      case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
187
		if (*(macroNeighbour->getElement()->getDOF(2)) == *(element->getDOF(0))) {
188
189
		  neighbourCoord[i][0] = macroNeighbour->coord[2];
		  neighbourCoord[i][1] = macroNeighbour->coord[1];
Thomas Witkowski's avatar
Thomas Witkowski committed
190
		} else if (*(macroNeighbour->getElement()->getDOF(2)) == *(element->getDOF(1))) {
191
192
		  neighbourCoord[i][0] = macroNeighbour->coord[0];
		  neighbourCoord[i][1] = macroNeighbour->coord[2];		 
193
194
195
196
		} else {
		  ERROR_EXIT("Should not happen!\n");
		}

197
198
199
		// 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.
200
201
		//		ERROR_EXIT("Should not happen!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
202
203
		break;

204
	      default:
Thomas Witkowski's avatar
Thomas Witkowski committed
205
		std::cout << "------------- Error --------------" << std::endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
206
		std::cout << "  Neighbour counter = " << i << "\n";
Thomas Witkowski's avatar
Thomas Witkowski committed
207
		std::cout << "  Element index     = " << element->getIndex() << "\n\n";
Thomas Witkowski's avatar
Thomas Witkowski committed
208
209
210
211
212
213
214
215
		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;
		  }
216
217
218
		  std::cout << "  OppVertex " << j << ": " 
			    << static_cast<int>(mel->getOppVertex(j)) 
			    << std::endl << std::endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
219
		}
220
221
		ERROR_EXIT("should not happen!\n");
		break;
222
223
224
	      }
	    }
	  } else {
225
226
227
228
229
230
231

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

232
	    if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
233
	      oppCoord[i] = macroNeighbour->coord[edgeNo];
234
	      neighbourCoord[i] = macroNeighbour->coord;	      
235
	    }
236
	  }
237
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
238
	  neighbour[i] = NULL;
239
        }
240
      }
241
    }
242
    
Thomas Witkowski's avatar
Thomas Witkowski committed
243
244
245
246
247
    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);      
248
249
250
251
252
253
254
255
    }
  }


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

256
  void ElInfo2d::fillElInfo(int ichild, const ElInfo *elInfoOld)
257
  {
258
    FUNCNAME("ElInfo::fillElInfo()");
259

Thomas Witkowski's avatar
Thomas Witkowski committed
260
    Element *elem = elInfoOld->element;
261
262
    Element *nb;

Thomas Witkowski's avatar
Thomas Witkowski committed
263
    Flag fill_flag = elInfoOld->fillFlag;
264

265
    TEST_EXIT_DBG(elem->getFirstChild())("no children?\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
266
    element = const_cast<Element*>((ichild == 0) ? 
267
268
				    elem->getFirstChild() : 
				    elem->getSecondChild());
Thomas Witkowski's avatar
Thomas Witkowski committed
269
    TEST_EXIT_DBG(element)("missing child %d?\n", ichild);
270

Thomas Witkowski's avatar
Thomas Witkowski committed
271
272
273
    macroElement  = elInfoOld->macroElement;
    fillFlag = fill_flag;
    parent = elem;
274
    level = elInfoOld->level + 1;
275
    iChild = ichild;
276

Thomas Witkowski's avatar
Thomas Witkowski committed
277
278
279
    if (fillFlag.isSet(Mesh::FILL_COORDS) || 
	fillFlag.isSet(Mesh::FILL_DET)    ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
280
      
281
      if (elem->isNewCoordSet())
Thomas Witkowski's avatar
Thomas Witkowski committed
282
	coord[2] = *(elem->getNewCoord());
283
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
284
	coord[2].setMidpoint(elInfoOld->coord[0], elInfoOld->coord[1]);      
285
286
      
      if (ichild == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
287
288
	coord[0] = elInfoOld->coord[2];
	coord[1] = elInfoOld->coord[0];
289
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
290
291
	coord[0] = elInfoOld->coord[1];
	coord[1] = elInfoOld->coord[2];
292
293
294
295
296
297
298
299
300
301
      }
    }

    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
302
	neighbour[2] = elInfoOld->neighbour[1];
303
	oppVertex[2] = elInfoOld->oppVertex[1];
304
	
Thomas Witkowski's avatar
Thomas Witkowski committed
305
306
	if (neighbour[2] && fill_opp_coords) {
	  oppCoord[2] = elInfoOld->oppCoord[1];
307
	  neighbourCoord[2] = elInfoOld->neighbourCoord[1];
308
	}
309
310
311
312
313
314
315
316
317
	
	
	// 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
318
	  neighbour[1] = elem->getSecondChild()->getSecondChild();
319
	  oppVertex[1] = 2;
320
321
	  
	  if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
322
323
324
325
326
327
328
329
            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];  
330
331
	  }
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
332
	  neighbour[1] = elem->getSecondChild();
333
	  oppVertex[1] = 0;
334
335

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

Thomas Witkowski's avatar
Thomas Witkowski committed
338
339
340
	    neighbourCoord[1][0] = elInfoOld->coord[1];
	    neighbourCoord[1][1] = elInfoOld->coord[2];
	    neighbourCoord[1][2] = coord[2];
341
342
343
344
	  }
	}


345
346
347
	// Calculation of the neighbour 0, its oppCoords and the
	// cooresponding oppVertex.
	
Thomas Witkowski's avatar
Thomas Witkowski committed
348
	nb = elInfoOld->neighbour[2];
349
	if (nb) {
350
	  TEST(elInfoOld->oppVertex[2] == 2)("invalid neighbour\n"); 
351
352
	  TEST_EXIT_DBG(nb->getFirstChild())("missing first child?\n");
	  TEST_EXIT_DBG(nb->getSecondChild())("missing second child?\n");
353
354
355
356
	 
	  nb = nb->getSecondChild();

	  if (nb->getFirstChild()) {
357
	    oppVertex[0] = 2;
358

359
	    if (fill_opp_coords) {
360
	      if (nb->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
361
		oppCoord[0] = *(nb->getNewCoord());
362
	      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
363
		oppCoord[0].setMidpoint(elInfoOld->neighbourCoord[2][1],
364
					 elInfoOld->neighbourCoord[2][2]);
365
	      }
366

367
368
369
	      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
370
	      neighbourCoord[0][2] = oppCoord[0];
371
372
373
374
	    }	   
 
	    nb = nb->getFirstChild();
	  } else {
375
	    oppVertex[0] = 1;
376

377
	    if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
378
	      oppCoord[0] = elInfoOld->oppCoord[2];    
379

380
381
382
383
	      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]);
384
385
	    }
	  }
386
387
	}
	
Thomas Witkowski's avatar
Thomas Witkowski committed
388
	neighbour[0] = nb;
389
390
391
392
      } else {   /* ichild == 1 */
	// Calculation of the neighbour 2, its oppCoords and the
	// cooresponding oppVertex.

Thomas Witkowski's avatar
Thomas Witkowski committed
393
	neighbour[2] = elInfoOld->neighbour[0];
394
	oppVertex[2] = elInfoOld->oppVertex[0];
395

Thomas Witkowski's avatar
Thomas Witkowski committed
396
397
	if (neighbour[2] && fill_opp_coords) {
	  oppCoord[2] = elInfoOld->oppCoord[0];
398
	  neighbourCoord[2] = elInfoOld->neighbourCoord[0];
399
400
401
402
403
404
405
	}
	

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

	if (elem->getFirstChild()->getFirstChild()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
406
	  neighbour[0] = elem->getFirstChild()->getFirstChild();
407
	  oppVertex[0] = 2;
408
409

	  if (fill_opp_coords) {
410
            if (elem->getFirstChild()->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
411
	      oppCoord[0] = *(elem->getFirstChild()->getNewCoord());
412
	    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
413
414
	      oppCoord[0].setMidpoint(elInfoOld->coord[0], 
				       elInfoOld->coord[2]);
415
	    }
416

Thomas Witkowski's avatar
Thomas Witkowski committed
417
418
419
	    neighbourCoord[0][0] = coord[2];
	    neighbourCoord[0][1] = coord[1];
	    neighbourCoord[0][2] = oppCoord[0];
420
	  }
421
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
422
	  neighbour[0] = elem->getFirstChild();
423
	  oppVertex[0] = 1;
424
425

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

Thomas Witkowski's avatar
Thomas Witkowski committed
428
429
430
	    neighbourCoord[0][0] = elInfoOld->coord[2];
	    neighbourCoord[0][1] = elInfoOld->coord[0];
	    neighbourCoord[0][2] = coord[2];
431
	  }
432
433
434
435
436
	}

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

Thomas Witkowski's avatar
Thomas Witkowski committed
437
	nb = elInfoOld->neighbour[2];
438
	if (nb) {
439
	  TEST(elInfoOld->oppVertex[2] == 2)("invalid neighbour\n"); 
440
441
442
	  TEST((nb = nb->getFirstChild()))("missing child?\n");

	  if (nb->getFirstChild()) {
443
	    oppVertex[1] = 2;
444

445
	    if (fill_opp_coords) {
446
	      if (nb->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
447
		oppCoord[1] = *(nb->getNewCoord());
448
	      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
449
		oppCoord[1].setMidpoint(elInfoOld->neighbourCoord[2][0],
450
					 elInfoOld->neighbourCoord[2][2]);
451
	      }
452

453
454
455
	      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
456
	      neighbourCoord[1][2] = oppCoord[1];
457
	    }
458
459
	    nb = nb->getSecondChild();

460
	  } else {
461
	    oppVertex[1] = 0;
462

463
	    if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
464
	      oppCoord[1] = elInfoOld->oppCoord[2];
465

466
467
468
469
	      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]);
470
471
472
	    }
	  }
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
473
	neighbour[1] = nb;
474
      } // if (ichild == 0) {} else
Thomas Witkowski's avatar
Thomas Witkowski committed
475
    } // if (fill_flag.isSet(Mesh::FILL_NEIGH) || fillFlag.isSet(Mesh::FILL_OPP_COORDS))
476
477
    

478
    if (fill_flag.isSet(Mesh::FILL_BOUND)) {
479
      if (elInfoOld->getBoundary(2))
Thomas Witkowski's avatar
Thomas Witkowski committed
480
	boundary[5] = elInfoOld->getBoundary(2);
481
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
482
	boundary[5] = INTERIOR;
483

484
      if (ichild == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
485
486
487
488
489
	boundary[3] = elInfoOld->getBoundary(5);
	boundary[4] = elInfoOld->getBoundary(3);
	boundary[0] = elInfoOld->getBoundary(2);
	boundary[1] = INTERIOR;
	boundary[2] = elInfoOld->getBoundary(1);
490
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
491
492
493
494
495
	boundary[3] = elInfoOld->getBoundary(4);
	boundary[4] = elInfoOld->getBoundary(5);
	boundary[0] = INTERIOR;
	boundary[1] = elInfoOld->getBoundary(2);
	boundary[2] = elInfoOld->getBoundary(0);
496
497
      }

498
499
      if (elInfoOld->getProjection(0) && 
	  elInfoOld->getProjection(0)->getType() == VOLUME_PROJECTION) {
500
	
Thomas Witkowski's avatar
Thomas Witkowski committed
501
	projection[0] = elInfoOld->getProjection(0);
502
503
      } else { // boundary projection
	if (ichild == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
504
505
506
	  projection[0] = elInfoOld->getProjection(2);
	  projection[1] = NULL;
	  projection[2] = elInfoOld->getProjection(1);
507
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
508
509
510
	  projection[0] = NULL;
	  projection[1] = elInfoOld->getProjection(2);
	  projection[2] = elInfoOld->getProjection(0);
511
512
513
514
515
	}
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
516
  double ElInfo2d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
517
  {
518
    FUNCNAME("ElInfo2d::calcGrdLambda()");
519
520

    testFlag(Mesh::FILL_COORDS);
521
522
  
    double adet = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
523
    int dim = mesh->getDim();
524

Thomas Witkowski's avatar
Thomas Witkowski committed
525
    for (int i = 0; i < dimOfWorld; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
526
527
      (*e1)[i] = coord[1][i] - coord[0][i];
      (*e2)[i] = coord[2][i] - coord[0][i];
528
529
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
530
    if (dimOfWorld == 2) {
531
      double sdet = (*e1)[0] * (*e2)[1] - (*e1)[1] * (*e2)[0];
532
533
534
535
      adet = abs(sdet);

      if (adet < 1.0E-25) {
	MSG("abs(det) = %f\n", adet);
536
	for (int i = 0; i <= dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
537
	  for (int j = 0; j < dimOfWorld; j++)
538
	    grd_lam[i][j] = 0.0;
539
540
      } else {
	double det1 = 1.0 / sdet;
541
542
543
544
545

	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
546
547
	grd_lam[0][0] = - grd_lam[1][0] - grd_lam[2][0];
	grd_lam[0][1] = - grd_lam[1][1] - grd_lam[2][1];
548
      }
549
550
    } else {  
      vectorProduct(*e1, *e2, *normal);
551

552
      adet = norm(normal);
553
554
555

      if (adet < 1.0E-15) {
	MSG("abs(det) = %lf\n", adet);
556
	for (int i = 0; i <= dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
557
	  for (int j = 0; j < dimOfWorld; j++)
558
559
	    grd_lam[i][j] = 0.0;
      } else {
560
561
	vectorProduct(*e2, *normal, grd_lam[1]);
	vectorProduct(*normal, *e1, grd_lam[2]);
562
      
563
	double adet2 = 1.0 / (adet * adet);
564

Thomas Witkowski's avatar
Thomas Witkowski committed
565
	for (int i = 0; i < dimOfWorld; i++) {
566
567
568
569
570
571
572
573
574
575
576
577
578
	  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;
  }

579
580

  const int ElInfo2d::worldToCoord(const WorldVector<double>& xy, 
581
582
				   DimVec<double>* lambda) const
  {
583
    FUNCNAME("ElInfo::worldToCoord()");
584

585
    TEST_EXIT_DBG(lambda)("lambda must not be NULL\n");
586

Thomas Witkowski's avatar
Thomas Witkowski committed
587
    DimVec<WorldVector<double> > edge(mesh->getDim(), NO_INIT);
588
    WorldVector<double> x; 
Thomas Witkowski's avatar
Thomas Witkowski committed
589
    static DimVec<double> vec(mesh->getDim(), NO_INIT);
590

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

593
    for (int j = 0; j < dimOfWorld; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
594
      double x0 = coord[dim][j];
595
      x[j] = xy[j] - x0;
596
      for (int i = 0; i < dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
597
	edge[i][j] = coord[i][j] - x0;
598
599
    }
  
600
601
602
    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]; 
603
604
605

    if (abs(det) < DBL_TOL) {
      ERROR("det = %le; abort\n", det);
606
607
      for (int i = 0; i <= dim; i++) 
	(*lambda)[i] = 1.0/dim;
608
609
610
611
612
613
614
615
616
      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;
617
    for (int i = 0; i <= dim; i++) {
618
619
620
621
622
623
624
625
626
627
628
629
      if ((*lambda)[i] < -1.E-5) {
	if ((*lambda)[i] < lmin) {
	  k = i;
	  lmin = (*lambda)[i];
	}
      }
    }

    return k;
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
630
  double ElInfo2d::getNormal(int side, WorldVector<double> &normal)
631
  {
632
    FUNCNAME("ElInfo::getNormal()");
633

634
635
    int i0 = (side + 1) % 3;
    int i1 = (side + 2) % 3;
636

Thomas Witkowski's avatar
Thomas Witkowski committed
637
    if (dimOfWorld == 2){
Thomas Witkowski's avatar
Thomas Witkowski committed
638
639
      normal[0] = coord[i1][1] - coord[i0][1];
      normal[1] = coord[i0][0] - coord[i1][0];
640
641
642
    } else { // dow == 3
      WorldVector<double> e0, e1,e2, elementNormal;

Thomas Witkowski's avatar
Thomas Witkowski committed
643
644
645
646
647
648
      e0 = coord[i1]; 
      e0 -= coord[i0];
      e1 = coord[i1]; 
      e1 -= coord[side];
      e2 = coord[i0]; 
      e2 -= coord[side];
649
650
651
652
653

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

654
    double det = norm(&normal);
655

656
    TEST_EXIT_DBG(det > 1.e-30)("det = 0 on face %d\n", side);
657
658
659

    normal *= 1.0 / det;
    
Thomas Witkowski's avatar
Thomas Witkowski committed
660
    return det;
661
662
  }

663

664
665
666
667
668
  /****************************************************************************/
  /*  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
669
  double ElInfo2d::getElementNormal(WorldVector<double> &elementNormal) const
670
  {
671
    FUNCNAME("ElInfo::getElementNormal()");
672

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

Thomas Witkowski's avatar
Thomas Witkowski committed
676
677
    WorldVector<double> e0 = coord[1] - coord[0];
    WorldVector<double> e1 = coord[2] - coord[0];
678
679
680

    vectorProduct(e0, e1, elementNormal);

681
    double det = norm(&elementNormal);
682

683
    TEST_EXIT_DBG(det > 1.e-30)("det = 0");
684
685
686

    elementNormal *= 1.0 / det;
    
Thomas Witkowski's avatar
Thomas Witkowski committed
687
    return det;
688
  }
689

690
691
692
693
694
695
696
697
698
699
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

  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");
  }
726
}