ElInfo3d.cc 18.6 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");

134
135
136
    std::vector<double> &e1 = tmpWorldVecs[0];
    std::vector<double> &e2 = tmpWorldVecs[1];
    std::vector<double> &e3 = tmpWorldVecs[2];
137
138
139

    testFlag(Mesh::FILL_COORDS);

140
    for (int i = 0; i < 3; i++) {
141
142
143
      e1[i] = coord[1][i] - coord[0][i];
      e2[i] = coord[2][i] - coord[0][i];
      e3[i] = coord[3][i] - coord[0][i];
144
145
    }

146
147
148
149
    double 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]);
150

151
    double adet = abs(det);
152

153
154
155
156
157
158
159
    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;
160
161
162
163
164
165
166
167
168
169
170
      /* (a_ij) = A^{-T} */

      grd_lam[1][0] = (e2[1] * e3[2] - e2[2] * e3[1]) * det;
      grd_lam[1][1] = (e2[2] * e3[0] - e2[0] * e3[2]) * det;
      grd_lam[1][2] = (e2[0] * e3[1] - e2[1] * e3[0]) * det;
      grd_lam[2][0] = (e1[2] * e3[1] - e1[1] * e3[2]) * det;
      grd_lam[2][1] = (e1[0] * e3[2] - e1[2] * e3[0]) * det;
      grd_lam[2][2] = (e1[1] * e3[0] - e1[0] * e3[1]) * det;
      grd_lam[3][0] = (e1[1] * e2[2] - e1[2] * e2[1]) * det;
      grd_lam[3][1] = (e1[2] * e2[0] - e1[0] * e2[2]) * det;
      grd_lam[3][2] = (e1[0] * e2[1] - e1[1] * e2[0]) * det;
171

Thomas Witkowski's avatar
Thomas Witkowski committed
172
173
174
      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];
175
    }
176
177
178
179

    return adet;
  }

180

181
182
183
  const int ElInfo3d::worldToCoord(const WorldVector<double>& xy,
				   DimVec<double>* lambda) const
  {
184
185
    FUNCNAME("ElInfo::worldToCoord()");

Thomas Witkowski's avatar
Thomas Witkowski committed
186
    DimVec<WorldVector<double> > edge(mesh->getDim(), NO_INIT);
187
188
189
    WorldVector<double> x;
    double  x0, det, det0, det1, det2;
  
Thomas Witkowski's avatar
Thomas Witkowski committed
190
    static DimVec<double> vec(mesh->getDim(), NO_INIT);
191

192
    TEST_EXIT_DBG(lambda)("lambda must not be NULL\n");
193

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

196
    TEST_EXIT_DBG(dim == dimOfWorld)("dim!=dimOfWorld not yet implemented\n");
197
198
199
200
201
202
203
    
    /*  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                 */

204
    for (int j = 0; j < dimOfWorld; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
205
      x0 = coord[dim][j];
206
      x[j] = xy[j] - x0;
207
208

      for (int i = 0; i < dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
209
	edge[i][j] = coord[i][j] - x0;
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
    }

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

240
      for (int i = 0; i <= dim; i++)
241
242
	(*lambda)[i] = 1.0 / dim;

243
244
245
246
247
248
249
250
251
252
      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;
253
254

    for (int i = 0; i <= dim; i++) {
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
      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()
  {
273
    FUNCNAME("ElInfo::update()");
274

Thomas Witkowski's avatar
Thomas Witkowski committed
275
276
    int neighbours = mesh->getGeo(NEIGH);
    int vertices = mesh->getGeo(VERTEX);
277
  
Thomas Witkowski's avatar
Thomas Witkowski committed
278
    if (fillFlag.isSet(Mesh::FILL_NEIGH) || fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {
279
      Tetrahedron *nb;
Thomas Witkowski's avatar
Thomas Witkowski committed
280
      Flag fill_opp_coords = fillFlag & Mesh::FILL_OPP_COORDS;
281
282
      
      for (int ineigh = 0; ineigh < neighbours; ineigh++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
283
	if ((nb = dynamic_cast<Tetrahedron*>(const_cast<Element*>(neighbour[ineigh])))) {
284
	  int ov = oppVertex[ineigh];
285
286
287
288
	  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
289
		if (element->getDOF(j) == nb->getDOF(1 - ov)) 
290
291
292
		  k = j;
	      
	      if (k == -1) {
293
		for (int j = 0; j < vertices; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
294
		  if (mesh->associated(element->getDOF(j, 0), nb->getDOF(1 - ov, 0)))
295
		    k = j;
296
	      }
297
	      TEST_EXIT_DBG(k >= 0)("neighbour dof not found\n");
298
299
	      
	      if (nb->isNewCoordSet())
Thomas Witkowski's avatar
Thomas Witkowski committed
300
		oppCoord[ineigh] = *(nb->getNewCoord());
301
	      else
Thomas Witkowski's avatar
Thomas Witkowski committed
302
		for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
303
		  oppCoord[ineigh][j] = (oppCoord[ineigh][j] + coord[k][j]) / 2;
304
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
305
	    neighbour[ineigh] = dynamic_cast<Tetrahedron*>(const_cast<Element*>(nb->getChild(1-ov)));
306
	    oppVertex[ineigh] = 3;
307
308
309
	  }
	}
      }
310
    }
311
312
313
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
314
  double ElInfo3d::getNormal(int face, WorldVector<double> &normal)
315
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
316
317
    FUNCNAME("ElInfo3d::getNormal()");

318
    double det = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
319

320
    WorldVector<double> e0, e1, e2;
321

Thomas Witkowski's avatar
Thomas Witkowski committed
322
    if (dimOfWorld == 3) {
323
324
325
      int i0 = (face + 1) % 4;
      int i1 = (face + 2) % 4;
      int i2 = (face + 3) % 4;
326

Thomas Witkowski's avatar
Thomas Witkowski committed
327
      for (int i = 0; i < dimOfWorld; i++) {
328
329
330
	e0[i] = coord[i1][i] - coord[i0][i];
	e1[i] = coord[i2][i] - coord[i0][i];
	e2[i] = coord[face][i] - coord[i0][i];
331
332
      }

333
      vectorProduct(e0, e1, normal);
334

335
      if ((e2 * normal) < 0.0)
Thomas Witkowski's avatar
Thomas Witkowski committed
336
	for (int i = 0; i < dimOfWorld; i++)
337
338
339
	  normal[i] = -normal[i];

      det = norm(&normal);
340
      TEST_EXIT_DBG(det > 1.e-30)("det = 0 on face %d\n", face);
341
342
343
344
345

      normal[0] /= det;
      normal[1] /= det;
      normal[2] /= det;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
346
      MSG("not implemented for DIM_OF_WORLD = %d in 3d\n", dimOfWorld);
347
348
    }

349
    return det;
350
351
352
  }


353
  void ElInfo3d::fillElInfo(int ichild, const ElInfo *elInfoOld)
354
  {
355
356
357
358
    FUNCNAME("ElInfo3d::fillElInfo()");

    int ochild = 0;             /* index of other child = 1-ichild */
    int *cv = NULL;             /* cv = child_vertex[el_type][ichild] */
359
    const int (*cvg)[4] = NULL;     /* cvg = child_vertex[el_type] */
360
    int *ce;                    /* ce = child_edge[el_type][ichild] */
361
    Element *nb, *nbk;
Thomas Witkowski's avatar
Thomas Witkowski committed
362
363
    Element *el_old = elInfoOld->element;
    Flag fillFlag_local = elInfoOld->fillFlag;
364
    DegreeOfFreedom *dof;
365
    int ov = -1;
366
    Mesh *mesh = elInfoOld->getMesh();
367

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

Thomas Witkowski's avatar
Thomas Witkowski committed
370
371
372
373
    element = const_cast<Element*>( el_old->getChild(ichild));
    macroElement = elInfoOld->macroElement;
    fillFlag = fillFlag_local;
    parent  = el_old;
374
    level = elInfoOld->level + 1;
375
376
    iChild = ichild;
    int el_type_local = (dynamic_cast<const ElInfo3d*>(elInfoOld))->getType();
377
    elType = (el_type_local + 1) % 3;
378

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

Thomas Witkowski's avatar
Thomas Witkowski committed
381
    if (fillFlag_local.isAnySet()) {
382
      cvg = Tetrahedron::childVertex[el_type_local];
383
      cv = const_cast<int*>(cvg[ichild]);
384
      ochild = 1 - ichild;
385
386
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
387
388
389
    if (fillFlag_local.isSet(Mesh::FILL_COORDS) || 
	fillFlag.isSet(Mesh::FILL_DET) ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
390
      for (int i = 0; i < 3; i++) {
391
	for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
392
	  coord[i][j] = elInfoOld->coord[cv[i]][j];
393
      }
394
      if (el_old->getNewCoord()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
395
	coord[3] = *(el_old->getNewCoord());
396
      } else {
397
	for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
398
	  coord[3][j] = (elInfoOld->coord[0][j] + elInfoOld->coord[1][j]) / 2;
399
400
      }
    }
401

Thomas Witkowski's avatar
Thomas Witkowski committed
402
403
    if (fillFlag_local.isSet(Mesh::FILL_NEIGH) || 
	fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {
404

Thomas Witkowski's avatar
Thomas Witkowski committed
405
406
      FixVec<Element*, NEIGH> *neigh_local = &neighbour;
      const FixVec<Element*, NEIGH> *neigh_old = &elInfoOld->neighbour;
407
408

      Flag fill_opp_coords;
Thomas Witkowski's avatar
Thomas Witkowski committed
409
      fill_opp_coords.setFlags(fillFlag_local & Mesh::FILL_OPP_COORDS);
410
411
412
413
414
415
416
417
418
419
      
      /*----- 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
420
	      oppCoord[0]= *(nb->getNewCoord());
421
422
	    } else {
	      int k = cvg[ochild][1];
423
	      for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
424
		oppCoord[0][j] = (elInfoOld->coord[ochild][j] + elInfoOld->coord[k][j]) / 2;
425
426
	    }
	  }
427
	  (*neigh_local)[0] = const_cast<Element*>( nb->getChild(1));
428
	  oppVertex[0] = 3;
429
	} else {
430
431
	  if (fill_opp_coords.isAnySet())
	    for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
432
	      oppCoord[0][j] = elInfoOld->coord[ochild][j];
433

434
	  (*neigh_local)[0] = nb;
435
	  oppVertex[0] = 0;
436
	}
437
438
439
440
      } else {
	ERROR_EXIT("no other child");
	(*neigh_local)[0] = NULL;
      }
441
442


443
444
445
446
      /*----- 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]]))) {
447
	  TEST_EXIT_DBG(nb->getChild(0))("nonconforming triangulation\n");
448
	  
449
450
451
452
453
	  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 */
454
	      dof = const_cast<int*>(nb->getDOF(elInfoOld->oppVertex[cv[i]])); 
455
456
457
458
459
460
	      
	      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
461
		      oppCoord[i] = *(nbk->getNewCoord());
462
		    } else {
463
		      for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
464
465
			oppCoord[i][j] = (elInfoOld->oppCoord[cv[i]][j] + 
					   elInfoOld->coord[ichild][j]) / 2;		      
466
467
468
		    }
		  }
		  (*neigh_local)[i] = nbk->getChild(0);
469
		  oppVertex[i] = 3;
470
		  break;
471
		}
472
473
474
475
476
477
	      } else {
		if (dof != nbk->getDOF(2)) { 
		  ov = -1; 
		  break; 
		}
		ov = 2;
478
479
	      }

480
481
	      if (fill_opp_coords.isAnySet())
		for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
482
		  oppCoord[i][j] = elInfoOld->oppCoord[cv[i]][j];
483

484
	      (*neigh_local)[i] = nbk;
485
	      oppVertex[i] = ov;
486
	      break;
487
	    }
488
489
490
491
492
493
494
495
496
497
	    
	  } /* 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 */
498
		dof = const_cast<int*>(nb->getDOF(elInfoOld->oppVertex[cv[i]])); 
499
500
501
502
503
504
505
		
		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
506
			oppCoord[i] = *(nbk->getNewCoord());
507
		      } else {
508
			for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
509
510
			  oppCoord[i][j] = (elInfoOld->oppCoord[cv[i]][j] + 
					    elInfoOld->coord[ichild][j]) / 2;
511
512
513
		      }
		    }
		    (*neigh_local)[i] = nbk->getChild(0);
514
		    oppVertex[i] = 3;
515
516
517
		    break;
		  }
		} else {
518
519
		  TEST_EXIT_DBG(dof == nbk->getDOF(2) || 
				mesh->associated(dof[0], nbk->getDOF(2, 0)))
520
521
522
		    ("opp_vertex not found\n");
		  ov = 2;
		}
523

524
525
		if (fill_opp_coords.isAnySet())
		  for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
526
		    oppCoord[i][j] = elInfoOld->oppCoord[cv[i]][j];
527

528
		(*neigh_local)[i] = nbk;
529
		oppVertex[i] = ov;
530
531
532
533
		break;
	      }
	      
	    } /* end for k */
534

535
	    TEST_EXIT_DBG(k < 2)("child not found with vertex\n");
536
537
538
	  }
	} else {
	  (*neigh_local)[i] = NULL;
539
	}
540
541
542
543
544
545
      }  /* end for i */
      
      
      /*----- nb[3] is old neighbour neigh_old[ochild] ------------------------*/
      
      if (((*neigh_local)[3] = (*neigh_old)[ochild])) {
546
	oppVertex[3] = elInfoOld->oppVertex[ochild];
547

548
549
	if (fill_opp_coords.isAnySet())
	  for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
550
	    oppCoord[3][j] = elInfoOld->oppCoord[ochild][j];
551
552
      }
    }
553

Thomas Witkowski's avatar
Thomas Witkowski committed
554
    if (fillFlag_local.isSet(Mesh::FILL_BOUND)) {
555
      for (int i = 0; i < 3; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
556
	boundary[10 + i] = elInfoOld->getBoundary(10 + cv[i]);
557
      
Thomas Witkowski's avatar
Thomas Witkowski committed
558
      boundary[13] = elInfoOld->getBoundary(4);
559
      
Thomas Witkowski's avatar
Thomas Witkowski committed
560
561
562
563
      boundary[0] = INTERIOR;
      boundary[1] = elInfoOld->getBoundary(cv[1]);
      boundary[2] = elInfoOld->getBoundary(cv[2]);
      boundary[3] = elInfoOld->getBoundary(ochild);
564
      
Thomas Witkowski's avatar
Thomas Witkowski committed
565
      int geoFace = mesh->getGeo(FACE);
566

567
      ce = const_cast<int*>(Tetrahedron::childEdge[el_type_local][ichild]);
568
      for (int iedge = 0; iedge < 4; iedge++)
Thomas Witkowski's avatar
Thomas Witkowski committed
569
570
571
	boundary[geoFace + iedge] = elInfoOld->getBoundary(geoFace + ce[iedge]);      
      for (int iedge = 4; iedge < 6; iedge++)
	boundary[geoFace + iedge] = elInfoOld->getBoundary(5 - cv[iedge - 3]);
572

573
574
      if (elInfoOld->getProjection(0) &&
	  elInfoOld->getProjection(0)->getType() == VOLUME_PROJECTION) {
575
	
Thomas Witkowski's avatar
Thomas Witkowski committed
576
	projection[0] = elInfoOld->getProjection(0);      
577
      } else { // boundary projection
Thomas Witkowski's avatar
Thomas Witkowski committed
578
579
580
581
	projection[0] = NULL;
	projection[1] = elInfoOld->getProjection(cv[1]);
	projection[2] = elInfoOld->getProjection(cv[2]);
	projection[3] = elInfoOld->getProjection(ochild);
582
	
583
	for (int iedge = 0; iedge < 4; iedge++)
Thomas Witkowski's avatar
Thomas Witkowski committed
584
	  projection[geoFace + iedge] = elInfoOld->getProjection(geoFace + ce[iedge]);
585
	for (int iedge = 4; iedge < 6; iedge++)
Thomas Witkowski's avatar
Thomas Witkowski committed
586
	  projection[geoFace + iedge] = elInfoOld->getProjection(5 - cv[iedge - 3]);
587
      }
588
    }
589

590
    
Thomas Witkowski's avatar
Thomas Witkowski committed
591
    if (fillFlag.isSet(Mesh::FILL_ORIENTATION)) {
592
593
      orientation = 
	(dynamic_cast<ElInfo3d*>(const_cast<ElInfo*>(elInfoOld)))->orientation 
594
595
596
597
	* Tetrahedron::childOrientation[el_type_local][ichild];
    }
  }

598

599
  mtl::dense2D<double>& ElInfo3d::getSubElemCoordsMat(int degree) const
Thomas Witkowski's avatar
Thomas Witkowski committed
600
  {
601
602
    FUNCNAME("ElInfo3d::getSubElemCoordsMat()");

603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
    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;  
	}
628
	break;
629
630
631
632
      default:
	ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
      }
    }
633
634

    return subElemMatrices[degree][refinementPath];
Thomas Witkowski's avatar
Thomas Witkowski committed
635
636
  }

637

638
  mtl::dense2D<double>& ElInfo3d::getSubElemGradCoordsMat(int degree) const
Thomas Witkowski's avatar
Thomas Witkowski committed
639
  {
640
    FUNCNAME("ElInfo3d::getSubElemGradCoordsMat()");
641
642

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

644
    return subElemGradMatrices[degree][refinementPath];
Thomas Witkowski's avatar
Thomas Witkowski committed
645
646
  }

647
}