ElInfo3d.cc 19.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#include "ElInfo3d.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
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
  double ElInfo3d::mat_d1_left_val[4][4] = {{1.0, 0.0, 0.0, 0.5}, 
					    {0.0, 0.0, 0.0, 0.5},
					    {0.0, 1.0, 0.0, 0.0},
					    {0.0, 0.0, 1.0, 0.0}};
  mtl::dense2D<double> ElInfo3d::mat_d1_left(mat_d1_left_val);

  
  double ElInfo3d::mat_d1_l0_right_val[4][4] = {{0.0, 0.0, 0.0, 0.5}, 
						{1.0, 0.0, 0.0, 0.5},
						{0.0, 0.0, 1.0, 0.0},
						{0.0, 1.0, 0.0, 0.0}};
  double ElInfo3d::mat_d1_l12_right_val[4][4] = {{0.0, 0.0, 0.0, 0.5}, 
						 {1.0, 0.0, 0.0, 0.5},
						 {0.0, 1.0, 0.0, 0.0},
						 {0.0, 0.0, 1.0, 0.0}};
  mtl::dense2D<double> ElInfo3d::mat_d1_l0_right(mat_d1_l0_right_val);
  mtl::dense2D<double> ElInfo3d::mat_d1_l12_right(mat_d1_l12_right_val);


36
37
  void ElInfo3d::fillMacroInfo(const MacroElement * mel)
  {
38
39
    FUNCNAME("ElInfo3d::fillMacroInfo()");
    Element *nb;
40
    MacroElement *mnb;
41
    Flag fill_opp_coords;
42

Thomas Witkowski's avatar
Thomas Witkowski committed
43
44
45
    macroElement  = const_cast<MacroElement*>( mel);
    element  = const_cast<Element*>( mel->getElement());
    parent = NULL;
46
    level = 0;
47
    elType = const_cast<MacroElement*>(mel)->getElType();
48

Thomas Witkowski's avatar
Thomas Witkowski committed
49
    int vertices = mesh->getGeo(VERTEX);
50

Thomas Witkowski's avatar
Thomas Witkowski committed
51
52
53
    if (fillFlag.isSet(Mesh::FILL_COORDS) || 
	fillFlag.isSet(Mesh::FILL_DET) ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
54
      for (int i = 0; i < vertices; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
55
	coord[i] = mel->coord[i];
56
57
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
60
61
    if (fillFlag.isSet(Mesh::FILL_OPP_COORDS) || 
        fillFlag.isSet(Mesh::FILL_NEIGH)) {
62

Thomas Witkowski's avatar
Thomas Witkowski committed
63
      fill_opp_coords.setFlags(fillFlag & Mesh::FILL_OPP_COORDS);
64
      for (int i = 0; i < neighbours; i++) {
65
	if ((mnb = const_cast<MacroElement*>( mel->getNeighbour(i)))) {
Thomas Witkowski's avatar
Thomas Witkowski committed
66
67
	  neighbour[i] = const_cast<Element*>( mel->getNeighbour(i)->getElement());
	  nb = const_cast<Element*>( neighbour[i]);
68
	  int k;
69
	  k = oppVertex[i] = mel->getOppVertex(i);
70
71

	  if (nb->getChild(0) && (k < 2)) {   /*make nb nearest element.*/
72
	    if (k == 1) {
Thomas Witkowski's avatar
Thomas Witkowski committed
73
74
	      neighbour[i]      = const_cast<Element*>( nb->getChild(0));
	      nb = const_cast<Element*>( neighbour[i]);
75
	    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
76
77
	      neighbour[i]      = const_cast<Element*>( nb->getChild(1));
	      nb = const_cast<Element*>( neighbour[i]);
78
	    }
79
	    k = oppVertex[i] = 3;
80
81
82
	    if (fill_opp_coords.isAnySet()) {
	      /* always edge between vertices 0 and 1 is bisected! */
	      if (mnb->getElement()->isNewCoordSet())
Thomas Witkowski's avatar
Thomas Witkowski committed
83
		oppCoord[i] = *(mnb->getElement()->getNewCoord());
84
	      else
Thomas Witkowski's avatar
Thomas Witkowski committed
85
		oppCoord[i] = (mnb->coord[0] + mnb->coord[1]) * 0.5;
86
87
	    }
	  } else {
88
	    if  (fill_opp_coords.isAnySet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
89
	      oppCoord[i] = mnb->coord[k];
90
91
	    }
	  }
92
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
93
	  neighbour[i] = NULL;
94
95
96
97
	}
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
98
99
100
    if (fillFlag.isSet(Mesh::FILL_BOUND)) {
      for (int i = 0; i < element->getGeo(BOUNDARY); i++)
	boundary[i] = mel->getBoundary(i);     
101

Thomas Witkowski's avatar
Thomas Witkowski committed
102
103
      for (int i = 0; i < element->getGeo(PROJECTION); i++)
	projection[i] = mel->getProjection(i);
104
    }
105

Thomas Witkowski's avatar
Thomas Witkowski committed
106
    if (fillFlag.isSet(Mesh::FILL_ORIENTATION)) {
107
108
109
      WorldVector<WorldVector<double> > a;
      double s;

110
111
      for (int i = 0; i < 3; i++) {
	a[i] = mel->coord[i + 1];
112
113
114
115
116
117
118
119
120
121
122
123
124
125
	a[i] -= mel->coord[0];
      }

      s = (a[0][1] * a[1][2] - a[0][2] * a[1][1]) * a[2][0]
	+ (a[0][2] * a[1][0] - a[0][0] * a[1][2]) * a[2][1]
	+ (a[0][0] * a[1][1] - a[0][1] * a[1][0]) * a[2][2];

      if (s >= 0)
	orientation = 1;
      else
	orientation = -1;
    }
  }

126

Thomas Witkowski's avatar
Thomas Witkowski committed
127
  double ElInfo3d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
128
  {
129
130
    FUNCNAME("ElInfo3d::calcGrdLambda()");

Thomas Witkowski's avatar
Thomas Witkowski committed
131
    TEST_EXIT_DBG(dimOfWorld == 3)
132
133
      ("dim != dim_of_world ! use parametric elements!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
134
135
136
137
    WorldVector<double> *e1 = &tmpWorldVecs[0];
    WorldVector<double> *e2 = &tmpWorldVecs[1];
    WorldVector<double> *e3 = &tmpWorldVecs[2];
    WorldVector<double> *v0 = &tmpWorldVecs[3];
Thomas Witkowski's avatar
Thomas Witkowski committed
138

139
140
    double det, adet;
    double a11, a12, a13, a21, a22, a23, a31, a32, a33;
141
142
143

    testFlag(Mesh::FILL_COORDS);

Thomas Witkowski's avatar
Thomas Witkowski committed
144
    *v0 = coord[0];
145
    for (int i = 0; i < 3; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
146
147
148
      (*e1)[i] = coord[1][i] - (*v0)[i];
      (*e2)[i] = coord[2][i] - (*v0)[i];
      (*e3)[i] = coord[3][i] - (*v0)[i];
149
150
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
151
152
153
    det = (*e1)[0] * ((*e2)[1] * (*e3)[2] - (*e2)[2] * (*e3)[1])
      - (*e1)[1] * ((*e2)[0] * (*e3)[2] - (*e2)[2] * (*e3)[0])
      + (*e1)[2] * ((*e2)[0] * (*e3)[1] - (*e2)[1] * (*e3)[0]);
154
155
156

    adet = abs(det);

157
158
159
160
161
162
163
    if (adet < 1.0E-25) {
      MSG("abs(det) = %f\n",adet);
      for (int i = 0; i < 4; i++)
	for (int j = 0; j < 3; j++)
	  grd_lam[i][j] = 0.0;
    } else {
      det = 1.0 / det;
Thomas Witkowski's avatar
Thomas Witkowski committed
164
165
166
167
168
169
170
171
172
      a11 = ((*e2)[1] * (*e3)[2] - (*e2)[2] * (*e3)[1]) * det;    /* (a_ij) = A^{-T} */
      a12 = ((*e2)[2] * (*e3)[0] - (*e2)[0] * (*e3)[2]) * det;
      a13 = ((*e2)[0] * (*e3)[1] - (*e2)[1] * (*e3)[0]) * det;
      a21 = ((*e1)[2] * (*e3)[1] - (*e1)[1] * (*e3)[2]) * det;
      a22 = ((*e1)[0] * (*e3)[2] - (*e1)[2] * (*e3)[0]) * det;
      a23 = ((*e1)[1] * (*e3)[0] - (*e1)[0] * (*e3)[1]) * det;
      a31 = ((*e1)[1] * (*e2)[2] - (*e1)[2] * (*e2)[1]) * det;
      a32 = ((*e1)[2] * (*e2)[0] - (*e1)[0] * (*e2)[2]) * det;
      a33 = ((*e1)[0] * (*e2)[1] - (*e1)[1] * (*e2)[0]) * det;
173
174
175
176
177
178
179
180
181
182
183

      grd_lam[1][0] = a11;
      grd_lam[1][1] = a12;
      grd_lam[1][2] = a13;
      grd_lam[2][0] = a21;
      grd_lam[2][1] = a22;
      grd_lam[2][2] = a23;
      grd_lam[3][0] = a31;
      grd_lam[3][1] = a32;
      grd_lam[3][2] = a33;

Thomas Witkowski's avatar
Thomas Witkowski committed
184
185
186
      grd_lam[0][0] = -grd_lam[1][0] - grd_lam[2][0] - grd_lam[3][0];
      grd_lam[0][1] = -grd_lam[1][1] - grd_lam[2][1] - grd_lam[3][1];
      grd_lam[0][2] = -grd_lam[1][2] - grd_lam[2][2] - grd_lam[3][2];
187
    }
188
189
190
191

    return adet;
  }

192

193
194
195
  const int ElInfo3d::worldToCoord(const WorldVector<double>& xy,
				   DimVec<double>* lambda) const
  {
196
197
    FUNCNAME("ElInfo::worldToCoord()");

Thomas Witkowski's avatar
Thomas Witkowski committed
198
    DimVec<WorldVector<double> > edge(mesh->getDim(), NO_INIT);
199
200
201
    WorldVector<double> x;
    double  x0, det, det0, det1, det2;
  
Thomas Witkowski's avatar
Thomas Witkowski committed
202
    static DimVec<double> vec(mesh->getDim(), NO_INIT);
203

204
    TEST_EXIT_DBG(lambda)("lambda must not be NULL\n");
205

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

208
    TEST_EXIT_DBG(dim == dimOfWorld)("dim!=dimOfWorld not yet implemented\n");
209
210
211
212
213
214
215
    
    /*  wir haben das gleichungssystem zu loesen: */
    /*       ( q1x q2x q3x)  (lambda1)     (qx)      */
    /*       ( q1y q2y q3y)  (lambda2)  =  (qy)      */
    /*       ( q1z q2z q3z)  (lambda3)     (qz)      */
    /*      mit qi=pi-p3, q=xy-p3                 */

216
    for (int j = 0; j < dimOfWorld; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
217
      x0 = coord[dim][j];
218
      x[j] = xy[j] - x0;
219
220

      for (int i = 0; i < dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
221
	edge[i][j] = coord[i][j] - x0;
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
    }

    det =  edge[0][0] * edge[1][1] * edge[2][2]
      + edge[0][1] * edge[1][2] * edge[2][0]
      + edge[0][2] * edge[1][0] * edge[2][1]
      - edge[0][2] * edge[1][1] * edge[2][0]
      - edge[0][0] * edge[1][2] * edge[2][1]
      - edge[0][1] * edge[1][0] * edge[2][2];
    det0 =       x[0] * edge[1][1] * edge[2][2]
      +       x[1] * edge[1][2] * edge[2][0]
      +       x[2] * edge[1][0] * edge[2][1]
      -       x[2] * edge[1][1] * edge[2][0]
      -       x[0] * edge[1][2] * edge[2][1]
      -       x[1] * edge[1][0] * edge[2][2];
    det1 = edge[0][0] *       x[1] * edge[2][2]
      + edge[0][1] *       x[2] * edge[2][0]
      + edge[0][2] *       x[0] * edge[2][1]
      - edge[0][2] *       x[1] * edge[2][0]
      - edge[0][0] *       x[2] * edge[2][1]
      - edge[0][1] *       x[0] * edge[2][2];
    det2 = edge[0][0] * edge[1][1] *       x[2]
      + edge[0][1] * edge[1][2] *       x[0]
      + edge[0][2] * edge[1][0] *       x[1]
      - edge[0][2] * edge[1][1] *       x[0]
      - edge[0][0] * edge[1][2] *       x[1]
      - edge[0][1] * edge[1][0] *       x[2];
  
    if (abs(det) < DBL_TOL) {
      ERROR("det = %le; abort\n", det);
251

252
      for (int i = 0; i <= dim; i++)
253
254
	(*lambda)[i] = 1.0 / dim;

255
256
257
258
259
260
261
262
263
264
      return 0;
    }

    (*lambda)[0] = det0 / det;
    (*lambda)[1] = det1 / det;
    (*lambda)[2] = det2 / det;
    (*lambda)[3] = 1.0 - (*lambda)[0] - (*lambda)[1] - (*lambda)[2];
  
    int k = -1;
    double lmin = 0.0;
265
266

    for (int i = 0; i <= dim; i++) {
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
      if ((*lambda)[i] < -1.E-5) {
	if ((*lambda)[i] < lmin) {
	  k = i;
	  lmin = (*lambda)[i];
	}
      }
    }

    return k;
  }


  /****************************************************************************/
  /*   update EL_INFO structure after refinement (of some neighbours)	    */
  /****************************************************************************/

  void ElInfo3d::update()
  {
285
    FUNCNAME("ElInfo::update()");
286

Thomas Witkowski's avatar
Thomas Witkowski committed
287
288
    int neighbours = mesh->getGeo(NEIGH);
    int vertices = mesh->getGeo(VERTEX);
289
  
Thomas Witkowski's avatar
Thomas Witkowski committed
290
    if (fillFlag.isSet(Mesh::FILL_NEIGH) || fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {
291
      Tetrahedron *nb;
Thomas Witkowski's avatar
Thomas Witkowski committed
292
      Flag fill_opp_coords = fillFlag & Mesh::FILL_OPP_COORDS;
293
294
      
      for (int ineigh = 0; ineigh < neighbours; ineigh++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
295
	if ((nb = dynamic_cast<Tetrahedron*>(const_cast<Element*>(neighbour[ineigh])))) {
296
	  int ov = oppVertex[ineigh];
297
298
299
300
	  if (ov < 2 && nb->getFirstChild()) {
	    if (fill_opp_coords != Flag(0)) {
	      int k = -1;
	      for (int j = 0; j < vertices; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
301
		if (element->getDOF(j) == nb->getDOF(1 - ov)) 
302
303
304
		  k = j;
	      
	      if (k == -1) {
305
		for (int j = 0; j < vertices; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
306
		  if (mesh->associated(element->getDOF(j, 0), nb->getDOF(1 - ov, 0)))
307
		    k = j;
308
	      }
309
	      TEST_EXIT_DBG(k >= 0)("neighbour dof not found\n");
310
311
	      
	      if (nb->isNewCoordSet())
Thomas Witkowski's avatar
Thomas Witkowski committed
312
		oppCoord[ineigh] = *(nb->getNewCoord());
313
	      else
Thomas Witkowski's avatar
Thomas Witkowski committed
314
		for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
315
		  oppCoord[ineigh][j] = (oppCoord[ineigh][j] + coord[k][j]) / 2;
316
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
317
	    neighbour[ineigh] = dynamic_cast<Tetrahedron*>(const_cast<Element*>(nb->getChild(1-ov)));
318
	    oppVertex[ineigh] = 3;
319
320
321
	  }
	}
      }
322
    }
323
324
325
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
326
  double ElInfo3d::getNormal(int face, WorldVector<double> &normal)
327
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
328
329
    FUNCNAME("ElInfo3d::getNormal()");

330
    double det = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
331

Thomas Witkowski's avatar
Thomas Witkowski committed
332
333
334
    WorldVector<double> *e0 = &tmpWorldVecs[0];
    WorldVector<double> *e1 = &tmpWorldVecs[1];
    WorldVector<double> *e2 = &tmpWorldVecs[2];
335

Thomas Witkowski's avatar
Thomas Witkowski committed
336
    if (dimOfWorld == 3) {
337
338
339
      int i0 = (face + 1) % 4;
      int i1 = (face + 2) % 4;
      int i2 = (face + 3) % 4;
340

Thomas Witkowski's avatar
Thomas Witkowski committed
341
      for (int i = 0; i < dimOfWorld; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
342
343
344
	(*e0)[i] = coord[i1][i] - coord[i0][i];
	(*e1)[i] = coord[i2][i] - coord[i0][i];
	(*e2)[i] = coord[face][i] - coord[i0][i];
345
346
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
347
      vectorProduct(*e0, *e1, normal);
348

Thomas Witkowski's avatar
Thomas Witkowski committed
349
      if ((*e2 * normal) < 0.0)
Thomas Witkowski's avatar
Thomas Witkowski committed
350
	for (int i = 0; i < dimOfWorld; i++)
351
352
353
	  normal[i] = -normal[i];

      det = norm(&normal);
354
      TEST_EXIT_DBG(det > 1.e-30)("det = 0 on face %d\n", face);
355
356
357
358
359

      normal[0] /= det;
      normal[1] /= det;
      normal[2] /= det;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
360
      MSG("not implemented for DIM_OF_WORLD = %d in 3d\n", dimOfWorld);
361
362
363
364
365
366
    }

    return(det);
  }


367
  void ElInfo3d::fillElInfo(int ichild, const ElInfo *elInfoOld)
368
  {
369
370
371
372
    FUNCNAME("ElInfo3d::fillElInfo()");

    int ochild = 0;             /* index of other child = 1-ichild */
    int *cv = NULL;             /* cv = child_vertex[el_type][ichild] */
373
    const int (*cvg)[4] = NULL;     /* cvg = child_vertex[el_type] */
374
    int *ce;                    /* ce = child_edge[el_type][ichild] */
375
    Element *nb, *nbk;
Thomas Witkowski's avatar
Thomas Witkowski committed
376
377
    Element *el_old = elInfoOld->element;
    Flag fillFlag_local = elInfoOld->fillFlag;
378
    DegreeOfFreedom *dof;
379
    int ov = -1;
380
    Mesh *mesh = elInfoOld->getMesh();
381

Thomas Witkowski's avatar
Thomas Witkowski committed
382
    TEST_EXIT_DBG(el_old->getChild(0))("missing child?\n"); 
383

Thomas Witkowski's avatar
Thomas Witkowski committed
384
385
386
387
    element = const_cast<Element*>( el_old->getChild(ichild));
    macroElement = elInfoOld->macroElement;
    fillFlag = fillFlag_local;
    parent  = el_old;
388
    level = elInfoOld->level + 1;
389
390
    iChild = ichild;
    int el_type_local = (dynamic_cast<const ElInfo3d*>(elInfoOld))->getType();
391
    elType = (el_type_local + 1) % 3;
392

Thomas Witkowski's avatar
Thomas Witkowski committed
393
    TEST_EXIT_DBG(element)("missing child %d?\n", ichild);
394

Thomas Witkowski's avatar
Thomas Witkowski committed
395
    if (fillFlag_local.isAnySet()) {
396
      cvg = Tetrahedron::childVertex[el_type_local];
397
      cv = const_cast<int*>(cvg[ichild]);
398
      ochild = 1 - ichild;
399
400
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
401
402
403
    if (fillFlag_local.isSet(Mesh::FILL_COORDS) || 
	fillFlag.isSet(Mesh::FILL_DET) ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
404
      for (int i = 0; i < 3; i++) {
405
	for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
406
	  coord[i][j] = elInfoOld->coord[cv[i]][j];
407
      }
408
      if (el_old->getNewCoord()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
409
	coord[3] = *(el_old->getNewCoord());
410
      } else {
411
	for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
412
	  coord[3][j] = (elInfoOld->coord[0][j] + elInfoOld->coord[1][j]) / 2;
413
414
      }
    }
415

Thomas Witkowski's avatar
Thomas Witkowski committed
416
417
    if (fillFlag_local.isSet(Mesh::FILL_NEIGH) || 
	fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {
418

Thomas Witkowski's avatar
Thomas Witkowski committed
419
420
      FixVec<Element*, NEIGH> *neigh_local = &neighbour;
      const FixVec<Element*, NEIGH> *neigh_old = &elInfoOld->neighbour;
421
422

      Flag fill_opp_coords;
Thomas Witkowski's avatar
Thomas Witkowski committed
423
      fill_opp_coords.setFlags(fillFlag_local & Mesh::FILL_OPP_COORDS);
424
425
426
427
428
429
430
431
432
433
      
      /*----- nb[0] is other child --------------------------------------------*/
      
      /*    if (nb = el_old->child[ochild])   old version  */
      if (el_old->getChild(0)  &&  
	  (nb = const_cast<Element*>( el_old->getChild(ochild)))) {
	
	if (nb->getChild(0)) {         /* go down one level for direct neighbour */
	  if (fill_opp_coords.isAnySet()) {
	    if (nb->getNewCoord()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
434
	      oppCoord[0]= *(nb->getNewCoord());
435
436
	    } else {
	      int k = cvg[ochild][1];
437
	      for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
438
		oppCoord[0][j] = (elInfoOld->coord[ochild][j] + elInfoOld->coord[k][j]) / 2;
439
440
	    }
	  }
441
	  (*neigh_local)[0] = const_cast<Element*>( nb->getChild(1));
442
	  oppVertex[0] = 3;
443
	} else {
444
445
	  if (fill_opp_coords.isAnySet())
	    for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
446
	      oppCoord[0][j] = elInfoOld->coord[ochild][j];
447

448
	  (*neigh_local)[0] = nb;
449
	  oppVertex[0] = 0;
450
	}
451
452
453
454
      } else {
	ERROR_EXIT("no other child");
	(*neigh_local)[0] = NULL;
      }
455
456


457
458
459
460
      /*----- nb[1],nb[2] are childs of old neighbours nb_old[cv[i]] ----------*/
      
      for (int i = 1; i < 3; i++) {
	if ((nb = const_cast<Element*>( (*neigh_old)[cv[i]]))) {
461
	  TEST_EXIT_DBG(nb->getChild(0))("nonconforming triangulation\n");
462
	  
463
464
465
466
467
	  int k;
	  for (k = 0; k < 2; k++) { /* look at both childs of old neighbour */
	    nbk = const_cast<Element*>( nb->getChild(k));
	    if (nbk->getDOF(0) == el_old->getDOF(ichild)) {
	      /* opp. vertex */
468
	      dof = const_cast<int*>(nb->getDOF(elInfoOld->oppVertex[cv[i]])); 
469
470
471
472
473
474
	      
	      if (dof == nbk->getDOF(1)) {
		ov = 1;
		if (nbk->getChild(0)) {
		  if (fill_opp_coords.isAnySet()) {
		    if (nbk->getNewCoord()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
475
		      oppCoord[i] = *(nbk->getNewCoord());
476
		    } else {
477
		      for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
478
479
			oppCoord[i][j] = (elInfoOld->oppCoord[cv[i]][j] + 
					   elInfoOld->coord[ichild][j]) / 2;		      
480
481
482
		    }
		  }
		  (*neigh_local)[i] = nbk->getChild(0);
483
		  oppVertex[i] = 3;
484
		  break;
485
		}
486
487
488
489
490
491
	      } else {
		if (dof != nbk->getDOF(2)) { 
		  ov = -1; 
		  break; 
		}
		ov = 2;
492
493
	      }

494
495
	      if (fill_opp_coords.isAnySet())
		for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
496
		  oppCoord[i][j] = elInfoOld->oppCoord[cv[i]][j];
497

498
	      (*neigh_local)[i] = nbk;
499
	      oppVertex[i] = ov;
500
	      break;
501
	    }
502
503
504
505
506
507
508
509
510
511
	    
	  } /* end for k */

	  // periodic ?
	  if (k == 2 || ov == -1) {
	    for (k = 0; k < 2; k++) {  /* look at both childs of old neighbour */
	      nbk = const_cast<Element*>( nb->getChild(k));
	      if (nbk->getDOF(0) == el_old->getDOF(ichild) ||
		  mesh->associated(nbk->getDOF(0, 0), el_old->getDOF(ichild, 0))) {
		/* opp. vertex */
512
		dof = const_cast<int*>(nb->getDOF(elInfoOld->oppVertex[cv[i]])); 
513
514
515
516
517
518
519
		
		if (dof == nbk->getDOF(1) || 
		    mesh->associated(dof[0], nbk->getDOF(1, 0))) {
		  ov = 1;
		  if (nbk->getChild(0)) {
		    if (fill_opp_coords.isAnySet()) {
		      if (nbk->getNewCoord()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
520
			oppCoord[i] = *(nbk->getNewCoord());
521
		      } else {
522
			for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
523
524
			  oppCoord[i][j] = (elInfoOld->oppCoord[cv[i]][j] + 
					    elInfoOld->coord[ichild][j]) / 2;
525
526
527
		      }
		    }
		    (*neigh_local)[i] = nbk->getChild(0);
528
		    oppVertex[i] = 3;
529
530
531
		    break;
		  }
		} else {
532
533
		  TEST_EXIT_DBG(dof == nbk->getDOF(2) || 
				mesh->associated(dof[0], nbk->getDOF(2, 0)))
534
535
536
		    ("opp_vertex not found\n");
		  ov = 2;
		}
537

538
539
		if (fill_opp_coords.isAnySet())
		  for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
540
		    oppCoord[i][j] = elInfoOld->oppCoord[cv[i]][j];
541

542
		(*neigh_local)[i] = nbk;
543
		oppVertex[i] = ov;
544
545
546
547
		break;
	      }
	      
	    } /* end for k */
548

549
	    TEST_EXIT_DBG(k < 2)("child not found with vertex\n");
550
551
552
	  }
	} else {
	  (*neigh_local)[i] = NULL;
553
	}
554
555
556
557
558
559
      }  /* end for i */
      
      
      /*----- nb[3] is old neighbour neigh_old[ochild] ------------------------*/
      
      if (((*neigh_local)[3] = (*neigh_old)[ochild])) {
560
	oppVertex[3] = elInfoOld->oppVertex[ochild];
561

562
563
	if (fill_opp_coords.isAnySet())
	  for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
564
	    oppCoord[3][j] = elInfoOld->oppCoord[ochild][j];
565
566
      }
    }
567

Thomas Witkowski's avatar
Thomas Witkowski committed
568
    if (fillFlag_local.isSet(Mesh::FILL_BOUND)) {
569
      for (int i = 0; i < 3; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
570
	boundary[10 + i] = elInfoOld->getBoundary(10 + cv[i]);
571
      
Thomas Witkowski's avatar
Thomas Witkowski committed
572
      boundary[13] = elInfoOld->getBoundary(4);
573
      
Thomas Witkowski's avatar
Thomas Witkowski committed
574
575
576
577
      boundary[0] = INTERIOR;
      boundary[1] = elInfoOld->getBoundary(cv[1]);
      boundary[2] = elInfoOld->getBoundary(cv[2]);
      boundary[3] = elInfoOld->getBoundary(ochild);
578
      
Thomas Witkowski's avatar
Thomas Witkowski committed
579
      int geoFace = mesh->getGeo(FACE);
580

581
      ce = const_cast<int*>(Tetrahedron::childEdge[el_type_local][ichild]);
582
      for (int iedge = 0; iedge < 4; iedge++)
Thomas Witkowski's avatar
Thomas Witkowski committed
583
584
585
	boundary[geoFace + iedge] = elInfoOld->getBoundary(geoFace + ce[iedge]);      
      for (int iedge = 4; iedge < 6; iedge++)
	boundary[geoFace + iedge] = elInfoOld->getBoundary(5 - cv[iedge - 3]);
586

587
588
      if (elInfoOld->getProjection(0) &&
	  elInfoOld->getProjection(0)->getType() == VOLUME_PROJECTION) {
589
	
Thomas Witkowski's avatar
Thomas Witkowski committed
590
	projection[0] = elInfoOld->getProjection(0);      
591
      } else { // boundary projection
Thomas Witkowski's avatar
Thomas Witkowski committed
592
593
594
595
	projection[0] = NULL;
	projection[1] = elInfoOld->getProjection(cv[1]);
	projection[2] = elInfoOld->getProjection(cv[2]);
	projection[3] = elInfoOld->getProjection(ochild);
596
	
597
	for (int iedge = 0; iedge < 4; iedge++)
Thomas Witkowski's avatar
Thomas Witkowski committed
598
	  projection[geoFace + iedge] = elInfoOld->getProjection(geoFace + ce[iedge]);
599
	for (int iedge = 4; iedge < 6; iedge++)
Thomas Witkowski's avatar
Thomas Witkowski committed
600
	  projection[geoFace + iedge] = elInfoOld->getProjection(5 - cv[iedge - 3]);
601
      }
602
    }
603

604
    
Thomas Witkowski's avatar
Thomas Witkowski committed
605
    if (fillFlag.isSet(Mesh::FILL_ORIENTATION)) {
606
607
      orientation = 
	(dynamic_cast<ElInfo3d*>(const_cast<ElInfo*>(elInfoOld)))->orientation 
608
609
610
611
	* Tetrahedron::childOrientation[el_type_local][ichild];
    }
  }

612

613
  mtl::dense2D<double>& ElInfo3d::getSubElemCoordsMat(int degree) const
Thomas Witkowski's avatar
Thomas Witkowski committed
614
  {
615
616
    FUNCNAME("ElInfo3d::getSubElemCoordsMat()");

617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
    using namespace mtl;

    if (subElemMatrices[degree].count(refinementPath) == 0) {
      switch (degree) {
      case 1:
	{
	  dense2D<double> mat(4, 4), tmpMat(4, 4);
	  mat = 1;
	  
	  for (int i = 0; i < refinementPathLength; i++) {
	    if (refinementPath & (1 << i)) {
	      if ((level + i) % 3 == 0)
		tmpMat = mat_d1_l0_right * mat;
	      else
		tmpMat = mat_d1_l12_right * mat;

	      mat = tmpMat;
	    } else  {
	      tmpMat = mat_d1_left * mat;
	      mat = tmpMat;
	    }
	  }

	  subElemMatrices[1][refinementPath] = mat;  
	}
      default:
	ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
      }
    }
646
647

    return subElemMatrices[degree][refinementPath];
Thomas Witkowski's avatar
Thomas Witkowski committed
648
649
  }

650

651
  mtl::dense2D<double>& ElInfo3d::getSubElemGradCoordsMat(int degree) const
Thomas Witkowski's avatar
Thomas Witkowski committed
652
  {
653
    FUNCNAME("ElInfo3d::getSubElemGradCoordsMat()");
654
655

    ERROR_EXIT("Not yet implemented!\n");
656

657
    return subElemGradMatrices[degree][refinementPath];
Thomas Witkowski's avatar
Thomas Witkowski committed
658
659
  }

660
}