ElInfo2d.cc 27.4 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
  double ElInfo2d::mat_d1_left_val[3][3] = {{0.0, 1.0, 0.5}, 
18
  					    {0.0, 0.0, 0.5},
19
  					    {1.0, 0.0, 0.0}}; 
20
21
  mtl::dense2D<double> ElInfo2d::mat_d1_left(mat_d1_left_val);

22
  
23
  double ElInfo2d::mat_d1_right_val[3][3] = {{0.0, 0.0, 0.5}, 
24
  					     {1.0, 0.0, 0.5},
25
  					     {0.0, 1.0, 0.0}}; 
26
27
  mtl::dense2D<double> ElInfo2d::mat_d1_right(mat_d1_right_val);

28

29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46

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


47

48
49
  double ElInfo2d::mat_d3_left_val[10][10] = {{0.0, 1.0, -0.0625, 0.3125, 0.0, 0.0, 0.0625, 0.0, 0.0, -0.0625},
					      {0.0, 0.0, -0.0625, 0.0625, 0.0, 0.0, 0.0625, 0.0, 0.0, 0.0625},
50
					      {1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
51
					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.0, 0.0, -0.125},
52
53
					      {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.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0},
54
55
56
					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.0, 1.0, 0.375},
					      {0.0, 0.0, 0.5625, 0.9375, 1.0, 0.0, -0.0625, 0.0, 0.0, 0.1875},
					      {0.0, 0.0, 0.5625, -0.3125, 0.0, 0.0, -0.0625, 0.0, 0.0, 0.0},
57
58
59
					      {0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.75}};
  mtl::dense2D<double> ElInfo2d::mat_d3_left(mat_d3_left_val);

60
61
  double ElInfo2d::mat_d3_right_val[10][10] = {{0.0, 0.0, -0.0625, 0.0625, 0.0, 0.0, 0.0625, 0.0, 0.0, 0.0625},
					       {1.0, 0.0, -0.0625, 0.0625, 0.0, 0.0, 0.3125, 0.0, 0.0, -0.0625},
62
					       {0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
63
					       {0.0, 0.0, 0.0, -0.25, 0.0, 0.0, 0.0, 1.0, 0.0, 0.375},
64
65
					       {0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
					       {0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
66
67
68
					       {0.0, 0.0, 0.0, -0.25, 0.0, 0.0, 0.0, 0.0, 0.0, -0.125},
					       {0.0, 0.0, 0.5625, -0.0625, 0.0, 0.0, -0.3125, 0.0, 0.0, -0.1875},
					       {0.0, 0.0, 0.5625, -0.0625, 0.0, 1.0, 0.9375, 0.0, 0.0, 0.1875},
69
70
71
72
73
					       {0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.0, 0.75}};
  mtl::dense2D<double> ElInfo2d::mat_d3_right(mat_d3_right_val);



74
75
  double ElInfo2d::mat_d4_left_val[15][15] = {{0.0, 1.0, 0.0, 2.734375e-01, 0.0, -3.90625e-02, 2.34375e-02, 0.0, -3.90625e-02, 0.0, 0.0, 0.0, 2.34375e-02, -3.90625e-02, 0.0},
					      {0.0, 0.0, 0.0, -3.90625e-02, 0.0, 2.34375e-02, 2.34375e-02, 0.0, -3.90625e-02, 0.0, 0.0, 0.0, -3.90625e-02, -3.90625e-02, 0.0},
76
					      {1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
77
78
79
80
81
82
83
84
85
86
87
88
					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -6.25e-02, 0.0, 1.875e-01, 0.0, 0.0, 0.0, 1.25e-01, 6.25e-02, 0.0},
					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -3.75e-01, 0.0, 0.0, 0.0, -1.25e-01, 0.0, 0.0},
					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.e-01, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.e-01, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -3.75e-01, 0.0, 1.0, 0.0, 3.75e-01, 0.0, 0.0},
					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -6.25e-02, 0.0, 1.875e-01, 0.0, 0.0, 1.0, -1.25e-01, 3.125e-01, 0.0},
					      {0.0, 0.0, 0.0, 1.09375e+00, 1.0, 4.6875e-01, -9.375e-02, 0.0, 3.125e-02, 0.0, 0.0, 0.0, -3.125e-02, 1.5625e-01, 0.0},
					      {0.0, 0.0, 1.0, -5.46875e-01, 0.0, 7.03125e-01, 1.40625e-01, 0.0, 1.5625e-02, 0.0, 0.0, 0.0, -4.6875e-02, -2.34375e-01, 0.0},
					      {0.0, 0.0, 0.0, 2.1875e-01, 0.0, -1.5625e-01, -9.375e-02, 0.0, 3.125e-02, 0.0, 0.0, 0.0, 9.375e-02, 1.5625e-01, 0.0},
					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.625e-01, 0.0, -1.875e-01, 0.0, 0.0, 0.0, 3.75e-01, 9.375e-01, 1.0},
					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.625e-01, 0.0, -1.875e-01, 0.0, 0.0, 0.0, -3.75e-01, -3.125e-01, 0.0},
					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 7.5e-01, 0.0, 0.0, 0.0, 7.5e-01, 0.0, 0.0}};
89
90
  mtl::dense2D<double> ElInfo2d::mat_d4_left(mat_d4_left_val);

91
92
  double ElInfo2d::mat_d4_right_val[15][15] = {{0.0, 0.0, 0.0, -3.90625e-02, 0.0, 2.34375e-02, 2.34375e-02, 0.0, -3.90625e-02, 0.0, 0.0, 0.0, -3.90625e-02, -3.90625e-02, 0.0},
					       {1.0, 0.0, 0.0, -3.90625e-02, 0.0, 2.34375e-02, -3.90625e-02, 0.0, 2.734375e-01, 0.0, 0.0, 0.0, -3.90625e-02, 2.34375e-02, 0.0},
93
					       {0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
94
95
96
97
98
99
100
101
102
103
104
105
					       {0.0, 0.0, 0.0, 1.875e-01, 0.0, -6.25e-02, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 3.125e-01, -1.25e-01, 0.0},
					       {0.0, 0.0, 0.0, -3.75e-01, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 3.75e-01, 0.0},
					       {0.0, 0.0, 0.0, 5.0e-01, 0.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, 5.0e-01, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
					       {0.0, 0.0, 0.0, -3.75e-01, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.25e-01, 0.0},
					       {0.0, 0.0, 0.0, 1.875e-01, 0.0, -6.25e-02, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 6.25e-02, 1.25e-01, 0.0},
					       {0.0, 0.0, 0.0, 3.125e-02, 0.0, -9.375e-02, -1.5625e-01, 0.0, 2.1875e-01, 0.0, 0.0, 0.0, 1.5625e-01, 9.375e-02, 0.0},
					       {0.0, 0.0, 1.0, 1.5625e-02, 0.0, 1.40625e-01, 7.03125e-01, 0.0, -5.46875e-01, 0.0, 0.0, 0.0, -2.34375e-01, -4.6875e-02, 0.0},
					       {0.0, 0.0, 0.0, 3.125e-02, 0.0, -9.375e-02, 4.6875e-01, 1.0, 1.09375e+00, 0.0, 0.0, 0.0, 1.5625e-01, -3.125e-02, 0.0},
					       {0.0, 0.0, 0.0, -1.875e-01, 0.0, 5.625e-01, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -3.125e-01, -3.75e-01, 0.0},
					       {0.0, 0.0, 0.0, -1.875e-01, 0.0, 5.625e-01, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 9.375e-01, 3.75e-01, 1.0},
					       {0.0, 0.0, 0.0, 7.5e-01, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 7.5e-01, 0.0}};
106
107
108
  mtl::dense2D<double> ElInfo2d::mat_d4_right(mat_d4_right_val);


109
110
  ElInfo2d::ElInfo2d(Mesh *aMesh) 
    : ElInfo(aMesh) 
111
  {}
112

113

114
  ElInfo2d::~ElInfo2d()
115
  {}
116

117

118
119
  void ElInfo2d::fillMacroInfo(const MacroElement * mel)
  {
120
121
    FUNCNAME("ElInfo::fillMacroInfo()");
 
Thomas Witkowski's avatar
Thomas Witkowski committed
122
123
124
    macroElement = const_cast<MacroElement*>(mel);
    element = const_cast<Element*>(mel->getElement());
    parent = NULL;
125
    level = 0;
126

Thomas Witkowski's avatar
Thomas Witkowski committed
127
128
129
    if (fillFlag.isSet(Mesh::FILL_COORDS) || 
	fillFlag.isSet(Mesh::FILL_DET)    ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
130

Thomas Witkowski's avatar
Thomas Witkowski committed
131
      int vertices = mesh->getGeo(VERTEX);
132
      for (int i = 0; i < vertices; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
133
	coord[i] = mel->coord[i];
134
135
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
138
139
    if (fillFlag.isSet(Mesh::FILL_OPP_COORDS) || 
	fillFlag.isSet(Mesh::FILL_NEIGH)) {
140

Thomas Witkowski's avatar
Thomas Witkowski committed
141
      bool fill_opp_coords = (fillFlag.isSet(Mesh::FILL_OPP_COORDS));
142
    
143
144
145
146
      for (int i = 0; i < neighbours; i++) {
	MacroElement *macroNeighbour = mel->getNeighbour(i);

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

150
	  int edgeNo = oppVertex[i] = mel->getOppVertex(i);
151
152


153
	  if (nb->getFirstChild() && edgeNo != 2) {  
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173

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


174
	    if (edgeNo == 0) {
175
176
177
178
179
180
181
182
183
184
185
186
	      // The situation is as follows:
	      //
	      //          -------
	      //          \    /|\
	      //           \  / | \
	      //            \/  |  \
	      //             \  |   \
	      //              \ |    \ 
	      //               \|     \
	      //                -------
	      //
	      //            nb     el
187
              //
188
189
190
191
	      // 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
192
	      nb = neighbour[i] = nb->getSecondChild();
193
	    } else {
194
195
	      // 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
196
	      nb = neighbour[i] = nb->getFirstChild();
197
198
	    }

199
200
	    // In both cases the opp vertex number is 2, as one can see in the 
	    // pictures above.
201
	    oppVertex[i] = 2;
202
203

	    if (fill_opp_coords) {
204
	      if (nb->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
205
		oppCoord[i] = *(nb->getNewCoord());
206
	      } else {
207
208
209
		// 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
210
		oppCoord[i] = (macroNeighbour->coord[0] + 
211
212
				macroNeighbour->coord[1]) * 0.5;
	      }
213
214
215
	      
	      switch (i) {
	      case 0:
216
217
218
219
		// The common edge is the edge 0 of this element.

		switch (edgeNo) {
		case 1:
220
221
		  neighbourCoord[i][0] = macroNeighbour->coord[2];
		  neighbourCoord[i][1] = macroNeighbour->coord[0];
222
223
		  break;
		case 0:		  
224
225
		  neighbourCoord[i][0] = macroNeighbour->coord[1];
		  neighbourCoord[i][1] = macroNeighbour->coord[2];
226
227
		  break;
		default:
Thomas Witkowski's avatar
Thomas Witkowski committed
228
		  ERROR_EXIT("Should not happen!\n");
229
		}
230
	
Thomas Witkowski's avatar
Thomas Witkowski committed
231
		neighbourCoord[i][2] = oppCoord[i];
232
233
234
		break;
		
	      case 1:
235
236
237
		// The commonedge is the edge 1 of this element.
		switch (edgeNo) {
		case 0:
238
239
		  neighbourCoord[i][0] = macroNeighbour->coord[1];
		  neighbourCoord[i][1] = macroNeighbour->coord[2];
240
241
		  break;
		case 1:
242
243
		  neighbourCoord[i][0] = macroNeighbour->coord[2];
		  neighbourCoord[i][1] = macroNeighbour->coord[0];
244
245
		  break;
		default:
Thomas Witkowski's avatar
Thomas Witkowski committed
246
		  ERROR_EXIT("Should not happen!\n");
247
248
		}
		
Thomas Witkowski's avatar
Thomas Witkowski committed
249
		neighbourCoord[i][2] = oppCoord[i];
250
251
		break;
		
Thomas Witkowski's avatar
Thomas Witkowski committed
252
	      case 2:
253
		if (*(macroNeighbour->getElement()->getDof(2)) == *(element->getDof(0))) {
254
255
		  neighbourCoord[i][0] = macroNeighbour->coord[2];
		  neighbourCoord[i][1] = macroNeighbour->coord[1];
256
		} else if (*(macroNeighbour->getElement()->getDof(2)) == *(element->getDof(1))) {
257
258
		  neighbourCoord[i][0] = macroNeighbour->coord[0];
		  neighbourCoord[i][1] = macroNeighbour->coord[2];		 
259
260
261
262
		} else {
		  ERROR_EXIT("Should not happen!\n");
		}

263
264
265
		// 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.
266
267
		//		ERROR_EXIT("Should not happen!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
268
269
		break;

270
	      default:
Thomas Witkowski's avatar
Thomas Witkowski committed
271
		std::cout << "------------- Error --------------" << std::endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
272
		std::cout << "  Neighbour counter = " << i << "\n";
Thomas Witkowski's avatar
Thomas Witkowski committed
273
		std::cout << "  Element index     = " << element->getIndex() << "\n\n";
Thomas Witkowski's avatar
Thomas Witkowski committed
274
275
276
277
278
279
280
281
		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;
		  }
282
283
284
		  std::cout << "  OppVertex " << j << ": " 
			    << static_cast<int>(mel->getOppVertex(j)) 
			    << std::endl << std::endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
285
		}
286
287
		ERROR_EXIT("should not happen!\n");
		break;
288
289
290
	      }
	    }
	  } else {
291
292
293
294
295
296
297

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

298
	    if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
299
	      oppCoord[i] = macroNeighbour->coord[edgeNo];
300
	      neighbourCoord[i] = macroNeighbour->coord;	      
301
	    }
302
	  }
303
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
304
	  neighbour[i] = NULL;
305
        }
306
      }
307
    }
308
    
Thomas Witkowski's avatar
Thomas Witkowski committed
309
310
311
312
313
    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);      
314
315
316
317
318
319
320
321
    }
  }


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

322
  void ElInfo2d::fillElInfo(int ichild, const ElInfo *elInfoOld)
323
  {
324
    FUNCNAME("ElInfo::fillElInfo()");
325

Thomas Witkowski's avatar
Thomas Witkowski committed
326
327
    Element *elem = elInfoOld->element;
    Flag fill_flag = elInfoOld->fillFlag;
328

329
    TEST_EXIT_DBG(elem->getFirstChild())("no children?\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
330
    element = const_cast<Element*>((ichild == 0) ? 
331
332
				    elem->getFirstChild() : 
				    elem->getSecondChild());
Thomas Witkowski's avatar
Thomas Witkowski committed
333
    TEST_EXIT_DBG(element)("missing child %d?\n", ichild);
334

Thomas Witkowski's avatar
Thomas Witkowski committed
335
336
337
    macroElement  = elInfoOld->macroElement;
    fillFlag = fill_flag;
    parent = elem;
338
    level = elInfoOld->level + 1;
339
    iChild = ichild;
340

Thomas Witkowski's avatar
Thomas Witkowski committed
341
342
343
    if (fillFlag.isSet(Mesh::FILL_COORDS) || 
	fillFlag.isSet(Mesh::FILL_DET)    ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
344
      
345
      if (elem->isNewCoordSet())
Thomas Witkowski's avatar
Thomas Witkowski committed
346
	coord[2] = *(elem->getNewCoord());
347
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
348
	coord[2].setMidpoint(elInfoOld->coord[0], elInfoOld->coord[1]);      
349
350
      
      if (ichild == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
351
352
	coord[0] = elInfoOld->coord[2];
	coord[1] = elInfoOld->coord[0];
353
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
354
355
	coord[0] = elInfoOld->coord[1];
	coord[1] = elInfoOld->coord[2];
356
357
358
359
360
361
362
363
364
365
      }
    }

    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
366
	neighbour[2] = elInfoOld->neighbour[1];
367
	oppVertex[2] = elInfoOld->oppVertex[1];
368
	
Thomas Witkowski's avatar
Thomas Witkowski committed
369
370
	if (neighbour[2] && fill_opp_coords) {
	  oppCoord[2] = elInfoOld->oppCoord[1];
371
	  neighbourCoord[2] = elInfoOld->neighbourCoord[1];
372
	}
373
374
375
376
377
378
379
380
381
	
	
	// 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
382
	  neighbour[1] = elem->getSecondChild()->getSecondChild();
383
	  oppVertex[1] = 2;
384
385
	  
	  if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
386
387
388
389
390
391
392
393
            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];  
394
395
	  }
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
396
	  neighbour[1] = elem->getSecondChild();
397
	  oppVertex[1] = 0;
398
399

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

Thomas Witkowski's avatar
Thomas Witkowski committed
402
403
404
	    neighbourCoord[1][0] = elInfoOld->coord[1];
	    neighbourCoord[1][1] = elInfoOld->coord[2];
	    neighbourCoord[1][2] = coord[2];
405
406
407
408
	  }
	}


409
410
411
	// Calculation of the neighbour 0, its oppCoords and the
	// cooresponding oppVertex.
	
412
	Element *nb = elInfoOld->neighbour[2];
413
	if (nb) {
414
415
416
417
418
419
420
421
422
423
424
	  TEST_EXIT_DBG(elInfoOld->oppVertex[2] == 2)
	    ("Fill child %d of element %d (mel %d): Invalid neighbour %d!\n",
	     ichild,
	     elInfoOld->getElement()->getIndex(), 
	     elInfoOld->getMacroElement()->getIndex(),
	     nb->getIndex());

	  TEST_EXIT_DBG(nb->getFirstChild())
	    ("Missing first child in element %d!\n", nb->getIndex());	  
	  TEST_EXIT_DBG(nb->getSecondChild())
	    ("Missing second child in element %d!\n", nb->getIndex());
425
426
427
428
	 
	  nb = nb->getSecondChild();

	  if (nb->getFirstChild()) {
429
	    oppVertex[0] = 2;
430

431
	    if (fill_opp_coords) {
432
	      if (nb->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
433
		oppCoord[0] = *(nb->getNewCoord());
434
	      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
435
		oppCoord[0].setMidpoint(elInfoOld->neighbourCoord[2][1],
436
					 elInfoOld->neighbourCoord[2][2]);
437
	      }
438

439
440
441
	      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
442
	      neighbourCoord[0][2] = oppCoord[0];
443
444
445
446
	    }	   
 
	    nb = nb->getFirstChild();
	  } else {
447
	    oppVertex[0] = 1;
448

449
	    if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
450
	      oppCoord[0] = elInfoOld->oppCoord[2];    
451

452
453
454
455
	      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]);
456
457
	    }
	  }
458
459
	}
	
Thomas Witkowski's avatar
Thomas Witkowski committed
460
	neighbour[0] = nb;
461
462
463
464
      } else {   /* ichild == 1 */
	// Calculation of the neighbour 2, its oppCoords and the
	// cooresponding oppVertex.

Thomas Witkowski's avatar
Thomas Witkowski committed
465
	neighbour[2] = elInfoOld->neighbour[0];
466
	oppVertex[2] = elInfoOld->oppVertex[0];
467

Thomas Witkowski's avatar
Thomas Witkowski committed
468
469
	if (neighbour[2] && fill_opp_coords) {
	  oppCoord[2] = elInfoOld->oppCoord[0];
470
	  neighbourCoord[2] = elInfoOld->neighbourCoord[0];
471
472
473
474
475
476
477
	}
	

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

	if (elem->getFirstChild()->getFirstChild()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
478
	  neighbour[0] = elem->getFirstChild()->getFirstChild();
479
	  oppVertex[0] = 2;
480
481

	  if (fill_opp_coords) {
482
            if (elem->getFirstChild()->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
483
	      oppCoord[0] = *(elem->getFirstChild()->getNewCoord());
484
	    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
485
486
	      oppCoord[0].setMidpoint(elInfoOld->coord[0], 
				       elInfoOld->coord[2]);
487
	    }
488

Thomas Witkowski's avatar
Thomas Witkowski committed
489
490
491
	    neighbourCoord[0][0] = coord[2];
	    neighbourCoord[0][1] = coord[1];
	    neighbourCoord[0][2] = oppCoord[0];
492
	  }
493
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
494
	  neighbour[0] = elem->getFirstChild();
495
	  oppVertex[0] = 1;
496
497

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

Thomas Witkowski's avatar
Thomas Witkowski committed
500
501
502
	    neighbourCoord[0][0] = elInfoOld->coord[2];
	    neighbourCoord[0][1] = elInfoOld->coord[0];
	    neighbourCoord[0][2] = coord[2];
503
	  }
504
505
506
507
508
	}

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

509
	Element *nb = elInfoOld->neighbour[2];
510
	if (nb) {
511
	  TEST(elInfoOld->oppVertex[2] == 2)("invalid neighbour\n"); 
512
513
514
	  TEST((nb = nb->getFirstChild()))("missing child?\n");

	  if (nb->getFirstChild()) {
515
	    oppVertex[1] = 2;
516

517
	    if (fill_opp_coords) {
518
	      if (nb->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
519
		oppCoord[1] = *(nb->getNewCoord());
520
	      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
521
		oppCoord[1].setMidpoint(elInfoOld->neighbourCoord[2][0],
522
					 elInfoOld->neighbourCoord[2][2]);
523
	      }
524

525
526
527
	      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
528
	      neighbourCoord[1][2] = oppCoord[1];
529
	    }
530
531
	    nb = nb->getSecondChild();

532
	  } else {
533
	    oppVertex[1] = 0;
534

535
	    if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
536
	      oppCoord[1] = elInfoOld->oppCoord[2];
537

538
539
540
541
	      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]);
542
543
544
	    }
	  }
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
545
	neighbour[1] = nb;
546
      } // if (ichild == 0) {} else
Thomas Witkowski's avatar
Thomas Witkowski committed
547
    } // if (fill_flag.isSet(Mesh::FILL_NEIGH) || fillFlag.isSet(Mesh::FILL_OPP_COORDS))
548
549
    

550
    if (fill_flag.isSet(Mesh::FILL_BOUND)) {
551
      if (elInfoOld->getBoundary(2))
Thomas Witkowski's avatar
Thomas Witkowski committed
552
	boundary[5] = elInfoOld->getBoundary(2);
553
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
554
	boundary[5] = INTERIOR;
555

556
      if (ichild == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
557
558
559
560
561
	boundary[3] = elInfoOld->getBoundary(5);
	boundary[4] = elInfoOld->getBoundary(3);
	boundary[0] = elInfoOld->getBoundary(2);
	boundary[1] = INTERIOR;
	boundary[2] = elInfoOld->getBoundary(1);
562
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
563
564
565
566
567
	boundary[3] = elInfoOld->getBoundary(4);
	boundary[4] = elInfoOld->getBoundary(5);
	boundary[0] = INTERIOR;
	boundary[1] = elInfoOld->getBoundary(2);
	boundary[2] = elInfoOld->getBoundary(0);
568
569
      }

570
571
      if (elInfoOld->getProjection(0) && 
	  elInfoOld->getProjection(0)->getType() == VOLUME_PROJECTION) {
572
	
Thomas Witkowski's avatar
Thomas Witkowski committed
573
	projection[0] = elInfoOld->getProjection(0);
574
575
      } else { // boundary projection
	if (ichild == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
576
577
578
	  projection[0] = elInfoOld->getProjection(2);
	  projection[1] = NULL;
	  projection[2] = elInfoOld->getProjection(1);
579
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
580
581
582
	  projection[0] = NULL;
	  projection[1] = elInfoOld->getProjection(2);
	  projection[2] = elInfoOld->getProjection(0);
583
584
585
586
587
	}
      }
    }
  }

588
  double ElInfo2d::calcGrdLambda(DimVec<WorldVector<double> >& grd)
589
  {
590
    FUNCNAME("ElInfo2d::calcGrdLambda()");
591
592

    testFlag(Mesh::FILL_COORDS);
593
594
  
    double adet = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
595
    int dim = mesh->getDim();
596

Thomas Witkowski's avatar
Thomas Witkowski committed
597
    for (int i = 0; i < dimOfWorld; i++) {
598
599
      e1[i] = coord[1][i] - coord[0][i];
      e2[i] = coord[2][i] - coord[0][i];
600
601
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
602
    if (dimOfWorld == 2) {
603
      double sdet = e1[0] * e2[1] - e1[1] * e2[0];
604
605
606
607
      adet = abs(sdet);

      if (adet < 1.0E-25) {
	MSG("abs(det) = %f\n", adet);
608
	for (int i = 0; i <= dim; i++)
609
	  grd[i].set(0.0);
610
611
      } else {
	double det1 = 1.0 / sdet;
612

613
614
615
616
617
618
	grd[1][0] = e2[1] * det1;  // a11: (a_ij) = A^{-T}
	grd[1][1] = -e2[0] * det1; // a21
	grd[2][0] = -e1[1] * det1; // a12
	grd[2][1] = e1[0] * det1;  // a22
	grd[0][0] = -grd[1][0] - grd[2][0];
	grd[0][1] = -grd[1][1] - grd[2][1];
619
      }
620
    } else {  
621
      vectorProduct(e1, e2, normal);
622

623
      adet = norm(normal);
624
625
626

      if (adet < 1.0E-15) {
	MSG("abs(det) = %lf\n", adet);
627
	for (int i = 0; i <= dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
628
	  for (int j = 0; j < dimOfWorld; j++)
629
	    grd[i][j] = 0.0;
630
      } else {
631
632
	vectorProduct(e2, normal, grd[1]);
	vectorProduct(normal, e1, grd[2]);
633
      
634
	double adet2 = 1.0 / (adet * adet);
635

Thomas Witkowski's avatar
Thomas Witkowski committed
636
	for (int i = 0; i < dimOfWorld; i++) {
637
638
	  grd[1][i] *= adet2;
	  grd[2][i] *= adet2;
639
640
	}

641
642
643
	grd[0][0] = -grd[1][0] - grd[2][0];
	grd[0][1] = -grd[1][1] - grd[2][1];
	grd[0][2] = -grd[1][2] - grd[2][2];
644
645
646
647
648
649
      }
    }

    return adet;
  }

650
651

  const int ElInfo2d::worldToCoord(const WorldVector<double>& xy, 
652
653
				   DimVec<double>* lambda) const
  {
654
    FUNCNAME("ElInfo::worldToCoord()");
655

656
    TEST_EXIT_DBG(lambda)("lambda must not be NULL\n");
657

Thomas Witkowski's avatar
Thomas Witkowski committed
658
    DimVec<WorldVector<double> > edge(mesh->getDim(), NO_INIT);
659
    WorldVector<double> x; 
Thomas Witkowski's avatar
Thomas Witkowski committed
660
    static DimVec<double> vec(mesh->getDim(), NO_INIT);
661

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

664
    for (int j = 0; j < dimOfWorld; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
665
      double x0 = coord[dim][j];
666
      x[j] = xy[j] - x0;
667
      for (int i = 0; i < dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
668
	edge[i][j] = coord[i][j] - x0;
669
670
    }
  
671
672
673
    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]; 
674
675
676

    if (abs(det) < DBL_TOL) {
      ERROR("det = %le; abort\n", det);
677
678
      for (int i = 0; i <= dim; i++) 
	(*lambda)[i] = 1.0/dim;
679
680
681
682
683
684
685
686
687
      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;
688
    for (int i = 0; i <= dim; i++) {
689
690
691
692
693
694
695
696
697
698
699
700
      if ((*lambda)[i] < -1.E-5) {
	if ((*lambda)[i] < lmin) {
	  k = i;
	  lmin = (*lambda)[i];
	}
      }
    }

    return k;
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
701
  double ElInfo2d::getNormal(int side, WorldVector<double> &normal)
702
  {
703
    FUNCNAME("ElInfo::getNormal()");
704

705
706
    int i0 = (side + 1) % 3;
    int i1 = (side + 2) % 3;
707

Thomas Witkowski's avatar
Thomas Witkowski committed
708
    if (dimOfWorld == 2){
Thomas Witkowski's avatar
Thomas Witkowski committed
709
710
      normal[0] = coord[i1][1] - coord[i0][1];
      normal[1] = coord[i0][0] - coord[i1][0];
711
    } else { // dow == 3
712
      WorldVector<double> e0, elementNormal;
713

Thomas Witkowski's avatar
Thomas Witkowski committed
714
715
716
717
718
719
      e0 = coord[i1]; 
      e0 -= coord[i0];
      e1 = coord[i1]; 
      e1 -= coord[side];
      e2 = coord[i0]; 
      e2 -= coord[side];
720
721
722
723
724

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

725
    double det = norm(&normal);
726

727
    TEST_EXIT_DBG(det > 1.e-30)("det = 0 on face %d\n", side);
728
729
730

    normal *= 1.0 / det;
    
Thomas Witkowski's avatar
Thomas Witkowski committed
731
    return det;
732
733
  }

734

735
736
737
738
739
  /****************************************************************************/
  /*  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
740
  double ElInfo2d::getElementNormal(WorldVector<double> &elementNormal) const
741
  {
742
    FUNCNAME("ElInfo::getElementNormal()");
743

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

Thomas Witkowski's avatar
Thomas Witkowski committed
747
748
    WorldVector<double> e0 = coord[1] - coord[0];
    WorldVector<double> e1 = coord[2] - coord[0];
749
750
751

    vectorProduct(e0, e1, elementNormal);

752
    double det = norm(&elementNormal);
753

754
    TEST_EXIT_DBG(det > 1.e-30)("det = 0");
755
756
757

    elementNormal *= 1.0 / det;
    
Thomas Witkowski's avatar
Thomas Witkowski committed
758
    return det;
759
  }
760

761

762
  mtl::dense2D<double>& ElInfo2d::getSubElemCoordsMat(int degree) const
763
764
765
766
767
  {
    FUNCNAME("ElInfo2d::getSubElemCoordsMat()");

    using namespace mtl;

768
769
770
771
    if (subElemMatrices[degree].count(refinementPath) == 0) {
      switch (degree) {
      case 1:
	{
772
773
774
	  dense2D<double> mat(3, 3), tmpMat(3, 3);
	  mat = 1;

775
776
777
778
779
780
781
782
	  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;
	    }
783
784
	  }

785
786
787
	  subElemMatrices[1][refinementPath] = mat;  
	}
	break;
788
789
      case 2:
	{
790
791
792
	  dense2D<double> mat(6, 6), tmpMat(6, 6);
	  mat = 1;

793
794
795
796
797
798
799
800
801
802
803
804
805
	  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;
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
      case 3:
	{
	  dense2D<double> mat(10, 10), tmpMat(10, 10);
	  mat = 1;

	  for (int i = 0; i < refinementPathLength; i++) {
	    if (refinementPath & (1 << i)) {
	      tmpMat = mat_d3_right * mat;
	      mat = tmpMat;
	    } else  {
	      tmpMat = mat_d3_left * mat;
	      mat = tmpMat;
	    }
	  }

	  subElemMatrices[3][refinementPath] = mat;  
	}
	break;
      case 4:
	{
	  dense2D<double> mat(15, 15), tmpMat(15, 15);
	  mat = 1;

	  for (int i = 0; i < refinementPathLength; i++) {
	    if (refinementPath & (1 << i)) {
	      tmpMat = mat_d4_right * mat;
	      mat = tmpMat;
	    } else  {
	      tmpMat = mat_d4_left * mat;
	      mat = tmpMat;
	    }
	  }

	  subElemMatrices[4][refinementPath] = mat;  
	}
	break;

843
844
      default:
	ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
845
846
      }
    }
847
848

    return subElemMatrices[degree][refinementPath];
849
850
851
  }


852
  mtl::dense2D<double>& ElInfo2d::getSubElemGradCoordsMat(int degree) const
853
  {
854
    FUNCNAME("ElInfo2d::getSubElemGradCoordsMat()");
855

856
857
858
859
    TEST_EXIT(degree == 1)("Not supported for basis functions with degree > 1!\n");

    using namespace mtl;

860
    if (subElemGradMatrices[degree].count(refinementPath) == 0) {
861
862
      dense2D<double> mat(3, 3), tmpMat(3, 3);
      mat = 1;
863
864
865
866
867
868
869
870
871
872
873

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

874
875

      for (int i = 0; i < refinementPathLength; i++)
876
877
878
879
880
881
882
	if (refinementPath & (1 << i)) {
	  tmpMat = mat_right * mat;
	  mat = tmpMat;
	} else  {
	  tmpMat = mat_left * mat;
	  mat = tmpMat;
	}
883
884
885
886
887

      subElemGradMatrices[1][refinementPath] = mat;
    }

    return subElemGradMatrices[degree][refinementPath];
888

889
  }
890
}