Traverse.cc 34 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
16
17
18

namespace AMDiS {

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

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

    if (stack_size < 1)
      enlargeTraverseStack();

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

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

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

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

    traverse_mel = NULL;
40
    stack_used = 0;
41

42
    return traverseNext(NULL);
43
44
45
46
  }

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
61
    if (traverse_fill_flag.isSet(Mesh::CALL_LEAF_EL)) {
62
      elinfo = traverseLeafElement();
Thomas Witkowski's avatar
Thomas Witkowski committed
63
    } else {
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
      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
79
    }
80

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

88
    return elinfo;
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
108
109
110
    if (stack_size > 0) {
      for (int i = stack_size; i < new_stack_size; i++) {
	elinfo_stack[i]->setFillFlag(elinfo_stack[0]->getFillFlag());
      }
    }
111
112
113
114
115

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

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

    stack_size = new_stack_size;    
  }

  
  ElInfo* TraverseStack::traverseLeafElement()
  {
128
    FUNCNAME("TraverseStack::traverseLeafElement()");
129

130
    Element *el = NULL;
131

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

135
136
137
138
139
140
      // Search for the first macro element that corresponds to the current thread
      // and that is not hidden. So increment as long as one of this conditions does
      // not hold and the current macro element is not the last one.
      while (((*currentMacro)->getIndex() % maxThreads_ != myThreadId_ ||
	      (*currentMacro)->getElement()->isHidden()) &&
	     currentMacro != traverse_mesh->endOfMacroElements()) {
141
142
143
      	currentMacro++;
      }

144
      if (currentMacro == traverse_mesh->endOfMacroElements()) {
145
	return NULL;
146
      }
147

148
149
150
151
      traverse_mel = *currentMacro;
      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
      info_stack[stack_used] = 0;
152

153
154
      el = elinfo_stack[stack_used]->getElement();
      if ((el == NULL) || (el->getFirstChild() == NULL)) {
155
	return elinfo_stack[stack_used];
156
157
158
159
160
      }
    } else {
      el = elinfo_stack[stack_used]->getElement();
      
      /* go up in tree until we can go down again */
161
162
163
164
165
      while ((stack_used > 0) && 
	     ((info_stack[stack_used] >= 2) || (el->getFirstChild() == NULL))) {
	stack_used--;
	el = elinfo_stack[stack_used]->getElement();
      }
166
167
168
      
      /* goto next macro element */
      if (stack_used < 1) {
169
170
	// Again, search for the next macro element that corresponds to the current
	// thread and that is not hidden.
171
172
	do {	
	  currentMacro++;
173
174
175
	} while (currentMacro != traverse_mesh->endOfMacroElements() && 		 
		 ((*currentMacro)->getIndex() % maxThreads_ != myThreadId_ ||
		  (*currentMacro)->getElement()->isHidden()));
176

177
	if (currentMacro == traverse_mesh->endOfMacroElements()) {
178
	  return NULL;
179
	}
180
181
	traverse_mel = *currentMacro;
	
182
183
184
	stack_used = 1;
	elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
	info_stack[stack_used] = 0;
185
	
186
187
	el = elinfo_stack[stack_used]->getElement();
	if ((el == NULL) || (el->getFirstChild() == NULL)) {
188
	  return (elinfo_stack[stack_used]);
189
190
	}
      }
191
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
192
   
193
    /* go down tree until leaf */
194
    while (el->getFirstChild()) {
195
      if (stack_used >= stack_size - 1) 
196
197
	enlargeTraverseStack();
      
198
199
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
230
231
232
233
234
235
236
237
238
239

      int visitedChildren = info_stack[stack_used];
      if (visitedChildren == 0) {
	// No children of this element have been visited.

	if (el->getFirstChild()->isHidden()) {
	  // The first child is hidden, so omit it.
	  info_stack[stack_used]++;
	  visitedChildren++;

	  // And check the second child.
	  if (el->getSecondChild()->isHidden()) {
	    // Also the second child is hidden, so omit it too and start the
	    // traverse procedure again to go up into the stack.

	    info_stack[stack_used]++;

	    return traverseLeafElement();
	  } else {
	    // Only the first child is hidden, so go down on the second one.
	    el = const_cast<Element*>(el->getSecondChild());
	  }
	} else {
	  // The first child is not hidden, so go down on it.
	  el = const_cast<Element*>(el->getFirstChild());
	}

      } else {
	// The first child has been visited before, so check the second one.
	if (el->getSecondChild()->isHidden()) {
	  // The second one is hidden, so omit it and start the traverse
	  // procedure again to go up into the stack;
	  info_stack[stack_used]++;
	  
	  return traverseLeafElement();
	} else {
	  // The secod child is not hidden, so go down on it.
	  el = const_cast<Element*>(el->getSecondChild());
	}
      }


240
      info_stack[stack_used]++;
241
242
243
      elinfo_stack[stack_used + 1]->fillElInfo(visitedChildren, 
					       elinfo_stack[stack_used]);
      info_stack[stack_used + 1] = 0;
244
      stack_used++;
245

246
      TEST_EXIT_DBG(stack_used < stack_size)
247
	("stack_size=%d too small, level=(%d,%d)\n",
248
	 stack_size, elinfo_stack[stack_used]->getLevel());     
249
    }
250

251
252
253
254
255
    return elinfo_stack[stack_used];
  }

  ElInfo* TraverseStack::traverseLeafElementLevel()
  {
256
    FUNCNAME("TraverseStack::traverseLeafElementLevel()");
257
    ERROR_EXIT("not yet implemented\n");
258

259
260
261
262
263
    return NULL;
  }

  ElInfo* TraverseStack::traverseElementLevel()
  {
264
    FUNCNAME("TraverseStack::traverseElementLevel()");
265
266
267
268
269
270

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

271
    return elInfo;
272
273
274
275
  }

  ElInfo* TraverseStack::traverseMultiGridLevel()
  {
276
    FUNCNAME("TraverseStack::traverseMultiGridLevel()");
277

278
279
280
281
282
283
284
285
286
287
288
289
290
291
    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;
      
      if ((elinfo_stack[stack_used]->getLevel() == traverse_level) ||
	  (elinfo_stack[stack_used]->getLevel() < traverse_level && 
	   elinfo_stack[stack_used]->getElement()->isLeaf()))
	return(elinfo_stack[stack_used]);
    }
292
  
293
    Element *el = elinfo_stack[stack_used]->getElement();
294
295

    /* go up in tree until we can go down again */
296
297
298
299
300
    while ((stack_used > 0) && 
	   ((info_stack[stack_used] >= 2) || (el->getFirstChild()==NULL))) {
      stack_used--;
      el = elinfo_stack[stack_used]->getElement();
    }
301
302
303
304
305


    /* goto next macro element */
    if (stack_used < 1) {
      currentMacro++;
306
307
      if (currentMacro == traverse_mesh->endOfMacroElements()) 
	return(NULL);
308

309
      traverse_mel = *currentMacro;
310
311
312
313
      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
      info_stack[stack_used] = 0;

314
315
316
      if ((elinfo_stack[stack_used]->getLevel() == traverse_level) ||
	  (elinfo_stack[stack_used]->getLevel() < traverse_level && 
	   elinfo_stack[stack_used]->getElement()->isLeaf()))
317
318
319
320
321
322
	return(elinfo_stack[stack_used]);
    }


    /* go down tree */

323
324
325
326
    if (stack_used >= stack_size - 1) 
      enlargeTraverseStack();

    int i = info_stack[stack_used];
327
328
329
330
331
    info_stack[stack_used]++;

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

332
    TEST_EXIT_DBG(stack_used < stack_size)
333
334
335
336
337
      ("stack_size=%d too small, level=%d\n",
       stack_size, elinfo_stack[stack_used]->getLevel());

    info_stack[stack_used] = 0;
  
338
339
340
    if ((elinfo_stack[stack_used]->getLevel() == traverse_level) ||
	(elinfo_stack[stack_used]->getLevel() < traverse_level && 
	 elinfo_stack[stack_used]->getElement()->isLeaf()))
341
342
343
344
345
346
347
      return(elinfo_stack[stack_used]);

    return traverseMultiGridLevel();
  }

  ElInfo* TraverseStack::traverseEveryElementPreorder()
  {
348
    FUNCNAME("TraverseStack::traverseEveryElementPreorder()");
349
350
351

    Element *el;

352
353
354
355
356
357
358
359
360
    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;
361

362
363
      return(elinfo_stack[stack_used]);
    }
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
  
    el = elinfo_stack[stack_used]->getElement();

    /* go up in tree until we can go down again */
    while((stack_used > 0) && 
	  ((info_stack[stack_used] >= 2) || (el->getFirstChild()==NULL)))
      {
	stack_used--;
	el = elinfo_stack[stack_used]->getElement();
      }


    /* goto next macro element */
    if (stack_used < 1) {
      currentMacro++;
379
380
      if (currentMacro == traverse_mesh->endOfMacroElements()) 
	return(NULL);
381
382
383
384
385
386
387
388
389
390
391
392
      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 */

393
    if (stack_used >= stack_size - 1) 
394
395
396
      enlargeTraverseStack();

    int i = info_stack[stack_used];
397
398
399
400
401
    info_stack[stack_used]++;

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

402
    TEST_EXIT_DBG(stack_used < stack_size)
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
      ("stack_size=%d too small, level=%d\n",
       stack_size, elinfo_stack[stack_used]->getLevel());

    info_stack[stack_used] = 0;
  
    return(elinfo_stack[stack_used]);
  }

  ElInfo* TraverseStack::traverseEveryElementInorder()
  {
    FUNCNAME("TraverseStack::traverseEveryElementInorder");
    ERROR_EXIT("not yet implemented\n");
    return NULL;
  }

  ElInfo* TraverseStack::traverseEveryElementPostorder()
  {
    FUNCNAME("TraverseStack::traverseEveryElementPostorder");

    Element *el;
    int i;

    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 */
  
      el = elinfo_stack[stack_used]->getElement();

      /* go up in tree until we can go down again */          /* postorder!!! */
      while((stack_used > 0) && 
	    ((info_stack[stack_used] >= 3) || (el->getFirstChild()==NULL)))
	{
	  stack_used--;
	  el = elinfo_stack[stack_used]->getElement();
	}


      /* goto next macro element */
      if (stack_used < 1) {
	currentMacro++;
	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); */
      }
    }
    /* go down tree */

    while(elinfo_stack[stack_used]->getElement()->getFirstChild()
	  && (info_stack[stack_used] < 2)) {
467
468
469
      if (stack_used >= stack_size-1) 
	enlargeTraverseStack();

470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
      i = info_stack[stack_used];
      info_stack[stack_used]++;
      elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
      stack_used++;
      info_stack[stack_used] = 0;
    }
  
    info_stack[stack_used]++;      /* postorder!!! */

    return(elinfo_stack[stack_used]);
  }

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

  ElInfo* TraverseStack::traverseNeighbour3d(ElInfo* elinfo_old, int neighbour)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
502
503
    FUNCNAME("TraverseStackLLtraverseNeighbour3d()");

504
    Element *el2;
505
506
507
508
509
510
511
512
513
    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}}};

514
    TEST_EXIT_DBG(stack_used > 0)("no current element\n");
515
516
517

    Parametric *parametric = traverse_mesh->getParametric();

518
519
    if (parametric) 
      elinfo_old = parametric->removeParametricInfo(elinfo_old);
520

521
    TEST_EXIT_DBG(elinfo_old == elinfo_stack[stack_used])
522
523
      ("invalid old elinfo\n");

524
    TEST_EXIT_DBG(elinfo_stack[stack_used]->getFillFlag().isSet(Mesh::FILL_NEIGH))
525
      ("FILL_NEIGH not set");
526
    Element *el = elinfo_stack[stack_used]->getElement();
527
528
529
530
531
    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)) {
532
	if (stack_used >= stack_size - 1)
533
534
	  enlargeTraverseStack();
	i = 1 - neighbour;
535
536
	elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
	info_stack[stack_used] = i + 1;
537
538
539
540
541
542
543
544
545
	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;
546
    save_stack_used = stack_used;
547
548
549

    nb = neighbour;

550
551
552
553
554
555
556
557
    while (stack_used > 1) { /* go up in tree until we can go down again */
      stack_used--;
      typ = dynamic_cast<ElInfo3d*>( elinfo_stack[stack_used])->getType();
      nb = coarse_nb[typ][info_stack[stack_used]][nb];
      if (nb == -1) 
	break;
      TEST_EXIT_DBG(nb >= 0)("invalid coarse_nb %d\n",nb);
    }
558
559
560

    /* save hierarchy information about current element */

561
562
    for (i = stack_used; i <= save_stack_used; i++) {
      save_info_stack[i] = info_stack[i];
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
      *(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);
      if (traverse_mel == NULL)  return(NULL);
    
      if ((nb < 2) && (save_stack_used > 1)) {
	stack2_used = 2;                /* go down one level in OLD hierarchy */
      }
      else {
	stack2_used = 1;
      }
      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
588
    } else {                                                /* goto other child */
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
      stack2_used = stack_used + 1;
      if (save_stack_used > stack2_used) {
	stack2_used++;                  /* go down one level in OLD hierarchy */
      }
      elinfo2 = save_elinfo_stack[stack2_used];
      el2 = elinfo2->getElement();

      if (stack_used >= stack_size-1)
	enlargeTraverseStack();
      i = 2 - info_stack[stack_used];
      info_stack[stack_used] = i+1;
      elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
      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();
	info_stack[stack_used] = 2-nb;
Thomas Witkowski's avatar
Thomas Witkowski committed
617
	elinfo_stack[stack_used+1]->fillElInfo(1 - nb, elinfo_stack[stack_used]);
618
619
620
621
622
623
624
625
	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 */
626
	TEST_EXIT_DBG(el->getChild(0))("invalid NEW refinement?\n");
627
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
655
656
657
658
659
660
661
662
663
664
665
666
667
668

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

	info_stack[stack_used] = i+1;
	if (stack_used >= stack_size-1)
	  enlargeTraverseStack();
	elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
	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)) &&
669
670
671
672
673
674
	      !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();
	  } 
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
	}

      }
      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",
	  neighbour, old_elinfo->getElement()->getIndex(), reinterpret_cast<void*>( old_elinfo->getElement()));
      MSG(" originally: neighbour %d of element %d at %p\n",
	  sav_neighbour, sav_index, reinterpret_cast<void*>( old_elinfo->getElement()));
      MSG(" got element %d at %p with opp_vertex %d neigh %d\n",
	  elinfo->getElement()->getIndex(), reinterpret_cast<void*>( elinfo->getElement()),
	  opp_vertex, (elinfo->getNeighbour(opp_vertex))?(elinfo->getNeighbour(opp_vertex))->getIndex():-1);
694
      TEST_EXIT_DBG(elinfo->getNeighbour(opp_vertex) == old_elinfo->getElement())
695
696
697
698
699
700
701
702
703
	("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;  /* ??? */

704
705
706
    if (elinfo) {
      if (parametric) 
	elinfo = parametric->addParametricInfo(elinfo);
707
708
709
710
711
712
713
      elinfo->fillDetGrdLambda();
    }
    return(elinfo);
  }

  ElInfo* TraverseStack::traverseNeighbour2d(ElInfo* elinfo_old, int neighbour)
  {
714
715
    FUNCNAME("TraverseStack::traverseNeighbour2d()");

716
    Element *el, *el2, *sav_el;
717
    ElInfo *old_elinfo, *elinfo, *elinfo2;
718
    int i, nb, opp_vertex, stack2_used;
719
720
721
722

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

Thomas Witkowski's avatar
Thomas Witkowski committed
723
    int sav_index, sav_neighbour = neighbour;
724

725
    TEST_EXIT_DBG(stack_used > 0)("no current element");
Thomas Witkowski's avatar
Thomas Witkowski committed
726
727
    TEST_EXIT_DBG(traverse_fill_flag.isSet(Mesh::CALL_LEAF_EL))
      ("invalid traverse_fill_flag");
728
729
730
731
732
733

    Parametric *parametric = traverse_mesh->getParametric();

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

734
    TEST_EXIT_DBG(elinfo_old == elinfo_stack[stack_used])("invalid old elinfo");
735
736

    elinfo_stack[stack_used]->testFlag(Mesh::FILL_NEIGH);
737
    el = elinfo_stack[stack_used]->getElement();
738
    sav_index = el->getIndex();
739
    sav_el = el;
740
741
742
743
744
745
746
747
748
749
750
751
752

    /* first, goto to leaf level, if necessary... */
    if (!(el->isLeaf()) && (neighbour < 2)) {
      if (stack_used >= static_cast<int>( elinfo_stack.size())-1) enlargeTraverseStack();
      i = 1 - neighbour;
      elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
      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
753
754
    save_stack_used = stack_used;
    for (i = 0; i<=stack_used; i++)
755
      save_info_stack[i] = info_stack[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
756
    for (i = 0; i<=stack_used; i++)
757
758
759
760
761
762
763
764
765
766
767
768
769
770
      (*(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  */
    /* element of the OLD hierarchy branch to the NEW branch                    */
    /****************************************************************************/

    nb = neighbour;

771
772
773
774
775
776
777
778
    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);
    }
779
780
781
782
783
784
785
786
787
788
789

    /****************************************************************************/
    /* 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 */
790
      } else {
791
792
793
	stack2_used = 1;
      }
      elinfo2 = save_elinfo_stack[stack2_used];
794
      el2 = elinfo2->getElement();
795
796
797

      i = traverse_mel->getOppVertex(nb);
      traverse_mel = traverse_mel->getNeighbour(nb);
798
799
      if (traverse_mel == NULL)
	return NULL;
800
801
802
803
804
805
      nb = i;
    
      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(const_cast<MacroElement*>( traverse_mel));
      info_stack[stack_used] = 0;
    
806
    } else {                                               /* goto other child */
807
808
809
810
811
812

      stack2_used = stack_used + 1;
      if (save_stack_used > stack2_used) {
	stack2_used++;               /* go down one level in OLD hierarchy */
      }
      elinfo2 = save_elinfo_stack[stack2_used];
813
      el2 = elinfo2->getElement();
814
815


816
      if (stack_used >= stack_size - 1) {
817
	enlargeTraverseStack();
818
      }
819
      i = 2 - info_stack[stack_used];
820
821
      info_stack[stack_used] = i + 1;
      elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
822
      stack_used++;
823
      nb = 1 - i;
824
825
826
827
828
829
830
831
832
    }

    /****************************************************************************/
    /* 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       */
    /* NEW hierarchy branch to the OLD branch.                                  */
    /****************************************************************************/

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

835
    while (el->getFirstChild()) {
836
837

      if (nb < 2) {   /* go down one level in hierarchy */
838
	if (stack_used >= stack_size - 1)
839
	  enlargeTraverseStack();
840
841
	elinfo_stack[stack_used + 1]->fillElInfo(1 - nb, elinfo_stack[stack_used]);
	info_stack[stack_used] = 2 - nb;
842
843
844
845
846
	stack_used++;
	nb = 2;
      }

      if (save_stack_used > stack2_used) { /* `refine' both el and el2 */
847
	TEST_EXIT_DBG(el->getFirstChild())("invalid NEW refinement?");
848
849
850

	/* use child i, neighbour of el2->child[nb-1] */
	i = 2 - save_info_stack[stack2_used];
851
	TEST_EXIT_DBG(i < 2)("invalid OLD refinement?");
852
853
	info_stack[stack_used] = i + 1;
	elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
854
855
856
857
	stack_used++;
	nb = i;

	elinfo = elinfo_stack[stack_used];
858
	el = elinfo->getElement();
859
860
861
862
863
864

	stack2_used++;
	if (save_stack_used > stack2_used) {
	  stack2_used++;                /* go down one level in OLD hierarchy */
	}
	elinfo2 = save_elinfo_stack[stack2_used];
865
	el2 = elinfo2->getElement();
866

867
      } else {   /* now we're done... */
868
	elinfo = elinfo_stack[stack_used];
869
	el = elinfo->getElement();
870
871
872
873
874
875
876
877
878
879
880
      }
    }


    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,
881
	  elinfo->getNeighbour(opp_vertex)?elinfo->getNeighbour(opp_vertex)->getIndex() : -1);
882
      TEST_EXIT_DBG(elinfo->getNeighbour(opp_vertex) == old_elinfo->getElement())
883
884
	("didn't succeed !?!?!?");
    }
885

886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
    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);
901
902
903
904

      elinfo->fillDetGrdLambda();
    }

905
    return elinfo;
906
907
908
909
  }

  void TraverseStack::update()
  {
910
    FUNCNAME("TraverseStack::update()");
911

912
    TEST_EXIT_DBG(traverse_mesh->getDim() == 3)("update only in 3d\n");
913

914
    for (int i = stack_used; i > 0; i--) {
915
      dynamic_cast<ElInfo3d*>(elinfo_stack[i])->update();
916
917
918
    }
  }

919
  void TraverseStack::getCoordsInElem(const ElInfo *upperElInfo, 
920
				      DimMat<double> *coords)
921
  {
922
923
    int levelDif = elinfo_stack[stack_used]->getLevel() - upperElInfo->getLevel();

924
925
    upperElInfo->getRefSimplexCoords(coords);
    
926
927
    for (int i = 1; i <= levelDif; i++) {
      upperElInfo->getSubElementCoords(coords,
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
				       elinfo_stack[stack_used - levelDif + i]->
				         getIChild());
    }
  }


  /****************************************************************************/
  /* common (static) traversal routines for 2d and 3d                         */
  /* to be #included in traverse_r.c                                          */
  /****************************************************************************/

  /* traverse hierarchical mesh,  call el_fct() for each/some element */

  /****************************************************************************/
  /*   recursive_traverse:                                		    */
  /*   -------------------                                		    */
  /*   recursively traverse the mesh hierarchy tree       		    */
  /*   call the routine traverse_info->el_fct(el_info) with    		    */
  /*    - all tree leaves                                 		    */
  /*    - all tree leaves at a special level              		    */
  /*    - all tree elements in pre-/in-/post-order        		    */
  /*   depending on the traverse_info->level variable     		    */
  /****************************************************************************/
  int Traverse::recursive(ElInfoStack *elInfoStack)
  {
    FUNCNAME("Traverse::recursive()");

    ElInfo *elinfo = elInfoStack->getCurrentElement();

    Element *el = elinfo->getElement();
    int mg_level, sum = 0;
    Parametric *parametric = mesh->getParametric();

    if (flag.isSet(Mesh::CALL_LEAF_EL)) {
      if (el->getFirstChild()) {
	ElInfo* elinfo_new = elInfoStack->getNextElement();
	elinfo_new->fillElInfo(0, elinfo);
	sum += recursive(elInfoStack);
	elinfo_new->fillElInfo(1, elinfo);
	sum += recursive(elInfoStack);
	elInfoStack->getBackElement();
      } else {
	if (el_fct != NULL) {
	  if (parametric) 
	    elinfo = parametric->addParametricInfo(elinfo);
	  elinfo->fillDetGrdLambda();
	  sum += el_fct(elinfo);
	  if (parametric) 
	    elinfo = parametric->removeParametricInfo(elinfo);
	}
      }
      return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0;
    }

    if (flag.isSet(Mesh::CALL_LEAF_EL_LEVEL)) {
      if (el->getFirstChild()) {
	if ((elinfo->getLevel() >= level)) {
	  return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0;
	}
	ElInfo* elinfo_new = elInfoStack->getNextElement();
	elinfo_new->fillElInfo(0, elinfo);
	sum += recursive(elInfoStack);
	elinfo->fillElInfo(1, elinfo);
	sum += recursive(elInfoStack);
	elInfoStack->getBackElement();
      } else {
	if ((elinfo->getLevel() == level) && (el_fct!=NULL)) {
	  if (parametric) 
	    elinfo = parametric->addParametricInfo(elinfo);
	  elinfo->fillDetGrdLambda();
	  sum += el_fct(elinfo);
	  if (parametric) 
	    elinfo = parametric->removeParametricInfo(elinfo);
	}
      }
      return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0;
    }

    if (flag.isSet(Mesh::CALL_EL_LEVEL)) {
      if (elinfo->getLevel() == level) {
	if (NULL != el_fct) {
	  if (parametric) 
	    elinfo = parametric->addParametricInfo(elinfo);
	  elinfo->fillDetGrdLambda();
	  sum += el_fct(elinfo); 
	  if (parametric) 
	    elinfo = parametric->removeParametricInfo(elinfo);
	}
      } else {
	if (elinfo->getLevel() > level){
	  return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0;
	} else if (el->getFirstChild()) {
	  ElInfo* elinfo_new = elInfoStack->getNextElement();
	  elinfo_new->fillElInfo(0, elinfo);
	  sum += recursive(elInfoStack);
	  elinfo_new->fillElInfo(1, elinfo);
	  sum += recursive(elInfoStack);
	  elInfoStack->getBackElement();
	}
      }

      return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0;
    }

    if (flag.isSet(Mesh::CALL_MG_LEVEL)) {

      mg_level = (elinfo->getLevel() + mesh->getDim() - 1) / mesh->getDim();

      if (mg_level > level)  {
	return 0;
      }

      if (!(el->getFirstChild())) {
	if (el_fct != NULL) {
	  if (parametric) 
	    elinfo = parametric->addParametricInfo(elinfo);
	  elinfo->fillDetGrdLambda();
	  sum += el_fct(elinfo);
	  if (parametric) 
	    elinfo = parametric->removeParametricInfo(elinfo);
	}
	return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0;
      }
    
      if ((mg_level == level) && ((elinfo->getLevel() % mesh->getDim()) == 0)) {
	if (el_fct != NULL) {
	  if (parametric) 
	    elinfo = parametric->addParametricInfo(elinfo);
	  elinfo->fillDetGrdLambda();
	  sum += el_fct(elinfo); 
	  if (parametric) 
	    elinfo = parametric->removeParametricInfo(elinfo);
	}
	return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0;
      }

      ElInfo* elinfo_new = elInfoStack->getNextElement();
 
      elinfo_new->fillElInfo(0, elinfo);
      sum += recursive(elInfoStack);

      elinfo_new->fillElInfo(1, elinfo);
      sum += recursive(elInfoStack);

      elInfoStack->getBackElement();

      return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0;   
    }

    if (flag.isSet(Mesh::CALL_EVERY_EL_PREORDER)) {
      if (el_fct != NULL) {
	if (parametric) 
	  elinfo = parametric->addParametricInfo(elinfo);
	elinfo->fillDetGrdLambda();
	sum += el_fct(elinfo);
	if (parametric) 
	  elinfo = parametric->removeParametricInfo(elinfo);
      }
1086
    }
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129

    if (el->getFirstChild()) {
      ElInfo* elinfo_new = elInfoStack->getNextElement();
   
      elinfo_new->fillElInfo(0, elinfo);
      sum += recursive(elInfoStack);

      if (flag.isSet(Mesh::CALL_EVERY_EL_INORDER)) 
	if (el_fct!=NULL) { 
	  if (parametric) 
	    elinfo = parametric->addParametricInfo(elinfo);
	  elinfo->fillDetGrdLambda();
	  sum += el_fct(elinfo);
	  if (parametric) 
	    elinfo = parametric->removeParametricInfo(elinfo);
	}
      elinfo_new->fillElInfo(1, elinfo);
      sum += recursive(elInfoStack);

      elInfoStack->getBackElement();
    } else {
      if (flag.isSet(Mesh::CALL_EVERY_EL_INORDER)) 
	if (el_fct != NULL) {
	  if (parametric) 
	    elinfo = parametric->addParametricInfo(elinfo);
	  elinfo->fillDetGrdLambda();
	  sum += el_fct(elinfo);
	  if (parametric) 
	    elinfo = parametric->removeParametricInfo(elinfo);
	}
    }

    if (flag.isSet(Mesh::CALL_EVERY_EL_POSTORDER)) 
      if (el_fct != NULL) {
	if (parametric) 
	  elinfo = parametric->addParametricInfo(elinfo);
	elinfo->fillDetGrdLambda();
	sum += el_fct(elinfo);
	if (parametric) 
	  elinfo = parametric->removeParametricInfo(elinfo);
      }

    return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0;
1130
  }
1131

1132
}