Liebe Gitlab-Nutzerin, lieber Gitlab-Nutzer,
es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Die Konten der externen Nutzer:innen sind über den Reiter "Standard" erreichbar.
Die Administratoren


Dear Gitlab user,
it is now possible to log in to our service using the ZIH login/LDAP. The accounts of external users can be accessed via the "Standard" tab.
The administrators

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

namespace AMDiS {

17
  double ElInfo2d::mat_d1_left_val[3][3] = {{0.0, 1.0, 0.5}, 
18
  					    {0.0, 0.0, 0.5},
19
  					    {1.0, 0.0, 0.0}}; 
20 21
  mtl::dense2D<double> ElInfo2d::mat_d1_left(mat_d1_left_val);

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

28

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

  double ElInfo2d::mat_d2_left_val[6][6] = {{0.0, 1.0, 0.0, 0.375, -0.125, 0.0}, 
					    {0.0, 0.0, 0.0, -0.125, -0.125, 0.0}, 
					    {1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
					    {0.0, 0.0, 0.0, 0.0, 0.5, 0.0},
					    {0.0, 0.0, 0.0, 0.0, 0.5, 1.0},
					    {0.0, 0.0, 1.0, 0.75, 0.25, 0.0}};
  mtl::dense2D<double> ElInfo2d::mat_d2_left(mat_d2_left_val);

  double ElInfo2d::mat_d2_right_val[6][6] = {{0.0, 0.0, 0.0, -0.125, -0.125, 0.0}, 
					     {1.0, 0.0, 0.0, -0.125, 0.375, 0.0}, 
					     {0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
					     {0.0, 0.0, 0.0, 0.5, 0.0, 1.0},
					     {0.0, 0.0, 0.0, 0.5, 0.0, 0.0},
					     {0.0, 0.0, 1.0, 0.25, 0.75, 0.0}};
  mtl::dense2D<double> ElInfo2d::mat_d2_right(mat_d2_right_val);


47

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

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



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

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


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

113

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

117

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

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

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

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

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

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

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

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


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

	    // Search for the next neighbour. In many cases, the neighbour element 
	    // may be refinemed in a way, such that there is no new vertex on the 
	    // common edge. This situation is shown in the following picture: 
	    //
	    //               /|\
	    //              / | \
	    //             /  |  \
	    //            /\  |   \
	    //           /  \ |    \ 
	    //          /    \|     \
	    //          -------------
	    //
	    //            nb     el
	    //
	    // Note that we know (because of the last if statement), that the 
	    // neighbour element has children and the common edge is not the 
	    // refinement edge, which has always the number 2, of our element.


174
	    if (edgeNo == 0) {
175 176 177 178 179 180 181 182 183 184 185 186
	      // The situation is as follows:
	      //
	      //          -------
	      //          \    /|\
	      //           \  / | \
	      //            \/  |  \
	      //             \  |   \
	      //              \ |    \ 
	      //               \|     \
	      //                -------
	      //
	      //            nb     el
187
              //
188 189 190 191
	      // That means, the edge 0 of the same level neighbour is the common
	      // edge, i.e., the direct neighbour is the second child of the same
	      // level neighbour.

Thomas Witkowski's avatar
Thomas Witkowski committed
192
	      nb = neighbour[i] = nb->getSecondChild();
193
	    } else {
194 195
	      // The situation is as shown in the picture above. So the next
	      // neighbour is the first child of the same level neighbour element.
Thomas Witkowski's avatar
Thomas Witkowski committed
196
	      nb = neighbour[i] = nb->getFirstChild();
197 198
	    }

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

	    if (fill_opp_coords) {
204
	      if (nb->isNewCoordSet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
205
		oppCoord[i] = *(nb->getNewCoord());
206
	      } else {
207 208 209
		// In both cases, that are shown in the pictures above, the opp
		// vertex of the neighbour edge is the midpoint of the vertex 0
		// and vertex 1 of the same level neighbour element.
Thomas Witkowski's avatar
Thomas Witkowski committed
210
		oppCoord[i] = (macroNeighbour->coord[0] + 
211 212
				macroNeighbour->coord[1]) * 0.5;
	      }
213 214 215
	      
	      switch (i) {
	      case 0:
216 217 218 219
		// The common edge is the edge 0 of this element.

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

263 264 265
		// I've deleted here some code, be I think that this case is not
		// possible. If an error occurs in this line, please check AMDiS
		// revision <= 476 at the same position.
266 267
		//		ERROR_EXIT("Should not happen!\n");

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

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

	    // In this case, we know that the common edge is the refinement edge.
	    // This makes everything much more simpler, because we know that the
	    // next neighbour is equal to the samel level neighbour. If the same
	    // level neighbour would be refinement, also this element must to be 
	    // refinement, because they share the refinement edge.

298
	    if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
299
	      oppCoord[i] = macroNeighbour->coord[edgeNo];
300
	      neighbourCoord[i] = macroNeighbour->coord;	      
301
	    }
302
	  }
303
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
304
	  neighbour[i] = NULL;
305
        }
306
      }
307
    }
308
    
Thomas Witkowski's avatar
Thomas Witkowski committed
309 310 311 312 313
    if (fillFlag.isSet(Mesh::FILL_BOUND)) {   
      for (int i = 0; i < element->getGeo(BOUNDARY); i++)
	boundary[i] = mel->getBoundary(i);
      for (int i = 0; i < element->getGeo(PROJECTION); i++)
	projection[i] = mel->getProjection(i);      
314 315 316 317 318 319 320 321
    }
  }


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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
366
	neighbour[2] = elInfoOld->neighbour[1];
367
	oppVertex[2] = elInfoOld->oppVertex[1];
368
	
Thomas Witkowski's avatar
Thomas Witkowski committed
369 370
	if (neighbour[2] && fill_opp_coords) {
	  oppCoord[2] = elInfoOld->oppCoord[1];
371
	  neighbourCoord[2] = elInfoOld->neighbourCoord[1];
372
	}
373 374 375 376 377 378 379 380 381
	
	
	// Calculation of the neighbour 1, its oppCoords and the
	// cooresponding oppVertex.
	
	if (elem->getFirstChild()  &&  
	    elem->getSecondChild()->getFirstChild()  &&  
	    elem->getSecondChild()->getFirstChild()) {
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
382
	  neighbour[1] = elem->getSecondChild()->getSecondChild();
383
	  oppVertex[1] = 2;
384 385
	  
	  if (fill_opp_coords) {
Thomas Witkowski's avatar
Thomas Witkowski committed
386 387 388 389 390 391 392 393
            if (elem->getSecondChild()->isNewCoordSet())
	      oppCoord[1] = *(elem->getSecondChild()->getNewCoord());
	    else
	      oppCoord[1].setMidpoint(elInfoOld->coord[1], elInfoOld->coord[2]);

	    neighbourCoord[1][0] = coord[0];
	    neighbourCoord[1][1] = coord[2];
	    neighbourCoord[1][2] = oppCoord[1];  
394 395
	  }
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
396
	  neighbour[1] = elem->getSecondChild();
397
	  oppVertex[1] = 0;
398 399

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

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


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

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

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

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

439 440 441
	      neighbourCoord[0][0].setMidpoint(elInfoOld->neighbourCoord[2][0],
						elInfoOld->neighbourCoord[2][1]);
	      neighbourCoord[0][1] = elInfoOld->neighbourCoord[2][1];
Thomas Witkowski's avatar
Thomas Witkowski committed
442
	      neighbourCoord[0][2] = oppCoord[0];
443 444 445 446
	    }	   
 
	    nb = nb->getFirstChild();
	  } else {
447
	    oppVertex[0] = 1;
448

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

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

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

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

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

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

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

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

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

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

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

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

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

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

525 526 527
	      neighbourCoord[1][0] = elInfoOld->neighbourCoord[2][0];
	      neighbourCoord[1][1].setMidpoint(elInfoOld->neighbourCoord[2][0],
					       elInfoOld->neighbourCoord[2][1]);
Thomas Witkowski's avatar
Thomas Witkowski committed
528
	      neighbourCoord[1][2] = oppCoord[1];
529
	    }
530 531
	    nb = nb->getSecondChild();

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

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

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

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

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

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

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

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

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

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

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

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

623
      adet = norm(normal);
624 625 626

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

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

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

    return adet;
  }

650 651

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

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

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

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

664
    for (int j = 0; j < dimOfWorld; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
665
      double x0 = coord[dim][j];
666
      x[j] = xy[j] - x0;
667
      for (int i = 0; i < dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
668
	edge[i][j] = coord[i][j] - x0;
669 670
    }
  
671 672 673
    double det  = edge[0][0] * edge[1][1] - edge[0][1] * edge[1][0]; 
    double det0 =       x[0] * edge[1][1] -       x[1] * edge[1][0]; 
    double det1 = edge[0][0] * x[1]       - edge[0][1] * x[0]; 
674 675 676

    if (abs(det) < DBL_TOL) {
      ERROR("det = %le; abort\n", det);
677 678
      for (int i = 0; i <= dim; i++) 
	(*lambda)[i] = 1.0/dim;
679 680 681 682 683 684 685 686 687
      return 0;
    }

    (*lambda)[0] = det0 / det;
    (*lambda)[1] = det1 / det;
    (*lambda)[2] = 1.0 - (*lambda)[0] - (*lambda)[1];

    int k = -1;
    double lmin = 0.0;
688
    for (int i = 0; i <= dim; i++) {
689 690 691 692 693 694 695 696 697 698 699 700
      if ((*lambda)[i] < -1.E-5) {
	if ((*lambda)[i] < lmin) {
	  k = i;
	  lmin = (*lambda)[i];
	}
      }
    }

    return k;
  }


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

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

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

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

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

725
    double det = norm(&normal);
726

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

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

734

735 736 737 738 739
  /****************************************************************************/
  /*  calculate the normal of the element for dim of world = dim + 1          */
  /*  return the absulute value of the determinant from the                   */
  /*  transformation to the reference element                                 */
  /****************************************************************************/
Thomas Witkowski's avatar
Thomas Witkowski committed
740
  double ElInfo2d::getElementNormal(WorldVector<double> &elementNormal) const
741
  {
742
    FUNCNAME("ElInfo::getElementNormal()");
743

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

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

    vectorProduct(e0, e1, elementNormal);

752
    double det = norm(&elementNormal);
753

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

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