ElInfo3d.cc 17.5 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
    FUNCNAME("ElInfo3d::fillElInfo()");

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

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

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

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

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

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

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


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

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

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

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

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

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

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

}