Traverse.cc 31.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
14
15
16
17
18
19
20
21
22
23
#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"
24
#include "Debug.h"
25
26
27

namespace AMDiS {

28
  ElInfo* TraverseStack::traverseFirst(Mesh *mesh, int level, Flag fill_flag)
29
  {
30
    FUNCNAME("TraverseStack::traverseFirst()");
31

32
33
34
    TEST_EXIT_DBG(fill_flag.isSet(Mesh::CALL_REVERSE_MODE) == false ||
		  fill_flag.isSet(Mesh::CALL_EVERY_EL_PREORDER))
      ("Not yet implemented!\n");
35

36
37
38
    traverse_mesh = mesh;
    traverse_level = level;
    traverse_fill_flag = fill_flag; 
39

40
    TEST_EXIT_DBG(mesh)("No mesh!\n");
41
    TEST_EXIT(traverse_mesh->getMacroElements().size() > 0)("Mesh is empty!\n");
42

43
44
45
46
47
    if (stack_size < 1)
      enlargeTraverseStack();

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

48
    for (int i = 0; i < stack_size; i++)
49
50
51
52
53
      elinfo_stack[i]->setFillFlag(fill_flag & FILL_ANY);

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

54
    if (fill_flag.isSet(Mesh::CALL_LEAF_EL_LEVEL)) {
55
      TEST_EXIT_DBG(level >= 0)("invalid level: %d\n", level);   
56
    }
57
58

    traverse_mel = NULL;
59
    stack_used = 0;
60

61
    return traverseNext(NULL);
62
63
  }

64

65
66
67
68
69
70
71
72
73
74
75
76
77
78
  ElInfo* TraverseStack::traverseFirstOneMacro(Mesh *mesh, int macroIndex, int level, 
					       Flag fill_flag)
  {
    FUNCNAME("TraverseStack::traverseFirstOneMacro()");

    TEST_EXIT_DBG(macroIndex >= 0)("Invalid macro element index!\n");
    TEST_EXIT_DBG(traverse_fill_flag.isSet(Mesh::CALL_MG_LEVEL) == false)
      ("Multigrid level traverse not supported for one macro element only!\n");

    limitedToMacroElement = macroIndex;
    return TraverseStack::traverseFirst(mesh, level, fill_flag);
  }


79
80
  ElInfo* TraverseStack::traverseNext(ElInfo* elinfo_old)
  {
81
    FUNCNAME("TraverseStack::traverseNext()");
82

83
    ElInfo *elinfo = NULL;
84
85
86
87
88
89
    Parametric *parametric = traverse_mesh->getParametric();

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

90
      TEST_EXIT_DBG(elinfo_old == elinfo_stack[stack_used])("invalid old elinfo\n");
91
    } else {
92
      TEST_EXIT_DBG(elinfo_old == NULL)("invalid old elinfo != nil\n");
93
94
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
95
    if (traverse_fill_flag.isSet(Mesh::CALL_LEAF_EL)) {
96
      elinfo = traverseLeafElement();
Thomas Witkowski's avatar
Thomas Witkowski committed
97
    } else {
98
99
100
101
102
103
104
      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
105
	if (traverse_fill_flag.isSet(Mesh::CALL_EVERY_EL_PREORDER)) {
106
	  elinfo = traverseEveryElementPreorder();
107
	} else if (traverse_fill_flag.isSet(Mesh::CALL_EVERY_EL_INORDER))
108
109
110
111
112
	  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
113
    }
114

115
    if (elinfo) {
116
      if (parametric)
117
118
119
120
	elinfo = parametric->addParametricInfo(elinfo);
      elinfo->fillDetGrdLambda();
    }

121
    return elinfo;
Thomas Witkowski's avatar
Thomas Witkowski committed
122
123
  }

124

125
126
  void TraverseStack::enlargeTraverseStack()
  {
127
128
129
    FUNCNAME("TraverseStack::enlargeTraverseStack()");

    int new_stack_size = stack_size + 10;
130
131

    elinfo_stack.resize(new_stack_size, NULL);
132
 
133
    // create new elinfos
134
135
    for (int i = stack_size; i < new_stack_size; i++) {
      TEST_EXIT_DBG(elinfo_stack[i] == NULL)("???\n");
136
137
138
      elinfo_stack[i] = traverse_mesh->createNewElInfo();
    }

139
140
    if (stack_size > 0)
      for (int i = stack_size; i < new_stack_size; i++)
141
	elinfo_stack[i]->setFillFlag(elinfo_stack[0]->getFillFlag());
142
143
144
145
146

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

    // create new elinfos
147
148
    for (int i = stack_size; i < new_stack_size; i++) {
      TEST_EXIT_DBG(save_elinfo_stack[i] == NULL)("???\n");
149
150
151
152
153
154
155
156
157
158
      save_elinfo_stack[i] = traverse_mesh->createNewElInfo();
    }
    save_info_stack.resize(new_stack_size);

    stack_size = new_stack_size;    
  }

  
  ElInfo* TraverseStack::traverseLeafElement()
  {
159
    FUNCNAME("TraverseStack::traverseLeafElement()");
160

161
    Element *el = NULL;
162

163
164
    if (stack_used == 0) {   /* first call */
      currentMacro = traverse_mesh->firstMacroElement();
165

166
167
168
169
170
171
172
173
174
175
176
177
178
      if (limitedToMacroElement >= 0) {
	while ((*currentMacro)->getIndex() != limitedToMacroElement &&
	       currentMacro != traverse_mesh->endOfMacroElements())
	  currentMacro++;

	TEST_EXIT_DBG(currentMacro != traverse_mesh->endOfMacroElements())
	  ("Coult not find macro element with index %d!\n", limitedToMacroElement);

      } else {
	while (((*currentMacro)->getIndex() % maxThreads != myThreadId) &&
	       currentMacro != traverse_mesh->endOfMacroElements())
	  currentMacro++;      
      }
179

180
      if (currentMacro == traverse_mesh->endOfMacroElements())
181
	return NULL;
182

183
184
185
186
      traverse_mel = *currentMacro;
      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
      info_stack[stack_used] = 0;
187

188
      el = elinfo_stack[stack_used]->getElement();
189
      if (el == NULL || el->getFirstChild() == NULL)
190
	return elinfo_stack[stack_used];
191
192
193
194
    } else {
      el = elinfo_stack[stack_used]->getElement();
      
      /* go up in tree until we can go down again */
195
196
197
198
199
      while ((stack_used > 0) && 
	     ((info_stack[stack_used] >= 2) || (el->getFirstChild() == NULL))) {
	stack_used--;
	el = elinfo_stack[stack_used]->getElement();
      }
200
201
202
      
      /* goto next macro element */
      if (stack_used < 1) {
203
204
205
	if (limitedToMacroElement >= 0)
	  return NULL;

206
207
	do {	
	  currentMacro++;
Thomas Witkowski's avatar
Thomas Witkowski committed
208
	} while ((currentMacro != traverse_mesh->endOfMacroElements()) && 
209
		 ((*currentMacro)->getIndex() % maxThreads != myThreadId));
210

211
	if (currentMacro == traverse_mesh->endOfMacroElements())
212
	  return NULL;
213
214

	traverse_mel = *currentMacro;	
215
216
	stack_used = 1;
	elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
217
	info_stack[stack_used] = 0;	
218
	el = elinfo_stack[stack_used]->getElement();
219
220

	if (el == NULL || el->getFirstChild() == NULL)
221
	  return elinfo_stack[stack_used];
222
      }
223
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
224
   
225
    /* go down tree until leaf */
226
    while (el->getFirstChild()) {
227
      if (stack_used >= stack_size - 1) 
228
229
	enlargeTraverseStack();
      
Thomas Witkowski's avatar
Thomas Witkowski committed
230
231
      int i = info_stack[stack_used];
      el = const_cast<Element*>(((i == 0) ? el->getFirstChild() : el->getSecondChild()));
232
      info_stack[stack_used]++;
233
      elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
234
      stack_used++;
235

236
      TEST_EXIT_DBG(stack_used < stack_size)
237
	("stack_size = %d too small, level = (%d, %d)\n",
Thomas Witkowski's avatar
Thomas Witkowski committed
238
239
240
	 stack_size, elinfo_stack[stack_used]->getLevel());
      
      info_stack[stack_used] = 0;
241
    }
242

243
244
245
    return elinfo_stack[stack_used];
  }

246

247
248
  ElInfo* TraverseStack::traverseLeafElementLevel()
  {
249
    FUNCNAME("TraverseStack::traverseLeafElementLevel()");
250

251
    ERROR_EXIT("not yet implemented\n");
252

253
254
255
    return NULL;
  }

256

257
258
  ElInfo* TraverseStack::traverseElementLevel()
  {
259
    FUNCNAME("TraverseStack::traverseElementLevel()");
260
261
262
263
264
265

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

266
    return elInfo;
267
268
  }

269

270
271
  ElInfo* TraverseStack::traverseMultiGridLevel()
  {
272
    FUNCNAME("TraverseStack::traverseMultiGridLevel()");
273

274
275
276
    if (stack_used == 0) {   /* first call */
      currentMacro = traverse_mesh->firstMacroElement();
      traverse_mel = *currentMacro;
277
278
      if (traverse_mel == NULL)  
	return NULL;
279
280
281
282
283
284
285
286
      
      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()))
287
	return elinfo_stack[stack_used];
288
    }
289
  
290
    Element *el = elinfo_stack[stack_used]->getElement();
291
292

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


    /* goto next macro element */
    if (stack_used < 1) {
      currentMacro++;
303
      if (currentMacro == traverse_mesh->endOfMacroElements()) 
304
	return NULL;
305

306
      traverse_mel = *currentMacro;
307
308
309
310
      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
      info_stack[stack_used] = 0;

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


    /* go down tree */

320
321
322
323
    if (stack_used >= stack_size - 1) 
      enlargeTraverseStack();

    int i = info_stack[stack_used];
324
325
    info_stack[stack_used]++;

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

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

    info_stack[stack_used] = 0;
  
335
336
337
    if ((elinfo_stack[stack_used]->getLevel() == traverse_level) ||
	(elinfo_stack[stack_used]->getLevel() < traverse_level && 
	 elinfo_stack[stack_used]->getElement()->isLeaf()))
338
      return elinfo_stack[stack_used];
339
340
341
342

    return traverseMultiGridLevel();
  }

343

344
345
  ElInfo* TraverseStack::traverseEveryElementPreorder()
  {
346
    FUNCNAME("TraverseStack::traverseEveryElementPreorder()");
347

348
349
    if (stack_used == 0) {   /* first call */
      currentMacro = traverse_mesh->firstMacroElement();
350
351
352
353
354
355
356
357
358
359

      if (limitedToMacroElement >= 0) {
	while ((*currentMacro)->getIndex() != limitedToMacroElement &&
	       currentMacro != traverse_mesh->endOfMacroElements())
	  currentMacro++;

	TEST_EXIT_DBG(currentMacro != traverse_mesh->endOfMacroElements())
	  ("Coult not find macro element with index %d!\n", limitedToMacroElement);
      }

360
361
362
363
364
365
366
      traverse_mel = *currentMacro;
      if (traverse_mel == NULL)  
	return NULL;
      
      stack_used = 1;
      elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
      info_stack[stack_used] = 0;
367

368
      return elinfo_stack[stack_used];
369
    }
370
  
371
    Element *el = elinfo_stack[stack_used]->getElement();
372
373

    /* go up in tree until we can go down again */
374
375
376
377
378
    while (stack_used > 0 && 
	   (info_stack[stack_used] >= 2 || el->getFirstChild() == NULL)) {
      stack_used--;
      el = elinfo_stack[stack_used]->getElement();
    }
379
380
381
382


    /* goto next macro element */
    if (stack_used < 1) {
383
384
385
      if (limitedToMacroElement >= 0)
	return NULL;

386
      currentMacro++;
387
      if (currentMacro == traverse_mesh->endOfMacroElements()) 
388
	return NULL;
389
390
391
392
393
394
      traverse_mel = *currentMacro;

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

395
      return elinfo_stack[stack_used];
396
397
398
399
400
    }


    /* go down tree */

401
    if (stack_used >= stack_size - 1) 
402
403
      enlargeTraverseStack();

404
    int fillIthChild = info_stack[stack_used];
405

406
    info_stack[stack_used]++;
407
408
    if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
      fillIthChild = 1 - fillIthChild;
409

410
    elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);
411

412
413
    stack_used++;

414
    TEST_EXIT_DBG(stack_used < stack_size)
415
      ("stack_size = %d too small, level = %d\n",
416
417
418
419
       stack_size, elinfo_stack[stack_used]->getLevel());

    info_stack[stack_used] = 0;
  
420
    return elinfo_stack[stack_used];
421
422
  }

423

424
425
426
427
428
429
430
  ElInfo* TraverseStack::traverseEveryElementInorder()
  {
    FUNCNAME("TraverseStack::traverseEveryElementInorder");
    ERROR_EXIT("not yet implemented\n");
    return NULL;
  }

431

432
433
  ElInfo* TraverseStack::traverseEveryElementPostorder()
  {
434
    FUNCNAME("TraverseStack::traverseEveryElementPostorder()");
435

436
437
    if (stack_used == 0) {   /* first call */
      currentMacro = traverse_mesh->firstMacroElement();
438
439
440
441
442
443
444
445
446
      if (limitedToMacroElement >= 0) {
	while ((*currentMacro)->getIndex() != limitedToMacroElement &&
	       currentMacro != traverse_mesh->endOfMacroElements())
	  currentMacro++;

	TEST_EXIT_DBG(currentMacro != traverse_mesh->endOfMacroElements())
	  ("Coult not find macro element with index %d!\n", limitedToMacroElement);
      }

447
448
449
450
451
452
453
454
455
456
      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 */
457
      Element *el = elinfo_stack[stack_used]->getElement();
458
459

      /* go up in tree until we can go down again */          /* postorder!!! */
460
461
      while (stack_used > 0 && 
	     (info_stack[stack_used] >= 3 || el->getFirstChild() == NULL)) {
462
463
464
	stack_used--;
	el = elinfo_stack[stack_used]->getElement();
      }
465
466
467
468


      /* goto next macro element */
      if (stack_used < 1) {
469
470
471
	if (limitedToMacroElement >= 0)
	  return NULL;

472
	currentMacro++;
473
474
	if (currentMacro == traverse_mesh->endOfMacroElements()) 
	  return NULL;
475
476
477
478
479
480
481
482
483
484
485
	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 */

486
487
    while (elinfo_stack[stack_used]->getElement()->getFirstChild() &&
	   info_stack[stack_used] < 2) {
488
      if (stack_used >= stack_size-1) 
489
490
	enlargeTraverseStack();

491
      int i = info_stack[stack_used];
492
      info_stack[stack_used]++;
493
      elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
494
495
496
497
498
499
      stack_used++;
      info_stack[stack_used] = 0;
    }
  
    info_stack[stack_used]++;      /* postorder!!! */

500
    return elinfo_stack[stack_used];
501
502
  }

503

504
505
506
507
508
509
510
511
512
513
514
515
516
  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;
517
518
    default:
      ERROR_EXIT("invalid dim\n");
519
520
521
522
    }
    return NULL;
  }

523

524
525
  ElInfo* TraverseStack::traverseNeighbour3d(ElInfo* elinfo_old, int neighbour)
  {
526
    FUNCNAME("TraverseStack::traverseNeighbour3d()");
Thomas Witkowski's avatar
Thomas Witkowski committed
527

528
    Element *el2;
529
530
531
    ElInfo *elinfo2;
    int stack2_used;
    int sav_neighbour = neighbour;
532

533
534
535
536
    // 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}}};
537

538
    TEST_EXIT_DBG(stack_used > 0)("no current element\n");
539
540

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

544
    TEST_EXIT_DBG(elinfo_old == elinfo_stack[stack_used])("invalid old elinfo\n");
545
    TEST_EXIT_DBG(elinfo_stack[stack_used]->getFillFlag().isSet(Mesh::FILL_NEIGH))
546
      ("FILL_NEIGH not set");
547

548
    Element *el = elinfo_stack[stack_used]->getElement();
549
    int sav_index = el->getIndex();
550

551
    // First, goto to leaf level, if necessary ...
552
    if ((traverse_fill_flag & Mesh::CALL_LEAF_EL).isAnySet()) {
553
      if (el->getChild(0) && neighbour < 2) {
554
	if (stack_used >= stack_size - 1)
555
	  enlargeTraverseStack();
556
557
	int i = 1 - neighbour;

558
	elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
Thomas Witkowski's avatar
Thomas Witkowski committed
559
560
561
 	if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
 	  info_stack[stack_used] = (i == 0 ? 2 : 1);
 	else
562
	  info_stack[stack_used] = i + 1;
563
564
565
566
567
568
569
570
	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;
571
    save_stack_used = stack_used;
572

573
574
575
    // === First phase (see 2D). ===

    int nb = neighbour;
576

577
578
    while (stack_used > 1) { /* go up in tree until we can go down again */
      stack_used--;
579
580
581
582
583
584
585
586
587
588
589
      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];
590

591
592
      if (nb == -1) 
	break;
593

594
      TEST_EXIT_DBG(nb >= 0)("Invalid coarse_nb %d!\n", nb);
595
    }
596

Thomas Witkowski's avatar
Thomas Witkowski committed
597
598
599
600
601
602
603
    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);

604

605
606
    if (nb >= 0) {                           
      // Go to macro element neighbour.
607

608
      int i = traverse_mel->getOppVertex(nb);
609

610
      traverse_mel = traverse_mel->getNeighbour(nb);
611
      if (traverse_mel == NULL)  
612
	return NULL;
613
    
614
      if (nb < 2 && save_stack_used > 1)
615
	stack2_used = 2;                /* go down one level in OLD hierarchy */
616
      else
617
	stack2_used = 1;
618
      
619
620
621
622
623
624
      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;
625
626
627
    } else {                                                
      // Goto other child.

628
      stack2_used = stack_used + 1;
629
      if (save_stack_used > stack2_used)
630
	stack2_used++;                  /* go down one level in OLD hierarchy */
631

632
633
634
      elinfo2 = save_elinfo_stack[stack2_used];
      el2 = elinfo2->getElement();

635
      if (stack_used >= stack_size - 1)
636
	enlargeTraverseStack();
637
638
639
640
641
642

      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;
643
      elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);
644
645
646
647
648
      stack_used++;
      info_stack[stack_used] = 0;
      nb = 0;
    }

649
650
651
652

    // === Second phase. ===

    ElInfo *elinfo = elinfo_stack[stack_used];
653
654
    el = elinfo->getElement();

655
656
657
658
659
    while (el->getChild(0)) {
      if (nb < 2) {                         
	// Go down one level in hierarchy.

	if (stack_used >= stack_size - 1)
660
	  enlargeTraverseStack();
661
662
663
664
665
666
667
668


	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;	
	}
669
670
671
672

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

673
674
675
676
677
678
679
	stack_used++;
	info_stack[stack_used] = 0;
	elinfo = elinfo_stack[stack_used];
	el = elinfo->getElement();
	nb = 3;
      }

680
681
682
      if (save_stack_used > stack2_used) { 
	// `refine' both el and el2.

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

685
	int i = 0;
686
	if (el->getDof(0) == el2->getDof(0))
687
	  i = save_info_stack[stack2_used] - 1;
688
	else if (el->getDof(1) == el2->getDof(0))
689
690
	  i = 2 - save_info_stack[stack2_used];
	else {
691
	  if (traverse_mesh->associated(el->getDof(0, 0), el2->getDof(0, 0)))
692
	    i = save_info_stack[stack2_used] - 1;
693
	  else if (traverse_mesh->associated(el->getDof(1, 0), el2->getDof(0, 0)))
694
695
	    i = 2 - save_info_stack[stack2_used];
	  else {
696
	    ERROR_EXIT("No common refinement edge! %d\n", traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE));
697
698
699
	  }
	}

700
701
702
703
	int testChild = i;
 	if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
 	  testChild = 1 - testChild;

704
	if (el->getChild(0) &&  
705
706
	    (el->getChild(testChild)->getDof(1) == el->getDof(nb) ||
	     traverse_mesh->associated(el->getChild(testChild)->getDof(1, 0), el->getDof(nb, 0))))
707
708
709
	  nb = 1;	
	else
	  nb = 2;	
710

711
712
	info_stack[stack_used] = i + 1;	

713
	if (stack_used >= stack_size - 1)
714
	  enlargeTraverseStack();
715
716
717
718
719
720

	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]);

721
722
723
724
725
726
727
728
729
	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();
730

731
	if (save_stack_used > stack2_used) {
732
733
734
735
	  const DegreeOfFreedom *dof = el2->getDof(1);
	  if (dof != el->getDof(1) && dof != el->getDof(2) &&
	      !traverse_mesh->associated(dof[0], el->getDof(1, 0)) &&
	      !traverse_mesh->associated(dof[0], el->getDof(2, 0))) {
736
737
738
739
740
	  
	    stack2_used++;               /* go down one level in OLD hierarchy */
	    elinfo2 = save_elinfo_stack[stack2_used];
	    el2 = elinfo2->getElement();
	  } 
741
	}
742
743
      } else {   
	// Now we're done...
744
745
746

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

748
749
750
751
752
753
754
	break;
      }
    }


    if (elinfo->getNeighbour(opp_vertex) != old_elinfo->getElement()) {
      MSG(" looking for neighbour %d of element %d at %p\n",
755
	  neighbour, old_elinfo->getElement()->getIndex(), reinterpret_cast<void*>(old_elinfo->getElement()));
756
      MSG(" originally: neighbour %d of element %d at %p\n",
757
	  sav_neighbour, sav_index, reinterpret_cast<void*>(old_elinfo->getElement()));
758
      MSG(" got element %d at %p with opp_vertex %d neigh %d\n",
759
	  elinfo->getElement()->getIndex(), reinterpret_cast<void*>(elinfo->getElement()),
760
	  opp_vertex, (elinfo->getNeighbour(opp_vertex))?(elinfo->getNeighbour(opp_vertex))->getIndex():-1);
761
      TEST_EXIT_DBG(elinfo->getNeighbour(opp_vertex) == old_elinfo->getElement())
762
763
764
765
766
767
768
769
770
	("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;  /* ??? */

771
772
773
    if (elinfo) {
      if (parametric) 
	elinfo = parametric->addParametricInfo(elinfo);
774
775
      elinfo->fillDetGrdLambda();
    }
776
777

    return elinfo;
778
779
  }

780

781
782
  ElInfo* TraverseStack::traverseNeighbour2d(ElInfo* elinfo_old, int neighbour)
  {
783
784
    FUNCNAME("TraverseStack::traverseNeighbour2d()");

785
786
787
788
    Triangle *el2;
    ElInfo *elinfo2;
    int stack2_used;
    int sav_neighbour = neighbour;
789

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

793
    TEST_EXIT_DBG(stack_used > 0)("no current element");
794
795
796
797
798

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

799
    TEST_EXIT_DBG(elinfo_old == elinfo_stack[stack_used])("invalid old elinfo");
800
801

    elinfo_stack[stack_used]->testFlag(Mesh::FILL_NEIGH);
802
803
804
805
    Triangle *el = 
      dynamic_cast<Triangle*>(const_cast<Element*>(elinfo_stack[stack_used]->getElement()));
    int sav_index = el->getIndex();
    Triangle *sav_el = el;
806
807

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

810
811
812
813
814
815
816
817
818
819
820
821
822
      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;
      
823
      stack_used++;
824
      neighbour = 2;     
825
826
827
828
    }

    /* save information about current element and its position in the tree */
    save_traverse_mel = traverse_mel;
829
830
    save_stack_used = stack_used;
    for (int i = 0; i <= stack_used; i++) {
831
832
      save_info_stack[i] = info_stack[i];
      (*(save_elinfo_stack[i])) = (*(elinfo_stack[i]));
833
834
835
    }
    ElInfo *old_elinfo = save_elinfo_stack[stack_used];
    int opp_vertex = old_elinfo->getOppVertex(neighbour);
836
837
838
839
840
841


    /****************************************************************************/
    /* 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
842
    /* element of the OLD hierarchy branch to the new branch                    */
843
844
    /****************************************************************************/

845
    int nb = neighbour;
846

847
848
    while (stack_used > 1) {
      stack_used--;
849
      int elIsIthChild = info_stack[stack_used];
850

851
852
853
854
855
856
857
858
859
      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];
860
      if (nb == -1) 
861
862
	break;     

863
864
      TEST_EXIT_DBG(nb >= 0)("invalid coarse_nb %d\n",nb);
    }
865

866

867
868
869
870
871
872
    /****************************************************************************/
    /* 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)                               */
    /****************************************************************************/

873
874
    if (nb >= 0) {                        
      // Go to macro element neighbour.
875

876
      if (nb < 2 && save_stack_used > 1)
877
	stack2_used = 2;           /* go down one level in OLD hierarchy */
878
      else
879
	stack2_used = 1;
880
      
881
      elinfo2 = save_elinfo_stack[stack2_used];
882
      el2 = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo2->getElement()));
883

884
      int i = traverse_mel->getOppVertex(nb);
885
      traverse_mel = traverse_mel->getNeighbour(nb);
886
887
      if (traverse_mel == NULL)
	return NULL;
888
      nb = i;
889

890
      stack_used = 1;
891
      elinfo_stack[stack_used]->fillMacroInfo(const_cast<MacroElement*>(traverse_mel));
892
      info_stack[stack_used] = 0;
893
894
    } else {                                               
      // Go to other child.
895
896

      stack2_used = stack_used + 1;
897
      if (save_stack_used > stack2_used)
898
	stack2_used++;               /* go down one level in OLD hierarchy */
899