ElInfo3d.cc 18 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
121
122
123
    int myRank = omp_get_thread_num();

    WorldVector<double> *e1 = &tmpWorldVecs[myRank][0];
    WorldVector<double> *e2 = &tmpWorldVecs[myRank][1];
    WorldVector<double> *e3 = &tmpWorldVecs[myRank][2];
    WorldVector<double> *v0 = &tmpWorldVecs[myRank][3];

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

    testFlag(Mesh::FILL_COORDS);

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

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

    adet = abs(det);

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

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

    return adet;
  }

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

182
183
184
185
186
187
    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);

188
    TEST_EXIT_DBG(lambda)("lambda must not be NULL\n");
189
190
191

    int dim = mesh_->getDim();

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

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

      for (int i = 0; i < dim; i++)
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
233
234
	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);
235
236
237
238
239

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

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

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

    int neighbours = mesh_->getGeo(NEIGH);
    int vertices = mesh_->getGeo(VERTEX);
  
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
    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;
		  }
294
295
		}
	      }
296
	      TEST_EXIT_DBG(k >= 0)("neighbour dof not found\n");
297
298
299
300
	      
	      if (nb->isNewCoordSet())
		oppCoord_[ineigh] = *(nb->getNewCoord());
	      else
Thomas Witkowski's avatar
Thomas Witkowski committed
301
		for (int j = 0; j < dimOfWorld; j++)
302
		  oppCoord_[ineigh][j] = (oppCoord_[ineigh][j] + coord_[k][j]) / 2;
303
	    }
304
305
	    neighbour_[ineigh] = dynamic_cast<Tetrahedron*>(const_cast<Element*>(nb->getChild(1-ov)));
	    oppVertex_[ineigh] = 3;
306
307
308
	  }
	}
      }
309
    }
310
311
312
  }


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

317
    double det = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
318
319
320
321
322
323

    int myRank = omp_get_thread_num();

    WorldVector<double> *e0 = &tmpWorldVecs[myRank][0];
    WorldVector<double> *e1 = &tmpWorldVecs[myRank][1];
    WorldVector<double> *e2 = &tmpWorldVecs[myRank][2];
324

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
336
      vectorProduct(*e0, *e1, normal);
337

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

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

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

    return(det);
  }




  void ElInfo3d::fillElInfo(int ichild, const ElInfo *elinfo_old)
  {
360
361
362
363
    FUNCNAME("ElInfo3d::fillElInfo()");

    int ochild = 0;             /* index of other child = 1-ichild */
    int *cv = NULL;             /* cv = child_vertex[el_type][ichild] */
364
    const int (*cvg)[4] = NULL;     /* cvg = child_vertex[el_type] */
365
    int *ce;                    /* ce = child_edge[el_type][ichild] */
366
367
368
369
370
    Element *nb, *nbk;
    const FixVec<Element*, NEIGH> *neigh_old;
    Element *el_old = elinfo_old->element_;
    Flag fillFlag__local = elinfo_old->fillFlag_;
    DegreeOfFreedom *dof;
371
    int ov = -1;
372
    FixVec<Element*, NEIGH> *neigh_local;
373
    Flag fill_opp_coords;
374
375
    Mesh *mesh = elinfo_old->getMesh();

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

378
379
    element_ = const_cast<Element*>( el_old->getChild(ichild));
    macroElement_ = elinfo_old->macroElement_;
380
    fillFlag_ = fillFlag__local;
381
382
    parent_  = el_old;
    level_ = elinfo_old->level_ + 1;
383
384
    int el_type_local = ( dynamic_cast<ElInfo3d*>(const_cast<ElInfo*>( elinfo_old)))->getType();
    el_type = (el_type_local + 1) % 3;
385

386
    TEST_EXIT_DBG(element_)("missing child %d?\n", ichild);
387
388
389
390

    if (fillFlag__local.isAnySet()) {
      cvg = Tetrahedron::childVertex[el_type_local];
      cv = const_cast<int*>( cvg[ichild]);
391
      ochild = 1 - ichild;
392
393
    }

394
395
396
397
    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
398
	for (int j = 0; j < dimOfWorld; j++) {
399
	  coord_[i][j] = elinfo_old->coord_[cv[i]][j];
400
401
	}
      }
402
403
404
      if (el_old->getNewCoord()) {
	coord_[3] = *(el_old->getNewCoord());
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
405
	for (int j = 0; j < dimOfWorld; j++) {
406
407
408
409
	  coord_[3][j] = (elinfo_old->coord_[0][j] + elinfo_old->coord_[1][j]) / 2;
	}
      }
    }
410

411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
    if (fillFlag__local.isSet(Mesh::FILL_NEIGH) || 
	fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {

      neigh_local = &neighbour_;
      neigh_old = &elinfo_old->neighbour_;
      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
430
	      for (int j = 0; j < dimOfWorld; j++) {
431
		oppCoord_[0][j] = (elinfo_old->coord_[ochild][j] + elinfo_old->coord_[k][j]) / 2;
432
433
434
	      }
	    }
	  }
435
436
437
438
	  (*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
439
	    for (int j = 0; j < dimOfWorld; j++) {
440
441
442
443
444
	      oppCoord_[0][j] = elinfo_old->coord_[ochild][j];
	    }
	  }
	  (*neigh_local)[0] = nb;
	  oppVertex_[0] = 0;
445
	}
446
447
448
449
      } else {
	ERROR_EXIT("no other child");
	(*neigh_local)[0] = NULL;
      }
450
451


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

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

534
		if (fill_opp_coords.isAnySet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
535
		  for (int j = 0; j < dimOfWorld; j++) {
536
537
538
539
540
541
542
543
544
		    oppCoord_[i][j] = elinfo_old->oppCoord_[cv[i]][j];
		  }
		}
		(*neigh_local)[i] = nbk;
		oppVertex_[i] = ov;
		break;
	      }
	      
	    } /* end for k */
545
	    TEST_EXIT_DBG(k < 2)("child not found with vertex\n");
546
547
548
	  }
	} else {
	  (*neigh_local)[i] = NULL;
549
	}
550
551
552
553
554
555
556
557
558
      }  /* 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
559
	  for (int j = 0; j < dimOfWorld; j++) {
560
561
562
563
564
	    oppCoord_[3][j] = elinfo_old->oppCoord_[ochild][j];
	  }
	}
      }
    }
565

566
567
568
569
570
571
572
573
574
575
576
577
    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);
      
578
579
      int geoFace = mesh_->getGeo(FACE);

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

589
590
591
592
593
594
595
596
597
598
599
      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
600
	  projection_[geoFace + iedge] = elinfo_old->getProjection(geoFace + ce[iedge]);
601
	}
602
	for (int iedge = 4; iedge < 6; iedge++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
603
	  projection_[geoFace + iedge] = elinfo_old->getProjection(5 - cv[iedge - 3]);
604
605
	}
      }
606
    }
607

608
    
609
610
611
612
613
614
615
616
    if (fillFlag_.isSet(Mesh::FILL_ORIENTATION)) {
      orientation =
	( dynamic_cast<ElInfo3d*>(const_cast<ElInfo*>(elinfo_old)))->orientation 
	* Tetrahedron::childOrientation[el_type_local][ichild];
    }
  }

}