ElInfo3d.cc 17.6 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
110
111
	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;
    }
  }

  double ElInfo3d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam) const
  {
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");

117
118
119
120
    WorldVector<double> e1, e2, e3;
    WorldVector<double> v0;
    double det, adet;
    double a11, a12, a13, a21, a22, a23, a31, a32, a33;
121
122
123
124

    testFlag(Mesh::FILL_COORDS);

    v0 = coord_[0];
125
    for (int i = 0; i < 3; i++) {
126
127
128
129
130
      e1[i] = coord_[1][i] - v0[i];
      e2[i] = coord_[2][i] - v0[i];
      e3[i] = coord_[3][i] - v0[i];
    }

131
132
133
    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]);
134
135
136

    adet = abs(det);

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

      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;

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

    return adet;
  }

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

177
178
179
180
181
182
    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);

183
    TEST_EXIT_DBG(lambda)("lambda must not be NULL\n");
184
185
186

    int dim = mesh_->getDim();

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

195
    for (int j = 0; j < dimOfWorld; j++) {
196
197
      x0 = coord_[dim][j];
      x[j] = xy[j] - x0;
198
199

      for (int i = 0; i < dim; i++)
200
201
202
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
	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);
230
231
232
233
234

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

235
236
237
238
239
240
241
242
243
244
      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;
245
246

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

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


  double ElInfo3d::getNormal(int face, WorldVector<double> &normal) const
  {
    FUNCNAME("ElInfo3d::getNormal");
    double det = 0.0;
    WorldVector<double> e0, e1, e2;

Thomas Witkowski's avatar
Thomas Witkowski committed
314
    if (dimOfWorld == 3) {
315
316
317
      int i0 = (face + 1) % 4;
      int i1 = (face + 2) % 4;
      int i2 = (face + 3) % 4;
318

Thomas Witkowski's avatar
Thomas Witkowski committed
319
      for (int i = 0; i < dimOfWorld; i++) {
320
321
322
323
324
325
326
327
	e0[i] = coord_[i1][i] - coord_[i0][i];
	e1[i] = coord_[i2][i] - coord_[i0][i];
	e2[i] = coord_[face][i] - coord_[i0][i];
      }

      vectorProduct(e0, e1, normal);

      if ((e2*normal) < 0.0)
Thomas Witkowski's avatar
Thomas Witkowski committed
328
	for (int i = 0; i < dimOfWorld; i++)
329
330
331
	  normal[i] = -normal[i];

      det = norm(&normal);
332
      TEST_EXIT_DBG(det > 1.e-30)("det = 0 on face %d\n", face);
333
334
335
336
337

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

    return(det);
  }




  void ElInfo3d::fillElInfo(int ichild, const ElInfo *elinfo_old)
  {
349
350
351
352
353
    FUNCNAME("ElInfo3d::fillElInfo()");

    int el_type_local = 0;      /* el_type in {0,1,2} */
    int ochild = 0;             /* index of other child = 1-ichild */
    int *cv = NULL;             /* cv = child_vertex[el_type][ichild] */
354
    const int (*cvg)[4] = NULL;     /* cvg = child_vertex[el_type] */
355
    int *ce;                    /* ce = child_edge[el_type][ichild] */
356
357
358
359
360
    Element *nb, *nbk;
    const FixVec<Element*, NEIGH> *neigh_old;
    Element *el_old = elinfo_old->element_;
    Flag fillFlag__local = elinfo_old->fillFlag_;
    DegreeOfFreedom *dof;
361
    int ov = -1;
362
    FixVec<Element*, NEIGH> *neigh_local;
363
    Flag fill_opp_coords;
364
365
    Mesh *mesh = elinfo_old->getMesh();

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

368
369
    element_ = const_cast<Element*>( el_old->getChild(ichild));
    macroElement_ = elinfo_old->macroElement_;
370
    fillFlag_ = fillFlag__local;
371
372
373
    parent_  = el_old;
    level_ = elinfo_old->level_ + 1;
    el_type = (( dynamic_cast<ElInfo3d*>(const_cast<ElInfo*>( elinfo_old)))->el_type + 1) % 3;
374

375
    TEST_EXIT_DBG(element_)("missing child %d?\n", ichild);
376
377
378
379
380

    if (fillFlag__local.isAnySet()) {
      el_type_local = ( dynamic_cast<ElInfo3d*>(const_cast<ElInfo*>( elinfo_old)))->getType();
      cvg = Tetrahedron::childVertex[el_type_local];
      cv = const_cast<int*>( cvg[ichild]);
381
      ochild = 1 - ichild;
382
383
    }

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

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


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

479
	      if (fill_opp_coords.isAnySet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
480
		for (int j = 0; j < dimOfWorld; j++) {
481
482
		  oppCoord_[i][j] = elinfo_old->oppCoord_[cv[i]][j];
		}
483
	      }
484
485
486
	      (*neigh_local)[i] = nbk;
	      oppVertex_[i] = ov;
	      break;
487
	    }
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
	    
	  } /* 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
508
			for (int j = 0; j < dimOfWorld; j++) {
509
510
511
512
513
514
515
516
517
			  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 {
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
		if (fill_opp_coords.isAnySet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
525
		  for (int j = 0; j < dimOfWorld; j++) {
526
527
528
529
530
531
532
533
534
		    oppCoord_[i][j] = elinfo_old->oppCoord_[cv[i]][j];
		  }
		}
		(*neigh_local)[i] = nbk;
		oppVertex_[i] = ov;
		break;
	      }
	      
	    } /* end for k */
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
546
547
548
      }  /* 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
549
	  for (int j = 0; j < dimOfWorld; j++) {
550
551
552
553
554
	    oppCoord_[3][j] = elinfo_old->oppCoord_[ochild][j];
	  }
	}
      }
    }
555

556
557
558
559
560
561
562
563
564
565
566
567
    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);
      
568
569
      int geoFace = mesh_->getGeo(FACE);

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

579
580
581
582
583
584
585
586
587
588
589
      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
590
	  projection_[geoFace + iedge] = elinfo_old->getProjection(geoFace + ce[iedge]);
591
	}
592
	for (int iedge = 4; iedge < 6; iedge++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
593
	  projection_[geoFace + iedge] = elinfo_old->getProjection(5 - cv[iedge - 3]);
594
595
	}
      }
596
    }
597

598
    
599
600
601
602
603
604
605
606
    if (fillFlag_.isSet(Mesh::FILL_ORIENTATION)) {
      orientation =
	( dynamic_cast<ElInfo3d*>(const_cast<ElInfo*>(elinfo_old)))->orientation 
	* Tetrahedron::childOrientation[el_type_local][ichild];
    }
  }

}