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
    info_stack[stack_used]++;
348
    elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
349
350
    stack_used++;

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

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

360

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

368

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

373
374
375
376
377
378
379
380
381
382
383
384
    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 */
385
      Element *el = elinfo_stack[stack_used]->getElement();
386
387

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


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

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

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

425
    return elinfo_stack[stack_used];
426
427
  }

428

429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
  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;
  }

447

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

452
    Element *el2;
453
454
455
456
457
458
459
460
461
    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}}};

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

    Parametric *parametric = traverse_mesh->getParametric();

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

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

472
    TEST_EXIT_DBG(elinfo_stack[stack_used]->getFillFlag().isSet(Mesh::FILL_NEIGH))
473
      ("FILL_NEIGH not set");
474
    Element *el = elinfo_stack[stack_used]->getElement();
475
476
477
478
479
    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)) {
480
	if (stack_used >= stack_size - 1)
481
482
	  enlargeTraverseStack();
	i = 1 - neighbour;
483
	elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
484
	info_stack[stack_used] = i + 1;
485
486
487
488
489
490
491
492
	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;
493
    save_stack_used = stack_used;
494
495
496

    nb = neighbour;

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

    /* save hierarchy information about current element */

508
509
    for (i = stack_used; i <= save_stack_used; i++) {
      save_info_stack[i] = info_stack[i];
510
511
512
513
514
515
516
517
518
      *(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);
519
      if (traverse_mel == NULL)  
520
	return NULL;
521
    
522
      if (nb < 2 && save_stack_used > 1)
523
	stack2_used = 2;                /* go down one level in OLD hierarchy */
524
      else
525
	stack2_used = 1;
526
      
527
528
529
530
531
532
      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
533
    } else {                                                /* goto other child */
534
      stack2_used = stack_used + 1;
535
      if (save_stack_used > stack2_used)
536
	stack2_used++;                  /* go down one level in OLD hierarchy */
537

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

541
      if (stack_used >= stack_size-1)
542
543
	enlargeTraverseStack();
      i = 2 - info_stack[stack_used];
544
545
      info_stack[stack_used] = i+1;
      elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
546
547
548
549
550
551
552
553
554
555
556
557
558
      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();
559
560
	info_stack[stack_used] = 2-nb;
	elinfo_stack[stack_used+1]->fillElInfo(1 - nb, elinfo_stack[stack_used]);
561
562
563
564
565
566
567
568
	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
569
	TEST_EXIT_DBG(el->getChild(0))("invalid new refinement?\n");
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594

	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;
	}

595
596
	info_stack[stack_used] = i+1;
	if (stack_used >= stack_size-1)
597
	  enlargeTraverseStack();
598
	elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
599
600
601
602
603
604
605
606
607
608
609
610
611
	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)) &&
612
613
614
615
616
617
	      !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();
	  } 
618
619
620
621
622
623
624
625
626
627
628
629
630
	}

      }
      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",
631
	  neighbour, old_elinfo->getElement()->getIndex(), reinterpret_cast<void*>( old_elinfo->getElement()));
632
      MSG(" originally: neighbour %d of element %d at %p\n",
633
	  sav_neighbour, sav_index, reinterpret_cast<void*>( old_elinfo->getElement()));
634
      MSG(" got element %d at %p with opp_vertex %d neigh %d\n",
635
	  elinfo->getElement()->getIndex(), reinterpret_cast<void*>( elinfo->getElement()),
636
	  opp_vertex, (elinfo->getNeighbour(opp_vertex))?(elinfo->getNeighbour(opp_vertex))->getIndex():-1);
637
      TEST_EXIT_DBG(elinfo->getNeighbour(opp_vertex) == old_elinfo->getElement())
638
639
640
641
642
643
644
645
646
	("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;  /* ??? */

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

    return elinfo;
654
655
  }

656

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

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

    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] */

668
    int sav_index, sav_neighbour = neighbour;
669

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

    Parametric *parametric = traverse_mesh->getParametric();

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

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

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

    /* first, goto to leaf level, if necessary... */
    if (!(el->isLeaf()) && (neighbour < 2)) {
687
      if (stack_used >= static_cast<int>( elinfo_stack.size())-1) enlargeTraverseStack();
688
      i = 1 - neighbour;
689
      elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
690
691
692
693
694
695
696
      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
697
698
    save_stack_used   = stack_used;
    for (i=0; i<=stack_used; i++)
699
      save_info_stack[i] = info_stack[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
700
    for (i=0; i<=stack_used; i++)
701
702
703
704
705
706
707
708
709
      (*(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
710
    /* element of the OLD hierarchy branch to the new branch                    */
711
712
713
714
    /****************************************************************************/

    nb = neighbour;

715
716
717
718
719
720
721
722
    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);
    }
723
724
725
726
727
728
729
730
731
732
733

    /****************************************************************************/
    /* 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 */
734
      } else {
735
736
737
	stack2_used = 1;
      }
      elinfo2 = save_elinfo_stack[stack2_used];
738
      el2 = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo2->getElement()));
739
740
741

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

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


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

    /****************************************************************************/
    /* 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
772
    /* new hierarchy branch to the OLD branch.                                  */
773
774
775
    /****************************************************************************/

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

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

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

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

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

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

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

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

    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
824
	  elinfo->getNeighbour(opp_vertex)?elinfo->getNeighbour(opp_vertex)->getIndex():-1);
825
      TEST_EXIT_DBG(elinfo->getNeighbour(opp_vertex) == old_elinfo->getElement())
826
827
	("didn't succeed !?!?!?");
    }
828

829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
    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);
844
845
846
847

      elinfo->fillDetGrdLambda();
    }

848
    return elinfo;
849
850
  }

851

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

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

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

863

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

869
870
871
872
873
874
    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);
875
  }
876

877
}