1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 /****************************************************************************/
4 /* */
5 /* File: er.c */
6 /* */
7 /* Purpose: extract rules realized in a multigrid */
8 /* */
9 /* Author: Henrik Rentz-Reichert */
10 /* Institut fuer Computeranwendungen III */
11 /* Universitaet Stuttgart */
12 /* Pfaffenwaldring 27 */
13 /* 70550 Stuttgart */
14 /* email: ug@ica3.uni-stuttgart.de */
15 /* */
16 /* History: 06.12.97 begin, ug version 3.7 */
17 /* */
18 /* Remarks: */
19 /* */
20 /****************************************************************************/
21
22 #undef __DEBUG_ER__
23 #undef __PATCH_UG_RULES__
24
25 /****************************************************************************/
26 /*
27 storage allocation stages:
28
29 Notation: M - Mark(Tmp)Mem
30 G - Get(Tmp)Mem
31 R - Release(Tmp)Mem
32
33
34 ------------------------+-+-+-+-+---+-+-+-+-+
35 |BOTTOM | | TOP |
36 ------------------------+-+-+-+-+---+-+-+-+-+
37 stack level |0|1|2|3| |0|1|2|3|
38 ------------------------+-+-+-+-+---+-+-+-+-+
39 (1) for hrules |M| | | | | | | | |
40 (2) hash table | | | | | | |M| | |
41 | | | | | | |G| | |
42 (3) interface rules | | | | | | | |M| |
43 | | | | | | | |G| |
44 (4) hrules |G| | | | | | | | |
45 (5) free (3) | | | | | | | |R| |
46 (6) further hrules |G| | | | | | | | |
47 (6) hrule table |G| | | | | | | | |
48 (7) free (2) | | | | | | |R| | |
49 (7) mrule table (rm+er) | | | | | |G| | | |
50 (9) free hrules |R| | | | | | | | |
51 ------------------------+-+-+-+-+---+-+-+-+-+
52 */
53 /****************************************************************************/
54
55 /****************************************************************************/
56 /* */
57 /* include files */
58 /* system include files */
59 /* application include files */
60 /* */
61 /****************************************************************************/
62
63 #include <config.h>
64 #include <cstring>
65
66 #include <algorithm>
67
68 /* gm */
69 #include "rm.h"
70 #include "mgio.h"
71 #include "ugm.h"
72 #include "cw.h"
73 #include "elements.h"
74
75 /* own header */
76 #include "er.h"
77
78
79 USING_UG_NAMESPACES
80
81 /****************************************************************************/
82 /* */
83 /* defines in the following order */
84 /* */
85 /* compile time constants defining static data size (i.e. arrays) */
86 /* other constants */
87 /* macros */
88 /* */
89 /****************************************************************************/
90
91 /* TODO (HRR 971207): increase REFINE_LEN by shifting REFINE-ce to flag cw */
92 #define MAX_HRID (1<<REFINE_LEN)
93 #define MAX_IFID MAX_HRID
94
95 #ifdef __PATCH_UG_RULES__
96 /* CAUTION: save/load ok but REFINE(e) will not be consistent with rm anymore */
97 #define UGMAXRULE(tag) ((MaxRules[tag]>0) ? 1 : 0) /* include NO_REFINEMENT from ug */
98 #define HAS_NO_RULE(e) (REFINE(e)>=UGMAXRULE(TAG(e)))
99 #else
100 #define UGMAXRULE(tag) (MaxRules[tag])
101 #define HAS_NO_RULE(e) (REFINE(e)==COPY && REFINECLASS(e)==GREEN_CLASS)
102 #endif
103
104 #define BEYOND_UG_RULES(e) (REFINE(e)>=UGMAXRULE(TAG(e)))
105
106 #define NOTDONE -1
107
108 /* for computing son paths */
109 enum NB_STATUS {
110
111 NB_DONE,
112 NB_NOTDONE,
113 NB_TOUCHED
114 };
115
116 /* debug levels of extract rule module */
117 enum ER_DBG {
118
119 ER_DBG_GENERAL = 1,
120 ER_DBG_RULES = 2,
121 ER_DBG_RULE_VERBOSE = 3,
122 ER_DBG_ELEM = 3
123 };
124
125 #define HASH_SIZE 1000
126 #define HASH_FACTOR .61803398874989 /* golden section: 0.5*(sqrt(5)-1) */
127 #define HASH_ADDRESS(k) floor(HASH_SIZE*(k*HASH_FACTOR-floor(k*HASH_FACTOR)))
128
129 #define ER_NSONS(p) ((p)->nsons)
130 #define ER_NCO(p,i) ((p)->nco[i])
131 #define ER_DCO(p,i) ((p)->dco[i])
132
133 #define HR_ID(p) ((p)->id)
134 #define HR_KEY(p) ((p)->key)
135 #define HR_TAG(p) ((p)->tag)
136 #define HR_NEXT(p) ((p)->next)
137 #define HR_ERULE(p) (&((p)->erule))
138 #define HR_NSONS(p) ((p)->erule.nsons)
139 #define HR_NCO(p,i) ((p)->erule.nco[i])
140 #define HR_DCO(p,i) ((p)->erule.dco[i])
141 #define HR_OCOPTR(p) ((p)->erule.dco+HR_NSONS(p))
142 #define HR_OCO(p,i) ((p)->erule.dco[HR_NSONS(p)+i])
143
144 /* debug macros */
145 #define ELEM_INFO(f,c,e,ns) f ": %sconsiderd elem%ld tag%d refine%ld nsons%d (actually %d)\n", \
146 (c) ? "not " : "", \
147 (long)ID(e), \
148 (int)TAG(e), \
149 (long)REFINE(e), \
150 (int)NSONS(e), \
151 (int)ns
152
153 #ifdef Debug
154 #define PD_ERR(l,x,e) {PRINTDEBUG(gm,l,x); e++; /* ASSERT(false);*/}
155 #define PD(x) {PrintDebug x;}
156 #else
157 #define PD_ERR(l,x,e) /* no debug */
158 #define PD(x) /* no debug */
159 #ifdef PrintDebug
160 #undef PrintDebug
161 #define PrintDebug printf
162 #endif
163 #endif
164
165 /****************************************************************************/
166 /* */
167 /* data structures used in this source file (exported data structures are */
168 /* in the corresponding include file!) */
169 /* */
170 /****************************************************************************/
171
172 typedef struct {
173
174 SHORT nsons; /* number of son elements */
175 SHORT nco[MAX_SONS]; /* number of son corners */
176 DOUBLE dco[MAX_SONS]; /* corners coded as digits to base */
177 /* MAX_REFINED_CORNERS_DIM */
178
179 } ERULE;
180
181 typedef INT HRID;
182
183 struct hashed_rule {
184
185 HRID id; /* rule id will continue tag-wise rm-rules */
186 DOUBLE key; /* hash key */
187 SHORT tag; /* tag of element */
188 struct hashed_rule *next; /* list of rules per hash entry */
189
190 ERULE erule; /* refinement rule */
191 /* will be allocated with additional doubles for ordered sons */
192 /* CAUTION: storage occupied may be < sizeof(HRULE) */
193 };
194
195 /* HRULE example for 3 sons:
196 id
197 key
198 tag
199 next
200 nsons // begin of ERULE
201 nco[MAX_SONS]
202 dco[3]
203 oco[3] // corners in ascending order
204 */
205
206 typedef struct hashed_rule HRULE;
207 typedef REFRULE URULE;
208
209 /****************************************************************************/
210 /* */
211 /* definition of exported global variables */
212 /* */
213 /****************************************************************************/
214
215
216
217 /****************************************************************************/
218 /* */
219 /* definition of variables global to this source file only (static!) */
220 /* */
221 /****************************************************************************/
222
223 static struct {
224
225 HEAP *heap; /* multigrid heap */
226
227 HRULE **hash_table; /* hash table */
228 HRULE **hrule[TAGS]; /* tables with hrules sorted by id */
229
230 long maxrule[TAGS]; /* max rule id per tag (rm+er) */
231 long maxrules; /* max rule id (rm+er) */
232 long nelem_inspected[TAGS]; /* elements getting er-rules per tag */
233 long nelems_inspected; /* elements getting er-rules */
234 long nelem_not_inspected[TAGS]; /* elements having rm-rules per tag */
235 long nelems_not_inspected; /* elements having rm-rules */
236
237 #ifdef ModelP
238 long if_elems; /* number of elements getting er-rules */
239 ERULE *interface_rules; /* table of interface rules */
240 #endif
241
242 } global; /* static globals of this module */
243
244 /****************************************************************************/
245 /* */
246 /* forward declarations of functions used before they are defined */
247 /* */
248 /****************************************************************************/
249
250
251 /****************************************************************************/
252 /*doctext_disabled
253 Corner2DCorners - map corner array to DOUBLE
254
255 SYNOPSIS:
256 static DOUBLE Corner2DCorners (INT n, SHORT corners[])
257
258 PAAMETERS:
259 . n - number of corners
260 . corners - array of corner IDs
261
262 DESCRIPTION:
263 Convert list of corner IDs into one DOUBLE value. Corner IDs are digits
264 in a MAX_REFINED_CORNERS_DIM base.
265
266 RETURN VALUE:
267 DOUBLE
268 .n coded coners
269 doctext_disabled*/
270 /****************************************************************************/
271
Corner2DCorners(INT n,SHORT corners[])272 static DOUBLE Corner2DCorners (INT n, SHORT corners[])
273 {
274 DOUBLE dco = corners[0];
275 int i;
276
277 for (i=1; i<n; i++)
278 {
279 dco *= MAX_REFINED_CORNERS_DIM;
280 dco += corners[i];
281 }
282 return (dco);
283 }
284
285 /****************************************************************************/
286 /*doctext_disabled
287 DCorners2Corners - convert DOUBLE into array of corner IDs
288
289 SYNOPSIS:
290 static void DCorners2Corners (INT n, DOUBLE dco, SHORT corners[])
291
292 PAAMETERS:
293 . n - number of corners
294 . dco - corners coded on DOUBLE
295 . corners - array of corner IDs
296
297 DESCRIPTION:
298 Convert DOUBLE into array of corner IDs. For encoding see 'Corner2DCorners'.
299
300 RETURN VALUE:
301 void
302 .n none
303
304 SEE_ASLO:
305 Corner2DCorners
306 doctext_disabled*/
307 /****************************************************************************/
308
DCorners2Corners(INT n,DOUBLE dco,SHORT corners[])309 static void DCorners2Corners (INT n, DOUBLE dco, SHORT corners[])
310 {
311 int i;
312
313 for (i=n-1; i>=0; i--)
314 {
315 DOUBLE x = floor(dco/((DOUBLE)MAX_REFINED_CORNERS_DIM));
316 corners[i] = (SHORT)(dco-x*MAX_REFINED_CORNERS_DIM);
317 dco = x;
318 }
319 return;
320 }
321
322 /****************************************************************************/
323 /*doctext_disabled
324 FillOrderedSons - fill array with ordered sons
325
326 SYNOPSIS:
327 static void FillOrderedSons (const ERULE *er, DOUBLE oco[])
328
329 PAAMETERS:
330 . er - rule
331 . oco - array of ordered sons
332
333 DESCRIPTION:
334 Fill array with sons in ascending order whose ordered corners are coded as DOUBLEs.
335
336 RETURN VALUE:
337 void
338 .n none
339 doctext_disabled*/
340 /****************************************************************************/
341
FillOrderedSons(const ERULE * er,DOUBLE oco[])342 static void FillOrderedSons (const ERULE *er, DOUBLE oco[])
343 {
344 SHORT corners[MAX_CORNERS_OF_ELEM_DIM];
345 int s;
346
347 /* sort corners of each son */
348 for (s=0; s<ER_NSONS(er); s++)
349 {
350 DCorners2Corners(ER_NCO(er,s),ER_DCO(er,s),corners);
351 std::sort(corners, corners + ER_NCO(er,s));
352 oco[s] = Corner2DCorners(ER_NCO(er,s),corners);
353 }
354
355 /* sort sons in ascending order coded on a DOUBLE (cf. `Corner2DCorners`) */
356 std::sort(oco, oco + ER_NSONS(er));
357 }
358
359 /****************************************************************************/
360 /*doctext_disabled
361 Hash_Init - allocate and initialize hash table
362
363 SYNOPSIS:
364 static INT Hash_Init (int MarkKey)
365
366 PAAMETERS:
367 . MarkKey - mark key for memory alocation
368
369 DESCRIPTION:
370 Allocate and initialize hash table
371
372 RETURN VALUE:
373 INT
374 .n 0: ok
375 .n 1: no memory
376 doctext_disabled*/
377 /****************************************************************************/
378
379 #if (defined __THREEDIM__) || (defined __DEBUG_ER__)
Hash_Init(int MarkKey)380 static INT Hash_Init (int MarkKey)
381 {
382 int i;
383
384 global.hash_table = (HRULE**) GetTmpMem(global.heap,HASH_SIZE*sizeof(HRULE*),MarkKey);
385 if (global.hash_table==NULL)
386 REP_ERR_RETURN(1);
387
388 for (i=0; i<HASH_SIZE; i++)
389 global.hash_table[i] = NULL;
390
391 return (0);
392 }
393 #endif
394
395 /****************************************************************************/
396 /*doctext_disabled
397 Hash_InsertRule - insert a rule into the hash table
398
399 SYNOPSIS:
400 static HRID Hash_InsertRule (INT etag, INT key, const ERULE *er, const DOUBLE oco[], HRULE **next_handle)
401
402 PAAMETERS:
403 . etag - element tag
404 . key - hash key
405 . er - rule
406 . oco - ordered corners
407 . next_handle - insert here into linked list of rules in hash table
408
409 DESCRIPTION:
410 Make an HRULE from an ERULE and insert it into the hash table.
411
412 RETURN VALUE:
413 HRID
414 .n >=0: ok
415 .n -1: no memory for HRULE
416 doctext_disabled*/
417 /****************************************************************************/
418
Hash_InsertRule(INT etag,INT key,const ERULE * er,const DOUBLE oco[],HRULE ** next_handle)419 static HRID Hash_InsertRule (INT etag, INT key, const ERULE *er, const DOUBLE oco[], HRULE **next_handle)
420 {
421 size_t size = sizeof(HRULE) /* full HRULE */
422 +sizeof(DOUBLE)*
423 (2*ER_NSONS(er) /* #DOUBLEs needed */
424 -MAX_SONS); /* #DOUBLEs at end of HRULE */
425 // I think that nullptr can replace the call to GetCurrentMultigrid, because that method
426 // has been returning nullptr for a long time.
427 // In other words: it seems the method we are in is never called.
428 HRULE *hr = (HRULE*) GetMemoryForObject(nullptr/*GetCurrentMultigrid()*/,size,MAOBJ);
429 HRID id = global.maxrule[etag]++;
430
431
432
433
434 if (hr==NULL)
435 REP_ERR_RETURN (-1);
436
437 /* insert in list */
438 HR_NEXT(hr) = *next_handle;
439 *next_handle = hr;
440
441 /* init */
442 HR_ID(hr) = id;
443 HR_KEY(hr) = key;
444 HR_TAG(hr) = etag;
445 memcpy(HR_ERULE(hr), er, sizeof(ERULE)+sizeof(DOUBLE)*(ER_NSONS(er)-MAX_SONS));
446 memcpy(HR_OCOPTR(hr), oco, sizeof(DOUBLE)*ER_NSONS(er));
447
448 return (id);
449 }
450
451 /****************************************************************************/
452 /*doctext_disabled
453 SonsAreEqual - compare sons of rule with array of ordered sons
454
455 SYNOPSIS:
456 static INT SonsAreEqual (INT nsons, const DOUBLE oco[], const HRULE *hr)
457
458 PAAMETERS:
459 . nsons - number of sons
460 . oco - array of ordered sons
461 . hr - rule
462
463 DESCRIPTION:
464 Compare sons of rule with array of ordered sons
465
466 RETURN VALUE:
467 INT
468 .n YES: sons are equal
469 .n NO: sons are not equal
470 doctext_disabled*/
471 /****************************************************************************/
472
SonsAreEqual(INT nsons,const DOUBLE oco[],const HRULE * hr)473 static INT SonsAreEqual (INT nsons, const DOUBLE oco[], const HRULE *hr)
474 {
475 if (nsons!=HR_NSONS(hr))
476 return (NO);
477 else
478 {
479 int s;
480 for (s=0; s<nsons; s++)
481 if (oco[s]!=HR_OCO(hr,s))
482 return (NO);
483 return (YES);
484 }
485 }
486
487 /****************************************************************************/
488 /*doctext_disabled
489 GetRuleID - insert in hash table if not found there and return rule id for element
490
491 SYNOPSIS:
492 static HRID GetRuleID (ELEMENT *elem, INT etag, const ERULE *er)
493
494 PAAMETERS:
495 . elem - element (only in debug mode)
496 . etag - element tag
497 . er - rule
498
499 DESCRIPTION:
500 Insert in hash table if not found there and return rule id for element.
501
502 RETURN VALUE:
503 HRID
504 .n rule id (starting after rm IDs per tag)
505 doctext_disabled*/
506 /****************************************************************************/
507
508 #if (defined __THREEDIM__) || (defined __DEBUG_ER__)
GetRuleID(ELEMENT * elem,INT etag,const ERULE * er)509 static HRID GetRuleID
510 (
511 #ifdef Debug
512 ELEMENT *elem,
513 #endif
514 INT etag,
515 const ERULE *er
516 )
517 {
518 DOUBLE key = 0;
519 HRULE *hr;
520 DOUBLE oco[MAX_SONS];
521 int s,h;
522
523 global.nelem_inspected[etag]++;
524 global.nelems_inspected++;
525
526 PRINTDEBUG(gm,ER_DBG_ELEM,(ELEM_INFO("GetRuleID",YES,elem,ER_NSONS(er))));
527
528 FillOrderedSons(er,oco);
529
530 for (s=0; s<ER_NSONS(er); s++)
531 key += oco[s];
532
533 /** \todo (HRR 971211): use tag also for key? */
534
535 h = HASH_ADDRESS(key);
536
537 hr = global.hash_table[h];
538 if (hr==NULL)
539 return (Hash_InsertRule(etag,key,er,oco,&(global.hash_table[h])));
540
541 for (; hr!=NULL; hr=HR_NEXT(hr))
542 if (key==HR_KEY(hr) && etag==HR_TAG(hr) && SonsAreEqual(ER_NSONS(er),oco,hr))
543 return (HR_ID(hr));
544 else if ((HR_NEXT(hr)==NULL) || (key<HR_KEY(HR_NEXT(hr))))
545 break;
546
547 return (Hash_InsertRule(etag,key,er,oco,&HR_NEXT(hr)));
548 }
549 #endif
550
551 /****************************************************************************/
552 /*doctext_disabled
553 RuleCompare - compare er-rule with existing rm-rule
554
555 SYNOPSIS:
556 static INT RuleCompare (int id, const URULE *ur, const ERULE *er)
557
558 PAAMETERS:
559 . id - element id
560 . ur - rm rule
561 . er - er rule
562
563 DESCRIPTION:
564 Compare er-rule with existing rm-rule (used in __PATCH_UG_RULES__ mode).
565
566 RETURN VALUE:
567 INT
568 .n YES: rules are equal
569 .n NO: rules are not equal
570 doctext_disabled*/
571 /****************************************************************************/
572 #ifdef __DEBUG_ER__
RuleCompare(int id,const URULE * ur,const ERULE * er)573 static INT RuleCompare (int id, const URULE *ur, const ERULE *er)
574 {
575 const int ns0 = NSONS_OF_RULE(ur);
576 const int ns1 = ER_NSONS(er);
577 int s0,s1;
578
579 if (ns0!=ns1)
580 {
581 PrintDebug("CompareRules: rules not equal: elem%d (ns0=%d != ns1=%d)\n",id,ns0,ns1);
582 return (NO);
583 }
584
585 /* compare sons */
586 for (s0=0; s0<ns0; s0++)
587 {
588 const SONDATA *son0 = SON_OF_RULE(ur,s0);
589 int nco0 = CORNERS_OF_TAG(SON_TAG(son0));
590
591 for (s1=0; s1<ns1; s1++)
592 {
593 int nco1 = ER_NCO(er,s1);
594
595 if (nco0==nco1)
596 {
597 int corners_match = true;
598 int i;
599
600 for (i=0; i<nco0; i++)
601 {
602 int co0 = SON_CORNER(son0,i);
603 SHORT corners[MAX_CORNERS_OF_ELEM_DIM];
604 int j;
605
606 DCorners2Corners(ER_NCO(er,s1),ER_DCO(er,s1),corners);
607
608 for (j=0; j<nco1; j++)
609 if (co0==corners[j])
610 break;
611 if (j>=nco1)
612 {
613 corners_match = false;
614 break;
615 }
616 }
617 if (corners_match)
618 break;
619 }
620 }
621 if (s1>=ns1)
622 {
623 PrintDebug("CompareRules: rules not equal: elem%d (son=%d of urule not in erule)\n",id,s0);
624 return (NO);
625 }
626 }
627 return (YES);
628 }
629 #endif
630 /****************************************************************************/
631 /*doctext_disabled
632 ExtractERule - extract rule from element
633
634 SYNOPSIS:
635 static INT ExtractERule (ELEMENT *elem, ERULE *er)
636
637 PAAMETERS:
638 . elem - element
639 . er - extracted rule
640
641 DESCRIPTION:
642 Extract rule from element (rule may be incomplete in parallel)
643
644 RETURN VALUE:
645 INT
646 .n 0: ok
647 .n >0: error
648 doctext_disabled*/
649 /****************************************************************************/
650
ExtractERule(ELEMENT * elem,ERULE * er)651 static INT ExtractERule (ELEMENT *elem, ERULE *er)
652 {
653 int nsons = NSONS(elem);
654 ELEMENT *sons[MAX_SONS];
655 NODE *nodes[MAX_REFINED_CORNERS_DIM];
656 int s;
657
658 if (GetNodeContext(elem,nodes))
659 REP_ERR_RETURN(1);
660
661 /* get sons with master priority */
662 if (GetAllSons(elem,sons))
663 REP_ERR_RETURN(1);
664
665 ER_NSONS(er) = 0;
666 for (s=0; s<nsons; s++)
667 {
668 ELEMENT *son = sons[s];
669 int coe = CORNERS_OF_ELEM(son);
670 SHORT corners[MAX_CORNERS_OF_ELEM_DIM];
671 int k,j;
672
673 if (EGHOST(son)) continue;
674
675 ER_NSONS(er)++;
676
677 /* corners and dcorners */
678 ER_NCO(er,s) = coe;
679 for (j=0; j<coe; j++)
680 {
681 NODE *node = CORNER(son,j);
682
683 for (k=0; k<MAX_REFINED_CORNERS_DIM; k++)
684 if (node==nodes[k])
685 break;
686 ASSERT(k<MAX_REFINED_CORNERS_DIM);
687 corners[j] = k;
688 }
689 ER_DCO(er,s) = Corner2DCorners(CORNERS_OF_ELEM(son),corners);
690 }
691
692 #ifdef __DEBUG_ER__
693 if (!(REFINE(elem)==COPY && NSONS(elem)>1))
694 {
695 URULE *ur = MARK2RULEADR(elem,REFINE(elem));
696
697 RuleCompare(ID(elem),ur,er);
698 }
699 #endif
700
701 return (0);
702 }
703
704 /****************************************************************************/
705 /*
706 functions for parallel mode only
707 */
708 /****************************************************************************/
709
710 #ifdef ModelP
711
712 /****************************************************************************/
713 /*doctext_disabled
714 CountIFElements - count interface elements having no rm-rule
715
716 SYNOPSIS:
717 static int CountIFElements (DDD_OBJ obj)
718
719 PAAMETERS:
720 . obj - element
721
722 DESCRIPTION:
723 Count interface elements having no rm-rule (masters and VH-hgosts).
724 Those elements are flagged true.
725
726 RETURN VALUE:
727 int
728 .n 0: ok
729
730 SEE ALSO:
731 ExtractInterfaceRules
732 doctext_disabled*/
733 /****************************************************************************/
734
735 #if (defined __THREEDIM__) || (defined __DEBUG_ER__)
CountIFElements(DDD::DDDContext &,DDD_OBJ obj)736 static int CountIFElements (DDD::DDDContext&, DDD_OBJ obj)
737 {
738 ELEMENT *elem = (ELEMENT*) obj;
739
740 if (HAS_NO_RULE(elem))
741 {
742 SETTHEFLAG(elem,true);
743
744 ASSERT(global.if_elems<MAX_IFID);
745
746 SETREFINE(elem,global.if_elems);
747 global.if_elems++;
748 }
749 else
750 {
751 SETTHEFLAG(elem,false);
752 }
753
754 return (0);
755 }
756
757 /****************************************************************************/
758 /*doctext_disabled
759 InitMasterRules - init rules of master elements
760
761 SYNOPSIS:
762 static int InitMasterRules (DDD_OBJ obj)
763
764 PAAMETERS:
765 . obj - element
766
767 DESCRIPTION:
768 Init rules of master elements at processor interface with sons of master.
769
770 RETURN VALUE:
771 int
772 .n 0: ok
773 .n >0: error
774
775 SEE ALSO:
776 ExtractInterfaceRules
777 doctext_disabled*/
778 /****************************************************************************/
779
InitMasterRules(DDD::DDDContext &,DDD_OBJ obj)780 static int InitMasterRules (DDD::DDDContext&, DDD_OBJ obj)
781 {
782 ELEMENT *elem = (ELEMENT*) obj;
783
784 if (THEFLAG(elem))
785 {
786 ERULE *er = global.interface_rules+REFINE(elem);
787 return (ExtractERule((ELEMENT*)obj,er));
788 }
789 else
790 return (0);
791 }
792
793 /****************************************************************************/
794 /*doctext_disabled
795 Gather_ERULE - extract rules from interface elements
796
797 SYNOPSIS:
798 static int Gather_ERULE (DDD_OBJ obj, void *data)
799
800 PAAMETERS:
801 . obj - element
802 . data - rule
803
804 DESCRIPTION:
805 Extract rules from interface elements (masters and VH-ghosts)
806
807 RETURN VALUE:
808 int
809 .n 0: ok
810 .n >0: error
811
812 SEE ALSO:
813 ExtractInterfaceRules
814 doctext_disabled*/
815 /****************************************************************************/
816
Gather_ERULE(DDD::DDDContext &,DDD_OBJ obj,void * data)817 static int Gather_ERULE (DDD::DDDContext&, DDD_OBJ obj, void *data)
818 {
819 ELEMENT *elem = (ELEMENT*) obj;
820
821 if (THEFLAG(elem))
822 return (ExtractERule(elem,(ERULE*)data));
823 else
824 return (0);
825 }
826
827 /****************************************************************************/
828 /*doctext_disabled
829 Scatter_ERULE - write completed rule for VH-ghosts
830
831 SYNOPSIS:
832 static int Scatter_ERULE (DDD_OBJ obj, void *data)
833
834 PAAMETERS:
835 . obj - element
836 . data - rule previously completed by master
837
838 DESCRIPTION:
839 Write completed rule for VH-ghosts.
840
841 RETURN VALUE:
842 int
843 .n 0: ok
844
845 SEE ALSO:
846 ExtractInterfaceRules
847 doctext_disabled*/
848 /****************************************************************************/
849
Scatter_ERULE(DDD::DDDContext &,DDD_OBJ obj,void * data)850 static int Scatter_ERULE (DDD::DDDContext&, DDD_OBJ obj, void *data)
851 {
852 ELEMENT *elem = (ELEMENT*)obj;
853
854 if (THEFLAG(elem))
855 {
856 ERULE *er = global.interface_rules+REFINE(elem);
857 memcpy(er,data,sizeof(ERULE));
858 }
859
860 return (0);
861 }
862
863 /****************************************************************************/
864 /*doctext_disabled
865 Scatter_partial_ERULE - add son-info from VH-ghosts to master rule
866
867 SYNOPSIS:
868 static int Scatter_partial_ERULE (DDD_OBJ obj, void *data)
869
870 PAAMETERS:
871 . obj - element
872 . data - rule
873
874 DESCRIPTION:
875 Add son-info from VH-ghosts to master rule.
876
877 RETURN VALUE:
878 int
879 .n 0: ok
880
881 SEE ALSO:
882 ExtractInterfaceRules
883 doctext_disabled*/
884 /****************************************************************************/
885
Scatter_partial_ERULE(DDD::DDDContext &,DDD_OBJ obj,void * data)886 static int Scatter_partial_ERULE (DDD::DDDContext&, DDD_OBJ obj, void *data)
887 {
888 ELEMENT *elem = (ELEMENT*)obj;
889
890 if (THEFLAG(elem))
891 {
892 ERULE *er_ma = global.interface_rules+REFINE(elem);
893 ERULE *er_gh = (ERULE*)data;
894 int s_ma = ER_NSONS(er_ma);
895 int s_gh;
896
897 /* put sons at and of rule */
898 for (s_gh=0; s_gh<ER_NSONS(er_gh); s_gh++)
899 {
900 ASSERT(s_ma<MAX_SONS);
901
902 ER_NCO(er_ma,s_ma) = ER_NCO(er_gh,s_gh);
903 ER_DCO(er_ma,s_ma) = ER_DCO(er_gh,s_gh);
904 s_ma++;
905 }
906 ER_NSONS(er_ma) = s_ma;
907 }
908
909 return (0);
910 }
911
912 /****************************************************************************/
913 /*doctext_disabled
914 ExtractInterfaceERule - extract rule from interface element
915
916 SYNOPSIS:
917 static int ExtractInterfaceERule (DDD_OBJ obj)
918
919 PAAMETERS:
920 . obj - element
921
922 DESCRIPTION:
923 Extract rule from interface element and get id from hash table (if there
924 is no rm-rule for this element).
925
926 RETURN VALUE:
927 int
928 .n 0: ok
929
930 SEE ALSO:
931 ExtractInterfaceRules
932 doctext_disabled*/
933 /****************************************************************************/
934
ExtractInterfaceERule(DDD::DDDContext &,DDD_OBJ obj)935 static int ExtractInterfaceERule (DDD::DDDContext&, DDD_OBJ obj)
936 {
937 ELEMENT *elem = (ELEMENT*)obj;
938
939 if (THEFLAG(elem))
940 {
941 ERULE *er = global.interface_rules+REFINE(elem);
942 int id = GetRuleID(
943 #ifdef Debug
944 elem,
945 #endif
946 TAG(elem),er);
947
948 ASSERT(id<MAX_HRID);
949
950 SETREFINE(elem,id);
951 }
952 else
953 {
954 global.nelem_not_inspected[TAG(elem)]++;
955 global.nelems_not_inspected++;
956
957 IFDEBUG(gm,ER_DBG_ELEM)
958 int nsons;
959 ELEMENT *sons[MAX_SONS];
960
961 if (GetSons(elem,sons))
962 REP_ERR_RETURN(1);
963 for (nsons=0; nsons<MAX_SONS; nsons++)
964 if (sons[nsons]==NULL)
965 break;
966 PrintDebug(ELEM_INFO("ExtractInterfaceERule",NO,elem,nsons));
967 ENDDEBUG
968 }
969
970 SETUSED(elem,true);
971
972 return (0);
973 }
974
975 /****************************************************************************/
976 /*doctext_disabled
977 ExtractInterfaceRules - extract all rules from interface elements without rm-rule
978
979 SYNOPSIS:
980 static INT ExtractInterfaceRules (MULTIGRID *mg)
981
982 PAAMETERS:
983 . mg - multigrid
984
985 DESCRIPTION:
986 Extract all rules from interface elements without rm-rule.
987
988 RETURN VALUE:
989 INT
990 .n 0: ok
991 .n >0: error in storage allocation
992 doctext_disabled*/
993 /****************************************************************************/
994
ExtractInterfaceRules(MULTIGRID * mg)995 static INT ExtractInterfaceRules (MULTIGRID *mg)
996 {
997 auto& context = mg->dddContext();
998 const auto& dddctrl = ddd_ctrl(context);
999
1000 int lev;
1001
1002 /* TODO (HRR 971211): don't include TOPLEVEL (no elem refined there) */
1003 for (lev=0; lev<=TOPLEVEL(mg); lev++)
1004 {
1005 GRID *grid = GRID_ON_LEVEL(mg,lev);
1006
1007 /* count interface master and vhghost elements */
1008 global.if_elems = 1;
1009 DDD_IFAExecLocal(context, dddctrl.ElementVHIF, GRID_ATTR(grid), CountIFElements);
1010
1011 if (global.if_elems>1)
1012 {
1013 INT MarkKey;
1014
1015 PRINTDEBUG(gm,ER_DBG_GENERAL,("ExtractInterfaceRules: interface elements to consider on level %d: %ld\n",lev,global.if_elems));
1016
1017 if (MarkTmpMem(global.heap,&MarkKey))
1018 REP_ERR_RETURN(1);
1019
1020 /* REMARK (HRR 971207): storage could be Marked and Released in pieces if necessary
1021 to allow BOTTOM-storage to grow */
1022 global.interface_rules = (ERULE*)GetTmpMem(global.heap,(global.if_elems+1)*sizeof(ERULE),MarkKey);
1023 if (global.interface_rules==NULL)
1024 REP_ERR_RETURN(1);
1025
1026 /* init rules of masters */
1027 DDD_IFAExecLocal(context, dddctrl.ElementIF, GRID_ATTR(grid), InitMasterRules);
1028
1029 /* communicate VHghosts --> master */
1030 DDD_IFAOneway(context, dddctrl.ElementVHIF, GRID_ATTR(grid), IF_BACKWARD, sizeof(ERULE),
1031 Gather_ERULE, Scatter_partial_ERULE);
1032
1033 /* communicate master --> VHghosts */
1034 DDD_IFAOneway(context, dddctrl.ElementVHIF, GRID_ATTR(grid), IF_FORWARD, sizeof(ERULE),
1035 Gather_ERULE, Scatter_ERULE);
1036
1037 /* extract rules from interface elements */
1038 DDD_IFAExecLocal(context, dddctrl.ElementVHIF, GRID_ATTR(grid), ExtractInterfaceERule);
1039
1040 IFDEBUG(gm,ER_DBG_GENERAL)
1041 long N_er = 0;
1042 int tag;
1043
1044 PrintDebug("ExtractInterfaceRules: number of refrules extracted from interface on level %d:\n",lev);
1045 for (tag=0; tag<TAGS; tag++)
1046 {
1047 long n_er = global.maxrule[tag] - UGMAXRULE(tag);
1048 N_er += n_er;
1049 PrintDebug("tag %d: %ld extracted rules\n",tag,n_er);
1050 }
1051 PrintDebug("total: %ld extracted rules\n",N_er);
1052 ENDDEBUG
1053
1054 if (ReleaseTmpMem(global.heap,MarkKey))
1055 REP_ERR_RETURN(1);
1056 }
1057 }
1058
1059 return (0);
1060 }
1061 #endif
1062 #endif /* ModelP */
1063
1064 /****************************************************************************/
1065 /*doctext_disabled
1066 ExtractRules - extract all realized rules that have no rm-rule
1067
1068 SYNOPSIS:
1069 static INT ExtractRules (MULTIGRID *mg)
1070
1071 PAAMETERS:
1072 . mg - multigrid
1073
1074 DESCRIPTION:
1075 Extract all realized rules of the mg that have no rm-rule.
1076
1077 RETURN VALUE:
1078 INT
1079 .n 0: ok
1080 .n >0: error
1081 doctext_disabled*/
1082 /****************************************************************************/
1083
1084 #if (defined __THREEDIM__) || (defined __DEBUG_ER__)
ExtractRules(MULTIGRID * mg)1085 static INT ExtractRules (MULTIGRID *mg)
1086 {
1087 ELEMENT *elem;
1088 ERULE er;
1089 HRID id;
1090 INT MarkKey;
1091 int l,tag,h,maxrules;
1092
1093 /* for hash table */
1094 if (MarkTmpMem(global.heap,&MarkKey))
1095 REP_ERR_RETURN(1);
1096
1097 if (Hash_Init(MarkKey))
1098 REP_ERR_RETURN(1);
1099
1100 global.nelems_inspected = global.nelems_not_inspected = 0;
1101 for (tag=0; tag<TAGS; tag++)
1102 global.nelem_inspected[tag] = global.nelem_not_inspected[tag] = 0;
1103
1104 /* loop elements and extract rules */
1105 #ifdef ModelP
1106 /* TODO (HRR 971211): don't include TOPLEVEL (no elem refined there) */
1107 for (l=0; l<=TOPLEVEL(mg); l++)
1108 for (elem=FIRSTELEMENT(GRID_ON_LEVEL(mg,l)); elem!=NULL; elem=SUCCE(elem))
1109 SETUSED(elem,false);
1110 if (ExtractInterfaceRules(mg))
1111 REP_ERR_RETURN(1);
1112 #endif
1113
1114 /* TODO (HRR 971211): don't include TOPLEVEL (no elem refined there) */
1115 for (l=0; l<=TOPLEVEL(mg); l++)
1116 for (elem=FIRSTELEMENT(GRID_ON_LEVEL(mg,l)); elem!=NULL; elem=SUCCE(elem))
1117 #ifdef ModelP
1118 if (!USED(elem))
1119 {
1120 #endif
1121 if (HAS_NO_RULE(elem))
1122 {
1123 if (ExtractERule(elem,&er))
1124 REP_ERR_RETURN(1);
1125 id = GetRuleID(
1126 #ifdef Debug
1127 elem,
1128 #endif
1129 TAG(elem),&er);
1130 ASSERT(id<MAX_HRID);
1131 SETREFINE(elem,id);
1132 }
1133 else
1134 {
1135 global.nelem_not_inspected[TAG(elem)]++;
1136 global.nelems_not_inspected++;
1137 IFDEBUG(gm,ER_DBG_ELEM)
1138 int nsons;
1139 ELEMENT *sons[MAX_SONS];
1140
1141 if (GetSons(elem,sons))
1142 REP_ERR_RETURN(1);
1143 for (nsons=0; nsons<MAX_SONS; nsons++)
1144 if (sons[nsons]==NULL)
1145 break;
1146 PrintDebug(ELEM_INFO("ExtractRules",NO,elem,nsons));
1147 if (NSONS(elem)!=nsons)
1148 PrintDebug("------ ERROR: NSONS!=nsons -------\n");
1149 ENDDEBUG
1150 }
1151 #ifdef ModelP
1152 }
1153 #endif
1154
1155 global.maxrules = maxrules = 0;
1156 for (tag=0; tag<TAGS; tag++)
1157 {
1158 global.maxrules += global.maxrule[tag];
1159 maxrules += global.maxrule[tag] - UGMAXRULE(tag);
1160 }
1161
1162 if (maxrules>0)
1163 {
1164 int n = 0;
1165 int max_list_len = 0;
1166
1167 /* make tables of subsequent IDs */
1168 // I think that nullptr can replace the call to GetCurrentMultigrid, because that method
1169 // has been returning nullptr for a long time.
1170 // In other words: it seems the method we are in is never called.
1171 global.hrule[0] = (HRULE**)
1172 GetMemoryForObject(nullptr/*GetCurrentMultigrid()*/,
1173 global.maxrules*sizeof(HRULE*),MAOBJ);
1174 if (global.hrule[0]==NULL)
1175 REP_ERR_RETURN(1);
1176 for (tag=1; tag<TAGS; tag++)
1177 global.hrule[tag] = global.hrule[tag-1]+global.maxrule[tag-1];
1178
1179 for (h=0; h<HASH_SIZE; h++)
1180 {
1181 int list_len = 0;
1182 HRULE *hr;
1183
1184 for (hr=global.hash_table[h]; hr!=NULL; hr=HR_NEXT(hr))
1185 {
1186 global.hrule[HR_TAG(hr)][HR_ID(hr)] = hr;
1187 list_len++;
1188 n++;
1189 }
1190 max_list_len = MAX(max_list_len,list_len);
1191 }
1192 ASSERT(maxrules==n);
1193
1194 PRINTDEBUG(gm,ER_DBG_GENERAL,("ExtractRules: max members of hash table entries is %d\n",max_list_len));
1195 }
1196 else
1197 PRINTDEBUG(gm,ER_DBG_GENERAL,("ExtractRules: realized rules were completely covered by rm rules\n"));
1198
1199 if (ReleaseTmpMem(global.heap,MarkKey))
1200 REP_ERR_RETURN(1);
1201
1202 return (0);
1203 }
1204 #endif
1205
1206 /****************************************************************************/
1207 /*doctext_disabled
1208 FindPathForNeighbours - recursive function computing path infos
1209
1210 SYNOPSIS:
1211 static void FindPathForNeighbours (MGIO_RR_RULE *rule, SHORT myID, SHORT Status[MAX_SONS])
1212
1213 PAAMETERS:
1214 . rule - mrule
1215 . myID - son ID
1216 . Status - status (see 'enum NB_STATUS' in er.c)
1217
1218 DESCRIPTION:
1219 Recursive function computing path infos for sons.
1220
1221 RETURN VALUE:
1222 void
1223 .n none
1224 doctext_disabled*/
1225 /****************************************************************************/
1226
FindPathForNeighbours(MGIO_RR_RULE * rule,SHORT myID,SHORT Status[MAX_SONS])1227 static void FindPathForNeighbours (MGIO_RR_RULE *rule, SHORT myID, SHORT Status[MAX_SONS])
1228 {
1229 SHORT i,nbID;
1230
1231 for (i=0; i<MAX_SIDES_OF_ELEM; i++)
1232 if (((nbID=rule->sons[myID].nb[i])<FATHER_SIDE_OFFSET) && (Status[nbID]==NB_NOTDONE))
1233 {
1234 INT nbPath;
1235 SHORT nbPathDepth;
1236
1237 /* copy myPath to nbPath */
1238 nbPath = rule->sons[myID].path;
1239
1240 /* complete nbPath */
1241 nbPathDepth = PATHDEPTH(nbPath);
1242 SETNEXTSIDE(nbPath,nbPathDepth,i);
1243 SETPATHDEPTH(nbPath,++nbPathDepth);
1244 Status[nbID] = NB_TOUCHED;
1245 }
1246
1247 /* recursive call for NB_TOUCHED sons */
1248 for (nbID=1; nbID<rule->nsons; nbID++)
1249 if (Status[nbID]==NB_TOUCHED)
1250 {
1251 Status[nbID] = NB_DONE;
1252 FindPathForNeighbours(rule,nbID,Status);
1253 }
1254
1255 return;
1256 }
1257
1258 /****************************************************************************/
1259 /*doctext_disabled
1260 FillSonPaths - fill son path components in mrule
1261
1262 SYNOPSIS:
1263 static void FillSonPaths (MGIO_RR_RULE *rule)
1264
1265 PAAMETERS:
1266 . rule - mrule
1267
1268 DESCRIPTION:
1269 Fill son path components in mrule by recursion (using 'FindPathForNeighbours').
1270
1271 RETURN VALUE:
1272 void
1273 .n none
1274 doctext_disabled*/
1275 /****************************************************************************/
1276 #ifdef __THREEDIM__
FillSonPaths(MGIO_RR_RULE * rule)1277 static void FillSonPaths (MGIO_RR_RULE *rule)
1278 {
1279 SHORT Status[MAX_SONS];
1280 int i;
1281
1282 /* TODO (HRR 971211): debug recursive path construction
1283 (taken from GenerateRules but not used by ugio) */
1284
1285 /* son 0 has trivial path */
1286 Status[0] = NB_DONE;
1287 rule->sons[0].path = 0;
1288 for (i=1; i<rule->nsons; i++)
1289 Status[i] = NB_NOTDONE;
1290
1291 /* start recursion with son 0 */
1292 FindPathForNeighbours(rule,0,Status);
1293
1294 return;
1295 }
1296 #endif
1297 /****************************************************************************/
1298 /*doctext_disabled
1299 GetFSidesOfCorners - fill array for father side the given side corners belong to
1300
1301 SYNOPSIS:
1302 static INT GetFSidesOfCorners (int tag, int n, SHORT corners[MAX_CORNERS_OF_SIDE], SHORT corner_on_side[MAX_CORNERS_OF_SIDE][MAX_SIDES_OF_ELEM])
1303
1304 PAAMETERS:
1305 . tag - father tag
1306 . n - number of side nodes
1307 . corners - side node local IDs
1308 . corner_on_side - YES/NO indicating corner lying on father side
1309
1310 DESCRIPTION:
1311 Fill array for father side the given side corners belong to.
1312
1313 RETURN VALUE:
1314 INT
1315 .n YES: ok, maybe there is a common father side
1316 .n NO: sorry, center node definitely not on a father side
1317 doctext_disabled*/
1318 /****************************************************************************/
1319
GetFSidesOfCorners(int tag,int n,SHORT corners[MAX_CORNERS_OF_SIDE],SHORT corner_on_side[MAX_CORNERS_OF_SIDE][MAX_SIDES_OF_ELEM])1320 static INT GetFSidesOfCorners (int tag, int n, SHORT corners[MAX_CORNERS_OF_SIDE], SHORT corner_on_side[MAX_CORNERS_OF_SIDE][MAX_SIDES_OF_ELEM])
1321 {
1322 int coe = CORNERS_OF_TAG(tag);
1323 int eoe = EDGES_OF_TAG(tag);
1324 int soe = SIDES_OF_TAG(tag);
1325 int co,side;
1326
1327 for (co=0; co<n; co++)
1328 for (side=0; side<soe; side++)
1329 corner_on_side[co][side] = false;
1330
1331 for (co=0; co<n; co++)
1332 if (corners[co]==coe+CENTER_NODE_INDEX_TAG(tag))
1333 {
1334 /* center node: can not be part of a father side */
1335 return (NO);
1336 }
1337 else if (corners[co]<coe)
1338 {
1339 /* father corner */
1340 int fco = corners[co];
1341
1342 for (side=0; side<soe; side++)
1343 if (CORNER_OF_SIDE_INV_TAG(tag,side,fco)>=0)
1344 corner_on_side[co][side] = true;
1345 }
1346 else if (corners[co]<(coe+eoe))
1347 {
1348 /* edge mid node */
1349 int ed = corners[co]-coe;
1350
1351 #ifdef __TWODIM__
1352 corner_on_side[co][ed] = true;
1353 #else
1354 int i;
1355 for (i=0; i<MAX_SIDES_OF_EDGE; i++)
1356 {
1357 int sd = SIDE_WITH_EDGE_TAG(tag,ed,i);
1358 if (sd>=0)
1359 corner_on_side[co][sd] = true;
1360 }
1361 #endif
1362 }
1363 #ifdef __THREEDIM__
1364 else if (corners[co]<(coe+eoe+soe))
1365 {
1366 /* side mid */
1367 int sd = corners[co]-(coe+eoe);
1368 corner_on_side[co][sd] = true;
1369 }
1370 #endif
1371 else
1372 ASSERT(false); /* Huh??? */
1373
1374 return (YES);
1375 }
1376
1377 /****************************************************************************/
1378 /*doctext_disabled
1379 GetCommonFSide - return common father element side iff
1380
1381 SYNOPSIS:
1382 static INT GetCommonFSide (int nco, int nsi, SHORT corner_on_side[MAX_CORNERS_OF_SIDE][MAX_SIDES_OF_ELEM])
1383
1384 PAAMETERS:
1385 . nco - number of corners on side
1386 . nsi - number of sides
1387 . corner_on_side - array filled by 'GetFSidesOfCorners'
1388
1389 DESCRIPTION:
1390 Return common father element side iff.
1391
1392 RETURN VALUE:
1393 INT
1394 .n father side if common for all side corners
1395 .n nsi else
1396 doctext_disabled*/
1397 /****************************************************************************/
1398
GetCommonFSide(int nco,int nsi,SHORT corner_on_side[MAX_CORNERS_OF_SIDE][MAX_SIDES_OF_ELEM])1399 static INT GetCommonFSide (int nco, int nsi, SHORT corner_on_side[MAX_CORNERS_OF_SIDE][MAX_SIDES_OF_ELEM])
1400 {
1401 int i,side;
1402
1403 for (side=0; side<nsi; side++)
1404 {
1405 for (i=0; i<nco; i++)
1406 if (!corner_on_side[i][side])
1407 break;
1408 if (i>=nco)
1409 return (side);
1410 }
1411 return (nsi);
1412 }
1413
1414 /****************************************************************************/
1415 /*doctext_disabled
1416 IsOnFatherSide - determine whether son side lies on father side
1417
1418 SYNOPSIS:
1419 static INT IsOnFatherSide (int tag, int nsco, SHORT sco[], SHORT *nb)
1420
1421 PAAMETERS:
1422 . tag - father tag
1423 . nsco - number of side corners
1424 . sco - side corner array
1425 . nb - fill nb component with FATHER_SIDE_OFFSET+fside iff common father side
1426
1427 DESCRIPTION:
1428 Determine whether son side lies on father side.
1429
1430 RETURN VALUE:
1431 INT
1432 .n YES: there is a common father side
1433 .n NO: there is no common father side
1434 doctext_disabled*/
1435 /****************************************************************************/
1436
IsOnFatherSide(int tag,int nsco,SHORT sco[],SHORT * nb)1437 static INT IsOnFatherSide (int tag, int nsco, SHORT sco[], SHORT *nb)
1438 {
1439 SHORT sco_on_side[MAX_CORNERS_OF_SIDE][MAX_SIDES_OF_ELEM];
1440
1441 if (GetFSidesOfCorners(tag,nsco,sco,sco_on_side))
1442 {
1443 int fside = GetCommonFSide(nsco,SIDES_OF_TAG(tag),sco_on_side);
1444
1445 if (fside<SIDES_OF_TAG(tag))
1446 {
1447 *nb = FATHER_SIDE_OFFSET+fside;
1448 return (YES);
1449 }
1450 }
1451 return (NO);
1452 }
1453
1454 /****************************************************************************/
1455 /*doctext_disabled
1456 SidesMatch - check whether two son sides match
1457
1458 SYNOPSIS:
1459 static INT SidesMatch (int nsco, SHORT sco0[], SHORT sco1[])
1460
1461 PAAMETERS:
1462 . nsco - number of side corners
1463 . sco0 - side corners of first son
1464 . sco1 - side corners of second son
1465
1466 DESCRIPTION:
1467 Check whether two son sides match.
1468
1469 RETURN VALUE:
1470 INT
1471 .n YES: son sides are mathing
1472 .n NO: son sides are not mathing
1473 doctext_disabled*/
1474 /****************************************************************************/
1475
SidesMatch(int nsco,SHORT sco0[],SHORT sco1[])1476 static INT SidesMatch (int nsco, SHORT sco0[], SHORT sco1[])
1477 {
1478 int i;
1479
1480 /* try each permutation of first with reverse order of second */
1481 for (i=0; i<nsco; i++)
1482 {
1483 int match = true;
1484 int j;
1485
1486 for (j=0; j<nsco; j++)
1487 if (sco0[(i+j)%nsco]!=sco1[nsco-j-1])
1488 {
1489 match = false;
1490 break;
1491 }
1492 if (match)
1493 return (YES);
1494 }
1495 return (NO);
1496 }
1497
1498 /****************************************************************************/
1499 /*doctext_disabled
1500 HRule2Mrule - convert rule to MGIO_RR_RULE
1501
1502 SYNOPSIS:
1503 static void HRule2Mrule (const HRULE *hr, MGIO_RR_RULE *mr)
1504
1505 PAAMETERS:
1506 . hr - hashed rule
1507 . mr - mgio rule
1508
1509 DESCRIPTION:
1510 Convert local rule type HRULE (hashed rule) to MGIO_RR_RULE.
1511
1512 RETURN VALUE:
1513 void
1514 .n
1515 doctext_disabled*/
1516 /****************************************************************************/
1517
HRule2Mrule(const HRULE * hr,MGIO_RR_RULE * mr)1518 static void HRule2Mrule (const HRULE *hr, MGIO_RR_RULE *mr)
1519 {
1520 int coe = CORNERS_OF_TAG(HR_TAG(hr));
1521 int s,s0;
1522
1523 /* extracted rules are always irregular */
1524 /* TODO (HRR 971208): is that really ok (actually not used)? */
1525 mr->rclass = GREEN_CLASS;
1526 mr->nsons = HR_NSONS(hr);
1527
1528 {int i; for (i=0; i<MGIO_MAX_NEW_CORNERS; i++)
1529 mr->pattern[i] = 0;}
1530
1531 /* son tags, corners, sonandnode, pattern */
1532 for (s=0; s<mr->nsons; s++)
1533 {
1534 struct mgio_sondata *sonData = &(mr->sons[s]);
1535 int nco = HR_NCO(hr,s);
1536 int j;
1537 int sd0;
1538
1539 for (sd0=0; sd0<MAX_SIDES_OF_ELEM; sd0++)
1540 sonData->nb[sd0] = NOTDONE;
1541
1542 sonData->tag = REF2TAG(nco);
1543 DCorners2Corners(nco,HR_DCO(hr,s),sonData->corners);
1544 for (j=0; j<nco; j++)
1545 {
1546 int newco = sonData->corners[j]-coe;
1547
1548 if (newco<0)
1549 continue;
1550
1551 /* TODO (HRR 971208): what about center node pattern for tets? */
1552 mr->pattern[newco] = 1;
1553 mr->sonandnode[newco][0] = s;
1554 mr->sonandnode[newco][1] = j;
1555 }
1556 }
1557
1558 /* son nb */
1559 for (s0=0; s0<mr->nsons; s0++)
1560 {
1561 struct mgio_sondata *sonData0 = mr->sons+s0;
1562 int sd0;
1563
1564 for (sd0=0; sd0<MAX_SIDES_OF_ELEM; sd0++)
1565 if (sonData0->nb[sd0] == NOTDONE)
1566 {
1567 int stag0 = sonData0->tag;
1568 int nsco0 = CORNERS_OF_SIDE_TAG(stag0,sd0);
1569 SHORT sco0[MAX_CORNERS_OF_SIDE];
1570
1571 {int j; for (j=0; j<nsco0; j++)
1572 sco0[j] = sonData0->corners[CORNER_OF_SIDE_TAG(stag0,sd0,j)];}
1573
1574 if (!IsOnFatherSide(HR_TAG(hr),nsco0,sco0,&(sonData0->nb[sd0])))
1575 {
1576 int found = NO;
1577 int s1;
1578
1579 /* find matching side of other son */
1580 for (s1=s0+1; s1<mr->nsons; s1++)
1581 {
1582 struct mgio_sondata *sonData1 = mr->sons+s1;
1583 int sd1;
1584
1585 for (sd1=0; sd1<MAX_SIDES_OF_ELEM; sd1++)
1586 {
1587 int stag1 = sonData1->tag;
1588 int nsco1 = CORNERS_OF_SIDE_TAG(stag1,sd1);
1589
1590 if (nsco0==nsco1)
1591 {
1592 SHORT sco1[MAX_CORNERS_OF_SIDE];
1593
1594 {int j; for (j=0; j<nsco1; j++)
1595 sco1[j] = sonData1->corners[CORNER_OF_SIDE_TAG(stag1,sd1,j)];}
1596
1597 if (SidesMatch(nsco0,sco0,sco1))
1598 {
1599 /* neighbour found */
1600 sonData0->nb[sd0] = s1;
1601 sonData1->nb[sd1] = s0;
1602
1603 found = YES;
1604 break;
1605 }
1606 }
1607 }
1608 if (found)
1609 break;
1610 }
1611 ASSERT(found);
1612 }
1613 }
1614 }
1615
1616 #ifdef __THREEDIM__
1617 /* 2D: not even in rm */
1618 /* son path */
1619 FillSonPaths(mr);
1620 #endif
1621
1622 return;
1623 }
1624
1625 /****************************************************************************/
1626 /*doctext_disabled
1627 URule2Mrule - convert rule manager rule to MGIO_RR_RULE
1628
1629 SYNOPSIS:
1630 static void URule2Mrule (const URULE *ur, MGIO_RR_RULE *mr)
1631
1632 PAAMETERS:
1633 . ur - rule manager rule
1634 . mr - mgio rule
1635
1636 DESCRIPTION:
1637 Convert rule manager rule to mgio rule.
1638
1639 RETURN VALUE:
1640 void
1641 .n
1642 doctext_disabled*/
1643 /****************************************************************************/
1644
URule2Mrule(const URULE * ur,MGIO_RR_RULE * mr)1645 static void URule2Mrule (const URULE *ur, MGIO_RR_RULE *mr)
1646 {
1647 int j,k;
1648
1649 mr->rclass = ur->rclass;
1650 mr->nsons = ur->nsons;
1651 for (j=0; j<MGIO_MAX_NEW_CORNERS; j++)
1652 mr->pattern[j] = ur->pattern[j];
1653 for (j=0; j<MGIO_MAX_NEW_CORNERS; j++)
1654 {
1655 mr->sonandnode[j][0] = ur->sonandnode[j][0];
1656 mr->sonandnode[j][1] = ur->sonandnode[j][1];
1657 }
1658 for (j=0; j<mr->nsons; j++)
1659 {
1660 struct mgio_sondata *sonData = &(mr->sons[j]);
1661
1662 sonData->tag = ur->sons[j].tag;
1663 for (k=0; k<MGIO_MAX_CORNERS_OF_ELEM; k++)
1664 sonData->corners[k] = ur->sons[j].corners[k];
1665 for (k=0; k<MGIO_MAX_SIDES_OF_ELEM; k++)
1666 sonData->nb[k] = ur->sons[j].nb[k];
1667 sonData->path = ur->sons[j].path;
1668 }
1669 return;
1670 }
1671
WriteDebugInfo(void)1672 static void WriteDebugInfo (void)
1673 {
1674 long N_rm=0,N_er=0;
1675 int tag;
1676
1677 /* number of rules (rm+er) */
1678 PrintDebug("------------- Write_RefRules statistics --------------\n");
1679 PrintDebug("number of refrules:\n");
1680 for (tag=0; tag<TAGS; tag++)
1681 {
1682 long n_rm = UGMAXRULE(tag);
1683 long n_er = global.maxrule[tag]-UGMAXRULE(tag);
1684 N_rm += n_rm;
1685 N_er += n_er;
1686 PrintDebug("tag %d: %3ld rm rules, %4ld extracted rules (elems inspected: %6ld yes %6ld no)\n",
1687 tag,
1688 n_rm,
1689 n_er,
1690 global.nelem_inspected[tag],
1691 global.nelem_not_inspected[tag]);
1692 }
1693 PrintDebug("total: %3ld rm rules, %4ld extracted rules (elems inspected: %6ld yes %6ld no)\n",
1694 N_rm,
1695 N_er,
1696 global.nelem_inspected[tag],
1697 global.nelem_not_inspected[tag]);
1698
1699 PrintDebug("------------------------------------------------------\n");
1700
1701 return;
1702 }
1703
GetOrderedSons(ELEMENT * theElement,MGIO_RR_RULE * theRule,NODE ** NodeContext,ELEMENT ** SonList,INT * nmax)1704 INT NS_DIM_PREFIX GetOrderedSons (ELEMENT *theElement, MGIO_RR_RULE *theRule, NODE **NodeContext, ELEMENT **SonList, INT *nmax)
1705 {
1706 INT i,j,k,l,found;
1707 ELEMENT *NonorderedSonList[MAX_SONS];
1708 NODE *theNode;
1709
1710 *nmax = 0;
1711
1712 if (GetAllSons(theElement,NonorderedSonList)) REP_ERR_RETURN(1);
1713 for (i=0; i<theRule->nsons; i++)
1714 {
1715 found=1;
1716 for (j=0; j<CORNERS_OF_TAG(theRule->sons[i].tag); j++)
1717 if (NodeContext[theRule->sons[i].corners[j]] == NULL)
1718 {
1719 found=0;
1720 break;
1721 }
1722 if (!found)
1723 {
1724 SonList[i] = NULL;
1725 continue;
1726 }
1727
1728 /* identify (hopefully) an element of SonList */
1729 for (j=0; NonorderedSonList[j]!=NULL; j++)
1730 {
1731 found=0;
1732 for (l=0; l<CORNERS_OF_TAG(theRule->sons[i].tag); l++)
1733 {
1734 theNode = NodeContext[theRule->sons[i].corners[l]];
1735 for (k=0; k<CORNERS_OF_ELEM(NonorderedSonList[j]); k++)
1736 if (theNode==CORNER(NonorderedSonList[j],k))
1737 {
1738 found++;
1739 break;
1740 }
1741 }
1742 if (found==CORNERS_OF_TAG(theRule->sons[i].tag))
1743 {
1744 SonList[i] = NonorderedSonList[j];
1745 *nmax = i+1;
1746 break;
1747 }
1748 else
1749 SonList[i] = NULL;
1750 }
1751 }
1752
1753 return (0);
1754 }
1755
CheckNBrelations(MGIO_RR_RULE * mr,int i,int tag)1756 static int CheckNBrelations (MGIO_RR_RULE *mr, int i, int tag)
1757 {
1758 int n = mr->nsons;
1759 int err=0;
1760 int s;
1761
1762 for (s=0; s<n; s++)
1763 {
1764 struct mgio_sondata *son = &(mr->sons[s]);
1765 int ns = SIDES_OF_TAG(son->tag);
1766 int i;
1767
1768 for (i=0; i<ns; i++)
1769 {
1770 int nb = son->nb[i];
1771 if (nb<FATHER_SIDE_OFFSET)
1772 {
1773 struct mgio_sondata *nbson = &(mr->sons[nb]);
1774 int nnbs = SIDES_OF_TAG(nbson->tag);
1775 int j;
1776
1777 for (j=0; j<nnbs; j++)
1778 if (nbson->nb[j]==s)
1779 break;
1780 if (j>=nnbs)
1781 {
1782 PrintDebug("ERROR: rule %d of %d, asym in son %d, nb %d\n",i,tag,s,i);
1783 err++;
1784 }
1785 }
1786 }
1787 }
1788 return (err);
1789 }
1790
CheckMRules(MULTIGRID * mg,INT RefRuleOffset[],MGIO_RR_RULE * mrules)1791 static void CheckMRules (MULTIGRID *mg, INT RefRuleOffset[], MGIO_RR_RULE *mrules)
1792 {
1793 ELEMENT *elem;
1794 int l;
1795 int max_path_depth = 0;
1796 int use_bug_in_rule = true;
1797 short *bug_in_rule[TAGS];
1798 INT MarkKey;
1799 int tg;
1800
1801 if (MarkTmpMem(global.heap,&MarkKey))
1802 use_bug_in_rule = false;
1803 else
1804 {
1805 bug_in_rule[0] = (short*) GetTmpMem(global.heap,global.maxrules*sizeof(short),MarkKey);
1806 if (bug_in_rule[0]==NULL)
1807 use_bug_in_rule = false;
1808 else
1809 {
1810 int t,i;
1811 for (t=0; t<TAGS; t++)
1812 {
1813 bug_in_rule[t] = bug_in_rule[0] + RefRuleOffset[t];
1814 for (i=0; i<global.maxrule[t]; i++)
1815 bug_in_rule[t][i] = false;
1816 }
1817 }
1818 }
1819
1820 /* check symmetry of nb relations */
1821 for (tg=0; tg<TAGS; tg++)
1822 {
1823 int i;
1824 for (i=0; i<global.maxrule[tg]; i++)
1825 if (CheckNBrelations(mrules+RefRuleOffset[tg]+i,i,tg))
1826 if (use_bug_in_rule)
1827 bug_in_rule[tg][i] = true;
1828 }
1829
1830 for (l=0; l<TOPLEVEL(mg); l++)
1831 for (elem=PFIRSTELEMENT(GRID_ON_LEVEL(mg,l)); elem!=NULL; elem=SUCCE(elem))
1832 if (NSONS(elem)>0)
1833 {
1834 int refi = REFINE(elem);
1835 int tag = TAG(elem);
1836 MGIO_RR_RULE *mr = mrules + RefRuleOffset[tag] + refi;
1837 long id = ID(elem);
1838 int nsons = mr->nsons;
1839 int coe = CORNERS_OF_TAG(tag);
1840 int soe = SIDES_OF_TAG(tag);
1841 NODE *nodes[MAX_REFINED_CORNERS_DIM];
1842 NODE **newnodes = nodes+coe;
1843 INT maxsonex = 0;
1844 ELEMENT *sons[MAX_SONS];
1845 int error = 0;
1846 int s,i;
1847
1848 /* check nsons */
1849 if (NSONS(elem)!=nsons)
1850 PD_ERR(ER_DBG_RULE_VERBOSE,("ERROR in CheckMRules, elem %ld: NSONS(elem)!=mr->nsons\n",id),error);
1851
1852 if (GetNodeContext(elem,nodes))
1853 ASSERT(false);
1854
1855 /* check pattern */
1856 for (i=0; i<MAX_NEW_CORNERS_DIM; i++)
1857 if (mr->pattern[i])
1858 {
1859 #ifndef ModelP
1860 if (newnodes[i]==NULL)
1861 PD_ERR(ER_DBG_RULE_VERBOSE,("ERROR in CheckMRules, elem %ld: pattern %d inconsistent (says ex but doesn't)\n",id,i),error);
1862 #endif
1863 }
1864 else
1865 {
1866 if (newnodes[i]!=NULL)
1867 PD_ERR(ER_DBG_RULE_VERBOSE,("ERROR in CheckMRules, elem %ld: pattern %d inconsistent (says nonex but does)\n",id,i),error);
1868 }
1869
1870 /* check sons */
1871 if (GetOrderedSons(elem,mr,nodes,sons,&maxsonex))
1872 ASSERT(false);
1873 #ifndef ModelP
1874 if (maxsonex!=nsons)
1875 PD_ERR(ER_DBG_RULE_VERBOSE,("ERROR in CheckMRules, elem %ld: wrong number of sons (%d vs %d)\n",id,maxsonex,nsons),error);
1876 #endif
1877 for (s=0; s<maxsonex; s++)
1878 {
1879 ELEMENT *son = sons[s];
1880 if (son!=NULL)
1881 {
1882 struct mgio_sondata *rson = &(mr->sons[s]);
1883 int nco = CORNERS_OF_ELEM(son);
1884 int nsi = SIDES_OF_ELEM(son);
1885 int co,si;
1886
1887 if (rson->tag!=TAG(son))
1888 PD_ERR(ER_DBG_RULE_VERBOSE,("ERROR in CheckMRules, elem %ld: wrong tag of son %d (%d vs %d)\n",id,s,rson->tag,TAG(son)),error);
1889
1890 /* check corners */
1891 for (co=0; co<nco; co++)
1892 if (CORNER(son,co)!=nodes[rson->corners[co]])
1893 PD_ERR(ER_DBG_RULE_VERBOSE,("ERROR in CheckMRules, elem %ld: corner %d of son %d inconsistent\n",id,co,s),error);
1894
1895 /* check neighbours */
1896 for (si=0; si<nsi; si++)
1897 {
1898 ELEMENT *nb = NBELEM(son,si);
1899 if (nb!=NULL)
1900 {
1901 ELEMENT *nbf = EFATHER(nb);
1902
1903 if (nbf==elem)
1904 {
1905 /* check inner side */
1906 #ifdef ModelP
1907 if (sons[rson->nb[si]]!=NULL)
1908 #endif
1909 if (nb!=sons[rson->nb[si]])
1910 PD_ERR(ER_DBG_RULE_VERBOSE,("ERROR in CheckMRules, elem %ld: inner nb %d of son %d inconsistent\n",id,si,s),error);
1911 }
1912 else if (nbf!=NULL)
1913 {
1914 int j;
1915
1916 /* check father side */
1917 for (j=0; j<soe; j++)
1918 if (NBELEM(elem,j)==nbf)
1919 break;
1920 if (j!=(rson->nb[si]-FATHER_SIDE_OFFSET))
1921 PD_ERR(ER_DBG_RULE_VERBOSE,("ERROR in CheckMRules, elem %ld: outer nb %d of son %d inconsistent\n",id,si,s),error);
1922 }
1923 }
1924 else
1925 {
1926 #ifndef ModelP
1927 if (rson->nb[si]<FATHER_SIDE_OFFSET)
1928 PD_ERR(ER_DBG_RULE_VERBOSE,("ERROR in CheckMRules, elem %ld: son %d has no nb but %d no fatherside\n",id,s,si),error);
1929 #endif
1930 }
1931 }
1932 #ifndef ModelP
1933 #ifdef __THREEDIM__
1934 /* check path
1935 NOT CONSISTENT in rm since not used (says Stefan, 971219)
1936 if (s>0)
1937 {
1938 ELEMENT *nson = sons[0];
1939 int path = rson->path;
1940 int pd = PATHDEPTH(path);
1941 int son_path_ok = false;
1942 int j;
1943
1944 max_path_depth = MAX(max_path_depth,pd);
1945
1946 for (j=0; j<pd; j++)
1947 {
1948 int ns = NEXTSIDE(path,j);
1949
1950 if (nson==NULL)
1951 {
1952 PD_ERR(ER_DBG_RULE_VERBOSE,("ERROR in CheckMRules, elem %ld: nson==NULL for son %d at %d\n",id,s,j),error);
1953 break;
1954 }
1955 if (ns>=SIDES_OF_ELEM(nson))
1956 {
1957 PD_ERR(ER_DBG_RULE_VERBOSE,("ERROR in CheckMRules, elem %ld: path of son %d at %d has invalid side %s\n",id,s,j,ns),error);
1958 break;
1959 }
1960 nson = NBELEM(nson,ns);
1961 }
1962 if (nson!=son)
1963 PD_ERR(ER_DBG_RULE_VERBOSE,("ERROR in CheckMRules, elem %ld: wrong path of son %d\n",id,s),error)
1964 else
1965 son_path_ok = true;
1966 if (!son_path_ok)
1967 some_path_wrong = true;
1968 }*/
1969 #endif
1970 #endif
1971 }
1972 }
1973 /*if (!some_path_wrong)
1974 PD_ERR(ER_DBG_RULE_VERBOSE,("FINE in CheckMRules (%c-rule), elem %ld: all paths ok\n",(refi<UGMAXRULE(tag))?'r':'e',id),error);*/
1975
1976 /* check sonandnode
1977 NOT CONSISTENT in rm since not used (says Stefan, 971219)
1978 for (i=0; i<MAX_NEW_CORNERS_DIM; i++)
1979 if (newnodes[i]!=NULL)
1980 if (sons[mr->sonandnode[i][0]]!=NULL)
1981 if (CORNER(sons[mr->sonandnode[i][0]],mr->sonandnode[i][1]) != newnodes[i])
1982 PD_ERR(ER_DBG_RULE_VERBOSE,("ERROR in CheckMRules, elem %ld: sonandnode %d inconsistent\n",id,i),error);*/
1983
1984 /* summary */
1985 if (error)
1986 {
1987 PD(("in total %3d ERRORS in CheckMRules (%c-rule %3d of %d), elem %ld\n",error,(refi<UGMAXRULE(tag)) ? 'r' : 'e',refi,tag,id));
1988 if (use_bug_in_rule)
1989 bug_in_rule[tag][refi] = true;
1990 }
1991 else
1992 PrintDebug("FINE in CheckMRules (%c-rule %4d of %d), elem %ld: everything ok\n",(refi<UGMAXRULE(tag)) ? 'r' : 'e',refi,tag,id);
1993 }
1994 PrintDebug("CheckMRules: max_path_depth %d\n",max_path_depth);
1995
1996 if (use_bug_in_rule)
1997 {
1998 int t,i;
1999
2000 PrintDebug("--------------- CheckMRules: rules with bugs ---------------\n");
2001
2002 for (t=0; t<TAGS; t++)
2003 for (i=0; i<global.maxrule[t]; i++)
2004 if (bug_in_rule[t][i])
2005 {
2006 PrintDebug("-- rule %4d of %d\n",i,t);
2007 if (i<UGMAXRULE(t))
2008 ShowRefRuleX(t,i,PrintDebug);
2009 else
2010 PrintDebug("sorry, no ShowRefRule for mgio-rules\n");
2011 }
2012
2013 PrintDebug("------------------------------------------------------------\n");
2014
2015 if (ReleaseTmpMem(global.heap,MarkKey))
2016 ASSERT(false);
2017 }
2018 else
2019 PrintDebug("sorry, no storage available for use_bug_in_rule\n");
2020 }
2021
2022 /****************************************************************************/
2023 /*D
2024 NEW_Write_RefRules - write refinement rules of multigrid to file
2025
2026 SYNOPSIS:
2027 INT NEW_Write_RefRules (MULTIGRID *mg, INT RefRuleOffset[], INT MarkKey, MGIO_RR_RULE **mrule_handle)
2028
2029 PAAMETERS:
2030 . mg - multigrid
2031 . RefRuleOffset - tag-wise offsets for refrules
2032 . MarkKey - mark key for temporary storage of ugio
2033 . mrule_handle - pointer to mgio rules allocated on temporary storage of ugio
2034
2035 DESCRIPTION:
2036 Write refinement rules of multigrid to file (rm rules as far as available, the remainder
2037 is covered by rules extracted by the er module).
2038
2039 RETURN VALUE:
2040 INT
2041 .n 0: ok
2042 .n >0: error
2043 D*/
2044 /****************************************************************************/
2045
NEW_Write_RefRules(MULTIGRID * mg,INT RefRuleOffset[],INT MarkKey,MGIO_RR_RULE ** mrule_handle)2046 INT NS_DIM_PREFIX NEW_Write_RefRules (MULTIGRID *mg, INT RefRuleOffset[], INT MarkKey, MGIO_RR_RULE **mrule_handle)
2047 {
2048 MGIO_RR_GENERAL rr_general;
2049 MGIO_RR_RULE *mrule;
2050 INT BotMarkKey;
2051 int tag;
2052
2053 if (mg==NULL)
2054 REP_ERR_RETURN(1);
2055
2056 global.heap = MGHEAP(mg);
2057 if (MarkTmpMem(global.heap,&BotMarkKey))
2058 REP_ERR_RETURN(1);
2059
2060 /* init rule counters (continue with last rule IDs of rm) */
2061 /* TODO (HRR 971207): this is important for matching existing rules after loading a grid (coarsening)
2062 but CAUTION: refine should never address a rule beyond rm rules! */
2063 for (tag=0; tag<TAGS; tag++)
2064 global.maxrule[tag] = UGMAXRULE(tag);
2065
2066 #if (defined __THREEDIM__) || (defined __DEBUG_ER__)
2067 /* incomplete refrule set: extract non-existing rules that are realized in mg */
2068 if (ExtractRules(mg))
2069 REP_ERR_RETURN(1);
2070 #endif
2071
2072 global.maxrules = 0;
2073 for (tag=0; tag<TAGS; tag++)
2074 global.maxrules += global.maxrule[tag];
2075
2076 /* write refrules general */
2077 RefRuleOffset[0] = 0;
2078 for (tag=0; tag<TAGS; tag++)
2079 {
2080 if (tag>0) RefRuleOffset[tag] = RefRuleOffset[tag-1] + global.maxrule[tag-1];
2081 rr_general.RefRuleOffset[tag] = RefRuleOffset[tag];
2082 }
2083 rr_general.nRules = global.maxrules;
2084 if (Write_RR_General(&rr_general))
2085 REP_ERR_RETURN(1);
2086
2087 /* allocate MGIO_RR_RULE table (will stay in scope beyond this module) */
2088 *mrule_handle = (MGIO_RR_RULE*) GetTmpMem(global.heap,global.maxrules*sizeof(MGIO_RR_RULE),MarkKey);
2089 if (*mrule_handle==NULL)
2090 REP_ERR_RETURN(1);
2091
2092 /* write refrules */
2093 mrule = *mrule_handle;
2094 for (tag=0; tag<TAGS; tag++)
2095 {
2096 int i;
2097
2098 /* rm rules */
2099 for (i=0; i<UGMAXRULE(tag); i++)
2100 {
2101 ASSERT(RefRules[tag]+i!=NULL);
2102 URule2Mrule(RefRules[tag]+i,mrule);
2103 mrule++;
2104 }
2105
2106 /* er rules */
2107 for (; i<global.maxrule[tag]; i++)
2108 {
2109 ASSERT(global.hrule[tag][i]!=NULL);
2110 HRule2Mrule(global.hrule[tag][i],mrule);
2111 mrule++;
2112 }
2113
2114 }
2115 Write_RR_Rules(global.maxrules,*mrule_handle);
2116
2117 /* free hrules and hrule table */
2118 if (ReleaseTmpMem(global.heap,BotMarkKey))
2119 REP_ERR_RETURN(1);
2120
2121 IFDEBUG(gm,ER_DBG_GENERAL)
2122 WriteDebugInfo();
2123 ENDDEBUG
2124
2125 IFDEBUG(gm,ER_DBG_RULES)
2126 CheckMRules(mg,RefRuleOffset,*mrule_handle);
2127 ENDDEBUG
2128
2129 return (0);
2130 }
2131
2132 /****************************************************************************/
2133 /*D
2134 ResetRefineTagsBeyondRuleManager - reset refine tags to what refine expects there
2135
2136 SYNOPSIS:
2137 INT ResetRefineTagsBeyondRuleManager (MULTIGRID *mg)
2138
2139 PARAMETERS:
2140 . mg - multigrid
2141
2142 DESCRIPTION:
2143 Reset refine tags to what refine expects there.
2144
2145 RETURN VALUE:
2146 INT
2147 .n 0: ok
2148 D*/
2149 /****************************************************************************/
2150
ResetRefineTagsBeyondRuleManager(MULTIGRID * mg)2151 INT NS_DIM_PREFIX ResetRefineTagsBeyondRuleManager (MULTIGRID *mg)
2152 {
2153 ELEMENT *elem;
2154 int l,n=0;
2155
2156 /** \todo (HRR 971211): don't include TOPLEVEL (no elem refined there) */
2157 for (l=0; l<=TOPLEVEL(mg); l++)
2158 for (elem=PFIRSTELEMENT(GRID_ON_LEVEL(mg,l)); elem!=NULL; elem=SUCCE(elem))
2159 if (BEYOND_UG_RULES(elem))
2160 {
2161 SETREFINE(elem,COPY);
2162 n++;
2163 }
2164 PRINTDEBUG(gm,ER_DBG_GENERAL,("ResetRefineTags: done (for %d elements)\n",n));
2165
2166 return (0);
2167 }
2168