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
    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;  
	}
642
	break;
643
644
645
646
      default:
	ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
      }
    }
647
648

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

651

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

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

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

661
}