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

namespace AMDiS {

17
  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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108

  double ElInfo2d::mat_d3_left_val[10][10] = {{0.0, 1.0, 0.0, 0.3125, 0.0, 0.0, 0.0625, 0.0, 0.0, 0.0},
					      {0.0, 0.0, 0.0, 0.0625, 0.0, 0.0, 0.0625, 0.0, 0.0, 0.0625},
					      {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, 0.0, 0.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.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0},
					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.375},
					      {0.0, 0.0, 0.5625, 0.9375, 1.0, 0.0, 0.0, 0.0, 0.0, 0.1875},
					      {0.0, 0.0, 0.5625, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
					      {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);

  double ElInfo2d::mat_d3_right_val[10][10] = {{0.0, 0.0, 0.0, 0.0625, 0.0, 0.0, 0.0625, 0.0, 0.0, 0.0625},
					       {1.0, 0.0, 0.0, 0.0625, 0.0, 0.0, 0.3125, 0.0, 0.0, 0.0},
					       {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, 0.0, 0.0, 1.0, 0.0, 0.375},
					       {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},
					       {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.5625, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
					       {0.0, 0.0, 0.5625, 0.0, 0.0, 1.0, 0.9375, 0.0, 0.0, 0.1875},
					       {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);



  double ElInfo2d::mat_d4_left_val[15][15] = {{0.0, 1.0, 0.0, 0.273437, 0.0, 0.0, 0.023437, 0.0, 0.0, 0.0, 0.0, 0.0, 0.023437, 0.0, 0.0},
					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.023437, 0.023437, 0.0, 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, 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.0, 0.0, 0.1875, 0.0, 0.0, 0.0, 0.125, 0.0625, 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.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.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, 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, 1.0, 0.0, 0.375, 0.0, 0.0},
					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.1875, 0.0, 0.0, 1.0, 0.0, 0.3125, 0.0},
					      {0.0, 0.0, 0.0, 1.093750, 1.0, 0.468750, 0.0, 0.0, 0.031250, 0.0, 0.0, 0.0, 0.0, 0.156250, 0.0},
					      {0.0, 0.0, 1.0, 0.0, 0.0, 0.703125, 0.140625, 0.0, 0.015625, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
					      {0.0, 0.0, 0.0, 0.218750, 0.0, 0.0, 0.0, 0.0, 0.031250, 0.0, 0.0, 0.0, 0.093750, 0.156250, 0.0},
					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5625, 0.0, 0.0, 0.0, 0.0, 0.0, 0.375, 0.9375, 1.0},
					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5625, 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.0, 1.0, 0.75, 0.0, 0.0, 0.0, 0.75, 0.0, 0.0}};  
  mtl::dense2D<double> ElInfo2d::mat_d4_left(mat_d4_left_val);

  double ElInfo2d::mat_d4_right_val[15][15] = {{0.0, 0.0, 0.0, 0.0, 0.0, 0.023437, 0.023437, 0.0, 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.023437, 0.0, 0.0, 0.273437, 0.0, 0.0, 0.0, 0.0, 0.023437, 0.0},
					       {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},
					       {0.0, 0.0, 0.0, 0.1875, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.3125, 0.0, 0.0},
					       {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.375, 0.0},
					       {0.0, 0.0, 0.0, 0.5, 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, 0.5, 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.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, 0.1875, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0625, 0.125, 0.0},
					       {0.0, 0.0, 0.0, 0.031250, 0.0, 0.0, 0.0, 0.0, 0.218750, 0.0, 0.0, 0.0, 0.156250, 0.093750, 0.0},
					       {0.0, 0.0, 1.0, 0.015625, 0.0, 0.140625, 0.703125, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
					       {0.0, 0.0, 0.0, 0.031250, 0.0, 0.0, 0.468750, 1.0, 1.093750, 0.0, 0.0, 0.0, 0.156250, 0.0, 0.0},
					       {0.0, 0.0, 0.0, 0.0, 0.0, 0.5625, 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.5625, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.9375, 0.375, 1.0},
					       {0.0, 0.0, 0.0, 0.75, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.75, 0.0}};  
  mtl::dense2D<double> ElInfo2d::mat_d4_right(mat_d4_right_val);


109
110
111
  ElInfo2d::ElInfo2d(Mesh *aMesh) 
    : ElInfo(aMesh) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
112
113
114
    e1 = new WorldVector<double>;
    e2 = new WorldVector<double>;
    normal = new WorldVector<double>;
115
116
  }

117

118
119
  ElInfo2d::~ElInfo2d()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
120
121
122
    delete e1;
    delete e2;
    delete normal;
123
124
  }

125

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

Thomas Witkowski's avatar
Thomas Witkowski committed
135
136
137
    if (fillFlag.isSet(Mesh::FILL_COORDS) || 
	fillFlag.isSet(Mesh::FILL_DET)    ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
138

Thomas Witkowski's avatar
Thomas Witkowski committed
139
      int vertices = mesh->getGeo(VERTEX);
140
      for (int i = 0; i < vertices; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
141
	coord[i] = mel->coord[i];
142
143
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
146
147
    if (fillFlag.isSet(Mesh::FILL_OPP_COORDS) || 
	fillFlag.isSet(Mesh::FILL_NEIGH)) {
148

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

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

158
	  int edgeNo = oppVertex[i] = mel->getOppVertex(i);
159
160


161
	  if (nb->getFirstChild() && edgeNo != 2) {  
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181

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


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

207
208
	    // In both cases the opp vertex number is 2, as one can see in the 
	    // pictures above.
209
	    oppVertex[i] = 2;
210
211

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

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

271
272
273
		// 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.
274
275
		//		ERROR_EXIT("Should not happen!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
276
277
		break;

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

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

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


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

330
  void ElInfo2d::fillElInfo(int ichild, const ElInfo *elInfoOld)
331
  {
332
    FUNCNAME("ElInfo::fillElInfo()");
333

Thomas Witkowski's avatar
Thomas Witkowski committed
334
335
    Element *elem = elInfoOld->element;
    Flag fill_flag = elInfoOld->fillFlag;
336

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

Thomas Witkowski's avatar
Thomas Witkowski committed
343
344
345
    macroElement  = elInfoOld->macroElement;
    fillFlag = fill_flag;
    parent = elem;
346
    level = elInfoOld->level + 1;
347
    iChild = ichild;
348

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
410
411
412
	    neighbourCoord[1][0] = elInfoOld->coord[1];
	    neighbourCoord[1][1] = elInfoOld->coord[2];
	    neighbourCoord[1][2] = coord[2];
413
414
415
416
	  }
	}


417
418
419
	// Calculation of the neighbour 0, its oppCoords and the
	// cooresponding oppVertex.
	
420
	Element *nb = elInfoOld->neighbour[2];
421
	if (nb) {
422
	  TEST(elInfoOld->oppVertex[2] == 2)("invalid neighbour\n"); 
423
424
	  TEST_EXIT_DBG(nb->getFirstChild())("missing first child?\n");
	  TEST_EXIT_DBG(nb->getSecondChild())("missing second child?\n");
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
	}
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
588
  double ElInfo2d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
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++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
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++)
Thomas Witkowski's avatar
Thomas Witkowski committed
609
	  for (int j = 0; j < dimOfWorld; j++)
610
	    grd_lam[i][j] = 0.0;
611
612
      } else {
	double det1 = 1.0 / sdet;
613
614
615
616
617

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

624
      adet = norm(normal);
625
626
627

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

Thomas Witkowski's avatar
Thomas Witkowski committed
637
	for (int i = 0; i < dimOfWorld; i++) {
638
639
640
641
642
643
644
645
646
647
648
649
650
	  grd_lam[1][i] *= adet2;
	  grd_lam[2][i] *= adet2;
	}

	grd_lam[0][0] = - grd_lam[1][0] - grd_lam[2][0];
	grd_lam[0][1] = - grd_lam[1][1] - grd_lam[2][1];
	grd_lam[0][2] = - grd_lam[1][2] - grd_lam[2][2];
      }
    }

    return adet;
  }

651
652

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

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

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

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

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

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

    return k;
  }


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

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

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

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

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

726
    double det = norm(&normal);
727

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

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

735

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

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

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

    vectorProduct(e0, e1, elementNormal);

753
    double det = norm(&elementNormal);
754

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

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

762

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

    using namespace mtl;

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

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

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

794
795
796
797
798
799
800
801
802
803
804
805
806
	  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;
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
843
      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;

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

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


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

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

    using namespace mtl;

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

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

875
876

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

      subElemGradMatrices[1][refinementPath] = mat;
    }

    return subElemGradMatrices[degree][refinementPath];
889

890
  }
891
}