Tetrahedron.cc 16.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
20
21


22
23
24
25
26
#include "Tetrahedron.h"
#include "DOFAdmin.h"
#include "Mesh.h"
#include "CoarseningManager.h"
#include "FixVec.h"
27
#include "ElementDofIterator.h"
28
#include "Global.h"
29
30
31
32
33
34

namespace AMDiS {

  const unsigned char Tetrahedron::nChildEdge[3][2][2] = {{{5,4},{4,5}},
							  {{5,4},{5,4}},
							  {{5,4},{5,4}}};
35

36
37
38
  const unsigned char Tetrahedron::nChildFace[3][2][2] = {{{1,2},{2,1}},
							  {{1,2},{1,2}},
							  {{1,2},{1,2}}};
39

40
41
42
  const int Tetrahedron::childVertex[3][2][4] = {{{0,2,3,4},{1,3,2,4}},
						 {{0,2,3,4},{1,2,3,4}},
						 {{0,2,3,4},{1,2,3,4}}};
43

44
45
46
  const int Tetrahedron::childEdge[3][2][6] = {{{1,2,0,5,5,4}, {4,3,0,5,5,4}},
					       {{1,2,0,5,4,5}, {3,4,0,5,4,5}},
					       {{1,2,0,5,4,5}, {3,4,0,5,4,5}}};
47

48
49
50
51
52
53
  const unsigned char Tetrahedron::adjacentChild[2][2] = {{0,1}, {1,0}};

  const signed char Tetrahedron::childOrientation[3][2] = {{1,1}, 
							   {1,-1}, 
							   {1,-1}};

54
55
56
57
  const unsigned char Tetrahedron::edgeOfDofs[4][4] =  {{255, 0, 1, 2},
							{0, 255, 3, 4},
							{1, 3, 255, 5},
							{2, 4, 5, 255}};
58
59
60
61
62
63
64

  const int Tetrahedron::vertexOfEdge[6][2] = {{0,1},
					       {0,2},
					       {0,3},
					       {1,2},
					       {1,3},
					       {2,3}};
65

66
67
68
69
70
71
72
73
74
75
76
77
  const int Tetrahedron::vertexOfFace[4][3] = {{1,2,3},
					       {0,2,3},
					       {0,1,3},
					       {0,1,2}};

  const int Tetrahedron::sideOfChild[3][2][4] = {{{-1, 3, 1, 2},  // type 0
						  {3, -1, 2, 1}},
						 {{-1, 3, 1, 2},  // type 1
						  {3, -1, 1, 2}},
						 {{-1, 3, 1, 2},  // type 2
						  {3, -1, 1, 2}}};

78
79
80
81
82
83
84
  const int Tetrahedron::edgeOfChild[3][2][6] = {{{2, 0, 1, -1, -1, 3},   // type 0
						  {2, -1, -1, 1, 0, 3}},
						 {{2, 0, 1, -1, -1, 3},   // type 1
						  {2, -1, -1, 0, 1, 3}},
						 {{2, 0, 1, -1, -1, 3},   // type 2
						  {2, -1, -1, 0, 1, 3}}};

85
86
87
88
89
90
91
92
93
94
95
96
  const int Tetrahedron::vertexOfParent[3][2][4] = {{{0, 2, 3, -1},  // type 0
						     {1, 3, 2, -1}},
						    {{0, 2, 3, -1},  // type 1
						     {1, 2, 3, -1}},
						    {{0, 2, 3, -1},  // type 2
						     {1, 2, 3, -1}}};

  const int Tetrahedron::edgeOfFace[4][3] = {{5, 4, 3},  // face 0
					     {5, 2, 1},  // face 1
					     {4, 2, 0},  // face 2
					     {3, 1, 0}}; // face 3

97

98
99
  bool Tetrahedron::hasSide(Element* sideElem) const
  {
100
101
    FUNCNAME("Tetrahedron::hasSide()");
    TEST_EXIT_DBG(sideElem->isTriangle())("called for sideElem-type != Triangle\n");
102
103
104
105
    ERROR_EXIT("not yet\n");
    return false;
  }

106

107
  int Tetrahedron::getVertexOfPosition(GeoIndex position,
108
109
				       int positionIndex,
				       int vertexIndex) const 
110
  {
111
    FUNCNAME("Triangle::getVertexOfPosition()");
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
    switch(position) {
    case VERTEX:
      return positionIndex;
      break;
    case EDGE:
      return vertexOfEdge[positionIndex][vertexIndex];
      break;
    case FACE:
      return vertexOfFace[positionIndex][vertexIndex];
      break;
    case CENTER:
      return vertexIndex;
      break;
    default:
      ERROR_EXIT("invalid position\n");
      return 0;
    }
  }

131

132
  void Tetrahedron::sortFaceIndices(int face, FixVec<int,WORLD> &vec) const
133
  {
134
    // TODO: REMOVE STATIC
135
    static MatrixOfFixVecs<FixVec<int,WORLD> > *sorted_3d = nullptr;
136

137
    if (sorted_3d == nullptr) {
138
      sorted_3d = new MatrixOfFixVecs<FixVec<int,WORLD> >(3, 4, 7, NO_INIT);
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
      (*sorted_3d)[0][0][0] = (*sorted_3d)[0][0][1] =
	(*sorted_3d)[0][0][2] = (*sorted_3d)[1][0][0] =
	(*sorted_3d)[1][0][1] = (*sorted_3d)[1][0][2] =
	(*sorted_3d)[1][1][0] = (*sorted_3d)[1][2][1] =
	(*sorted_3d)[1][3][0] = (*sorted_3d)[1][4][2] =
	(*sorted_3d)[1][5][1] = (*sorted_3d)[1][6][2] =
	(*sorted_3d)[2][0][0] = (*sorted_3d)[2][0][1] =
	(*sorted_3d)[2][0][2] = (*sorted_3d)[2][1][0] =
	(*sorted_3d)[2][2][1] = (*sorted_3d)[2][3][0] =
	(*sorted_3d)[2][4][2] = (*sorted_3d)[2][5][1] =
	(*sorted_3d)[2][6][2] = (*sorted_3d)[3][0][0] =
	(*sorted_3d)[3][0][1] = (*sorted_3d)[3][0][2] =
	(*sorted_3d)[3][1][0] = (*sorted_3d)[3][2][1] =
	(*sorted_3d)[3][3][0] = (*sorted_3d)[3][4][2] =
	(*sorted_3d)[3][5][1] = (*sorted_3d)[3][6][2] = 0;

      (*sorted_3d)[0][1][0] = (*sorted_3d)[0][2][1] =
	(*sorted_3d)[0][3][0] = (*sorted_3d)[0][4][2] =
	(*sorted_3d)[0][5][1] = (*sorted_3d)[0][6][2] =
	(*sorted_3d)[2][1][2] = (*sorted_3d)[2][2][0] =
	(*sorted_3d)[2][3][1] = (*sorted_3d)[2][4][1] =
	(*sorted_3d)[2][5][2] = (*sorted_3d)[2][6][0] =
	(*sorted_3d)[3][1][2] = (*sorted_3d)[3][2][0] =
	(*sorted_3d)[3][3][1] = (*sorted_3d)[3][4][1] =
	(*sorted_3d)[3][5][2] = (*sorted_3d)[3][6][0] = 1;

      (*sorted_3d)[0][1][2] = (*sorted_3d)[0][2][0] =
	(*sorted_3d)[0][3][1] = (*sorted_3d)[0][4][1] =
	(*sorted_3d)[0][5][2] = (*sorted_3d)[0][6][0] =
	(*sorted_3d)[1][1][2] = (*sorted_3d)[1][2][0] =
	(*sorted_3d)[1][3][1] = (*sorted_3d)[1][4][1] =
	(*sorted_3d)[1][5][2] = (*sorted_3d)[1][6][0] =
	(*sorted_3d)[3][1][1] = (*sorted_3d)[3][2][2] =
	(*sorted_3d)[3][3][2] = (*sorted_3d)[3][4][0] =
	(*sorted_3d)[3][5][0] = (*sorted_3d)[3][6][1] = 2;

      (*sorted_3d)[0][1][1] = (*sorted_3d)[0][2][2] =
	(*sorted_3d)[0][3][2] = (*sorted_3d)[0][4][0] =
	(*sorted_3d)[0][5][0] = (*sorted_3d)[0][6][1] =
	(*sorted_3d)[1][1][1] = (*sorted_3d)[1][2][2] =
	(*sorted_3d)[1][3][2] = (*sorted_3d)[1][4][0] =
	(*sorted_3d)[1][5][0] = (*sorted_3d)[1][6][1] =
	(*sorted_3d)[2][1][1] = (*sorted_3d)[2][2][2] =
	(*sorted_3d)[2][3][2] = (*sorted_3d)[2][4][0] =
	(*sorted_3d)[2][5][0] = (*sorted_3d)[2][6][1] = 3;
185
186
    }

187
    const int *vof = vertexOfFace[face];
188
    int no = 0;
189
190
191
192
193
194
195
    if (dof[vof[0]][0] < dof[vof[1]][0])
      no++;
    if (dof[vof[1]][0] < dof[vof[2]][0])
      no += 2;
    if (dof[vof[2]][0] < dof[vof[0]][0])
      no += 4;

196
197
    TEST_EXIT_DBG(no >= 1 && no <= 6)
      ("Cannot sort face indices of element %d at face %d\n", getIndex(), face);
198

199
    vec = (*sorted_3d)[face][no];
200
201
  }

202

203
  void Tetrahedron::getNodeDofs(const FiniteElemSpace* feSpace, 
204
				BoundaryObject bound,
205
206
				DofContainer& dofs,
				bool baseDofPtr) const
207
  {
208
    FUNCNAME("Tetrahedron::getNodeDofs()");
209

210
    switch (bound.subObj) {
Thomas Witkowski's avatar
Thomas Witkowski committed
211
    case VERTEX:
212
      {
213
214
215
216
217
218
	if (baseDofPtr) {
	  dofs.push_back(dof[bound.ithObj]);
	} else {
	  int n0 = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
	  dofs.push_back(&(dof[bound.ithObj][n0]));
	}
219
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
220
      break;
221
    case EDGE:
222
      getNodeDofsAtEdge(feSpace, bound, dofs, baseDofPtr);
223
224
      break;
    case FACE:
225
      getNodeDofsAtFace(feSpace, bound, dofs, baseDofPtr);
226
227
228
229
230
231
232
      break;      
    default:
      ERROR_EXIT("Should not happen!\n");
    }
  }


233
  void Tetrahedron::getNodeDofsAtFace(const FiniteElemSpace* feSpace, 
234
				      BoundaryObject bound,
235
236
				      DofContainer& dofs,
				      bool baseDofPtr) const
237
  {
238
    FUNCNAME("Tetrahedron::getNodeDofsAtFace()");
239

240
    if (!child[0]){
241
      return;
242
	}
243
244

    switch (bound.ithObj) {
245
246
    case 0:
      {
247
248
	BoundaryObject nextBound = bound;
	prepareNextBound(nextBound, 1);
249
	child[1]->getNodeDofs(feSpace, nextBound, dofs, baseDofPtr);
250
251
252
253
      }
      break;
    case 1:
      {
254
255
	BoundaryObject nextBound = bound;
	prepareNextBound(nextBound, 0);
256
	child[0]->getNodeDofs(feSpace, nextBound, dofs, baseDofPtr);
257
258
259
260
261
262
      }
      break;

    case 2:
    case 3:
      {
Thomas Witkowski's avatar
Thomas Witkowski committed
263
	int n0 = (baseDofPtr ? 0 : feSpace->getAdmin()->getNumberOfPreDofs(VERTEX));
264
265
266
267
268
269
270
	BoundaryObject nextBound0 = bound, nextBound1 = bound;
	prepareNextBound(nextBound0, 0);
	prepareNextBound(nextBound1, 1);
	nextBound0.reverseMode = false;
	nextBound1.reverseMode = false;
	
	bool addDof = true;
271
272
	for (ExcludeList::iterator it = bound.excludedSubstructures.begin(); 
	     it != bound.excludedSubstructures.end(); ++it)
273
	  if (it->first == EDGE && it->second == 0)
274
	    addDof = false;	
275
276
	
	if (bound.reverseMode) {
277
	  child[1]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr);
278
279
	  
	  if (addDof)
280
	    dofs.push_back(&(child[0]->getDof()[3][n0]));
281
	  
282
	  child[0]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr);
283
	} else {
284
	  child[0]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr);
285
286
	  
	  if (addDof)
287
	    dofs.push_back(&(child[0]->getDof()[3][n0]));
288
	  
289
	  child[1]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr);
290
291
292
293
	}
      }
      break;
    default:
294
      ERROR_EXIT("Should never happen: %d\n", bound.ithObj);
295
    }
296
297
  }

298

299
  void Tetrahedron::getNodeDofsAtEdge(const FiniteElemSpace* feSpace, 
300
				      BoundaryObject bound,
301
302
				      DofContainer& dofs,
				      bool baseDofPtr) const
303
  {
304
    FUNCNAME_DBG("Tetrahedron::getNodeDofsAtEdge()");
305

306
307
308
309

    if (!child[0]){
      return;
	}
310

Thomas Witkowski's avatar
Thomas Witkowski committed
311
    int n0 = (baseDofPtr ? 0 : feSpace->getAdmin()->getNumberOfPreDofs(VERTEX));
312
313
314
315
    BoundaryObject nextBound0 = bound, nextBound1 = bound;
    prepareNextBound(nextBound0, 0);
    prepareNextBound(nextBound1, 1);

316
317
318
319
320
321
    switch (bound.ithObj) {
    case 0:      
      nextBound0.reverseMode = false;
      nextBound1.reverseMode = false;

      if (bound.reverseMode) {
322
323
324
325
			child[1]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr);
			dofs.push_back(&(child[0]->getDof()[3][n0]));
			child[0]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr);

326
      } else {
327
328
329
330
			child[0]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr);
			dofs.push_back(&(child[0]->getDof()[3][n0]));
			child[1]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr);

331
332
333
334
      }

      break;
    case 5:
335
336
      TEST_EXIT_DBG(nextBound0.ithObj == nextBound1.ithObj)
	("Should not happen!\n");
337
338

      if (nextBound0.ithObj != -1)
339
	child[0]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr);
340
341
342
343
344
345
      break;
    default:
      TEST_EXIT_DBG(nextBound0.ithObj == -1 || nextBound1.ithObj == -1)
	("This should not happen!\n");
      
      if (nextBound0.ithObj != -1)
346
	child[0]->getNodeDofs(feSpace, nextBound0, dofs, baseDofPtr);
347
348

      if (nextBound1.ithObj != -1)
349
	child[1]->getNodeDofs(feSpace, nextBound1, dofs, baseDofPtr);
350
    }
351

352
353
354
  }


355
  void Tetrahedron::getHigherOrderDofs(const FiniteElemSpace* feSpace,
356
				       BoundaryObject bound,
357
				       DofContainer& dofs,
358
359
				       bool baseDofPtr,
				       vector<GeoIndex>* dofGeoIndex) const
360
  {
361
    FUNCNAME("Tetrahedron::getHigherOrderDofs()");
362

363
364
365
366
367
    switch (bound.subObj) {
    case VERTEX:
      return;
      break;
    case EDGE:
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
      {
	// === Create boundary information objects for children elements. ===

	BoundaryObject nextBound0 = bound;
	prepareNextBound(nextBound0, 0);

	BoundaryObject nextBound1 = bound;
	prepareNextBound(nextBound1, 1);

	// === Check for boundary on children elements. ===

	if ((nextBound0.ithObj >= 0 || nextBound1.ithObj >= 0) &&  child[0]) {
	  // So, the edge is contained in at least on of the children and the
	  // element is also refined. Then we have go down further in refinement
	  // hierarchie.

384
385
386
387
388
389
	 	// Additonal Check if Boundary Edge is NOT Refinement Edge of BoundaryObject
	 	// AND the complete Edge is part of both childs
		// therefore only recursion on ONE child is necessary
	  if ( (bound.ithObj > 0) && (nextBound0.ithObj >= 0 && nextBound1.ithObj >= 0) ){
			child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs,
						   baseDofPtr, dofGeoIndex);
390
	  } else {
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
	  	if (bound.reverseMode) {
	    	if (nextBound1.ithObj >= 0)
	    	  child[1]->getHigherOrderDofs(feSpace, nextBound1, dofs, 
						   baseDofPtr, dofGeoIndex);
	    	if (nextBound0.ithObj >= 0)
	    	  child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs,
						   baseDofPtr, dofGeoIndex);
		  } else {
		    if (nextBound0.ithObj >= 0)
		      child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs,
						   baseDofPtr, dofGeoIndex);
		    if (nextBound1.ithObj >= 0)
		      child[1]->getHigherOrderDofs(feSpace, nextBound1, dofs,
						   baseDofPtr, dofGeoIndex);
		  }
406
	  }
407

408
409
410
411
412
413
414
	} else {
	  // Either the edge is not contained in further refined children, or
	  // the element is not refined further on this edge. Then we can get
	  // all the DOFs on this edge.

	  ElementDofIterator elDofIter(feSpace, true);
	  elDofIter.reset(this);
415
416
417
418

	  if (baseDofPtr) {
	    do {
	      if (elDofIter.getCurrentPos() == 1 && 
Thomas Witkowski's avatar
Thomas Witkowski committed
419
		  elDofIter.getCurrentElementPos() == bound.ithObj) {
420
		dofs.push_back(elDofIter.getBaseDof());	
421
		if (dofGeoIndex != nullptr)
Thomas Witkowski's avatar
Thomas Witkowski committed
422
423
		  dofGeoIndex->push_back(EDGE);
	      }
424
425
426
427
	    } while (elDofIter.nextStrict());
	  } else {
	    do {
	      if (elDofIter.getCurrentPos() == 1 && 
Thomas Witkowski's avatar
Thomas Witkowski committed
428
		  elDofIter.getCurrentElementPos() == bound.ithObj) {
429
430
				dofs.push_back(elDofIter.getDofPtr());	

431
				if (dofGeoIndex != nullptr)
432
				  dofGeoIndex->push_back(EDGE);
Thomas Witkowski's avatar
Thomas Witkowski committed
433
	      }
434
435
	    } while (elDofIter.next());
	  }
436
437
438
	}
      }

439
440
441
442
      break;
    case FACE:
      {      
	if (child[0]) {
443
444
445
446
447
	  BoundaryObject nextBound0 = bound, nextBound1 = bound;
	  prepareNextBound(nextBound0, 0);
	  prepareNextBound(nextBound1, 1);

	  TEST_EXIT_DBG(nextBound0.ithObj != -1 || nextBound1.ithObj != 1)
448
	    ("No new face for child elements!\n");
449
450

	  TEST_EXIT_DBG(!bound.reverseMode)("Not yet implemented!\n");
451
	  
452
453
	  if (nextBound0.ithObj != -1)
	    child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs,
454
					 baseDofPtr, dofGeoIndex);
455
	  
456
457
	  if (nextBound1.ithObj != -1)
	    child[1]->getHigherOrderDofs(feSpace, nextBound1, dofs,
458
					 baseDofPtr, dofGeoIndex);
459
	} else { 
Thomas Witkowski's avatar
Thomas Witkowski committed
460

461
462
463
	  // === On leaf elements iterate over all DOFs and collect the  === 
	  // === right ones.                                             ===

464
465
	  ElementDofIterator elDofIter(feSpace, true);
	  elDofIter.reset(this);
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503

	  bool next = true;
	  do {
	    bool addDof = false;

	    switch (elDofIter.getCurrentPos()) {
	    case 0:
	      // VERTEX is never a higher order DOF
	      addDof = false;
	      break;
	    case 1:
	      // On EDGE, first check whether this edge is part of the FACE and
	      // then check if this EDGE is not inside of excludedSubstructures.

	      {
		int localEdge = elDofIter.getCurrentElementPos();
		if (edgeOfFace[bound.ithObj][0] == localEdge ||
		    edgeOfFace[bound.ithObj][1] == localEdge ||
		    edgeOfFace[bound.ithObj][2] == localEdge) {
		  addDof = true;

		  for (ExcludeList::iterator it = bound.excludedSubstructures.begin(); 
		       it != bound.excludedSubstructures.end(); ++it)
		    if (it->first == EDGE && it->second == localEdge)
		      addDof = false;	
		} else {
		  addDof = false;
		}
	      }

	      break;
	    case 2:
	      // On FACE, check for the right one.
	      if (elDofIter.getCurrentElementPos() == bound.ithObj)
		addDof = true;
	      else
		addDof = false;
	      break;
504
505
506
507
		case 3:
		  // On CENTER, dof can not be part on internal boundary
		addDof = false;
		  break;
508
509
510
511
512
513
514
	    default:
	      ERROR_EXIT("Should not happen!\n");
	    }


	    if (addDof) {
	      if (baseDofPtr)
515
			dofs.push_back(elDofIter.getBaseDof());	
516
	      else
517
			dofs.push_back(elDofIter.getDofPtr());	
518
	      
519
	      if (dofGeoIndex != nullptr)
520
521
			dofGeoIndex->push_back(elDofIter.getPosIndex());

522
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
523

524
525
	    next = (baseDofPtr ? elDofIter.nextStrict() : elDofIter.next());	    
	  } while (next);
526
	}
527
      }
528
      break;
529
    default:
530
      ERROR_EXIT("Should not happen!\n");
531
    }
532

533
534
  }

535
536
537

  void Tetrahedron::getSubBoundary(BoundaryObject bound, 
				   vector<BoundaryObject> &subBound) const
538
  {}
539
540
541
542
543
544
545
546
547


  void Tetrahedron::prepareNextBound(BoundaryObject &bound, int ithChild) const
  {
    FUNCNAME("Tetrahedron::prepareNextBound()");

    switch (bound.subObj) {
    case FACE:      
      for (ExcludeList::iterator it = bound.excludedSubstructures.begin(); 
548
   	   it != bound.excludedSubstructures.end(); ++it)
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
	if (it->first == EDGE && it->second != -1)
	  it->second = edgeOfChild[bound.elType][ithChild][it->second];	
      
      bound.ithObj = sideOfChild[bound.elType][ithChild][bound.ithObj];
      bound.elType = (bound.elType + 1) % 3;

      break;
    case EDGE:
      bound.ithObj = edgeOfChild[bound.elType][ithChild][bound.ithObj];
      bound.elType = (bound.elType + 1) % 3;
      break;
    default:
      ERROR_EXIT("Should not happen!\n");
    }
  }

565
}