ElInfo3d.cc 19.4 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
	coord_[i] = mel->coord[i];
    }

    int neighbours = mesh_->getGeo(NEIGH);

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

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

	  if (nb->getChild(0) && (k < 2)) {   /*make nb nearest element.*/
53
	    if (k == 1) {
54
55
56
57
58
59
	      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]);
	    }
60
	    k = oppVertex[i] = 3;
61
62
63
64
65
66
67
68
	    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 {
69
70
71
72
	    if  (fill_opp_coords.isAnySet()) {
	      oppCoord_[i] = mnb->coord[k];
	    }
	  }
73
	} else {
74
75
76
77
78
	  neighbour_[i] = NULL;
	}
      }
    }

79
    if (fillFlag_.isSet(Mesh::FILL_BOUND)) {
80
81
82
83
      for (int i = 0; i < element_->getGeo(BOUNDARY); i++)
	boundary_[i] = mel->getBoundary(i);     

      for (int i = 0; i < element_->getGeo(PROJECTION); i++)
84
85
	projection_[i] = mel->getProjection(i);
    }
86
87
88
89
90

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

91
92
      for (int i = 0; i < 3; i++) {
	a[i] = mel->coord[i + 1];
93
94
95
96
97
98
99
100
101
102
103
104
105
106
	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
107
  double ElInfo3d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
108
  {
109
110
    FUNCNAME("ElInfo3d::calcGrdLambda()");

Thomas Witkowski's avatar
Thomas Witkowski committed
111
    TEST_EXIT_DBG(dimOfWorld == 3)
112
113
      ("dim != dim_of_world ! use parametric elements!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
114
115
116
117
    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
118

119
120
    double det, adet;
    double a11, a12, a13, a21, a22, a23, a31, a32, a33;
121
122
123

    testFlag(Mesh::FILL_COORDS);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
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
    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
144
145
146
147
148
149
150
151
152
      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;
153
154
155
156
157
158
159
160
161
162
163

      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
164
165
166
      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];
167
    }
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
      for (int i = 0; i <= dim; i++)
232
233
	(*lambda)[i] = 1.0 / dim;

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

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

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


Thomas Witkowski's avatar
Thomas Witkowski committed
305
  double ElInfo3d::getNormal(int face, WorldVector<double> &normal)
306
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
307
308
    FUNCNAME("ElInfo3d::getNormal()");

309
    double det = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
310

Thomas Witkowski's avatar
Thomas Witkowski committed
311
312
313
    WorldVector<double> *e0 = &tmpWorldVecs[0];
    WorldVector<double> *e1 = &tmpWorldVecs[1];
    WorldVector<double> *e2 = &tmpWorldVecs[2];
314

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
326
      vectorProduct(*e0, *e1, normal);
327

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

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

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

    return(det);
  }




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

    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
    Element *nb, *nbk;
357
358
    Element *el_old = elInfoOld->element_;
    Flag fillFlag__local = elInfoOld->fillFlag_;
359
    DegreeOfFreedom *dof;
360
    int ov = -1;
361
    Mesh *mesh = elInfoOld->getMesh();
362

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

365
    element_ = const_cast<Element*>( el_old->getChild(ichild));
366
    macroElement_ = elInfoOld->macroElement_;
367
    fillFlag_ = fillFlag__local;
368
    parent_  = el_old;
369
    level = elInfoOld->level + 1;
370
371
    iChild = ichild;
    int el_type_local = (dynamic_cast<const ElInfo3d*>(elInfoOld))->getType();
372
    el_type = (el_type_local + 1) % 3;
373

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

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

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

397
398
399
    if (fillFlag__local.isSet(Mesh::FILL_NEIGH) || 
	fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {

400
      FixVec<Element*, NEIGH> *neigh_local = &neighbour_;
401
      const FixVec<Element*, NEIGH> *neigh_old = &elInfoOld->neighbour_;
402
403

      Flag fill_opp_coords;
404
405
406
407
408
409
410
411
412
413
414
415
416
417
      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];
418
	      for (int j = 0; j < dimOfWorld; j++)
419
		oppCoord_[0][j] = (elInfoOld->coord_[ochild][j] + elInfoOld->coord_[k][j]) / 2;
420
421
	    }
	  }
422
	  (*neigh_local)[0] = const_cast<Element*>( nb->getChild(1));
423
	  oppVertex[0] = 3;
424
	} else {
425
426
	  if (fill_opp_coords.isAnySet())
	    for (int j = 0; j < dimOfWorld; j++)
427
	      oppCoord_[0][j] = elInfoOld->coord_[ochild][j];
428

429
	  (*neigh_local)[0] = nb;
430
	  oppVertex[0] = 0;
431
	}
432
433
434
435
      } else {
	ERROR_EXIT("no other child");
	(*neigh_local)[0] = NULL;
      }
436
437


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

475
476
	      if (fill_opp_coords.isAnySet())
		for (int j = 0; j < dimOfWorld; j++)
477
		  oppCoord_[i][j] = elInfoOld->oppCoord_[cv[i]][j];
478

479
	      (*neigh_local)[i] = nbk;
480
	      oppVertex[i] = ov;
481
	      break;
482
	    }
483
484
485
486
487
488
489
490
491
492
	    
	  } /* 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 */
493
		dof = const_cast<int*>(nb->getDOF(elInfoOld->oppVertex[cv[i]])); 
494
495
496
497
498
499
500
501
502
		
		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 {
503
504
505
			for (int j = 0; j < dimOfWorld; j++)
			  oppCoord_[i][j] = (elInfoOld->oppCoord_[cv[i]][j] + 
					     elInfoOld->coord_[ichild][j]) / 2;
506
507
508
		      }
		    }
		    (*neigh_local)[i] = nbk->getChild(0);
509
		    oppVertex[i] = 3;
510
511
512
		    break;
		  }
		} else {
513
514
		  TEST_EXIT_DBG(dof == nbk->getDOF(2) || 
				mesh->associated(dof[0], nbk->getDOF(2, 0)))
515
516
517
		    ("opp_vertex not found\n");
		  ov = 2;
		}
518

519
520
		if (fill_opp_coords.isAnySet())
		  for (int j = 0; j < dimOfWorld; j++)
521
		    oppCoord_[i][j] = elInfoOld->oppCoord_[cv[i]][j];
522

523
		(*neigh_local)[i] = nbk;
524
		oppVertex[i] = ov;
525
526
527
528
		break;
	      }
	      
	    } /* end for k */
529

530
	    TEST_EXIT_DBG(k < 2)("child not found with vertex\n");
531
532
533
	  }
	} else {
	  (*neigh_local)[i] = NULL;
534
	}
535
536
537
538
539
540
      }  /* end for i */
      
      
      /*----- nb[3] is old neighbour neigh_old[ochild] ------------------------*/
      
      if (((*neigh_local)[3] = (*neigh_old)[ochild])) {
541
	oppVertex[3] = elInfoOld->oppVertex[ochild];
542

543
544
	if (fill_opp_coords.isAnySet())
	  for (int j = 0; j < dimOfWorld; j++)
545
	    oppCoord_[3][j] = elInfoOld->oppCoord_[ochild][j];
546
547
      }
    }
548

549
    if (fillFlag__local.isSet(Mesh::FILL_BOUND)) {
550
      for (int i = 0; i < 3; i++)
551
	boundary_[10 + i] = elInfoOld->getBoundary(10 + cv[i]);
552
      
553
      boundary_[13] = elInfoOld->getBoundary(4);
554
555
      
      boundary_[0] = INTERIOR;
556
557
558
      boundary_[1] = elInfoOld->getBoundary(cv[1]);
      boundary_[2] = elInfoOld->getBoundary(cv[2]);
      boundary_[3] = elInfoOld->getBoundary(ochild);
559
      
560
561
      int geoFace = mesh_->getGeo(FACE);

562
      ce = const_cast<int*>(Tetrahedron::childEdge[el_type_local][ichild]);
563
564
      for (int iedge = 0; iedge < 4; iedge++)
	boundary_[geoFace + iedge] = elInfoOld->getBoundary(geoFace + ce[iedge]);      
565
566
      for (int iedge = 4; iedge < 6; iedge++) {
	int i = 5 - cv[iedge - 3];                /* old vertex opposite new edge */
567
	boundary_[geoFace + iedge] = elInfoOld->getBoundary(i);
568
      }
569

570
571
      if (elInfoOld->getProjection(0) &&
	  elInfoOld->getProjection(0)->getType() == VOLUME_PROJECTION) {
572
	
573
	projection_[0] = elInfoOld->getProjection(0);      
574
575
      } else { // boundary projection
	projection_[0] = NULL;
576
577
578
	projection_[1] = elInfoOld->getProjection(cv[1]);
	projection_[2] = elInfoOld->getProjection(cv[2]);
	projection_[3] = elInfoOld->getProjection(ochild);
579
	
580
	for (int iedge = 0; iedge < 4; iedge++)
581
	  projection_[geoFace + iedge] = elInfoOld->getProjection(geoFace + ce[iedge]);
582
	for (int iedge = 4; iedge < 6; iedge++)
583
	  projection_[geoFace + iedge] = elInfoOld->getProjection(5 - cv[iedge - 3]);
584
      }
585
    }
586

587
    
588
    if (fillFlag_.isSet(Mesh::FILL_ORIENTATION)) {
589
590
      orientation = 
	(dynamic_cast<ElInfo3d*>(const_cast<ElInfo*>(elInfoOld)))->orientation 
591
592
593
594
	* Tetrahedron::childOrientation[el_type_local][ichild];
    }
  }

595
#if 0
596
  void ElInfo3d::getRefSimplexCoords(const BasisFunction *basisFcts,
597
				     mtl::dense2D<double>& coords) const
Thomas Witkowski's avatar
Thomas Witkowski committed
598
  {
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
//     (*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;
Thomas Witkowski's avatar
Thomas Witkowski committed
618
619
  }

620
621
  void ElInfo3d::getSubElementCoords(const BasisFunction *basisFcts,
				     int iChild,
622
				     mtl::dense2D<double>& coords) const
Thomas Witkowski's avatar
Thomas Witkowski committed
623
624
625
626
627
  {
    FUNCNAME("ElInfo3d::getSubElementCoords()");

    ERROR_EXIT("Not yet\n");

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
//     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;
Thomas Witkowski's avatar
Thomas Witkowski committed
655
656
  }

657
#endif
658
}