Traverse.cc 29.2 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
#include "Debug.h"
14
15
16

namespace AMDiS {

17
  ElInfo* TraverseStack::traverseFirst(Mesh *mesh, int level, Flag fill_flag)
18
  {
19
    FUNCNAME("TraverseStack::traverseFirst()");
20
21
22
23
    
    TEST_EXIT_DBG(fill_flag.isSet(Mesh::CALL_REVERSE_MODE) == false ||
		  fill_flag.isSet(Mesh::CALL_EVERY_EL_PREORDER))
      ("Not yet implemented!\n");
24

25
26
27
    traverse_mesh = mesh;
    traverse_level = level;
    traverse_fill_flag = fill_flag; 
28
29
30
31
32
33

    if (stack_size < 1)
      enlargeTraverseStack();

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

34
    for (int i = 0; i < stack_size; i++)
35
36
37
38
39
      elinfo_stack[i]->setFillFlag(fill_flag & FILL_ANY);

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

40
    if (fill_flag.isSet(Mesh::CALL_LEAF_EL_LEVEL)) {
41
      TEST_EXIT_DBG(level >= 0)("invalid level: %d\n", level);   
42
    }
43
44

    traverse_mel = NULL;
45
    stack_used = 0;
46

47
    return traverseNext(NULL);
48
49
  }

50

51
52
  ElInfo* TraverseStack::traverseNext(ElInfo* elinfo_old)
  {
53
    FUNCNAME("TraverseStack::traverseNext()");
54

55
    ElInfo *elinfo = NULL;
56
57
58
59
60
61
    Parametric *parametric = traverse_mesh->getParametric();

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

62
      TEST_EXIT_DBG(elinfo_old == elinfo_stack[stack_used])("invalid old elinfo\n");
63
    } else {
64
      TEST_EXIT_DBG(elinfo_old == NULL)("invalid old elinfo != nil\n");
65
66
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
67
    if (traverse_fill_flag.isSet(Mesh::CALL_LEAF_EL)) {
68
      elinfo = traverseLeafElement();
Thomas Witkowski's avatar
Thomas Witkowski committed
69
    } else {
70
71
72
73
74
75
76
      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
77
	if (traverse_fill_flag.isSet(Mesh::CALL_EVERY_EL_PREORDER)) {
78
	  elinfo = traverseEveryElementPreorder();
79
	} else if (traverse_fill_flag.isSet(Mesh::CALL_EVERY_EL_INORDER))
80
81
82
83
84
	  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
85
    }
86

87
    if (elinfo) {
88
      if (parametric)
89
90
91
92
	elinfo = parametric->addParametricInfo(elinfo);
      elinfo->fillDetGrdLambda();
    }

93
    return elinfo;
Thomas Witkowski's avatar
Thomas Witkowski committed
94
95
  }

96

97
98
  void TraverseStack::enlargeTraverseStack()
  {
99
100
101
    FUNCNAME("TraverseStack::enlargeTraverseStack()");

    int new_stack_size = stack_size + 10;
102
103

    elinfo_stack.resize(new_stack_size, NULL);
104
 
105
    // create new elinfos
106
107
    for (int i = stack_size; i < new_stack_size; i++) {
      TEST_EXIT_DBG(elinfo_stack[i] == NULL)("???\n");
108
109
110
      elinfo_stack[i] = traverse_mesh->createNewElInfo();
    }

111
112
    if (stack_size > 0)
      for (int i = stack_size; i < new_stack_size; i++)
113
	elinfo_stack[i]->setFillFlag(elinfo_stack[0]->getFillFlag());
114
115
116
117
118

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

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

    stack_size = new_stack_size;    
  }

  
  ElInfo* TraverseStack::traverseLeafElement()
  {
131
    FUNCNAME("TraverseStack::traverseLeafElement()");
132

133
    Element *el = NULL;
134

135
136
    if (stack_used == 0) {   /* first call */
      currentMacro = traverse_mesh->firstMacroElement();
137

138
      while (((*currentMacro)->getIndex() % maxThreads != myThreadId) &&
139
140
      	     currentMacro != traverse_mesh->endOfMacroElements())
      	currentMacro++;      
141

142
      if (currentMacro == traverse_mesh->endOfMacroElements())
143
	return NULL;
144

145
146
147
148
      traverse_mel = *currentMacro;
      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
      info_stack[stack_used] = 0;
149

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

170
	if (currentMacro == traverse_mesh->endOfMacroElements())
171
	  return NULL;
172
173

	traverse_mel = *currentMacro;	
174
175
	stack_used = 1;
	elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
176
	info_stack[stack_used] = 0;	
177
	el = elinfo_stack[stack_used]->getElement();
178
179

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

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

202
203
204
    return elinfo_stack[stack_used];
  }

205

206
207
  ElInfo* TraverseStack::traverseLeafElementLevel()
  {
208
    FUNCNAME("TraverseStack::traverseLeafElementLevel()");
209

210
    ERROR_EXIT("not yet implemented\n");
211

212
213
214
    return NULL;
  }

215

216
217
  ElInfo* TraverseStack::traverseElementLevel()
  {
218
    FUNCNAME("TraverseStack::traverseElementLevel()");
219
220
221
222
223
224

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

225
    return elInfo;
226
227
  }

228

229
230
  ElInfo* TraverseStack::traverseMultiGridLevel()
  {
231
    FUNCNAME("TraverseStack::traverseMultiGridLevel()");
232

233
234
235
    if (stack_used == 0) {   /* first call */
      currentMacro = traverse_mesh->firstMacroElement();
      traverse_mel = *currentMacro;
236
237
      if (traverse_mel == NULL)  
	return NULL;
238
239
240
241
242
243
244
245
      
      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()))
246
	return elinfo_stack[stack_used];
247
    }
248
  
249
    Element *el = elinfo_stack[stack_used]->getElement();
250
251

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


    /* goto next macro element */
    if (stack_used < 1) {
      currentMacro++;
262
      if (currentMacro == traverse_mesh->endOfMacroElements()) 
263
	return NULL;
264

265
      traverse_mel = *currentMacro;
266
267
268
269
      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
      info_stack[stack_used] = 0;

270
271
272
      if ((elinfo_stack[stack_used]->getLevel() == traverse_level) ||
	  (elinfo_stack[stack_used]->getLevel() < traverse_level && 
	   elinfo_stack[stack_used]->getElement()->isLeaf()))
273
	return elinfo_stack[stack_used];
274
275
276
277
278
    }


    /* go down tree */

279
280
281
282
    if (stack_used >= stack_size - 1) 
      enlargeTraverseStack();

    int i = info_stack[stack_used];
283
284
    info_stack[stack_used]++;

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

288
    TEST_EXIT_DBG(stack_used < stack_size)
289
290
291
292
293
      ("stack_size=%d too small, level=%d\n",
       stack_size, elinfo_stack[stack_used]->getLevel());

    info_stack[stack_used] = 0;
  
294
295
296
    if ((elinfo_stack[stack_used]->getLevel() == traverse_level) ||
	(elinfo_stack[stack_used]->getLevel() < traverse_level && 
	 elinfo_stack[stack_used]->getElement()->isLeaf()))
297
      return elinfo_stack[stack_used];
298
299
300
301

    return traverseMultiGridLevel();
  }

302

303
304
  ElInfo* TraverseStack::traverseEveryElementPreorder()
  {
305
    FUNCNAME("TraverseStack::traverseEveryElementPreorder()");
306

307
308
309
310
311
312
313
314
315
    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;
316

317
      return elinfo_stack[stack_used];
318
    }
319
  
320
    Element *el = elinfo_stack[stack_used]->getElement();
321
322

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


    /* goto next macro element */
    if (stack_used < 1) {
      currentMacro++;
333
      if (currentMacro == traverse_mesh->endOfMacroElements()) 
334
	return NULL;
335
336
337
338
339
340
      traverse_mel = *currentMacro;

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

341
      return elinfo_stack[stack_used];
342
343
344
345
346
    }


    /* go down tree */

347
    if (stack_used >= stack_size - 1) 
348
349
      enlargeTraverseStack();

350
    int fillIthChild = info_stack[stack_used];
351

352
    info_stack[stack_used]++;
353
354
    if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
      fillIthChild = 1 - fillIthChild;
355

356
    elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);
357

358
359
    stack_used++;

360
    TEST_EXIT_DBG(stack_used < stack_size)
361
      ("stack_size = %d too small, level = %d\n",
362
363
364
365
       stack_size, elinfo_stack[stack_used]->getLevel());

    info_stack[stack_used] = 0;
  
366
    return elinfo_stack[stack_used];
367
368
  }

369

370
371
372
373
374
375
376
  ElInfo* TraverseStack::traverseEveryElementInorder()
  {
    FUNCNAME("TraverseStack::traverseEveryElementInorder");
    ERROR_EXIT("not yet implemented\n");
    return NULL;
  }

377

378
379
  ElInfo* TraverseStack::traverseEveryElementPostorder()
  {
380
    FUNCNAME("TraverseStack::traverseEveryElementPostorder()");
381

382
383
384
385
386
387
388
389
390
391
392
393
    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 */
394
      Element *el = elinfo_stack[stack_used]->getElement();
395
396

      /* go up in tree until we can go down again */          /* postorder!!! */
397
398
      while (stack_used > 0 && 
	     (info_stack[stack_used] >= 3 || el->getFirstChild() == NULL)) {
399
400
401
	stack_used--;
	el = elinfo_stack[stack_used]->getElement();
      }
402
403
404
405
406


      /* goto next macro element */
      if (stack_used < 1) {
	currentMacro++;
407
408
	if (currentMacro == traverse_mesh->endOfMacroElements()) 
	  return NULL;
409
410
411
412
413
414
415
416
417
418
419
	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 */

420
421
    while (elinfo_stack[stack_used]->getElement()->getFirstChild() &&
	   info_stack[stack_used] < 2) {
422
      if (stack_used >= stack_size-1) 
423
424
	enlargeTraverseStack();

425
      int i = info_stack[stack_used];
426
      info_stack[stack_used]++;
427
      elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
428
429
430
431
432
433
      stack_used++;
      info_stack[stack_used] = 0;
    }
  
    info_stack[stack_used]++;      /* postorder!!! */

434
    return elinfo_stack[stack_used];
435
436
  }

437

438
439
440
441
442
443
444
445
446
447
448
449
450
  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;
451
452
    default:
      ERROR_EXIT("invalid dim\n");
453
454
455
456
    }
    return NULL;
  }

457

458
459
  ElInfo* TraverseStack::traverseNeighbour3d(ElInfo* elinfo_old, int neighbour)
  {
460
    FUNCNAME("TraverseStack::traverseNeighbour3d()");
Thomas Witkowski's avatar
Thomas Witkowski committed
461

462
    Element *el2;
463
464
465
    ElInfo *elinfo2;
    int stack2_used;
    int sav_neighbour = neighbour;
466

467
468
469
470
    // father.neigh[coarse_nb[i][j]] == child[i - 1].neigh[j]
    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}}};
471

472
    TEST_EXIT_DBG(stack_used > 0)("no current element\n");
473
474

    Parametric *parametric = traverse_mesh->getParametric();
475
476
    if (parametric) 
      elinfo_old = parametric->removeParametricInfo(elinfo_old);
477

478
    TEST_EXIT_DBG(elinfo_old == elinfo_stack[stack_used])("invalid old elinfo\n");
479
    TEST_EXIT_DBG(elinfo_stack[stack_used]->getFillFlag().isSet(Mesh::FILL_NEIGH))
480
      ("FILL_NEIGH not set");
481

482
    Element *el = elinfo_stack[stack_used]->getElement();
483
    int sav_index = el->getIndex();
484

485
    // First, goto to leaf level, if necessary ...
486
    if ((traverse_fill_flag & Mesh::CALL_LEAF_EL).isAnySet()) {
487
      if (el->getChild(0) && neighbour < 2) {
488
	if (stack_used >= stack_size - 1)
489
	  enlargeTraverseStack();
490
491
	int i = 1 - neighbour;

492
	elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
Thomas Witkowski's avatar
Thomas Witkowski committed
493
494
495
 	if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
 	  info_stack[stack_used] = (i == 0 ? 2 : 1);
 	else
496
	  info_stack[stack_used] = i + 1;
497
498
499
500
501
502
503
504
	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;
505
    save_stack_used = stack_used;
506

507
508
509
    // === First phase (see 2D). ===

    int nb = neighbour;
510

511
512
    while (stack_used > 1) { /* go up in tree until we can go down again */
      stack_used--;
513
514
515
516
517
518
519
520
521
522
523
      int typ = elinfo_stack[stack_used]->getType();
      int elIsIthChild = info_stack[stack_used];
      if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE) && elIsIthChild != 0) 
	elIsIthChild = (elIsIthChild == 1 ? 2 : 1);

      TEST_EXIT_DBG(!elinfo_stack[stack_used + 1]->getParent() ||
	  elinfo_stack[stack_used + 1]->getParent()->getChild(elIsIthChild - 1) == 
	  elinfo_stack[stack_used + 1]->getElement())
	("Should not happen!\n");

      nb = coarse_nb[typ][elIsIthChild][nb];
524

525
526
      if (nb == -1) 
	break;
527

528
      TEST_EXIT_DBG(nb >= 0)("Invalid coarse_nb %d!\n", nb);
529
    }
530

Thomas Witkowski's avatar
Thomas Witkowski committed
531
532
533
534
535
536
537
    for (int i = stack_used; i <= save_stack_used; i++) {
      save_info_stack[i] = info_stack[i];
      *(save_elinfo_stack[i]) = *(elinfo_stack[i]);
    }
    ElInfo *old_elinfo = save_elinfo_stack[save_stack_used];
    int opp_vertex = old_elinfo->getOppVertex(neighbour);

538

539
540
    if (nb >= 0) {                           
      // Go to macro element neighbour.
541

542
      int i = traverse_mel->getOppVertex(nb);
543

544
      traverse_mel = traverse_mel->getNeighbour(nb);
545
      if (traverse_mel == NULL)  
546
	return NULL;
547
    
548
      if (nb < 2 && save_stack_used > 1)
549
	stack2_used = 2;                /* go down one level in OLD hierarchy */
550
      else
551
	stack2_used = 1;
552
      
553
554
555
556
557
558
      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;
559
560
561
    } else {                                                
      // Goto other child.

562
      stack2_used = stack_used + 1;
563
      if (save_stack_used > stack2_used)
564
	stack2_used++;                  /* go down one level in OLD hierarchy */
565

566
567
568
      elinfo2 = save_elinfo_stack[stack2_used];
      el2 = elinfo2->getElement();

569
      if (stack_used >= stack_size - 1)
570
	enlargeTraverseStack();
571
572
573
574
575
576

      int i = 2 - info_stack[stack_used];
      info_stack[stack_used] = i + 1;
      int fillIthChild = i;
      if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
	fillIthChild = 1 - fillIthChild;
577
      elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);
578
579
580
581
582
      stack_used++;
      info_stack[stack_used] = 0;
      nb = 0;
    }

583
584
585
586

    // === Second phase. ===

    ElInfo *elinfo = elinfo_stack[stack_used];
587
588
    el = elinfo->getElement();

589
590
591
592
593
    while (el->getChild(0)) {
      if (nb < 2) {                         
	// Go down one level in hierarchy.

	if (stack_used >= stack_size - 1)
594
	  enlargeTraverseStack();
595
596
597
598
599
600
601
602


	if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE)) {
	  int t = 2 - nb;
	  info_stack[stack_used] = (t == 2 ? 1 : 2);
	} else {	
	  info_stack[stack_used] = 2 - nb;	
	}
603
604
605
606

	int fillIthChild = 1 - nb;
	elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);

607
608
609
610
611
612
613
	stack_used++;
	info_stack[stack_used] = 0;
	elinfo = elinfo_stack[stack_used];
	el = elinfo->getElement();
	nb = 3;
      }

614
615
616
      if (save_stack_used > stack2_used) { 
	// `refine' both el and el2.

Thomas Witkowski's avatar
Thomas Witkowski committed
617
	TEST_EXIT_DBG(el->getChild(0))("invalid new refinement?\n");
618

619
	int i = 0;
620
621
622
623
624
625
626
627
628
629
	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 {
630
	    ERROR_EXIT("No common refinement edge! %d\n", traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE));
631
632
633
	  }
	}

634
635
636
637
	int testChild = i;
 	if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
 	  testChild = 1 - testChild;

638
	if (el->getChild(0) &&  
639
640
	    (el->getChild(testChild)->getDOF(1) == el->getDOF(nb) ||
	     traverse_mesh->associated(el->getChild(testChild)->getDOF(1, 0), el->getDOF(nb, 0))))
641
642
643
	  nb = 1;	
	else
	  nb = 2;	
644

645
646
	info_stack[stack_used] = i + 1;	

647
	if (stack_used >= stack_size - 1)
648
	  enlargeTraverseStack();
649
650
651
652
653
654

	int fillIthChild = i;
	if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
	  fillIthChild = 1 - fillIthChild;
	elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);

655
656
657
658
659
660
661
662
663
	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();
664

665
	if (save_stack_used > stack2_used) {
666
	  const DegreeOfFreedom *dof = el2->getDOF(1);
667
668
	  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
      } else {   
	// Now we're done...
678
679
680

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

682
683
684
685
686
687
688
	break;
      }
    }


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

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

    return elinfo;
712
713
  }

714

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

719
720
721
722
    Triangle *el2;
    ElInfo *elinfo2;
    int stack2_used;
    int sav_neighbour = neighbour;
723

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

727
    TEST_EXIT_DBG(stack_used > 0)("no current element");
728
729
730
731
732

    Parametric *parametric = traverse_mesh->getParametric();
    if (parametric) 
      elinfo_old = parametric->removeParametricInfo(elinfo_old);

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

    elinfo_stack[stack_used]->testFlag(Mesh::FILL_NEIGH);
736
737
738
739
    Triangle *el = 
      dynamic_cast<Triangle*>(const_cast<Element*>(elinfo_stack[stack_used]->getElement()));
    int sav_index = el->getIndex();
    Triangle *sav_el = el;
740
741

    /* first, goto to leaf level, if necessary... */
742
    if (!(el->isLeaf()) && neighbour < 2) {
743

744
745
746
747
748
749
750
751
752
753
754
755
756
      if (stack_used >= static_cast<int>(elinfo_stack.size()) - 1) 
	enlargeTraverseStack();

      // If we should search for neighbour 0, take second child, if for
      // neighbour 1, take the first child.
      int i = 1 - neighbour;

      elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
      if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
	info_stack[stack_used] = (i == 0 ? 2 : 1);
      else
	info_stack[stack_used] = i + 1;
      
757
      stack_used++;
758
      neighbour = 2;     
759
760
761
762
    }

    /* save information about current element and its position in the tree */
    save_traverse_mel = traverse_mel;
763
764
    save_stack_used = stack_used;
    for (int i = 0; i <= stack_used; i++) {
765
766
      save_info_stack[i] = info_stack[i];
      (*(save_elinfo_stack[i])) = (*(elinfo_stack[i]));
767
768
769
    }
    ElInfo *old_elinfo = save_elinfo_stack[stack_used];
    int opp_vertex = old_elinfo->getOppVertex(neighbour);
770
771
772
773
774
775


    /****************************************************************************/
    /* 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
776
    /* element of the OLD hierarchy branch to the new branch                    */
777
778
    /****************************************************************************/

779
    int nb = neighbour;
780

781
782
    while (stack_used > 1) {
      stack_used--;
783
      int elIsIthChild = info_stack[stack_used];
784

785
786
787
788
789
790
791
792
793
      if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE) && elIsIthChild != 0) 
	elIsIthChild = (elIsIthChild == 1 ? 2 : 1);

      TEST_EXIT_DBG(!elinfo_stack[stack_used + 1]->getParent() ||
	  elinfo_stack[stack_used + 1]->getParent()->getChild(elIsIthChild - 1) == 
	  elinfo_stack[stack_used + 1]->getElement())
	("Should not happen!\n");

      nb = coarse_nb[elIsIthChild][nb];
794
      if (nb == -1) 
795
796
	break;     

797
798
      TEST_EXIT_DBG(nb >= 0)("invalid coarse_nb %d\n",nb);
    }
799

800

801
802
803
804
805
806
    /****************************************************************************/
    /* 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)                               */
    /****************************************************************************/

807
808
    if (nb >= 0) {                        
      // Go to macro element neighbour.
809

810
      if (nb < 2 && save_stack_used > 1)
811
	stack2_used = 2;           /* go down one level in OLD hierarchy */
812
      else
813
	stack2_used = 1;
814
      
815
      elinfo2 = save_elinfo_stack[stack2_used];
816
      el2 = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo2->getElement()));
817

818
      int i = traverse_mel->getOppVertex(nb);
819
      traverse_mel = traverse_mel->getNeighbour(nb);
820
821
      if (traverse_mel == NULL)
	return NULL;
822
      nb = i;
823

824
      stack_used = 1;
825
      elinfo_stack[stack_used]->fillMacroInfo(const_cast<MacroElement*>(traverse_mel));
826
      info_stack[stack_used] = 0;
827
828
    } else {                                               
      // Go to other child.
829
830

      stack2_used = stack_used + 1;
831
      if (save_stack_used > stack2_used)
832
	stack2_used++;               /* go down one level in OLD hierarchy */
833
      
834
      elinfo2 = save_elinfo_stack[stack2_used];
835
      el2 = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo2->getElement()));
836

837
      if (stack_used >= stack_size - 1)
838
	enlargeTraverseStack();
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854

      TEST_EXIT_DBG(info_stack[stack_used] == 1 || info_stack[stack_used] == 2)
	("Should not happen!\n");

      int fillIthChild = -1;
      if (info_stack[stack_used] == 1) {
	info_stack[stack_used] = 2;
	fillIthChild = 1;
	nb = 0;
      } else {
	info_stack[stack_used] = 1;
	fillIthChild = 0;
	nb = 1;
      }

      if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE)) {
855
	fillIthChild = 1 - fillIthChild;
856
857
	nb = 1 - nb;
      }
858
      elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);
859
860
861
862
863
864
      stack_used++;
    }

    /****************************************************************************/
    /* 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
865
    /* new hierarchy branch to the OLD branch.                                  */
866
867
    /****************************************************************************/

868
869
    ElInfo *elinfo = elinfo_stack[stack_used];
    el = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo->getElement()));
870

871
    while (el->getFirstChild()) {
872
873
      if (nb < 2) {   
	// Go down one level in hierarchy.
874

875
	if (stack_used >= stack_size - 1)
876
	  enlargeTraverseStack();
877

878
879
880
881
882
883
884
	if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE)) {
	  int t = 2 - nb;
	  info_stack[stack_used] = (t == 2 ? 1 : 2);
	} else {	
	  info_stack[stack_used] = 2 - nb;	
	}

885
886
887
	int fillIthChild = 1 - nb;
	elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);

888
889
890
891
	stack_used++;
	nb = 2;