ElInfo2d.cc 27.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#include "ElInfo2d.h"
#include "BasisFunction.h"
#include "Element.h"
#include "Line.h"
#include "Triangle.h"
#include "Tetrahedron.h"
#include "FiniteElemSpace.h"
#include "Flag.h"
#include "MacroElement.h"
#include "Mesh.h"
#include "Global.h"
#include "FixVec.h"
#include "DOFVector.h"

namespace AMDiS {

17
  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:
Thomas Witkowski's avatar
Thomas Witkowski committed
253
		if (*(macroNeighbour->getElement()->getDOF(2)) == *(element->getDOF(0))) {
254
255
		  neighbourCoord[i][0] = macroNeighbour->coord[2];
		  neighbourCoord[i][1] = macroNeighbour->coord[1];
Thomas Witkowski's avatar
Thomas Witkowski committed
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
	  TEST(elInfoOld->oppVertex[2] == 2)("invalid neighbour\n"); 
415
416
	  TEST_EXIT_DBG(nb->getFirstChild())("missing first child?\n");
	  TEST_EXIT_DBG(nb->getSecondChild())("missing second child?\n");
417
418
419
420
	 
	  nb = nb->getSecondChild();

	  if (nb->getFirstChild()) {
421
	    oppVertex[0] = 2;
422

423
	    if (fill_opp_coords) {
424
	      if (nb->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
425
		oppCoord[0] = *(nb->getNewCoord());
426
	      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
427
		oppCoord[0].setMidpoint(elInfoOld->neighbourCoord[2][1],
428
					 elInfoOld->neighbourCoord[2][2]);
429
	      }
430

431
432
433
	      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
434
	      neighbourCoord[0][2] = oppCoord[0];
435
436
437
438
	    }	   
 
	    nb = nb->getFirstChild();
	  } else {
439
	    oppVertex[0] = 1;
440

441
	    if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
442
	      oppCoord[0] = elInfoOld->oppCoord[2];    
443

444
445
446
447
	      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]);
448
449
	    }
	  }
450
451
	}
	
Thomas Witkowski's avatar
Thomas Witkowski committed
452
	neighbour[0] = nb;
453
454
455
456
      } else {   /* ichild == 1 */
	// Calculation of the neighbour 2, its oppCoords and the
	// cooresponding oppVertex.

Thomas Witkowski's avatar
Thomas Witkowski committed
457
	neighbour[2] = elInfoOld->neighbour[0];
458
	oppVertex[2] = elInfoOld->oppVertex[0];
459

Thomas Witkowski's avatar
Thomas Witkowski committed
460
461
	if (neighbour[2] && fill_opp_coords) {
	  oppCoord[2] = elInfoOld->oppCoord[0];
462
	  neighbourCoord[2] = elInfoOld->neighbourCoord[0];
463
464
465
466
467
468
469
	}
	

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

	if (elem->getFirstChild()->getFirstChild()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
470
	  neighbour[0] = elem->getFirstChild()->getFirstChild();
471
	  oppVertex[0] = 2;
472
473

	  if (fill_opp_coords) {
474
            if (elem->getFirstChild()->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
475
	      oppCoord[0] = *(elem->getFirstChild()->getNewCoord());
476
	    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
477
478
	      oppCoord[0].setMidpoint(elInfoOld->coord[0], 
				       elInfoOld->coord[2]);
479
	    }
480

Thomas Witkowski's avatar
Thomas Witkowski committed
481
482
483
	    neighbourCoord[0][0] = coord[2];
	    neighbourCoord[0][1] = coord[1];
	    neighbourCoord[0][2] = oppCoord[0];
484
	  }
485
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
486
	  neighbour[0] = elem->getFirstChild();
487
	  oppVertex[0] = 1;
488
489

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

Thomas Witkowski's avatar
Thomas Witkowski committed
492
493
494
	    neighbourCoord[0][0] = elInfoOld->coord[2];
	    neighbourCoord[0][1] = elInfoOld->coord[0];
	    neighbourCoord[0][2] = coord[2];
495
	  }
496
497
498
499
500
	}

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

501
	Element *nb = elInfoOld->neighbour[2];
502
	if (nb) {
503
	  TEST(elInfoOld->oppVertex[2] == 2)("invalid neighbour\n"); 
504
505
506
	  TEST((nb = nb->getFirstChild()))("missing child?\n");

	  if (nb->getFirstChild()) {
507
	    oppVertex[1] = 2;
508

509
	    if (fill_opp_coords) {
510
	      if (nb->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
511
		oppCoord[1] = *(nb->getNewCoord());
512
	      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
513
		oppCoord[1].setMidpoint(elInfoOld->neighbourCoord[2][0],
514
					 elInfoOld->neighbourCoord[2][2]);
515
	      }
516

517
518
519
	      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
520
	      neighbourCoord[1][2] = oppCoord[1];
521
	    }
522
523
	    nb = nb->getSecondChild();

524
	  } else {
525
	    oppVertex[1] = 0;
526

527
	    if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
528
	      oppCoord[1] = elInfoOld->oppCoord[2];
529

530
531
532
533
	      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]);
534
535
536
	    }
	  }
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
537
	neighbour[1] = nb;
538
      } // if (ichild == 0) {} else
Thomas Witkowski's avatar
Thomas Witkowski committed
539
    } // if (fill_flag.isSet(Mesh::FILL_NEIGH) || fillFlag.isSet(Mesh::FILL_OPP_COORDS))
540
541
    

542
    if (fill_flag.isSet(Mesh::FILL_BOUND)) {
543
      if (elInfoOld->getBoundary(2))
Thomas Witkowski's avatar
Thomas Witkowski committed
544
	boundary[5] = elInfoOld->getBoundary(2);
545
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
546
	boundary[5] = INTERIOR;
547

548
      if (ichild == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
549
550
551
552
553
	boundary[3] = elInfoOld->getBoundary(5);
	boundary[4] = elInfoOld->getBoundary(3);
	boundary[0] = elInfoOld->getBoundary(2);
	boundary[1] = INTERIOR;
	boundary[2] = elInfoOld->getBoundary(1);
554
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
555
556
557
558
559
	boundary[3] = elInfoOld->getBoundary(4);
	boundary[4] = elInfoOld->getBoundary(5);
	boundary[0] = INTERIOR;
	boundary[1] = elInfoOld->getBoundary(2);
	boundary[2] = elInfoOld->getBoundary(0);
560
561
      }

562
563
      if (elInfoOld->getProjection(0) && 
	  elInfoOld->getProjection(0)->getType() == VOLUME_PROJECTION) {
564
	
Thomas Witkowski's avatar
Thomas Witkowski committed
565
	projection[0] = elInfoOld->getProjection(0);
566
567
      } else { // boundary projection
	if (ichild == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
568
569
570
	  projection[0] = elInfoOld->getProjection(2);
	  projection[1] = NULL;
	  projection[2] = elInfoOld->getProjection(1);
571
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
572
573
574
	  projection[0] = NULL;
	  projection[1] = elInfoOld->getProjection(2);
	  projection[2] = elInfoOld->getProjection(0);
575
576
577
578
579
	}
      }
    }
  }

580
  double ElInfo2d::calcGrdLambda(DimVec<WorldVector<double> >& grd)
581
  {
582
    FUNCNAME("ElInfo2d::calcGrdLambda()");
583
584

    testFlag(Mesh::FILL_COORDS);
585
586
  
    double adet = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
587
    int dim = mesh->getDim();
588

Thomas Witkowski's avatar
Thomas Witkowski committed
589
    for (int i = 0; i < dimOfWorld; i++) {
590
591
      e1[i] = coord[1][i] - coord[0][i];
      e2[i] = coord[2][i] - coord[0][i];
592
593
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
594
    if (dimOfWorld == 2) {
595
      double sdet = e1[0] * e2[1] - e1[1] * e2[0];
596
597
598
599
      adet = abs(sdet);

      if (adet < 1.0E-25) {
	MSG("abs(det) = %f\n", adet);
600
	for (int i = 0; i <= dim; i++)
601
	  grd[i].set(0.0);
602
603
      } else {
	double det1 = 1.0 / sdet;
604

605
606
607
608
609
610
	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];
611
      }
612
    } else {  
613
      vectorProduct(e1, e2, normal);
614

615
      adet = norm(normal);
616
617
618

      if (adet < 1.0E-15) {
	MSG("abs(det) = %lf\n", adet);
619
	for (int i = 0; i <= dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
620
	  for (int j = 0; j < dimOfWorld; j++)
621
	    grd[i][j] = 0.0;
622
      } else {
623
624
	vectorProduct(e2, normal, grd[1]);
	vectorProduct(normal, e1, grd[2]);
625
      
626
	double adet2 = 1.0 / (adet * adet);
627

Thomas Witkowski's avatar
Thomas Witkowski committed
628
	for (int i = 0; i < dimOfWorld; i++) {
629
630
	  grd[1][i] *= adet2;
	  grd[2][i] *= adet2;
631
632
	}

633
634
635
	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];
636
637
638
639
640
641
      }
    }

    return adet;
  }

642
643

  const int ElInfo2d::worldToCoord(const WorldVector<double>& xy, 
644
645
				   DimVec<double>* lambda) const
  {
646
    FUNCNAME("ElInfo::worldToCoord()");
647

648
    TEST_EXIT_DBG(lambda)("lambda must not be NULL\n");
649

Thomas Witkowski's avatar
Thomas Witkowski committed
650
    DimVec<WorldVector<double> > edge(mesh->getDim(), NO_INIT);
651
    WorldVector<double> x; 
Thomas Witkowski's avatar
Thomas Witkowski committed
652
    static DimVec<double> vec(mesh->getDim(), NO_INIT);
653

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

656
    for (int j = 0; j < dimOfWorld; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
657
      double x0 = coord[dim][j];
658
      x[j] = xy[j] - x0;
659
      for (int i = 0; i < dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
660
	edge[i][j] = coord[i][j] - x0;
661
662
    }
  
663
664
665
    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]; 
666
667
668

    if (abs(det) < DBL_TOL) {
      ERROR("det = %le; abort\n", det);
669
670
      for (int i = 0; i <= dim; i++) 
	(*lambda)[i] = 1.0/dim;
671
672
673
674
675
676
677
678
679
      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;
680
    for (int i = 0; i <= dim; i++) {
681
682
683
684
685
686
687
688
689
690
691
692
      if ((*lambda)[i] < -1.E-5) {
	if ((*lambda)[i] < lmin) {
	  k = i;
	  lmin = (*lambda)[i];
	}
      }
    }

    return k;
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
693
  double ElInfo2d::getNormal(int side, WorldVector<double> &normal)
694
  {
695
    FUNCNAME("ElInfo::getNormal()");
696

697
698
    int i0 = (side + 1) % 3;
    int i1 = (side + 2) % 3;
699

Thomas Witkowski's avatar
Thomas Witkowski committed
700
    if (dimOfWorld == 2){
Thomas Witkowski's avatar
Thomas Witkowski committed
701
702
      normal[0] = coord[i1][1] - coord[i0][1];
      normal[1] = coord[i0][0] - coord[i1][0];
703
    } else { // dow == 3
704
      WorldVector<double> e0, elementNormal;
705

Thomas Witkowski's avatar
Thomas Witkowski committed
706
707
708
709
710
711
      e0 = coord[i1]; 
      e0 -= coord[i0];
      e1 = coord[i1]; 
      e1 -= coord[side];
      e2 = coord[i0]; 
      e2 -= coord[side];
712
713
714
715
716

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

717
    double det = norm(&normal);
718

719
    TEST_EXIT_DBG(det > 1.e-30)("det = 0 on face %d\n", side);
720
721
722

    normal *= 1.0 / det;
    
Thomas Witkowski's avatar
Thomas Witkowski committed
723
    return det;
724
725
  }

726

727
728
729
730
731
  /****************************************************************************/
  /*  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
732
  double ElInfo2d::getElementNormal(WorldVector<double> &elementNormal) const
733
  {
734
    FUNCNAME("ElInfo::getElementNormal()");
735

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

Thomas Witkowski's avatar
Thomas Witkowski committed
739
740
    WorldVector<double> e0 = coord[1] - coord[0];
    WorldVector<double> e1 = coord[2] - coord[0];
741
742
743

    vectorProduct(e0, e1, elementNormal);

744
    double det = norm(&elementNormal);
745

746
    TEST_EXIT_DBG(det > 1.e-30)("det = 0");
747
748
749

    elementNormal *= 1.0 / det;
    
Thomas Witkowski's avatar
Thomas Witkowski committed
750
    return det;
751
  }
752

753

754
  mtl::dense2D<double>& ElInfo2d::getSubElemCoordsMat(int degree) const
755
756
757
758
759
  {
    FUNCNAME("ElInfo2d::getSubElemCoordsMat()");

    using namespace mtl;

760
761
762
763
    if (subElemMatrices[degree].count(refinementPath) == 0) {
      switch (degree) {
      case 1:
	{
764
765
766
	  dense2D<double> mat(3, 3), tmpMat(3, 3);
	  mat = 1;

767
768
769
770
771
772
773
774
	  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;
	    }
775
776
	  }

777
778
779
	  subElemMatrices[1][refinementPath] = mat;  
	}
	break;
780
781
      case 2:
	{
782
783
784
	  dense2D<double> mat(6, 6), tmpMat(6, 6);
	  mat = 1;

785
786
787
788
789
790
791
792
793
794
795
796
797
	  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;
798
799
800
801
802
803
804
805
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
      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;

835
836
      default:
	ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
837
838
      }
    }
839
840

    return subElemMatrices[degree][refinementPath];
841
842
843
  }


844
  mtl::dense2D<double>& ElInfo2d::getSubElemGradCoordsMat(int degree) const
845
  {
846
    FUNCNAME("ElInfo2d::getSubElemGradCoordsMat()");
847

848
849
850
851
    TEST_EXIT(degree == 1)("Not supported for basis functions with degree > 1!\n");

    using namespace mtl;

852
    if (subElemGradMatrices[degree].count(refinementPath) == 0) {
853
854
      dense2D<double> mat(3, 3), tmpMat(3, 3);
      mat = 1;
855
856
857
858
859
860
861
862
863
864
865

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

866
867

      for (int i = 0; i < refinementPathLength; i++)
868
869
870
871
872
873
874
	if (refinementPath & (1 << i)) {
	  tmpMat = mat_right * mat;
	  mat = tmpMat;
	} else  {
	  tmpMat = mat_left * mat;
	  mat = tmpMat;
	}
875
876
877
878
879

      subElemGradMatrices[1][refinementPath] = mat;
    }

    return subElemGradMatrices[degree][refinementPath];
880

881
  }
882
}