Traverse.cc 32.7 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
36
    if (fill_flag.isSet(Mesh::CALL_LEAF_EL_LEVEL))
      TEST_EXIT_DBG(level >= 0)("invalid level: %d\n", level);   
37
38

    traverse_mel = NULL;
39
    stack_used = 0;
40

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

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

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

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

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

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

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

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


Thomas Witkowski's avatar
Thomas Witkowski committed
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
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
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290

  /****************************************************************************/
  /* 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);
      }
    }

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

291
292
  void TraverseStack::enlargeTraverseStack()
  {
293
294
295
    FUNCNAME("TraverseStack::enlargeTraverseStack()");

    int new_stack_size = stack_size + 10;
296
297

    elinfo_stack.resize(new_stack_size, NULL);
298
 
299
    // create new elinfos
300
301
    for (int i = stack_size; i < new_stack_size; i++) {
      TEST_EXIT_DBG(elinfo_stack[i] == NULL)("???\n");
302
303
304
      elinfo_stack[i] = traverse_mesh->createNewElInfo();
    }

305
306
307
308
309
    if (stack_size > 0) {
      for (int i = stack_size; i < new_stack_size; i++) {
	elinfo_stack[i]->setFillFlag(elinfo_stack[0]->getFillFlag());
      }
    }
310
311
312
313
314

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

    // create new elinfos
315
316
    for (int i = stack_size; i < new_stack_size; i++) {
      TEST_EXIT_DBG(save_elinfo_stack[i] == NULL)("???\n");
317
318
319
320
321
322
323
324
325
326
      save_elinfo_stack[i] = traverse_mesh->createNewElInfo();
    }
    save_info_stack.resize(new_stack_size);

    stack_size = new_stack_size;    
  }

  
  ElInfo* TraverseStack::traverseLeafElement()
  {
327
    FUNCNAME("TraverseStack::traverseLeafElement()");
328

329
    Element *el = NULL;
330

331
332
    if (stack_used == 0) {   /* first call */
      currentMacro = traverse_mesh->firstMacroElement();
333

334
      while (((*currentMacro)->getIndex() % maxThreads != myThreadId) &&
Thomas Witkowski's avatar
Thomas Witkowski committed
335
      	     (currentMacro != traverse_mesh->endOfMacroElements())) {
336
337
338
      	currentMacro++;
      }

339
      if (currentMacro == traverse_mesh->endOfMacroElements())
340
	return NULL;
341

342
343
344
345
      traverse_mel = *currentMacro;
      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
      info_stack[stack_used] = 0;
346

347
      el = elinfo_stack[stack_used]->getElement();
348
349
      if (el == NULL || el->getFirstChild() == NULL)
	return (elinfo_stack[stack_used]);      
350
351
352
353
    } else {
      el = elinfo_stack[stack_used]->getElement();
      
      /* go up in tree until we can go down again */
354
355
356
357
358
      while ((stack_used > 0) && 
	     ((info_stack[stack_used] >= 2) || (el->getFirstChild() == NULL))) {
	stack_used--;
	el = elinfo_stack[stack_used]->getElement();
      }
359
360
361
      
      /* goto next macro element */
      if (stack_used < 1) {
362
363
	do {	
	  currentMacro++;
Thomas Witkowski's avatar
Thomas Witkowski committed
364
	} while ((currentMacro != traverse_mesh->endOfMacroElements()) && 
365
		 ((*currentMacro)->getIndex() % maxThreads != myThreadId));
366

367
	if (currentMacro == traverse_mesh->endOfMacroElements())
368
	  return NULL;
369
370

	traverse_mel = *currentMacro;	
371
372
	stack_used = 1;
	elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
373
	info_stack[stack_used] = 0;	
374
	el = elinfo_stack[stack_used]->getElement();
375
376
377

	if (el == NULL || el->getFirstChild() == NULL)
	  return (elinfo_stack[stack_used]);	
378
      }
379
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
380
   
381
    /* go down tree until leaf */
382
    while (el->getFirstChild()) {
383
      if (stack_used >= stack_size - 1) 
384
385
	enlargeTraverseStack();
      
Thomas Witkowski's avatar
Thomas Witkowski committed
386
387
      int i = info_stack[stack_used];
      el = const_cast<Element*>(((i == 0) ? el->getFirstChild() : el->getSecondChild()));
388
      info_stack[stack_used]++;
Thomas Witkowski's avatar
Thomas Witkowski committed
389
      elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
390
      stack_used++;
391

392
      TEST_EXIT_DBG(stack_used < stack_size)
393
	("stack_size=%d too small, level=(%d,%d)\n",
Thomas Witkowski's avatar
Thomas Witkowski committed
394
395
396
	 stack_size, elinfo_stack[stack_used]->getLevel());
      
      info_stack[stack_used] = 0;
397
    }
398

399
400
401
402
403
    return elinfo_stack[stack_used];
  }

  ElInfo* TraverseStack::traverseLeafElementLevel()
  {
404
    FUNCNAME("TraverseStack::traverseLeafElementLevel()");
405
    ERROR_EXIT("not yet implemented\n");
406

407
408
409
410
411
    return NULL;
  }

  ElInfo* TraverseStack::traverseElementLevel()
  {
412
    FUNCNAME("TraverseStack::traverseElementLevel()");
413
414
415
416
417
418

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

419
    return elInfo;
420
421
422
423
  }

  ElInfo* TraverseStack::traverseMultiGridLevel()
  {
424
    FUNCNAME("TraverseStack::traverseMultiGridLevel()");
425

426
427
428
429
430
431
432
433
434
435
436
437
438
439
    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]);
    }
440
  
441
    Element *el = elinfo_stack[stack_used]->getElement();
442
443

    /* go up in tree until we can go down again */
444
445
446
447
448
    while ((stack_used > 0) && 
	   ((info_stack[stack_used] >= 2) || (el->getFirstChild()==NULL))) {
      stack_used--;
      el = elinfo_stack[stack_used]->getElement();
    }
449
450
451
452
453


    /* goto next macro element */
    if (stack_used < 1) {
      currentMacro++;
454
455
      if (currentMacro == traverse_mesh->endOfMacroElements()) 
	return(NULL);
456

457
      traverse_mel = *currentMacro;
458
459
460
461
      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
      info_stack[stack_used] = 0;

462
463
464
      if ((elinfo_stack[stack_used]->getLevel() == traverse_level) ||
	  (elinfo_stack[stack_used]->getLevel() < traverse_level && 
	   elinfo_stack[stack_used]->getElement()->isLeaf()))
465
466
467
468
469
470
	return(elinfo_stack[stack_used]);
    }


    /* go down tree */

471
472
473
474
    if (stack_used >= stack_size - 1) 
      enlargeTraverseStack();

    int i = info_stack[stack_used];
475
476
477
478
479
    info_stack[stack_used]++;

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

480
    TEST_EXIT_DBG(stack_used < stack_size)
481
482
483
484
485
      ("stack_size=%d too small, level=%d\n",
       stack_size, elinfo_stack[stack_used]->getLevel());

    info_stack[stack_used] = 0;
  
486
487
488
    if ((elinfo_stack[stack_used]->getLevel() == traverse_level) ||
	(elinfo_stack[stack_used]->getLevel() < traverse_level && 
	 elinfo_stack[stack_used]->getElement()->isLeaf()))
489
490
491
492
493
494
495
      return(elinfo_stack[stack_used]);

    return traverseMultiGridLevel();
  }

  ElInfo* TraverseStack::traverseEveryElementPreorder()
  {
496
    FUNCNAME("TraverseStack::traverseEveryElementPreorder()");
497
498
499

    Element *el;

500
501
502
503
504
505
506
507
508
    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;
509

510
511
      return(elinfo_stack[stack_used]);
    }
512
513
514
515
  
    el = elinfo_stack[stack_used]->getElement();

    /* go up in tree until we can go down again */
516
517
518
519
520
    while (stack_used > 0 && 
	   (info_stack[stack_used] >= 2 || el->getFirstChild() == NULL)) {
      stack_used--;
      el = elinfo_stack[stack_used]->getElement();
    }
521
522
523
524
525


    /* goto next macro element */
    if (stack_used < 1) {
      currentMacro++;
526
527
      if (currentMacro == traverse_mesh->endOfMacroElements()) 
	return(NULL);
528
529
530
531
532
533
534
535
536
537
538
539
      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 */

540
    if (stack_used >= stack_size - 1) 
541
542
543
      enlargeTraverseStack();

    int i = info_stack[stack_used];
544
545
546
547
548
    info_stack[stack_used]++;

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

549
    TEST_EXIT_DBG(stack_used < stack_size)
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
      ("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;

572
573
574
575
576
577
578
579
580
581
582
583
    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 */
584
585
586
587
  
      el = elinfo_stack[stack_used]->getElement();

      /* go up in tree until we can go down again */          /* postorder!!! */
588
589
590
591
592
      while ((stack_used > 0) && 
	     ((info_stack[stack_used] >= 3) || (el->getFirstChild()==NULL))) {
	stack_used--;
	el = elinfo_stack[stack_used]->getElement();
      }
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609


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

610
611
    while (elinfo_stack[stack_used]->getElement()->getFirstChild()
	   && (info_stack[stack_used] < 2)) {
612
613
614
      if (stack_used >= stack_size-1) 
	enlargeTraverseStack();

615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
      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
647
648
    FUNCNAME("TraverseStackLLtraverseNeighbour3d()");

649
    Element *el2;
650
651
652
653
654
655
656
657
658
    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}}};

659
    TEST_EXIT_DBG(stack_used > 0)("no current element\n");
660
661
662

    Parametric *parametric = traverse_mesh->getParametric();

663
664
    if (parametric) 
      elinfo_old = parametric->removeParametricInfo(elinfo_old);
665

666
    TEST_EXIT_DBG(elinfo_old == elinfo_stack[stack_used])
667
668
      ("invalid old elinfo\n");

669
    TEST_EXIT_DBG(elinfo_stack[stack_used]->getFillFlag().isSet(Mesh::FILL_NEIGH))
670
      ("FILL_NEIGH not set");
671
    Element *el = elinfo_stack[stack_used]->getElement();
672
673
674
675
676
    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)) {
677
	if (stack_used >= stack_size - 1)
678
679
	  enlargeTraverseStack();
	i = 1 - neighbour;
680
681
	elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
	info_stack[stack_used] = i + 1;
682
683
684
685
686
687
688
689
690
	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;
691
    save_stack_used = stack_used;
692
693
694

    nb = neighbour;

695
696
697
698
699
700
701
702
    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);
    }
703
704
705

    /* save hierarchy information about current element */

706
707
    for (i = stack_used; i <= save_stack_used; i++) {
      save_info_stack[i] = info_stack[i];
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
      *(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
733
    } else {                                                /* goto other child */
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
      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
762
	elinfo_stack[stack_used+1]->fillElInfo(1 - nb, elinfo_stack[stack_used]);
763
764
765
766
767
768
769
770
	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
771
	TEST_EXIT_DBG(el->getChild(0))("invalid new refinement?\n");
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813

	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)) &&
814
815
816
817
818
819
	      !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();
	  } 
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
	}

      }
      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);
839
      TEST_EXIT_DBG(elinfo->getNeighbour(opp_vertex) == old_elinfo->getElement())
840
841
842
843
844
845
846
847
848
	("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;  /* ??? */

849
850
851
    if (elinfo) {
      if (parametric) 
	elinfo = parametric->addParametricInfo(elinfo);
852
853
854
855
856
857
858
      elinfo->fillDetGrdLambda();
    }
    return(elinfo);
  }

  ElInfo* TraverseStack::traverseNeighbour2d(ElInfo* elinfo_old, int neighbour)
  {
859
860
    FUNCNAME("TraverseStack::traverseNeighbour2d()");

Thomas Witkowski's avatar
Thomas Witkowski committed
861
    Triangle *el, *el2, *sav_el;
862
    ElInfo *old_elinfo, *elinfo, *elinfo2;
863
    int i, nb, opp_vertex, stack2_used;
864
865
866
867

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

868
    int sav_index, sav_neighbour = neighbour;
869

870
    TEST_EXIT_DBG(stack_used > 0)("no current element");
Thomas Witkowski's avatar
Thomas Witkowski committed
871
    TEST_EXIT_DBG(traverse_fill_flag.isSet(Mesh::CALL_LEAF_EL))("invalid traverse_fill_flag");
872
873
874
875
876
877

    Parametric *parametric = traverse_mesh->getParametric();

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

878
    TEST_EXIT_DBG(elinfo_old == elinfo_stack[stack_used])("invalid old elinfo");
879
880

    elinfo_stack[stack_used]->testFlag(Mesh::FILL_NEIGH);
Thomas Witkowski's avatar
Thomas Witkowski committed
881
    el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo_stack[stack_used]->getElement()));
882
    sav_index = el->getIndex();
883
    sav_el = el;
884
885
886
887
888
889
890
891
892
893
894
895
896

    /* 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
897
898
    save_stack_used   = stack_used;
    for (i=0; i<=stack_used; i++)
899
      save_info_stack[i] = info_stack[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
900
    for (i=0; i<=stack_used; i++)
901
902
903
904
905
906
907
908
909
      (*(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
910
    /* element of the OLD hierarchy branch to the new branch                    */
911
912
913
914
    /****************************************************************************/

    nb = neighbour;

915
916
917
918
919
920
921
922
    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);
    }
923
924
925
926
927
928
929
930
931
932
933

    /****************************************************************************/
    /* 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 */
934
      } else {
935
936
937
	stack2_used = 1;
      }
      elinfo2 = save_elinfo_stack[stack2_used];
Thomas Witkowski's avatar
Thomas Witkowski committed
938
      el2 = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo2->getElement()));
939
940
941

      i = traverse_mel->getOppVertex(nb);
      traverse_mel = traverse_mel->getNeighbour(nb);
942
943
      if (traverse_mel == NULL)
	return NULL;
944
945
946
947
948
949
      nb = i;
    
      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(const_cast<MacroElement*>( traverse_mel));
      info_stack[stack_used] = 0;
    
950
    } else {                                               /* goto other child */
951
952
953
954
955
956

      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];
Thomas Witkowski's avatar
Thomas Witkowski committed
957
      el2 = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo2->getElement()));
958
959


Thomas Witkowski's avatar
Thomas Witkowski committed
960
      if (stack_used >= stack_size-1) {
961
	enlargeTraverseStack();
962
      }
963
      i = 2 - info_stack[stack_used];
Thomas Witkowski's avatar
Thomas Witkowski committed
964
965
      info_stack[stack_used] = i+1;
      elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
966
      stack_used++;
Thomas Witkowski's avatar
Thomas Witkowski committed
967
      nb = 1-i;
968
969
970
971
972
    }

    /****************************************************************************/
    /* 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
973
    /* new hierarchy branch to the OLD branch.                                  */
974
975
976
    /****************************************************************************/

    elinfo = elinfo_stack[stack_used];
Thomas Witkowski's avatar
Thomas Witkowski committed
977
    el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo->getElement()));
978

979
    while (el->getFirstChild()) {
980
981

      if (nb < 2) {   /* go down one level in hierarchy */
Thomas Witkowski's avatar
Thomas Witkowski committed
982
	if (stack_used >= stack_size-1)
983
	  enlargeTraverseStack();
Thomas Witkowski's avatar
Thomas Witkowski committed
984
985
	elinfo_stack[stack_used+1]->fillElInfo(1-nb, elinfo_stack[stack_used]);
	info_stack[stack_used] = 2-nb;
986
987
988
989
990
	stack_used++;
	nb = 2;
      }

      if (save_stack_used > stack2_used) { /* `refine' both el and el2 */
Thomas Witkowski's avatar
Thomas Witkowski committed
991
	TEST_EXIT_DBG(el->getFirstChild())("invalid new refinement?");
992
993
994

	/* use child i, neighbour of el2->child[nb-1] */
	i = 2 - save_info_stack[stack2_used];
995
	TEST_EXIT_DBG(i < 2)("invalid OLD refinement?");
Thomas Witkowski's avatar
Thomas Witkowski committed
996
997
	info_stack[stack_used] = i+1;
	elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
998
999
1000
1001
	stack_used++;
	nb = i;

	elinfo = elinfo_stack[stack_used];
Thomas Witkowski's avatar
Thomas Witkowski committed
1002
	el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo->getElement()));
1003
1004
1005
1006
1007
1008

	stack2_used++;
	if (save_stack_used > stack2_used) {
	  stack2_used++;                /* go down one level in OLD hierarchy */
	}
	elinfo2 = save_elinfo_stack[stack2_used];
Thomas Witkowski's avatar
Thomas Witkowski committed
1009
	el2 = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo2->getElement()));
1010

Thomas Witkowski's avatar
Thomas Witkowski committed
1011
1012
      }
      else {   /* now we're done... */
1013
	elinfo = elinfo_stack[stack_used];
Thomas Witkowski's avatar
Thomas Witkowski committed
1014
	el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo->getElement()));
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
      }
    }


    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
1026
	  elinfo->getNeighbour(opp_vertex)?elinfo->getNeighbour(opp_vertex)->getIndex():-1);
1027
      TEST_EXIT_DBG(elinfo->getNeighbour(opp_vertex) == old_elinfo->getElement())
1028
1029
	("didn't succeed !?!?!?");
    }
1030

1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
    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);
1046
1047
1048
1049

      elinfo->fillDetGrdLambda();
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1050
    return(elinfo);
1051
1052
1053
1054
  }

  void TraverseStack::update()
  {
1055
    FUNCNAME("TraverseStack::update()");
1056

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

1059
    for (int i = stack_used; i > 0; i--)
1060
      dynamic_cast<ElInfo3d*>(elinfo_stack[i])->update();
1061
1062
  }

1063
  void TraverseStack::fillRefinementPath(ElInfo& elInfo, const ElInfo& upperElInfo)
1064
  {
1065
1066
    int levelDif = elinfo_stack[stack_used]->getLevel() - upperElInfo.getLevel();
    unsigned long rPath = 0;
1067

1068
1069
1070
1071
1072
1073
    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);
1074
1075
  }

1076
}