ElInfo3d.cc 19.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#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 {

  void ElInfo3d::fillMacroInfo(const MacroElement * mel)
  {
19
20
    FUNCNAME("ElInfo3d::fillMacroInfo()");
    Element *nb;
21
    MacroElement *mnb;
22
    Flag fill_opp_coords;
23
24

    macroElement_  = const_cast<MacroElement*>( mel);
25
26
    element_  = const_cast<Element*>( mel->getElement());
    parent_ = NULL;
27
    level = 0;
28
    el_type = const_cast<MacroElement*>(mel)->getElType();
29
30
31
32
33
34

    int vertices = mesh_->getGeo(VERTEX);

    if (fillFlag_.isSet(Mesh::FILL_COORDS) || 
	fillFlag_.isSet(Mesh::FILL_DET) ||
	fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
35
      for (int i = 0; i < vertices; i++) {
36
37
38
39
40
41
	coord_[i] = mel->coord[i];
      }
    }

    int neighbours = mesh_->getGeo(NEIGH);

42
43
    if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS) || 
        fillFlag_.isSet(Mesh::FILL_NEIGH)) {
44
45

      fill_opp_coords.setFlags(fillFlag_ & Mesh::FILL_OPP_COORDS);
46
      for (int i = 0; i < neighbours; i++) {
47
48
49
	if ((mnb = const_cast<MacroElement*>( mel->getNeighbour(i)))) {
	  neighbour_[i] = const_cast<Element*>( mel->getNeighbour(i)->getElement());
	  nb = const_cast<Element*>( neighbour_[i]);
50
51
	  int k;
	  k = oppVertex_[i] = mel->getOppVertex(i);
52
53

	  if (nb->getChild(0) && (k < 2)) {   /*make nb nearest element.*/
54
	    if (k == 1) {
55
56
57
58
59
60
	      neighbour_[i]      = const_cast<Element*>( nb->getChild(0));
	      nb = const_cast<Element*>( neighbour_[i]);
	    } else {
	      neighbour_[i]      = const_cast<Element*>( nb->getChild(1));
	      nb = const_cast<Element*>( neighbour_[i]);
	    }
61
62
63
64
65
66
67
68
69
	    k = oppVertex_[i] = 3;
	    if (fill_opp_coords.isAnySet()) {
	      /* always edge between vertices 0 and 1 is bisected! */
	      if (mnb->getElement()->isNewCoordSet())
		oppCoord_[i] = *(mnb->getElement()->getNewCoord());
	      else
		oppCoord_[i] = (mnb->coord[0] + mnb->coord[1]) * 0.5;
	    }
	  } else {
70
71
72
73
	    if  (fill_opp_coords.isAnySet()) {
	      oppCoord_[i] = mnb->coord[k];
	    }
	  }
74
	} else {
75
76
77
78
79
	  neighbour_[i] = NULL;
	}
      }
    }

80
81
82
83
84
85
86
    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);
87
      }
88
    }
89
90
91
92
93

    if (fillFlag_.isSet(Mesh::FILL_ORIENTATION)) {
      WorldVector<WorldVector<double> > a;
      double s;

94
95
      for (int i = 0; i < 3; i++) {
	a[i] = mel->coord[i + 1];
96
97
98
99
100
101
102
103
104
105
106
107
108
109
	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;
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
110
  double ElInfo3d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
111
  {
112
113
    FUNCNAME("ElInfo3d::calcGrdLambda()");

Thomas Witkowski's avatar
Thomas Witkowski committed
114
    TEST_EXIT_DBG(dimOfWorld == 3)
115
116
      ("dim != dim_of_world ! use parametric elements!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
117
118
119
120
    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
121

122
123
    double det, adet;
    double a11, a12, a13, a21, a22, a23, a31, a32, a33;
124
125
126

    testFlag(Mesh::FILL_COORDS);

Thomas Witkowski's avatar
Thomas Witkowski committed
127
    *v0 = coord_[0];
128
    for (int i = 0; i < 3; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
129
130
131
      (*e1)[i] = coord_[1][i] - (*v0)[i];
      (*e2)[i] = coord_[2][i] - (*v0)[i];
      (*e3)[i] = coord_[3][i] - (*v0)[i];
132
133
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
134
135
136
    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]);
137
138
139

    adet = abs(det);

140
141
142
143
144
145
146
    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
147
148
149
150
151
152
153
154
155
      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;
156
157
158
159
160
161
162
163
164
165
166

      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
167
168
169
      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];
170
    }
171
172
173
174
175
176
177

    return adet;
  }

  const int ElInfo3d::worldToCoord(const WorldVector<double>& xy,
				   DimVec<double>* lambda) const
  {
178
179
    FUNCNAME("ElInfo::worldToCoord()");

180
181
182
183
184
185
    DimVec<WorldVector<double> > edge(mesh_->getDim(), NO_INIT);
    WorldVector<double> x;
    double  x0, det, det0, det1, det2;
  
    static DimVec<double> vec(mesh_->getDim(), NO_INIT);

186
    TEST_EXIT_DBG(lambda)("lambda must not be NULL\n");
187
188
189

    int dim = mesh_->getDim();

190
    TEST_EXIT_DBG(dim == dimOfWorld)("dim!=dimOfWorld not yet implemented\n");
191
192
193
194
195
196
197
    
    /*  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                 */

198
    for (int j = 0; j < dimOfWorld; j++) {
199
200
      x0 = coord_[dim][j];
      x[j] = xy[j] - x0;
201
202

      for (int i = 0; i < dim; i++)
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
	edge[i][j] = coord_[i][j] - x0;
    }

    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);
233
234
235
236
237

      for (int i = 0; i <= dim; i++) {
	(*lambda)[i] = 1.0 / dim;
      }

238
239
240
241
242
243
244
245
246
247
      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;
248
249

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

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


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

315
    double det = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
316

Thomas Witkowski's avatar
Thomas Witkowski committed
317
318
319
    WorldVector<double> *e0 = &tmpWorldVecs[0];
    WorldVector<double> *e1 = &tmpWorldVecs[1];
    WorldVector<double> *e2 = &tmpWorldVecs[2];
320

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
332
      vectorProduct(*e0, *e1, normal);
333

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

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

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

    return(det);
  }




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

    int ochild = 0;             /* index of other child = 1-ichild */
    int *cv = NULL;             /* cv = child_vertex[el_type][ichild] */
360
    const int (*cvg)[4] = NULL;     /* cvg = child_vertex[el_type] */
361
    int *ce;                    /* ce = child_edge[el_type][ichild] */
362
    Element *nb, *nbk;
363
364
    Element *el_old = elInfoOld->element_;
    Flag fillFlag__local = elInfoOld->fillFlag_;
365
    DegreeOfFreedom *dof;
366
    int ov = -1;
367
    Mesh *mesh = elInfoOld->getMesh();
368

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

371
    element_ = const_cast<Element*>( el_old->getChild(ichild));
372
    macroElement_ = elInfoOld->macroElement_;
373
    fillFlag_ = fillFlag__local;
374
    parent_  = el_old;
375
    level = elInfoOld->level + 1;
376
377
    iChild = ichild;
    int el_type_local = (dynamic_cast<const ElInfo3d*>(elInfoOld))->getType();
378
    el_type = (el_type_local + 1) % 3;
379

380
    TEST_EXIT_DBG(element_)("missing child %d?\n", ichild);
381
382
383

    if (fillFlag__local.isAnySet()) {
      cvg = Tetrahedron::childVertex[el_type_local];
384
      cv = const_cast<int*>(cvg[ichild]);
385
      ochild = 1 - ichild;
386
387
    }

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

405
406
407
    if (fillFlag__local.isSet(Mesh::FILL_NEIGH) || 
	fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {

408
      FixVec<Element*, NEIGH> *neigh_local = &neighbour_;
409
      const FixVec<Element*, NEIGH> *neigh_old = &elInfoOld->neighbour_;
410
411

      Flag fill_opp_coords;
412
413
414
415
416
417
418
419
420
421
422
423
424
425
      fill_opp_coords.setFlags(fillFlag__local & Mesh::FILL_OPP_COORDS);
      
      /*----- 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()) {
	      oppCoord_[0]= *(nb->getNewCoord());
	    } else {
	      int k = cvg[ochild][1];
Thomas Witkowski's avatar
Thomas Witkowski committed
426
	      for (int j = 0; j < dimOfWorld; j++) {
427
		oppCoord_[0][j] = (elInfoOld->coord_[ochild][j] + elInfoOld->coord_[k][j]) / 2;
428
429
430
	      }
	    }
	  }
431
432
433
434
	  (*neigh_local)[0] = const_cast<Element*>( nb->getChild(1));
	  oppVertex_[0] = 3;
	} else {
	  if (fill_opp_coords.isAnySet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
435
	    for (int j = 0; j < dimOfWorld; j++) {
436
	      oppCoord_[0][j] = elInfoOld->coord_[ochild][j];
437
438
439
440
	    }
	  }
	  (*neigh_local)[0] = nb;
	  oppVertex_[0] = 0;
441
	}
442
443
444
445
      } else {
	ERROR_EXIT("no other child");
	(*neigh_local)[0] = NULL;
      }
446
447


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

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

530
		if (fill_opp_coords.isAnySet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
531
		  for (int j = 0; j < dimOfWorld; j++) {
532
		    oppCoord_[i][j] = elInfoOld->oppCoord_[cv[i]][j];
533
534
535
536
537
538
539
540
		  }
		}
		(*neigh_local)[i] = nbk;
		oppVertex_[i] = ov;
		break;
	      }
	      
	    } /* end for k */
541
	    TEST_EXIT_DBG(k < 2)("child not found with vertex\n");
542
543
544
	  }
	} else {
	  (*neigh_local)[i] = NULL;
545
	}
546
547
548
549
550
551
      }  /* end for i */
      
      
      /*----- nb[3] is old neighbour neigh_old[ochild] ------------------------*/
      
      if (((*neigh_local)[3] = (*neigh_old)[ochild])) {
552
	oppVertex_[3] = elInfoOld->oppVertex_[ochild];
553
554

	if (fill_opp_coords.isAnySet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
555
	  for (int j = 0; j < dimOfWorld; j++) {
556
	    oppCoord_[3][j] = elInfoOld->oppCoord_[ochild][j];
557
558
559
560
	  }
	}
      }
    }
561

562
563
    if (fillFlag__local.isSet(Mesh::FILL_BOUND)) {
      for (int i = 0; i < 3; i++) {
564
	boundary_[10 + i] = elInfoOld->getBoundary(10 + cv[i]);
565
566
      }
      
567
      boundary_[13] = elInfoOld->getBoundary(4);
568
569
      
      boundary_[0] = INTERIOR;
570
571
572
      boundary_[1] = elInfoOld->getBoundary(cv[1]);
      boundary_[2] = elInfoOld->getBoundary(cv[2]);
      boundary_[3] = elInfoOld->getBoundary(ochild);
573
      
574
575
      int geoFace = mesh_->getGeo(FACE);

576
577
      ce = const_cast<int*>(Tetrahedron::childEdge[el_type_local][ichild]);
      for (int iedge = 0; iedge < 4; iedge++) {
578
	boundary_[geoFace + iedge] = elInfoOld->getBoundary(geoFace + ce[iedge]);
579
580
581
      }
      for (int iedge = 4; iedge < 6; iedge++) {
	int i = 5 - cv[iedge - 3];                /* old vertex opposite new edge */
582
	boundary_[geoFace + iedge] = elInfoOld->getBoundary(i);
583
      }
584

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
612
613
614
615
616
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
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
  void ElInfo3d::getRefSimplexCoords(DimMat<double> *coords) const
  {
    (*coords)[0][0] = 1.0;
    (*coords)[1][0] = 0.0;
    (*coords)[2][0] = 0.0;
    (*coords)[3][0] = 0.0;

    (*coords)[0][1] = 0.0;
    (*coords)[1][1] = 1.0;
    (*coords)[2][1] = 0.0;
    (*coords)[3][1] = 0.0;

    (*coords)[0][2] = 0.0;
    (*coords)[1][2] = 0.0;
    (*coords)[2][2] = 1.0;
    (*coords)[3][2] = 0.0;

    (*coords)[0][3] = 0.0;
    (*coords)[1][3] = 0.0;
    (*coords)[2][3] = 0.0;
    (*coords)[3][3] = 1.0;
  }

  void ElInfo3d::getSubElementCoords(DimMat<double> *coords,
				     int iChild) const
  {
    FUNCNAME("ElInfo3d::getSubElementCoords()");

    ERROR_EXIT("Not yet\n");

    double c0 = ((*coords)[0][0] + (*coords)[0][1]) * 0.5;
    double c1 = ((*coords)[1][0] + (*coords)[1][1]) * 0.5;
    double c2 = ((*coords)[2][0] + (*coords)[2][1]) * 0.5;
    double c3 = ((*coords)[3][0] + (*coords)[3][1]) * 0.5;

    if (iChild == 0) {
      (*coords)[0][1] = (*coords)[0][0];
      (*coords)[1][1] = (*coords)[1][0];
      (*coords)[2][1] = (*coords)[2][0];

      (*coords)[0][0] = (*coords)[0][2];
      (*coords)[1][0] = (*coords)[1][2];
      (*coords)[2][0] = (*coords)[2][2];
    } else {
      (*coords)[0][1] = (*coords)[0][2];
      (*coords)[1][1] = (*coords)[1][2];
      (*coords)[2][1] = (*coords)[2][2];

      (*coords)[0][0] = (*coords)[0][1];
      (*coords)[1][0] = (*coords)[1][1];
      (*coords)[2][0] = (*coords)[2][1];      
    }

    (*coords)[0][3] = c0;
    (*coords)[1][3] = c1;
    (*coords)[2][3] = c2;
    (*coords)[3][3] = c3;
  }

671
}