ElInfo3d.cc 17.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#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)
  {
    FUNCNAME("ElInfo3d::fillMacroInfo");
    Element      *nb;
    MacroElement *mnb;
    Flag    fill_opp_coords;

    macroElement_  = const_cast<MacroElement*>( mel);
25
26
27
28
    element_  = const_cast<Element*>( mel->getElement());
    parent_ = NULL;
    level_ = 0;
    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
354
355
    }

    return(det);
  }




  void ElInfo3d::fillElInfo(int ichild, const ElInfo *elinfo_old)
  {
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
363
364
365
    Element *nb, *nbk;
    Element *el_old = elinfo_old->element_;
    Flag fillFlag__local = elinfo_old->fillFlag_;
    DegreeOfFreedom *dof;
366
    int ov = -1;
367
368
    Mesh *mesh = elinfo_old->getMesh();

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

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

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

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

387
388
389
390
    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
391
	for (int j = 0; j < dimOfWorld; j++) {
392
	  coord_[i][j] = elinfo_old->coord_[cv[i]][j];
393
394
	}
      }
395
396
397
      if (el_old->getNewCoord()) {
	coord_[3] = *(el_old->getNewCoord());
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
398
	for (int j = 0; j < dimOfWorld; j++) {
399
400
401
402
	  coord_[3][j] = (elinfo_old->coord_[0][j] + elinfo_old->coord_[1][j]) / 2;
	}
      }
    }
403

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

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

      Flag fill_opp_coords;
411
412
413
414
415
416
417
418
419
420
421
422
423
424
      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
425
	      for (int j = 0; j < dimOfWorld; j++) {
426
		oppCoord_[0][j] = (elinfo_old->coord_[ochild][j] + elinfo_old->coord_[k][j]) / 2;
427
428
429
	      }
	    }
	  }
430
431
432
433
	  (*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
434
	    for (int j = 0; j < dimOfWorld; j++) {
435
436
437
438
439
	      oppCoord_[0][j] = elinfo_old->coord_[ochild][j];
	    }
	  }
	  (*neigh_local)[0] = nb;
	  oppVertex_[0] = 0;
440
	}
441
442
443
444
      } else {
	ERROR_EXIT("no other child");
	(*neigh_local)[0] = NULL;
      }
445
446


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

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

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

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

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

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

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

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

}