ElInfo2d.cc 22.3 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
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61


  double ElInfo2d::mat_d2_val[6][6] = {{1.0, 0.0, 0.0, 0.0, 0.0, 0.0}, 
				       {0.0, 1.0, 0.0, 0.0, 0.0, 0.0}, 
				       {0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
				       {0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
				       {0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
				       {0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
  mtl::dense2D<double> ElInfo2d::mat_d2(mat_d2_val);

  double ElInfo2d::mat_d2_left_val[6][6] = {{0.0, 1.0, 0.0, 0.375, -0.125, 0.0}, 
					    {0.0, 0.0, 0.0, -0.125, -0.125, 0.0}, 
					    {1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
					    {0.0, 0.0, 0.0, 0.0, 0.5, 0.0},
					    {0.0, 0.0, 0.0, 0.0, 0.5, 1.0},
					    {0.0, 0.0, 1.0, 0.75, 0.25, 0.0}};
  mtl::dense2D<double> ElInfo2d::mat_d2_left(mat_d2_left_val);

  double ElInfo2d::mat_d2_right_val[6][6] = {{0.0, 0.0, 0.0, -0.125, -0.125, 0.0}, 
					     {1.0, 0.0, 0.0, -0.125, 0.375, 0.0}, 
					     {0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
					     {0.0, 0.0, 0.0, 0.5, 0.0, 1.0},
					     {0.0, 0.0, 0.0, 0.5, 0.0, 0.0},
					     {0.0, 0.0, 1.0, 0.25, 0.75, 0.0}};
  mtl::dense2D<double> ElInfo2d::mat_d2_right(mat_d2_right_val);


62
63
64
  ElInfo2d::ElInfo2d(Mesh *aMesh) 
    : ElInfo(aMesh) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
65
66
67
    e1 = new WorldVector<double>;
    e2 = new WorldVector<double>;
    normal = new WorldVector<double>;
68
69
  }

70

71
72
  ElInfo2d::~ElInfo2d()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
73
74
75
    delete e1;
    delete e2;
    delete normal;
76
77
  }

78

79
80
  void ElInfo2d::fillMacroInfo(const MacroElement * mel)
  {
81
82
    FUNCNAME("ElInfo::fillMacroInfo()");
 
Thomas Witkowski's avatar
Thomas Witkowski committed
83
84
85
    macroElement = const_cast<MacroElement*>(mel);
    element = const_cast<Element*>(mel->getElement());
    parent = NULL;
86
    level = 0;
87

Thomas Witkowski's avatar
Thomas Witkowski committed
88
89
90
    if (fillFlag.isSet(Mesh::FILL_COORDS) || 
	fillFlag.isSet(Mesh::FILL_DET)    ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
91

Thomas Witkowski's avatar
Thomas Witkowski committed
92
      int vertices = mesh->getGeo(VERTEX);
93
      for (int i = 0; i < vertices; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
94
	coord[i] = mel->coord[i];
95
96
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
99
100
    if (fillFlag.isSet(Mesh::FILL_OPP_COORDS) || 
	fillFlag.isSet(Mesh::FILL_NEIGH)) {
101

Thomas Witkowski's avatar
Thomas Witkowski committed
102
      bool fill_opp_coords = (fillFlag.isSet(Mesh::FILL_OPP_COORDS));
103
    
104
105
106
107
      for (int i = 0; i < neighbours; i++) {
	MacroElement *macroNeighbour = mel->getNeighbour(i);

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

111
	  int edgeNo = oppVertex[i] = mel->getOppVertex(i);
112
113


114
	  if (nb->getFirstChild() && edgeNo != 2) {  
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134

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


135
	    if (edgeNo == 0) {
136
137
138
139
140
141
142
143
144
145
146
147
	      // The situation is as follows:
	      //
	      //          -------
	      //          \    /|\
	      //           \  / | \
	      //            \/  |  \
	      //             \  |   \
	      //              \ |    \ 
	      //               \|     \
	      //                -------
	      //
	      //            nb     el
148
              //
149
150
151
152
	      // 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
153
	      nb = neighbour[i] = nb->getSecondChild();
154
	    } else {
155
156
	      // 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
157
	      nb = neighbour[i] = nb->getFirstChild();
158
159
	    }

160
161
	    // In both cases the opp vertex number is 2, as one can see in the 
	    // pictures above.
162
	    oppVertex[i] = 2;
163
164

	    if (fill_opp_coords) {
165
	      if (nb->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
166
		oppCoord[i] = *(nb->getNewCoord());
167
	      } else {
168
169
170
		// 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
171
		oppCoord[i] = (macroNeighbour->coord[0] + 
172
173
				macroNeighbour->coord[1]) * 0.5;
	      }
174
175
176
	      
	      switch (i) {
	      case 0:
177
178
179
180
		// The common edge is the edge 0 of this element.

		switch (edgeNo) {
		case 1:
181
182
		  neighbourCoord[i][0] = macroNeighbour->coord[2];
		  neighbourCoord[i][1] = macroNeighbour->coord[0];
183
184
		  break;
		case 0:		  
185
186
		  neighbourCoord[i][0] = macroNeighbour->coord[1];
		  neighbourCoord[i][1] = macroNeighbour->coord[2];
187
188
		  break;
		default:
Thomas Witkowski's avatar
Thomas Witkowski committed
189
		  ERROR_EXIT("Should not happen!\n");
190
		}
191
	
Thomas Witkowski's avatar
Thomas Witkowski committed
192
		neighbourCoord[i][2] = oppCoord[i];
193
194
195
		break;
		
	      case 1:
196
197
198
		// The commonedge is the edge 1 of this element.
		switch (edgeNo) {
		case 0:
199
200
		  neighbourCoord[i][0] = macroNeighbour->coord[1];
		  neighbourCoord[i][1] = macroNeighbour->coord[2];
201
202
		  break;
		case 1:
203
204
		  neighbourCoord[i][0] = macroNeighbour->coord[2];
		  neighbourCoord[i][1] = macroNeighbour->coord[0];
205
206
		  break;
		default:
Thomas Witkowski's avatar
Thomas Witkowski committed
207
		  ERROR_EXIT("Should not happen!\n");
208
209
		}
		
Thomas Witkowski's avatar
Thomas Witkowski committed
210
		neighbourCoord[i][2] = oppCoord[i];
211
212
		break;
		
Thomas Witkowski's avatar
Thomas Witkowski committed
213
	      case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
214
		if (*(macroNeighbour->getElement()->getDOF(2)) == *(element->getDOF(0))) {
215
216
		  neighbourCoord[i][0] = macroNeighbour->coord[2];
		  neighbourCoord[i][1] = macroNeighbour->coord[1];
Thomas Witkowski's avatar
Thomas Witkowski committed
217
		} else if (*(macroNeighbour->getElement()->getDOF(2)) == *(element->getDOF(1))) {
218
219
		  neighbourCoord[i][0] = macroNeighbour->coord[0];
		  neighbourCoord[i][1] = macroNeighbour->coord[2];		 
220
221
222
223
		} else {
		  ERROR_EXIT("Should not happen!\n");
		}

224
225
226
		// 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.
227
228
		//		ERROR_EXIT("Should not happen!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
229
230
		break;

231
	      default:
Thomas Witkowski's avatar
Thomas Witkowski committed
232
		std::cout << "------------- Error --------------" << std::endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
233
		std::cout << "  Neighbour counter = " << i << "\n";
Thomas Witkowski's avatar
Thomas Witkowski committed
234
		std::cout << "  Element index     = " << element->getIndex() << "\n\n";
Thomas Witkowski's avatar
Thomas Witkowski committed
235
236
237
238
239
240
241
242
		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;
		  }
243
244
245
		  std::cout << "  OppVertex " << j << ": " 
			    << static_cast<int>(mel->getOppVertex(j)) 
			    << std::endl << std::endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
246
		}
247
248
		ERROR_EXIT("should not happen!\n");
		break;
249
250
251
	      }
	    }
	  } else {
252
253
254
255
256
257
258

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

259
	    if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
260
	      oppCoord[i] = macroNeighbour->coord[edgeNo];
261
	      neighbourCoord[i] = macroNeighbour->coord;	      
262
	    }
263
	  }
264
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
265
	  neighbour[i] = NULL;
266
        }
267
      }
268
    }
269
    
Thomas Witkowski's avatar
Thomas Witkowski committed
270
271
272
273
274
    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);      
275
276
277
278
279
280
281
282
    }
  }


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

283
  void ElInfo2d::fillElInfo(int ichild, const ElInfo *elInfoOld)
284
  {
285
    FUNCNAME("ElInfo::fillElInfo()");
286

Thomas Witkowski's avatar
Thomas Witkowski committed
287
    Element *elem = elInfoOld->element;
288
289
    Element *nb;

Thomas Witkowski's avatar
Thomas Witkowski committed
290
    Flag fill_flag = elInfoOld->fillFlag;
291

292
    TEST_EXIT_DBG(elem->getFirstChild())("no children?\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
293
    element = const_cast<Element*>((ichild == 0) ? 
294
295
				    elem->getFirstChild() : 
				    elem->getSecondChild());
Thomas Witkowski's avatar
Thomas Witkowski committed
296
    TEST_EXIT_DBG(element)("missing child %d?\n", ichild);
297

Thomas Witkowski's avatar
Thomas Witkowski committed
298
299
300
    macroElement  = elInfoOld->macroElement;
    fillFlag = fill_flag;
    parent = elem;
301
    level = elInfoOld->level + 1;
302
    iChild = ichild;
303

Thomas Witkowski's avatar
Thomas Witkowski committed
304
305
306
    if (fillFlag.isSet(Mesh::FILL_COORDS) || 
	fillFlag.isSet(Mesh::FILL_DET)    ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
307
      
308
      if (elem->isNewCoordSet())
Thomas Witkowski's avatar
Thomas Witkowski committed
309
	coord[2] = *(elem->getNewCoord());
310
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
311
	coord[2].setMidpoint(elInfoOld->coord[0], elInfoOld->coord[1]);      
312
313
      
      if (ichild == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
314
315
	coord[0] = elInfoOld->coord[2];
	coord[1] = elInfoOld->coord[0];
316
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
317
318
	coord[0] = elInfoOld->coord[1];
	coord[1] = elInfoOld->coord[2];
319
320
321
322
323
324
325
326
327
328
      }
    }

    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
329
	neighbour[2] = elInfoOld->neighbour[1];
330
	oppVertex[2] = elInfoOld->oppVertex[1];
331
	
Thomas Witkowski's avatar
Thomas Witkowski committed
332
333
	if (neighbour[2] && fill_opp_coords) {
	  oppCoord[2] = elInfoOld->oppCoord[1];
334
	  neighbourCoord[2] = elInfoOld->neighbourCoord[1];
335
	}
336
337
338
339
340
341
342
343
344
	
	
	// 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
345
	  neighbour[1] = elem->getSecondChild()->getSecondChild();
346
	  oppVertex[1] = 2;
347
348
	  
	  if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
349
350
351
352
353
354
355
356
            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];  
357
358
	  }
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
359
	  neighbour[1] = elem->getSecondChild();
360
	  oppVertex[1] = 0;
361
362

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

Thomas Witkowski's avatar
Thomas Witkowski committed
365
366
367
	    neighbourCoord[1][0] = elInfoOld->coord[1];
	    neighbourCoord[1][1] = elInfoOld->coord[2];
	    neighbourCoord[1][2] = coord[2];
368
369
370
371
	  }
	}


372
373
374
	// Calculation of the neighbour 0, its oppCoords and the
	// cooresponding oppVertex.
	
Thomas Witkowski's avatar
Thomas Witkowski committed
375
	nb = elInfoOld->neighbour[2];
376
	if (nb) {
377
	  TEST(elInfoOld->oppVertex[2] == 2)("invalid neighbour\n"); 
378
379
	  TEST_EXIT_DBG(nb->getFirstChild())("missing first child?\n");
	  TEST_EXIT_DBG(nb->getSecondChild())("missing second child?\n");
380
381
382
383
	 
	  nb = nb->getSecondChild();

	  if (nb->getFirstChild()) {
384
	    oppVertex[0] = 2;
385

386
	    if (fill_opp_coords) {
387
	      if (nb->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
388
		oppCoord[0] = *(nb->getNewCoord());
389
	      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
390
		oppCoord[0].setMidpoint(elInfoOld->neighbourCoord[2][1],
391
					 elInfoOld->neighbourCoord[2][2]);
392
	      }
393

394
395
396
	      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
397
	      neighbourCoord[0][2] = oppCoord[0];
398
399
400
401
	    }	   
 
	    nb = nb->getFirstChild();
	  } else {
402
	    oppVertex[0] = 1;
403

404
	    if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
405
	      oppCoord[0] = elInfoOld->oppCoord[2];    
406

407
408
409
410
	      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]);
411
412
	    }
	  }
413
414
	}
	
Thomas Witkowski's avatar
Thomas Witkowski committed
415
	neighbour[0] = nb;
416
417
418
419
      } else {   /* ichild == 1 */
	// Calculation of the neighbour 2, its oppCoords and the
	// cooresponding oppVertex.

Thomas Witkowski's avatar
Thomas Witkowski committed
420
	neighbour[2] = elInfoOld->neighbour[0];
421
	oppVertex[2] = elInfoOld->oppVertex[0];
422

Thomas Witkowski's avatar
Thomas Witkowski committed
423
424
	if (neighbour[2] && fill_opp_coords) {
	  oppCoord[2] = elInfoOld->oppCoord[0];
425
	  neighbourCoord[2] = elInfoOld->neighbourCoord[0];
426
427
428
429
430
431
432
	}
	

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

	if (elem->getFirstChild()->getFirstChild()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
433
	  neighbour[0] = elem->getFirstChild()->getFirstChild();
434
	  oppVertex[0] = 2;
435
436

	  if (fill_opp_coords) {
437
            if (elem->getFirstChild()->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
438
	      oppCoord[0] = *(elem->getFirstChild()->getNewCoord());
439
	    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
440
441
	      oppCoord[0].setMidpoint(elInfoOld->coord[0], 
				       elInfoOld->coord[2]);
442
	    }
443

Thomas Witkowski's avatar
Thomas Witkowski committed
444
445
446
	    neighbourCoord[0][0] = coord[2];
	    neighbourCoord[0][1] = coord[1];
	    neighbourCoord[0][2] = oppCoord[0];
447
	  }
448
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
449
	  neighbour[0] = elem->getFirstChild();
450
	  oppVertex[0] = 1;
451
452

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

Thomas Witkowski's avatar
Thomas Witkowski committed
455
456
457
	    neighbourCoord[0][0] = elInfoOld->coord[2];
	    neighbourCoord[0][1] = elInfoOld->coord[0];
	    neighbourCoord[0][2] = coord[2];
458
	  }
459
460
461
462
463
	}

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

Thomas Witkowski's avatar
Thomas Witkowski committed
464
	nb = elInfoOld->neighbour[2];
465
	if (nb) {
466
	  TEST(elInfoOld->oppVertex[2] == 2)("invalid neighbour\n"); 
467
468
469
	  TEST((nb = nb->getFirstChild()))("missing child?\n");

	  if (nb->getFirstChild()) {
470
	    oppVertex[1] = 2;
471

472
	    if (fill_opp_coords) {
473
	      if (nb->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
474
		oppCoord[1] = *(nb->getNewCoord());
475
	      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
476
		oppCoord[1].setMidpoint(elInfoOld->neighbourCoord[2][0],
477
					 elInfoOld->neighbourCoord[2][2]);
478
	      }
479

480
481
482
	      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
483
	      neighbourCoord[1][2] = oppCoord[1];
484
	    }
485
486
	    nb = nb->getSecondChild();

487
	  } else {
488
	    oppVertex[1] = 0;
489

490
	    if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
491
	      oppCoord[1] = elInfoOld->oppCoord[2];
492

493
494
495
496
	      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]);
497
498
499
	    }
	  }
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
500
	neighbour[1] = nb;
501
      } // if (ichild == 0) {} else
Thomas Witkowski's avatar
Thomas Witkowski committed
502
    } // if (fill_flag.isSet(Mesh::FILL_NEIGH) || fillFlag.isSet(Mesh::FILL_OPP_COORDS))
503
504
    

505
    if (fill_flag.isSet(Mesh::FILL_BOUND)) {
506
      if (elInfoOld->getBoundary(2))
Thomas Witkowski's avatar
Thomas Witkowski committed
507
	boundary[5] = elInfoOld->getBoundary(2);
508
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
509
	boundary[5] = INTERIOR;
510

511
      if (ichild == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
512
513
514
515
516
	boundary[3] = elInfoOld->getBoundary(5);
	boundary[4] = elInfoOld->getBoundary(3);
	boundary[0] = elInfoOld->getBoundary(2);
	boundary[1] = INTERIOR;
	boundary[2] = elInfoOld->getBoundary(1);
517
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
518
519
520
521
522
	boundary[3] = elInfoOld->getBoundary(4);
	boundary[4] = elInfoOld->getBoundary(5);
	boundary[0] = INTERIOR;
	boundary[1] = elInfoOld->getBoundary(2);
	boundary[2] = elInfoOld->getBoundary(0);
523
524
      }

525
526
      if (elInfoOld->getProjection(0) && 
	  elInfoOld->getProjection(0)->getType() == VOLUME_PROJECTION) {
527
	
Thomas Witkowski's avatar
Thomas Witkowski committed
528
	projection[0] = elInfoOld->getProjection(0);
529
530
      } else { // boundary projection
	if (ichild == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
531
532
533
	  projection[0] = elInfoOld->getProjection(2);
	  projection[1] = NULL;
	  projection[2] = elInfoOld->getProjection(1);
534
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
535
536
537
	  projection[0] = NULL;
	  projection[1] = elInfoOld->getProjection(2);
	  projection[2] = elInfoOld->getProjection(0);
538
539
540
541
542
	}
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
543
  double ElInfo2d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
544
  {
545
    FUNCNAME("ElInfo2d::calcGrdLambda()");
546
547

    testFlag(Mesh::FILL_COORDS);
548
549
  
    double adet = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
550
    int dim = mesh->getDim();
551

Thomas Witkowski's avatar
Thomas Witkowski committed
552
    for (int i = 0; i < dimOfWorld; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
553
554
      (*e1)[i] = coord[1][i] - coord[0][i];
      (*e2)[i] = coord[2][i] - coord[0][i];
555
556
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
557
    if (dimOfWorld == 2) {
558
      double sdet = (*e1)[0] * (*e2)[1] - (*e1)[1] * (*e2)[0];
559
560
561
562
      adet = abs(sdet);

      if (adet < 1.0E-25) {
	MSG("abs(det) = %f\n", adet);
563
	for (int i = 0; i <= dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
564
	  for (int j = 0; j < dimOfWorld; j++)
565
	    grd_lam[i][j] = 0.0;
566
567
      } else {
	double det1 = 1.0 / sdet;
568
569
570
571
572

	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
573
574
	grd_lam[0][0] = - grd_lam[1][0] - grd_lam[2][0];
	grd_lam[0][1] = - grd_lam[1][1] - grd_lam[2][1];
575
      }
576
577
    } else {  
      vectorProduct(*e1, *e2, *normal);
578

579
      adet = norm(normal);
580
581
582

      if (adet < 1.0E-15) {
	MSG("abs(det) = %lf\n", adet);
583
	for (int i = 0; i <= dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
584
	  for (int j = 0; j < dimOfWorld; j++)
585
586
	    grd_lam[i][j] = 0.0;
      } else {
587
588
	vectorProduct(*e2, *normal, grd_lam[1]);
	vectorProduct(*normal, *e1, grd_lam[2]);
589
      
590
	double adet2 = 1.0 / (adet * adet);
591

Thomas Witkowski's avatar
Thomas Witkowski committed
592
	for (int i = 0; i < dimOfWorld; i++) {
593
594
595
596
597
598
599
600
601
602
603
604
605
	  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;
  }

606
607

  const int ElInfo2d::worldToCoord(const WorldVector<double>& xy, 
608
609
				   DimVec<double>* lambda) const
  {
610
    FUNCNAME("ElInfo::worldToCoord()");
611

612
    TEST_EXIT_DBG(lambda)("lambda must not be NULL\n");
613

Thomas Witkowski's avatar
Thomas Witkowski committed
614
    DimVec<WorldVector<double> > edge(mesh->getDim(), NO_INIT);
615
    WorldVector<double> x; 
Thomas Witkowski's avatar
Thomas Witkowski committed
616
    static DimVec<double> vec(mesh->getDim(), NO_INIT);
617

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

620
    for (int j = 0; j < dimOfWorld; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
621
      double x0 = coord[dim][j];
622
      x[j] = xy[j] - x0;
623
      for (int i = 0; i < dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
624
	edge[i][j] = coord[i][j] - x0;
625
626
    }
  
627
628
629
    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]; 
630
631
632

    if (abs(det) < DBL_TOL) {
      ERROR("det = %le; abort\n", det);
633
634
      for (int i = 0; i <= dim; i++) 
	(*lambda)[i] = 1.0/dim;
635
636
637
638
639
640
641
642
643
      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;
644
    for (int i = 0; i <= dim; i++) {
645
646
647
648
649
650
651
652
653
654
655
656
      if ((*lambda)[i] < -1.E-5) {
	if ((*lambda)[i] < lmin) {
	  k = i;
	  lmin = (*lambda)[i];
	}
      }
    }

    return k;
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
657
  double ElInfo2d::getNormal(int side, WorldVector<double> &normal)
658
  {
659
    FUNCNAME("ElInfo::getNormal()");
660

661
662
    int i0 = (side + 1) % 3;
    int i1 = (side + 2) % 3;
663

Thomas Witkowski's avatar
Thomas Witkowski committed
664
    if (dimOfWorld == 2){
Thomas Witkowski's avatar
Thomas Witkowski committed
665
666
      normal[0] = coord[i1][1] - coord[i0][1];
      normal[1] = coord[i0][0] - coord[i1][0];
667
668
669
    } else { // dow == 3
      WorldVector<double> e0, e1,e2, elementNormal;

Thomas Witkowski's avatar
Thomas Witkowski committed
670
671
672
673
674
675
      e0 = coord[i1]; 
      e0 -= coord[i0];
      e1 = coord[i1]; 
      e1 -= coord[side];
      e2 = coord[i0]; 
      e2 -= coord[side];
676
677
678
679
680

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

681
    double det = norm(&normal);
682

683
    TEST_EXIT_DBG(det > 1.e-30)("det = 0 on face %d\n", side);
684
685
686

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

690

691
692
693
694
695
  /****************************************************************************/
  /*  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
696
  double ElInfo2d::getElementNormal(WorldVector<double> &elementNormal) const
697
  {
698
    FUNCNAME("ElInfo::getElementNormal()");
699

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

Thomas Witkowski's avatar
Thomas Witkowski committed
703
704
    WorldVector<double> e0 = coord[1] - coord[0];
    WorldVector<double> e1 = coord[2] - coord[0];
705
706
707

    vectorProduct(e0, e1, elementNormal);

708
    double det = norm(&elementNormal);
709

710
    TEST_EXIT_DBG(det > 1.e-30)("det = 0");
711
712
713

    elementNormal *= 1.0 / det;
    
Thomas Witkowski's avatar
Thomas Witkowski committed
714
    return det;
715
  }
716

717

718
  mtl::dense2D<double>& ElInfo2d::getSubElemCoordsMat(int degree) const
719
720
721
722
723
  {
    FUNCNAME("ElInfo2d::getSubElemCoordsMat()");

    using namespace mtl;

724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
    if (subElemMatrices[degree].count(refinementPath) == 0) {
      switch (degree) {
      case 1:
	{
	  dense2D<double> 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;
	    }
739
740
	  }

741
742
743
	  subElemMatrices[1][refinementPath] = mat;  
	}
	break;
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
      case 2:
	{
	  dense2D<double> mat(mat_d2);
	  dense2D<double> tmpMat(num_rows(mat), num_rows(mat));
	  
	  for (int i = 0; i < refinementPathLength; i++) {
	    if (refinementPath & (1 << i)) {
	      tmpMat = mat_d2_right * mat;
	      mat = tmpMat;
	    } else  {
	      tmpMat = mat_d2_left * mat;
	      mat = tmpMat;
	    }
	  }

	  subElemMatrices[2][refinementPath] = mat;  
	}
	break;
762
763
      default:
	ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
764
765
      }
    }
766
767

    return subElemMatrices[degree][refinementPath];
768
769
770
  }


771
  mtl::dense2D<double>& ElInfo2d::getSubElemGradCoordsMat(int degree) const
772
  {
773
    FUNCNAME("ElInfo2d::getSubElemGradCoordsMat()");
774

775
776
777
778
    TEST_EXIT(degree == 1)("Not supported for basis functions with degree > 1!\n");

    using namespace mtl;

779
    if (subElemGradMatrices[degree].count(refinementPath) == 0) {
780
      dense2D<double> mat(mat_d1);
781
782
783
784
785
786
787
788
789
790
791
792
      dense2D<double> tmpMat(num_rows(mat), num_rows(mat));

      double test_left[3][3] = {{0.0, 0.0, 0.5},
				{-0.5, -0.5, 0.0},
				{1.0, 0.0, 0.0}};
      double test_right[3][3] = {{0.0, 0.0, 0.5},
				 {0.5, -0.5, 0.0},
				 {0.0, 1.0, 0.0}};
      
      mtl::dense2D<double> mat_left(test_left);
      mtl::dense2D<double> mat_right(test_right);

793
794

      for (int i = 0; i < refinementPathLength; i++)
795
796
797
798
799
800
801
	if (refinementPath & (1 << i)) {
	  tmpMat = mat_right * mat;
	  mat = tmpMat;
	} else  {
	  tmpMat = mat_left * mat;
	  mat = tmpMat;
	}
802
803
804
805
806

      subElemGradMatrices[1][refinementPath] = mat;
    }

    return subElemGradMatrices[degree][refinementPath];
807

808
  }
809
}