Traverse.cc 26.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
#include "Traverse.h"
#include "Element.h"
#include "Line.h"
#include "Triangle.h"
#include "Tetrahedron.h"
#include "ElInfo.h"
#include "Mesh.h"
#include "Flag.h"
#include "MacroElement.h"
#include "FixVec.h"
#include "Parametric.h"
12
#include "OpenMP.h"
13
14
15

namespace AMDiS {

16
  ElInfo* TraverseStack::traverseFirst(Mesh *mesh, int level, Flag fill_flag)
17
  {
18
    FUNCNAME("TraverseStack::traverseFirst()");
19

20
21
22
    traverse_mesh = mesh;
    traverse_level = level;
    traverse_fill_flag = fill_flag; 
23
24
25
26
27
28

    if (stack_size < 1)
      enlargeTraverseStack();

    Flag FILL_ANY = Mesh::getFillAnyFlag(mesh->getDim());

29
    for (int i = 0; i < stack_size; i++)
30
31
32
33
34
      elinfo_stack[i]->setFillFlag(fill_flag & FILL_ANY);

    elinfo_stack[0]->setMesh(mesh);
    elinfo_stack[1]->setMesh(mesh);

35
    if (fill_flag.isSet(Mesh::CALL_LEAF_EL_LEVEL)) {
36
      TEST_EXIT_DBG(level >= 0)("invalid level: %d\n", level);   
37
    }
38
39

    traverse_mel = NULL;
40
    stack_used = 0;
41

42
    return traverseNext(NULL);
43
44
  }

45

46
47
  ElInfo* TraverseStack::traverseNext(ElInfo* elinfo_old)
  {
48
    FUNCNAME("TraverseStack::traverseNext()");
49

50
    ElInfo *elinfo = NULL;
51
52
53
54
55
56
    Parametric *parametric = traverse_mesh->getParametric();

    if (stack_used) {
      if (parametric) 
	elinfo_old = parametric->removeParametricInfo(elinfo_old);

57
      TEST_EXIT_DBG(elinfo_old == elinfo_stack[stack_used])("invalid old elinfo\n");
58
    } else {
59
      TEST_EXIT_DBG(elinfo_old == NULL)("invalid old elinfo != nil\n");
60
61
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
62
    if (traverse_fill_flag.isSet(Mesh::CALL_LEAF_EL)) {
63
      elinfo = traverseLeafElement();
Thomas Witkowski's avatar
Thomas Witkowski committed
64
    } else {
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
      if (traverse_fill_flag.isSet(Mesh::CALL_LEAF_EL_LEVEL))
	elinfo = traverseLeafElementLevel();
      else if (traverse_fill_flag.isSet(Mesh::CALL_EL_LEVEL))
	elinfo = traverseElementLevel();
      else if (traverse_fill_flag.isSet(Mesh::CALL_MG_LEVEL))
	elinfo = traverseMultiGridLevel();
      else
	if (traverse_fill_flag.isSet(Mesh::CALL_EVERY_EL_PREORDER))
	  elinfo = traverseEveryElementPreorder();
	else if (traverse_fill_flag.isSet(Mesh::CALL_EVERY_EL_INORDER))
	  elinfo = traverseEveryElementInorder();
	else if (traverse_fill_flag.isSet(Mesh::CALL_EVERY_EL_POSTORDER))
	  elinfo = traverseEveryElementPostorder();
	else
	  ERROR_EXIT("invalid traverse_flag\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
80
    }
81

82
    if (elinfo) {
83
      if (parametric)
84
85
86
87
	elinfo = parametric->addParametricInfo(elinfo);
      elinfo->fillDetGrdLambda();
    }

88
    return elinfo;
Thomas Witkowski's avatar
Thomas Witkowski committed
89
90
  }

91

92
93
  void TraverseStack::enlargeTraverseStack()
  {
94
95
96
    FUNCNAME("TraverseStack::enlargeTraverseStack()");

    int new_stack_size = stack_size + 10;
97
98

    elinfo_stack.resize(new_stack_size, NULL);
99
 
100
    // create new elinfos
101
102
    for (int i = stack_size; i < new_stack_size; i++) {
      TEST_EXIT_DBG(elinfo_stack[i] == NULL)("???\n");
103
104
105
      elinfo_stack[i] = traverse_mesh->createNewElInfo();
    }

106
107
    if (stack_size > 0)
      for (int i = stack_size; i < new_stack_size; i++)
108
	elinfo_stack[i]->setFillFlag(elinfo_stack[0]->getFillFlag());
109
110
111
112
113

    info_stack.resize(new_stack_size);
    save_elinfo_stack.resize(new_stack_size, NULL);

    // create new elinfos
114
115
    for (int i = stack_size; i < new_stack_size; i++) {
      TEST_EXIT_DBG(save_elinfo_stack[i] == NULL)("???\n");
116
117
118
119
120
121
122
123
124
125
      save_elinfo_stack[i] = traverse_mesh->createNewElInfo();
    }
    save_info_stack.resize(new_stack_size);

    stack_size = new_stack_size;    
  }

  
  ElInfo* TraverseStack::traverseLeafElement()
  {
126
    FUNCNAME("TraverseStack::traverseLeafElement()");
127

128
    Element *el = NULL;
129

130
131
    if (stack_used == 0) {   /* first call */
      currentMacro = traverse_mesh->firstMacroElement();
132

133
      while (((*currentMacro)->getIndex() % maxThreads != myThreadId) &&
Thomas Witkowski's avatar
Thomas Witkowski committed
134
      	     (currentMacro != traverse_mesh->endOfMacroElements())) {
135
136
137
      	currentMacro++;
      }

138
      if (currentMacro == traverse_mesh->endOfMacroElements())
139
	return NULL;
140

141
142
143
144
      traverse_mel = *currentMacro;
      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
      info_stack[stack_used] = 0;
145

146
      el = elinfo_stack[stack_used]->getElement();
147
      if (el == NULL || el->getFirstChild() == NULL)
148
	return elinfo_stack[stack_used];
149
150
151
152
    } else {
      el = elinfo_stack[stack_used]->getElement();
      
      /* go up in tree until we can go down again */
153
154
155
156
157
      while ((stack_used > 0) && 
	     ((info_stack[stack_used] >= 2) || (el->getFirstChild() == NULL))) {
	stack_used--;
	el = elinfo_stack[stack_used]->getElement();
      }
158
159
160
      
      /* goto next macro element */
      if (stack_used < 1) {
161
162
	do {	
	  currentMacro++;
Thomas Witkowski's avatar
Thomas Witkowski committed
163
	} while ((currentMacro != traverse_mesh->endOfMacroElements()) && 
164
		 ((*currentMacro)->getIndex() % maxThreads != myThreadId));
165

166
	if (currentMacro == traverse_mesh->endOfMacroElements())
167
	  return NULL;
168
169

	traverse_mel = *currentMacro;	
170
171
	stack_used = 1;
	elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
172
	info_stack[stack_used] = 0;	
173
	el = elinfo_stack[stack_used]->getElement();
174
175

	if (el == NULL || el->getFirstChild() == NULL)
176
	  return elinfo_stack[stack_used];
177
      }
178
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
179
   
180
    /* go down tree until leaf */
181
    while (el->getFirstChild()) {
182
      if (stack_used >= stack_size - 1) 
183
184
	enlargeTraverseStack();
      
Thomas Witkowski's avatar
Thomas Witkowski committed
185
186
      int i = info_stack[stack_used];
      el = const_cast<Element*>(((i == 0) ? el->getFirstChild() : el->getSecondChild()));
187
      info_stack[stack_used]++;
188
      elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
189
      stack_used++;
190

191
      TEST_EXIT_DBG(stack_used < stack_size)
192
	("stack_size=%d too small, level=(%d,%d)\n",
Thomas Witkowski's avatar
Thomas Witkowski committed
193
194
195
	 stack_size, elinfo_stack[stack_used]->getLevel());
      
      info_stack[stack_used] = 0;
196
    }
197

198
199
200
    return elinfo_stack[stack_used];
  }

201

202
203
  ElInfo* TraverseStack::traverseLeafElementLevel()
  {
204
    FUNCNAME("TraverseStack::traverseLeafElementLevel()");
205

206
    ERROR_EXIT("not yet implemented\n");
207

208
209
210
    return NULL;
  }

211

212
213
  ElInfo* TraverseStack::traverseElementLevel()
  {
214
    FUNCNAME("TraverseStack::traverseElementLevel()");
215
216
217
218
219
220

    ElInfo *elInfo;
    do {
      elInfo = traverseEveryElementPreorder();
    } while (elInfo != NULL && elInfo->getLevel() != traverse_level);

221
    return elInfo;
222
223
  }

224

225
226
  ElInfo* TraverseStack::traverseMultiGridLevel()
  {
227
    FUNCNAME("TraverseStack::traverseMultiGridLevel()");
228

229
230
231
    if (stack_used == 0) {   /* first call */
      currentMacro = traverse_mesh->firstMacroElement();
      traverse_mel = *currentMacro;
232
233
      if (traverse_mel == NULL)  
	return NULL;
234
235
236
237
238
239
240
241
      
      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
      info_stack[stack_used] = 0;
      
      if ((elinfo_stack[stack_used]->getLevel() == traverse_level) ||
	  (elinfo_stack[stack_used]->getLevel() < traverse_level && 
	   elinfo_stack[stack_used]->getElement()->isLeaf()))
242
	return elinfo_stack[stack_used];
243
    }
244
  
245
    Element *el = elinfo_stack[stack_used]->getElement();
246
247

    /* go up in tree until we can go down again */
248
249
250
251
252
    while ((stack_used > 0) && 
	   ((info_stack[stack_used] >= 2) || (el->getFirstChild()==NULL))) {
      stack_used--;
      el = elinfo_stack[stack_used]->getElement();
    }
253
254
255
256
257


    /* goto next macro element */
    if (stack_used < 1) {
      currentMacro++;
258
      if (currentMacro == traverse_mesh->endOfMacroElements()) 
259
	return NULL;
260

261
      traverse_mel = *currentMacro;
262
263
264
265
      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
      info_stack[stack_used] = 0;

266
267
268
      if ((elinfo_stack[stack_used]->getLevel() == traverse_level) ||
	  (elinfo_stack[stack_used]->getLevel() < traverse_level && 
	   elinfo_stack[stack_used]->getElement()->isLeaf()))
269
	return elinfo_stack[stack_used];
270
271
272
273
274
    }


    /* go down tree */

275
276
277
278
    if (stack_used >= stack_size - 1) 
      enlargeTraverseStack();

    int i = info_stack[stack_used];
279
280
    info_stack[stack_used]++;

281
    elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
282
283
    stack_used++;

284
    TEST_EXIT_DBG(stack_used < stack_size)
285
286
287
288
289
      ("stack_size=%d too small, level=%d\n",
       stack_size, elinfo_stack[stack_used]->getLevel());

    info_stack[stack_used] = 0;
  
290
291
292
    if ((elinfo_stack[stack_used]->getLevel() == traverse_level) ||
	(elinfo_stack[stack_used]->getLevel() < traverse_level && 
	 elinfo_stack[stack_used]->getElement()->isLeaf()))
293
      return elinfo_stack[stack_used];
294
295
296
297

    return traverseMultiGridLevel();
  }

298

299
300
  ElInfo* TraverseStack::traverseEveryElementPreorder()
  {
301
    FUNCNAME("TraverseStack::traverseEveryElementPreorder()");
302

303
304
305
306
307
308
309
310
311
    if (stack_used == 0) {   /* first call */
      currentMacro = traverse_mesh->firstMacroElement();
      traverse_mel = *currentMacro;
      if (traverse_mel == NULL)  
	return NULL;
      
      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
      info_stack[stack_used] = 0;
312

313
      return elinfo_stack[stack_used];
314
    }
315
  
316
    Element *el = elinfo_stack[stack_used]->getElement();
317
318

    /* go up in tree until we can go down again */
319
320
321
322
323
    while (stack_used > 0 && 
	   (info_stack[stack_used] >= 2 || el->getFirstChild() == NULL)) {
      stack_used--;
      el = elinfo_stack[stack_used]->getElement();
    }
324
325
326
327
328


    /* goto next macro element */
    if (stack_used < 1) {
      currentMacro++;
329
      if (currentMacro == traverse_mesh->endOfMacroElements()) 
330
	return NULL;
331
332
333
334
335
336
      traverse_mel = *currentMacro;

      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
      info_stack[stack_used] = 0;

337
      return elinfo_stack[stack_used];
338
339
340
341
342
    }


    /* go down tree */

343
    if (stack_used >= stack_size - 1) 
344
345
346
      enlargeTraverseStack();

    int i = info_stack[stack_used];
347
348
    info_stack[stack_used]++;

349
    elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
350
351
    stack_used++;

352
    TEST_EXIT_DBG(stack_used < stack_size)
353
354
355
356
357
      ("stack_size=%d too small, level=%d\n",
       stack_size, elinfo_stack[stack_used]->getLevel());

    info_stack[stack_used] = 0;
  
358
    return elinfo_stack[stack_used];
359
360
  }

361

362
363
364
365
366
367
368
  ElInfo* TraverseStack::traverseEveryElementInorder()
  {
    FUNCNAME("TraverseStack::traverseEveryElementInorder");
    ERROR_EXIT("not yet implemented\n");
    return NULL;
  }

369

370
371
  ElInfo* TraverseStack::traverseEveryElementPostorder()
  {
372
    FUNCNAME("TraverseStack::traverseEveryElementPostorder()");
373

374
375
376
377
378
379
380
381
382
383
384
385
    if (stack_used == 0) {   /* first call */
      currentMacro = traverse_mesh->firstMacroElement();
      if (currentMacro == traverse_mesh->endOfMacroElements()) 
	return NULL;
      traverse_mel = *currentMacro;
      
      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
      info_stack[stack_used] = 0;
      
      //return(elinfo_stack[stack_used]);
    } else { /* don't go up on first call */
386
      Element *el = elinfo_stack[stack_used]->getElement();
387
388

      /* go up in tree until we can go down again */          /* postorder!!! */
389
390
      while (stack_used > 0 && 
	     (info_stack[stack_used] >= 3 || el->getFirstChild() == NULL)) {
391
392
393
	stack_used--;
	el = elinfo_stack[stack_used]->getElement();
      }
394
395
396
397
398


      /* goto next macro element */
      if (stack_used < 1) {
	currentMacro++;
399
400
	if (currentMacro == traverse_mesh->endOfMacroElements()) 
	  return NULL;
401
402
403
404
405
406
407
408
409
410
411
	traverse_mel = *currentMacro;

	stack_used = 1;
	elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
	info_stack[stack_used] = 0;

	/*    return(elinfo_stack+stack_used); */
      }
    }
    /* go down tree */

412
413
    while (elinfo_stack[stack_used]->getElement()->getFirstChild() &&
	   info_stack[stack_used] < 2) {
414
      if (stack_used >= stack_size-1) 
415
416
	enlargeTraverseStack();

417
      int i = info_stack[stack_used];
418
      info_stack[stack_used]++;
419
      elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
420
421
422
423
424
425
      stack_used++;
      info_stack[stack_used] = 0;
    }
  
    info_stack[stack_used]++;      /* postorder!!! */

426
    return elinfo_stack[stack_used];
427
428
  }

429

430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
  ElInfo* TraverseStack::traverseNeighbour(ElInfo* elinfo_old, int neighbour)
  {
    int dim = elinfo_old->getMesh()->getDim();
    switch(dim) {
    case 1:
      ERROR_EXIT("invalid dim\n");
      break;
    case 2:
      return traverseNeighbour2d(elinfo_old, neighbour);
      break;
    case 3:
      return traverseNeighbour3d(elinfo_old, neighbour);
      break;
    default: ERROR_EXIT("invalid dim\n");
    }
    return NULL;
  }

448

449
450
  ElInfo* TraverseStack::traverseNeighbour3d(ElInfo* elinfo_old, int neighbour)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
451
452
    FUNCNAME("TraverseStackLLtraverseNeighbour3d()");

453
    Element *el2;
454
455
456
457
458
459
460
461
462
    ElInfo *old_elinfo, *elinfo, *elinfo2;
    const DegreeOfFreedom *dof;
    int i, nb, opp_vertex, stack2_used, typ = 0;
    int sav_index, sav_neighbour = neighbour;

    static int coarse_nb[3][3][4] = {{{-2,-2,-2,-2}, {-1,2,3,1}, {-1,3,2,0}},
				     {{-2,-2,-2,-2}, {-1,2,3,1}, {-1,2,3,0}},
				     {{-2,-2,-2,-2}, {-1,2,3,1}, {-1,2,3,0}}};

463
    TEST_EXIT_DBG(stack_used > 0)("no current element\n");
464
465
466

    Parametric *parametric = traverse_mesh->getParametric();

467
468
    if (parametric) 
      elinfo_old = parametric->removeParametricInfo(elinfo_old);
469

470
    TEST_EXIT_DBG(elinfo_old == elinfo_stack[stack_used])
471
472
      ("invalid old elinfo\n");

473
    TEST_EXIT_DBG(elinfo_stack[stack_used]->getFillFlag().isSet(Mesh::FILL_NEIGH))
474
      ("FILL_NEIGH not set");
475
    Element *el = elinfo_stack[stack_used]->getElement();
476
477
478
479
480
    sav_index = el->getIndex();

    /* first, goto to leaf level, if necessary... */
    if ((traverse_fill_flag & Mesh::CALL_LEAF_EL).isAnySet()) {
      if ((el->getChild(0)) && (neighbour < 2)) {
481
	if (stack_used >= stack_size - 1)
482
483
	  enlargeTraverseStack();
	i = 1 - neighbour;
484
	elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
485
	info_stack[stack_used] = i + 1;
486
487
488
489
490
491
492
493
	stack_used++;
	info_stack[stack_used] = 0;
	neighbour = 3;
      }
    }

    /* save information about current element and its position in the tree */
    save_traverse_mel = traverse_mel;
494
    save_stack_used = stack_used;
495
496
497

    nb = neighbour;

498
499
    while (stack_used > 1) { /* go up in tree until we can go down again */
      stack_used--;
500
      typ = elinfo_stack[stack_used]->getType();
501
502
503
504
505
      nb = coarse_nb[typ][info_stack[stack_used]][nb];
      if (nb == -1) 
	break;
      TEST_EXIT_DBG(nb >= 0)("invalid coarse_nb %d\n",nb);
    }
506
507
508

    /* save hierarchy information about current element */

509
510
    for (i = stack_used; i <= save_stack_used; i++) {
      save_info_stack[i] = info_stack[i];
511
512
513
514
515
516
517
518
519
      *(save_elinfo_stack[i]) = *(elinfo_stack[i]);
    }

    old_elinfo = save_elinfo_stack[save_stack_used];
    opp_vertex = old_elinfo->getOppVertex(neighbour);

    if (nb >= 0) {                           /* go to macro element neighbour */
      i = traverse_mel->getOppVertex(nb);
      traverse_mel = traverse_mel->getNeighbour(nb);
520
      if (traverse_mel == NULL)  
521
	return NULL;
522
    
523
      if (nb < 2 && save_stack_used > 1)
524
	stack2_used = 2;                /* go down one level in OLD hierarchy */
525
      else
526
	stack2_used = 1;
527
      
528
529
530
531
532
533
      elinfo2 = save_elinfo_stack[stack2_used];
      el2 = elinfo2->getElement();
      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
      info_stack[stack_used] = 0;
      nb = i;
Thomas Witkowski's avatar
Thomas Witkowski committed
534
    } else {                                                /* goto other child */
535
      stack2_used = stack_used + 1;
536
      if (save_stack_used > stack2_used)
537
	stack2_used++;                  /* go down one level in OLD hierarchy */
538

539
540
541
      elinfo2 = save_elinfo_stack[stack2_used];
      el2 = elinfo2->getElement();

542
      if (stack_used >= stack_size-1)
543
544
	enlargeTraverseStack();
      i = 2 - info_stack[stack_used];
545
546
      info_stack[stack_used] = i+1;
      elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
547
548
549
550
551
552
553
554
555
556
557
558
559
      stack_used++;
      info_stack[stack_used] = 0;
      nb = 0;
    }

    elinfo = elinfo_stack[stack_used];
    el = elinfo->getElement();

    while(el->getChild(0)) {
 
      if (nb < 2) {                         /* go down one level in hierarchy */
	if (stack_used >= stack_size-1)
	  enlargeTraverseStack();
560
561
	info_stack[stack_used] = 2-nb;
	elinfo_stack[stack_used+1]->fillElInfo(1 - nb, elinfo_stack[stack_used]);
562
563
564
565
566
567
568
569
	stack_used++;
	info_stack[stack_used] = 0;
	elinfo = elinfo_stack[stack_used];
	el = elinfo->getElement();
	nb = 3;
      }

      if (save_stack_used > stack2_used) { /* `refine' both el and el2 */
Thomas Witkowski's avatar
Thomas Witkowski committed
570
	TEST_EXIT_DBG(el->getChild(0))("invalid new refinement?\n");
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595

	if (el->getDOF(0) == el2->getDOF(0))
	  i = save_info_stack[stack2_used] - 1;
	else if (el->getDOF(1) == el2->getDOF(0))
	  i = 2 - save_info_stack[stack2_used];
	else {
	  if (traverse_mesh->associated(el->getDOF(0, 0), el2->getDOF(0, 0)))
	    i = save_info_stack[stack2_used] - 1;
	  else if (traverse_mesh->associated(el->getDOF(1, 0), el2->getDOF(0, 0)))
	    i = 2 - save_info_stack[stack2_used];
	  else {
	    ERROR_EXIT("no common refinement edge\n");
	  }
	}

	if (el->getChild(0)  &&  
	    (el->getChild(i)->getDOF(1) == el->getDOF(nb) ||
	     traverse_mesh->associated(el->getChild(i)->getDOF(1, 0), el->getDOF(nb, 0)))) 
	  {   /*  el->child[0]  Kuni  22.08.96  */
	    nb = 1;
	  }
	else {
	  nb = 2;
	}

596
597
	info_stack[stack_used] = i+1;
	if (stack_used >= stack_size-1)
598
	  enlargeTraverseStack();
599
	elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
600
601
602
603
604
605
606
607
608
609
610
611
612
	stack_used++;
	info_stack[stack_used] = 0;

	elinfo = elinfo_stack[stack_used];
	el = elinfo->getElement();

	stack2_used++;
	elinfo2 = save_elinfo_stack[stack2_used];
	el2 = elinfo2->getElement();
	if (save_stack_used > stack2_used) {
	  dof = el2->getDOF(1);
	  if (dof != el->getDOF(1) && dof != el->getDOF(2) &&
	      !traverse_mesh->associated(dof[0], el->getDOF(1, 0)) &&
613
614
615
616
617
618
	      !traverse_mesh->associated(dof[0], el->getDOF(2, 0))) {
	  
	    stack2_used++;               /* go down one level in OLD hierarchy */
	    elinfo2 = save_elinfo_stack[stack2_used];
	    el2 = elinfo2->getElement();
	  } 
619
620
621
622
623
624
625
626
627
628
629
630
631
	}

      }
      else {   /* now we're done... */
	elinfo = elinfo_stack[stack_used];
	el = elinfo->getElement();
	break;
      }
    }


    if (elinfo->getNeighbour(opp_vertex) != old_elinfo->getElement()) {
      MSG(" looking for neighbour %d of element %d at %p\n",
632
	  neighbour, old_elinfo->getElement()->getIndex(), reinterpret_cast<void*>( old_elinfo->getElement()));
633
      MSG(" originally: neighbour %d of element %d at %p\n",
634
	  sav_neighbour, sav_index, reinterpret_cast<void*>( old_elinfo->getElement()));
635
      MSG(" got element %d at %p with opp_vertex %d neigh %d\n",
636
	  elinfo->getElement()->getIndex(), reinterpret_cast<void*>( elinfo->getElement()),
637
	  opp_vertex, (elinfo->getNeighbour(opp_vertex))?(elinfo->getNeighbour(opp_vertex))->getIndex():-1);
638
      TEST_EXIT_DBG(elinfo->getNeighbour(opp_vertex) == old_elinfo->getElement())
639
640
641
642
643
644
645
646
647
	("didn't succeed !?!?!?\n");
    }


    if ((traverse_fill_flag & Mesh::CALL_EVERY_EL_POSTORDER).isAnySet())
      info_stack[stack_used] = 3;
    else if ((traverse_fill_flag & Mesh::CALL_EVERY_EL_INORDER).isAnySet())
      info_stack[stack_used] = 1;  /* ??? */

648
649
650
    if (elinfo) {
      if (parametric) 
	elinfo = parametric->addParametricInfo(elinfo);
651
652
      elinfo->fillDetGrdLambda();
    }
653
654

    return elinfo;
655
656
  }

657

658
659
  ElInfo* TraverseStack::traverseNeighbour2d(ElInfo* elinfo_old, int neighbour)
  {
660
661
    FUNCNAME("TraverseStack::traverseNeighbour2d()");

Thomas Witkowski's avatar
Thomas Witkowski committed
662
    Triangle *el, *el2, *sav_el;
663
    ElInfo *old_elinfo, *elinfo, *elinfo2;
664
    int i, nb, opp_vertex, stack2_used;
665
666
667
668

    static int coarse_nb[3][3] = {{-2,-2,-2}, {2,-1,1}, {-1,2,0}};
    /* father.neigh[coarse_nb[i][j]] == child[i-1].neigh[j] */

669
    int sav_index, sav_neighbour = neighbour;
670

671
    TEST_EXIT_DBG(stack_used > 0)("no current element");
Thomas Witkowski's avatar
Thomas Witkowski committed
672
    TEST_EXIT_DBG(traverse_fill_flag.isSet(Mesh::CALL_LEAF_EL))("invalid traverse_fill_flag");
673
674
675
676
677
678

    Parametric *parametric = traverse_mesh->getParametric();

    if (parametric) 
      elinfo_old = parametric->removeParametricInfo(elinfo_old);

679
    TEST_EXIT_DBG(elinfo_old == elinfo_stack[stack_used])("invalid old elinfo");
680
681

    elinfo_stack[stack_used]->testFlag(Mesh::FILL_NEIGH);
682
    el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo_stack[stack_used]->getElement()));
683
    sav_index = el->getIndex();
684
    sav_el = el;
685
686
687

    /* first, goto to leaf level, if necessary... */
    if (!(el->isLeaf()) && (neighbour < 2)) {
688
      if (stack_used >= static_cast<int>( elinfo_stack.size())-1) enlargeTraverseStack();
689
      i = 1 - neighbour;
690
      elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
691
692
693
694
695
696
697
      info_stack[stack_used] = i + 1;
      stack_used++;
      neighbour = 2;
    }

    /* save information about current element and its position in the tree */
    save_traverse_mel = traverse_mel;
Thomas Witkowski's avatar
Thomas Witkowski committed
698
699
    save_stack_used   = stack_used;
    for (i=0; i<=stack_used; i++)
700
      save_info_stack[i] = info_stack[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
701
    for (i=0; i<=stack_used; i++)
702
703
704
705
706
707
708
709
710
      (*(save_elinfo_stack[i])) = (*(elinfo_stack[i]));
    old_elinfo = save_elinfo_stack[stack_used];
    opp_vertex = old_elinfo->getOppVertex(neighbour);


    /****************************************************************************/
    /* First phase: go up in tree until we can go down again.                   */
    /*                                                                          */
    /* During this first phase, nb is the neighbour index which points from an  */
Thomas Witkowski's avatar
Thomas Witkowski committed
711
    /* element of the OLD hierarchy branch to the new branch                    */
712
713
714
715
    /****************************************************************************/

    nb = neighbour;

716
717
718
719
720
721
722
723
    while (stack_used > 1) {
      stack_used--;
      nb = coarse_nb[info_stack[stack_used]][nb];
      if (nb == -1) 
	break;
      
      TEST_EXIT_DBG(nb >= 0)("invalid coarse_nb %d\n",nb);
    }
724
725
726
727
728
729
730
731
732
733
734

    /****************************************************************************/
    /* Now, goto neighbouring element at the local hierarchy entry              */
    /* This is either a macro element neighbour or the other child of parent.   */
    /* initialize nb for second phase (see below)                               */
    /****************************************************************************/

    if (nb >= 0) {                        /* go to macro element neighbour */

      if ((nb < 2) && (save_stack_used > 1)) {
	stack2_used = 2;           /* go down one level in OLD hierarchy */
735
      } else {
736
737
738
	stack2_used = 1;
      }
      elinfo2 = save_elinfo_stack[stack2_used];
739
      el2 = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo2->getElement()));
740
741
742

      i = traverse_mel->getOppVertex(nb);
      traverse_mel = traverse_mel->getNeighbour(nb);
743
744
      if (traverse_mel == NULL)
	return NULL;
745
746
747
      nb = i;
    
      stack_used = 1;
748
      elinfo_stack[stack_used]->fillMacroInfo(const_cast<MacroElement*>( traverse_mel));
749
      info_stack[stack_used] = 0;
750
    } else {                                               /* goto other child */
751
752

      stack2_used = stack_used + 1;
753
      if (save_stack_used > stack2_used)
754
	stack2_used++;               /* go down one level in OLD hierarchy */
755
      
756
      elinfo2 = save_elinfo_stack[stack2_used];
757
      el2 = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo2->getElement()));
758
759


760
      if (stack_used >= stack_size - 1)
761
	enlargeTraverseStack();
762
      
763
      i = 2 - info_stack[stack_used];
764
765
      info_stack[stack_used] = i + 1;
      elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
766
      stack_used++;
767
      nb = 1-i;
768
769
770
771
772
    }

    /****************************************************************************/
    /* Second phase: go down in a new hierarchy branch until leaf level.        */
    /* Now, nb is the neighbour index which points from an element of the       */
Thomas Witkowski's avatar
Thomas Witkowski committed
773
    /* new hierarchy branch to the OLD branch.                                  */
774
775
776
    /****************************************************************************/

    elinfo = elinfo_stack[stack_used];
777
    el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo->getElement()));
778

779
    while (el->getFirstChild()) {
780
781

      if (nb < 2) {   /* go down one level in hierarchy */
782
	if (stack_used >= stack_size-1)
783
	  enlargeTraverseStack();
784
785
	elinfo_stack[stack_used+1]->fillElInfo(1-nb, elinfo_stack[stack_used]);
	info_stack[stack_used] = 2-nb;
786
787
788
789
790
	stack_used++;
	nb = 2;
      }

      if (save_stack_used > stack2_used) { /* `refine' both el and el2 */
Thomas Witkowski's avatar
Thomas Witkowski committed
791
	TEST_EXIT_DBG(el->getFirstChild())("invalid new refinement?");
792
793
794

	/* use child i, neighbour of el2->child[nb-1] */
	i = 2 - save_info_stack[stack2_used];
795
	TEST_EXIT_DBG(i < 2)("invalid OLD refinement?");
796
797
	info_stack[stack_used] = i+1;
	elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
798
799
800
801
	stack_used++;
	nb = i;

	elinfo = elinfo_stack[stack_used];
802
	el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo->getElement()));
803
804

	stack2_used++;
805
806
807
	if (save_stack_used > stack2_used) {
	  stack2_used++;                /* go down one level in OLD hierarchy */
	}
808
	elinfo2 = save_elinfo_stack[stack2_used];
809
	el2 = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo2->getElement()));
810

811
812
      }
      else {   /* now we're done... */
813
	elinfo = elinfo_stack[stack_used];
814
	el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo->getElement()));
815
816
817
818
819
820
821
822
823
824
      }
    }

    if (elinfo->getNeighbour(opp_vertex) != old_elinfo->getElement()) {
      MSG(" looking for neighbour %d of element %d at %8X\n",
	  neighbour, old_elinfo->getElement()->getIndex(), old_elinfo->getElement());
      MSG(" originally: neighbour %d of element %d at %8X\n",
	  sav_neighbour, sav_index, sav_el);
      MSG(" got element %d at %8X with opp_vertex %d neigh %d\n",
	  elinfo->getElement()->getIndex(), elinfo->getElement(), opp_vertex,
Thomas Witkowski's avatar
Thomas Witkowski committed
825
	  elinfo->getNeighbour(opp_vertex)?elinfo->getNeighbour(opp_vertex)->getIndex():-1);
826
      TEST_EXIT_DBG(elinfo->getNeighbour(opp_vertex) == old_elinfo->getElement())
827
828
	("didn't succeed !?!?!?");
    }
829

830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
    if (elinfo->getElement()->getFirstChild()) {
      MSG(" looking for neighbour %d of element %d at %8X\n",
	  neighbour, old_elinfo->getElement()->getIndex(), old_elinfo->getElement());
      MSG(" originally: neighbour %d of element %d at %8X\n",
	  sav_neighbour, sav_index, sav_el);
      MSG(" got element %d at %8X with opp_vertex %d neigh %d\n",
	  elinfo->getElement()->getIndex(), elinfo->getElement(), opp_vertex,
	  elinfo->getNeighbour(opp_vertex)->getIndex());
      MSG("got no leaf element\n");
      WAIT_REALLY;
    }

    if (elinfo) {
      if (parametric) 
	elinfo = parametric->addParametricInfo(elinfo);
845
846
847
848

      elinfo->fillDetGrdLambda();
    }

849
    return elinfo;
850
851
  }

852

853
854
  void TraverseStack::update()
  {
855
    FUNCNAME("TraverseStack::update()");
856

857
858
    TEST_EXIT_DBG(traverse_mesh->getDim() == 3)
      ("Update only in 3d, mesh is d = %d\n", traverse_mesh->getDim());
859

860
    for (int i = stack_used; i > 0; i--)
861
      dynamic_cast<ElInfo3d*>(elinfo_stack[i])->update();
862
863
  }

864

865
  void TraverseStack::fillRefinementPath(ElInfo& elInfo, const ElInfo& upperElInfo)
866
  {
867
868
    int levelDif = elinfo_stack[stack_used]->getLevel() - upperElInfo.getLevel();
    unsigned long rPath = 0;
869

870
871
872
873
874
875
    for (int i = 1; i <= levelDif; i++) 
      if (elinfo_stack[stack_used - levelDif + i]->getIChild())
	rPath = rPath | (1 << (i - 1));

    elInfo.setRefinementPath(rPath);
    elInfo.setRefinementPathLength(levelDif);
876
  }
877

878
}