ElInfo3d.cc 17.9 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

Thomas Witkowski's avatar
Thomas Witkowski committed
24
25
26
    macroElement  = const_cast<MacroElement*>( mel);
    element  = const_cast<Element*>( mel->getElement());
    parent = NULL;
27
    level = 0;
28
    elType = const_cast<MacroElement*>(mel)->getElType();
29

Thomas Witkowski's avatar
Thomas Witkowski committed
30
    int vertices = mesh->getGeo(VERTEX);
31

Thomas Witkowski's avatar
Thomas Witkowski committed
32
33
34
    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++)
Thomas Witkowski's avatar
Thomas Witkowski committed
36
	coord[i] = mel->coord[i];
37
38
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
39
    int neighbours = mesh->getGeo(NEIGH);
40

Thomas Witkowski's avatar
Thomas Witkowski committed
41
42
    if (fillFlag.isSet(Mesh::FILL_OPP_COORDS) || 
        fillFlag.isSet(Mesh::FILL_NEIGH)) {
43

Thomas Witkowski's avatar
Thomas Witkowski committed
44
      fill_opp_coords.setFlags(fillFlag & Mesh::FILL_OPP_COORDS);
45
      for (int i = 0; i < neighbours; i++) {
46
	if ((mnb = const_cast<MacroElement*>( mel->getNeighbour(i)))) {
Thomas Witkowski's avatar
Thomas Witkowski committed
47
48
	  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) {
Thomas Witkowski's avatar
Thomas Witkowski committed
54
55
	      neighbour[i]      = const_cast<Element*>( nb->getChild(0));
	      nb = const_cast<Element*>( neighbour[i]);
56
	    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
57
58
	      neighbour[i]      = const_cast<Element*>( nb->getChild(1));
	      nb = const_cast<Element*>( neighbour[i]);
59
	    }
60
	    k = oppVertex[i] = 3;
61
62
63
	    if (fill_opp_coords.isAnySet()) {
	      /* always edge between vertices 0 and 1 is bisected! */
	      if (mnb->getElement()->isNewCoordSet())
Thomas Witkowski's avatar
Thomas Witkowski committed
64
		oppCoord[i] = *(mnb->getElement()->getNewCoord());
65
	      else
Thomas Witkowski's avatar
Thomas Witkowski committed
66
		oppCoord[i] = (mnb->coord[0] + mnb->coord[1]) * 0.5;
67
68
	    }
	  } else {
69
	    if  (fill_opp_coords.isAnySet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
70
	      oppCoord[i] = mnb->coord[k];
71
72
	    }
	  }
73
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
74
	  neighbour[i] = NULL;
75
76
77
78
	}
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
79
80
81
    if (fillFlag.isSet(Mesh::FILL_BOUND)) {
      for (int i = 0; i < element->getGeo(BOUNDARY); i++)
	boundary[i] = mel->getBoundary(i);     
82

Thomas Witkowski's avatar
Thomas Witkowski committed
83
84
      for (int i = 0; i < element->getGeo(PROJECTION); i++)
	projection[i] = mel->getProjection(i);
85
    }
86

Thomas Witkowski's avatar
Thomas Witkowski committed
87
    if (fillFlag.isSet(Mesh::FILL_ORIENTATION)) {
88
89
90
      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()");

Thomas Witkowski's avatar
Thomas Witkowski committed
177
    DimVec<WorldVector<double> > edge(mesh->getDim(), NO_INIT);
178
179
180
    WorldVector<double> x;
    double  x0, det, det0, det1, det2;
  
Thomas Witkowski's avatar
Thomas Witkowski committed
181
    static DimVec<double> vec(mesh->getDim(), NO_INIT);
182

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

Thomas Witkowski's avatar
Thomas Witkowski committed
185
    int dim = mesh->getDim();
186

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++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
196
      x0 = coord[dim][j];
197
      x[j] = xy[j] - x0;
198
199

      for (int i = 0; i < dim; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
200
	edge[i][j] = coord[i][j] - x0;
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
    }

    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

Thomas Witkowski's avatar
Thomas Witkowski committed
266
267
    int neighbours = mesh->getGeo(NEIGH);
    int vertices = mesh->getGeo(VERTEX);
268
  
Thomas Witkowski's avatar
Thomas Witkowski committed
269
    if (fillFlag.isSet(Mesh::FILL_NEIGH) || fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {
270
      Tetrahedron *nb;
Thomas Witkowski's avatar
Thomas Witkowski committed
271
      Flag fill_opp_coords = fillFlag & Mesh::FILL_OPP_COORDS;
272
273
      
      for (int ineigh = 0; ineigh < neighbours; ineigh++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
274
	if ((nb = dynamic_cast<Tetrahedron*>(const_cast<Element*>(neighbour[ineigh])))) {
275
	  int ov = oppVertex[ineigh];
276
277
278
279
	  if (ov < 2 && nb->getFirstChild()) {
	    if (fill_opp_coords != Flag(0)) {
	      int k = -1;
	      for (int j = 0; j < vertices; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
280
		if (element->getDOF(j) == nb->getDOF(1 - ov)) 
281
282
283
		  k = j;
	      
	      if (k == -1) {
284
		for (int j = 0; j < vertices; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
285
		  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
	      
	      if (nb->isNewCoordSet())
Thomas Witkowski's avatar
Thomas Witkowski committed
291
		oppCoord[ineigh] = *(nb->getNewCoord());
292
	      else
Thomas Witkowski's avatar
Thomas Witkowski committed
293
		for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
294
		  oppCoord[ineigh][j] = (oppCoord[ineigh][j] + coord[k][j]) / 2;
295
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
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;
Thomas Witkowski's avatar
Thomas Witkowski committed
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

Thomas Witkowski's avatar
Thomas Witkowski committed
365
366
367
368
    element = const_cast<Element*>( el_old->getChild(ichild));
    macroElement = elInfoOld->macroElement;
    fillFlag = fillFlag_local;
    parent  = el_old;
369
    level = elInfoOld->level + 1;
370
371
    iChild = ichild;
    int el_type_local = (dynamic_cast<const ElInfo3d*>(elInfoOld))->getType();
372
    elType = (el_type_local + 1) % 3;
373

Thomas Witkowski's avatar
Thomas Witkowski committed
374
    TEST_EXIT_DBG(element)("missing child %d?\n", ichild);
375

Thomas Witkowski's avatar
Thomas Witkowski committed
376
    if (fillFlag_local.isAnySet()) {
377
      cvg = Tetrahedron::childVertex[el_type_local];
378
      cv = const_cast<int*>(cvg[ichild]);
379
      ochild = 1 - ichild;
380
381
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
382
383
384
    if (fillFlag_local.isSet(Mesh::FILL_COORDS) || 
	fillFlag.isSet(Mesh::FILL_DET) ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
385
      for (int i = 0; i < 3; i++) {
386
	for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
387
	  coord[i][j] = elInfoOld->coord[cv[i]][j];
388
      }
389
      if (el_old->getNewCoord()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
390
	coord[3] = *(el_old->getNewCoord());
391
      } else {
392
	for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
393
	  coord[3][j] = (elInfoOld->coord[0][j] + elInfoOld->coord[1][j]) / 2;
394
395
      }
    }
396

Thomas Witkowski's avatar
Thomas Witkowski committed
397
398
    if (fillFlag_local.isSet(Mesh::FILL_NEIGH) || 
	fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {
399

Thomas Witkowski's avatar
Thomas Witkowski committed
400
401
      FixVec<Element*, NEIGH> *neigh_local = &neighbour;
      const FixVec<Element*, NEIGH> *neigh_old = &elInfoOld->neighbour;
402
403

      Flag fill_opp_coords;
Thomas Witkowski's avatar
Thomas Witkowski committed
404
      fill_opp_coords.setFlags(fillFlag_local & Mesh::FILL_OPP_COORDS);
405
406
407
408
409
410
411
412
413
414
      
      /*----- 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()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
415
	      oppCoord[0]= *(nb->getNewCoord());
416
417
	    } else {
	      int k = cvg[ochild][1];
418
	      for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
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++)
Thomas Witkowski's avatar
Thomas Witkowski committed
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
	      
	      if (dof == nbk->getDOF(1)) {
		ov = 1;
		if (nbk->getChild(0)) {
		  if (fill_opp_coords.isAnySet()) {
		    if (nbk->getNewCoord()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
456
		      oppCoord[i] = *(nbk->getNewCoord());
457
		    } else {
458
		      for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
459
460
			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++)
Thomas Witkowski's avatar
Thomas Witkowski committed
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
		
		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()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
501
			oppCoord[i] = *(nbk->getNewCoord());
502
		      } else {
503
			for (int j = 0; j < dimOfWorld; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
504
505
			  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++)
Thomas Witkowski's avatar
Thomas Witkowski committed
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++)
Thomas Witkowski's avatar
Thomas Witkowski committed
545
	    oppCoord[3][j] = elInfoOld->oppCoord[ochild][j];
546
547
      }
    }
548

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

562
      ce = const_cast<int*>(Tetrahedron::childEdge[el_type_local][ichild]);
563
      for (int iedge = 0; iedge < 4; iedge++)
Thomas Witkowski's avatar
Thomas Witkowski committed
564
565
566
	boundary[geoFace + iedge] = elInfoOld->getBoundary(geoFace + ce[iedge]);      
      for (int iedge = 4; iedge < 6; iedge++)
	boundary[geoFace + iedge] = elInfoOld->getBoundary(5 - cv[iedge - 3]);
567

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

585
    
Thomas Witkowski's avatar
Thomas Witkowski committed
586
    if (fillFlag.isSet(Mesh::FILL_ORIENTATION)) {
587
588
      orientation = 
	(dynamic_cast<ElInfo3d*>(const_cast<ElInfo*>(elInfoOld)))->orientation 
589
590
591
592
	* Tetrahedron::childOrientation[el_type_local][ichild];
    }
  }

593

594
  mtl::dense2D<double>& ElInfo3d::getSubElemCoordsMat(int degree) const
Thomas Witkowski's avatar
Thomas Witkowski committed
595
  {
596
597
598
    FUNCNAME("ElInfo3d::getSubElemCoordsMat()");

    ERROR_EXIT("Not yet implemented!\n");
599
600

    return subElemMatrices[degree][refinementPath];
Thomas Witkowski's avatar
Thomas Witkowski committed
601
602
  }

603

604
  mtl::dense2D<double>& ElInfo3d::getSubElemCoordsMat_so(int degree) const
Thomas Witkowski's avatar
Thomas Witkowski committed
605
  {
606
607
608
    FUNCNAME("ElInfo3d::getSubElemCoordsMat_so()");

    ERROR_EXIT("Not yet implemented!\n");
609
610

    return subElemMatrices[degree][refinementPath];
Thomas Witkowski's avatar
Thomas Witkowski committed
611
612
  }

613
}