1 /****************************************************************************/
2 /* */
3 /* File: refine.c */
4 /* */
5 /* Purpose: unstructured grid adaption using a general element concept */
6 /* (dimension independent for 2/3D) */
7 /* */
8 /* Author: Stefan Lang */
9 /* Institut fuer Computeranwendungen III */
10 /* Universitaet Stuttgart */
11 /* Pfaffenwaldring 27 */
12 /* 70569 Stuttgart */
13 /* email: ug@ica3.uni-stuttgart.de */
14 /* */
15 /* History: 08.08.95 begin serial algorithm, ug version 3.0 */
16 /* */
17 /* Remarks: - level 0 grid consists of red elements only */
18 /* - the only restriction in the element hierarchy is that */
19 /* green or yellow elements might not have sons of class */
20 /* green or red */
21 /* - the rule set for refinement consists of regular (red) */
22 /* and irregular rules; regular rules create red elements */
23 /* while irregular rules result in green elements */
24 /* (green elements are needed for the closure of the grid, */
25 /* yellow elements, which are from copy rules, save */
26 /* the numerical properties of the solver and are handsome for */
27 /* the discretisation */
28 /* - if the rule set for the red rules is not complete for build-*/
29 /* up a consistent red refined region the FIFO might be used */
30 /* for some (hopefully not too much) iterations to find a */
31 /* consistent one */
32 /* - in 2D: exists a complete rule set for grids of triangles */
33 /* and quadrilaterals exclusively */
34 /* - in 3D: exists a complete rule set for tetrahedrons */
35 /* and we assume after some analysation a complete set */
36 /* of rules described by an algorithm for hexahedrons */
37 /* - for mixed element types in arbitrary dimension no */
38 /* rule set for the closure exists */
39 /* - BEFORE refinement we assume a situation where the error */
40 /* estimator has detected and marked the leaf elements for */
41 /* further refinement */
42 /* - AFTER refinement all elements are refined by a rule in way */
43 /* that no hanging nodes remain (this is the default mode) */
44 /* or with hanging nodes (in the hanging node mode) */
45 /* if you use inconsistent red refinement, you need to tell */
46 /* the algorithm explicitly to use the FIFO */
47 /* */
48 /****************************************************************************/
49
50 /****************************************************************************/
51 /* */
52 /* include files */
53 /* system include files */
54 /* application include files */
55 /* */
56 /****************************************************************************/
57
58 /* standard C library */
59 #include <config.h>
60 #include <cassert>
61 #include <cmath>
62 #include <cstdio>
63 #include <cstdlib>
64 #include <cstring>
65
66 #include <algorithm>
67
68 /* low module */
69 #include <dune/uggrid/low/debug.h>
70 #include <dune/uggrid/low/heaps.h>
71 #include <dune/uggrid/low/misc.h>
72 #include <dune/uggrid/low/ugtimer.h>
73 #include <dune/uggrid/low/ugtypes.h>
74
75 /* dev module */
76 #include <dune/uggrid/ugdevices.h>
77
78 /* gm module */
79 #include "algebra.h"
80 #include "evm.h"
81 #include "gm.h"
82 #include "cw.h"
83 #include "refine.h"
84 #include "elements.h"
85 #include "rm.h"
86 #include "ugm.h"
87 #include "mgheapmgr.h"
88
89 /* parallel modules */
90 #ifdef ModelP
91 #include <dune/uggrid/parallel/ppif/ppif.h>
92 #include <dune/uggrid/parallel/ddd/include/ddd.h>
93 #include <dune/uggrid/parallel/dddif/debugger.h>
94 #include <dune/uggrid/parallel/dddif/identify.h>
95 #include <dune/uggrid/parallel/dddif/parallel.h>
96 #include "pargm.h"
97 #endif
98
99 #ifdef Debug
100 #include <dune/uggrid/low/architecture.h>
101 #endif
102
103 #include <dune/uggrid/parallel/ppif/ppifcontext.hh>
104
105 USING_UG_NAMESPACES
106 using namespace PPIF;
107
108 /****************************************************************************/
109 /* */
110 /* defines in the following order */
111 /* */
112 /* compile time constants defining static data size (i.e. arrays) */
113 /* other constants */
114 /* macros */
115 /* */
116 /****************************************************************************/
117
118 /* undefine if overlap should be only updated where needed */
119 /* This does not work since the connection of the overlap needs the
120 fatherelements on both sides (ghost and master sons) and this is
121 not ensured.
122 #define UPDATE_FULLOVERLAP
123 */
124
125 /* define for use of DDD_ConsCheck() */
126 /*
127 #define DDD_CONSCHECK(context) if (DDD_ConsCheck(context)) assert(0);
128 */
129 #define DDD_CONSCHECK(context)
130
131 #define MINVNCLASS 2 /* determines copies, dep. on discr. ! */
132
133 /* defines for side matching of elements 8 bits: */
134 /* _ _ _ _ (4 bits for corner of one element) _ _ _ _ (4 bits for the other)*/
135
136 #define LINEPOINTS 51 /* 0011 0011 */
137 #define TRIPOINTS 119 /* 0111 0111 */
138 #define QUADPOINTS 255 /* 1111 1111 */
139
140 #define MAX_GREEN_SONS 32 /* max num of sons for green refinement */
141
142 /* define to control element which is responsible */
143 /* equal side refinement of neighboring elements */
144 #ifdef ModelP
145 #define _EID_ EGID
146 #define _ID_ GID
147 #else
148 #define _EID_ ID
149 #define _ID_ ID
150 #endif
151
152 #define EDGE_IN_PATTERN(p,i) (((p)[(i)]) & 0x1)
153 #define SIDE_IN_PATTERN(e,p,i) (((p)[EDGES_OF_ELEM(e)+(i)]) & 0x1)
154 #define EDGE_IN_PAT(p,i) (((p)>>(i)) & 0x1)
155 #define SIDE_IN_PAT(p,i) (((p)>>(i)) & 0x1)
156
157 #define MARK_BISECT_EDGE(r,i) ((r)->pattern[(i)]==1)
158
159 #define REF_TYPE_CHANGES(e) ((REFINE(e)!=MARK(e)) || \
160 (REFINECLASS(e)!=MARKCLASS(e)))
161 #define MARKED(e) (MARK(e)!=NO_REFINEMENT)
162
163 /* green marked elements were NEWGREEN is true are refined without rule */
164 #ifdef DUNE_UGGRID_TET_RULESET
165 #define NEWGREEN(e) (TAG(e)==HEXAHEDRON || TAG(e)== PRISM || \
166 TAG(e)==PYRAMID)
167 #else
168 #define NEWGREEN(e) (TAG(e)==HEXAHEDRON || TAG(e)== PRISM || \
169 TAG(e)==PYRAMID || TAG(e)== TETRAHEDRON)
170 #endif
171
172 /* marked elem with new green refinement (without rule, only 3D) */
173 #ifdef __ANISOTROPIC__
174 #define MARKED_NEW_GREEN(elem) \
175 (DIM==3 && ((NEWGREEN(elem) && MARKCLASS(elem)==GREEN_CLASS) || \
176 ((TAG(elem)==PRISM && MARKCLASS(elem)==RED_CLASS && USED(elem)==1))))
177 #else
178 #define MARKED_NEW_GREEN(elem) \
179 (DIM==3 && NEWGREEN(elem) && MARKCLASS(elem)==GREEN_CLASS)
180 #endif
181
182 /** \brief Refined elem with new green refinement (without rule, only 3D) */
183 #define REFINED_NEW_GREEN(elem) \
184 (DIM==3 && NEWGREEN(elem) && REFINECLASS(elem)==GREEN_CLASS)
185
186 /** \brief Macro to test whether element changes refinement */
187 #define REFINEMENT_CHANGES(elem) \
188 (REF_TYPE_CHANGES(elem) || \
189 (MARKED_NEW_GREEN(elem) && \
190 (REFINECLASS(elem)!=GREEN_CLASS || (REFINECLASS(elem)==GREEN_CLASS \
191 && USED(elem)==1))))
192
193 /** @name Macros for storing sparse data needed in ExchangeClosureInfo() */
194 /*@{*/
195 #define MARKCLASSDATA_SHIFT 20
196 #define GETMARKCLASSDATA(elem,dataadr) \
197 (*(dataadr)) = (*(dataadr)) | ((MARKCLASS(elem))<<MARKCLASSDATA_SHIFT)
198 #define SETMARKCLASSDATA(elem,data) \
199 SETMARKCLASS((elem),((data)>>MARKCLASSDATA_SHIFT)&((1<<MARKCLASS_LEN)-1))
200 /*@}*/
201
202 /** @name Macros for storing sparse data needed in ExchangeClosureInfo() */
203 /*@{*/
204 #define MARKDATA_SHIFT 22
205 #define GETMARKDATA(elem,dataadr) \
206 (*(dataadr)) = (*(dataadr)) | ((MARK(elem))<<MARKDATA_SHIFT)
207 #define SETMARKDATA(elem,data) \
208 SETMARK((elem),MARK(elem)|(((data)>>MARKDATA_SHIFT)&((1<<MARK_LEN)-1)))
209 /*@}*/
210
211 /** @name Macros for storing sparse data needed in ExchangeClosureInfo() */
212 /*@{*/
213 #define COARSENDATA_SHIFT 19
214 #define GETCOARSENDATA(elem,dataadr) \
215 (*(dataadr)) = (*(dataadr)) | ((COARSEN(elem))<<COARSENDATA_SHIFT)
216 #define SETCOARSENDATA(elem,data) \
217 SETCOARSEN((elem),((data)>>COARSENDATA_SHIFT)&((1<<COARSEN_LEN)-1))
218 /*@}*/
219
220 /** @name Macros to get and set the edgepattern on an element */
221 /*@{*/
222 #define GetEdgeInfo(elem,patadr,macro) \
223 { \
224 INT i,pattern; \
225 EDGE *theEdge; \
226 \
227 pattern = 0; \
228 for (i=EDGES_OF_ELEM(elem)-1; i>=0; i--) \
229 { \
230 theEdge=GetEdge(CORNER_OF_EDGE_PTR((elem),i,0), \
231 CORNER_OF_EDGE_PTR((elem),i,1)); \
232 ASSERT(theEdge!=NULL); \
233 \
234 pattern = (pattern<<1) | macro(theEdge); \
235 } \
236 \
237 (*(patadr)) |= pattern; \
238 }
239
240 #define SetEdgeInfo(elem,pat,macro,OP) \
241 { \
242 INT i,pattern; \
243 EDGE *theEdge; \
244 \
245 pattern = (pat); \
246 for (i=0; i<EDGES_OF_ELEM(elem); i++) \
247 { \
248 theEdge = GetEdge(CORNER_OF_EDGE_PTR((elem),i,0), \
249 CORNER_OF_EDGE_PTR((elem),i,1)); \
250 ASSERT(theEdge!=NULL); \
251 \
252 SET ## macro(theEdge, macro(theEdge) OP (pattern&0x1)); \
253 pattern >>= 1; \
254 } \
255 }
256 /*@}*/
257
258 #define REFINE_CONTEXT_LIST(d,context) \
259 IFDEBUG(gm,2) \
260 { \
261 INT i; \
262 \
263 UserWrite(" UpdateContext is :\n"); \
264 for(i=0; i<MAX_CORNERS_OF_ELEM+MAX_NEW_CORNERS_DIM; i++) \
265 UserWriteF(" %3d",i); \
266 UserWrite("\n"); \
267 for(i=0; i<MAX_CORNERS_OF_ELEM+MAX_NEW_CORNERS_DIM; i++) \
268 if ((context)[i] != NULL) \
269 UserWriteF(" %3d",ID((context)[i])); \
270 else \
271 UserWriteF(" "); \
272 UserWrite("\n"); \
273 } \
274 ENDDEBUG
275
276 #ifndef STAT_OUT
277 #undef NEW_TIMER
278 #undef DEL_TIMER
279 #undef START_TIMER
280 #undef STOP_TIMER
281 #undef DIFF_TIMER
282 #undef SUM_TIMER
283 #undef EVAL_TIMER
284
285 #define NEW_TIMER(n)
286 #define DEL_TIMER(n)
287 #define START_TIMER(n)
288 #define STOP_TIMER(n)
289 #define DIFF_TIMER(n)
290 #define SUM_TIMER(n)
291 #define EVAL_TIMER(n)
292 #endif
293
294 /****************************************************************************/
295 /* */
296 /* data structures used in this source file (exported data structures are */
297 /* in the corresponding include file!) */
298 /* */
299 /****************************************************************************/
300
301 typedef NODE *ELEMENTCONTEXT[MAX_CORNERS_OF_ELEM+MAX_NEW_CORNERS_DIM];
302
303
304 /****************************************************************************/
305 /* */
306 /* definition of exported global variables */
307 /* */
308 /****************************************************************************/
309
310 /* information used by the estimator and refine*/
311 REFINEINFO NS_DIM_PREFIX refine_info;
312
313 #ifdef ModelP
314 /* control words for identiftication of new nodes and edges */
315 INT NS_DIM_PREFIX ce_NEW_NIDENT;
316 INT NS_DIM_PREFIX ce_NEW_EDIDENT;
317 #endif
318
319 /****************************************************************************/
320 /* */
321 /* definition of variables global to this source file only (static!) */
322 /* */
323 /****************************************************************************/
324
325 /** \brief type of refinement */
326 static int rFlag=GM_REFINE_TRULY_LOCAL;
327
328 /** \brief refine with hanging nodes? */
329 static int hFlag=0;
330
331 /** \brief use fifo? 0=no 1=yes */
332 static int fifoFlag=0;
333
334 /** \brief fifo loop counter */
335 static int first;
336
337 /** \brief first element in fifo work list */
338 static ELEMENT *fifo_first=NULL;
339
340 /** \brief last element in fifo work list */
341 static ELEMENT *fifo_last=NULL;
342
343 /** \brief first element in fifo insertlist */
344 static ELEMENT *fifo_insertfirst=NULL;
345
346 /** \brief last element in fifo insertlist */
347 static ELEMENT *fifo_insertlast=NULL;
348
349 /** \brief first element to consider for next loop */
350 static ELEMENT *firstElement=NULL;
351
352 /** \brief Counter for green refinements doesn't need to be updated */
353 static INT No_Green_Update;
354
355 /** \brief green refined element counter */
356 static INT Green_Marks;
357
358 /** \brief 0/1: do/do not parallel part */
359 static INT refine_seq = 0;
360
361 /** \brief counter for FIFO loops */
362 static INT fifoloop = 0;
363
364 /** \brief count of adapted elements */
365 static INT total_adapted = 0;
366
367 #ifdef STAT_OUT
368 /* timing variables */
369 static int adapt_timer,closure_timer,gridadapt_timer,gridadapti_timer;
370 static int gridadaptl_timer,ident_timer,overlap_timer,gridcons_timer;
371 static int algebra_timer;
372 #endif
373
374 #ifdef DUNE_UGGRID_TET_RULESET
375 /* determine number of edge from reduced (i.e. restricted to one side) edgepattern */
376 /* if there are two edges marked for bisection, if not deliver -1. If the edge- */
377 /* is not reduced (i.e. marked edges lying on more than one side) deliver -2 */
378 static INT TriSectionEdge[64][2] =
379 { {-1,-1},{-1,-1},{-1,-1},{ 1, 0},{-1,-1},{ 0, 2},{ 2, 1},{-1,-1},
380 {-1,-1},{ 3, 0},{-2,-2},{-2,-2},{ 2, 3},{-2,-2},{-2,-2},{-2,-2},
381 {-1,-1},{ 0, 4},{ 4, 1},{-2,-2},{-2,-2},{-2,-2},{-2,-2},{-2,-2},
382 { 4, 3},{-1,-1},{-2,-2},{-2,-2},{-2,-2},{-2,-2},{-2,-2},{-2,-2},
383 {-1,-1},{-2,-2},{ 1, 5},{-2,-2},{ 5, 2},{-2,-2},{-2,-2},{-2,-2},
384 { 3, 5},{-2,-2},{-2,-2},{-2,-2},{-1,-1},{-2,-2},{-2,-2},{-2,-2},
385 { 5, 4},{-2,-2},{-1,-1},{-2,-2},{-2,-2},{-2,-2},{-2,-2},{-2,-2},
386 {-2,-2},{-2,-2},{-2,-2},{-2,-2},{-2,-2},{-2,-2},{-2,-2},{-2,-2}};
387
388 /* the indices of the edges of each side */
389 static INT CondensedEdgeOfSide[4] = {0x07,0x32,0x2C,0x19};
390 #endif
391
392 /****************************************************************************/
393 /* */
394 /* forward declarations of functions used before they are defined */
395 /* */
396 /****************************************************************************/
397
398 #ifdef ModelP
399 static void CheckConsistency (MULTIGRID *theMG, INT level ,INT debugstart, INT gmlevel, int *check);
400 #endif
401
402 /****************************************************************************/
403 /** \brief Fill refineinfo structure
404
405 This function fills the refineinfo structure
406
407 \return <ul>
408 <li> GM_OK if ok </li>
409 <li> GM_ERROR if an error occurs </li>
410 </ul>
411 */
412 /****************************************************************************/
413
SetRefineInfo(MULTIGRID * theMG)414 INT NS_DIM_PREFIX SetRefineInfo (MULTIGRID *theMG)
415 {
416 if (MultiGridStatus(theMG,1,0,0,0) != GM_OK) return(GM_ERROR);
417
418 return(GM_OK);
419 }
420
421
422 /****************************************************************************/
423 /** \brief Test entries of refineinfo structure
424
425 This function tests entries of refineinfo structure
426
427 \return <ul>
428 .n GM_OK if MG can be refined
429 .n GM_ERROR if MG refinement will lead to heap overflow
430 */
431 /****************************************************************************/
432
TestRefineInfo(MULTIGRID * theMG)433 INT NS_DIM_PREFIX TestRefineInfo (MULTIGRID *theMG)
434 {
435 if (PREDNEW0(REFINEINFO(theMG)) > PREDMAX(REFINEINFO(theMG)))
436 return(GM_ERROR);
437 else
438 return(GM_OK);
439 }
440
441 /****************************************************************************/
442 /** \brief Drop marks from leafelements to first regular
443
444 This function drops marks from leafelements to first regular, and resets
445 marks on all elements above (important for restrict marks)
446
447 \return <ul>
448 .n GM_OK if ok
449 .n GM_ERROR if an error occurs
450 */
451 /****************************************************************************/
452
DropMarks(MULTIGRID * theMG)453 static INT DropMarks (MULTIGRID *theMG)
454 {
455 INT k, Mark;
456 GRID *theGrid;
457 ELEMENT *theElement, *FatherElement;
458
459 return(GM_OK);
460
461 for (k=TOPLEVEL(theMG); k>0; k--)
462 {
463 theGrid = GRID_ON_LEVEL(theMG,k);
464 for (theElement=FIRSTELEMENT(theGrid); theElement!=NULL;
465 theElement=SUCCE(theElement))
466 if ((MARKCLASS(theElement) == RED_CLASS) &&
467 (ECLASS(theElement) != RED_CLASS))
468 {
469 Mark = MARK(theElement);
470 /** \todo marks must be changed if element type changes */
471 if (TAG(theElement)!=HEXAHEDRON &&
472 TAG(EFATHER(theElement))==HEXAHEDRON) Mark = HEXA_RED;
473 if (TAG(theElement)!=PYRAMID &&
474 TAG(EFATHER(theElement))==PYRAMID) Mark = PYR_RED;
475 FatherElement = theElement;
476
477 SETMARK(FatherElement,NO_REFINEMENT);
478 SETMARKCLASS(FatherElement,NO_CLASS);
479 FatherElement = EFATHER(FatherElement);
480
481 SETMARK(FatherElement,Mark);
482 SETMARKCLASS(FatherElement,RED_CLASS);
483 }
484 }
485 return(GM_OK);
486 }
487
488 /* Functions for realizing the (parallel) closure FIFO */
489
490
491 /****************************************************************************/
492 /** \brief Function for realizing the (parallel) closure FIFO
493 */
494 /****************************************************************************/
495
InitClosureFIFO(void)496 static INT InitClosureFIFO (void)
497 {
498 fifo_first=fifo_last=fifo_insertfirst=fifo_insertlast=NULL;
499 first = 1;
500 fifoloop = 0;
501 if (0) UserWriteF("Using FIFO: loop %d\n",fifoloop);
502
503 return (GM_OK);
504 }
505
506
507 /****************************************************************************/
508 /** \brief Function for realizing the (parallel) closure FIFO
509
510 \param theGrid - pointer to grid structure
511 . theElement
512 . thePattern
513 . NewPattern
514
515 DESCRIPTION:
516
517 */
518 /****************************************************************************/
519
UpdateFIFOLists(GRID * theGrid,ELEMENT * theElement,INT thePattern,INT NewPattern)520 static INT UpdateFIFOLists (GRID *theGrid, ELEMENT *theElement, INT thePattern, INT NewPattern)
521 {
522 #ifdef __TWODIM__
523 INT j;
524 EDGE *theEdge;
525 ELEMENT *NbElement;
526 #endif
527
528 if (MARKCLASS(theElement)==RED_CLASS && thePattern!=NewPattern)
529 {
530 #ifdef __TWODIM__
531 for (j=0; j<EDGES_OF_ELEM(theElement); j++)
532 {
533 if (EDGE_IN_PAT(thePattern,j)==0 &&
534 EDGE_IN_PAT(NewPattern,j))
535 {
536
537 theEdge=GetEdge(CORNER_OF_EDGE_PTR(theElement,j,0),
538 CORNER_OF_EDGE_PTR(theElement,j,1));
539 ASSERT(theEdge != NULL);
540
541 SETPATTERN(theEdge,1);
542
543 /* boundary case */
544 if (SIDE_ON_BND(theElement,j)) continue;
545
546 /* add the element sharing this edge to fifo_queue */
547 NbElement = NBELEM(theElement,j);
548
549 if (NbElement==NULL) continue;
550
551 PRINTDEBUG(gm,1,(" ADDING to FIFO: NBID=%d\n",
552 ID(NbElement)))
553
554 /* unlink element from element list */
555 if (PREDE(NbElement) != NULL)
556 SUCCE(PREDE(NbElement)) = SUCCE(NbElement);
557 if (SUCCE(NbElement) != NULL)
558 PREDE(SUCCE(NbElement)) = PREDE(NbElement);
559 if (FIRSTELEMENT(theGrid) == NbElement)
560 FIRSTELEMENT(theGrid) = SUCCE(NbElement);
561
562 SUCCE(NbElement) = PREDE(NbElement) = NULL;
563 /* insert into fifo */
564 if (fifo_insertfirst == NULL)
565 {
566 fifo_insertfirst = fifo_insertlast = NbElement;
567 }
568 else
569 {
570 SUCCE(fifo_insertlast) = NbElement;
571 PREDE(NbElement) = fifo_insertlast;
572 fifo_insertlast = NbElement;
573 }
574 }
575
576 if (EDGE_IN_PAT(thePattern,j) &&
577 EDGE_IN_PAT(NewPattern,j)==0)
578 {
579
580 UserWriteF("UpdateFIFOLists(): ERROR EID=%d in fifo "
581 "thePattern=%d has edge=%d refined but "
582 "NewPattern=%d NOT!\n",
583 ID(theElement),thePattern,j,NewPattern);
584 RETURN(-1);
585 }
586 }
587 #endif
588 #ifdef __THREEDIM__
589 UserWriteF("UpdateFIFOLists(): ERROR fifo for 3D NOT implemented!\n");
590 ASSERT(0);
591 #endif
592 }
593
594 return(GM_OK);
595 }
596
597
598 /****************************************************************************/
599 /** \brief Function for realizing the (parallel) closure FIFO
600
601 \param theGrid - pointer to grid structure
602
603 */
604 /****************************************************************************/
605
UpdateClosureFIFO(GRID * theGrid)606 static INT UpdateClosureFIFO (GRID *theGrid)
607 {
608 ELEMENT *theElement;
609
610 /* insert fifo work list into elementlist */
611 for (theElement=fifo_last; theElement!=NULL;
612 theElement=PREDE(theElement))
613 {
614 SUCCE(theElement) = FIRSTELEMENT(theGrid);
615 PREDE(FIRSTELEMENT(theGrid)) = theElement;
616 FIRSTELEMENT(theGrid) = theElement;
617 }
618
619 PREDE(FIRSTELEMENT(theGrid)) = NULL;
620
621 if (fifo_insertfirst != NULL)
622 {
623 /* append fifo insert list to fifo work list */
624 firstElement = fifo_first = fifo_insertfirst;
625 fifo_last = fifo_insertlast;
626
627 IFDEBUG(gm,2)
628 UserWriteF(" FIFO Queue:");
629 for (theElement=fifo_first; theElement!=NULL;
630 theElement=SUCCE(theElement))
631 UserWriteF(" %d\n", ID(theElement));
632 ENDDEBUG
633
634 fifo_insertfirst = fifo_insertlast = NULL;
635 first = 0;
636 fifoloop++;
637 UserWriteF(" loop %d",fifoloop);
638 return(1);
639 }
640
641 return(0);
642 }
643
644
645 /****************************************************************************/
646 /*
647 ManageParallelFIFO - function for realizing the (parallel) closure FIFO
648
649 SYNOPSIS:
650 static INT ManageParallelFIFO (ELEMENT *firstElement);
651
652 PARAMETERS:
653 \param firstElement
654
655 DESCRIPTION:
656
657 \return <ul>
658 INT
659 */
660 /****************************************************************************/
661
ManageParallelFIFO(const PPIF::PPIFContext & context,ELEMENT * firstElement)662 static INT ManageParallelFIFO (const PPIF::PPIFContext& context, ELEMENT *firstElement)
663 {
664 #if defined(FIFO) && defined(ModelP)
665 ELEMENT *theElement;
666
667 if (context.procs() == 1) return(0);
668
669 do
670 {
671 /* exchange FIFO flag and PATTERN from slaves to master */
672 IF_FIF0AndPat_S2M(GLEVEL(grid));
673
674 /* add all master elements of horizontal interface to FIFO */
675 for (theElement=firstElement; theElement!=NULL;
676 theElement=SUCCE(theElement))
677 {
678 if (IS_HOR_MASTER(theElement) && FIFO(theElement))
679 {
680 AddToFIFIO(theElement);
681 SETFIFO(theElement,0);
682 }
683 }
684
685 /* check condition for termination of pattern adaptation */
686 AllFIFOsEmpty = CheckGlobalFIFOStatus(fifo);
687 }
688 while (fifo==NULL && AllFIFOsEmpty==1);
689 #else
690 return (0);
691 #endif
692 }
693
694 /****************************************************************************/
695 /** \brief Changes refinement of element
696
697 \param theElement - pointer to element
698
699 This function returns a boolean value to indicate, whether an element
700 changes its refinement (1) or not (0).
701
702 \return <ul>
703 .n 0 if element will not change refinement
704 .n 1 if it will
705 */
706 /****************************************************************************/
707
Refinement_Changes(ELEMENT * theElement)708 INT NS_DIM_PREFIX Refinement_Changes (ELEMENT *theElement)
709 {
710 return(REFINEMENT_CHANGES(theElement));
711 }
712
713 /****************************************************************************/
714 /*
715 GridClosure - compute closure for next level
716
717 SYNOPSIS:
718 static INT PrepareGridClosure (GRID *theGrid);
719
720 PARAMETERS:
721 \param theGrid - pointer to grid structure
722
723 DESCRIPTION:
724 This function computes the closure for next level. A closure can only be
725 determined if the rule set for the used elements is complete. This means
726 that for all side and edge patterns possible for an element type exists a
727 rule which closes the element. In this case a FIFO for computing the closure
728 is not needed any more and the closure can be computed in one step.
729
730 \return <ul>
731 INT
732 .n >0 if elements will be refined
733 .n 0 if no elements will be refined
734 .n -1 if an error occured
735 */
736 /****************************************************************************/
737
PrepareGridClosure(GRID * theGrid)738 static INT PrepareGridClosure (GRID *theGrid)
739 {
740 INT j;
741 ELEMENT *theElement;
742 EDGE *theEdge;
743
744 /* reset USED flag of elements and PATTERN and */
745 /* ADDPATTERN flag on the edges */
746 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL;
747 theElement=SUCCE(theElement))
748 {
749 SETUSED(theElement,0);
750 if (EGHOST(theElement))
751 {
752 SETCOARSEN(theElement,0);
753 SETMARK(theElement,NO_REFINEMENT);
754 SETMARKCLASS(theElement,0);
755 }
756
757 for (j=0; j<EDGES_OF_ELEM(theElement); j++)
758 {
759 theEdge=GetEdge(CORNER_OF_EDGE_PTR(theElement,j,0),
760 CORNER_OF_EDGE_PTR(theElement,j,1));
761 ASSERT(theEdge != NULL);
762
763 SETPATTERN(theEdge,0);
764 SETADDPATTERN(theEdge,1); /* needed in RestrictMarks() */
765 }
766 }
767
768 return(GM_OK);
769 }
770
771 #ifdef ModelP
772
773
774 /****************************************************************************/
775 /*
776 Gather_ElementClosureInfo -
777
778 SYNOPSIS:
779 static int Gather_ElementClosureInfo (DDD_OBJ obj, void *data);
780
781 PARAMETERS:
782 . obj
783 . data
784
785 DESCRIPTION:
786
787 \return <ul>
788 int
789 */
790 /****************************************************************************/
791
Gather_ElementClosureInfo(DDD::DDDContext &,DDD_OBJ obj,void * data,DDD_PROC proc,DDD_PRIO prio)792 static int Gather_ElementClosureInfo (DDD::DDDContext&, DDD_OBJ obj, void *data, DDD_PROC proc, DDD_PRIO prio)
793 {
794 INT refinedata;
795 ELEMENT *theElement = (ELEMENT *)obj;
796
797 PRINTDEBUG(gm,1,("Gather_ElementClosureInfo(): e=" EID_FMTX "\n",
798 EID_PRTX(theElement)))
799
800 refinedata = 0;
801
802 #ifdef __TWODIM__
803 GetEdgeInfo(theElement,&refinedata,PATTERN);
804 #endif
805
806 /* mark and sidepattern have same control word positions */
807 /* if this changes sidepattern must be sent separately */
808 GETMARKDATA(theElement,&refinedata);
809 GETMARKCLASSDATA(theElement,&refinedata);
810 GETCOARSENDATA(theElement,&refinedata);
811 ((INT *)data)[0] = refinedata;
812
813 PRINTDEBUG(gm,1,("Gather_ElementClosureInfo(): refinedata=%08x "
814 "sidepattern=%d markclass=%d mark=%d coarse=%d\n",refinedata,
815 SIDEPATTERN(theElement),MARKCLASS(theElement),MARK(theElement),COARSEN(theElement)))
816
817 return(GM_OK);
818 }
819
820
821
822 /****************************************************************************/
823 /*
824 Scatter_ElementClosureInfo -
825
826 SYNOPSIS:
827 static int Scatter_ElementClosureInfo (DDD_OBJ obj, void *data);
828
829 PARAMETERS:
830 . obj
831 . data
832
833 DESCRIPTION:
834
835 \return <ul>
836 int
837 */
838 /****************************************************************************/
839
Scatter_ElementClosureInfo(DDD::DDDContext &,DDD_OBJ obj,void * data,DDD_PROC proc,DDD_PRIO prio)840 static int Scatter_ElementClosureInfo (DDD::DDDContext&, DDD_OBJ obj, void *data, DDD_PROC proc, DDD_PRIO prio)
841 {
842 INT refinedata;
843 ELEMENT *theElement = (ELEMENT *)obj;
844
845 PRINTDEBUG(gm,1,("Scatter_ElementClosureInfo(): e=" EID_FMTX "\n",
846 EID_PRTX(theElement)))
847
848 refinedata = ((INT *)data)[0];
849
850 #ifdef __TWODIM__
851 SetEdgeInfo(theElement,refinedata,PATTERN,|);
852 #endif
853
854 /* mark and sidepattern have same control word positions */
855 /* if this changes sidepattern must be sent separately */
856 SETMARKDATA(theElement,refinedata);
857
858 if (EMASTER(theElement)) return(GM_OK);
859 if (EGHOST(theElement) && EGHOSTPRIO(prio)) return(GM_OK);
860
861 SETMARKCLASSDATA(theElement,refinedata);
862 SETCOARSENDATA(theElement,refinedata);
863
864 PRINTDEBUG(gm,1,("Scatter_ElementClosureInfo(): refinedata=%08x "
865 "sidepattern=%d markclass=%d mark=%d coarse=%d\n",refinedata,
866 SIDEPATTERN(theElement),MARKCLASS(theElement),MARK(theElement),COARSEN(theElement)))
867
868 return(GM_OK);
869 }
870
ExchangeElementClosureInfo(GRID * theGrid)871 static INT ExchangeElementClosureInfo (GRID *theGrid)
872 {
873 auto& context = theGrid->dddContext();
874 const auto& dddctrl = ddd_ctrl(context);
875
876 /* exchange information of elements to compute closure */
877 DDD_IFAOnewayX(context,
878 dddctrl.ElementSymmVHIF,GRID_ATTR(theGrid),IF_FORWARD,sizeof(INT),
879 Gather_ElementClosureInfo, Scatter_ElementClosureInfo);
880
881 return(GM_OK);
882 }
883
Gather_ElementRefine(DDD::DDDContext &,DDD_OBJ obj,void * data,DDD_PROC proc,DDD_PRIO prio)884 static int Gather_ElementRefine (DDD::DDDContext&, DDD_OBJ obj, void *data, DDD_PROC proc, DDD_PRIO prio)
885 {
886 ELEMENT *theElement = (ELEMENT *)obj;
887
888 PRINTDEBUG(gm,1,("Gather_ElementRefine(): e=" EID_FMTX "\n",
889 EID_PRTX(theElement)))
890
891 ((INT *)data)[0] = MARKCLASS(theElement);
892 ((INT *)data)[1] = MARK(theElement);
893
894 return(GM_OK);
895 }
896
Scatter_ElementRefine(DDD::DDDContext &,DDD_OBJ obj,void * data,DDD_PROC proc,DDD_PRIO prio)897 static int Scatter_ElementRefine (DDD::DDDContext&, DDD_OBJ obj, void *data, DDD_PROC proc, DDD_PRIO prio)
898 {
899 ELEMENT *theElement = (ELEMENT *)obj;
900
901 PRINTDEBUG(gm,1,("Scatter_ElementClosureInfo(): e=" EID_FMTX "\n",
902 EID_PRTX(theElement)))
903
904 if (EMASTER(theElement)) return(GM_OK);
905 if (EGHOST(theElement) && EGHOSTPRIO(prio)) return(GM_OK);
906
907 SETMARKCLASS(theElement,((INT *)data)[0]);
908 SETMARK(theElement,((INT *)data)[1]);
909
910 return(GM_OK);
911 }
912
ExchangeElementRefine(GRID * theGrid)913 static INT ExchangeElementRefine (GRID *theGrid)
914 {
915 auto& context = theGrid->dddContext();
916 const auto& dddctrl = ddd_ctrl(context);
917
918 /* exchange information of elements to compute closure */
919 DDD_IFAOnewayX(context,
920 dddctrl.ElementSymmVHIF,GRID_ATTR(theGrid),IF_FORWARD,2*sizeof(INT),
921 Gather_ElementRefine, Scatter_ElementRefine);
922
923 return(GM_OK);
924 }
925
926 #ifdef __THREEDIM__
927
928 /****************************************************************************/
929 /*
930 Gather_EdgeClosureInfo -
931
932 SYNOPSIS:
933 static int Gather_EdgeClosureInfo (DDD_OBJ obj, void *data);
934
935 PARAMETERS:
936 . obj
937 . data
938
939 DESCRIPTION:
940
941 \return <ul>
942 int
943 */
944 /****************************************************************************/
945
Gather_EdgeClosureInfo(DDD::DDDContext &,DDD_OBJ obj,void * data)946 static int Gather_EdgeClosureInfo (DDD::DDDContext&, DDD_OBJ obj, void *data)
947 {
948 INT pattern;
949 EDGE *theEdge = (EDGE *)obj;
950
951 PRINTDEBUG(gm,1,("Gather_EdgeClosureInfo(): e=" ID_FMTX "pattern=%d \n",
952 ID_PRTX(theEdge),PATTERN(theEdge)))
953
954 pattern = PATTERN(theEdge);
955
956 ((INT *)data)[0] = pattern;
957
958 return(GM_OK);
959 }
960
961
962
963 /****************************************************************************/
964 /*
965 Scatter_EdgeClosureInfo -
966
967 SYNOPSIS:
968 static int Scatter_EdgeClosureInfo (DDD_OBJ obj, void *data);
969
970 PARAMETERS:
971 . obj
972 . data
973
974 DESCRIPTION:
975
976 \return <ul>
977 int
978 */
979 /****************************************************************************/
980
Scatter_EdgeClosureInfo(DDD::DDDContext &,DDD_OBJ obj,void * data)981 static int Scatter_EdgeClosureInfo (DDD::DDDContext&, DDD_OBJ obj, void *data)
982 {
983 INT pattern;
984 EDGE *theEdge = (EDGE *)obj;
985
986 pattern = MAX(PATTERN(theEdge),((INT *)data)[0]);
987
988 PRINTDEBUG(gm,1,("Gather_EdgeClosureInfo(): e=" ID_FMTX "pattern=%d \n",
989 ID_PRTX(theEdge),pattern))
990
991 SETPATTERN(theEdge,pattern);
992
993 return(GM_OK);
994 }
995
ExchangeEdgeClosureInfo(GRID * theGrid)996 INT ExchangeEdgeClosureInfo (GRID *theGrid)
997 {
998 auto& context = theGrid->dddContext();
999 const auto& dddctrl = ddd_ctrl(context);
1000
1001 /* exchange information of edges to compute closure */
1002 DDD_IFAOneway(context,
1003 dddctrl.EdgeVHIF,GRID_ATTR(theGrid),IF_FORWARD,sizeof(INT),
1004 Gather_EdgeClosureInfo, Scatter_EdgeClosureInfo);
1005
1006 return(GM_OK);
1007 }
1008 #endif
1009
1010 /****************************************************************************/
1011 /*
1012 ExchangeClosureInfo -
1013
1014 SYNOPSIS:
1015 static INT ExchangeClosureInfo (GRID *theGrid);
1016
1017 PARAMETERS:
1018 . theGrid
1019
1020 DESCRIPTION:
1021
1022 \return <ul>
1023 INT
1024 */
1025 /****************************************************************************/
1026
ExchangeClosureInfo(GRID * theGrid)1027 static INT ExchangeClosureInfo (GRID *theGrid)
1028 {
1029
1030 /* exchange information of elements to compute closure */
1031 if (ExchangeElementClosureInfo(theGrid) != GM_OK) RETURN(GM_ERROR);
1032
1033 #ifdef __THREEDIM__
1034 /* exchange information of edges to compute closure */
1035 if (ExchangeEdgeClosureInfo(theGrid) != GM_OK) RETURN(GM_ERROR);
1036 #endif
1037
1038 return(GM_OK);
1039 }
1040 #endif
1041
1042
1043 /****************************************************************************/
1044 /*
1045 ComputePatterns -
1046
1047 SYNOPSIS:
1048 static INT ComputePatterns (GRID *theGrid);
1049
1050 PARAMETERS:
1051 . theGrid - pointer to grid structure
1052
1053 DESCRIPTION:
1054
1055 \return <ul>
1056 INT
1057 */
1058 /****************************************************************************/
1059
ComputePatterns(GRID * theGrid)1060 static INT ComputePatterns (GRID *theGrid)
1061 {
1062 SHORT *thePattern;
1063 INT i,Mark;
1064 ELEMENT *theElement;
1065 EDGE *theEdge;
1066
1067 /* ComputePatterns works only on master elements */
1068 /* since ghost elements have no information from */
1069 /* RestrictMarks() up to this time and this may */
1070 /* leed to inconsistency while coarsening */
1071
1072 /* reset EDGE/SIDEPATTERN in elements */
1073 /* set SIDEPATTERN in elements */
1074 /* set PATTERN on the edges */
1075 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL;
1076 theElement=SUCCE(theElement))
1077 {
1078 #ifdef ModelP
1079 if (EGHOST(theElement))
1080 {
1081 #ifdef __THREEDIM__
1082 SETSIDEPATTERN(theElement,0);
1083 #endif
1084 continue;
1085 }
1086 #endif
1087
1088 if (MARKCLASS(theElement)==RED_CLASS)
1089 {
1090 Mark = MARK(theElement);
1091 thePattern = MARK2PATTERN(theElement,Mark);
1092
1093 for (i=0; i<EDGES_OF_ELEM(theElement); i++)
1094 if (EDGE_IN_PATTERN(thePattern,i))
1095 {
1096 theEdge=GetEdge(CORNER_OF_EDGE_PTR(theElement,i,0),
1097 CORNER_OF_EDGE_PTR(theElement,i,1));
1098
1099 ASSERT(theEdge != NULL);
1100
1101 SETPATTERN(theEdge,1);
1102 }
1103
1104 #ifdef __THREEDIM__
1105 /* SIDEPATTERN must be reset here for master elements, */
1106 /* because it overlaps with MARK (980217 s.l.) */
1107 SETSIDEPATTERN(theElement,0);
1108 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
1109 {
1110 #ifdef DUNE_UGGRID_TET_RULESET
1111 if (CORNERS_OF_SIDE(theElement,i)==4)
1112 {
1113 #endif
1114 /* set SIDEPATTERN if side has node */
1115 if(SIDE_IN_PATTERN(theElement,thePattern,i))
1116 SETSIDEPATTERN(theElement,
1117 SIDEPATTERN(theElement) | 1<<i);
1118 #ifdef DUNE_UGGRID_TET_RULESET
1119 }
1120 #endif
1121 }
1122 #endif
1123 }
1124 else
1125 {
1126 #ifdef __THREEDIM__
1127 /* SIDEPATTERN must be reset here for master elements, */
1128 /* because it overlaps with MARK (980217 s.l.) */
1129 SETSIDEPATTERN(theElement,0);
1130 #endif
1131 SETMARKCLASS(theElement,NO_CLASS);
1132 }
1133 }
1134
1135 return(GM_OK);
1136 }
1137
1138 #ifdef __THREEDIM__
1139 #ifdef DUNE_UGGRID_TET_RULESET
1140
1141 /****************************************************************************/
1142 /*
1143 CorrectTetrahedronSidePattern -
1144
1145 SYNOPSIS:
1146 static INT CorrectTetrahedronSidePattern (ELEMENT *theElement, INT i, ELEMENT *theNeighbor, INT j);
1147
1148 PARAMETERS:
1149 . theElement
1150 . i
1151 . theNeighbor
1152 . j
1153
1154 DESCRIPTION:
1155
1156 \return <ul>
1157 INT
1158 */
1159 /****************************************************************************/
1160
CorrectTetrahedronSidePattern(ELEMENT * theElement,INT i,ELEMENT * theNeighbor,INT j)1161 static INT CorrectTetrahedronSidePattern (ELEMENT *theElement, INT i, ELEMENT *theNeighbor, INT j)
1162 {
1163 INT k;
1164 INT theEdgeNum,theEdgePattern=0;
1165 INT NbEdgeNum,NbEdgePattern,NbSidePattern,NbSideMask;
1166 EDGE *theEdge,*NbEdge;
1167
1168 if (TAG(theElement)==PYRAMID || TAG(theElement)==PRISM)
1169 return(GM_OK);
1170
1171 for (k=EDGES_OF_ELEM(theElement)-1; k>=0; k--)
1172 {
1173 theEdge=GetEdge(CORNER_OF_EDGE_PTR(theElement,k,0),
1174 CORNER_OF_EDGE_PTR(theElement,k,1));
1175 ASSERT(theEdge!=NULL);
1176
1177 theEdgePattern = (theEdgePattern<<1) | PATTERN(theEdge);
1178 }
1179
1180 /* because SIDEPATTERN is set to zero, */
1181 /* I choose TriSectionEdge[0] */
1182 theEdgeNum = TriSectionEdge[theEdgePattern
1183 &CondensedEdgeOfSide[i]][0];
1184
1185 if (theEdgeNum == -2) RETURN(-1);
1186
1187 if (theEdgeNum == -1) return(GM_OK);
1188
1189 switch (TAG(theNeighbor))
1190 {
1191
1192 case TETRAHEDRON :
1193
1194 NbEdgePattern = 0;
1195
1196 for (k=0; k<EDGES_OF_ELEM(theNeighbor); k++)
1197 {
1198 NbEdge=GetEdge(CORNER_OF_EDGE_PTR(theNeighbor,k,0),
1199 CORNER_OF_EDGE_PTR(theNeighbor,k,1));
1200 ASSERT(NbEdge!=NULL);
1201 NbEdgePattern = NbEdgePattern | (PATTERN(NbEdge)<<k);
1202 }
1203
1204 NbEdgeNum = TriSectionEdge[NbEdgePattern
1205 &CondensedEdgeOfSide[j]][0];
1206
1207 if (NbEdgeNum == -2 || NbEdgeNum == -1)
1208 RETURN(-1);
1209
1210 if (!(CORNER_OF_EDGE_PTR(theElement,theEdgeNum,0) ==
1211 CORNER_OF_EDGE_PTR(theNeighbor,NbEdgeNum,0) &&
1212 CORNER_OF_EDGE_PTR(theElement,theEdgeNum,1) ==
1213 CORNER_OF_EDGE_PTR(theNeighbor,NbEdgeNum,1) )
1214 &&
1215 !(CORNER_OF_EDGE_PTR(theElement,theEdgeNum,0) ==
1216 CORNER_OF_EDGE_PTR(theNeighbor,NbEdgeNum,1) &&
1217 CORNER_OF_EDGE_PTR(theElement,theEdgeNum,1) ==
1218 CORNER_OF_EDGE_PTR(theNeighbor,NbEdgeNum,0) ) )
1219 {
1220 NbSidePattern = SIDEPATTERN(theNeighbor);
1221 NbSideMask = (1<<j);
1222
1223 if ( NbSidePattern & NbSideMask )
1224 {
1225 NbSidePattern &= ~NbSideMask;
1226 #ifdef ModelP
1227 /* in this case ExchangeSidePatterns() fails */
1228 /* does it occur ? (980217 s.l.) */
1229 assert(0);
1230 #endif
1231 }
1232 else
1233 NbSidePattern |= NbSideMask;
1234
1235 PRINTDEBUG(gm,1,("CorrectTetrahedronSidePattern(): nb=" EID_FMTX
1236 " new nbsidepattern=%d\n",EID_PRTX(theNeighbor),NbSidePattern));
1237 SETSIDEPATTERN(theNeighbor,NbSidePattern);
1238 }
1239 break;
1240
1241 case PYRAMID :
1242 case PRISM :
1243 {
1244 INT trisectionedge=-1;
1245
1246 for (k=0; k<CORNERS_OF_SIDE(theNeighbor,j); k++)
1247 {
1248 INT edge;
1249
1250 edge = EDGE_OF_SIDE(theElement,j,k);
1251
1252 NbEdge=GetEdge(CORNER_OF_EDGE_PTR(theNeighbor,edge,0),
1253 CORNER_OF_EDGE_PTR(theNeighbor,edge,1));
1254 ASSERT(NbEdge!=NULL);
1255
1256 if (PATTERN(NbEdge) && (edge>trisectionedge))
1257 trisectionedge = edge;
1258 }
1259 assert(trisectionedge != -1);
1260
1261 if (theEdgeNum != trisectionedge)
1262 SETSIDEPATTERN(theNeighbor,
1263 SIDEPATTERN(theNeighbor)|(1<<j));
1264
1265 break;
1266 }
1267
1268 default :
1269 ASSERT(0);
1270 }
1271
1272 return(GM_OK);
1273 }
1274 #endif
1275
1276
1277 /****************************************************************************/
1278 /*
1279 CorrectElementSidePattern -
1280
1281 SYNOPSIS:
1282 static INT CorrectElementSidePattern (ELEMENT *theElement, ELEMENT *theNeighbor, INT i);
1283
1284 PARAMETERS:
1285 . theElement
1286 . theNeighbor
1287 . i
1288
1289 DESCRIPTION:
1290
1291 \return <ul>
1292 INT
1293 */
1294 /****************************************************************************/
1295
CorrectElementSidePattern(ELEMENT * theElement,ELEMENT * theNeighbor,INT i)1296 static INT CorrectElementSidePattern (ELEMENT *theElement, ELEMENT *theNeighbor, INT i)
1297 {
1298 INT j;
1299
1300 #ifdef ModelP
1301 /* this case should not occur */
1302 if (theNeighbor == NULL)
1303 {
1304 ASSERT(EGHOST(theElement));
1305 UserWriteF("CorrectElementSidePattern(): error elem=" EID_FMTX " nb[%d]="
1306 EID_FMTX " nb=" EID_FMTX
1307 "\n", EID_PRTX(theElement),i,EID_PRTX(NBELEM(theElement,i)),
1308 EID_PRTX(theNeighbor));
1309 return(GM_OK);
1310 }
1311 #endif
1312
1313 /* search neighbors side */
1314 for (j=0; j<SIDES_OF_ELEM(theNeighbor); j++)
1315 if (NBELEM(theNeighbor,j) == theElement)
1316 break;
1317 #ifdef ModelP
1318 if (j >= SIDES_OF_ELEM(theNeighbor))
1319 {
1320 if (!(EGHOST(theElement) && EGHOST(theNeighbor)))
1321 {
1322 UserWriteF("CorrectElementSidePattern(): ERROR nbelem not found elem=%08x/"
1323 EID_FMTX " nb=%08x/" EID_FMTX "\n",
1324 theElement,EID_PRTX(theElement),theNeighbor,EID_PRTX(theNeighbor));
1325 }
1326 ASSERT(EGHOST(theElement) && EGHOST(theNeighbor));
1327 return(GM_OK);
1328 }
1329 #else
1330 ASSERT(j<SIDES_OF_ELEM(theNeighbor));
1331 #endif
1332
1333 /* side is triangle or quadrilateral */
1334 switch (CORNERS_OF_SIDE(theElement,i))
1335 {
1336 case 3 :
1337 #ifdef DUNE_UGGRID_TET_RULESET
1338 /* handle case with 2 edges of the side refined */
1339 if (CorrectTetrahedronSidePattern(theElement,i,theNeighbor,j) != GM_OK)
1340 RETURN(GM_ERROR);
1341 #endif
1342 break;
1343
1344 case 4 :
1345 /* if side of one of the neighboring elements has a */
1346 /* sidenode, then both need a sidenode */
1347 if (SIDE_IN_PAT(SIDEPATTERN(theElement),i))
1348 {
1349 SETSIDEPATTERN(theNeighbor,
1350 SIDEPATTERN(theNeighbor) | (1<<j));
1351 }
1352 else if (SIDE_IN_PAT(SIDEPATTERN(theNeighbor),j))
1353 {
1354 SETSIDEPATTERN(theElement,
1355 SIDEPATTERN(theElement) | (1<<i));
1356 }
1357 break;
1358
1359 default :
1360 ASSERT(0);
1361 }
1362
1363 return(GM_OK);
1364 }
1365
1366
1367 /****************************************************************************/
1368 /*
1369 SetElementSidePatterns -
1370
1371 SYNOPSIS:
1372 static INT SetElementSidePatterns (GRID *theGrid, ELEMENT *firstElement);
1373
1374 PARAMETERS:
1375 . theGrid - pointer to grid structure
1376 . firstElement
1377
1378 DESCRIPTION:
1379
1380 \return <ul>
1381 INT
1382 */
1383 /****************************************************************************/
1384
SetElementSidePatterns(GRID * theGrid,ELEMENT * firstElement)1385 static INT SetElementSidePatterns (GRID *theGrid, ELEMENT *firstElement)
1386 {
1387 INT i;
1388 ELEMENT *theElement,*theNeighbor;
1389
1390 /* set pattern (edge and side) on the elements */
1391 for (theElement=firstElement; theElement!=NULL;
1392 theElement=SUCCE(theElement))
1393 {
1394 /* make edgepattern consistent with pattern of edges */
1395 SETUSED(theElement,1);
1396
1397 #ifndef __ANISOTROPIC__
1398 /** \todo change this for red refinement of pyramids */
1399 if (DIM==3 && TAG(theElement)==PYRAMID) continue;
1400 #endif
1401
1402 /* make sidepattern consistent with neighbors */
1403 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
1404 {
1405 theNeighbor = NBELEM(theElement,i);
1406 if (theNeighbor == NULL) continue;
1407
1408 /* only one of the neighboring elements does corrections */
1409 /* determine element for side correction by (g)id */
1410 if (_EID_(theElement) < _EID_(theNeighbor)) continue;
1411
1412 /* determine element for side correction by used flag */
1413 /** \todo delete this:
1414 if (USED(theNeighbor)) continue;
1415 */
1416
1417 /* edgepatterns from theElement and theNeighbor are in final state */
1418
1419 if (CorrectElementSidePattern(theElement,theNeighbor,i) != GM_OK) RETURN(GM_ERROR);
1420 }
1421 }
1422
1423 return(GM_OK);
1424 }
1425 #endif
1426
1427
1428
1429 /****************************************************************************/
1430 /*
1431 SetElementRules -
1432
1433 SYNOPSIS:
1434 static INT SetElementRules (GRID *theGrid, ELEMENT *firstElement, INT *cnt);
1435
1436 PARAMETERS:
1437 . theGrid - pointer to grid structure
1438 . firstElement
1439 . cnt
1440
1441 DESCRIPTION:
1442
1443 \return <ul>
1444 INT
1445 */
1446 /****************************************************************************/
1447
SetElementRules(GRID * theGrid,ELEMENT * firstElement,INT * cnt)1448 static INT SetElementRules (GRID *theGrid, ELEMENT *firstElement, INT *cnt)
1449 {
1450 INT Mark,NewPattern;
1451 INT thePattern,theEdgePattern,theSidePattern=0;
1452 ELEMENT *theElement;
1453
1454 [[maybe_unused]] const int me = theGrid->ppifContext().me();
1455
1456 /* set refinement rules from edge- and sidepattern */
1457 (*cnt) = 0;
1458 for (theElement=firstElement; theElement!=NULL;
1459 theElement=SUCCE(theElement))
1460 {
1461 theEdgePattern = 0;
1462
1463 /* compute element pattern */
1464 GetEdgeInfo(theElement,&theEdgePattern,PATTERN);
1465
1466 #ifdef __TWODIM__
1467 thePattern = theEdgePattern;
1468 PRINTDEBUG(gm,2,(PFMT "SetElementRules(): e=" EID_FMTX " edgepattern=%d\n",
1469 me,EID_PRTX(theElement),theEdgePattern));
1470 #endif
1471 #ifdef __THREEDIM__
1472 theSidePattern = SIDEPATTERN(theElement);
1473 thePattern = theSidePattern<<EDGES_OF_ELEM(theElement) | theEdgePattern;
1474 PRINTDEBUG(gm,2,(PFMT "SetElementRules(): e=" EID_FMTX
1475 " edgepattern=%03x sidepattern=%02x\n",
1476 me,EID_PRTX(theElement),theEdgePattern,theSidePattern));
1477 #endif
1478
1479 /* get Mark from pattern */
1480 Mark = PATTERN2MARK(theElement,thePattern);
1481
1482 /* treat Mark according to mode */
1483 if (fifoFlag)
1484 {
1485 /* directed refinement */
1486 if (Mark == -1 && MARKCLASS(theElement)==RED_CLASS)
1487 {
1488 /* there is no rule for this pattern, switch to red */
1489 Mark = RED;
1490 }
1491 else
1492 ASSERT(Mark != -1);
1493 }
1494 else if (hFlag==0 && MARKCLASS(theElement)!=RED_CLASS)
1495 {
1496 /* refinement with hanging nodes */
1497 Mark = NO_REFINEMENT;
1498 }
1499 else
1500 {
1501 /* refinement with closure (default) */
1502 #ifdef __ANISOTROPIC__
1503 if (MARKCLASS(theElement)==RED_CLASS && TAG(theElement)==PRISM)
1504 {
1505 ASSERT(USED(theElement)==1);
1506 if (Mark==-1)
1507 {
1508 ASSERT(TAG(theElement)==PRISM);
1509 /* to implement the anisotropic case for other elements */
1510 /* and anisotropic refinements the initial anisotropic */
1511 /* rule is needed here. (980316 s.l.) */
1512 Mark = PRI_QUADSECT;
1513 }
1514 else
1515 SETUSED(theElement,0);
1516 }
1517 #endif
1518 ASSERT(Mark != -1);
1519
1520 /* switch green class to red class? */
1521 if (MARKCLASS(theElement)!=RED_CLASS &&
1522 SWITCHCLASS(CLASS_OF_RULE(MARK2RULEADR(theElement,Mark))))
1523 {
1524 IFDEBUG(gm,1)
1525 UserWriteF(" Switching MARKCLASS=%d for MARK=%d of EID=%d "
1526 "to RED_CLASS\n",
1527 MARKCLASS(theElement),Mark,ID(theElement));
1528 ENDDEBUG
1529 SETMARKCLASS(theElement,RED_CLASS);
1530 }
1531 }
1532
1533 REFINE_ELEMENT_LIST(1,theElement,"");
1534
1535 #ifdef __THREEDIM__
1536 /* choose best tet_red rule according to (*theFullRefRule)() */
1537 if (TAG(theElement)==TETRAHEDRON && MARKCLASS(theElement)==RED_CLASS)
1538 {
1539 #ifndef DUNE_UGGRID_TET_RULESET
1540 if ((Mark==TET_RED || Mark==TET_RED_0_5 ||
1541 Mark==TET_RED_1_3))
1542 #endif
1543 {
1544 PRINTDEBUG(gm,5,("FullRefRule() call with mark=%d\n",Mark))
1545
1546 Mark = (*theFullRefRule)(theElement);
1547 assert( Mark==FULL_REFRULE_0_5 ||
1548 Mark==FULL_REFRULE_1_3 ||
1549 Mark==FULL_REFRULE_2_4);
1550 }
1551 }
1552 #endif
1553
1554 /* get new pattern from mark */
1555 NewPattern = MARK2PAT(theElement,Mark);
1556 IFDEBUG(gm,2)
1557 UserWriteF(" thePattern=%d EdgePattern=%d SidePattern=%d NewPattern=%d Mark=%d\n",
1558 thePattern,theEdgePattern,theSidePattern,NewPattern,Mark);
1559 ENDDEBUG
1560
1561
1562 if (fifoFlag)
1563 {
1564 if (UpdateFIFOLists(theGrid,theElement,thePattern,NewPattern) != GM_OK) return(GM_OK);
1565 }
1566
1567 if (Mark) (*cnt)++;
1568 SETMARK(theElement,Mark);
1569 }
1570
1571 return(GM_OK);
1572 }
1573
1574
1575 #ifdef ModelP
1576
1577 /****************************************************************************/
1578 /*
1579 Gather_AddEdgePattern -
1580
1581 SYNOPSIS:
1582 static int Gather_AddEdgePattern (DDD_OBJ obj, void *data);
1583
1584 PARAMETERS:
1585 . obj
1586 . data
1587
1588 DESCRIPTION:
1589
1590 \return <ul>
1591 int
1592 */
1593 /****************************************************************************/
1594
Gather_AddEdgePattern(DDD::DDDContext &,DDD_OBJ obj,void * data)1595 static int Gather_AddEdgePattern (DDD::DDDContext&, DDD_OBJ obj, void *data)
1596 {
1597 #ifdef __TWODIM__
1598 INT pat;
1599 ELEMENT *theElement = (ELEMENT *)obj;
1600
1601 pat = 0;
1602 GetEdgeInfo(theElement,&pat,ADDPATTERN);
1603
1604 ((INT *)data)[0] = pat;
1605
1606 PRINTDEBUG(gm,4,("Gather_AddEdgePattern(): elem=" EID_FMTX "pat=%08x\n",
1607 EID_PRTX(theElement),pat));
1608 return(GM_OK);
1609 #endif
1610
1611 #ifdef __THREEDIM__
1612 INT addpattern;
1613 EDGE *theEdge = (EDGE *)obj;
1614
1615 addpattern = ADDPATTERN(theEdge);
1616
1617 ((INT *)data)[0] = addpattern;
1618
1619 PRINTDEBUG(gm,4,("Gather_AddEdgePattern(): edge=" ID_FMTX "pat=%08x\n",
1620 ID_PRTX(theEdge),addpattern));
1621 return(GM_OK);
1622 #endif
1623 }
1624
1625
1626 /****************************************************************************/
1627 /*
1628 Scatter_AddEdgePattern -
1629
1630 SYNOPSIS:
1631 static int Scatter_AddEdgePattern (DDD_OBJ obj, void *data);
1632
1633 PARAMETERS:
1634 . obj
1635 . data
1636
1637 DESCRIPTION:
1638
1639 \return <ul>
1640 int
1641 */
1642 /****************************************************************************/
1643
Scatter_AddEdgePattern(DDD::DDDContext &,DDD_OBJ obj,void * data)1644 static int Scatter_AddEdgePattern (DDD::DDDContext&, DDD_OBJ obj, void *data)
1645 {
1646 #ifdef __TWODIM__
1647 INT pat;
1648 ELEMENT *theElement = (ELEMENT *)obj;
1649
1650 /** \todo (HRR 971207): output after SetEdgeInfo (pat not init)? */
1651 PRINTDEBUG(gm,4,("Scatter_AddEdgePattern(): elem=" EID_FMTX "pat=%08x\n",
1652 EID_PRTX(theElement),pat));
1653
1654 pat = ((INT *)data)[0];
1655 SetEdgeInfo(theElement,pat,ADDPATTERN,&);
1656
1657 return(GM_OK);
1658 #endif
1659
1660 #ifdef __THREEDIM__
1661 INT addpattern;
1662 EDGE *theEdge = (EDGE *)obj;
1663
1664 addpattern = MIN(ADDPATTERN(theEdge),((INT *)data)[0]);
1665
1666 PRINTDEBUG(gm,4,("Gather_AddEdgePattern(): edge=" ID_FMTX "pat=%08x\n",
1667 ID_PRTX(theEdge),addpattern));
1668
1669 SETADDPATTERN(theEdge,addpattern);
1670
1671 return(GM_OK);
1672 #endif
1673 }
1674
1675
1676 /****************************************************************************/
1677 /*
1678 ExchangeAddPatterns -
1679
1680 SYNOPSIS:
1681 static INT ExchangeAddPatterns (GRID *theGrid);
1682
1683 PARAMETERS:
1684 . theGrid - pointer to grid structure
1685
1686 DESCRIPTION:
1687
1688 \return <ul>
1689 INT
1690 */
1691 /****************************************************************************/
1692
ExchangeAddPatterns(GRID * theGrid)1693 static INT ExchangeAddPatterns (GRID *theGrid)
1694 {
1695 auto& context = theGrid->dddContext();
1696 const auto& dddctrl = ddd_ctrl(context);
1697
1698 /* exchange addpatterns of edges */
1699 #ifdef __TWODIM__
1700 DDD_IFAOneway(context,
1701 dddctrl.ElementVHIF,GRID_ATTR(theGrid),IF_FORWARD,sizeof(INT),
1702 Gather_AddEdgePattern, Scatter_AddEdgePattern);
1703 #endif
1704 #ifdef __THREEDIM__
1705 DDD_IFAOneway(context,
1706 dddctrl.EdgeVHIF,GRID_ATTR(theGrid),IF_FORWARD,sizeof(INT),
1707 Gather_AddEdgePattern, Scatter_AddEdgePattern);
1708 #endif
1709
1710 return(GM_OK);
1711 }
1712 #endif
1713
1714
1715 /****************************************************************************/
1716 /*
1717 SetAddPatterns -
1718
1719 SYNOPSIS:
1720 static INT SetAddPatterns (GRID *theGrid);
1721
1722 PARAMETERS:
1723 . theGrid - pointer to grid structure
1724
1725 DESCRIPTION:
1726
1727 \return <ul>
1728 INT
1729 */
1730 /****************************************************************************/
1731
SetAddPatterns(GRID * theGrid)1732 static INT SetAddPatterns (GRID *theGrid)
1733 {
1734 INT j;
1735 ELEMENT *theElement;
1736 EDGE *theEdge;
1737
1738 /* set additional pattern on the edges */
1739 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL;
1740 theElement=SUCCE(theElement))
1741 {
1742 if (MARKCLASS(theElement)!=RED_CLASS) continue;
1743
1744 REFINE_ELEMENT_LIST(1,theElement,"SetAddPatterns(): addpattern=0");
1745
1746 for (j=0; j<EDGES_OF_ELEM(theElement); j++)
1747 {
1748 /* no green elements for this edge if there is no edge node */
1749 if (!NODE_OF_RULE(theElement,MARK(theElement),j))
1750 continue;
1751
1752 theEdge=GetEdge(CORNER_OF_EDGE_PTR(theElement,j,0),
1753 CORNER_OF_EDGE_PTR(theElement,j,1));
1754 ASSERT(theEdge != NULL);
1755
1756 /* ADDPATTERN is now set to 0 for all edges of red elements */
1757 SETADDPATTERN(theEdge,0);
1758 }
1759 }
1760
1761 #ifdef ModelP
1762 if (ExchangeAddPatterns(theGrid)) RETURN(GM_FATAL);
1763 #endif
1764
1765 return(GM_OK);
1766 }
1767
1768
1769 /****************************************************************************/
1770 /*
1771 BuildGreenClosure -
1772
1773 SYNOPSIS:
1774 static INT BuildGreenClosure (GRID *theGrid);
1775
1776 PARAMETERS:
1777 . theGrid - pointer to grid structure
1778
1779 DESCRIPTION:
1780
1781 \return <ul>
1782 INT
1783 */
1784 /****************************************************************************/
1785
BuildGreenClosure(GRID * theGrid)1786 static INT BuildGreenClosure (GRID *theGrid)
1787 {
1788 INT i;
1789 ELEMENT *theElement;
1790 EDGE *theEdge;
1791 #ifdef __THREEDIM__
1792 INT j;
1793 #endif
1794
1795 /* build a green covering around the red elements */
1796 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL;
1797 theElement=SUCCE(theElement))
1798 {
1799 #ifdef __ANISOTROPIC__
1800 if (MARKCLASS(theElement)==RED_CLASS &&
1801 !(TAG(theElement)==PRISM && MARK(theElement)==PRI_QUADSECT)) continue;
1802
1803 ASSERT(MARKCLASS(theElement)!=RED_CLASS ||
1804 (MARKCLASS(theElement)==RED_CLASS && TAG(theElement)==PRISM
1805 && MARK(theElement)==PRI_QUADSECT));
1806 #else
1807 if (MARKCLASS(theElement)==RED_CLASS) continue;
1808 #endif
1809
1810 SETUPDATE_GREEN(theElement,0);
1811
1812 /* if edge node exists element needs to be green */
1813 for (i=0; i<EDGES_OF_ELEM(theElement); i++)
1814 {
1815 theEdge=GetEdge(CORNER_OF_EDGE_PTR(theElement,i,0),
1816 CORNER_OF_EDGE_PTR(theElement,i,1));
1817 ASSERT(theEdge != NULL);
1818
1819 /* if edge is refined this will be a green element */
1820 if (ADDPATTERN(theEdge) == 0)
1821 {
1822 /* for pyramids, prisms and hexhedra Patterns2Rules returns 0 */
1823 /* for non red elements, because there is no complete rule set */
1824 /* switch to mark COPY, because COPY rule refines no edges */
1825 #ifdef DUNE_UGGRID_TET_RULESET
1826 if (DIM==3 && TAG(theElement)!=TETRAHEDRON)
1827 #else
1828 if (DIM==3)
1829 #endif
1830 {
1831 /* set to no-empty rule, e.g. COPY rule */
1832 #ifdef __ANISOTROPIC__
1833 if (MARKCLASS(theElement) != RED_CLASS)
1834 #endif
1835 SETMARK(theElement,COPY);
1836
1837 /* no existing edge node renew green refinement */
1838 if (MIDNODE(theEdge)==NULL)
1839 {
1840 SETUPDATE_GREEN(theElement,1);
1841 }
1842 }
1843 /* tetrahedra in 3D and 2D elements have a complete rule set */
1844 else if (MARK(theElement) == NO_REFINEMENT)
1845 {
1846 IFDEBUG(gm,2)
1847 UserWriteF(" ERROR: green tetrahedron with no rule! "
1848 "EID=%d TAG=%d "
1849 "REFINECLASS=%d REFINE=%d MARKCLASS=%d MARK=%d\n",
1850 ID(theElement),TAG(theElement),REFINECLASS(theElement),
1851 REFINE(theElement),MARKCLASS(theElement),
1852 MARK(theElement));
1853 ENDDEBUG
1854 }
1855
1856 #ifdef __ANISOTROPIC__
1857 if (MARKCLASS(theElement) != RED_CLASS)
1858 #endif
1859 SETMARKCLASS(theElement,GREEN_CLASS);
1860 }
1861 else
1862 {
1863 /* existing edge node is deleted */
1864 /* renew green refinement if element will be a green one */
1865 if (MIDNODE(theEdge)!=NULL)
1866 SETUPDATE_GREEN(theElement,1);
1867 }
1868 }
1869
1870 #ifdef __THREEDIM__
1871 /* if side node exists element needs to be green */
1872 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
1873 {
1874 ELEMENT *theNeighbor;
1875
1876 theNeighbor = NBELEM(theElement,i);
1877
1878 if (theNeighbor==NULL) continue;
1879
1880 for (j=0; j<SIDES_OF_ELEM(theNeighbor); j++)
1881 if (NBELEM(theNeighbor,j) == theElement)
1882 break;
1883
1884 #ifdef ModelP
1885 if (j >= SIDES_OF_ELEM(theNeighbor))
1886 {
1887 ASSERT(EGHOST(theElement) && EGHOST(theNeighbor));
1888 continue;
1889 }
1890 #else
1891 ASSERT(j<SIDES_OF_ELEM(theNeighbor));
1892 #endif
1893
1894 if (NODE_OF_RULE(theNeighbor,MARK(theNeighbor),
1895 EDGES_OF_ELEM(theNeighbor)+j))
1896 {
1897 #ifdef DUNE_UGGRID_TET_RULESET
1898 if (TAG(theNeighbor)==TETRAHEDRON)
1899 printf("ERROR: no side nodes for tetrahedra! side=%d\n",j);
1900 #endif
1901 #ifdef __ANISOTROPIC__
1902 if (MARKCLASS(theElement) != RED_CLASS)
1903 #endif
1904 SETMARKCLASS(theElement,GREEN_CLASS);
1905 }
1906
1907
1908 /* side node change? */
1909 if ((!NODE_OF_RULE(theNeighbor,REFINE(theNeighbor),
1910 EDGES_OF_ELEM(theNeighbor)+j) &&
1911 NODE_OF_RULE(theNeighbor,MARK(theNeighbor),
1912 EDGES_OF_ELEM(theNeighbor)+j)) ||
1913 (NODE_OF_RULE(theNeighbor,REFINE(theNeighbor),
1914 EDGES_OF_ELEM(theNeighbor)+j) &&
1915 !NODE_OF_RULE(theNeighbor,MARK(theNeighbor),
1916 EDGES_OF_ELEM(theNeighbor)+j)))
1917 {
1918 SETUPDATE_GREEN(theElement,1);
1919 }
1920 }
1921 #endif
1922
1923
1924 #ifndef ModelP
1925 /* if element is green before refinement and will be green after */
1926 /* refinement and nothing changes -> reset USED flag */
1927 /* in parallel case: one communication to determine the minimum */
1928 /* over all copies of an green element would be needed */
1929 if (REFINECLASS(theElement)==GREEN_CLASS &&
1930 MARKCLASS(theElement)==GREEN_CLASS && UPDATE_GREEN(theElement)==0)
1931 {
1932 /* do not renew green refinement */
1933 SETUSED(theElement,0);
1934 }
1935 #ifdef __ANISOTROPIC__
1936 if (MARKCLASS(theElement)==RED_CLASS && UPDATE_GREEN(theElement)==0)
1937 {
1938 ASSERT(TAG(theElement)==PRISM && MARK(theElement)==PRI_QUADSECT);
1939 SETUSED(theElement,0);
1940 }
1941 #endif
1942 #endif
1943 }
1944
1945 return(GM_OK);
1946 }
1947
1948 #if defined(ModelP) && defined(Debug)
1949
1950 #define CEIL(n) ((n)+((ALIGNMENT-((n)&(ALIGNMENT-1)))&(ALIGNMENT-1)))
1951
1952 /****************************************************************************/
1953 /*
1954 Gather_ElementInfo -
1955
1956 SYNOPSIS:
1957 static int Gather_ElementInfo (DDD_OBJ obj, void *Data);
1958
1959 PARAMETERS:
1960 . obj
1961 . Data
1962
1963 DESCRIPTION:
1964
1965 \return <ul>
1966 int
1967 */
1968 /****************************************************************************/
1969
Gather_ElementInfo(DDD::DDDContext &,DDD_OBJ obj,void * Data)1970 static int Gather_ElementInfo (DDD::DDDContext&, DDD_OBJ obj, void *Data)
1971 {
1972 INT epat,eaddpat;
1973 ELEMENT *theElement = (ELEMENT *)obj;
1974 char *data = (char *)Data;
1975
1976 PRINTDEBUG(gm,4,("Gather_ElementInfo(): elem=" EID_FMTX "\n",
1977 EID_PRTX(theElement)));
1978
1979 memcpy(data,theElement,sizeof(struct generic_element));
1980 data += CEIL(sizeof(struct generic_element));
1981
1982 epat = 0;
1983 GetEdgeInfo(theElement,&epat,PATTERN);
1984 ((int *)data)[0] = epat;
1985 data += sizeof(INT);
1986
1987 eaddpat = 0;
1988 GetEdgeInfo(theElement,&eaddpat,ADDPATTERN);
1989 ((int *)data)[0] = eaddpat;
1990
1991 return(GM_OK);
1992 }
1993
1994 /* this macro compares control word macros of two elements */
1995 #define COMPARE_MACRO(elem0,elem1,macro,print) \
1996 if (macro(elem0) != macro(elem1)) \
1997 { \
1998 print("e=" EID_FMTX " macro=%s differs value0=%d value1=%d \n", \
1999 EID_PRTX(elem0),# macro,macro(elem0),macro(elem1)); \
2000 assert(0); \
2001 }
2002
2003 #ifdef DUNE_UGGRID_TET_RULESET
2004 #define COMPARE_MACROX(elem0,elem1,macro,print) \
2005 { \
2006 INT _mark0,_mark1,_pat0,_pat1; \
2007 \
2008 _mark0 = macro(elem0); _mark1 = macro(elem1); \
2009 _pat0 = MARK2PAT((elem0),_mark0); _pat1 = MARK2PAT((elem1),_mark1); \
2010 if ((_pat0 & ((1<<10)-1)) != (_pat1 & ((1<<10)-1))) \
2011 COMPARE_MACRO(elem0,elem1,macro,print) \
2012 }
2013 #else
2014 #define COMPARE_MACROX(elem0,elem1,macro,print) \
2015 COMPARE_MACRO(elem0,elem1,macro,print)
2016 #endif
2017
2018 /* this macro compares two values */
2019 #define COMPARE_VALUE(elem0,val0,val1,string,print) \
2020 if ((val0) != (val1)) \
2021 { \
2022 (print)("e=" EID_FMTX " %s differs value0=%d value1=%d \n", \
2023 EID_PRTX(elem0),(string),(val0),(val1)); \
2024 assert(0); \
2025 }
2026
2027
2028 /****************************************************************************/
2029 /*
2030 Scatter_ElementInfo -
2031
2032 SYNOPSIS:
2033 static int Scatter_ElementInfo (DDD_OBJ obj, void *Data);
2034
2035 PARAMETERS:
2036 . obj
2037 . Data
2038
2039 DESCRIPTION:
2040
2041 \return <ul>
2042 int
2043 */
2044 /****************************************************************************/
2045
Scatter_ElementInfo(DDD::DDDContext &,DDD_OBJ obj,void * Data)2046 static int Scatter_ElementInfo (DDD::DDDContext&, DDD_OBJ obj, void *Data)
2047 {
2048 INT epat,mpat,eaddpat,maddpat;
2049 ELEMENT *theElement = (ELEMENT *)obj;
2050 struct generic_element ge;
2051 ELEMENT *theMaster = (ELEMENT *)≥
2052 char *data = (char *)Data;
2053
2054 memcpy(theMaster,data,sizeof(struct generic_element));
2055 data += CEIL(sizeof(struct generic_element));
2056
2057 PRINTDEBUG(gm,4,("Scatter_ElementInfo(): Comparing elem="
2058 EID_FMTX " master=" EID_FMTX "\n",EID_PRTX(theElement),
2059 EID_PRTX(theMaster)));
2060
2061 /* now compare the control entries of master with it local copy */
2062 COMPARE_MACRO(theElement,theMaster,REFINECLASS,PrintDebug)
2063 COMPARE_MACRO(theElement,theMaster,MARKCLASS,PrintDebug)
2064 COMPARE_MACROX(theElement,theMaster,REFINE,PrintDebug)
2065 COMPARE_MACROX(theElement,theMaster,MARK,PrintDebug)
2066 COMPARE_MACRO(theElement,theMaster,COARSEN,PrintDebug)
2067 COMPARE_MACRO(theElement,theMaster,USED,PrintDebug)
2068 /* #ifndef DUNE_UGGRID_TET_RULESET */
2069 #ifdef __THREEDIM__
2070 COMPARE_MACRO(theElement,theMaster,SIDEPATTERN,PrintDebug)
2071 #endif
2072 /* #endif */
2073
2074 epat = 0;
2075 GetEdgeInfo(theElement,&epat,PATTERN);
2076 mpat = ((INT *)data)[0];
2077 data += sizeof(INT);
2078 COMPARE_VALUE(theElement,epat,mpat,"EdgePattern",PrintDebug)
2079
2080 eaddpat = 0;
2081 GetEdgeInfo(theElement,&eaddpat,ADDPATTERN);
2082 maddpat = ((INT *)data)[0];
2083 COMPARE_VALUE(theElement,eaddpat,maddpat,"EdgeAddPattern",PrintDebug)
2084
2085 return(GM_OK);
2086 }
2087
CheckElementInfo(GRID * theGrid)2088 static INT CheckElementInfo (GRID *theGrid)
2089 {
2090 auto& context = theGrid->dddContext();
2091 const auto& dddctrl = ddd_ctrl(context);
2092
2093 /* exchange element info */
2094 DDD_IFAOneway(context,
2095 dddctrl.ElementVHIF,GRID_ATTR(theGrid),IF_FORWARD,
2096 CEIL(sizeof(struct generic_element))+2*sizeof(INT),
2097 Gather_ElementInfo, Scatter_ElementInfo);
2098
2099 return(GM_OK);
2100 }
2101 #endif
2102
2103
2104 /****************************************************************************/
2105 /* */
2106 /* Function: GridClosure */
2107 /* */
2108 /* Purpose: compute closure for next level. A closure can only be */
2109 /* determined if the rule set for the used elements is complete. */
2110 /* This means that for all side and edge patterns possible for */
2111 /* an element type exists a rule which closes the element. */
2112 /* In this case a FIFO for computing the closure is not needed */
2113 /* any more and the closure can be computed in one step. */
2114 /* */
2115 /* Param: GRID *theGrid: pointer to grid structure */
2116 /* */
2117 /* return: INT >0: elements will be refined */
2118 /* INT 0: no elements will be refined */
2119 /* INT -1: an error occured */
2120 /* */
2121 /****************************************************************************/
2122 /****************************************************************************/
2123 /** \brief Compute closure for next level
2124
2125 \param theGrid - pointer to grid structure
2126
2127 This function computes closure for next level. A closure can only be
2128 determined if the rule set for the used elements is complete. This means
2129 that for all side and edge patterns possible for an element type exists
2130 a rule which closes the element. In this case a FIFO for computing the
2131 closure is not needed any more and the closure can be computed in one step.
2132
2133 \return <ul>
2134 .n >0 elements will be refined
2135 .n =0 no elements will be refined
2136 .n =-1 an error accured
2137 */
2138 /****************************************************************************/
2139
GridClosure(GRID * theGrid)2140 static int GridClosure (GRID *theGrid)
2141 {
2142 INT cnt;
2143
2144 /* initialize used control word entries */
2145 if (PrepareGridClosure(theGrid) != GM_OK) RETURN(GM_ERROR);
2146
2147 /* compute pattern on edges and elements */
2148 if (ComputePatterns(theGrid) != GM_OK) RETURN(GM_ERROR);
2149
2150 firstElement = PFIRSTELEMENT(theGrid);
2151
2152 if (fifoFlag)
2153 if (InitClosureFIFO() != GM_OK) return(GM_OK);
2154
2155 /* fifo loop */
2156 do
2157 {
2158 #ifdef __THREEDIM__
2159 #if defined(ModelP) && defined(DUNE_UGGRID_TET_RULESET)
2160 /* edge pattern is needed consistently in CorrectTetrahedronSidePattern() */
2161 if (!refine_seq)
2162 {
2163 if (ExchangeEdgeClosureInfo(theGrid) != GM_OK) return(GM_ERROR);
2164 }
2165 #endif
2166
2167 /* set side patterns on the elements */
2168 if (SetElementSidePatterns(theGrid,firstElement) != GM_OK) RETURN(GM_ERROR);
2169 #endif
2170
2171 #ifdef ModelP
2172 if (ExchangeClosureInfo(theGrid) != GM_OK) RETURN(GM_ERROR);
2173 #endif
2174
2175 /* set rules on the elements */
2176 if (SetElementRules(theGrid,firstElement,&cnt) != GM_OK) RETURN(GM_ERROR);
2177
2178 }
2179 /* exit only if fifo not active or fifo queue */
2180 /* empty or all processor have finished closure */
2181 while (fifoFlag && UpdateClosureFIFO(theGrid) &&
2182 ManageParallelFIFO(theGrid->ppifContext(), firstElement));
2183
2184 /* set patterns on all edges of red elements */
2185 if (SetAddPatterns(theGrid) != GM_OK) RETURN(GM_ERROR);
2186
2187 /* build the closure around the red elements */
2188 if (BuildGreenClosure(theGrid) != GM_OK) RETURN(GM_ERROR);
2189
2190 #if defined(Debug) && defined(ModelP)
2191 if (CheckElementInfo(theGrid)) RETURN(GM_ERROR);
2192 #endif
2193
2194 return(cnt);
2195 }
2196
2197
2198
2199 /****************************************************************************/
2200 /*
2201 GetNeighborSons - fill SonList for theElement
2202
2203 SYNOPSIS:
2204 static INT GetNeighborSons (ELEMENT *theElement, ELEMENT *theSon,
2205 ELEMENT *SonList[MAX_SONS], int count, int nsons);
2206
2207 PARAMETERS:
2208 . theElement - father element
2209 . theSon - currently visited son
2210 . SonList[MAX_SONS] - the list of sons
2211 . count - son count
2212 . sons - number of sons
2213
2214 DESCRIPTION:
2215 This function fills SonList for theElement with a breadth first search
2216
2217 \return <ul>
2218 INT
2219 .n new son count
2220 */
2221 /****************************************************************************/
2222
GetNeighborSons(ELEMENT * theElement,ELEMENT * theSon,ELEMENT * SonList[MAX_SONS],int count,int nsons)2223 static INT GetNeighborSons (ELEMENT *theElement, ELEMENT *theSon,
2224 ELEMENT *SonList[MAX_SONS], int count, int nsons)
2225 {
2226 ELEMENT *NbElement;
2227 int i,j, startson, stopson;
2228
2229 startson = count;
2230
2231 for (i=0; i<SIDES_OF_ELEM(theSon); i++)
2232 {
2233 NbElement = NBELEM(theSon,i);
2234 if (NbElement == NULL) continue;
2235 if (EFATHER(NbElement) == theElement)
2236 {
2237 /* is NbElement already in list */
2238 for (j=0; j<count; j++)
2239 if (SonList[j] == NbElement)
2240 break;
2241 if (j==count && count<nsons)
2242 SonList[count++] = NbElement;
2243 }
2244 }
2245 if (count == nsons) return(count);
2246
2247 stopson = count;
2248 for (i=startson; i<stopson; i++)
2249 {
2250 if (count<nsons) count = GetNeighborSons(theElement,SonList[i],
2251 SonList,count,nsons);
2252 else return(count);
2253 }
2254 return(count);
2255 }
2256
2257
2258 #ifdef ModelP
2259
2260
2261 /****************************************************************************/
2262 /*
2263 GetSons - fill SonList for theElement
2264
2265 SYNOPSIS:
2266 INT GetAllSons (ELEMENT *theElement, ELEMENT *SonList[MAX_SONS]);
2267
2268 PARAMETERS:
2269 . theElement
2270 . SonList[MAX_SONS]
2271
2272 DESCRIPTION:
2273 This function fills SonList for theElement
2274
2275 \return <ul>
2276 INT
2277 .n 0 if ok
2278 .n 1 if an error occurs
2279 */
2280 /****************************************************************************/
2281
GetAllSons(const ELEMENT * theElement,ELEMENT * SonList[MAX_SONS])2282 INT NS_DIM_PREFIX GetAllSons (const ELEMENT *theElement, ELEMENT *SonList[MAX_SONS])
2283 {
2284 ELEMENT *son;
2285 int SonID,i;
2286
2287 ASSERT(theElement != NULL);
2288
2289 for (SonID=0; SonID<MAX_SONS; SonID++)
2290 SonList[SonID] = NULL;
2291
2292 if (NSONS(theElement) == 0) return(GM_OK);
2293
2294 SonID = 0;
2295
2296 for (i=0; i<2; i++)
2297 {
2298 if (i == 0)
2299 son = SON(theElement,PRIO2INDEX(PrioMaster));
2300 else
2301 son = SON(theElement,PRIO2INDEX(PrioHGhost));
2302
2303 if (son == NULL)
2304 continue;
2305 else
2306 SonList[SonID++] = son;
2307
2308 while (SUCCE(son) != NULL)
2309 {
2310 if (EFATHER(SUCCE(son)) == theElement
2311 && PRIO2INDEX(EPRIO(son))==PRIO2INDEX(EPRIO(SUCCE(son)))
2312 )
2313 {
2314 SonList[SonID++] = SUCCE(son);
2315 son = SUCCE(son);
2316 ASSERT(SonID <= MAX_SONS);
2317 }
2318 else
2319 break;
2320 }
2321 }
2322
2323 return(GM_OK);
2324 }
2325 #endif
2326
2327
2328 /****************************************************************************/
2329 /*
2330 GetSons -
2331
2332 SYNOPSIS:
2333 INT GetSons (ELEMENT *theElement, ELEMENT *SonList[MAX_SONS]);
2334
2335 PARAMETERS:
2336 . theElement
2337 . SonList[MAX_SONS]
2338
2339 DESCRIPTION:
2340
2341 \return <ul>
2342 INT
2343 */
2344 /****************************************************************************/
2345
GetSons(const ELEMENT * theElement,ELEMENT * SonList[MAX_SONS])2346 INT NS_DIM_PREFIX GetSons (const ELEMENT *theElement, ELEMENT *SonList[MAX_SONS])
2347 {
2348 int SonID;
2349 ELEMENT *son;
2350
2351 if (theElement==NULL) RETURN(GM_ERROR);
2352
2353 for (SonID=0; SonID<MAX_SONS; SonID++)
2354 SonList[SonID] = NULL;
2355
2356 if (NSONS(theElement) == 0) return(GM_OK);
2357
2358 SonID = 0;
2359 SonList[SonID++] = son = SON(theElement,PRIO2INDEX(PrioMaster));
2360
2361 if (son == NULL) return(GM_OK);
2362
2363 while (SUCCE(son) != NULL)
2364 {
2365 if (EFATHER(SUCCE(son)) == theElement
2366 #ifdef ModelP
2367 && PRIO2INDEX(EPRIO(son))==PRIO2INDEX(EPRIO(SUCCE(son)))
2368 #endif
2369 )
2370 {
2371 SonList[SonID++] = SUCCE(son);
2372 son = SUCCE(son);
2373 ASSERT(SonID <= MAX_SONS);
2374 }
2375 else
2376 break;
2377
2378 }
2379
2380 return(GM_OK);
2381 }
2382
2383 /****************************************************************************/
2384 /* */
2385 /* Function: RestrictElementMark */
2386 /* */
2387 /* Purpose: restrict refinement marks of an element whose sons are */
2388 /* further marked for refinement */
2389 /* */
2390 /* Param: ELEMENT *theElement: pointer to the element */
2391 /* */
2392 /* return: INT: =0 ok */
2393 /* >0 error */
2394 /* */
2395 /****************************************************************************/
2396 /****************************************************************************/
2397 /*
2398 RestrictElementMark - restrict refinement marks
2399
2400 SYNOPSIS:
2401 static INT RestrictElementMark(ELEMENT *theElement);
2402
2403 PARAMETERS:
2404 . theElement - pointer to the element
2405
2406 DESCRIPTION:
2407 This function restricts refinement marks of an element whose sons are
2408 further marked for refinement
2409
2410 \return <ul>
2411 INT
2412 .n =0 - ok
2413 .n >0 - error
2414 */
2415 /****************************************************************************/
2416
RestrictElementMark(ELEMENT * theElement)2417 static INT RestrictElementMark(ELEMENT *theElement)
2418 {
2419 #ifdef __THREEDIM__
2420 #ifdef DUNE_UGGRID_TET_RULESET
2421 EDGE *theEdge;
2422 int j,Rule,Pattern;
2423 #endif
2424 #endif
2425
2426 if (MARKCLASS(theElement)==RED_CLASS)
2427 {
2428 /** \todo this mark is from DropMarks()! */
2429 /* theElement is marked from outside */
2430 /** \todo edit this for new element type or */
2431 /* for different restrictions */
2432 switch (TAG(theElement))
2433 {
2434 #ifdef __TWODIM__
2435 case TRIANGLE :
2436 SETMARK(theElement,T_RED);
2437 break;
2438 case QUADRILATERAL :
2439 SETMARK(theElement,Q_RED);
2440 break;
2441 #endif
2442 #ifdef __THREEDIM__
2443 case TETRAHEDRON :
2444 #ifdef DUNE_UGGRID_TET_RULESET
2445 if (MARK(theElement)!=RED)
2446 /** \todo Is REFINE always as red rule available? */
2447 SETMARK(theElement,REFINE(theElement));
2448 #else
2449 SETMARK(theElement,TET_RED);
2450 #endif
2451 break;
2452 case PYRAMID :
2453 SETMARK(theElement,PYR_RED);
2454 break;
2455 case PRISM :
2456 SETMARK(theElement,PRI_RED);
2457 break;
2458 case HEXAHEDRON :
2459 SETMARK(theElement,HEXA_RED);
2460 break;
2461 #endif
2462
2463 default :
2464 ASSERT(0);
2465 }
2466 }
2467 else
2468 {
2469 /** \todo edit this for new element type or for different restrictions */
2470 switch (TAG(theElement))
2471 {
2472 #ifdef __TWODIM__
2473 case TRIANGLE :
2474 SETMARK(theElement,T_RED);
2475 break;
2476
2477 case QUADRILATERAL :
2478 SETMARK(theElement,Q_RED);
2479 break;
2480 #endif
2481 #ifdef __THREEDIM__
2482 case TETRAHEDRON :
2483 #ifdef DUNE_UGGRID_TET_RULESET
2484 /* theElement is not marked from outside, */
2485 /* so find a reg. rule being consistent */
2486 /* with those neighbors of all sons of */
2487 /* theElement which are marked for refine.*/
2488 /* this choice will make sure these marks */
2489 /* will not be distroyed. */
2490 Pattern = RULE2PAT(theElement,
2491 REFINE(theElement));
2492 for (j=0; j<EDGES_OF_ELEM(theElement); j++)
2493 {
2494 theEdge=GetEdge(CORNER_OF_EDGE_PTR(theElement,j,0),
2495 CORNER_OF_EDGE_PTR(theElement,j,1));
2496 ASSERT(theEdge != NULL);
2497
2498 /** \todo What's on when MIDNODE exists?? */
2499 if (MIDNODE(theEdge)==NULL)
2500 {
2501 theEdge=GetEdge(
2502 SONNODE(CORNER_OF_EDGE_PTR(theElement,j,0)),
2503 SONNODE(CORNER_OF_EDGE_PTR(theElement,j,1)));
2504 ASSERT(theEdge != NULL);
2505
2506 /** \todo Is ADDPATTERN needed for fitting with other green elements?? */
2507 if (ADDPATTERN(theEdge))
2508 Pattern |= (1<<j);
2509 PRINTDEBUG(gm,4,("RestrictElementMark(): modified Pattern=%d bisects now edge=%d too\n",Pattern,j))
2510 }
2511 }
2512 Rule = PATTERN2RULE(theElement,Pattern);
2513 SETMARK(theElement,RULE2MARK(theElement,Rule));
2514 /** \todo delete this old code
2515 SETMARK(theElement,PATTERN2MARK(theElement,Pattern)); */
2516 /** \todo this would be the quick fix
2517 SETMARK(theElement,FULL_REFRULE); */
2518 #else
2519 SETMARK(theElement,TET_RED);
2520 #endif
2521 break;
2522 case PYRAMID :
2523 SETMARK(theElement,PYR_RED);
2524 break;
2525 case PRISM :
2526 #ifdef __ANISOTROPIC__
2527 SETMARK(theElement,PRI_QUADSECT);
2528 #else
2529 SETMARK(theElement,PRI_RED);
2530 #endif
2531 break;
2532 case HEXAHEDRON :
2533 SETMARK(theElement,HEXA_RED);
2534 break;
2535 #endif
2536 default :
2537 ASSERT(0);
2538 }
2539 SETMARKCLASS(theElement,RED_CLASS);
2540 }
2541
2542 return(GM_OK);
2543 }
2544
2545 /****************************************************************************/
2546 /*
2547 RestrictMarks - restrict refinement marks when going down
2548
2549 SYNOPSIS:
2550 static INT RestrictMarks (GRID *theGrid);
2551
2552 PARAMETERS:
2553 . theGrid - pointer to grid structure
2554
2555 DESCRIPTION:
2556 This function restricts refinement marks when going down
2557
2558 \return <ul>
2559 INT
2560 .n 0 if ok
2561 .n >0 if an error occurs
2562 */
2563 /****************************************************************************/
2564
RestrictMarks(GRID * theGrid)2565 static INT RestrictMarks (GRID *theGrid)
2566 {
2567 ELEMENT *theElement,*SonList[MAX_SONS];
2568 int i,flag;
2569
2570 for (theElement=FIRSTELEMENT(theGrid); theElement!=NULL;
2571 theElement=SUCCE(theElement))
2572 {
2573 if (GetSons(theElement,SonList)!=GM_OK) RETURN(GM_ERROR);
2574
2575 if (hFlag)
2576 {
2577 if (
2578 /* if element is not refined anyway, */
2579 /* then there are no restrictions to apply */
2580 REFINE(theElement) == NO_REFINEMENT ||
2581
2582 /* irregular elements are marked by estimator, */
2583 /* because they are leaf elements */
2584 ECLASS(theElement) == YELLOW_CLASS ||
2585 ECLASS(theElement) == GREEN_CLASS ||
2586
2587 /* regular elements with YELLOW_CLASS copies are */
2588 /* marked by estimator, because the marks are dropped */
2589 REFINECLASS(theElement) == YELLOW_CLASS
2590 )
2591 {
2592 continue;
2593 }
2594
2595 /* regular elements with GREEN_CLASS refinement */
2596 /* go to no refinement or red refinement */
2597 if (REFINECLASS(theElement)==GREEN_CLASS)
2598 {
2599 for (i=0; i<NSONS(theElement); i++)
2600 {
2601 #ifdef ModelP
2602 if (SonList[i] == NULL) break;
2603 #endif
2604
2605 /* Is the son marked for further refinement */
2606 if (MARK(SonList[i])>NO_REFINEMENT)
2607 {
2608 if (RestrictElementMark(theElement)) RETURN(GM_ERROR);
2609
2610 /* this must be done only once for each element */
2611 break;
2612 }
2613 }
2614 continue;
2615 }
2616
2617 /* regular elements with regular refinement are */
2618 /* the only ones to coarsen */
2619 if (REFINECLASS(theElement) == RED_CLASS)
2620 {
2621 #ifndef __ANISOTROPIC__
2622 SETMARK(theElement,REFINE(theElement));
2623 SETMARKCLASS(theElement,REFINECLASS(theElement));
2624 #else
2625 ASSERT(MARK(theElement)>=1);
2626 #endif
2627 }
2628 }
2629
2630 #ifdef ModelP
2631 /* if no (or not all) sons are found by GetSons() on */
2632 /* this proc then coarsening is not allowed */
2633 if (REFINECLASS(theElement)==RED_CLASS &&
2634 SonList[0]==NULL) continue;
2635 #endif
2636
2637 flag = 0;
2638 for (i=0; SonList[i]!=NULL; i++)
2639 {
2640 /* if not all sons are marked no unrefinement is possible */
2641 if (!COARSEN(SonList[i]) || REFINECLASS(SonList[i])==RED_CLASS)
2642 {
2643 flag = 1;
2644 break;
2645 }
2646 }
2647
2648 if (flag) continue;
2649
2650 /* preserve regular refinement marks */
2651 if (hFlag==0 && SonList[0]==NULL) continue;
2652
2653 /* remove refinement */
2654 SETMARK(theElement,NO_REFINEMENT);
2655 SETMARKCLASS(theElement,NO_CLASS);
2656 SETCOARSEN(theElement,1);
2657 }
2658
2659 return(GM_OK);
2660 }
2661
2662 /****************************************************************************/
2663 /*
2664 ComputeCopies - determine copy elements from node classes
2665
2666 SYNOPSIS:
2667 static int ComputeCopies (GRID *theGrid);
2668
2669 PARAMETERS:
2670 . theGrid - pointer to grid structure
2671
2672 DESCRIPTION:
2673 This function determines copy elements from node classes
2674
2675 \return <ul>
2676 int
2677 .n GM_OK if ok
2678 */
2679 /****************************************************************************/
2680
ComputeCopies(GRID * theGrid)2681 static int ComputeCopies (GRID *theGrid)
2682 {
2683 ELEMENT *theElement;
2684 int flag;
2685 int cnt = 0;
2686 PRINTDEBUG(gm,1,("ComputeCopies on level %d\n",GLEVEL(theGrid)));
2687
2688 [[maybe_unused]] const int me = theGrid->ppifContext().me();
2689
2690 /* set class of all dofs on next level to 0 */
2691 ClearNextNodeClasses(theGrid);
2692
2693 /* seed dofs of regularly and irregularly refined elements to 3 */
2694 flag = 0;
2695 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL;
2696 theElement=SUCCE(theElement))
2697 {
2698 if (MARK(theElement)!=NO_REFINEMENT &&
2699 (MARKCLASS(theElement)==RED_CLASS ||
2700 MARKCLASS(theElement)==GREEN_CLASS))
2701 {
2702 SeedNextNodeClasses(theElement);
2703 flag=1; /* there is at least one element to be refined */
2704 }
2705 }
2706
2707 /* copy all option or neighborhood */
2708 if (rFlag==GM_COPY_ALL)
2709 {
2710 #ifdef ModelP
2711 flag = UG_GlobalMaxINT(theGrid->ppifContext(), flag);
2712 #endif
2713 if (flag)
2714 for (theElement=FIRSTELEMENT(theGrid); theElement!=NULL;
2715 theElement=SUCCE(theElement))
2716 {
2717 SeedNextNodeClasses(theElement);
2718 }
2719 }
2720 else
2721 {
2722 PropagateNextNodeClasses(theGrid);
2723 }
2724
2725 /* an element is copied if it has a dof of class 2 and higher */
2726 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL;
2727 theElement=SUCCE(theElement))
2728 {
2729 INT maxclass = 0;
2730
2731 if ((MARK(theElement)==NO_REFINEMENT)&&
2732 ((maxclass=MaxNextNodeClass(theElement))>=MINVNCLASS))
2733 {
2734 PRINTDEBUG(gm,1,(PFMT "ComputeCopies(): level=%d e=" EID_FMTX "yellow marked\n",
2735 me,LEVEL(theElement),EID_PRTX(theElement)));
2736 SETMARK(theElement,COPY);
2737 SETMARKCLASS(theElement,YELLOW_CLASS);
2738 cnt++;
2739 }
2740 else
2741 {
2742 PRINTDEBUG(gm,1,(PFMT "ComputeCopies(): level=%d e=" EID_FMTX "not yellow marked"
2743 " mark=%d maxclass=%d\n",
2744 me,LEVEL(theElement),EID_PRTX(theElement),MARK(theElement),maxclass));
2745 }
2746 }
2747
2748 return(cnt);
2749 }
2750
2751 /****************************************************************************/
2752 /*
2753 CheckElementContextConsistency - check NTYPE flags of nodes
2754
2755 SYNOPSIS:
2756 static void CheckElementContextConsistency(ELEMENT *theElement);
2757
2758 PARAMETERS:
2759 . theElement - element to check
2760
2761 DESCRIPTION:
2762 This function checks NTYPE flags of nodes in elementcontextt with the sons
2763
2764 \return <ul>
2765 void
2766 */
2767 /****************************************************************************/
2768 #ifdef Debug
CheckElementContextConsistency(ELEMENT * theElement,ELEMENTCONTEXT theElementContext)2769 static void CheckElementContextConsistency(ELEMENT *theElement,
2770 ELEMENTCONTEXT theElementContext)
2771 {
2772 int i;
2773 [[maybe_unused]] int errorflag = 0;
2774 int errortype[MAX_CORNERS_OF_ELEM+MAX_NEW_CORNERS_DIM];
2775 int correcttype[MAX_CORNERS_OF_ELEM+MAX_NEW_CORNERS_DIM];
2776
2777 for (i=0; i<MAX_CORNERS_OF_ELEM+MAX_NEW_CORNERS_DIM; i++)
2778 errortype[i] = correcttype[i] = -1;
2779
2780
2781 /* check corner nodes */
2782 for (i=0; i<CORNERS_OF_ELEM(theElement); i++)
2783 if (theElementContext[i] != NULL)
2784 if(!CORNERTYPE(theElementContext[i]))
2785 {
2786 errortype[i] = NTYPE(theElementContext[i]);
2787 correcttype[i] = CORNER_NODE;
2788 }
2789
2790 /* check mid nodes */
2791 for (i=CORNERS_OF_ELEM(theElement);
2792 i<CORNERS_OF_ELEM(theElement)+EDGES_OF_ELEM(theElement);
2793 i++)
2794 {
2795 if (theElementContext[i] != NULL)
2796 if(NTYPE(theElementContext[i]) != MID_NODE)
2797 {
2798 errortype[i] = NTYPE(theElementContext[i]);
2799 correcttype[i] = MID_NODE;
2800 }
2801 }
2802
2803 #ifdef __THREEDIM__
2804 /* check side nodes */
2805 for (i=CORNERS_OF_ELEM(theElement)+EDGES_OF_ELEM(theElement);
2806 i<CORNERS_OF_ELEM(theElement)+EDGES_OF_ELEM(theElement)+
2807 SIDES_OF_ELEM(theElement);
2808 i++)
2809 {
2810 if (theElementContext[i] != NULL)
2811 if(NTYPE(theElementContext[i]) != SIDE_NODE)
2812 {
2813 errortype[i] = NTYPE(theElementContext[i]);
2814 correcttype[i] = SIDE_NODE;
2815 }
2816 }
2817
2818 #endif
2819
2820 /* check center node */
2821 i = CORNERS_OF_ELEM(theElement)+CENTER_NODE_INDEX(theElement);
2822 if (theElementContext[i] != NULL)
2823 if(NTYPE(theElementContext[i]) != CENTER_NODE)
2824 {
2825 errortype[i] = NTYPE(theElementContext[i]);
2826 correcttype[i] = CENTER_NODE;
2827 }
2828
2829 for (i=0; i<MAX_CORNERS_OF_ELEM+MAX_NEW_CORNERS_DIM; i++)
2830 if (errortype[i] != -1)
2831 {
2832 printf("ERROR: TAG=%d NTYPE(CONTEXT(i=%d)=%d should be %d\n",
2833 TAG(theElement),i,errortype[i],correcttype[i]);
2834 fflush(stdout);
2835 errorflag = 1;
2836 }
2837
2838 ASSERT(errorflag == 0);
2839 }
2840 #endif
2841
2842 /****************************************************************************/
2843 /*
2844 UpdateContext - assemble references
2845
2846 SYNOPSIS:
2847 static int UpdateContext (GRID *theGrid, ELEMENT *theElement, NODE **theElementContext);
2848
2849 PARAMETERS:
2850 . theGrid - grid level of the sons of theElement
2851 . theElement - element to refine
2852 . theContext - context structure to update
2853
2854 DESCRIPTION:
2855 This function assembles references to objects which interact with the sons
2856 of the given element, i.e. objects are allocated, kept or deleted as
2857 indicated by MARK (i) corner nodes (ii) nodes at midpoints of edges
2858
2859 \return <ul>
2860 int
2861 .n 0 - ok
2862 .n 1 - fatal memory error
2863 */
2864 /****************************************************************************/
2865
UpdateContext(GRID * theGrid,ELEMENT * theElement,NODE ** theElementContext)2866 static int UpdateContext (GRID *theGrid, ELEMENT *theElement, NODE **theElementContext)
2867 {
2868 bool toCreate;
2869 #ifdef __THREEDIM__
2870 ELEMENT *theNeighbor; /* neighbor and a son of current elem. */
2871 EDGE *fatherEdge;
2872 INT j;
2873 #endif
2874
2875 /* reset context to NULL */
2876 for(INT i=0; i<MAX_CORNERS_OF_ELEM+MAX_NEW_CORNERS_DIM; i++)
2877 theElementContext[i] = NULL;
2878
2879 /* is element to refine */
2880 if (!MARKED(theElement)) return(GM_OK);
2881
2882 const INT Mark = MARK(theElement);
2883
2884 /* allocate corner nodes if necessary */
2885 for (INT i=0; i<CORNERS_OF_ELEM(theElement); i++)
2886 {
2887 NODE* theNode = CORNER(theElement,i);
2888 if (SONNODE(theNode)==NULL)
2889 {
2890 SONNODE(theNode) = CreateSonNode(theGrid,theNode);
2891 if (SONNODE(theNode)==NULL)
2892 RETURN(GM_FATAL);
2893 #ifdef IDENT_ONLY_NEW
2894 SETNEW_NIDENT(SONNODE(theNode),1);
2895 #endif
2896 }
2897 theElementContext[i] = SONNODE(theNode);
2898 }
2899
2900 /* allocate edge midpoint nodes */
2901 /* nodes on refined edges */
2902 NODE** MidNodes = theElementContext+CORNERS_OF_ELEM(theElement);
2903 for (INT i=0; i<EDGES_OF_ELEM(theElement); i++)
2904 {
2905 const INT Corner0 = CORNER_OF_EDGE(theElement,i,0);
2906 const INT Corner1 = CORNER_OF_EDGE(theElement,i,1);
2907
2908 bool toBisect = false;
2909
2910 if (MARKED_NEW_GREEN(theElement))
2911 {
2912 EDGE* theEdge = GetEdge(CORNER(theElement,Corner0),
2913 CORNER(theElement,Corner1));
2914 ASSERT(theEdge != NULL);
2915
2916 if (ADDPATTERN(theEdge) == 0)
2917 {
2918 toBisect = true;
2919 MidNodes[i] = MIDNODE(theEdge);
2920 }
2921 }
2922 #ifndef __ANISOTROPIC__
2923 else
2924 #endif
2925 if (NODE_OF_RULE(theElement,Mark,i))
2926 {
2927 toBisect = true;
2928 }
2929
2930 IFDEBUG(gm,2)
2931 if (MidNodes[i] == NULL)
2932 UserWriteF("\n MidNodes[%d]: toBisect=%d ID(Corner0)=%d "
2933 "ID(Corner1)=%d",
2934 i,toBisect,ID(CORNER(theElement,Corner0)),
2935 ID(CORNER(theElement,Corner1)));
2936 else
2937 UserWriteF("\n MidNodes[%d]: toBisect=%d ID(Corner0)=%d "
2938 "ID(Corner1)=%d"
2939 " ID(MidNode)=%d",i,toBisect,ID(CORNER(theElement,Corner0)),
2940 ID(CORNER(theElement,Corner1)),ID(MidNodes[i]));
2941 ENDDEBUG
2942
2943 if (toBisect)
2944 {
2945 /* we need a midpoint node */
2946 if (MidNodes[i]!=NULL) continue;
2947 NODE* Node0 = CORNER(theElement,Corner0);
2948 NODE* Node1 = CORNER(theElement,Corner1);
2949 EDGE* theEdge = GetEdge(Node0,Node1);
2950 if (theEdge == nullptr)
2951 RETURN(GM_FATAL);
2952 MidNodes[i] = MIDNODE(theEdge);
2953 /* MidNodes[i] = GetMidNode(theElement,i); */
2954 if (MidNodes[i] == NULL)
2955 {
2956 MidNodes[i] = CreateMidNode(theGrid,theElement,NULL,i);
2957 if (MidNodes[i]==NULL)
2958 RETURN(GM_FATAL);
2959 #ifdef IDENT_ONLY_NEW
2960 SETNEW_NIDENT(MidNodes[i],1);
2961 #endif
2962 IFDEBUG(gm,2)
2963 UserWriteF(" created ID(MidNode)=%d for edge=%d",ID(MidNodes[i]),i);
2964 ENDDEBUG
2965 }
2966 assert(MidNodes[i]!=NULL);
2967 }
2968 }
2969
2970 IFDEBUG(gm,2)
2971 UserWriteF("\n");
2972 ENDDEBUG
2973
2974 #ifdef __THREEDIM__
2975 /* nodes on refined sides */
2976 NODE** SideNodes = theElementContext+CORNERS_OF_ELEM(theElement)+
2977 EDGES_OF_ELEM(theElement);
2978 for (INT i=0; i<SIDES_OF_ELEM(theElement); i++)
2979 {
2980 /* no side nodes for triangular sides yet */
2981 #ifdef DUNE_UGGRID_TET_RULESET
2982 if (CORNERS_OF_SIDE(theElement,i) == 3) continue;
2983 #endif
2984
2985 toCreate = false;
2986 /* is side node needed */
2987 if (MARKED_NEW_GREEN(theElement))
2988 {
2989 theNeighbor = NBELEM(theElement,i);
2990
2991 if (theNeighbor!=NULL)
2992 {
2993 if (MARKCLASS(theNeighbor)!=GREEN_CLASS &&
2994 MARKCLASS(theNeighbor)!=YELLOW_CLASS)
2995 {
2996
2997 for (j=0; j<SIDES_OF_ELEM(theNeighbor); j++)
2998 {
2999 if (NBELEM(theNeighbor,j) == theElement) break;
3000 }
3001 ASSERT(j<SIDES_OF_ELEM(theNeighbor));
3002 if (NODE_OF_RULE(theNeighbor,MARK(theNeighbor),
3003 EDGES_OF_ELEM(theNeighbor)+j))
3004 toCreate = true;
3005 }
3006 }
3007 }
3008 #ifndef __ANISOTROPIC__
3009 else
3010 #endif
3011 if (NODE_OF_RULE(theElement,Mark,EDGES_OF_ELEM(theElement)+i))
3012 {
3013 toCreate = true;
3014 }
3015
3016 IFDEBUG(gm,2)
3017 if (SideNodes[i] == NULL)
3018 UserWriteF(" SideNode[%d]: create=%d old=%x",
3019 i,toCreate,SideNodes[i]);
3020 else
3021 UserWriteF(" SideNode[%d]: create=%d old=%x oldID=%d",i,toCreate,
3022 SideNodes[i],ID(SideNodes[i]));
3023 if (SideNodes[i] != NULL)
3024 if (START(SideNodes[i])!=NULL)
3025 {
3026 LINK *sidelink;
3027 UserWriteF("\n NO_OF_ELEM of EDGES:");
3028 for (sidelink=START(SideNodes[i]);
3029 sidelink!=NULL;
3030 sidelink=NEXT(sidelink))
3031 {
3032 UserWriteF(" NO=%d NodeTo=%d",NO_OF_ELEM(MYEDGE(sidelink)),
3033 ID(NBNODE(sidelink)));
3034 }
3035 UserWrite("\n");
3036 }
3037 ENDDEBUG
3038
3039 if (toCreate)
3040 {
3041
3042 theNeighbor = NBELEM(theElement,i);
3043
3044 IFDEBUG(gm,1)
3045 if (theNeighbor != NULL)
3046 {
3047 IFDEBUG(gm,3)
3048 UserWriteF(" ID(theNeighbor)=%d nbadr=%x:\n",
3049 ID(theNeighbor),theNeighbor);
3050 ENDDEBUG
3051 }
3052 else
3053 {
3054 /* this must be a boundary side */
3055 ASSERT(SIDE_ON_BND(theElement,i));
3056 }
3057 ENDDEBUG
3058
3059 /*
3060 if (MARKED_NEW_GREEN(theNeighbor) && USED(theNeighbor)==0)
3061 */
3062 if (theNeighbor !=NULL)
3063 {
3064
3065 IFDEBUG(gm,3)
3066 UserWriteF(" Searching for side node already allocated:\n");
3067 ENDDEBUG
3068
3069 /* check for side node */
3070 SideNodes[i] = GetSideNode(theElement,i);
3071 }
3072
3073 if (SideNodes[i] == NULL)
3074 {
3075
3076 /* allocate the sidenode */
3077 if ((SideNodes[i] = CreateSideNode(theGrid,theElement,NULL,i)) == NULL)
3078 RETURN(GM_FATAL);
3079
3080 #ifdef IDENT_ONLY_NEW
3081 SETNEW_NIDENT(SideNodes[i],1);
3082 #endif
3083 }
3084
3085 IFDEBUG(gm,0)
3086 ASSERT(SideNodes[i]!=NULL);
3087 for (j=0; j<EDGES_OF_SIDE(theElement,i); j++)
3088 {
3089
3090 fatherEdge = GetEdge(CORNER_OF_EDGE_PTR(theElement,EDGE_OF_SIDE(theElement,i,j),0),
3091 CORNER_OF_EDGE_PTR(theElement,EDGE_OF_SIDE(theElement,i,j),1));
3092
3093 [[maybe_unused]] NODE* Node0 = MIDNODE(fatherEdge);
3094
3095 /* if side node exists all mid nodes must exist */
3096 ASSERT(Node0 != NULL);
3097 }
3098 ENDDEBUG
3099 }
3100
3101 IFDEBUG(gm,2)
3102 if (SideNodes[i] != NULL)
3103 UserWriteF(" new=%x newID=%d\n",SideNodes[i],ID(SideNodes[i]));
3104 else
3105 UserWriteF(" new=%x\n",SideNodes[i]);
3106 ENDDEBUG
3107 }
3108 #endif
3109
3110 /* allocate center node */
3111 NODE** CenterNode = MidNodes+CENTER_NODE_INDEX(theElement);
3112 CenterNode[0] = NULL;
3113
3114 toCreate = false;
3115 if (CenterNode[0] == NULL)
3116 {
3117
3118 if (MARKED_NEW_GREEN(theElement))
3119 {
3120 toCreate = true;
3121 }
3122 #ifndef __ANISOTROPIC__
3123 else
3124 #endif
3125 if (NODE_OF_RULE(theElement,Mark,CENTER_NODE_INDEX(theElement)))
3126 {
3127 toCreate = true;
3128 }
3129 }
3130
3131 IFDEBUG(gm,2)
3132 if (CenterNode[0] == NULL)
3133 UserWriteF(" CenterNode: create=%d old=%x",toCreate,CenterNode[0]);
3134 else
3135 UserWriteF(" CenterNode: create=%d old=%x oldID=%d",toCreate,
3136 CenterNode[0],ID(CenterNode[0]));
3137 ENDDEBUG
3138
3139 if (toCreate)
3140 {
3141 if ((CenterNode[0] = CreateCenterNode(theGrid,theElement,NULL)) == NULL)
3142 RETURN(GM_FATAL);
3143 }
3144
3145 #ifdef ModelP
3146 /* mark nodes as needed */
3147 for(INT i=0; i<MAX_CORNERS_OF_ELEM+MAX_NEW_CORNERS_DIM; i++)
3148 {
3149 if (theElementContext[i] != NULL) SETUSED(theElementContext[i],1);
3150 }
3151 #endif
3152
3153 IFDEBUG(gm,2)
3154 if (CenterNode[0] != NULL)
3155 UserWriteF(" new=%x newID=%d\n",CenterNode[0],ID(CenterNode[0]));
3156 else
3157 UserWriteF(" new=%x\n",CenterNode[0]);
3158 ENDDEBUG
3159
3160 return(GM_OK);
3161 }
3162
3163
3164 /****************************************************************************/
3165 /*
3166 UnrefineElement - remove previous refinement
3167
3168 SYNOPSIS:
3169 static INT UnrefineElement (GRID *theGrid, ELEMENT *theElement);
3170
3171 PARAMETERS:
3172 . theGrid: grid level of sons of theElement
3173 . theElement: element to refine
3174
3175 DESCRIPTION:
3176 This function removes previous refinement for an element and all sonelement
3177 recursively, deletes: (i) all connections, (ii) all interior nodes and edges
3178 are deleted, (iii) sons are deleted and references to sons reset to NULL.
3179
3180 \return <ul>
3181 INT
3182 */
3183 /****************************************************************************/
3184
UnrefineElement(GRID * theGrid,ELEMENT * theElement)3185 static INT UnrefineElement (GRID *theGrid, ELEMENT *theElement)
3186 {
3187 int s;
3188 ELEMENT *theSon,*SonList[MAX_SONS];
3189 [[maybe_unused]] const int me = theGrid->ppifContext().me();
3190
3191 /* something to do ? */
3192 if ((REFINE(theElement)==NO_REFINEMENT)||(theGrid==NULL)) return(GM_OK);
3193
3194 if (GetAllSons(theElement,SonList)!=GM_OK) RETURN(GM_FATAL);
3195
3196 for (s=0; SonList[s]!=NULL; s++)
3197 {
3198 theSon = SonList[s];
3199 SETMARK(theSon,NO_REFINEMENT);
3200 if (IS_REFINED(theSon))
3201 {
3202 if (UnrefineElement(UPGRID(theGrid),theSon)) RETURN(GM_FATAL);
3203 }
3204 }
3205
3206 /* remove connections in neighborhood of sons */
3207 for (s=0; SonList[s]!=NULL; s++)
3208 {
3209 DisposeConnectionsInNeighborhood(theGrid,SonList[s]);
3210 }
3211
3212 /* remove son elements */
3213 #ifndef ModelP
3214 IFDEBUG(gm,1)
3215 if (!REFINED_NEW_GREEN(theElement))
3216 {
3217 if (NSONS(theElement) != NSONS_OF_RULE(MARK2RULEADR(theElement,
3218 REFINE(theElement))))
3219 {
3220 UserWriteF("ERROR: NSONS=%d but rule.sons=%d\n",NSONS(theElement),
3221 NSONS_OF_RULE(MARK2RULEADR(theElement,REFINE(theElement))));
3222 }
3223 }
3224 ENDDEBUG
3225 #endif
3226
3227 for (s=0; SonList[s]!=NULL; s++)
3228 {
3229 /** \todo delete special debug */
3230 /* if (ID(SonList[s])==11644) { RETURN(GM_FATAL); } */
3231 PRINTDEBUG(gm,1,(PFMT "UnrefineElement(): DisposeElement[%d]="
3232 EID_FMTX "\n",me,s,EID_PRTX(SonList[s])));
3233
3234 if (DisposeElement(theGrid,SonList[s],true)!=0) RETURN(GM_FATAL);
3235 }
3236
3237 return (GM_OK);
3238 }
3239
3240
3241 START_UGDIM_NAMESPACE
3242
3243 struct compare_record
3244 {
3245 /** \brief Element to connect */
3246 ELEMENT *elem;
3247
3248 /** \brief Side of element to connect */
3249 INT side;
3250
3251 /** \brief Number of nodes of side */
3252 INT nodes;
3253
3254 /** \brief Number of nodes in descending order */
3255 NODE *nodeptr[4];
3256 };
3257 typedef struct compare_record COMPARE_RECORD;
3258
3259 END_UGDIM_NAMESPACE
3260
GetSonSideNodes(const ELEMENT * theElement,INT side,INT * nodes,NODE * SideNodes[MAX_SIDE_NODES],INT ioflag)3261 INT NS_DIM_PREFIX GetSonSideNodes (const ELEMENT *theElement, INT side, INT *nodes,
3262 NODE *SideNodes[MAX_SIDE_NODES], INT ioflag)
3263 {
3264 INT i,ncorners,nedges;
3265
3266 ncorners = CORNERS_OF_SIDE(theElement,side);
3267 nedges = EDGES_OF_SIDE(theElement,side);
3268 (*nodes) = 0;
3269
3270
3271 /* reset pointers */
3272 for (i=0; i<MAX_SIDE_NODES; i++)
3273 {
3274 SideNodes[i] = NULL;
3275 }
3276
3277 /* determine corner nodes */
3278 for (i=0; i<ncorners; i++)
3279 {
3280 SideNodes[i] = SONNODE(CORNER_OF_SIDE_PTR(theElement,side,i));
3281 #ifndef ModelP
3282 assert(SideNodes[i]!=NULL);
3283 #endif
3284 if (!ioflag)
3285 assert(SideNodes[i]==NULL || CORNERTYPE(SideNodes[i]));
3286 (*nodes)++;
3287 }
3288
3289 /* determine mid nodes */
3290 for (i=0; i<nedges; i++)
3291 {
3292 /** \todo delete
3293 #ifdef __TWODIM__
3294 ASSERT(OBJT(NFATHER(SideNodes[i])) == NDOBJ);
3295 ASSERT(OBJT(NFATHER(SideNodes[i+1])) == NDOBJ);
3296 theEdge = GetEdge((NODE *)NFATHER(SideNodes[i]),
3297 (NODE *)NFATHER(SideNodes[i+1]));
3298 #endif
3299 #ifdef __THREEDIM__
3300 ASSERT(OBJT(NFATHER(SideNodes[i])) == NDOBJ);
3301 ASSERT(OBJT(NFATHER(SideNodes[(i+1)%nedges])) == NDOBJ);
3302 theEdge = GetEdge((NODE *)NFATHER(SideNodes[i]),
3303 (NODE *)NFATHER(SideNodes[(i+1)%nedges]));
3304 #endif
3305 assert(theEdge != NULL);
3306
3307 IFDEBUG(gm,4)
3308 UserWriteF("theEdge=%x midnode=%x\n",theEdge,MIDNODE(theEdge));
3309 ENDDEBUG
3310
3311 if (MIDNODE(theEdge) != NULL)
3312 {
3313 SideNodes[ncorners+i] = MIDNODE(theEdge);
3314 assert(NTYPE(MIDNODE(theEdge)) == MID_NODE);
3315 (*nodes)++;
3316 }
3317 */
3318
3319 SideNodes[ncorners+i] = GetMidNode(theElement,EDGE_OF_SIDE(theElement,side,i));
3320 if (SideNodes[ncorners+i] != NULL)
3321 {
3322 assert(NTYPE(SideNodes[ncorners+i]) == MID_NODE);
3323 (*nodes)++;
3324 }
3325 }
3326
3327 #ifdef __THREEDIM__
3328 /* determine side node */
3329 {
3330 NODE *theNode;
3331
3332 theNode = GetSideNode(theElement,side);
3333 if (theNode != NULL)
3334 {
3335 (*nodes)++;
3336 }
3337
3338 SideNodes[ncorners+nedges] = theNode;
3339
3340 IFDEBUG(gm,4)
3341 UserWriteF("sidenode=%x\n",theNode);
3342 ENDDEBUG
3343 }
3344 #endif
3345
3346 IFDEBUG(gm,2)
3347 UserWriteF("GetSonSideNodes\n");
3348 for (i=0; i<MAX_SIDE_NODES; i++) UserWriteF(" %5d",i);
3349 UserWriteF("\n");
3350 for (i=0; i<MAX_SIDE_NODES; i++)
3351 if (SideNodes[i]!=NULL) UserWriteF(" %5d",ID(SideNodes[i]));
3352 UserWriteF("\n");
3353 ENDDEBUG
3354
3355 return(GM_OK);
3356 }
3357
3358
3359 /****************************************************************************/
3360 /*
3361 compare_node -
3362
3363 SYNOPSIS:
3364 static bool compare_node (const NODE* a, const NODE* b)
3365
3366 PARAMETERS:
3367 . a
3368 . b
3369
3370 DESCRIPTION:
3371
3372 \return <ul>
3373 bool
3374 */
3375 /****************************************************************************/
3376
compare_node(const NODE * a,const NODE * b)3377 static bool compare_node(const NODE* a, const NODE* b)
3378 {
3379 return a > b;
3380 }
3381
3382 /****************************************************************************/
3383 /** \brief Get the sons of an element side
3384
3385 \param[in] theElement Input element
3386 \param[in] side Input element side number
3387 \param[out] Sons_of_Side Number of topological sons of the element side
3388 \param[in,out] SonList[MAX_SONS] Output elements
3389 \param[out] SonSides Output element side numbers
3390 \param[in] NeedSons If this is false, the correct list of sons is expected to be
3391 provided in SonList. If not it is recomputed.
3392 \param[in] ioflag An obsolete debugging flag
3393 \param[in] useRefineClass
3394
3395 For a given side of an element, this routine computes all element sides
3396 on the next finer grid level which are topological sons of the input
3397 element side.
3398 */
3399 /****************************************************************************/
3400
Get_Sons_of_ElementSide(const ELEMENT * theElement,INT side,INT * Sons_of_Side,ELEMENT * SonList[MAX_SONS],INT * SonSides,INT NeedSons,INT ioflag,INT useRefineClass)3401 INT NS_DIM_PREFIX Get_Sons_of_ElementSide (const ELEMENT *theElement, INT side, INT *Sons_of_Side,
3402 ELEMENT *SonList[MAX_SONS], INT *SonSides,
3403 INT NeedSons, INT ioflag,
3404 INT useRefineClass)
3405 {
3406 INT i,j,nsons;
3407 enum MarkClass markclass;
3408
3409 /* reset soncount */
3410 *Sons_of_Side = 0;
3411 nsons = 0;
3412
3413 /* get sons of element */
3414 if (NeedSons)
3415 if (GetAllSons(theElement,SonList) != GM_OK) RETURN(GM_FATAL);
3416
3417 IFDEBUG(gm,2)
3418 UserWriteF(" Get_Sons_of_ElementSide():"
3419 " id=%d tag=%d, refineclass=%d markclass=%d refine=%d mark=%d coarse=%d"
3420 " used=%d nsons=%d side=%d needsons=%d\n",
3421 ID(theElement),TAG(theElement),REFINECLASS(theElement),MARKCLASS(theElement),
3422 REFINE(theElement),MARK(theElement),COARSEN(theElement),
3423 USED(theElement),NSONS(theElement),side,NeedSons);
3424 for (i=0; SonList[i]!=NULL; i++)
3425 UserWriteF(" son[%d]=" EID_FMTX "\n",i,EID_PRTX(SonList[i]));
3426 ENDDEBUG
3427
3428 #ifdef __TWODIM__
3429 markclass = RED_CLASS;
3430 #endif
3431 #ifdef __THREEDIM__
3432 /* The following line used to read: markclass = (enum MarkClass) MARKCLASS(theElement);
3433 This works well within the UG grid refinement context. However, now I want to use this
3434 method within the DUNE UGGridLeafIntersectionIterator. The problem is that the user
3435 may have randomly marked elements before calling the iterator. In that case the
3436 mark classes may not be set the way this method expects it. Hence we allow the option
3437 to use REFINECLASS instead. To be absolutely certain we don't break existing code
3438 we keep the old behaviour as the default.*/
3439 markclass = (enum MarkClass) ((useRefineClass) ? REFINECLASS(theElement) : MARKCLASS(theElement));
3440 #endif
3441
3442 /** \todo quick fix */
3443 #ifdef ModelP
3444 if (EHGHOST(theElement))
3445 markclass = GREEN_CLASS;
3446 #endif
3447
3448 /* select sons to connect */
3449 switch (markclass)
3450 {
3451
3452 case YELLOW_CLASS :
3453 {
3454 *Sons_of_Side = 1;
3455 SonSides[0] = side;
3456 break;
3457 }
3458
3459 case GREEN_CLASS :
3460 case RED_CLASS :
3461 {
3462 /* determine sonnodes of side */
3463 NODE *SideNodes[MAX_SIDE_NODES];
3464 INT corner[MAX_CORNERS_OF_SIDE];
3465 INT n,nodes;
3466
3467 /* determine nodes of sons on side of element */
3468 GetSonSideNodes(theElement,side,&nodes,SideNodes,ioflag);
3469
3470 /* sort side nodes in descending adress order */
3471 std::sort(SideNodes, SideNodes + MAX_SIDE_NODES, compare_node);
3472
3473 IFDEBUG(gm,3)
3474 UserWriteF("After qsort:\n");
3475 for (i=0; i<MAX_SIDE_NODES; i++) UserWriteF(" %8d",i);
3476 UserWriteF("\n");
3477 for (i=0; i<MAX_SIDE_NODES; i++)
3478 if (SideNodes[i]!=NULL) UserWriteF(" %x",SideNodes[i]);
3479 else UserWriteF(" %8d",0);
3480 UserWriteF("\n");
3481 ENDDEBUG
3482
3483 /* determine sonnode on side */
3484 /* for (i=0; i<NSONS(theElement); i++) */
3485 for (i=0; SonList[i]!=NULL; i++)
3486 {
3487 n = 0;
3488
3489 for (j=0; j<MAX_CORNERS_OF_SIDE; j++)
3490 corner[j] = -1;
3491
3492 IFDEBUG(gm,4)
3493 UserWriteF("son=%d\n",i);
3494 ENDDEBUG
3495
3496 /* soncorners on side */
3497 for (j=0; j<CORNERS_OF_ELEM(SonList[i]); j++)
3498 {
3499 NODE *nd = CORNER(SonList[i],j);
3500 if (std::binary_search(SideNodes, SideNodes + nodes, nd, compare_node))
3501 {
3502 corner[n] = j;
3503 n++;
3504 }
3505 }
3506 assert(n<5);
3507
3508 IFDEBUG(gm,4)
3509 UserWriteF("\n nodes on side n=%d:",n);
3510 for (j=0; j<MAX_CORNERS_OF_SIDE; j++)
3511 UserWriteF(" %d",corner[j]);
3512 ENDDEBUG
3513
3514
3515 IFDEBUG(gm,0)
3516 if (n==3) assert(TAG(SonList[i])!=HEXAHEDRON);
3517 if (n==4) assert(TAG(SonList[i])!=TETRAHEDRON);
3518 ENDDEBUG
3519
3520 /* sonside on side */
3521 #ifdef __TWODIM__
3522 assert(n<=2);
3523 if (n==2)
3524 {
3525 if (corner[0]+1 == corner[1])
3526 SonSides[nsons] = corner[0];
3527 else
3528 {
3529 /** \todo find proper assert
3530 assert(corner[1] == CORNERS_OF_ELEM(theElement)-1); */
3531 SonSides[nsons] = corner[1];
3532 }
3533 SonList[nsons] = SonList[i];
3534 nsons++;
3535 }
3536 #endif
3537 #ifdef __THREEDIM__
3538 if (n==3 || n==4)
3539 {
3540 INT edge0,edge1,sonside,side0,side1;
3541
3542 /* determine side number */
3543 edge0 = edge1 = -1;
3544 edge0 = EDGE_WITH_CORNERS(SonList[i],corner[0],corner[1]);
3545 edge1 = EDGE_WITH_CORNERS(SonList[i],corner[1],corner[2]);
3546 /* corners are not stored in local side numbering, */
3547 /* therefore corner[x]-corner[y] might be the diagonal */
3548 if (n==4 && edge0==-1)
3549 edge0 = EDGE_WITH_CORNERS(SonList[i],corner[0],
3550 corner[3]);
3551 if (n==4 && edge1==-1)
3552 edge1 = EDGE_WITH_CORNERS(SonList[i],corner[1],
3553 corner[3]);
3554 assert(edge0!=-1 && edge1!=-1);
3555
3556 sonside = -1;
3557 for (side0=0; side0<MAX_SIDES_OF_EDGE; side0++)
3558 {
3559 for (side1=0; side1<MAX_SIDES_OF_EDGE; side1++)
3560 {
3561 IFDEBUG(gm,5)
3562 UserWriteF("edge0=%d side0=%d SIDE_WITH_EDGE=%d\n",
3563 edge0, side0,
3564 SIDE_WITH_EDGE(SonList[i],edge0,side0));
3565 UserWriteF("edge1=%d side1=%d SIDE_WITH_EDGE=%d\n",
3566 edge1, side1,
3567 SIDE_WITH_EDGE(SonList[i],edge1,side1));
3568 ENDDEBUG
3569 if (SIDE_WITH_EDGE(SonList[i],edge0,side0) ==
3570 SIDE_WITH_EDGE(SonList[i],edge1,side1))
3571 {
3572 sonside = SIDE_WITH_EDGE(SonList[i],edge0,side0);
3573 break;
3574 }
3575 }
3576 if (sonside != -1) break;
3577 }
3578 assert(sonside != -1);
3579 IFDEBUG(gm,4)
3580 UserWriteF(" son[%d]=%x with sonside=%d on eside=%d\n",
3581 i,SonList[i],sonside,side);
3582 ENDDEBUG
3583
3584 IFDEBUG(gm,3)
3585 INT k;
3586 ELEMENT *Nb;
3587
3588 for (k=0; k<SIDES_OF_ELEM(SonList[i]); k++)
3589 {
3590 Nb = NBELEM(SonList[i],k);
3591 if (Nb!=NULL)
3592 {
3593 INT j;
3594 for (j=0; j<SIDES_OF_ELEM(Nb); j++)
3595 {
3596 if (NBELEM(Nb,j)==SonList[i]) break;
3597 }
3598 if (j<SIDES_OF_ELEM(Nb))
3599 UserWriteF(" sonside=%d has backptr to son "
3600 "Nb=%x Nbside=%d\n",k,Nb,j);
3601 }
3602 }
3603 ENDDEBUG
3604
3605 ASSERT(CORNERS_OF_SIDE(SonList[i],sonside) == n);
3606
3607 SonSides[nsons] = sonside;
3608 SonList[nsons] = SonList[i];
3609 nsons++;
3610 }
3611 #endif
3612 }
3613 #ifndef ModelP
3614 assert(nsons>0 && nsons<6);
3615 #endif
3616
3617 IFDEBUG(gm,3)
3618 UserWriteF(" nsons on side=%d\n",nsons);
3619 ENDDEBUG
3620
3621 *Sons_of_Side = nsons;
3622 break;
3623 }
3624
3625 /* old style
3626 This case would work for the sequential case if SonList is ordered according to
3627 rule (e.g. introduce GetOrderedSonList()). The upper case, used now in each
3628 situation, is a factor of 2
3629 slower for uniform refinement!! (s.l. 981016)
3630 #ifndef ModelP
3631 case RED_CLASS:
3632 #endif
3633 */
3634 {
3635 SONDATA *sondata;
3636
3637 for (i=0; SonList[i]!=NULL; i++)
3638 {
3639 sondata = SON_OF_RULE(MARK2RULEADR(theElement,
3640 MARK(theElement)),i);
3641
3642 for (j=0; j<SIDES_OF_ELEM(SonList[i]); j++)
3643 if (SON_NB(sondata,j) == FATHER_SIDE_OFFSET+side)
3644 {
3645 SonSides[nsons] = j;
3646 SonList[nsons] = SonList[i];
3647 nsons ++;
3648 }
3649 }
3650 *Sons_of_Side = nsons;
3651 break;
3652 }
3653
3654 default :
3655 RETURN(GM_FATAL);
3656 }
3657
3658 #ifdef ModelP
3659 IFDEBUG(gm,4)
3660 UserWriteF("Sons_of_Side=%d\n",*Sons_of_Side);
3661 for (i=0; i<*Sons_of_Side; i++)
3662 UserWriteF("son[%d]=" EID_FMTX " sonside[%d]=%d\n",i,
3663 EID_PRTX(SonList[i]),i,SonSides[i]);
3664 ENDDEBUG
3665 #endif
3666
3667 for (i=*Sons_of_Side; i<MAX_SONS; i++)
3668 SonList[i] = NULL;
3669
3670 return(GM_OK);
3671 }
3672
3673
3674 /****************************************************************************/
3675 /*
3676 Sort_Node_Ptr -
3677
3678 SYNOPSIS:
3679 static INT Sort_Node_Ptr (INT n,NODE **nodes);
3680
3681 PARAMETERS:
3682 . n
3683 . nodes
3684
3685 DESCRIPTION:
3686
3687 \return <ul>
3688 INT
3689 */
3690 /****************************************************************************/
3691
Sort_Node_Ptr(INT n,NODE ** nodes)3692 static INT Sort_Node_Ptr (INT n,NODE **nodes)
3693 {
3694 NODE* nd;
3695 INT i,j,max;
3696
3697 max = 0;
3698
3699 switch (n)
3700 {
3701
3702 #ifdef __TWODIM__
3703 case 2 :
3704 #endif
3705 #ifdef __THREEDIM__
3706 case 3 :
3707 case 4 :
3708 #endif
3709 for (i=0; i<n; i++)
3710 {
3711 max = i;
3712 for (j=i+1; j<n; j++)
3713 if (nodes[max]<nodes[j]) max = j;
3714 if (i != max)
3715 {
3716 nd = nodes[i];
3717 nodes[i] = nodes[max];
3718 nodes[max] = nd;
3719 }
3720 }
3721 break;
3722
3723 default :
3724 RETURN(GM_FATAL);
3725 }
3726
3727 return(GM_OK);
3728 }
3729
3730 /****************************************************************************/
3731 /*
3732 Fill_Comp_Table -
3733
3734 SYNOPSIS:
3735 static INT Fill_Comp_Table (COMPARE_RECORD **SortTable, COMPARE_RECORD *Table, INT nelems, ELEMENT **Elements, INT *Sides);
3736
3737 PARAMETERS:
3738 . SortTable
3739 . Table
3740 . nelems
3741 . Elements
3742 . Sides
3743
3744 DESCRIPTION:
3745
3746 \return <ul>
3747 INT
3748 */
3749 /****************************************************************************/
3750
Fill_Comp_Table(COMPARE_RECORD ** SortTable,COMPARE_RECORD * Table,INT nelems,ELEMENT ** Elements,INT * Sides)3751 static INT Fill_Comp_Table (COMPARE_RECORD **SortTable, COMPARE_RECORD *Table, INT nelems,
3752 ELEMENT **Elements, INT *Sides)
3753 {
3754 COMPARE_RECORD *Entry;
3755 INT i,j;
3756
3757 for (i=0; i<nelems; i++)
3758 {
3759 SortTable[i] = Table+i;
3760 Entry = Table+i;
3761 Entry->elem = Elements[i];
3762 Entry->side = Sides[i];
3763 Entry->nodes = CORNERS_OF_SIDE(Entry->elem,Entry->side);
3764 for (j=0; j<CORNERS_OF_SIDE(Entry->elem,Entry->side); j++)
3765 Entry->nodeptr[j] = CORNER_OF_SIDE_PTR(Entry->elem,Entry->side,j);
3766 if (Sort_Node_Ptr(Entry->nodes,Entry->nodeptr)!=GM_OK) RETURN(GM_FATAL);
3767 }
3768
3769 return(GM_OK);
3770 }
3771
3772
3773 /****************************************************************************/
3774 /*
3775 compare_nodes -
3776
3777 SYNOPSIS:
3778 static bool compare_nodes(const COMPARE_RECORD* a, const COMPARE_RECORD* b)
3779
3780 PARAMETERS:
3781 . a
3782 . b
3783
3784 DESCRIPTION:
3785
3786 \return <ul>
3787 bool
3788 */
3789 /****************************************************************************/
3790
compare_nodes(const COMPARE_RECORD * a,const COMPARE_RECORD * b)3791 static bool compare_nodes(const COMPARE_RECORD* a, const COMPARE_RECORD* b)
3792 {
3793 const int n = (a->nodes == 4 && b->nodes == 4) ? 4 : 3;
3794
3795 for (int i = 0; i < n; ++i) {
3796 if (a->nodeptr[i] > b->nodeptr[i])
3797 return true;
3798 else if (a->nodeptr[i] < b->nodeptr[i])
3799 return false;
3800 }
3801
3802 return false;
3803 }
3804
3805 /****************************************************************************/
3806 /*
3807 Connect_Sons_of_ElementSide -
3808
3809 SYNOPSIS:
3810 INT Connect_Sons_of_ElementSide (GRID *theGrid, ELEMENT *theElement, INT side, INT Sons_of_Side, ELEMENT **Sons_of_Side_List, INT *SonSides, INT ioflag);
3811
3812 PARAMETERS:
3813 . theGrid
3814 . theElement
3815 . side
3816 . Sons_of_Side
3817 . Sons_of_Side_List
3818 . SonSides
3819 . ioflag
3820
3821 DESCRIPTION:
3822
3823 \return <ul>
3824 INT
3825 */
3826 /****************************************************************************/
3827
Connect_Sons_of_ElementSide(GRID * theGrid,ELEMENT * theElement,INT side,INT Sons_of_Side,ELEMENT ** Sons_of_Side_List,INT * SonSides,INT ioflag)3828 INT NS_DIM_PREFIX Connect_Sons_of_ElementSide (GRID *theGrid, ELEMENT *theElement, INT side, INT Sons_of_Side, ELEMENT **Sons_of_Side_List, INT *SonSides, INT ioflag)
3829 {
3830 COMPARE_RECORD ElemSonTable[MAX_SONS];
3831 COMPARE_RECORD NbSonTable[MAX_SONS];
3832 COMPARE_RECORD *ElemSortTable[MAX_SONS];
3833 COMPARE_RECORD *NbSortTable[MAX_SONS];
3834
3835 ELEMENT *theNeighbor;
3836 ELEMENT *Sons_of_NbSide_List[MAX_SONS];
3837 INT nbside,Sons_of_NbSide,NbSonSides[MAX_SONS];
3838 INT i,j,k;
3839
3840 IFDEBUG(gm,2)
3841 UserWriteF("Connect_Sons_of_ElementSide: ID(elem)=%d side=%d "
3842 "Sons_of_Side=%d\n",ID(theElement),side,Sons_of_Side);
3843 REFINE_ELEMENT_LIST(0,theElement,"theElement:");
3844 ENDDEBUG
3845
3846 if (Sons_of_Side <= 0) return(GM_OK);
3847
3848 /* connect to boundary */
3849 if (OBJT(theElement)==BEOBJ && SIDE_ON_BND(theElement,side))
3850 {
3851 /** \todo connect change test */
3852
3853 for (i=0; i<Sons_of_Side; i++)
3854 {
3855
3856 assert(OBJT(Sons_of_Side_List[i])==BEOBJ);
3857 if (CreateSonElementSide(theGrid,theElement,side,
3858 Sons_of_Side_List[i],SonSides[i]) != GM_OK)
3859 {
3860 return(GM_FATAL);
3861 }
3862 }
3863
3864 /* internal boundaries not connected */
3865 /* return(GM_OK); */
3866 }
3867
3868 /* connect to neighbor element */
3869 theNeighbor = NBELEM(theElement,side);
3870
3871 if (theNeighbor==NULL) return(GM_OK);
3872
3873 /* master elements only connect to master elements */
3874 /* ghost elements connect to ghost and master elements */
3875 #ifdef ModelP
3876 if (!ioflag && EMASTER(theElement) && EHGHOST(theNeighbor))
3877 return(GM_OK);
3878 #endif
3879
3880 /* only yellow elements may have no neighbors */
3881 if (MARKCLASS(theNeighbor)==NO_CLASS)
3882 {
3883
3884 if (hFlag) assert(MARKCLASS(theElement)==YELLOW_CLASS);
3885
3886 return(GM_OK);
3887 }
3888
3889 if (REFINEMENT_CHANGES(theNeighbor)) return(GM_OK);
3890
3891 /* determine corresponding side of neighbor */
3892 for (nbside=0; nbside<SIDES_OF_ELEM(theNeighbor); nbside++)
3893 if (NBELEM(theNeighbor,nbside) == theElement) break;
3894 assert(nbside<SIDES_OF_ELEM(theNeighbor));
3895
3896 /* get sons of neighbor to connect */
3897 Get_Sons_of_ElementSide(theNeighbor,nbside,&Sons_of_NbSide,
3898 Sons_of_NbSide_List,NbSonSides,1,ioflag);
3899
3900 #ifdef ModelP
3901 if (!ioflag)
3902 {
3903 #endif
3904 /* match exactly */
3905 if (1) ASSERT(Sons_of_Side == Sons_of_NbSide && Sons_of_NbSide>0
3906 && Sons_of_NbSide<6);
3907 else
3908 if (!(Sons_of_Side == Sons_of_NbSide && Sons_of_NbSide>0
3909 && Sons_of_NbSide<6))
3910 {
3911 INT i;
3912 ELEMENT *elem=NULL;
3913 ELEMENT *SonList[MAX_SONS];
3914
3915 REFINE_ELEMENT_LIST(0,theElement,"theElement:");
3916 REFINE_ELEMENT_LIST(0,theNeighbor,"theNeighbor:");
3917 UserWriteF("elem=" EID_FMTX " nb=" EID_FMTX
3918 "Sons_of_Side=%d Sons_of_NbSide=%d\n",EID_PRTX(theElement),
3919 EID_PRTX(theNeighbor),Sons_of_Side,Sons_of_NbSide);
3920 fflush(stdout);
3921 GetAllSons(theElement,SonList);
3922 for (i=0; SonList[i]!=NULL; i++)
3923 {
3924 REFINE_ELEMENT_LIST(0,SonList[i],"son:");
3925 }
3926 GetAllSons(theNeighbor,SonList);
3927 for (i=0; SonList[i]!=NULL; i++)
3928 {
3929 REFINE_ELEMENT_LIST(0,SonList[i],"nbson:");
3930 }
3931
3932 /* sig bus error to see stack in totalview */
3933 SET_NBELEM(elem,0,NULL);
3934 assert(0);
3935
3936 }
3937 #ifdef ModelP
3938 }
3939 #endif
3940
3941 IFDEBUG(gm,2)
3942 UserWriteF("Connect_Sons_of_ElementSide: NBID(elem)=%d side=%d "
3943 "Sons_of_Side=%d\n",ID(theNeighbor),nbside,Sons_of_NbSide);
3944 REFINE_ELEMENT_LIST(0,theNeighbor,"theNeighbor:");
3945 ENDDEBUG
3946
3947 /* fill sort and comparison tables */
3948 Fill_Comp_Table(ElemSortTable,ElemSonTable,Sons_of_Side,Sons_of_Side_List,
3949 SonSides);
3950 Fill_Comp_Table(NbSortTable,NbSonTable,Sons_of_NbSide,Sons_of_NbSide_List,
3951 NbSonSides);
3952
3953 IFDEBUG(gm,5)
3954 INT i,j;
3955
3956 if (!ioflag)
3957 {
3958 UserWriteF("BEFORE qsort\n");
3959
3960 /* test whether all entries are corresponding */
3961 for (i=0; i<Sons_of_Side; i++)
3962 {
3963 COMPARE_RECORD *Entry, *NbEntry;
3964
3965 Entry = ElemSortTable[i];
3966 NbEntry = NbSortTable[i];
3967
3968 if (Entry->nodes != NbEntry->nodes)
3969 UserWriteF("Connect_Sons_of_ElementSide(): LIST Sorttables[%d]"
3970 " eNodes=%d nbNodes=%d\n",
3971 i,Entry->nodes,NbEntry->nodes);
3972 for (j=0; j<Entry->nodes; j++)
3973 UserWriteF("Connect_Sons_of_ElementSide(): LIST Sorttables[%d][%d]"
3974 " eNodePtr=%d/%8x/%d nbNodePtr=%d/%8x/%d\n",
3975 i,j,
3976 ID(Entry->nodeptr[j]),Entry->nodeptr[j],NTYPE(Entry->nodeptr[j]),
3977 ID(NbEntry->nodeptr[j]),NbEntry->nodeptr[j],NTYPE(NbEntry->nodeptr[j]));
3978 UserWriteF("\n");
3979 }
3980 UserWriteF("\n\n");
3981 }
3982 ENDDEBUG
3983
3984 /* qsort the tables using nodeptrs */
3985 std::sort(ElemSortTable, ElemSortTable + Sons_of_Side, compare_nodes);
3986 std::sort(NbSortTable, NbSortTable + Sons_of_NbSide, compare_nodes);
3987
3988 #ifdef ModelP
3989 if (!ioflag && Sons_of_NbSide!=Sons_of_Side) ASSERT(0);
3990 #endif
3991
3992 #ifdef Debug
3993 if (!ioflag)
3994 /* check whether both sort table match exactly */
3995 for (i=0; i<Sons_of_Side; i++)
3996 {
3997 COMPARE_RECORD *Entry, *NbEntry;
3998 INT j;
3999
4000 Entry = ElemSortTable[i];
4001 NbEntry = NbSortTable[i];
4002
4003 if (Entry->nodes != NbEntry->nodes)
4004 {
4005 printf("Connect_Sons_of_ElementSide(): ERROR Sorttables[%d]"\
4006 " eNodes=%d nbNodes=%d\n",i,Entry->nodes,NbEntry->nodes);
4007 assert(0);
4008 }
4009 for (j=0; j<Entry->nodes; j++)
4010 if (Entry->nodeptr[j] != NbEntry->nodeptr[j])
4011 {
4012 printf("Connect_Sons_of_ElementSide(): "
4013 "ERROR Sorttables[%d][%d]"\
4014 " eNodePtr=%p nbNodePtr=%p\n",
4015 i,j,Entry->nodeptr[j],NbEntry->nodeptr[j]);
4016 /* assert(0);*/
4017 }
4018 }
4019 #endif
4020
4021 IFDEBUG(gm,4)
4022 INT i;
4023 if (!ioflag)
4024 {
4025 UserWriteF("After qsort\n");
4026
4027 /* test whether all entries are corresponding */
4028 UserWriteF("SORTTABLELIST:\n");
4029 for (i=0; i<Sons_of_Side; i++)
4030 {
4031 COMPARE_RECORD *Entry, *NbEntry;
4032
4033 Entry = ElemSortTable[i];
4034 NbEntry = NbSortTable[i];
4035
4036 UserWriteF("EAdr=%x side=%d realNbAdr=%x NbAdr=%x nbside=%x "
4037 "realNbAdr=%x\n",
4038 Entry->elem, Entry->side, NBELEM(Entry->elem,Entry->side),
4039 NbEntry->elem, NbEntry->side,NBELEM(NbEntry->elem,NbEntry->side));
4040 }
4041
4042 for (i=0; i<Sons_of_Side; i++)
4043 {
4044 COMPARE_RECORD *Entry, *NbEntry;
4045
4046 Entry = ElemSortTable[i];
4047 NbEntry = NbSortTable[i];
4048
4049 if (NBELEM(Entry->elem,Entry->side)!=NbEntry->elem)
4050 {
4051 UserWriteF("NOTEQUAL for i=%d elem=%x: elemrealnb=%x "
4052 "elemsortnb=%x\n",
4053 i,Entry->elem,NBELEM(Entry->elem,Entry->side),NbEntry->elem);
4054 REFINE_ELEMENT_LIST(0,theElement,"theElement:");
4055 REFINE_ELEMENT_LIST(0,theNeighbor,"theNeighbor:");
4056 }
4057 if (NBELEM(NbEntry->elem,NbEntry->side)!=Entry->elem)
4058 {
4059 UserWriteF("NOTEQUAL for i=%d nb=%x: nbrealnb=%x nbsortnb=%x\n",
4060 i,NbEntry->elem,NBELEM(NbEntry->elem,NbEntry->side),
4061 Entry->elem);
4062 REFINE_ELEMENT_LIST(0,theElement,"theE:");
4063 REFINE_ELEMENT_LIST(0,theNeighbor,"theN:");
4064 }
4065 }
4066 UserWriteF("\n\n");
4067 }
4068 ENDDEBUG
4069
4070 /* set neighborship relations */
4071 if (ioflag)
4072 {
4073 COMPARE_RECORD *Entry, *NbEntry;
4074
4075 for (i=0; i<Sons_of_Side; i++)
4076 {
4077 Entry = ElemSortTable[i];
4078 for (k=0; k<Sons_of_NbSide; k++)
4079 {
4080 NbEntry = NbSortTable[k];
4081
4082 if (Entry->nodes != NbEntry->nodes) continue;
4083 for (j=0; j<Entry->nodes; j++)
4084 if (Entry->nodeptr[j] != NbEntry->nodeptr[j])
4085 break;
4086 if (j == Entry->nodes)
4087 {
4088 SET_NBELEM(ElemSortTable[i]->elem,ElemSortTable[i]->side,
4089 NbSortTable[k]->elem);
4090 SET_NBELEM(NbSortTable[k]->elem,NbSortTable[k]->side,
4091 ElemSortTable[i]->elem);
4092 }
4093 }
4094 }
4095 }
4096 else
4097 /* all entires need to match exactly */
4098 for (i=0; i<Sons_of_Side; i++)
4099 {
4100 SET_NBELEM(ElemSortTable[i]->elem,ElemSortTable[i]->side,
4101 NbSortTable[i]->elem);
4102 SET_NBELEM(NbSortTable[i]->elem,NbSortTable[i]->side,
4103 ElemSortTable[i]->elem);
4104 #ifdef __THREEDIM__
4105 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,SIDEVEC))
4106 if (DisposeDoubledSideVector(theGrid,ElemSortTable[i]->elem,
4107 ElemSortTable[i]->side,
4108 NbSortTable[i]->elem,
4109 NbSortTable[i]->side))
4110 RETURN(GM_FATAL);
4111 #endif
4112 }
4113
4114 return(GM_OK);
4115 }
4116
4117 /****************************************************************************/
4118 /** \brief Copy an element
4119
4120 \param theGrid - grid level of sons of theElement
4121 \param theElement - element to refine
4122 \param theContext - nodes needed for new elements
4123
4124 This function copies an element, (i) corner nodes are already allocated, (iv) create son and set references to sons
4125
4126 \return <ul>
4127 <li> 0 - ok </li>
4128 <li> 1 - fatal memory error </li>
4129 </ul>
4130 */
4131 /****************************************************************************/
RefineElementYellow(GRID * theGrid,ELEMENT * theElement,NODE ** theContext)4132 static INT RefineElementYellow (GRID *theGrid, ELEMENT *theElement, NODE **theContext)
4133 {
4134 INT i;
4135 bool boundaryelement = false;
4136 const int me = theGrid->ppifContext().me();
4137
4138 /* check for boundary */
4139 if (OBJT(theElement) == BEOBJ)
4140 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
4141 {
4142 /* at the boundary */
4143 if (SIDE_ON_BND(theElement,i))
4144 {
4145 boundaryelement = true;
4146 break;
4147 }
4148 }
4149
4150 #ifdef Debug
4151 /* check son nodes validity */
4152 for (i=0; i<CORNERS_OF_ELEM(theElement); i++)
4153 {
4154 assert (theContext[i] != NULL);
4155 }
4156 #endif
4157
4158 /* create son */
4159 const INT theSonType = boundaryelement ? BEOBJ : IEOBJ;
4160 ELEMENT* theSon = CreateElement(theGrid,TAG(theElement), theSonType,
4161 theContext,theElement,1);
4162 if (theSon==NULL) RETURN(GM_ERROR);
4163 SETECLASS(theSon,MARKCLASS(theElement));
4164
4165 /* connect son */
4166 IFDEBUG(gm,2)
4167 UserWriteF(PFMT "CONNECTING elem=" EID_FMTX "\n",me,EID_PRTX(theElement));
4168 ENDDEBUG
4169 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
4170 {
4171 INT j,Sons_of_Side;
4172 ELEMENT *Sons_of_Side_List[MAX_SONS];
4173 INT SonSides[MAX_SIDE_NODES];
4174
4175 IFDEBUG(gm,2)
4176 UserWriteF(PFMT " CONNECT side=%i of elem=" EID_FMTX "\n",me,i,EID_PRTX(theElement));
4177 ENDDEBUG
4178
4179 for (j=0; j<MAX_SONS; j++)
4180 Sons_of_Side_List[j] = NULL;
4181
4182 Sons_of_Side = 1;
4183 Sons_of_Side_List[0] = theSon;
4184 SonSides[0] = i;
4185
4186 if (Connect_Sons_of_ElementSide(theGrid,theElement,i,Sons_of_Side,
4187 Sons_of_Side_List,SonSides,0)!=GM_OK) RETURN(GM_FATAL);
4188
4189 #ifdef ModelP
4190 if (Identify_Objects_of_ElementSide(theGrid,theElement,i)) RETURN(GM_FATAL);
4191 #endif
4192 }
4193
4194 return(GM_OK);
4195 }
4196
4197
4198 /****************************************************************************/
4199 /** \brief Refine an element without context
4200
4201 \param theGrid - grid level of sons of theElement
4202 \param theElement - element to refine
4203 \param theContext - nodes needed for new elements
4204
4205 This function refines an element without context,
4206 (i) corner and midnodes are already allocated,
4207 (ii) edges between corner and midnodes are ok,
4208 (iii) create interior nodes and edges,
4209 (iv) create sons and set references to sons.
4210
4211 \return <ul>
4212 INT
4213 .n 0 - ok
4214 .n 1 - fatal memory error
4215 */
4216 /****************************************************************************/
RefineElementGreen(GRID * theGrid,ELEMENT * theElement,NODE ** theContext)4217 static int RefineElementGreen (GRID *theGrid, ELEMENT *theElement, NODE **theContext)
4218 {
4219 struct greensondata
4220 {
4221 /** \brief Element type */
4222 short tag;
4223 /** \brief Boundary element: yes (=1) or no (=0) */
4224 short bdy;
4225 NODE *corners[MAX_CORNERS_OF_ELEM];
4226 int nb[MAX_SIDES_OF_ELEM];
4227 ELEMENT *theSon;
4228 };
4229 typedef struct greensondata GREENSONDATA;
4230
4231 GREENSONDATA sons[MAX_GREEN_SONS];
4232
4233 NODE *theNode, *theNode1;
4234 NODE *theSideNodes[8];
4235 NODE *ElementNodes[MAX_CORNERS_OF_ELEM];
4236 int i,j,k,l,m,n,s;
4237 int nelem,nedges,node0;
4238 bool bdy,found;
4239 int edge, sides[4], side0, side1;
4240 int tetNode0, tetNode1, tetNode2, tetEdge0, tetEdge1, tetEdge2,
4241 tetSideNode0Node1, tetSideNode0Node2, tetSideNode1Node2,
4242 pyrNode0, pyrNode1, pyrNode2, pyrNode3,
4243 pyrEdge0, pyrEdge1, pyrEdge2, pyrEdge3,
4244 pyrSideNode0Node1, pyrSideNode1Node2, pyrSideNode2Node3,
4245 pyrSideNode0Node3;
4246 int elementsSide0[5], elementsSide1[5];
4247
4248 IFDEBUG(gm,1)
4249 UserWriteF("RefineElementGreen(): ELEMENT ID=%d\n",ID(theElement));
4250 ENDDEBUG
4251
4252 /* init son data array */
4253 for (i=0; i<MAX_GREEN_SONS; i++)
4254 {
4255 sons[i].tag = -1;
4256 sons[i].bdy = -1;
4257 for (j=0; j<MAX_CORNERS_OF_ELEM; j++) sons[i].corners[j] = NULL;
4258 for (j=0; j<MAX_SIDES_OF_ELEM; j++) sons[i].nb[j] = -1;
4259 sons[i].theSon = NULL;
4260 }
4261
4262 IFDEBUG(gm,2)
4263 UserWriteF(" Element ID=%d actual CONTEXT is:\n",ID(theElement));
4264 for (i=0; i<MAX_CORNERS_OF_ELEM+MAX_NEW_CORNERS_DIM; i++)
4265 UserWriteF(" %3d",i);
4266 UserWrite("\n");
4267 for (i=0; i<MAX_CORNERS_OF_ELEM+MAX_NEW_CORNERS_DIM; i++)
4268 if (theContext[i] != NULL)
4269 UserWriteF(" %3d",ID(theContext[i]));
4270 else
4271 UserWriteF(" ");
4272 UserWrite("\n");
4273 ENDDEBUG
4274 IFDEBUG(gm,3)
4275 for (i=0; i<MAX_CORNERS_OF_ELEM+MAX_NEW_CORNERS_DIM; i++)
4276 {
4277 if (theContext[i] == NULL) continue;
4278 if (NDOBJ != OBJT(theContext[i]))
4279 UserWriteF(" ERROR NO NDOBJ(5) OBJT(i=%d)=%d ID=%d adr=%x\n",\
4280 i,OBJT(theContext[i]),ID(theContext[i]),theContext[i]);
4281 }
4282 ENDDEBUG
4283
4284 /* init indices for son elements */
4285 /* outer side for tetrahedra is side 0 */
4286 tetNode0 = CORNER_OF_SIDE_TAG(TETRAHEDRON,0,0);
4287 tetNode1 = CORNER_OF_SIDE_TAG(TETRAHEDRON,0,1);
4288 tetNode2 = CORNER_OF_SIDE_TAG(TETRAHEDRON,0,2);
4289
4290 tetEdge0 = EDGE_OF_SIDE_TAG(TETRAHEDRON,0,0);
4291 tetEdge1 = EDGE_OF_SIDE_TAG(TETRAHEDRON,0,1);
4292 tetEdge2 = EDGE_OF_SIDE_TAG(TETRAHEDRON,0,2);
4293
4294 tetSideNode0Node1 = SIDE_WITH_EDGE_TAG(TETRAHEDRON,tetEdge0,0);
4295 if (tetSideNode0Node1 == 0)
4296 tetSideNode0Node1 = SIDE_WITH_EDGE_TAG(TETRAHEDRON,tetEdge0,1);
4297
4298 tetSideNode1Node2 = SIDE_WITH_EDGE_TAG(TETRAHEDRON,tetEdge1,0);
4299 if (tetSideNode1Node2 == 0)
4300 tetSideNode1Node2 = SIDE_WITH_EDGE_TAG(TETRAHEDRON,tetEdge1,1);
4301
4302 tetSideNode0Node2 = SIDE_WITH_EDGE_TAG(TETRAHEDRON,tetEdge2,0);
4303 if (tetSideNode0Node2 == 0)
4304 tetSideNode0Node2 = SIDE_WITH_EDGE_TAG(TETRAHEDRON,tetEdge2,1);
4305
4306 /* outer side for pyramid has 4 corners */
4307 for (i=0; i<SIDES_OF_TAG(PYRAMID); i++)
4308 if (CORNERS_OF_SIDE_TAG(PYRAMID,i) == 4)
4309 break;
4310 pyrNode0 = CORNER_OF_SIDE_TAG(PYRAMID,i,0);
4311 pyrNode1 = CORNER_OF_SIDE_TAG(PYRAMID,i,1);
4312 pyrNode2 = CORNER_OF_SIDE_TAG(PYRAMID,i,2);
4313 pyrNode3 = CORNER_OF_SIDE_TAG(PYRAMID,i,3);
4314
4315 pyrEdge0 = EDGE_OF_SIDE_TAG(PYRAMID,i,0);
4316 pyrEdge1 = EDGE_OF_SIDE_TAG(PYRAMID,i,1);
4317 pyrEdge2 = EDGE_OF_SIDE_TAG(PYRAMID,i,2);
4318 pyrEdge3 = EDGE_OF_SIDE_TAG(PYRAMID,i,3);
4319
4320 pyrSideNode0Node1 = SIDE_WITH_EDGE_TAG(PYRAMID,pyrEdge0,1);
4321 if (pyrSideNode0Node1 == i)
4322 pyrSideNode0Node1 = SIDE_WITH_EDGE_TAG(PYRAMID,pyrEdge0,0);
4323
4324 pyrSideNode1Node2 = SIDE_WITH_EDGE_TAG(PYRAMID,pyrEdge1,1);
4325 if (pyrSideNode1Node2 == i)
4326 pyrSideNode1Node2 = SIDE_WITH_EDGE_TAG(PYRAMID,pyrEdge1,0);
4327
4328 pyrSideNode2Node3 = SIDE_WITH_EDGE_TAG(PYRAMID,pyrEdge2,1);
4329 if (pyrSideNode2Node3 == i)
4330 pyrSideNode2Node3 = SIDE_WITH_EDGE_TAG(PYRAMID,pyrEdge2,0);
4331
4332 pyrSideNode0Node3 = SIDE_WITH_EDGE_TAG(PYRAMID,pyrEdge3,1);
4333 if (pyrSideNode0Node3 == i)
4334 pyrSideNode0Node3 = SIDE_WITH_EDGE_TAG(PYRAMID,pyrEdge3,0);
4335
4336 /* create edges on inner of sides, create son elements and connect them */
4337 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
4338 {
4339 theNode = theContext[CORNERS_OF_ELEM(theElement)+
4340 EDGES_OF_ELEM(theElement)+i];
4341 nedges = EDGES_OF_SIDE(theElement,i);
4342
4343 bdy = (OBJT(theElement) == BEOBJ && SIDE_ON_BND(theElement,i));
4344
4345 nelem = 5*i; /* A face in 3d gets subdivided into at most 5 (yes, 5) parts */
4346 for (j=nelem; j<(nelem+5); j++)
4347 sons[j].bdy = bdy;
4348
4349 k = 0;
4350 for (j=0; j<EDGES_OF_SIDE(theElement,i); j++)
4351 {
4352 edge = EDGE_OF_SIDE(theElement,i,j);
4353 for (l=0; l<MAX_SIDES_OF_ELEM; l++)
4354 if (SIDE_WITH_EDGE(theElement,edge,l) != i)
4355 {
4356 sides[k++] = SIDE_WITH_EDGE(theElement,edge,l)+
4357 MAX_GREEN_SONS;
4358 break;
4359 }
4360 ASSERT(l<2);
4361 }
4362
4363 k = 0;
4364 for (j=0; j<nedges; j++)
4365 {
4366 theSideNodes[2*j] = theContext[CORNER_OF_SIDE(theElement,i,j)];
4367 theSideNodes[2*j+1] = theContext[CORNERS_OF_ELEM(theElement)+
4368 EDGE_OF_SIDE(theElement,i,j)];
4369 if (theSideNodes[2*j+1] != NULL) k++;
4370 }
4371
4372 IFDEBUG(gm,2)
4373 UserWriteF(" SIDE %d has %d nodes and sidenode=%x\n",i,k,theNode);
4374 ENDDEBUG
4375
4376 switch (CORNERS_OF_SIDE(theElement,i))
4377 {
4378
4379 case 4 :
4380
4381 /* theNode points to a potential new side node */
4382 if (theNode == NULL)
4383 {
4384 switch (k) /* The number of nodes on the edges of this side */
4385 {
4386 case 0 :
4387 sons[nelem].tag = PYRAMID;
4388 sons[nelem].corners[pyrNode0] = theSideNodes[0];
4389 sons[nelem].corners[pyrNode1] = theSideNodes[2];
4390 sons[nelem].corners[pyrNode2] = theSideNodes[4];
4391 sons[nelem].corners[pyrNode3] = theSideNodes[6];
4392
4393 sons[nelem].nb[pyrSideNode0Node1] = sides[0];
4394 sons[nelem].nb[pyrSideNode1Node2] = sides[1];
4395 sons[nelem].nb[pyrSideNode2Node3] = sides[2];
4396 sons[nelem].nb[pyrSideNode0Node3] = sides[3];
4397 nelem++;
4398
4399 break;
4400
4401 case 1 :
4402 for (j=0; j<nedges; j++)
4403 {
4404 node0 = 2*j+1;
4405 if (theSideNodes[node0] != NULL)
4406 {
4407
4408 /* define the son corners and inner side relations */
4409 sons[nelem].tag = TETRAHEDRON;
4410 sons[nelem].corners[tetNode0] = theSideNodes[node0];
4411 sons[nelem].corners[tetNode1] = theSideNodes[(node0+1)%(2*nedges)];
4412 sons[nelem].corners[tetNode2] = theSideNodes[(node0+3)%(2*nedges)];
4413
4414 sons[nelem].nb[tetSideNode0Node1] = sides[j];
4415 sons[nelem].nb[tetSideNode1Node2] = sides[(j+1)%nedges];
4416 sons[nelem].nb[tetSideNode0Node2] = nelem+2;
4417 nelem++;
4418
4419 sons[nelem].tag = TETRAHEDRON;
4420 sons[nelem].corners[tetNode0] = theSideNodes[node0];
4421 sons[nelem].corners[tetNode1] = theSideNodes[(node0+5)%(2*nedges)];
4422 sons[nelem].corners[tetNode2] = theSideNodes[(node0+7)%(2*nedges)];
4423
4424 sons[nelem].nb[tetSideNode0Node1] = nelem+1;
4425 sons[nelem].nb[tetSideNode1Node2] = sides[(j+3)%nedges];
4426 sons[nelem].nb[tetSideNode0Node2] = sides[j];
4427 nelem++;
4428
4429 sons[nelem].tag = TETRAHEDRON;
4430 sons[nelem].corners[tetNode0] = theSideNodes[node0];
4431 sons[nelem].corners[tetNode1] = theSideNodes[(node0+3)%(2*nedges)];
4432 sons[nelem].corners[tetNode2] = theSideNodes[(node0+5)%(2*nedges)];
4433
4434 sons[nelem].nb[tetSideNode0Node1] = nelem-2;
4435 sons[nelem].nb[tetSideNode1Node2] = sides[(j+2)%nedges];
4436 sons[nelem].nb[tetSideNode0Node2] = nelem-1;
4437 nelem++;
4438
4439 break;
4440 }
4441 }
4442 break;
4443
4444 case 2 :
4445 /* two cases: sidenodes are not on neighboring edges OR are on neighboring edges */
4446 for (j=0; j<nedges; j++)
4447 {
4448 node0 = 2*j+1;
4449 if (theSideNodes[node0] != NULL)
4450 break;
4451 }
4452 if (theSideNodes[(node0+6)%(2*nedges)] != NULL)
4453 {
4454 node0 = (node0+6)%(2*nedges);
4455 j = (j+3)%nedges;
4456 }
4457 if (theSideNodes[(node0+4)%(2*nedges)] == NULL)
4458 {
4459
4460 sons[nelem].tag = TETRAHEDRON;
4461 sons[nelem].corners[tetNode0] = theSideNodes[node0];
4462 sons[nelem].corners[tetNode1] = theSideNodes[(node0+1)%(2*nedges)];
4463 sons[nelem].corners[tetNode2] = theSideNodes[(node0+2)%(2*nedges)];
4464
4465 sons[nelem].nb[tetSideNode0Node1] = sides[(j)%nedges];
4466 sons[nelem].nb[tetSideNode1Node2] = sides[(j+1)%nedges];
4467 sons[nelem].nb[tetSideNode0Node2] = nelem+3;
4468 nelem++;
4469
4470 sons[nelem].tag = TETRAHEDRON;
4471 sons[nelem].corners[tetNode0] = theSideNodes[node0];
4472 sons[nelem].corners[tetNode1] = theSideNodes[(node0+5)%(2*nedges)];
4473 sons[nelem].corners[tetNode2] = theSideNodes[(node0+7)%(2*nedges)];
4474
4475 sons[nelem].nb[tetSideNode0Node1] = nelem+2;
4476 sons[nelem].nb[tetSideNode1Node2] = sides[(j+3)%nedges];
4477 sons[nelem].nb[tetSideNode0Node2] = sides[(j)%nedges];
4478 nelem++;
4479
4480 sons[nelem].tag = TETRAHEDRON;
4481 sons[nelem].corners[tetNode0] = theSideNodes[(node0+2)%(2*nedges)];
4482 sons[nelem].corners[tetNode1] = theSideNodes[(node0+3)%(2*nedges)];
4483 sons[nelem].corners[tetNode2] = theSideNodes[(node0+5)%(2*nedges)];
4484
4485 sons[nelem].nb[tetSideNode0Node1] = sides[(j+1)%nedges];
4486 sons[nelem].nb[tetSideNode1Node2] = sides[(j+2)%nedges];
4487 sons[nelem].nb[tetSideNode0Node2] = nelem+1;
4488 nelem++;
4489
4490 sons[nelem].tag = TETRAHEDRON;
4491 sons[nelem].corners[tetNode0] = theSideNodes[node0];
4492 sons[nelem].corners[tetNode1] = theSideNodes[(node0+2)%(2*nedges)];
4493 sons[nelem].corners[tetNode2] = theSideNodes[(node0+5)%(2*nedges)];
4494
4495 sons[nelem].nb[tetSideNode0Node1] = nelem-3;
4496 sons[nelem].nb[tetSideNode1Node2] = nelem-1;
4497 sons[nelem].nb[tetSideNode0Node2] = nelem-2;
4498 nelem++;
4499 }
4500 else
4501 {
4502 sons[nelem].tag = PYRAMID;
4503 sons[nelem].corners[pyrNode0] = theSideNodes[node0];
4504 sons[nelem].corners[pyrNode1] = theSideNodes[(node0+1)%(2*nedges)];
4505 sons[nelem].corners[pyrNode2] = theSideNodes[(node0+3)%(2*nedges)];
4506 sons[nelem].corners[pyrNode3] = theSideNodes[(node0+4)%(2*nedges)];
4507
4508 sons[nelem].nb[pyrSideNode0Node1] = sides[(j)%nedges];
4509 sons[nelem].nb[pyrSideNode1Node2] = sides[(j+1)%nedges];
4510 sons[nelem].nb[pyrSideNode2Node3] = sides[(j+2)%nedges];
4511 sons[nelem].nb[pyrSideNode0Node3] = nelem+1;
4512 nelem++;
4513
4514 sons[nelem].tag = PYRAMID;
4515 sons[nelem].corners[pyrNode0] = theSideNodes[(node0+4)%(2*nedges)];
4516 sons[nelem].corners[pyrNode1] = theSideNodes[(node0+5)%(2*nedges)];
4517 sons[nelem].corners[pyrNode2] = theSideNodes[(node0+7)%(2*nedges)];
4518 sons[nelem].corners[pyrNode3] = theSideNodes[(node0+8)%(2*nedges)];
4519
4520 sons[nelem].nb[pyrSideNode0Node1] = sides[(j+2)%nedges];
4521 sons[nelem].nb[pyrSideNode1Node2] = sides[(j+3)%nedges];
4522 sons[nelem].nb[pyrSideNode2Node3] = sides[(j)%nedges];
4523 sons[nelem].nb[pyrSideNode0Node3] = nelem-1;
4524 nelem++;
4525 }
4526 break;
4527
4528 case 3 :
4529 for (j=0; j<nedges; j++)
4530 {
4531 node0 = 2*j+1;
4532 if (theSideNodes[node0] == NULL)
4533 break;
4534 }
4535
4536 sons[nelem].tag = PYRAMID;
4537 sons[nelem].corners[pyrNode0] = theSideNodes[(node0+1)%(2*nedges)];
4538 sons[nelem].corners[pyrNode1] = theSideNodes[(node0+2)%(2*nedges)];
4539 sons[nelem].corners[pyrNode2] = theSideNodes[(node0+6)%(2*nedges)];
4540 sons[nelem].corners[pyrNode3] = theSideNodes[(node0+7)%(2*nedges)];
4541
4542 sons[nelem].nb[pyrSideNode0Node1] = sides[(j+1)%nedges];
4543 sons[nelem].nb[pyrSideNode1Node2] = nelem+3;
4544 sons[nelem].nb[pyrSideNode2Node3] = sides[(j+3)%nedges];
4545 sons[nelem].nb[pyrSideNode0Node3] = sides[(j)%nedges];
4546 nelem++;
4547
4548 sons[nelem].tag = TETRAHEDRON;
4549 sons[nelem].corners[tetNode0] = theSideNodes[(node0+2)%(2*nedges)];
4550 sons[nelem].corners[tetNode1] = theSideNodes[(node0+3)%(2*nedges)];
4551 sons[nelem].corners[tetNode2] = theSideNodes[(node0+4)%(2*nedges)];
4552
4553 sons[nelem].nb[tetSideNode0Node1] = sides[(j+1)%nedges];
4554 sons[nelem].nb[tetSideNode1Node2] = sides[(j+2)%nedges];
4555 sons[nelem].nb[tetSideNode0Node2] = nelem+2;
4556 nelem++;
4557
4558 sons[nelem].tag = TETRAHEDRON;
4559 sons[nelem].corners[tetNode0] = theSideNodes[(node0+4)%(2*nedges)];
4560 sons[nelem].corners[tetNode1] = theSideNodes[(node0+5)%(2*nedges)];
4561 sons[nelem].corners[tetNode2] = theSideNodes[(node0+6)%(2*nedges)];
4562
4563 sons[nelem].nb[tetSideNode0Node1] = sides[(j+2)%nedges];
4564 sons[nelem].nb[tetSideNode1Node2] = sides[(j+3)%nedges];
4565 sons[nelem].nb[tetSideNode0Node2] = nelem+1;
4566 nelem++;
4567
4568 sons[nelem].tag = TETRAHEDRON;
4569 sons[nelem].corners[tetNode0] = theSideNodes[(node0+2)%(2*nedges)];
4570 sons[nelem].corners[tetNode1] = theSideNodes[(node0+4)%(2*nedges)];
4571 sons[nelem].corners[tetNode2] = theSideNodes[(node0+6)%(2*nedges)];
4572
4573 sons[nelem].nb[tetSideNode0Node1] = nelem-2;
4574 sons[nelem].nb[tetSideNode1Node2] = nelem-1;
4575 sons[nelem].nb[tetSideNode0Node2] = nelem-3;
4576 nelem++;
4577
4578 break;
4579
4580 case 4 :
4581 for (j=0; j<nedges; j++)
4582 {
4583 node0 = 2*j+1;
4584
4585 sons[nelem].tag = TETRAHEDRON;
4586 sons[nelem].corners[tetNode0] = theSideNodes[node0];
4587 sons[nelem].corners[tetNode1] = theSideNodes[(node0+1)%(2*nedges)];
4588 sons[nelem].corners[tetNode2] = theSideNodes[(node0+2)%(2*nedges)];
4589
4590 sons[nelem].nb[tetSideNode0Node1] = sides[(j)%nedges];
4591 sons[nelem].nb[tetSideNode1Node2] = sides[(j+1)%nedges];
4592 sons[nelem].nb[tetSideNode0Node2] = nelem+(nedges-j);
4593 nelem++;
4594 }
4595
4596 sons[nelem].tag = PYRAMID;
4597 sons[nelem].corners[pyrNode0] = theSideNodes[1];
4598 sons[nelem].corners[pyrNode1] = theSideNodes[3];
4599 sons[nelem].corners[pyrNode2] = theSideNodes[5];
4600 sons[nelem].corners[pyrNode3] = theSideNodes[7];
4601
4602 sons[nelem].nb[pyrSideNode0Node1] = nelem-4;
4603 sons[nelem].nb[pyrSideNode1Node2] = nelem-3;
4604 sons[nelem].nb[pyrSideNode2Node3] = nelem-2;
4605 sons[nelem].nb[pyrSideNode0Node3] = nelem-1;
4606 nelem++;
4607 break;
4608
4609 default :
4610 RETURN(GM_FATAL);
4611 }
4612 }
4613 else /* theNode != NULL */
4614 {
4615 /* create the four side edges */
4616 for (j=0; j<nedges; j++)
4617 {
4618 node0 = 2*j+1;
4619 if (theSideNodes[node0] == NULL) break;
4620
4621 sons[nelem].tag = PYRAMID;
4622 sons[nelem].corners[pyrNode0] = theSideNodes[node0%(2*nedges)];
4623 sons[nelem].corners[pyrNode1] = theSideNodes[(node0+1)%(2*nedges)];
4624 sons[nelem].corners[pyrNode2] = theSideNodes[(node0+2)%(2*nedges)];
4625 sons[nelem].corners[pyrNode3] = theNode;
4626
4627 sons[nelem].nb[pyrSideNode0Node1] = sides[(j)%nedges];
4628 sons[nelem].nb[pyrSideNode1Node2] = sides[(j+1)%nedges];
4629 if (j == 3)
4630 sons[nelem].nb[pyrSideNode2Node3] = nelem-3;
4631 else
4632 sons[nelem].nb[pyrSideNode2Node3] = nelem+1;
4633 if (j == 0)
4634 sons[nelem].nb[pyrSideNode0Node3] = nelem+3;
4635 else
4636 sons[nelem].nb[pyrSideNode0Node3] = nelem-1;
4637 nelem++;
4638
4639 }
4640 ASSERT(j==4);
4641 }
4642
4643 break;
4644
4645 case 3 :
4646 if (theNode == NULL)
4647 {
4648 switch (k)
4649 {
4650 case 0 :
4651 sons[nelem].tag = TETRAHEDRON;
4652 sons[nelem].corners[tetNode0] = theSideNodes[0];
4653 sons[nelem].corners[tetNode1] = theSideNodes[2];
4654 sons[nelem].corners[tetNode2] = theSideNodes[4];
4655
4656 sons[nelem].nb[tetSideNode0Node1] = sides[0];
4657 sons[nelem].nb[tetSideNode1Node2] = sides[1];
4658 sons[nelem].nb[tetSideNode0Node2] = sides[2];
4659 nelem++;
4660
4661 break;
4662
4663 case 1 :
4664 for (j=0; j<nedges; j++)
4665 {
4666 node0 = 2*j+1;
4667 if (theSideNodes[node0] != NULL)
4668 {
4669
4670 /* define the son corners and inner side relations */
4671 sons[nelem].tag = TETRAHEDRON;
4672 sons[nelem].corners[tetNode0] = theSideNodes[node0];
4673 sons[nelem].corners[tetNode1] = theSideNodes[(node0+1)%(2*nedges)];
4674 sons[nelem].corners[tetNode2] = theSideNodes[(node0+3)%(2*nedges)];
4675
4676 sons[nelem].nb[tetSideNode0Node1] = sides[j];
4677 sons[nelem].nb[tetSideNode1Node2] = sides[(j+1)%nedges];
4678 sons[nelem].nb[tetSideNode0Node2] = nelem+1;
4679 nelem++;
4680
4681 sons[nelem].tag = TETRAHEDRON;
4682 sons[nelem].corners[tetNode0] = theSideNodes[node0];
4683 sons[nelem].corners[tetNode1] = theSideNodes[(node0+3)%(2*nedges)];
4684 sons[nelem].corners[tetNode2] = theSideNodes[(node0+5)%(2*nedges)];
4685
4686 sons[nelem].nb[tetSideNode0Node1] = nelem-1;
4687 sons[nelem].nb[tetSideNode1Node2] = sides[(j+2)%nedges];
4688 sons[nelem].nb[tetSideNode0Node2] = sides[j];
4689 nelem++;
4690
4691 break;
4692 }
4693 }
4694 break;
4695
4696 case 2 :
4697 {
4698 INT node,k;
4699 /*INT maxedge=-1;*/
4700 #ifdef ModelP
4701 unsigned int maxid = 0;
4702 #else
4703 INT maxid = -1;
4704 #endif
4705
4706 node0 = -1;
4707 for (k=0; k<nedges; k++)
4708 {
4709 node = (2*k+3)%(2*nedges);
4710 if (theSideNodes[node] == NULL)
4711 {
4712 node0 = 2*k+1;
4713 j = k;
4714 }
4715 /*
4716 if (EDGE_OF_SIDE(theElement,i,k)>maxedge)
4717 maxedge = 2*k+1;
4718 */
4719 /* neighboring elements need to refine in the same way */
4720 /* in parallel case all copies of the elements also */
4721 if (theSideNodes[2*k+1]!=NULL && _ID_(theSideNodes[2*k+1])>maxid)
4722 maxid = _ID_(theSideNodes[2*k+1]);
4723 }
4724 assert(maxid != -1);
4725 assert(node0 != -1);
4726
4727 /* if (node0 == maxedge && ((SIDEPATTERN(theElement)&(1<<i)) == 0)) */
4728 if (_ID_(theSideNodes[node0]) == maxid)
4729 {
4730
4731 sons[nelem].tag = TETRAHEDRON;
4732 sons[nelem].corners[tetNode0] = theSideNodes[node0];
4733 sons[nelem].corners[tetNode1] = theSideNodes[(node0+1)%(2*nedges)];
4734 sons[nelem].corners[tetNode2] = theSideNodes[(node0+3)%(2*nedges)];
4735
4736 sons[nelem].nb[tetSideNode0Node1] = sides[j];
4737 sons[nelem].nb[tetSideNode1Node2] = sides[(j+1)%nedges];
4738 sons[nelem].nb[tetSideNode0Node2] = nelem+2;
4739 nelem++;
4740
4741 sons[nelem].tag = TETRAHEDRON;
4742 sons[nelem].corners[tetNode0] = theSideNodes[node0];
4743 sons[nelem].corners[tetNode1] = theSideNodes[(node0+4)%(2*nedges)];
4744 sons[nelem].corners[tetNode2] = theSideNodes[(node0+5)%(2*nedges)];
4745
4746 sons[nelem].nb[tetSideNode0Node1] = nelem+1;
4747 sons[nelem].nb[tetSideNode1Node2] = sides[(j+2)%nedges];
4748 sons[nelem].nb[tetSideNode0Node2] = sides[j];
4749 nelem++;
4750
4751 sons[nelem].tag = TETRAHEDRON;
4752 sons[nelem].corners[tetNode0] = theSideNodes[node0];
4753 sons[nelem].corners[tetNode1] = theSideNodes[(node0+3)%(2*nedges)];
4754 sons[nelem].corners[tetNode2] = theSideNodes[(node0+4)%(2*nedges)];
4755
4756 sons[nelem].nb[tetSideNode0Node1] = nelem-2;
4757 sons[nelem].nb[tetSideNode1Node2] = sides[(j+2)%nedges];
4758 sons[nelem].nb[tetSideNode0Node2] = nelem-1;
4759 nelem++;
4760
4761 }
4762 else
4763 {
4764 sons[nelem].tag = TETRAHEDRON;
4765 sons[nelem].corners[tetNode0] = theSideNodes[node0];
4766 sons[nelem].corners[tetNode1] = theSideNodes[(node0+1)%(2*nedges)];
4767 sons[nelem].corners[tetNode2] = theSideNodes[(node0+4)%(2*nedges)];
4768
4769 sons[nelem].nb[tetSideNode0Node1] = sides[j];
4770 sons[nelem].nb[tetSideNode1Node2] = nelem+1;
4771 sons[nelem].nb[tetSideNode0Node2] = nelem+2;
4772 nelem++;
4773
4774 sons[nelem].tag = TETRAHEDRON;
4775 sons[nelem].corners[tetNode0] = theSideNodes[(node0+4)%(2*nedges)];
4776 sons[nelem].corners[tetNode1] = theSideNodes[(node0+1)%(2*nedges)];
4777 sons[nelem].corners[tetNode2] = theSideNodes[(node0+3)%(2*nedges)];
4778
4779 sons[nelem].nb[tetSideNode0Node1] = nelem-1;
4780 sons[nelem].nb[tetSideNode1Node2] = sides[(j+1)%nedges];
4781 sons[nelem].nb[tetSideNode0Node2] = sides[(j+2)%nedges];
4782 nelem++;
4783
4784 sons[nelem].tag = TETRAHEDRON;
4785 sons[nelem].corners[tetNode0] = theSideNodes[node0];
4786 sons[nelem].corners[tetNode1] = theSideNodes[(node0+4)%(2*nedges)];
4787 sons[nelem].corners[tetNode2] = theSideNodes[(node0+5)%(2*nedges)];
4788
4789 sons[nelem].nb[tetSideNode0Node1] = nelem-2;
4790 sons[nelem].nb[tetSideNode1Node2] = sides[(j+2)%nedges];
4791 sons[nelem].nb[tetSideNode0Node2] = sides[j];
4792 nelem++;
4793
4794 }
4795
4796 break;
4797 }
4798 case 3 :
4799 j = 0;
4800 node0 = 1;
4801
4802 sons[nelem].tag = TETRAHEDRON;
4803 sons[nelem].corners[tetNode0] = theSideNodes[node0];
4804 sons[nelem].corners[tetNode1] = theSideNodes[(node0+1)%(2*nedges)];
4805 sons[nelem].corners[tetNode2] = theSideNodes[(node0+2)%(2*nedges)];
4806
4807 sons[nelem].nb[tetSideNode0Node1] = sides[j];
4808 sons[nelem].nb[tetSideNode1Node2] = sides[(j+1)%nedges];
4809 sons[nelem].nb[tetSideNode0Node2] = nelem+3;
4810 nelem++;
4811
4812 sons[nelem].tag = TETRAHEDRON;
4813 sons[nelem].corners[tetNode0] = theSideNodes[node0];
4814 sons[nelem].corners[tetNode1] = theSideNodes[(node0+4)%(2*nedges)];
4815 sons[nelem].corners[tetNode2] = theSideNodes[(node0+5)%(2*nedges)];
4816
4817 sons[nelem].nb[tetSideNode0Node1] = nelem+2;
4818 sons[nelem].nb[tetSideNode1Node2] = sides[(j+2)%nedges];
4819 sons[nelem].nb[tetSideNode0Node2] = sides[j];
4820 nelem++;
4821
4822 sons[nelem].tag = TETRAHEDRON;
4823 sons[nelem].corners[tetNode0] = theSideNodes[node0+2];
4824 sons[nelem].corners[tetNode1] = theSideNodes[(node0+3)%(2*nedges)];
4825 sons[nelem].corners[tetNode2] = theSideNodes[(node0+4)%(2*nedges)];
4826
4827 sons[nelem].nb[tetSideNode0Node1] = sides[(j+1)%nedges];
4828 sons[nelem].nb[tetSideNode1Node2] = sides[(j+2)%nedges];
4829 sons[nelem].nb[tetSideNode0Node2] = nelem+1;
4830 nelem++;
4831
4832 sons[nelem].tag = TETRAHEDRON;
4833 sons[nelem].corners[tetNode0] = theSideNodes[node0];
4834 sons[nelem].corners[tetNode1] = theSideNodes[(node0+2)%(2*nedges)];
4835 sons[nelem].corners[tetNode2] = theSideNodes[(node0+4)%(2*nedges)];
4836
4837 sons[nelem].nb[tetSideNode0Node1] = nelem-3;
4838 sons[nelem].nb[tetSideNode1Node2] = nelem-1;
4839 sons[nelem].nb[tetSideNode0Node2] = nelem-2;
4840 nelem++;
4841
4842 break;
4843
4844 default :
4845 assert(0);
4846 break;
4847 }
4848 }
4849
4850 else /* theNode != NULL */
4851 {
4852 /* create the four side edges */
4853 for (j=0; j<nedges; j++)
4854 {
4855 node0 = 2*j+1;
4856 if (theSideNodes[node0] == NULL) break;
4857
4858 sons[nelem].tag = PYRAMID;
4859 sons[nelem].corners[pyrNode0] = theSideNodes[node0%(2*nedges)];
4860 sons[nelem].corners[pyrNode1] = theSideNodes[(node0+1)%(2*nedges)];
4861 sons[nelem].corners[pyrNode2] = theSideNodes[(node0+2)%(2*nedges)];
4862 sons[nelem].corners[pyrNode3] = theNode;
4863
4864 sons[nelem].nb[pyrSideNode0Node1] = sides[(j)%nedges];
4865 sons[nelem].nb[pyrSideNode1Node2] = sides[(j+1)%nedges];
4866 if (j == 2)
4867 sons[nelem].nb[pyrSideNode2Node3] = nelem-2;
4868 else
4869 sons[nelem].nb[pyrSideNode2Node3] = nelem+1;
4870 if (j == 0)
4871 sons[nelem].nb[pyrSideNode0Node3] = nelem+2;
4872 else
4873 sons[nelem].nb[pyrSideNode0Node3] = nelem-1;
4874 nelem++;
4875
4876 }
4877 ASSERT(j==nedges);
4878 }
4879
4880 break;
4881
4882 default : /* Side with neither 3 nor 4 vertices found */
4883 assert(0);
4884 break;
4885 }
4886 }
4887
4888 /* connect elements over edges */
4889 for (i=0; i<EDGES_OF_ELEM(theElement); i++)
4890 {
4891 side0 = SIDE_WITH_EDGE(theElement,i,0);
4892 side1 = SIDE_WITH_EDGE(theElement,i,1);
4893
4894 if (theContext[i+CORNERS_OF_ELEM(theElement)] == NULL) /* No new node in the middle of this edge */
4895 {
4896 /* two elements share this edge */
4897
4898 /* get son elements for this edge */
4899 found = false;
4900 for (j=side0*5; j<(side0*5+5); j++)
4901 {
4902 for (k=0; k<MAX_SIDES_OF_ELEM; k++)
4903 if ((sons[j].nb[k]-MAX_GREEN_SONS)==side1)
4904 {
4905 found = true;
4906 break;
4907 }
4908 if (found) break;
4909 }
4910 ASSERT(j<side0*5+5);
4911
4912 found = false;
4913 for (l=side1*5; l<side1*5+5; l++)
4914 {
4915 for (m=0; m<MAX_SIDES_OF_ELEM; m++)
4916 if ((sons[l].nb[m]-MAX_GREEN_SONS)==side0)
4917 {
4918 found = true;
4919 break;
4920 }
4921 if (found) break;
4922 }
4923 ASSERT(l<side1*5+5);
4924
4925 sons[j].nb[k] = l;
4926 sons[l].nb[m] = j;
4927 }
4928 else
4929 {
4930 /* four elements share this edge */
4931
4932 /* get son elements for this edge */
4933 l = 0;
4934 for (j=side0*5; j<(side0*5+5); j++)
4935 {
4936 for (k=0; k<MAX_SIDES_OF_ELEM; k++)
4937 if ((sons[j].nb[k]-MAX_GREEN_SONS)==side1)
4938 elementsSide0[l++] = j;
4939 }
4940 ASSERT(l==2);
4941
4942 l = 0;
4943 for (j=side1*5; j<(side1*5+5); j++)
4944 {
4945 for (m=0; m<MAX_SIDES_OF_ELEM; m++)
4946 if ((sons[j].nb[m]-MAX_GREEN_SONS)==side0)
4947 elementsSide1[l++] = j;
4948 }
4949 ASSERT(l==2);
4950
4951 /* determine neighboring elements */
4952 for (j=0; j<CORNERS_OF_EDGE; j++)
4953 {
4954 theNode1 = theContext[CORNER_OF_EDGE(theElement,i,j)];
4955 found = false;
4956 for (l=0; l<2; l++)
4957 {
4958 for (k=0; k<MAX_CORNERS_OF_ELEM; k++)
4959 {
4960 if (theNode1 == sons[elementsSide0[l]].corners[k])
4961 {
4962 found = true;
4963 break;
4964 }
4965 }
4966 if (found) break;
4967 }
4968 ASSERT(k<MAX_CORNERS_OF_ELEM);
4969 ASSERT(l<2);
4970
4971 found = false;
4972 for (m=0; m<2; m++)
4973 {
4974 for (k=0; k<MAX_CORNERS_OF_ELEM; k++)
4975 {
4976 if (theNode1 == sons[elementsSide1[m]].corners[k])
4977 {
4978 found = true;
4979 break;
4980 }
4981 }
4982 if (found) break;
4983 }
4984 ASSERT(k<MAX_CORNERS_OF_ELEM);
4985 ASSERT(m<2);
4986
4987 /* init neighbor field */
4988 for (k=0; k<MAX_SIDES_OF_ELEM; k++)
4989 if ((sons[elementsSide0[l]].nb[k]-MAX_GREEN_SONS)==side1)
4990 break;
4991 ASSERT(k<MAX_SIDES_OF_ELEM);
4992 sons[elementsSide0[l]].nb[k] = elementsSide1[m];
4993
4994 for (k=0; k<MAX_SIDES_OF_ELEM; k++)
4995 if ((sons[elementsSide1[m]].nb[k]-MAX_GREEN_SONS)==side0)
4996 break;
4997 ASSERT(k<MAX_SIDES_OF_ELEM);
4998 sons[elementsSide1[m]].nb[k] = elementsSide0[l];
4999 }
5000 }
5001 }
5002
5003 /* create son elements */
5004 IFDEBUG(gm,1)
5005 UserWriteF(" Creating SON elements for element ID=%d:\n",ID(theElement));
5006 ENDDEBUG
5007 n = 0;
5008 for (i=0; i<MAX_GREEN_SONS; i++)
5009 {
5010 if (sons[i].tag >= 0)
5011 {
5012
5013 IFDEBUG(gm,2)
5014 if (i%5 == 0)
5015 UserWriteF(" SIDE %d:\n",i/5);
5016 ENDDEBUG
5017
5018 l = 0;
5019 for (j=0; j<CORNERS_OF_TAG(sons[i].tag); j++)
5020 {
5021 if (sons[i].corners[j] == NULL)
5022 {
5023 sons[i].corners[j] = theContext[CORNERS_OF_ELEM(theElement)+
5024 CENTER_NODE_INDEX(theElement)];
5025 l++;
5026 }
5027 ElementNodes[j] = sons[i].corners[j];
5028 }
5029 ASSERT(l == 1);
5030
5031 if (sons[i].bdy == 1)
5032 sons[i].theSon = CreateElement(theGrid,sons[i].tag,BEOBJ,
5033 ElementNodes,theElement,1);
5034 else
5035 sons[i].theSon = CreateElement(theGrid,sons[i].tag,IEOBJ,
5036 ElementNodes,theElement,1);
5037 if (sons[i].theSon==NULL) RETURN(GM_FATAL);
5038
5039 IFDEBUG(gm,0)
5040 for (j=0; j<CORNERS_OF_ELEM(sons[i].theSon); j++)
5041 for (m=0; m<CORNERS_OF_ELEM(sons[i].theSon); m++)
5042 if (sons[i].corners[j] == NULL || sons[i].corners[m] == NULL)
5043 {
5044 if ((m!=j) && (sons[i].corners[j] == sons[i].corners[m]))
5045 UserWriteF(" ERROR: son %d has equivalent corners "
5046 "%d=%d adr=%x adr=%x\n",
5047 n,j,m,sons[i].corners[j],sons[i].corners[m]);
5048 }
5049 else
5050 if ((m!=j) && (sons[i].corners[j] == sons[i].corners[m] ||
5051 (_ID_(sons[i].corners[j]) == _ID_(sons[i].corners[m]))))
5052 UserWriteF(" ERROR: son %d has equivalent corners "
5053 "%d=%d ID=%d ID=%d adr=%x adr=%x\n",
5054 n,j,m,_ID_(sons[i].corners[j]),_ID_(sons[i].corners[m]),
5055 sons[i].corners[j],sons[i].corners[m]);
5056 ENDDEBUG
5057
5058 IFDEBUG(gm,2)
5059 UserWriteF(" SONS[i=%d] ID=%d: CORNERS ",i,ID(sons[i].theSon));
5060 for (j=0; j<CORNERS_OF_ELEM(sons[i].theSon); j++)
5061 {
5062 if (sons[i].corners[j] != NULL)
5063 UserWriteF(" %d",_ID_(sons[i].corners[j]));
5064 }
5065 UserWriteF("\n");
5066 ENDDEBUG
5067
5068 #ifdef __ANISOTROPIC__
5069 if (MARK(theElement) != COPY)
5070 SETECLASS(sons[i].theSon,RED_CLASS);
5071 else
5072 SETECLASS(sons[i].theSon,GREEN_CLASS);
5073 #else
5074 SETECLASS(sons[i].theSon,GREEN_CLASS);
5075 #endif
5076 if (i == 0) SET_SON(theElement,0,sons[i].theSon);
5077 for (s=0; s<SIDES_OF_ELEM(sons[i].theSon); s++)
5078 SET_NBELEM(sons[i].theSon,s,NULL);
5079
5080 n++;
5081 }
5082 }
5083 IFDEBUG(gm,1)
5084 UserWriteF(" n=%d sons created NSONS=%d\n",n,NSONS(theElement));
5085 ENDDEBUG
5086
5087 /* translate neighbor information */
5088 for (i=0; i<MAX_GREEN_SONS; i++)
5089 {
5090 if (sons[i].tag >= 0) /* valid son entry */
5091 {
5092 l = 0;
5093 IFDEBUG(gm,0)
5094 for (j=0; j<SIDES_OF_ELEM(sons[i].theSon); j++)
5095 {
5096 for (m=0; m<SIDES_OF_ELEM(sons[i].theSon); m++)
5097 {
5098 if (sons[i].nb[j] == sons[i].nb[m] && (m!=j))
5099 UserWriteF(" ERROR: son %d has equivalent "
5100 "neighbors %d=%d NB=%d\n",n,j,m,sons[i].nb[m]);
5101 }
5102 }
5103 ENDDEBUG
5104 for (j=0; j<SIDES_OF_ELEM(sons[i].theSon); j++)
5105 {
5106 if (sons[i].nb[j] != -1)
5107 SET_NBELEM(sons[i].theSon,j,sons[sons[i].nb[j]].theSon);
5108 else
5109 l++;
5110 }
5111 /* l counts the number of element sides without a neighboring element.
5112 * Since all elements are pyramids/tetrahedra with exactly one vertex in the interior,
5113 * this value must be one. */
5114 ASSERT(l == 1);
5115 }
5116 }
5117
5118 #ifdef __THREEDIM__
5119 /* If there are side vectors for the elements, then the CreateElement methods called above have
5120 * allocated one side vector for each new element face. Therefore, for each face shared by
5121 * two elements there are now two side vectors, even though there should be only one (shared).
5122 * The following loops gets rid of the second redundant side vector. The redundant side
5123 * vectors shared with elements in the rest of the grid are treated in the method Connect_Sons_of_ElementSide,
5124 * called further below.
5125 */
5126 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,SIDEVEC))
5127 {
5128 for (i=0; i<MAX_GREEN_SONS; i++)
5129 {
5130 if (sons[i].tag < 0) /* empty son entry */
5131 continue;
5132
5133 for (j=0; j<SIDES_OF_ELEM(sons[i].theSon); j++)
5134 {
5135 if (sons[i].nb[j] != -1) /* We have neighbor on the j-th face */
5136 {
5137 if (sons[i].nb[j] <= i) /* Visit every element pair only once */
5138 continue;
5139
5140 ELEMENT* theNeighbor = sons[sons[i].nb[j]].theSon;
5141 /* what neighbor are we for the neighbor? */
5142 for (l=0; l<SIDES_OF_ELEM(theNeighbor); l++)
5143 if (sons[sons[i].nb[j]].nb[l] == i)
5144 break;
5145 ASSERT(l < SIDES_OF_ELEM(theNeighbor));
5146
5147 DisposeDoubledSideVector(theGrid, sons[i].theSon, j, theNeighbor, l);
5148 }
5149 }
5150 }
5151 }
5152 #endif
5153
5154 /* connect sons over outer sides */
5155 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
5156 {
5157 INT j,k,Sons_of_Side;
5158 ELEMENT *Sons_of_Side_List[MAX_SONS];
5159 INT SonSides[MAX_SIDE_NODES];
5160
5161 Sons_of_Side = 0;
5162
5163 for (j=0; j<MAX_SONS; j++)
5164 Sons_of_Side_List[j] = NULL;
5165
5166 for (j=0; j<5; j++)
5167 {
5168 if (sons[i*5+j].tag < 0) break;
5169 Sons_of_Side_List[j] = sons[i*5+j].theSon;
5170 Sons_of_Side++;
5171 SonSides[j] = 0;
5172 if (sons[i*5+j].tag == PYRAMID)
5173 {
5174 for (k=0; k<SIDES_OF_TAG(PYRAMID); k++)
5175 if (CORNERS_OF_SIDE_TAG(PYRAMID,k) == 4)
5176 break;
5177 SonSides[j] = k;
5178 }
5179 }
5180 assert(Sons_of_Side>0 && Sons_of_Side<6);
5181
5182 if (Connect_Sons_of_ElementSide(theGrid,theElement,i,Sons_of_Side,
5183 Sons_of_Side_List,SonSides,0)!=GM_OK) RETURN(GM_FATAL);
5184
5185 #ifdef ModelP
5186 if (Identify_Objects_of_ElementSide(theGrid,theElement,i)) RETURN(GM_FATAL);
5187 #endif
5188 }
5189
5190 return(GM_OK);
5191 }
5192
5193
5194 /****************************************************************************/
5195 /*
5196 RefineElementRed - refine an element in the given context
5197
5198 SYNOPSIS:
5199 static int RefineElementRed (GRID *theGrid, ELEMENT *theElement, NODE **theElementContext);
5200
5201 PARAMETERS:
5202 \param theGrid - grid level of sons of theElement
5203 \param theElement - element to refine
5204 \param theContext - current context of element
5205
5206 DESCRIPTION:
5207 This function refines an element in the given context,
5208 (i) corner and midnodes are already allocated,
5209 (ii) edges between corner and midnodes are ok,
5210 (iii) create interior nodes and edges,
5211 (iv) create sons and set references to sons.
5212
5213 \return <ul>
5214 INT
5215 .n GM_OK - ok
5216 .n GM_FATAL - fatal memory error
5217 */
5218 /****************************************************************************/
RefineElementRed(GRID * theGrid,ELEMENT * theElement,NODE ** theElementContext)5219 static int RefineElementRed (GRID *theGrid, ELEMENT *theElement, NODE **theElementContext)
5220 {
5221 const int me = theGrid->ppifContext().me();
5222
5223 /* is something to do ? */
5224 if (!MARKED(theElement)) return(GM_OK);
5225
5226 ELEMENT *SonList[MAX_SONS];
5227 for (INT i=0; i<MAX_SONS; i++) SonList[i] = NULL;
5228
5229 REFRULE const* rule = MARK2RULEADR(theElement,MARK(theElement));
5230
5231 /* create elements */
5232 for (INT s=0; s<NSONS_OF_RULE(rule); s++)
5233 {
5234 bool boundaryelement = false;
5235 /** \todo how can boundary detection be generalized */
5236 if (OBJT(theElement) == BEOBJ)
5237 {
5238 for (INT i=0; i<SIDES_OF_TAG(SON_TAG_OF_RULE(rule,s)); i++)
5239 {
5240 const INT side = SON_NB_OF_RULE(rule,s,i);
5241 /* exterior side */
5242 if (side >= FATHER_SIDE_OFFSET)
5243 {
5244 /* at the boundary */
5245 if (SIDE_ON_BND(theElement,side-FATHER_SIDE_OFFSET))
5246 {
5247 boundaryelement = true;
5248 break;
5249 }
5250 }
5251 }
5252 }
5253
5254 NODE* ElementNodes[MAX_CORNERS_OF_ELEM];
5255 for (INT i=0; i<CORNERS_OF_TAG(SON_TAG_OF_RULE(rule,s)); i++)
5256 {
5257 ASSERT(theElementContext[SON_CORNER_OF_RULE(rule,s,i)]!=NULL);
5258 ElementNodes[i] = theElementContext[SON_CORNER_OF_RULE(rule,s,i)];
5259 }
5260
5261 const INT theSonType = boundaryelement ? BEOBJ : IEOBJ;
5262 ELEMENT* theSon = CreateElement(theGrid,SON_TAG_OF_RULE(rule,s),theSonType,
5263 ElementNodes,theElement,1);
5264 if (theSon==NULL) RETURN(GM_ERROR);
5265
5266 /* fill in son data */
5267 SonList[s] = theSon;
5268 SETECLASS(theSon,MARKCLASS(theElement));
5269 }
5270
5271 /* connect elements */
5272 for (INT s=0; s<NSONS_OF_RULE(rule); s++)
5273 {
5274 SONDATA const* sdata = SON_OF_RULE(rule,s);
5275 for (INT i=0; i<SIDES_OF_ELEM(SonList[s]); i++)
5276 {
5277 SET_NBELEM(SonList[s],i,NULL);
5278
5279 const INT side = SON_NB(sdata,i);
5280 /* an interior face */
5281 if (side < FATHER_SIDE_OFFSET)
5282 {
5283 SET_NBELEM(SonList[s],i,SonList[side]);
5284
5285 IFDEBUG(gm,3)
5286 UserWriteF("elid=%3d: side:",ID(SonList[s]));
5287 for (INT p=0; p<CORNERS_OF_SIDE(SonList[s],i); p++)
5288 UserWriteF(" %2d",ID(CORNER_OF_SIDE_PTR(SonList[s],i,p)));
5289 UserWriteF(" INSIDE of father");
5290 UserWriteF("\nnbid=%3d: side:",ID(SonList[side]));
5291 {
5292 int ss,qq,pp,f,pts;
5293 f=0;
5294 for (ss=0; ss<SIDES_OF_ELEM(SonList[side]); ss++)
5295 {
5296 pts = 0;
5297 for (pp=0; pp<CORNERS_OF_SIDE(SonList[s],i); pp++)
5298 for (qq=0; qq<CORNERS_OF_SIDE(SonList[side],ss); qq++)
5299 if (CORNER_OF_SIDE_PTR(SonList[s],i,pp) ==
5300 CORNER_OF_SIDE_PTR(SonList[side],ss,qq))
5301 {
5302 pts |= ((1<<pp) | (16<<qq));
5303 break;
5304 }
5305 switch (pts)
5306 {
5307 #ifdef __TWODIM__
5308 case (LINEPOINTS) :
5309 #endif
5310 #ifdef __THREEDIM__
5311 case (TRIPOINTS) :
5312 case (QUADPOINTS) :
5313 #endif
5314 if (pts==TRIPOINTS &&
5315 CORNERS_OF_SIDE(SonList[s],i)==4)
5316 {
5317 PrintErrorMessage('E',"RefineElement",
5318 "quad side with 3 equal nodes");
5319 RETURN(GM_FATAL);
5320 }
5321 f=1;
5322 break;
5323 default :
5324 break;
5325 }
5326 if (f) break;
5327 }
5328 ASSERT(f==1);
5329 for (pp=0; pp<CORNERS_OF_SIDE(SonList[side],ss); pp++)
5330 UserWriteF(" %2d",ID(CORNER_OF_SIDE_PTR(SonList[side],ss,pp)));
5331 }
5332 UserWriteF("\n\n");
5333 ENDDEBUG
5334
5335 ASSERT(SonList[side]!=NULL);
5336
5337 /* dispose doubled side vectors if */
5338 #ifdef __THREEDIM__
5339 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,SIDEVEC))
5340 {
5341 INT l;
5342 for (l=0; l<SIDES_OF_ELEM(SonList[side]); l++)
5343 if (NBELEM(SonList[side],l)==SonList[s])
5344 break;
5345
5346 if (l<SIDES_OF_ELEM(SonList[side]))
5347 {
5348 /* assert consistency of rule set */
5349 ASSERT(SON_NB_OF_RULE(rule,side,l)==s);
5350 ASSERT(SON_NB_OF_RULE(rule,s,i)==side);
5351 ASSERT(NBELEM(SonList[s],i)==SonList[side]
5352 && NBELEM(SonList[side],l)==SonList[s]);
5353 if (DisposeDoubledSideVector(theGrid,SonList[s],i,
5354 SonList[side],l)) RETURN(GM_FATAL);
5355 }
5356 }
5357 #endif
5358 continue;
5359 }
5360 }
5361 }
5362
5363 IFDEBUG(gm,2)
5364 UserWriteF(PFMT "CONNECTING elem=" EID_FMTX "\n",me,EID_PRTX(theElement));
5365 ENDDEBUG
5366 for (INT i=0; i<SIDES_OF_ELEM(theElement); i++)
5367 {
5368 INT Sons_of_Side;
5369 ELEMENT *Sons_of_Side_List[MAX_SONS];
5370 INT SonSides[MAX_SIDE_NODES];
5371
5372 IFDEBUG(gm,2)
5373 UserWriteF(PFMT " CONNECT side=%i of elem=" EID_FMTX "\n",me,i,EID_PRTX(theElement));
5374 ENDDEBUG
5375
5376 for (INT j=0; j<MAX_SONS; j++)
5377 Sons_of_Side_List[j] = NULL;
5378
5379 for (INT j=0; j<NSONS_OF_RULE(rule); j++)
5380 Sons_of_Side_List[j] = SonList[j];
5381
5382 if (Get_Sons_of_ElementSide(theElement,i,&Sons_of_Side,
5383 Sons_of_Side_List,SonSides,0,0)!=GM_OK) RETURN(GM_FATAL);
5384
5385 if (Connect_Sons_of_ElementSide(theGrid,theElement,i,Sons_of_Side,
5386 Sons_of_Side_List,SonSides,0)!=GM_OK) RETURN(GM_FATAL);
5387
5388 #ifdef ModelP
5389 if (Identify_Objects_of_ElementSide(theGrid,theElement,i)) RETURN(GM_FATAL);
5390 #endif
5391 }
5392
5393 return(GM_OK);
5394 }
5395
5396
5397 /****************************************************************************/
5398 /*
5399 - refine an element
5400
5401 SYNOPSIS:
5402 static INT RefineElement (GRID *UpGrid, ELEMENT *theElement,NODE** theNodeContext);
5403
5404 PARAMETERS:
5405 \param UpGrid - grid level to refine
5406 \param theElement - element to refine
5407 \param theNodeContext - nodecontext for refinement
5408
5409 DESCRIPTION:
5410 This function refines an element
5411 \return <ul>
5412 INT
5413 .n GM_OK - ok
5414 .n GM_FATAL - fatal memory error
5415 */
5416 /****************************************************************************/
5417
RefineElement(GRID * UpGrid,ELEMENT * theElement,NODE ** theNodeContext)5418 static INT RefineElement (GRID *UpGrid, ELEMENT *theElement,NODE** theNodeContext)
5419 {
5420 switch (MARKCLASS(theElement))
5421 {
5422 case (YELLOW_CLASS) :
5423 if (RefineElementYellow(UpGrid,theElement,theNodeContext)!=GM_OK)
5424 RETURN(GM_FATAL);
5425 break;
5426
5427 case (GREEN_CLASS) :
5428 #ifdef __ANISOTROPIC__
5429 case (RED_CLASS) :
5430 #endif
5431 if (MARKED_NEW_GREEN(theElement))
5432 {
5433 /* elements with incomplete rules set */
5434 if (RefineElementGreen(UpGrid,theElement,theNodeContext) != GM_OK)
5435 RETURN(GM_FATAL);
5436 }
5437 else
5438 {
5439 /* elements with complete rules set */
5440 if (RefineElementRed(UpGrid,theElement,theNodeContext)!=GM_OK)
5441 RETURN(GM_FATAL);
5442 }
5443 break;
5444
5445 #ifndef __ANISOTROPIC__
5446 case (RED_CLASS) :
5447 if (RefineElementRed(UpGrid,theElement,theNodeContext)!=GM_OK)
5448 RETURN(GM_FATAL);
5449 break;
5450 #endif
5451
5452 default :
5453 RETURN(GM_FATAL);
5454 }
5455
5456 return(GM_OK);
5457 }
5458
5459
5460 /****************************************************************************/
5461 /*
5462 AdaptGrid - adapt one level of the multigrid
5463
5464 SYNOPSIS:
5465 static int AdaptGrid (GRID *theGrid, INT *nadapted)
5466
5467 PARAMETERS:
5468 \param theGrid - grid level to refine
5469 \param nadapted - number elements have been changed
5470
5471 DESCRIPTION:
5472 This function refines one level of the grid
5473
5474 \return <ul>
5475 INT
5476 .n GM_OK - ok
5477 .n GM_FATAL - fatal memory error
5478 */
5479 /****************************************************************************/
5480
5481 #ifdef ModelP
AdaptLocalGrid(GRID * theGrid,INT * nadapted)5482 static int AdaptLocalGrid (GRID *theGrid, INT *nadapted)
5483 #else
5484 static int AdaptGrid (GRID *theGrid, INT *nadapted)
5485 #endif
5486 {
5487 INT modified = 0;
5488 ELEMENT *theElement;
5489 ELEMENT *NextElement;
5490 ELEMENTCONTEXT theContext;
5491 GRID *UpGrid;
5492 const int me = theGrid->ppifContext().me();
5493 #ifdef ModelP
5494 auto& dddContext = theGrid->dddContext();
5495 #endif
5496
5497 UpGrid = UPGRID(theGrid);
5498 if (UpGrid==NULL)
5499 RETURN(GM_FATAL);
5500
5501 REFINE_GRID_LIST(1,MYMG(theGrid),GLEVEL(theGrid),("AdaptGrid(%d):\n",GLEVEL(theGrid)),"");
5502
5503 #ifdef IDENT_ONLY_NEW
5504 /* reset ident flags for old objects */
5505 {
5506 INT i;
5507 NODE *theNode;
5508 EDGE *theEdge;
5509
5510 for (theNode=PFIRSTNODE(UpGrid); theNode!=NULL; theNode=SUCCN(theNode))
5511 SETNEW_NIDENT(theNode,0);
5512
5513 for (theElement=PFIRSTELEMENT(UpGrid); theElement!=NULL; theElement=SUCCE(theElement))
5514 {
5515 for (i=0; i<EDGES_OF_ELEM(theElement); i++)
5516 {
5517 theEdge = GetEdge(CORNER_OF_EDGE_PTR(theElement,i,0),
5518 CORNER_OF_EDGE_PTR(theElement,i,1));
5519 SETNEW_EDIDENT(theEdge,0);
5520 }
5521 }
5522 }
5523 #endif
5524
5525 /* refine elements */
5526 /* ModelP: first loop over master elems, then loop over ghost elems */
5527 /* this assures that no unnecessary disposures of objects are done */
5528 /* which may cause trouble during identification (s.l. 9803020 */
5529 for (theElement=FIRSTELEMENT(theGrid); theElement!=NULL;
5530 theElement=NextElement)
5531 {
5532 NextElement = SUCCE(theElement);
5533 #ifdef ModelP
5534 /* loop over master elems first, then over ghost elems */
5535 if (NextElement == NULL) NextElement=PFIRSTELEMENT(theGrid);
5536 if (NextElement == FIRSTELEMENT(theGrid)) NextElement = NULL;
5537 #endif
5538
5539 #ifdef ModelP
5540 /* reset update overlap flag */
5541 SETTHEFLAG(theElement,0);
5542 #endif
5543
5544 /* do not change PrioVGhost elements */
5545 if (EVGHOST(theElement)) continue;
5546
5547 if (REFINEMENT_CHANGES(theElement))
5548 {
5549 #ifdef ModelP
5550 {
5551 /* check for valid load balancing */
5552 int *proclist = EPROCLIST(dddContext, theElement);
5553 proclist += 2;
5554 while (*proclist != -1)
5555 {
5556 if (*(proclist+1)!=PrioMaster && (*(proclist+1)!=PrioHGhost))
5557 {
5558 UserWriteF(PFMT "ERROR invalid load balancing: element="
5559 EID_FMTX " has copies of type=%d\n",
5560 me,EID_PRTX(theElement),*(proclist+1));
5561 REFINE_ELEMENT_LIST(0,theElement,"ERROR element: ");
5562 assert(0);
5563 }
5564 proclist += 2;
5565 }
5566 }
5567 #endif
5568
5569 if (hFlag==0 && MARKCLASS(theElement)!=RED_CLASS)
5570 {
5571 /* remove copy marks */
5572 SETMARK(theElement,NO_REFINEMENT);
5573 SETMARKCLASS(theElement,NO_CLASS);
5574 /* continue; */
5575 }
5576
5577 REFINE_ELEMENT_LIST(1,theElement,"REFINING element: ");
5578
5579 if (UnrefineElement(UpGrid,theElement))
5580 RETURN(GM_FATAL);
5581
5582 #ifdef ModelP
5583 /* dispose hghost elements with EFATHER==NULL */
5584 /** \todo how handle this situation? */
5585 /* situation possibly some elements to be coarsened are */
5586 /* disconnected from their fathers (970109 s.l.) */
5587 if (0)
5588 if (EHGHOST(theElement) && COARSEN(theElement))
5589 {
5590 if (LEVEL(theElement)>0 && EFATHER(theElement)==NULL)
5591 {
5592 DisposeElement(theGrid,theElement,true);
5593 continue;
5594 }
5595 }
5596 #endif
5597
5598 if (EMASTER(theElement))
5599 {
5600 if (UpdateContext(UpGrid,theElement,theContext)!=0)
5601 RETURN(GM_FATAL);
5602
5603 REFINE_CONTEXT_LIST(2,theContext);
5604
5605 #ifdef Debug
5606 CheckElementContextConsistency(theElement,theContext);
5607 #endif
5608
5609 /* is something to do ? */
5610 if (MARKED(theElement))
5611 if (RefineElement(UpGrid,theElement,theContext))
5612 RETURN(GM_FATAL);
5613 }
5614
5615 /* refine and refineclass flag */
5616 SETREFINE(theElement,MARK(theElement));
5617 SETREFINECLASS(theElement,MARKCLASS(theElement));
5618 SETUSED(theElement,0);
5619
5620 #ifdef ModelP
5621 /* set update overlap flag */
5622 SETTHEFLAG(theElement,1);
5623 #endif
5624
5625 /* this grid is modified */
5626 modified++;
5627 }
5628 else
5629 {
5630 #ifdef ModelP
5631 /* dispose hghost elements with EFATHER==NULL */
5632 /** \todo how handle this situation? */
5633 /* situation possibly some elements to be coarsened are */
5634 /* disconnected from their fathers (970109 s.l.) */
5635 if (0)
5636 if (EHGHOST(theElement) && COARSEN(theElement))
5637 {
5638 if (LEVEL(theElement)>0 && EFATHER(theElement)==NULL)
5639 {
5640 DisposeElement(theGrid,theElement,true);
5641 continue;
5642 }
5643 }
5644 #endif
5645
5646 #ifdef __ANISOTROPIC__
5647 if (USED(theElement)==0 && MARKCLASS(theElement)==GREEN_CLASS)
5648 #else
5649 if (USED(theElement)==0)
5650 #endif
5651 {
5652 /* count not updated green refinements */
5653 No_Green_Update++;
5654 }
5655 }
5656
5657 /* count green marks */
5658 if (MARKCLASS(theElement) == GREEN_CLASS) Green_Marks++;
5659
5660 /* reset coarse flag */
5661 SETCOARSEN(theElement,0);
5662 }
5663
5664 if (UG_GlobalMaxINT(theGrid->ppifContext(), modified))
5665 {
5666 /* reset (multi)grid status */
5667 SETGLOBALGSTATUS(UpGrid);
5668 RESETMGSTATUS(MYMG(UpGrid));
5669 }
5670 REFINE_GRID_LIST(1,MYMG(theGrid),GLEVEL(theGrid),
5671 ("END AdaptGrid(%d):\n",GLEVEL(theGrid)),"");
5672
5673 *nadapted = modified;
5674
5675 return(GM_OK);
5676 }
5677
5678 #ifdef ModelP
AdaptGrid(GRID * theGrid,INT toplevel,INT level,INT newlevel,INT * nadapted)5679 static int AdaptGrid (GRID *theGrid, INT toplevel, INT level, INT newlevel, INT *nadapted)
5680 {
5681 GRID *FinerGrid = UPGRID(theGrid);
5682
5683 START_TIMER(gridadapti_timer)
5684
5685 #ifdef UPDATE_FULLOVERLAP
5686 DDD_XferBegin(theGrid->dddContext());
5687 {
5688 ELEMENT *theElement;
5689 for (theElement=PFIRSTELEMENT(FinerGrid);
5690 theElement!=NULL;
5691 theElement=SUCCE(theElement))
5692 {
5693 if (EPRIO(theElement) == PrioHGhost) DisposeElement(FinerGrid,theElement,true);
5694 }
5695 }
5696 DDD_XferEnd(theGrid->dddContext());
5697 #endif
5698
5699 DDD_IdentifyBegin(theGrid->dddContext());
5700 SET_IDENT_MODE(IDENT_ON);
5701 DDD_XferBegin(theGrid->dddContext());
5702
5703 DDD_CONSCHECK(theGrid->dddContext());
5704
5705 /* now really manipulate the next finer level */
5706 START_TIMER(gridadaptl_timer)
5707
5708 #ifdef DDDOBJMGR
5709 DDD_ObjMgrBegin();
5710 #endif
5711 if (level<toplevel || newlevel)
5712 if (AdaptLocalGrid(theGrid,nadapted)!=GM_OK)
5713 RETURN(GM_FATAL);
5714 #ifdef DDDOBJMGR
5715 DDD_ObjMgrEnd();
5716 #endif
5717
5718 DDD_XferEnd(theGrid->dddContext());
5719
5720 SUM_TIMER(gridadaptl_timer)
5721
5722 DDD_CONSCHECK(theGrid->dddContext());
5723
5724 {
5725 int check=1;
5726 int debugstart=3;
5727 #ifdef Debug
5728 int gmlevel=Debuggm;
5729 #else
5730 int gmlevel=0;
5731 #endif
5732
5733
5734 if (IDENT_IN_STEPS)
5735 {
5736 DDD_IdentifyEnd(theGrid->dddContext());
5737 }
5738
5739 /* if no grid adaption has occured adapt next level */
5740 *nadapted = UG_GlobalSumINT(theGrid->ppifContext(), *nadapted);
5741 if (*nadapted == 0)
5742 {
5743 if (!IDENT_IN_STEPS)
5744 {
5745 SET_IDENT_MODE(IDENT_OFF);
5746 DDD_IdentifyEnd(theGrid->dddContext());
5747 }
5748
5749 SUM_TIMER(gridadapti_timer)
5750
5751 return(GM_OK);
5752 }
5753
5754 /* DDD_JoinBegin(); */
5755 if (IDENT_IN_STEPS)
5756 DDD_IdentifyBegin(theGrid->dddContext());
5757
5758 DDD_CONSCHECK(theGrid->dddContext());
5759
5760 START_TIMER(ident_timer)
5761
5762 if (Identify_SonObjects(theGrid)) RETURN(GM_FATAL);
5763
5764 SET_IDENT_MODE(IDENT_OFF);
5765 DDD_IdentifyEnd(theGrid->dddContext());
5766
5767 SUM_TIMER(ident_timer)
5768 /* DDD_JoinEnd(); */
5769
5770
5771 DDD_CONSCHECK(theGrid->dddContext());
5772
5773 if (level<toplevel || newlevel)
5774 {
5775 START_TIMER(overlap_timer)
5776 DDD_XferBegin(theGrid->dddContext());
5777 if (0) /* delete sine this is already done in */
5778 /* ConstructConsistentGrid() (s.l. 980522) */
5779 if (SetGridBorderPriorities(theGrid)) RETURN(GM_FATAL);
5780 if (UpdateGridOverlap(theGrid)) RETURN(GM_FATAL);
5781
5782
5783 DDD_XferEnd(theGrid->dddContext());
5784
5785 DDD_CONSCHECK(theGrid->dddContext());
5786
5787 DDD_XferBegin(theGrid->dddContext());
5788 if (ConnectGridOverlap(theGrid)) RETURN(GM_FATAL);
5789 DDD_XferEnd(theGrid->dddContext());
5790
5791 DDD_CONSCHECK(theGrid->dddContext());
5792
5793 /* this is needed due to special cases while coarsening */
5794 /* sample scene: a ghost element is needed as overlap */
5795 /* for two master elements, one of the master elements */
5796 /* is coarsened, then the prio of nodes of the ghost */
5797 /* element must eventually be downgraded from master */
5798 /* to ghost prio (s.l. 971020) */
5799 /* this is done as postprocessing step, since this needs */
5800 /* 2 XferBegin/Ends here per modified gridlevel (s.l. 980522) */
5801 /*
5802 #ifndef NEW_GRIDCONS_STYLE
5803 ConstructConsistentGrid(FinerGrid);
5804 #endif
5805 */
5806 SUM_TIMER(overlap_timer)
5807 }
5808
5809 DDD_CONSCHECK(theGrid->dddContext());
5810
5811
5812 CheckConsistency(MYMG(theGrid),(INT)level,(INT)debugstart,(INT)gmlevel,&check);
5813
5814
5815 }
5816
5817 if (0) CheckGrid(FinerGrid,1,0,1,1);
5818
5819 SUM_TIMER(gridadapti_timer)
5820
5821 return(GM_OK);
5822 }
5823 #endif
5824
5825
5826 #ifdef ModelP
5827 /* parameters for CheckGrid() */
5828 #define GHOSTS 1
5829 #define GEOM 1
5830 #define ALG 0
5831 #define LIST 1
5832 #define IF 1
5833
5834
5835 /****************************************************************************/
5836 /*
5837 CheckConsistency -
5838
5839 SYNOPSIS:
5840 void CheckConsistency (MULTIGRID *theMG, INT level ,INT debugstart, INT gmlevel, INT *check);
5841
5842 PARAMETERS:
5843 . theMG
5844 . level
5845 . debugstart
5846 . gmlevel
5847 . check
5848
5849 */
5850 /****************************************************************************/
5851
CheckConsistency(MULTIGRID * theMG,INT level,INT debugstart,INT gmlevel,int * check)5852 static void CheckConsistency (MULTIGRID *theMG, INT level ,INT debugstart, INT gmlevel, int *check)
5853 {
5854 GRID *theGrid = GRID_ON_LEVEL(theMG,level);
5855
5856 IFDEBUG(gm,debugstart)
5857 printf(PFMT "AdaptMultiGrid(): %d. ConsCheck() on level=%d\n",theMG->ppifContext().me(),(*check)++,level);
5858 #ifdef Debug
5859 Debuggm = GHOSTS;
5860 #endif
5861 CheckGrid(theGrid,GEOM,ALG,LIST,IF);
5862 #ifdef Debug
5863 Debuggm=gmlevel;
5864 #endif
5865 if (DDD_ConsCheck(theMG->dddContext()) > 0) buggy(theMG);
5866 ENDDEBUG
5867 }
5868 #endif
5869
5870
CheckMultiGrid(MULTIGRID * theMG)5871 static INT CheckMultiGrid (MULTIGRID *theMG)
5872 {
5873 INT level;
5874
5875 UserWriteF("CheckMultiGrid() begin\n");
5876
5877 for (level=0; level<=TOPLEVEL(theMG); level++)
5878 #ifdef ModelP
5879 CheckGrid(GRID_ON_LEVEL(theMG,level),1,1,1,1);
5880 #else
5881 CheckGrid(GRID_ON_LEVEL(theMG,level),1,1,1);
5882 #endif
5883
5884 UserWriteF("CheckMultiGrid() end\n");
5885
5886 return(0);
5887 }
5888
5889
5890 #ifdef STAT_OUT
Manage_Adapt_Timer(int alloc)5891 void NS_DIM_PREFIX Manage_Adapt_Timer (int alloc)
5892 {
5893 if (alloc)
5894 {
5895 NEW_TIMER(adapt_timer)
5896 NEW_TIMER(closure_timer)
5897 NEW_TIMER(gridadapt_timer)
5898 NEW_TIMER(gridadapti_timer)
5899 NEW_TIMER(gridadaptl_timer)
5900 NEW_TIMER(ident_timer)
5901 NEW_TIMER(overlap_timer)
5902 NEW_TIMER(gridcons_timer)
5903 NEW_TIMER(algebra_timer)
5904 }
5905 else
5906 {
5907 DEL_TIMER(adapt_timer)
5908 DEL_TIMER(closure_timer)
5909 DEL_TIMER(gridadapt_timer)
5910 DEL_TIMER(gridadapti_timer)
5911 DEL_TIMER(gridadaptl_timer)
5912 DEL_TIMER(ident_timer)
5913 DEL_TIMER(overlap_timer)
5914 DEL_TIMER(gridcons_timer)
5915 DEL_TIMER(algebra_timer)
5916 }
5917 }
5918
Print_Adapt_Timer(const MULTIGRID * theMG,int total_adapted)5919 void NS_DIM_PREFIX Print_Adapt_Timer (const MULTIGRID* theMG, int total_adapted)
5920 {
5921 UserWriteF("ADAPT: total_adapted=%d t_adapt=%.2f: t_closure=%.2f t_gridadapt=%.2f t_gridadapti=%.2f "
5922 "t_gridadaptl=%.2f t_overlap=%.2f t_ident=%.2f t_gridcons=%.2f t_algebra=%.2f\n",
5923 total_adapted,EVAL_TIMER(adapt_timer),EVAL_TIMER(closure_timer),EVAL_TIMER(gridadapt_timer),
5924 EVAL_TIMER(gridadapti_timer),EVAL_TIMER(gridadaptl_timer),EVAL_TIMER(overlap_timer),
5925 EVAL_TIMER(ident_timer),EVAL_TIMER(gridcons_timer),EVAL_TIMER(algebra_timer));
5926 UserWriteF("ADAPTMAX: total_adapted=%d t_adapt=%.2f: t_closure=%.2f t_gridadapt=%.2f "
5927 "t_gridadapti=%.2f "
5928 "t_gridadaptl=%.2f t_overlap=%.2f t_ident=%.2f t_gridcons=%.2f t_algebra=%.2f\n",
5929 total_adapted,UG_GlobalMaxDOUBLE(theMG->ppifContext(), EVAL_TIMER(adapt_timer)),
5930 UG_GlobalMaxDOUBLE(theMG->ppifContext(), EVAL_TIMER(closure_timer)),
5931 UG_GlobalMaxDOUBLE(theMG->ppifContext(), EVAL_TIMER(gridadapt_timer)),UG_GlobalMaxDOUBLE(theMG->ppifContext(), EVAL_TIMER(gridadapti_timer)),
5932 UG_GlobalMaxDOUBLE(theMG->ppifContext(), EVAL_TIMER(gridadaptl_timer)),
5933 UG_GlobalMaxDOUBLE(theMG->ppifContext(), EVAL_TIMER(overlap_timer)),
5934 UG_GlobalMaxDOUBLE(theMG->ppifContext(), EVAL_TIMER(ident_timer)),UG_GlobalMaxDOUBLE(theMG->ppifContext(), EVAL_TIMER(gridcons_timer)),
5935 UG_GlobalMaxDOUBLE(theMG->ppifContext(), EVAL_TIMER(algebra_timer)));
5936 }
5937 #endif
5938
PreProcessAdaptMultiGrid(MULTIGRID * theMG)5939 static INT PreProcessAdaptMultiGrid(MULTIGRID *theMG)
5940 {
5941 if (DisposeBottomHeapTmpMemory(theMG)) REP_ERR_RETURN(1);
5942
5943 /* The matrices for the calculation are removed, to remember the
5944 recalculating the MGSTATUS is set to 1 */
5945
5946 MGSTATUS(theMG) = 1 ;
5947
5948 return(0);
5949 }
5950
PostProcessAdaptMultiGrid(MULTIGRID * theMG)5951 static INT PostProcessAdaptMultiGrid(MULTIGRID *theMG)
5952 {
5953 START_TIMER(algebra_timer)
5954 if (CreateAlgebra(theMG)) REP_ERR_RETURN(1);
5955 SUM_TIMER(algebra_timer)
5956
5957 REFINE_MULTIGRID_LIST(1,theMG,"END AdaptMultiGrid():\n","","");
5958
5959 /*
5960 if (hFlag)
5961 UserWriteF(" Number of green refinements not updated: "
5962 "%d (%d green marks)\n",No_Green_Update,Green_Marks);
5963 */
5964
5965 /* increment step count */
5966 SETREFINESTEP(REFINEINFO(theMG),REFINESTEP(REFINEINFO(theMG))+1);
5967
5968 SUM_TIMER(adapt_timer)
5969 /*
5970 CheckMultiGrid(theMG);
5971 */
5972
5973 #ifdef STAT_OUT
5974 Print_Adapt_Timer(theMG, total_adapted);
5975 Manage_Adapt_Timer(0);
5976 #endif
5977
5978 return(0);
5979 }
5980
5981 /****************************************************************************/
5982 /** \brief Adapt whole multigrid structure
5983
5984 \param theMG - multigrid to refine
5985 \param flag - flag for switching between different yellow closures
5986
5987 This function refines whole multigrid structure
5988
5989 \return <ul>
5990 <li> 0 - ok
5991 <li> 1 - out of memory, but data structure as before
5992 <li> 2 - fatal memory error, data structure corrupted
5993 </ul>
5994 */
5995 /****************************************************************************/
5996
AdaptMultiGrid(MULTIGRID * theMG,INT flag,INT seq,INT mgtest)5997 INT NS_DIM_PREFIX AdaptMultiGrid (MULTIGRID *theMG, INT flag, INT seq, INT mgtest)
5998 {
5999 INT level,toplevel,nrefined,nadapted;
6000 INT newlevel;
6001 NODE *theNode;
6002 GRID *theGrid, *FinerGrid;
6003 ELEMENT *theElement;
6004
6005 /* check necessary condition */
6006 if (!MG_COARSE_FIXED(theMG))
6007 return (GM_COARSE_NOT_FIXED);
6008
6009 if (PreProcessAdaptMultiGrid(theMG)) REP_ERR_RETURN(1);
6010
6011 #ifdef ModelP
6012 {
6013 /* check and restrict partitioning of elements */
6014 if (CheckPartitioning(theMG))
6015 {
6016 /* Each call to RestrictPartitioning fixes the partitionings of children with
6017 * respect to their fathers. To fix also fix the partitionings of the grandchildren,
6018 * the method has to be called again. To be on the save side we call it once for
6019 * every level.
6020 * If I omit the loop (the way it was until 11.4.2014 I get assertion failures here
6021 * when mixing load balancing and adaptive refinement. I am not really sure whether
6022 * the loop is the correct fix, or whether it just papers over the problem.
6023 * Anyway, no crashes for now.
6024 */
6025 for (level=0; level<TOPLEVEL(theMG); level++)
6026 if (RestrictPartitioning(theMG)) RETURN(GM_FATAL);
6027 if (CheckPartitioning(theMG)) assert(0);
6028 }
6029 }
6030 #endif
6031
6032 #ifdef STAT_OUT
6033 Manage_Adapt_Timer(1);
6034 #endif
6035
6036 START_TIMER(adapt_timer)
6037
6038 /* set up information in refine_info */
6039 #ifndef ModelP
6040 if (TOPLEVEL(theMG) == 0)
6041 #else
6042 if (UG_GlobalMaxINT(theMG->ppifContext(), TOPLEVEL(theMG)) == 0)
6043 #endif
6044 {
6045 SETREFINESTEP(REFINEINFO(theMG),0);
6046 }
6047
6048 /* set info for refinement prediction */
6049 SetRefineInfo(theMG);
6050 /* evaluate prediction */
6051 if (mgtest)
6052 {
6053 UserWriteF("refinetest: predicted_new0=%9.0f predicted_new1=%9.0f"
6054 " predicted_max=%9.0f\n",
6055 PREDNEW0(REFINEINFO(theMG)), PREDNEW1(REFINEINFO(theMG)),
6056 PREDMAX(REFINEINFO(theMG)));
6057
6058 /* refinement possible or not */
6059 if (TestRefineInfo(theMG) != GM_OK)
6060 {
6061 UserWriteF("Too much marked elements: "
6062 "Number of marked elements would cause heap overflow\n");
6063 return(GM_ERROR);
6064 }
6065 }
6066
6067 /* set flags for different modes */
6068 rFlag=flag & 0x03; /* copy local or all */
6069 hFlag=!((flag>>2)&0x1); /* use hanging nodes */
6070 fifoFlag=(flag>>3)&0x1; /* use fifo */
6071
6072 refine_seq = seq;
6073
6074 No_Green_Update=0;
6075 Green_Marks=0;
6076
6077 /* drop marks to regular elements */
6078 if (hFlag)
6079 if (DropMarks(theMG)) RETURN(GM_ERROR);
6080
6081 /* prepare algebra (set internal flags correctly) */
6082 START_TIMER(algebra_timer)
6083
6084 PrepareAlgebraModification(theMG);
6085
6086 SUM_TIMER(algebra_timer)
6087
6088 toplevel = TOPLEVEL(theMG);
6089
6090 REFINE_MULTIGRID_LIST(1,theMG,"AdaptMultiGrid()","","")
6091
6092 /* compute modification of coarser levels from above */
6093 START_TIMER(closure_timer)
6094
6095 for (level=toplevel; level>0; level--)
6096 {
6097 theGrid = GRID_ON_LEVEL(theMG,level);
6098
6099 if (hFlag)
6100 {
6101 PRINTDEBUG(gm,1,("Begin GridClosure(%d,down):\n",level))
6102
6103 if ((nrefined = GridClosure(GRID_ON_LEVEL(theMG,level)))<0)
6104 {
6105 PrintErrorMessage('E',"AdaptMultiGrid","error in GridClosure");
6106 RETURN(GM_ERROR);
6107 }
6108
6109 REFINE_GRID_LIST(1,theMG,level,("End GridClosure(%d,down):\n",level),"");
6110 }
6111 #ifdef ModelP
6112 else
6113 {
6114 ExchangeElementRefine(theGrid);
6115 }
6116 #endif
6117
6118 /* restrict marks on next lower grid level */
6119 if (RestrictMarks(GRID_ON_LEVEL(theMG,level-1))!=GM_OK) RETURN(GM_ERROR);
6120
6121 REFINE_GRID_LIST(1,theMG,level-1,("End RestrictMarks(%d,down):\n",level),"");
6122 }
6123
6124 SUM_TIMER(closure_timer)
6125
6126
6127 #ifdef ModelP
6128 IdentifyInit(theMG);
6129 #endif
6130
6131 newlevel = 0;
6132 for (level=0; level<=toplevel; level++)
6133 {
6134 theGrid = GRID_ON_LEVEL(theMG,level);
6135 if (level<toplevel) FinerGrid = GRID_ON_LEVEL(theMG,level+1);else FinerGrid = NULL;
6136
6137 START_TIMER(closure_timer)
6138
6139 /* reset MODIFIED flags for grid and nodes */
6140 SETMODIFIED(theGrid,0);
6141 for (theNode=FIRSTNODE(theGrid); theNode!=NULL; theNode=SUCCN(theNode)) SETMODIFIED(theNode,0);
6142
6143 if (hFlag)
6144 {
6145 /* leave only regular marks */
6146 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL; theElement=SUCCE(theElement))
6147 {
6148 if ((ECLASS(theElement)==RED_CLASS) && MARKCLASS(theElement)==RED_CLASS) continue;
6149 SETMARK(theElement,NO_REFINEMENT);
6150 }
6151
6152 PRINTDEBUG(gm,1,("Begin GridClosure(%d,up):\n",level));
6153
6154 /* determine regular and irregular elements on next level */
6155 if ((nrefined = GridClosure(theGrid))<0)
6156 {
6157 PrintErrorMessage('E',"AdaptMultiGrid","error in 2. GridClosure");
6158 RETURN(GM_ERROR);
6159 }
6160
6161 REFINE_GRID_LIST(1,theMG,level,("End GridClosure(%d,up):\n",level),"");
6162 }
6163 #ifdef ModelP
6164 else
6165 {
6166 ExchangeElementRefine(theGrid);
6167 }
6168 #endif
6169
6170 nrefined += ComputeCopies(theGrid);
6171
6172 if (hFlag)
6173 {
6174 /* dispose connections that may be changed on next level, this is determined */
6175 /* by the neighborhood of elements were MARK != REFINE. */
6176 /* This will leave some flags where to rebuild connections later */
6177 if (level<toplevel)
6178 {
6179 for (theElement=FIRSTELEMENT(FinerGrid); theElement!=NULL; theElement=SUCCE(theElement))
6180 {
6181 ELEMENT *theFather = EFATHER(theElement);
6182 ASSERT(theFather != NULL);
6183 /*
6184 if (REFINE(EFATHER(theElement))!=MARK(EFATHER(theElement)))
6185 */
6186 if (REFINEMENT_CHANGES(theFather))
6187 {
6188 ASSERT(EMASTER(theFather));
6189 if (DisposeConnectionsInNeighborhood(FinerGrid,theElement)!=GM_OK)
6190 RETURN(GM_FATAL);
6191 }
6192 }
6193 }
6194 }
6195
6196 /** \todo bug fix to force new level creation */
6197 if (!hFlag)
6198 {
6199 /* set this variable>0 */
6200 nrefined = 1;
6201 }
6202
6203 /* create a new grid level, if at least one element is refined on finest level */
6204 if (nrefined>0 && level==toplevel) newlevel = 1;
6205 #ifdef ModelP
6206 newlevel = UG_GlobalMaxINT(theMG->ppifContext(), newlevel);
6207 #endif
6208 if (newlevel)
6209 {
6210 if (CreateNewLevel(theMG)==NULL)
6211 RETURN(GM_FATAL);
6212 FinerGrid = GRID_ON_LEVEL(theMG,toplevel+1);
6213 }
6214
6215 PRINTDEBUG(gm,1,(PFMT "AdaptMultiGrid(): toplevel=%d nrefined=%d newlevel=%d\n",
6216 me,toplevel,nrefined,newlevel));
6217
6218 SUM_TIMER(closure_timer)
6219
6220 /* now really manipulate the next finer level */
6221 START_TIMER(gridadapt_timer)
6222
6223 nadapted = 0;
6224
6225 if (level<toplevel || newlevel)
6226 #ifndef ModelP
6227 if (AdaptGrid(theGrid,&nadapted)!=GM_OK)
6228 RETURN(GM_FATAL);
6229 #else
6230 if (AdaptGrid(theGrid,toplevel,level,newlevel,&nadapted)!=GM_OK)
6231 RETURN(GM_FATAL);
6232 #endif
6233
6234 SUM_TIMER(gridadapt_timer)
6235
6236 /* if no grid adaption has occured adapt next level */
6237 if (nadapted == 0) continue;
6238
6239 total_adapted += nadapted;
6240
6241 if (level<toplevel || newlevel)
6242 {
6243 START_TIMER(algebra_timer)
6244
6245 /* and compute the vector classes on the new (or changed) level */
6246 ClearNodeClasses(FinerGrid);
6247
6248 for (theElement=FIRSTELEMENT(FinerGrid); theElement!=NULL; theElement=SUCCE(theElement))
6249 if (ECLASS(theElement)>=GREEN_CLASS || (rFlag==GM_COPY_ALL)) {
6250 SeedNodeClasses(theElement);
6251 }
6252
6253 PropagateNodeClasses(FinerGrid);
6254
6255 SUM_TIMER(algebra_timer)
6256 }
6257 }
6258
6259 #ifdef ModelP
6260 IdentifyExit();
6261
6262 /* now repair inconsistencies */
6263 /* former done on each grid level (s.l. 980522) */
6264 START_TIMER(gridcons_timer);
6265
6266 ConstructConsistentMultiGrid(theMG);
6267
6268 SUM_TIMER(gridcons_timer)
6269 #endif
6270
6271 DisposeTopLevel(theMG);
6272 if (TOPLEVEL(theMG) > 0) DisposeTopLevel(theMG);
6273 CURRENTLEVEL(theMG) = TOPLEVEL(theMG);
6274
6275 if (PostProcessAdaptMultiGrid(theMG)) REP_ERR_RETURN(1);
6276
6277 return(GM_OK);
6278 }
6279