ElInfo2d.cc 27.6 KB
Newer Older
1 2 3 4 5 6 7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9 10 11 12 13 14 15 16 17
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
18
 *
19
 ******************************************************************************/
20 21


22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
#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 {

38
  double ElInfo2d::mat_d1_left_val[3][3] = {{0.0, 1.0, 0.5},
39
  					    {0.0, 0.0, 0.5},
40
  					    {1.0, 0.0, 0.0}};
41 42
  mtl::dense2D<double> ElInfo2d::mat_d1_left(mat_d1_left_val);

43 44

  double ElInfo2d::mat_d1_right_val[3][3] = {{0.0, 0.0, 0.5},
45
  					     {1.0, 0.0, 0.5},
46
  					     {0.0, 1.0, 0.0}};
47 48
  mtl::dense2D<double> ElInfo2d::mat_d1_right(mat_d1_right_val);

49

50

51 52
  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},
53 54 55 56 57 58
					    {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);

59 60
  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},
61 62 63 64 65 66 67
					     {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);


68

69 70
  double ElInfo2d::mat_d3_left_val[10][10] = {{0.0,  1.0, -6.25e-02,  3.125e-01,  0.0,  0.0,  6.25e-02,  0.0,  0.0, -6.25e-02},
					      {0.0,  0.0, -6.25e-02,  6.25e-02,  0.0,  0.0,  6.25e-02,  0.0,  0.0, 6.25e-02},
71
					      {1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
72
					      {0.0,  0.0,  0.0,  0.0,  0.0,  0.0, -2.5e-01,  0.0,  0.0, -0.125},
73 74
					      {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},
75 76 77 78
					      {0.0,  0.0,  0.0,  0.0,  0.0,  0.0, -2.5e-01,  0.0,  1.0, 0.375},
					      {0.0,  0.0,  5.625e-01,  9.375e-01,  1.0,  0.0, -6.25e-02,  0.0,  0.0, 1.875e-01},
					      {0.0,  0.0,  5.625e-01, -3.125e-01,  0.0,  0.0, -6.25e-02,  0.0,  0.0, -1.875e-01},
					      {0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 7.5e-01}};
79 80
  mtl::dense2D<double> ElInfo2d::mat_d3_left(mat_d3_left_val);

81 82 83 84 85 86 87 88 89 90
  double ElInfo2d::mat_d3_right_val[10][10] = {{0.0,  0.0, -6.25e-02,  6.25e-02,  0.0,  0.0,  6.25e-02,  0.0,  0.0, 6.25e-02},
					      {1.0,  0.0, -6.25e-02,  6.25e-02,  0.0,  0.0,  3.125e-01,  0.0,  0.0, -6.25e-02},
					      {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, -2.5e-01,  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, -2.5e-01,  0.0,  0.0,  0.0,  0.0,  0.0, -0.125},
					      {0.0,  0.0,  5.625e-01, -6.25e-02,  0.0,  0.0, -3.125e-01,  0.0,  0.0, -1.875e-01},
					      {0.0,  0.0,  5.625e-01, -6.25e-02,  0.0,  1.0,  9.375e-01,  0.0,  0.0, 1.875e-01},
					      {0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.0, 7.5e-01}};
91 92 93 94
  mtl::dense2D<double> ElInfo2d::mat_d3_right(mat_d3_right_val);



95 96
  double ElInfo2d::mat_d4_left_val[15][15] = {{0.0,  1.0,  0.0,  2.734375e-01,  0.0, -3.906250e-02,  2.343750e-02,  0.0, -3.906250e-02,  0.0,  0.0,  0.0,  2.343750e-02, -3.906250e-02, 0.0},
					      {0.0,  0.0,  0.0, -3.906250e-02,  0.0,  2.343750e-02,  2.343750e-02,  0.0, -3.906250e-02,  0.0,  0.0,  0.0, -3.906250e-02, -3.906250e-02, 0.0},
97
					      {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},
98 99 100 101 102 103 104 105 106 107 108
					      {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,  0.125,  6.25e-02, 0.0},
					      {0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0, -0.375,  0.0,  0.0,  0.0, -0.125,  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.375,  0.0,  1.0,  0.0,  0.375,  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, -0.125,  3.125e-01, 0.0},
					      {0.0,  0.0,  0.0,  1.093750e+00,  1.0,  4.687500e-01, -9.375e-02,  0.0,  3.125e-02,  0.0,  0.0,  0.0, -3.125e-02,  1.562500e-01, 0.0},
					      {0.0,  0.0,  1.0, -5.468750e-01,  0.0,  7.031250e-01,  1.406250e-01,  0.0,  1.562500e-02,  0.0,  0.0,  0.0, -4.687500e-02, -2.343750e-01, 0.0},
					      {0.0,  0.0,  0.0,  2.187500e-01,  0.0, -1.562500e-01, -9.375e-02,  0.0,  3.125e-02,  0.0,  0.0,  0.0,  9.375e-02,  1.562500e-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,  0.375,  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, -0.375, -3.125e-01, 0.0},
109
					      {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}};
110 111
  mtl::dense2D<double> ElInfo2d::mat_d4_left(mat_d4_left_val);

112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
  double ElInfo2d::mat_d4_right_val[15][15] = {{0.0,  0.0,  0.0, -3.906250e-02,  0.0,  2.343750e-02,  2.343750e-02,  0.0, -3.906250e-02,  0.0,  0.0,  0.0, -3.906250e-02, -3.906250e-02, 0.0},
						{1.0,  0.0,  0.0, -3.906250e-02,  0.0,  2.343750e-02, -3.906250e-02,  0.0,  2.734375e-01,  0.0,  0.0,  0.0, -3.906250e-02,  2.343750e-02, 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,  1.875e-01,  0.0, -6.25e-02,  0.0,  0.0,  0.0,  1.0,  0.0,  0.0,  3.125e-01, -0.125, 0.0},
						{0.0,  0.0,  0.0, -0.375,  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.375,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0, -0.125, 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,  0.125, 0.0},
						{0.0,  0.0,  0.0,  3.125e-02,  0.0, -9.375e-02, -1.562500e-01,  0.0,  2.187500e-01,  0.0,  0.0,  0.0,  1.562500e-01,  9.375e-02, 0.0},
						{0.0,  0.0,  1.0,  1.562500e-02,  0.0,  1.406250e-01,  7.031250e-01,  0.0, -5.468750e-01,  0.0,  0.0,  0.0, -2.343750e-01, -4.687500e-02, 0.0},
						{0.0,  0.0,  0.0,  3.125e-02,  0.0, -9.375e-02,  4.687500e-01,  1.0,  1.093750e+00,  0.0,  0.0,  0.0,  1.562500e-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, -0.375, 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,  0.375, 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}};
127 128 129
  mtl::dense2D<double> ElInfo2d::mat_d4_right(mat_d4_right_val);


130 131
  ElInfo2d::ElInfo2d(Mesh *aMesh)
    : ElInfo(aMesh)
132
  {}
133

134

135
  ElInfo2d::~ElInfo2d()
136
  {}
137

138

139 140
  void ElInfo2d::fillMacroInfo(const MacroElement * mel)
  {
141
    FUNCNAME("ElInfo::fillMacroInfo()");
142

Thomas Witkowski's avatar
Thomas Witkowski committed
143 144
    macroElement = const_cast<MacroElement*>(mel);
    element = const_cast<Element*>(mel->getElement());
145
    parent = NULL;
146
    level = 0;
147

148
    if (fillFlag.isSet(Mesh::FILL_COORDS) ||
Thomas Witkowski's avatar
Thomas Witkowski committed
149 150
	fillFlag.isSet(Mesh::FILL_DET)    ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
151

Thomas Witkowski's avatar
Thomas Witkowski committed
152
      int vertices = mesh->getGeo(VERTEX);
153
      for (int i = 0; i < vertices; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
154
	coord[i] = mel->coord[i];
155 156
    }

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

159
    if (fillFlag.isSet(Mesh::FILL_OPP_COORDS) ||
Thomas Witkowski's avatar
Thomas Witkowski committed
160
	fillFlag.isSet(Mesh::FILL_NEIGH)) {
161

Thomas Witkowski's avatar
Thomas Witkowski committed
162
      bool fill_opp_coords = (fillFlag.isSet(Mesh::FILL_OPP_COORDS));
163

164 165 166 167
      for (int i = 0; i < neighbours; i++) {
	MacroElement *macroNeighbour = mel->getNeighbour(i);

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

171
	  int edgeNo = oppVertex[i] = mel->getOppVertex(i);
172 173


174
	  if (nb->getFirstChild() && edgeNo != 2) {
175
            /*
176 177 178
	    * 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:
179 180 181 182 183
	    *
	    *               /|\
	    *              / | \
	    *             /  |  \
	    *            /\  |   \
184
	    *           /  \ |    \
185 186 187 188 189
	    *          /    \|     \
	    *          -------------
	    *
	    *            nb     el
	    *
190 191
	    * Note that we know (because of the last if statement), that the
	    * neighbour element has children and the common edge is not the
192 193
	    * refinement edge, which has always the number 2, of our element.
	    */
194

195
	    if (edgeNo == 0) {
196
	      /*
197 198 199 200 201 202 203
	      * The situation is as follows:
	      *
	      *          -------
	      *          \    /|\
	      *           \  / | \
	      *            \/  |  \
	      *             \  |   \
204
	      *              \ |    \
205 206 207 208 209 210 211 212 213
	      *               \|     \
	      *                -------
	      *
	      *            nb     el
              *
	      * 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
214
	      nb = neighbour[i] = nb->getSecondChild();
215
	    } else {
216 217
	      // 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
218
	      nb = neighbour[i] = nb->getFirstChild();
219 220
	    }

221
	    // In both cases the opp vertex number is 2, as one can see in the
222
	    // pictures above.
223
	    oppVertex[i] = 2;
224 225

	    if (fill_opp_coords) {
226
	      if (nb->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
227
		oppCoord[i] = *(nb->getNewCoord());
228
	      } else {
229 230 231
		// 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.
232
		oppCoord[i] = (macroNeighbour->coord[0] +
233 234
				macroNeighbour->coord[1]) * 0.5;
	      }
235

236 237
	      switch (i) {
	      case 0:
238 239 240 241
		// The common edge is the edge 0 of this element.

		switch (edgeNo) {
		case 1:
242 243
		  neighbourCoord[i][0] = macroNeighbour->coord[2];
		  neighbourCoord[i][1] = macroNeighbour->coord[0];
244
		  break;
245
		case 0:
246 247
		  neighbourCoord[i][0] = macroNeighbour->coord[1];
		  neighbourCoord[i][1] = macroNeighbour->coord[2];
248 249
		  break;
		default:
Thomas Witkowski's avatar
Thomas Witkowski committed
250
		  ERROR_EXIT("Should not happen!\n");
251
		}
252

Thomas Witkowski's avatar
Thomas Witkowski committed
253
		neighbourCoord[i][2] = oppCoord[i];
254
		break;
255

256
	      case 1:
Thomas Witkowski's avatar
a  
Thomas Witkowski committed
257
		// The common edge is the edge 1 of this element.
258 259
		switch (edgeNo) {
		case 0:
260 261
		  neighbourCoord[i][0] = macroNeighbour->coord[1];
		  neighbourCoord[i][1] = macroNeighbour->coord[2];
262 263
		  break;
		case 1:
264 265
		  neighbourCoord[i][0] = macroNeighbour->coord[2];
		  neighbourCoord[i][1] = macroNeighbour->coord[0];
266 267
		  break;
		default:
Thomas Witkowski's avatar
Thomas Witkowski committed
268
		  ERROR_EXIT("Should not happen!\n");
269
		}
270

Thomas Witkowski's avatar
Thomas Witkowski committed
271
		neighbourCoord[i][2] = oppCoord[i];
272
		break;
273

Thomas Witkowski's avatar
Thomas Witkowski committed
274
	      case 2:
275
		if (*(macroNeighbour->getElement()->getDof(2)) == *(element->getDof(0))) {
276 277
		  neighbourCoord[i][0] = macroNeighbour->coord[2];
		  neighbourCoord[i][1] = macroNeighbour->coord[1];
278
		} else if (*(macroNeighbour->getElement()->getDof(2)) == *(element->getDof(1))) {
279
		  neighbourCoord[i][0] = macroNeighbour->coord[0];
280
		  neighbourCoord[i][1] = macroNeighbour->coord[2];
281
		} else {
282
		  ERROR_EXIT("Should not happen! Non-conforming AMDiS-mesh? Periodic mesh with corrected index-circles?\n");
283 284
		}

285 286 287
		// 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.
288 289
		//		ERROR_EXIT("Should not happen!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
290 291
		break;

292
	      default:
Thomas Witkowski's avatar
Thomas Witkowski committed
293
		std::cout << "------------- Error --------------" << std::endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
294
		std::cout << "  Neighbour counter = " << i << "\n";
Thomas Witkowski's avatar
Thomas Witkowski committed
295
		std::cout << "  Element index     = " << element->getIndex() << "\n\n";
Thomas Witkowski's avatar
Thomas Witkowski committed
296 297
		for (int j = 0; j < neighbours; j++) {
		  if (mel->getNeighbour(j)) {
298 299
		    std::cout << "  Neighbour " << j << ": "
			      << mel->getNeighbour(j)->getElement()->getIndex()
Thomas Witkowski's avatar
Thomas Witkowski committed
300 301 302 303
			      << std::endl;
		  } else {
		    std::cout << "  Neighbour " << j << ": not existing" << std::endl;
		  }
304 305
		  std::cout << "  OppVertex " << j << ": "
			    << static_cast<int>(mel->getOppVertex(j))
306
			    << std::endl << std::endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
307
		}
308 309
		ERROR_EXIT("should not happen!\n");
		break;
310 311 312
	      }
	    }
	  } else {
313 314 315 316

	    // 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
317
	    // level neighbour would be refinement, also this element must to be
318 319
	    // refinement, because they share the refinement edge.

320
	    if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
321
	      oppCoord[i] = macroNeighbour->coord[edgeNo];
322
	      neighbourCoord[i] = macroNeighbour->coord;
323
	    }
324
	  }
325
	} else {
326
	  neighbour[i] = NULL;
327
        }
328
      }
329
    }
330 331

    if (fillFlag.isSet(Mesh::FILL_BOUND)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
332 333 334
      for (int i = 0; i < element->getGeo(BOUNDARY); i++)
	boundary[i] = mel->getBoundary(i);
      for (int i = 0; i < element->getGeo(PROJECTION); i++)
335
	projection[i] = mel->getProjection(i);
336 337 338 339 340 341 342 343
    }
  }


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

344
  void ElInfo2d::fillElInfo(int ichild, const ElInfo *elInfoOld)
345
  {
346
    FUNCNAME("ElInfo::fillElInfo()");
347

Thomas Witkowski's avatar
Thomas Witkowski committed
348 349
    Element *elem = elInfoOld->element;
    Flag fill_flag = elInfoOld->fillFlag;
350

351
    TEST_EXIT_DBG(elem->getFirstChild())("no children?\n");
352 353
    element = const_cast<Element*>((ichild == 0) ?
				    elem->getFirstChild() :
354
				    elem->getSecondChild());
Thomas Witkowski's avatar
Thomas Witkowski committed
355
    TEST_EXIT_DBG(element)("missing child %d?\n", ichild);
356

Thomas Witkowski's avatar
Thomas Witkowski committed
357 358 359
    macroElement  = elInfoOld->macroElement;
    fillFlag = fill_flag;
    parent = elem;
360
    level = elInfoOld->level + 1;
361
    iChild = ichild;
362

363
    if (fillFlag.isSet(Mesh::FILL_COORDS) ||
Thomas Witkowski's avatar
Thomas Witkowski committed
364 365
	fillFlag.isSet(Mesh::FILL_DET)    ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
366

367
      if (elem->isNewCoordSet())
Thomas Witkowski's avatar
Thomas Witkowski committed
368
	coord[2] = *(elem->getNewCoord());
369
      else
370 371
	coord[2].setMidpoint(elInfoOld->coord[0], elInfoOld->coord[1]);

372
      if (ichild == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
373 374
	coord[0] = elInfoOld->coord[2];
	coord[1] = elInfoOld->coord[0];
375
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
376 377
	coord[0] = elInfoOld->coord[1];
	coord[1] = elInfoOld->coord[2];
378 379 380 381 382
      }
    }

    bool fill_opp_coords = (fill_flag.isSet(Mesh::FILL_OPP_COORDS));

383
    if (fill_flag.isSet(Mesh::FILL_NEIGH) || fill_opp_coords) {
384 385 386 387
      if (ichild == 0) {
	// Calculation of the neighbour 2, its oppCoords and the
	// cooresponding oppVertex.

Thomas Witkowski's avatar
Thomas Witkowski committed
388
	neighbour[2] = elInfoOld->neighbour[1];
389
	oppVertex[2] = elInfoOld->oppVertex[1];
390

Thomas Witkowski's avatar
Thomas Witkowski committed
391 392
	if (neighbour[2] && fill_opp_coords) {
	  oppCoord[2] = elInfoOld->oppCoord[1];
393
	  neighbourCoord[2] = elInfoOld->neighbourCoord[1];
394
	}
395 396


397 398
	// Calculation of the neighbour 1, its oppCoords and the
	// cooresponding oppVertex.
399 400 401

	if (elem->getFirstChild()  &&
	    elem->getSecondChild()->getFirstChild()  &&
402
	    elem->getSecondChild()->getFirstChild()) {
403

Thomas Witkowski's avatar
Thomas Witkowski committed
404
	  neighbour[1] = elem->getSecondChild()->getSecondChild();
405
	  oppVertex[1] = 2;
406

407
	  if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
408 409 410 411 412 413 414
            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];
415
	    neighbourCoord[1][2] = oppCoord[1];
416 417
	  }
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
418
	  neighbour[1] = elem->getSecondChild();
419
	  oppVertex[1] = 0;
420 421

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

Thomas Witkowski's avatar
Thomas Witkowski committed
424 425 426
	    neighbourCoord[1][0] = elInfoOld->coord[1];
	    neighbourCoord[1][1] = elInfoOld->coord[2];
	    neighbourCoord[1][2] = coord[2];
427 428 429 430
	  }
	}


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

434
	Element *nb = elInfoOld->neighbour[2];
435
	if (nb) {
436 437 438
	  TEST_EXIT_DBG(elInfoOld->oppVertex[2] == 2)
	    ("Fill child %d of element %d (mel %d): Invalid neighbour %d!\n",
	     ichild,
439
	     elInfoOld->getElement()->getIndex(),
440 441 442 443
	     elInfoOld->getMacroElement()->getIndex(),
	     nb->getIndex());

	  TEST_EXIT_DBG(nb->getFirstChild())
444
	    ("Missing first child in element %d!\n", nb->getIndex());
445 446
	  TEST_EXIT_DBG(nb->getSecondChild())
	    ("Missing second child in element %d!\n", nb->getIndex());
447

448 449 450
	  nb = nb->getSecondChild();

	  if (nb->getFirstChild()) {
451
	    oppVertex[0] = 2;
452

453
	    if (fill_opp_coords) {
454
	      if (nb->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
455
		oppCoord[0] = *(nb->getNewCoord());
456
	      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
457
		oppCoord[0].setMidpoint(elInfoOld->neighbourCoord[2][1],
458
					 elInfoOld->neighbourCoord[2][2]);
459
	      }
460

461 462 463
	      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
464
	      neighbourCoord[0][2] = oppCoord[0];
465 466
	    }

467 468
	    nb = nb->getFirstChild();
	  } else {
469
	    oppVertex[0] = 1;
470

471
	    if (fill_opp_coords) {
472
	      oppCoord[0] = elInfoOld->oppCoord[2];
473

474 475 476 477
	      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]);
478 479
	    }
	  }
480
	}
481

Thomas Witkowski's avatar
Thomas Witkowski committed
482
	neighbour[0] = nb;
483 484 485 486
      } else {   /* ichild == 1 */
	// Calculation of the neighbour 2, its oppCoords and the
	// cooresponding oppVertex.

Thomas Witkowski's avatar
Thomas Witkowski committed
487
	neighbour[2] = elInfoOld->neighbour[0];
488
	oppVertex[2] = elInfoOld->oppVertex[0];
489

Thomas Witkowski's avatar
Thomas Witkowski committed
490 491
	if (neighbour[2] && fill_opp_coords) {
	  oppCoord[2] = elInfoOld->oppCoord[0];
492
	  neighbourCoord[2] = elInfoOld->neighbourCoord[0];
493
	}
494

495 496 497 498 499

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

	if (elem->getFirstChild()->getFirstChild()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
500
	  neighbour[0] = elem->getFirstChild()->getFirstChild();
501
	  oppVertex[0] = 2;
502 503

	  if (fill_opp_coords) {
504
            if (elem->getFirstChild()->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
505
	      oppCoord[0] = *(elem->getFirstChild()->getNewCoord());
506
	    } else {
507
	      oppCoord[0].setMidpoint(elInfoOld->coord[0],
Thomas Witkowski's avatar
Thomas Witkowski committed
508
				       elInfoOld->coord[2]);
509
	    }
510

Thomas Witkowski's avatar
Thomas Witkowski committed
511 512 513
	    neighbourCoord[0][0] = coord[2];
	    neighbourCoord[0][1] = coord[1];
	    neighbourCoord[0][2] = oppCoord[0];
514
	  }
515
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
516
	  neighbour[0] = elem->getFirstChild();
517
	  oppVertex[0] = 1;
518 519

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

Thomas Witkowski's avatar
Thomas Witkowski committed
522 523 524
	    neighbourCoord[0][0] = elInfoOld->coord[2];
	    neighbourCoord[0][1] = elInfoOld->coord[0];
	    neighbourCoord[0][2] = coord[2];
525
	  }
526 527 528 529 530
	}

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

531
	Element *nb = elInfoOld->neighbour[2];
532
	if (nb) {
533
	  TEST(elInfoOld->oppVertex[2] == 2)("invalid neighbour\n");
534 535 536
	  TEST((nb = nb->getFirstChild()))("missing child?\n");

	  if (nb->getFirstChild()) {
537
	    oppVertex[1] = 2;
538

539
	    if (fill_opp_coords) {
540
	      if (nb->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
541
		oppCoord[1] = *(nb->getNewCoord());
542
	      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
543
		oppCoord[1].setMidpoint(elInfoOld->neighbourCoord[2][0],
544
					 elInfoOld->neighbourCoord[2][2]);
545
	      }
546

547 548 549
	      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
550
	      neighbourCoord[1][2] = oppCoord[1];
551
	    }
552 553
	    nb = nb->getSecondChild();

554
	  } else {
555
	    oppVertex[1] = 0;
556

557
	    if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
558
	      oppCoord[1] = elInfoOld->oppCoord[2];
559

560
	      neighbourCoord[1][0] = elInfoOld->neighbourCoord[2][2];
561 562 563
	      neighbourCoord[1][1] = elInfoOld->neighbourCoord[2][0];
	      neighbourCoord[1][2].setMidpoint(elInfoOld->neighbourCoord[2][0],
					       elInfoOld->neighbourCoord[2][1]);
564 565 566
	    }
	  }
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
567
	neighbour[1] = nb;
568
      } // if (ichild == 0) {} else
Thomas Witkowski's avatar
Thomas Witkowski committed
569
    } // if (fill_flag.isSet(Mesh::FILL_NEIGH) || fillFlag.isSet(Mesh::FILL_OPP_COORDS))
570

571

572
    if (fill_flag.isSet(Mesh::FILL_BOUND)) {
573
      if (elInfoOld->getBoundary(2))
Thomas Witkowski's avatar
Thomas Witkowski committed
574
	boundary[5] = elInfoOld->getBoundary(2);
575
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
576
	boundary[5] = INTERIOR;
577

578
      if (ichild == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
579 580 581 582 583
	boundary[3] = elInfoOld->getBoundary(5);
	boundary[4] = elInfoOld->getBoundary(3);
	boundary[0] = elInfoOld->getBoundary(2);
	boundary[1] = INTERIOR;
	boundary[2] = elInfoOld->getBoundary(1);
584
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
585 586 587 588 589
	boundary[3] = elInfoOld->getBoundary(4);
	boundary[4] = elInfoOld->getBoundary(5);
	boundary[0] = INTERIOR;
	boundary[1] = elInfoOld->getBoundary(2);
	boundary[2] = elInfoOld->getBoundary(0);
590 591
      }

592
      if (elInfoOld->getProjection(0) &&
593
	  elInfoOld->getProjection(0)->getType() == VOLUME_PROJECTION) {
594

Thomas Witkowski's avatar
Thomas Witkowski committed
595
	projection[0] = elInfoOld->getProjection(0);
596 597
      } else { // boundary projection
	if (ichild == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
598
	  projection[0] = elInfoOld->getProjection(2);
599
	  projection[1] = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
600
	  projection[2] = elInfoOld->getProjection(1);
601
	} else {
602
	  projection[0] = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
603 604
	  projection[1] = elInfoOld->getProjection(2);
	  projection[2] = elInfoOld->getProjection(0);
605 606 607 608 609
	}
      }
    }
  }

610
  double ElInfo2d::calcGrdLambda(DimVec<WorldVector<double> >& grd)
611
  {
612
    FUNCNAME("ElInfo2d::calcGrdLambda()");
613 614

    testFlag(Mesh::FILL_COORDS);
615

616
    double adet = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
617
    int dim = mesh->getDim();
618

Thomas Witkowski's avatar
Thomas Witkowski committed
619
    for (int i = 0; i < dimOfWorld; i++) {
620 621
      e1[i] = coord[1][i] - coord[0][i];
      e2[i] = coord[2][i] - coord[0][i];
622 623
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
624
    if (dimOfWorld == 2) {
625
      double sdet = e1[0] * e2[1] - e1[1] * e2[0];
626 627 628 629
      adet = abs(sdet);

      if (adet < 1.0E-25) {
	MSG("abs(det) = %f\n", adet);
630
	for (int i = 0; i <= dim; i++)
631
	  grd[i].set(0.0);
632 633
      } else {
	double det1 = 1.0 / sdet;
634

635 636 637 638 639 640
	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];
641
      }
642
    } else {
643
      vectorProduct(e1, e2, normal);
644

645
      adet = norm(normal);
646 647 648

      if (adet < 1.0E-15) {
	MSG("abs(det) = %lf\n", adet);
649
	for (int i = 0; i <= dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
650
	  for (int j = 0; j < dimOfWorld; j++)
651
	    grd[i][j] = 0.0;
652
      } else {
653 654
	vectorProduct(e2, normal, grd[1]);
	vectorProduct(normal, e1, grd[2]);
655

656
	double adet2 = 1.0 / (adet * adet);
657

Thomas Witkowski's avatar
Thomas Witkowski committed
658
	for (int i = 0; i < dimOfWorld; i++) {
659 660
	  grd[1][i] *= adet2;
	  grd[2][i] *= adet2;
661 662
	}

663 664 665
	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];
666 667 668 669 670 671
      }
    }

    return adet;
  }

672

673
  int ElInfo2d::worldToCoord(const WorldVector<double>& xy,
674 675
				   DimVec<double>* lambda) const
  {
676
    FUNCNAME("ElInfo::worldToCoord()");
677

678
    TEST_EXIT_DBG(lambda)("lambda must not be NULL\n");
679

Thomas Witkowski's avatar
Thomas Witkowski committed
680
    DimVec<WorldVector<double> > edge(mesh->getDim(), NO_INIT);
681
    WorldVector<double> x;
682

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

685
    for (int j = 0; j < dimOfWorld; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
686
      double x0 = coord[dim][j];
687
      x[j] = xy[j] - x0;
688
      for (int i = 0; i < dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
689
	edge[i][j] = coord[i][j] - x0;
690
    }
691

692 693 694
    double det  = edge<