ElInfo2d.cc 28 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
#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 {

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

34
  
35
  double ElInfo2d::mat_d1_right_val[3][3] = {{0.0, 0.0, 0.5}, 
36
  					     {1.0, 0.0, 0.5},
37
  					     {0.0, 1.0, 0.0}}; 
38
39
  mtl::dense2D<double> ElInfo2d::mat_d1_right(mat_d1_right_val);

40

41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58

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


59

60
61
  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},
62
					      {1.0, 0.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.0, 0.0, 0.0, -0.25, 0.0, 0.0, -0.125},
64
65
					      {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},
66
67
68
					      {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},
69
70
71
					      {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);

72
73
  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},
74
					       {0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
75
					       {0.0, 0.0, 0.0, -0.25, 0.0, 0.0, 0.0, 1.0, 0.0, 0.375},
76
77
					       {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},
78
79
80
					       {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},
81
82
83
84
85
					       {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);



86
87
  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},
88
					      {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},
89
90
91
92
93
94
95
96
97
98
99
100
					      {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}};
101
102
  mtl::dense2D<double> ElInfo2d::mat_d4_left(mat_d4_left_val);

103
104
  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},
105
					       {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},
106
107
108
109
110
111
112
113
114
115
116
117
					       {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}};
118
119
120
  mtl::dense2D<double> ElInfo2d::mat_d4_right(mat_d4_right_val);


121
122
  ElInfo2d::ElInfo2d(Mesh *aMesh) 
    : ElInfo(aMesh) 
123
  {}
124

125

126
  ElInfo2d::~ElInfo2d()
127
  {}
128

129

130
131
  void ElInfo2d::fillMacroInfo(const MacroElement * mel)
  {
132
133
    FUNCNAME("ElInfo::fillMacroInfo()");
 
Thomas Witkowski's avatar
Thomas Witkowski committed
134
135
136
    macroElement = const_cast<MacroElement*>(mel);
    element = const_cast<Element*>(mel->getElement());
    parent = NULL;
137
    level = 0;
138

Thomas Witkowski's avatar
Thomas Witkowski committed
139
140
141
    if (fillFlag.isSet(Mesh::FILL_COORDS) || 
	fillFlag.isSet(Mesh::FILL_DET)    ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
142

Thomas Witkowski's avatar
Thomas Witkowski committed
143
      int vertices = mesh->getGeo(VERTEX);
144
      for (int i = 0; i < vertices; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
145
	coord[i] = mel->coord[i];
146
147
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
150
151
    if (fillFlag.isSet(Mesh::FILL_OPP_COORDS) || 
	fillFlag.isSet(Mesh::FILL_NEIGH)) {
152

Thomas Witkowski's avatar
Thomas Witkowski committed
153
      bool fill_opp_coords = (fillFlag.isSet(Mesh::FILL_OPP_COORDS));
154
    
155
156
157
158
      for (int i = 0; i < neighbours; i++) {
	MacroElement *macroNeighbour = mel->getNeighbour(i);

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

162
	  int edgeNo = oppVertex[i] = mel->getOppVertex(i);
163
164


165
	  if (nb->getFirstChild() && edgeNo != 2) {  
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185

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


186
	    if (edgeNo == 0) {
187
188
189
190
191
192
193
194
195
196
197
198
	      // The situation is as follows:
	      //
	      //          -------
	      //          \    /|\
	      //           \  / | \
	      //            \/  |  \
	      //             \  |   \
	      //              \ |    \ 
	      //               \|     \
	      //                -------
	      //
	      //            nb     el
199
              //
200
201
202
203
	      // 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
204
	      nb = neighbour[i] = nb->getSecondChild();
205
	    } else {
206
207
	      // 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
208
	      nb = neighbour[i] = nb->getFirstChild();
209
210
	    }

211
212
	    // In both cases the opp vertex number is 2, as one can see in the 
	    // pictures above.
213
	    oppVertex[i] = 2;
214
215

	    if (fill_opp_coords) {
216
	      if (nb->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
217
		oppCoord[i] = *(nb->getNewCoord());
218
	      } else {
219
220
221
		// 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
222
		oppCoord[i] = (macroNeighbour->coord[0] + 
223
224
				macroNeighbour->coord[1]) * 0.5;
	      }
225
226
227
	      
	      switch (i) {
	      case 0:
228
229
230
231
		// The common edge is the edge 0 of this element.

		switch (edgeNo) {
		case 1:
232
233
		  neighbourCoord[i][0] = macroNeighbour->coord[2];
		  neighbourCoord[i][1] = macroNeighbour->coord[0];
234
235
		  break;
		case 0:		  
236
237
		  neighbourCoord[i][0] = macroNeighbour->coord[1];
		  neighbourCoord[i][1] = macroNeighbour->coord[2];
238
239
		  break;
		default:
Thomas Witkowski's avatar
Thomas Witkowski committed
240
		  ERROR_EXIT("Should not happen!\n");
241
		}
242
	
Thomas Witkowski's avatar
Thomas Witkowski committed
243
		neighbourCoord[i][2] = oppCoord[i];
244
245
246
		break;
		
	      case 1:
247
248
249
		// The commonedge is the edge 1 of this element.
		switch (edgeNo) {
		case 0:
250
251
		  neighbourCoord[i][0] = macroNeighbour->coord[1];
		  neighbourCoord[i][1] = macroNeighbour->coord[2];
252
253
		  break;
		case 1:
254
255
		  neighbourCoord[i][0] = macroNeighbour->coord[2];
		  neighbourCoord[i][1] = macroNeighbour->coord[0];
256
257
		  break;
		default:
Thomas Witkowski's avatar
Thomas Witkowski committed
258
		  ERROR_EXIT("Should not happen!\n");
259
260
		}
		
Thomas Witkowski's avatar
Thomas Witkowski committed
261
		neighbourCoord[i][2] = oppCoord[i];
262
263
		break;
		
Thomas Witkowski's avatar
Thomas Witkowski committed
264
	      case 2:
265
		if (*(macroNeighbour->getElement()->getDof(2)) == *(element->getDof(0))) {
266
267
		  neighbourCoord[i][0] = macroNeighbour->coord[2];
		  neighbourCoord[i][1] = macroNeighbour->coord[1];
268
		} else if (*(macroNeighbour->getElement()->getDof(2)) == *(element->getDof(1))) {
269
270
		  neighbourCoord[i][0] = macroNeighbour->coord[0];
		  neighbourCoord[i][1] = macroNeighbour->coord[2];		 
271
272
273
274
		} else {
		  ERROR_EXIT("Should not happen!\n");
		}

275
276
277
		// 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.
278
279
		//		ERROR_EXIT("Should not happen!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
280
281
		break;

282
	      default:
Thomas Witkowski's avatar
Thomas Witkowski committed
283
		std::cout << "------------- Error --------------" << std::endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
284
		std::cout << "  Neighbour counter = " << i << "\n";
Thomas Witkowski's avatar
Thomas Witkowski committed
285
		std::cout << "  Element index     = " << element->getIndex() << "\n\n";
Thomas Witkowski's avatar
Thomas Witkowski committed
286
287
288
289
290
291
292
293
		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;
		  }
294
295
296
		  std::cout << "  OppVertex " << j << ": " 
			    << static_cast<int>(mel->getOppVertex(j)) 
			    << std::endl << std::endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
297
		}
298
299
		ERROR_EXIT("should not happen!\n");
		break;
300
301
302
	      }
	    }
	  } else {
303
304
305
306
307
308
309

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

310
	    if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
311
	      oppCoord[i] = macroNeighbour->coord[edgeNo];
312
	      neighbourCoord[i] = macroNeighbour->coord;	      
313
	    }
314
	  }
315
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
316
	  neighbour[i] = NULL;
317
        }
318
      }
319
    }
320
    
Thomas Witkowski's avatar
Thomas Witkowski committed
321
322
323
324
325
    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);      
326
327
328
329
330
331
332
333
    }
  }


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

334
  void ElInfo2d::fillElInfo(int ichild, const ElInfo *elInfoOld)
335
  {
336
    FUNCNAME("ElInfo::fillElInfo()");
337

Thomas Witkowski's avatar
Thomas Witkowski committed
338
339
    Element *elem = elInfoOld->element;
    Flag fill_flag = elInfoOld->fillFlag;
340

341
    TEST_EXIT_DBG(elem->getFirstChild())("no children?\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
342
    element = const_cast<Element*>((ichild == 0) ? 
343
344
				    elem->getFirstChild() : 
				    elem->getSecondChild());
Thomas Witkowski's avatar
Thomas Witkowski committed
345
    TEST_EXIT_DBG(element)("missing child %d?\n", ichild);
346

Thomas Witkowski's avatar
Thomas Witkowski committed
347
348
349
    macroElement  = elInfoOld->macroElement;
    fillFlag = fill_flag;
    parent = elem;
350
    level = elInfoOld->level + 1;
351
    iChild = ichild;
352

Thomas Witkowski's avatar
Thomas Witkowski committed
353
354
355
    if (fillFlag.isSet(Mesh::FILL_COORDS) || 
	fillFlag.isSet(Mesh::FILL_DET)    ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
356
      
357
      if (elem->isNewCoordSet())
Thomas Witkowski's avatar
Thomas Witkowski committed
358
	coord[2] = *(elem->getNewCoord());
359
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
360
	coord[2].setMidpoint(elInfoOld->coord[0], elInfoOld->coord[1]);      
361
362
      
      if (ichild == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
363
364
	coord[0] = elInfoOld->coord[2];
	coord[1] = elInfoOld->coord[0];
365
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
366
367
	coord[0] = elInfoOld->coord[1];
	coord[1] = elInfoOld->coord[2];
368
369
370
371
372
373
374
375
376
377
      }
    }

    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
378
	neighbour[2] = elInfoOld->neighbour[1];
379
	oppVertex[2] = elInfoOld->oppVertex[1];
380
	
Thomas Witkowski's avatar
Thomas Witkowski committed
381
382
	if (neighbour[2] && fill_opp_coords) {
	  oppCoord[2] = elInfoOld->oppCoord[1];
383
	  neighbourCoord[2] = elInfoOld->neighbourCoord[1];
384
	}
385
386
387
388
389
390
391
392
393
	
	
	// 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
394
	  neighbour[1] = elem->getSecondChild()->getSecondChild();
395
	  oppVertex[1] = 2;
396
397
	  
	  if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
398
399
400
401
402
403
404
405
            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];  
406
407
	  }
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
408
	  neighbour[1] = elem->getSecondChild();
409
	  oppVertex[1] = 0;
410
411

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

Thomas Witkowski's avatar
Thomas Witkowski committed
414
415
416
	    neighbourCoord[1][0] = elInfoOld->coord[1];
	    neighbourCoord[1][1] = elInfoOld->coord[2];
	    neighbourCoord[1][2] = coord[2];
417
418
419
420
	  }
	}


421
422
423
	// Calculation of the neighbour 0, its oppCoords and the
	// cooresponding oppVertex.
	
424
	Element *nb = elInfoOld->neighbour[2];
425
	if (nb) {
426
427
428
429
430
431
432
433
434
435
436
	  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());
437
438
439
440
	 
	  nb = nb->getSecondChild();

	  if (nb->getFirstChild()) {
441
	    oppVertex[0] = 2;
442

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

451
452
453
	      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
454
	      neighbourCoord[0][2] = oppCoord[0];
455
456
457
458
	    }	   
 
	    nb = nb->getFirstChild();
	  } else {
459
	    oppVertex[0] = 1;
460

461
	    if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
462
	      oppCoord[0] = elInfoOld->oppCoord[2];    
463

464
465
466
467
	      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]);
468
469
	    }
	  }
470
471
	}
	
Thomas Witkowski's avatar
Thomas Witkowski committed
472
	neighbour[0] = nb;
473
474
475
476
      } else {   /* ichild == 1 */
	// Calculation of the neighbour 2, its oppCoords and the
	// cooresponding oppVertex.

Thomas Witkowski's avatar
Thomas Witkowski committed
477
	neighbour[2] = elInfoOld->neighbour[0];
478
	oppVertex[2] = elInfoOld->oppVertex[0];
479

Thomas Witkowski's avatar
Thomas Witkowski committed
480
481
	if (neighbour[2] && fill_opp_coords) {
	  oppCoord[2] = elInfoOld->oppCoord[0];
482
	  neighbourCoord[2] = elInfoOld->neighbourCoord[0];
483
484
485
486
487
488
489
	}
	

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

	if (elem->getFirstChild()->getFirstChild()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
490
	  neighbour[0] = elem->getFirstChild()->getFirstChild();
491
	  oppVertex[0] = 2;
492
493

	  if (fill_opp_coords) {
494
            if (elem->getFirstChild()->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
495
	      oppCoord[0] = *(elem->getFirstChild()->getNewCoord());
496
	    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
497
498
	      oppCoord[0].setMidpoint(elInfoOld->coord[0], 
				       elInfoOld->coord[2]);
499
	    }
500

Thomas Witkowski's avatar
Thomas Witkowski committed
501
502
503
	    neighbourCoord[0][0] = coord[2];
	    neighbourCoord[0][1] = coord[1];
	    neighbourCoord[0][2] = oppCoord[0];
504
	  }
505
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
506
	  neighbour[0] = elem->getFirstChild();
507
	  oppVertex[0] = 1;
508
509

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

Thomas Witkowski's avatar
Thomas Witkowski committed
512
513
514
	    neighbourCoord[0][0] = elInfoOld->coord[2];
	    neighbourCoord[0][1] = elInfoOld->coord[0];
	    neighbourCoord[0][2] = coord[2];
515
	  }
516
517
518
519
520
	}

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

521
	Element *nb = elInfoOld->neighbour[2];
522
	if (nb) {
523
	  TEST(elInfoOld->oppVertex[2] == 2)("invalid neighbour\n"); 
524
525
526
	  TEST((nb = nb->getFirstChild()))("missing child?\n");

	  if (nb->getFirstChild()) {
527
	    oppVertex[1] = 2;
528

529
	    if (fill_opp_coords) {
530
	      if (nb->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
531
		oppCoord[1] = *(nb->getNewCoord());
532
	      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
533
		oppCoord[1].setMidpoint(elInfoOld->neighbourCoord[2][0],
534
					 elInfoOld->neighbourCoord[2][2]);
535
	      }
536

537
538
539
	      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
540
	      neighbourCoord[1][2] = oppCoord[1];
541
	    }
542
543
	    nb = nb->getSecondChild();

544
	  } else {
545
	    oppVertex[1] = 0;
546

547
	    if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
548
	      oppCoord[1] = elInfoOld->oppCoord[2];
549

550
551
552
553
	      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]);
554
555
556
	    }
	  }
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
557
	neighbour[1] = nb;
558
      } // if (ichild == 0) {} else
Thomas Witkowski's avatar
Thomas Witkowski committed
559
    } // if (fill_flag.isSet(Mesh::FILL_NEIGH) || fillFlag.isSet(Mesh::FILL_OPP_COORDS))
560
561
    

562
    if (fill_flag.isSet(Mesh::FILL_BOUND)) {
563
      if (elInfoOld->getBoundary(2))
Thomas Witkowski's avatar
Thomas Witkowski committed
564
	boundary[5] = elInfoOld->getBoundary(2);
565
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
566
	boundary[5] = INTERIOR;
567

568
      if (ichild == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
569
570
571
572
573
	boundary[3] = elInfoOld->getBoundary(5);
	boundary[4] = elInfoOld->getBoundary(3);
	boundary[0] = elInfoOld->getBoundary(2);
	boundary[1] = INTERIOR;
	boundary[2] = elInfoOld->getBoundary(1);
574
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
575
576
577
578
579
	boundary[3] = elInfoOld->getBoundary(4);
	boundary[4] = elInfoOld->getBoundary(5);
	boundary[0] = INTERIOR;
	boundary[1] = elInfoOld->getBoundary(2);
	boundary[2] = elInfoOld->getBoundary(0);
580
581
      }

582
583
      if (elInfoOld->getProjection(0) && 
	  elInfoOld->getProjection(0)->getType() == VOLUME_PROJECTION) {
584
	
Thomas Witkowski's avatar
Thomas Witkowski committed
585
	projection[0] = elInfoOld->getProjection(0);
586
587
      } else { // boundary projection
	if (ichild == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
588
589
590
	  projection[0] = elInfoOld->getProjection(2);
	  projection[1] = NULL;
	  projection[2] = elInfoOld->getProjection(1);
591
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
592
593
594
	  projection[0] = NULL;
	  projection[1] = elInfoOld->getProjection(2);
	  projection[2] = elInfoOld->getProjection(0);
595
596
597
598
599
	}
      }
    }
  }

600
  double ElInfo2d::calcGrdLambda(DimVec<WorldVector<double> >& grd)
601
  {
602
    FUNCNAME("ElInfo2d::calcGrdLambda()");
603
604

    testFlag(Mesh::FILL_COORDS);
605
606
  
    double adet = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
607
    int dim = mesh->getDim();
608

Thomas Witkowski's avatar
Thomas Witkowski committed
609
    for (int i = 0; i < dimOfWorld; i++) {
610
611
      e1[i] = coord[1][i] - coord[0][i];
      e2[i] = coord[2][i] - coord[0][i];
612
613
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
614
    if (dimOfWorld == 2) {
615
      double sdet = e1[0] * e2[1] - e1[1] * e2[0];
616
617
618
619
      adet = abs(sdet);

      if (adet < 1.0E-25) {
	MSG("abs(det) = %f\n", adet);
620
	for (int i = 0; i <= dim; i++)
621
	  grd[i].set(0.0);
622
623
      } else {
	double det1 = 1.0 / sdet;
624

625
626
627
628
629
630
	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];
631
      }
632
    } else {  
633
      vectorProduct(e1, e2, normal);
634

635
      adet = norm(normal);
636
637
638

      if (adet < 1.0E-15) {
	MSG("abs(det) = %lf\n", adet);
639
	for (int i = 0; i <= dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
640
	  for (int j = 0; j < dimOfWorld; j++)
641
	    grd[i][j] = 0.0;
642
      } else {
643
644
	vectorProduct(e2, normal, grd[1]);
	vectorProduct(normal, e1, grd[2]);
645
      
646
	double adet2 = 1.0 / (adet * adet);
647

Thomas Witkowski's avatar
Thomas Witkowski committed
648
	for (int i = 0; i < dimOfWorld; i++) {
649
650
	  grd[1][i] *= adet2;
	  grd[2][i] *= adet2;
651
652
	}

653
654
655
	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];
656
657
658
659
660
661
      }
    }

    return adet;
  }

662
663

  const int ElInfo2d::worldToCoord(const WorldVector<double>& xy, 
664
665
				   DimVec<double>* lambda) const
  {
666
    FUNCNAME("ElInfo::worldToCoord()");
667

668
    TEST_EXIT_DBG(lambda)("lambda must not be NULL\n");
669

Thomas Witkowski's avatar
Thomas Witkowski committed
670
    DimVec<WorldVector<double> > edge(mesh->getDim(), NO_INIT);
671
    WorldVector<double> x; 
Thomas Witkowski's avatar
Thomas Witkowski committed
672
    static DimVec<double> vec(mesh->getDim(), NO_INIT);
673

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

676
    for (int j = 0; j < dimOfWorld; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
677
      double x0 = coord[dim][j];
678
      x[j] = xy[j] - x0;
679
      for (int i = 0; i < dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
680
	edge[i][j] = coord[i][j] - x0;
681
    }
682

683
684
685
    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]; 
686
687
688

    if (abs(det) < DBL_TOL) {
      ERROR("det = %le; abort\n", det);
689
      for (int i = 0; i <= dim; i++) 
690
	(*lambda)[i] = 1.0 / dim;
691
692
693
694
695
696
697
698
699
      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;
700
    for (int i = 0; i <= dim; i++) {
701
      if ((*lambda)[i] < -1.e-5) {
702
703
704
705
706
707
708
709
710
711
712
	if ((*lambda)[i] < lmin) {
	  k = i;
	  lmin = (*lambda)[i];
	}
      }
    }

    return k;
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
713
  double ElInfo2d::getNormal(int side, WorldVector<double> &normal)
714
  {
715
    FUNCNAME("ElInfo::getNormal()");
716

717
718
    int i0 = (side + 1) % 3;
    int i1 = (side + 2) % 3;
719

Thomas Witkowski's avatar
Thomas Witkowski committed
720
    if (dimOfWorld == 2){
Thomas Witkowski's avatar
Thomas Witkowski committed
721
722
      normal[0] = coord[i1][1] - coord[i0][1];
      normal[1] = coord[i0][0] - coord[i1][0];
723
    } else { // dow == 3
724
      WorldVector<double> e0, elementNormal;
725

Thomas Witkowski's avatar
Thomas Witkowski committed
726
727
728
729
730
731
      e0 = coord[i1]; 
      e0 -= coord[i0];
      e1 = coord[i1]; 
      e1 -= coord[side];
      e2 = coord[i0]; 
      e2 -= coord[side];
732
733
734
735
736

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

737
    double det = norm(&normal);
738

739
    TEST_EXIT_DBG(det > 1.e-30)("det = 0 on face %d\n", side);
740
741
742

    normal *= 1.0 / det;
    
Thomas Witkowski's avatar
Thomas Witkowski committed
743
    return det;
744
745
  }

746

747
748
749
750
751
  /****************************************************************************/
  /*  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
752
  double ElInfo2d::getElementNormal(WorldVector<double> &elementNormal) const
753
  {
754
    FUNCNAME("ElInfo::getElementNormal()");
755

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

Thomas Witkowski's avatar
Thomas Witkowski committed
759
760
    WorldVector<double> e0 = coord[1] - coord[0];
    WorldVector<double> e1 = coord[2] - coord[0];
761
762
763

    vectorProduct(e0, e1, elementNormal);

764
    double det = norm(&elementNormal);
765

766
    TEST_EXIT_DBG(det > 1.e-30)("det = 0");
767
768
769

    elementNormal *= 1.0 / det;
    
Thomas Witkowski's avatar
Thomas Witkowski committed
770
    return det;
771
  }
772

773

774
  mtl::dense2D<double>& ElInfo2d::getSubElemCoordsMat(int degree) const
775
776
777
778
779
  {
    FUNCNAME("ElInfo2d::getSubElemCoordsMat()");

    using namespace mtl;

Thomas Witkowski's avatar
Thomas Witkowski committed
780
    if (subElemMatrices[degree].count(std::make_pair(refinementPathLength, refinementPath)) == 0) {
781
782
783
      switch (degree) {
      case 1:
	{
784
785
786
	  dense2D<double> mat(3, 3), tmpMat(3, 3);
	  mat = 1;

787
788
	  for (int i = 0; i < refinementPathLength; i++) {
	    if (refinementPath & (1 << i)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
789
	      tmpMat = mat * mat_d1_right;
790
791
	      mat = tmpMat;
	    } else  {
Thomas Witkowski's avatar
Thomas Witkowski committed
792
	      tmpMat = mat * mat_d1_left;
793
794
	      mat = tmpMat;
	    }
795
796
	  }

Thomas Witkowski's avatar
Thomas Witkowski committed
797
	  subElemMatrices[1][std::make_pair(refinementPathLength, refinementPath)] = mat;  
798
799
	}
	break;
800
801
      case 2:
	{
802
803
804
	  dense2D<double> mat(6, 6), tmpMat(6, 6);
	  mat = 1;

805
806
	  for (int i = 0; i < refinementPathLength; i++) {
	    if (refinementPath & (1 << i)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
807
	      tmpMat = mat * mat_d2_right;
808
809
	      mat = tmpMat;
	    } else  {
Thomas Witkowski's avatar
Thomas Witkowski committed
810
	      tmpMat = mat * mat_d2_left;
811
812
813
814
	      mat = tmpMat;
	    }
	  }

Thomas Witkowski's avatar
Thomas Witkowski committed
815
	  subElemMatrices[2][std::make_pair(refinementPathLength, refinementPath)] = mat;
816
817
	}
	break;
818
819
820
821
822
823
824
      case 3:
	{
	  dense2D<double> mat(10, 10), tmpMat(10, 10);
	  mat = 1;

	  for (int i = 0; i < refinementPathLength; i++) {
	    if (refinementPath & (1 << i)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
825
	      tmpMat = mat * mat_d3_right;
826
827
	      mat = tmpMat;
	    } else  {
Thomas Witkowski's avatar
Thomas Witkowski committed
828
	      tmpMat = mat * mat_d3_left;
829
830
831
832
	      mat = tmpMat;
	    }
	  }

Thomas Witkowski's avatar
Thomas Witkowski committed
833
	  subElemMatrices[3][std::make_pair(refinementPathLength, refinementPath)] = mat;
834
835
836
837
838
839
840
841
842
	}
	break;
      case 4:
	{
	  dense2D<double> mat(15, 15), tmpMat(15, 15);
	  mat = 1;

	  for (int i = 0; i < refinementPathLength; i++) {
	    if (refinementPath & (1 << i)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
843
	      tmpMat = mat * mat_d4_right;
844
845
	      mat = tmpMat;
	    } else  {
Thomas Witkowski's avatar
Thomas Witkowski committed
846
	      tmpMat = mat * mat_d4_left;
847
848
849
850
	      mat = tmpMat;
	    }
	  }

Thomas Witkowski's avatar
Thomas Witkowski committed
851
	  subElemMatrices[4][std::make_pair(refinementPathLength, refinementPath)] = mat;
852
853
854
	}
	break;

855
856
      default:
	ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
857
858
      }
    }
859

Thomas Witkowski's avatar
Thomas Witkowski committed
860
    return subElemMatrices[degree][std::make_pair(refinementPathLength, refinementPath)];
861
862
863
  }


864
  mtl::dense2D<double>& ElInfo2d::getSubElemGradCoordsMat(int degree) const
865
  {
866
    FUNCNAME("ElInfo2d::getSubElemGradCoordsMat()");
867

868
869
870
871
    TEST_EXIT(degree == 1)("Not supported for basis functions with degree > 1!\n");

    using namespace mtl;

Thomas Witkowski's avatar
Thomas Witkowski committed
872
    if (subElemGradMatrices[degree].count(std::make_pair(refinementPathLength, refinementPath)) == 0) {
873
874
      dense2D<double> mat(3, 3), tmpMat(3, 3);
      mat = 1;
875
876
877
878
879
880
881
882
883
884
885

      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<