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: ugio.c */
6 /* */
7 /* Purpose: ug's grid input/output */
8 /* */
9 /* Author: Peter Bastian */
10 /* Interdisziplinaeres Zentrum fuer Wissenschaftliches Rechnen */
11 /* Universitaet Heidelberg */
12 /* Im Neuenheimer Feld 368 */
13 /* 6900 Heidelberg */
14 /* email: ug@ica3.uni-stuttgart.de */
15 /* */
16 /* History: 29.01.92 begin, ug version 2.0 */
17 /* */
18 /* Remarks: */
19 /* */
20 /****************************************************************************/
21
22 #define EXTRACT_RULES 1 /* 1: use er.c, 0: use old version */
23
24 /****************************************************************************/
25 /* */
26 /* include files */
27 /* system include files */
28 /* application include files */
29 /* */
30 /****************************************************************************/
31
32 #include <config.h>
33 #include <cstdlib>
34 #include <cstdio>
35 #include <cstring>
36 #include <string>
37 #include <cmath>
38 #include <climits>
39 #include <ctime>
40
41 #include <dune/uggrid/parallel/ppif/ppifcontext.hh>
42
43 #include <dune/uggrid/low/architecture.h>
44 #include <dune/uggrid/low/bio.h>
45 #include <dune/uggrid/low/debug.h>
46 #include <dune/uggrid/low/fifo.h>
47 #include <dune/uggrid/low/fileopen.h>
48 #include <dune/uggrid/low/heaps.h>
49 #include <dune/uggrid/low/misc.h>
50 #include <dune/uggrid/low/ugstruct.h>
51 #include <dune/uggrid/low/ugtypes.h>
52
53 #include <dune/uggrid/ugdevices.h>
54 #ifdef ModelP
55 #include <dune/uggrid/parallel/dddif/parallel.h>
56 #endif
57 #include "gm.h"
58 #include "algebra.h"
59 #include "ugm.h"
60 #include "ugio.h"
61 #include "elements.h"
62 #include "shapes.h"
63 #include "rm.h"
64 #include "mgio.h"
65 #include "mgheapmgr.h"
66
67 /* include refine because of macros accessed */
68 #include "refine.h"
69 #include "rm.h"
70
71 #include "er.h"
72 #include "cw.h"
73
74
75 USING_UG_NAMESPACE
76 USING_UGDIM_NAMESPACE
77 using namespace PPIF;
78
79 /****************************************************************************/
80 /* */
81 /* defines in the following order */
82 /* */
83 /* compile time constants defining static data size (i.e. arrays) */
84 /* other constants */
85 /* macros */
86 /* */
87 /****************************************************************************/
88
89 /* define LOCAL_FILE_SYSTEM if each processor has a local (own) filesystem
90 #define LOCAL_FILE_SYSTEM
91 */
92
93 /* define OPTIMIZED_IO to use optimized configuration */
94 #ifdef ModelP
95 #define OPTIMIZED_IO
96 #endif
97
98 #define BUFFERSIZE 512 /* size of the general purpose text buff*/
99
100 #define HEADER_FMT "# grid on level 0 for %s\n# saved %s\n# %s\n# %s\n"
101 #define BN_HEADER_FMT "\n# boundary nodes\n"
102 #define BN_FMT "bn %d"
103 #define IN_HEADER_FMT "\n# inner nodes\n"
104 #define IN_FMT "in "
105 #define IE_HEADER_FMT "\n# elements\n"
106 #define IE_FMT "ie "
107 #define EOL_FMT ";\n"
108 #define EOF_FMT "# end of file\n"
109
110 /* size of integer list storing processor-ids of copies of objects,sufficient at least for one cg_element or refinement */
111 #define ELEMPROCLISTSIZE 2000 /* n_elem + n_edge + n_node */
112 /* n_elem = 6+30 (6 h-ghosts, 30 v-ghosts) */
113 /* n_edge = 12*30 (12 max_edges_of_elem) (30 probably 30 elements per edge) */
114 /* n_node = 100*8 (8 max_node_of_elem) 100 probably elements per node) */
115 /*#define PROCLISTSIZE (ELEMPROCLISTSIZE*MAX_SONS * 2) size not enough and change to static variable; Christian Wrobel 980128 */
116 #define PROCLISTSIZE_VALUE (ELEMPROCLISTSIZE*MAX_SONS * MAX(13,(int)(2.0+log((double)procs))))
117 #define PROCLISTSIZE proc_list_size
118
119 /* orphan condition for elements */
120 #define EORPHAN(e) (EFATHER(e)==NULL || THEFLAG(e))
121
122 /* define this for reconstruction of vertical node pointers */
123 #define __ConnectVerticalOverlap__
124
125 /****************************************************************************/
126 /* */
127 /* data structures used in this source file (exported data structures are */
128 /* in the corresponding include file!) */
129 /* */
130 /****************************************************************************/
131
132 /****************************************************************************/
133 /* */
134 /* definition of exported global variables */
135 /* */
136 /****************************************************************************/
137
138 /****************************************************************************/
139 /* */
140 /* definition of variables global to this source file only (static!) */
141 /* */
142 /****************************************************************************/
143
144 static INT gridpaths_set=false;
145 static MGIO_RR_RULE *rr_rules;
146 static INT rr_rule_offsets[TAGS];
147 static unsigned short *ProcList, *ActProcListPos;
148 static int foid,non;
149 static NODE **nid_n; /* mapping: orphan-node-id --> orphan node */
150 static NODE **vid_n; /* mapping: orphan-vertex-id --> lowest orphan node */
151 static INT nov; /* number of orphan vertices */
152 static ELEMENT **eid_e;
153 static int nparfiles; /* nb of parallel files */
154 #ifdef __MGIO_PE_INFO__
155 static int npe_info; /* partitioning information? */
156 #endif
157 static int proc_list_size = -1; /* hold the computed value for PROCLISTSIZE; initialized with crazy dummy */
158
159 /****************************************************************************/
160 /* */
161 /* forward declarations of functions used before they are defined */
162 /* */
163 /****************************************************************************/
164
165 static INT WriteElementParInfo (GRID *theGrid, ELEMENT *theElement, MGIO_PARINFO *pinfo);
166
167
168 /************************************************************************************/
169 /* */
170 /* Enumeration of nodes from left to right corresponding to 0...n-1 */
171 /* ---------------------------------------------------------------- */
172 /* */
173 /* orphan nodes */
174 /* |~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~| */
175 /* N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N */
176 /* |~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~|~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~| */
177 /* ghost nodes master nodes */
178 /* */
179 /* */
180 /* */
181 /* Mapping should exist for all orphan nodes to vertices. 'foid' is id of first */
182 /* orphan node, 'non' number of orphan nodes. (MGIO_PAR) it looks (consequently): */
183 /* */
184 /* */
185 /* orphan nodes */
186 /* |~~~~~~~~~~~~~~~~~~~~| */
187 /* N N N N N N N N N N N N N N N N N N N */
188 /* |~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~| */
189 /* master nodes */
190 /* */
191 /* */
192 /************************************************************************************/
193
RenumberNodes(MULTIGRID * theMG,INT * foid,INT * non)194 static INT RenumberNodes (MULTIGRID *theMG, INT *foid, INT *non)
195 {
196 INT i,nid;
197 NODE *theNode;
198
199 nid=0;
200 if (theMG->ppifContext().procs() == 1)
201 {
202 /* ids for all nodes */
203 for (theNode=FIRSTNODE(GRID_ON_LEVEL(theMG,0)); theNode!=NULL; theNode=SUCCN(theNode))
204 {
205 ID(theNode) = ID(MYVERTEX(theNode));
206 nid = MAX(nid,ID(theNode));
207 }
208 nid++;
209 non[0] = nid;
210 for (i=1; i<=TOPLEVEL(theMG); i++)
211 for (theNode=FIRSTNODE(GRID_ON_LEVEL(theMG,i)); theNode!=NULL; theNode=SUCCN(theNode))
212 ID(theNode) = nid++;
213 foid[0] = 0;
214 }
215 else
216 {
217 /* ids for other nodes */
218 for (i=0; i<=TOPLEVEL(theMG); i++)
219 for (theNode=PFIRSTNODE(GRID_ON_LEVEL(theMG,i)); theNode!=NULL; theNode=SUCCN(theNode))
220 if (GHOST(theNode) && !USED(theNode))
221 ID(theNode) = nid++;
222 foid[0] = nid;
223
224 for (i=0; i<=TOPLEVEL(theMG); i++)
225 for (theNode=PFIRSTNODE(GRID_ON_LEVEL(theMG,i)); theNode!=NULL; theNode=SUCCN(theNode))
226 if (GHOST(theNode) && USED(theNode))
227 ID(theNode) = nid++;
228
229 /* ids for master nodes */
230 for (i=0; i<=TOPLEVEL(theMG); i++)
231 for (theNode=FIRSTNODE(GRID_ON_LEVEL(theMG,i)); theNode!=NULL; theNode=SUCCN(theNode))
232 if (MASTER(theNode) && USED(theNode))
233 ID(theNode) = nid++;
234 non[0] = nid-foid[0];
235
236 /* ids for master nodes */
237 for (i=0; i<=TOPLEVEL(theMG); i++)
238 for (theNode=FIRSTNODE(GRID_ON_LEVEL(theMG,i)); theNode!=NULL; theNode=SUCCN(theNode))
239 if (MASTER(theNode) && !USED(theNode))
240 ID(theNode) = nid++;
241
242 }
243
244 return (0);
245 }
246
SaveSurfaceGrid(MULTIGRID * theMG,FILE * stream)247 static INT SaveSurfaceGrid (MULTIGRID *theMG, FILE *stream)
248 {
249 NODE *theNode;
250 ELEMENT *theElement;
251 VERTEX *theVertex;
252 DOUBLE *global;
253 char buffer[BUFFERSIZE];
254 INT i,id,move,l,tl,part;
255
256 tl = CURRENTLEVEL(theMG);
257 for (l=0; l<= tl; l++)
258 for (theElement = FIRSTELEMENT(GRID_ON_LEVEL(theMG,l));
259 theElement != NULL; theElement = SUCCE(theElement))
260 if (NSONS(theElement) == 0)
261 for (i=0; i<CORNERS_OF_ELEM(theElement); i++)
262 ID(MYVERTEX(CORNER(theElement,i))) = 0;
263
264 /* find all boundary nodes witch are no corner nodes */
265 fprintf(stream,BN_HEADER_FMT);
266 id = 0;
267 for (theNode=FIRSTNODE(GRID_ON_LEVEL(theMG,0));
268 theNode!= NULL; theNode=SUCCN(theNode)) {
269 theVertex = MYVERTEX(theNode);
270 if (OBJT(theVertex) == IVOBJ)
271 continue;
272 if (BNDP_BndPDesc(V_BNDP(theVertex),&move,&part))
273 RETURN(1);
274 if (move == 0)
275 ID(theVertex) = id++;
276 }
277 for (l=0; l<= tl; l++)
278 for (theElement = FIRSTELEMENT(GRID_ON_LEVEL(theMG,l));
279 theElement != NULL; theElement = SUCCE(theElement))
280 if (NSONS(theElement) == 0)
281 for (i=0; i<CORNERS_OF_ELEM(theElement); i++) {
282 theVertex = MYVERTEX(CORNER(theElement,i));
283 if (OBJT(theVertex) == IVOBJ)
284 continue;
285 /* skip corner points */
286 if (BNDP_BndPDesc(V_BNDP(theVertex),&move,&part))
287 RETURN(1);
288 if (move == 0)
289 continue;
290 if (ID(theVertex) > 0)
291 continue;
292 ID(theVertex) = id++;
293 if (BNDP_SaveInsertedBndP(V_BNDP(theVertex),
294 buffer,BUFFERSIZE))
295 RETURN(1);
296 fprintf(stream,"%s",buffer);
297 fprintf(stream,EOL_FMT);
298 }
299 /* find all inner nodes */
300 fprintf(stream,IN_HEADER_FMT);
301 for (l=0; l<= tl; l++)
302 for (theElement = FIRSTELEMENT(GRID_ON_LEVEL(theMG,l));
303 theElement != NULL; theElement = SUCCE(theElement))
304 if (NSONS(theElement) == 0)
305 for (i=0; i<CORNERS_OF_ELEM(theElement); i++) {
306 theVertex = MYVERTEX(CORNER(theElement,i));
307 if (OBJT(theVertex) == BVOBJ)
308 continue;
309 if (ID(theVertex) > 0)
310 continue;
311 global = CVECT(theVertex);
312 fprintf(stream,IN_FMT);
313 for (i=0; i<DIM; i++)
314 fprintf(stream," %f",global[i]);
315 fprintf(stream,EOL_FMT);
316 ID(theVertex) = id++;
317 }
318
319 /* elements */
320 fprintf(stream,IE_HEADER_FMT);
321 for (l=0; l<= tl; l++)
322 for (theElement = FIRSTELEMENT(GRID_ON_LEVEL(theMG,l));
323 theElement != NULL; theElement = SUCCE(theElement))
324 if (NSONS(theElement) == 0) {
325 fprintf(stream,IE_FMT);
326 for (i=0; i<CORNERS_OF_ELEM(theElement); i++)
327 fprintf(stream," %d",ID(MYVERTEX(CORNER(theElement,i))));
328 fprintf(stream,EOL_FMT);
329 }
330
331 /* trailer */
332 fprintf(stream,EOF_FMT);
333 fclose(stream);
334 return(GM_OK);
335 }
336
SaveMultiGrid_SCR(MULTIGRID * theMG,const char * name,const char * comment)337 static INT SaveMultiGrid_SCR (MULTIGRID *theMG, const char *name, const char *comment)
338 {
339 FILE *stream;
340 GRID *theGrid;
341 NODE *theNode;
342 ELEMENT *theElement;
343 VERTEX *theVertex;
344 DOUBLE *global;
345 time_t Time;
346 const char *fmt;
347 char buffer[BUFFERSIZE];
348 BVP_DESC theBVPDesc;
349 INT i,id,move,part;
350
351 if (gridpaths_set)
352 /* this way grids are stored to path[0] */
353 stream = FileOpenUsingSearchPaths(name,"w","gridpaths");
354 else
355 stream = fileopen(name,"w");
356 if (stream==NULL)
357 {
358 PrintErrorMessage('E',"SaveMultiGrid","cannot open file");
359 RETURN(GM_FILEOPEN_ERROR);
360 }
361
362 /* get BVPDesc */
363 if (BVP_SetBVPDesc(MG_BVP(theMG),&theBVPDesc))
364 RETURN (GM_ERROR);
365
366 /* get time */
367 fmt = "%a %b %d %H:%M:%S %Y";
368 time(&Time);
369 strftime(buffer,BUFFERSIZE,fmt,localtime(&Time));
370
371 /* write header */
372 fprintf(stream,HEADER_FMT,BVPD_NAME(&theBVPDesc),buffer,name,comment);
373
374 if (TOPLEVEL(theMG) > 0)
375 return(SaveSurfaceGrid(theMG,stream));
376
377 theGrid = GRID_ON_LEVEL(theMG,0);
378
379 /* find all boundary nodes witch are no corner nodes */
380 fprintf(stream,BN_HEADER_FMT);
381 id = 0;
382 for (theNode=FIRSTNODE(theGrid); theNode!= NULL; theNode=SUCCN(theNode))
383 {
384 theVertex = MYVERTEX(theNode);
385 if (OBJT(theVertex) == IVOBJ)
386 continue;
387 if (BNDP_BndPDesc(V_BNDP(theVertex),&move,&part))
388 RETURN(1);
389 if (move == 0)
390 ID(theNode) = id++;
391 }
392 for (theNode=FIRSTNODE(theGrid); theNode!= NULL; theNode=SUCCN(theNode))
393 {
394 theVertex = MYVERTEX(theNode);
395 if (OBJT(theVertex) == IVOBJ)
396 continue;
397 /* skip corner points */
398 if (BNDP_BndPDesc(V_BNDP(theVertex),&move,&part))
399 RETURN(1);
400 if (move == 0)
401 continue;
402 if (BNDP_SaveInsertedBndP(V_BNDP(theVertex),buffer,BUFFERSIZE))
403 RETURN(1);
404 fprintf(stream,"%s",buffer);
405 fprintf(stream,EOL_FMT);
406 ID(theNode) = id++;
407 }
408 /* find all inner nodes */
409 fprintf(stream,IN_HEADER_FMT);
410 for (theNode=FIRSTNODE(theGrid); theNode!= NULL; theNode=SUCCN(theNode))
411 {
412 theVertex = MYVERTEX(theNode);
413 if (OBJT(theVertex) == BVOBJ)
414 continue;
415 global = CVECT(theVertex);
416 fprintf(stream,IN_FMT);
417 for (i=0; i<DIM; i++)
418 fprintf(stream," %f",global[i]);
419 fprintf(stream,EOL_FMT);
420 ID(theNode) = id++;
421 }
422 if (id != NN(theGrid))
423 RETURN(1);
424
425 /* elements */
426 fprintf(stream,IE_HEADER_FMT);
427 for (theElement=FIRSTELEMENT(theGrid); theElement!= NULL;
428 theElement=SUCCE(theElement))
429 {
430 fprintf(stream,IE_FMT);
431 for (i=0; i<CORNERS_OF_ELEM(theElement); i++)
432 fprintf(stream," %d",ID(CORNER(theElement,i)));
433 fprintf(stream,EOL_FMT);
434 }
435
436 /* trailer */
437 fprintf(stream,EOF_FMT);
438 fclose(stream);
439 return(GM_OK);
440 }
441
MarkAsOrphan(GRID * theGrid,ELEMENT * theElement)442 static void MarkAsOrphan( GRID *theGrid, ELEMENT *theElement)
443 {
444 ELEMENT *succe = SUCCE(theElement);
445 ELEMENT *theFather = EFATHER(theElement);
446
447 assert(EGHOST(theElement) || LEVEL(theElement)==0);
448
449 if (theFather != NULL)
450 {
451 int index = PRIO2INDEX(EPRIO(theElement));
452 ELEMENT *Next = NULL;
453 ASSERT(index!=-1 && index<2);
454
455 /* this is an orphan */
456 SETTHEFLAG(theElement,1);
457 if (0)
458 {
459 if (SON(theFather,index) == theElement)
460 {
461 if (succe != NULL)
462 {
463 if (EFATHER(succe)==theFather)
464 if (EPRIO(succe) == EPRIO(theElement))
465 {
466 Next = succe;
467 }
468 }
469 SET_SON(theFather,index,Next);
470 }
471
472 SETNSONS(theFather,NSONS(theFather)-1);
473 SET_EFATHER(theElement,NULL);
474 GRID_UNLINK_ELEMENT(theGrid,theElement);
475 GRID_LINK_ELEMENT(theGrid,theElement,EPRIO(theElement));
476 }
477 PRINTDEBUG(gm,1,("OrphanCons(): new orphan elem=" EID_FMTX "\n",
478 EID_PRTX(theElement)));
479 }
480 }
OrphanCons(MULTIGRID * theMG)481 static INT OrphanCons(MULTIGRID *theMG)
482 {
483 INT i,j,error,orphan;
484 GRID *theGrid;
485 ELEMENT *theElement, *el_father, *nb_el, *nb_el_father;
486 NODE *theNode,*FatherNode;
487 EDGE *theEdge;
488
489 #if defined(ModelP) && defined(__ConnectVerticalOverlap__)
490 if (ConnectVerticalOverlap(theMG)) assert(0);
491 #endif
492
493 error = 0;
494 for (i=0; i<=TOPLEVEL(theMG); i++)
495 {
496 theGrid = GRID_ON_LEVEL(theMG,i);
497 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL; theElement=SUCCE(theElement))
498 {
499 SETTHEFLAG(theElement,0);
500 orphan = 0;
501 for (j=0; j<CORNERS_OF_ELEM(theElement); j++)
502 {
503 theNode = CORNER(theElement,j);
504 switch (NTYPE(theNode))
505 {
506 case (CORNER_NODE) :
507 FatherNode = (NODE *)NFATHER(theNode);
508 if (FatherNode == NULL)
509 {
510 if (EGHOST(theElement)) orphan = 1;
511 else if (LEVEL(theElement) != 0)
512 {
513 error++;
514 PRINTDEBUG(gm,1,("OrphanCons(): ERROR: elem=" EID_FMTX " cornernode=" ID_FMTX " is orphannode\n",
515 EID_PRTX(theElement),ID_PRTX(theNode)));
516 }
517 }
518 else
519 assert (SONNODE(FatherNode) == theNode);
520 break;
521 case (MID_NODE) :
522 theEdge = (EDGE *)NFATHER(theNode);
523 if (theEdge == NULL)
524 {
525 if (EGHOST(theElement)) orphan = 1;
526 else if (LEVEL(theElement) != 0)
527 {
528 error++;
529 PRINTDEBUG(gm,1,("OrphanCons(): ERROR: elem=" EID_FMTX " midnode=" ID_FMTX " is orphannode\n",
530 EID_PRTX(theElement),ID_PRTX(theNode)));
531 }
532 }
533 else
534 assert (MIDNODE(theEdge) == theNode);
535 break;
536 default :
537 break;
538 }
539 }
540
541 /* there exists a second case in which an artificial orphanisation is necessary:
542 if we have 2 neigboring elements, whose fathers are different and are neighbors
543 but have no neighbor-pointers to each other (for examples if both fathers are
544 ghosts since 2 ghosts do not have to have neighbor-pointers to each other).
545 In this situation the 2 children will not get their neighbor-pointers to
546 each other, which is an serious error if a child is a master (masters must
547 have always pointers to all neighbors).
548
549 The idea to solve this problem is: orphanise suitable elements explicitley
550 to force as well the storage of the necessary information as the correct
551 reconstruction while loading.
552
553 TODO: solve the following listed drawbacks. Christian Wrobel 980313.
554
555 //The old idea was:
556 // A preliminary proceeding to do this is: orphanise both children. This increases
557 // the number of orphans more than necessary, but it is safe.
558 //But this is impossible if a child is a master element (masters are not allowed to be
559 //orpahns).
560 //The new idea:
561
562 Orphanise fathers if they are VGhosts or VHGhosts. This situation is handled already
563 correctly.
564
565 In contradiction the requirement to the fathers may be too weak. The faulty
566 situation could perhaps be: also for the fathers holds the same case as for
567 the current children (they have no neighboriung pointers to each other) and
568 thus we rely us for clearing an uncertain case on the already cleared same
569 case for the fathers. This would be perhaps correct if we reconstruct the
570 grid level-wise and ensure the consistancy for this level before proceeding
571 to the next one; but this is not done in ug.
572 */
573 if( !orphan && EMASTER(theElement) ) /* only a master element must know its neighbors */
574 {
575 el_father = EFATHER(theElement);
576 if( el_father != NULL )
577 {
578 if( EMASTER(el_father) )
579 {
580 /* the father will know all his neighbors because he is a master.
581 Will this be true in all situations??? Refer to the above TODO. */
582 }
583 else
584 {
585 ASSERT(EVGHOST(el_father)); /* it is impossible that a master el has a HGHOST-father,
586 must be V[H]GHOST */
587
588 for( j=0; j<SIDES_OF_ELEM(theElement); j++ )
589 {
590 nb_el = NBELEM(theElement,j);
591 if (nb_el == NULL)
592 continue;
593 nb_el_father = EFATHER(nb_el);
594 if( nb_el_father == el_father )
595 continue; /* the elements have the same father; thus they know
596 each other and can set their neighbor-pointers */
597
598 if( nb_el_father==NULL )
599 {
600 /* the neighbor nb_el will be orphanised because of father==NULL
601 and theElement will know an orphan-neighbor */
602 }
603 else if(EMASTER(nb_el_father))
604 {
605 /* the nb-father will know all his neighbors because he is a master.
606 Will this be true in all situations??? Refer to the above TODO. */
607 }
608 else
609 {
610 /* the fathers will perhaps have no neighbor-pointers because they
611 are both ghosts (for which neighbor-pointers isn't mandatory).
612 Orphanise the V[H]Ghost-father to force saving the necessary
613 information. Pure HGhost fathers aren't orphanized. */
614 PRINTDEBUG(gm,1,("OrphanCons(): orphanize father " EID_FMTX " of element " EID_FMTX "\n", EID_PRTX(el_father), EID_PRTX(theElement)));
615 PRINTDEBUG(gm,1,("OrphanCons(): orphanize father " EID_FMTX " of neighbor element " EID_FMTX "\n", EID_PRTX(nb_el_father), EID_PRTX(nb_el)));
616 MarkAsOrphan(theGrid,el_father);
617 /*if(EVGHOST(nb_el_father)) HGHOST-fathers must be orphanized too! */
618 MarkAsOrphan(theGrid,nb_el_father);
619 }
620
621 }
622 /* no need for additional artificial orphanisation */
623 }
624 }
625 else
626 {
627 /* el_father==NULL thus will be treated automatically as orphan
628 by the following code */
629 }
630 }
631 else
632 {
633 /* the element is already orphanised */
634 }
635
636 if (orphan)
637 MarkAsOrphan(theGrid,theElement);
638 }
639 }
640
641 return(error);
642 }
643
644 /****************************************************************************/
645 /** \brief Recalculate ids in the current order
646
647 \param theMG - structure to renumber
648
649 This function recalculates ids in the current order.
650
651 \return 0 if ok
652 */
653 /****************************************************************************/
654
RenumberMultiGrid(MULTIGRID * theMG,INT * nboe,INT * nioe,INT * nbov,INT * niov,NODE *** vid_n,INT * foid,INT * non,INT MarkKey)655 INT NS_DIM_PREFIX RenumberMultiGrid (MULTIGRID *theMG, INT *nboe, INT *nioe, INT *nbov, INT *niov, NODE ***vid_n, INT *foid, INT *non, INT MarkKey)
656 {
657 NODE *theNode;
658 ELEMENT *theElement;
659 INT i,n_ioe,n_boe,vid,n_iov,n_bov,eid,lfoid,lnon,j;
660
661 /* if called from shell, call OrphanCons first */
662 if (nboe==NULL && nioe==NULL && nbov==NULL && niov==NULL && vid_n==NULL && foid==NULL && non==NULL) {
663 i=OrphanCons(theMG);
664 if (i!=0) REP_ERR_RETURN (1);
665 }
666
667 /* init used-flags */
668 for (i=0; i<=TOPLEVEL(theMG); i++)
669 for (theNode=PFIRSTNODE(GRID_ON_LEVEL(theMG,i)); theNode!=NULL; theNode=SUCCN(theNode))
670 {
671 SETUSED(theNode,0);
672 SETUSED(MYVERTEX(theNode),0);
673 SETTHEFLAG(MYVERTEX(theNode),0);
674 }
675
676 /* renumber elements and set orphan-flags for nodes and vertices */
677 eid=n_ioe=n_boe=0;
678 for (i=0; i<=TOPLEVEL(theMG); i++)
679 for (theElement=PFIRSTELEMENT(GRID_ON_LEVEL(theMG,i)); theElement!=NULL; theElement=SUCCE(theElement))
680 if (EFATHER(theElement)==NULL || THEFLAG(theElement)==1)
681 {
682 ID(theElement) = eid++;
683 if (OBJT(theElement)==BEOBJ) n_boe++;
684 else n_ioe++;
685 for (j=0; j<CORNERS_OF_ELEM(theElement); j++)
686 {
687 SETUSED(CORNER(theElement,j),1);
688 SETUSED(MYVERTEX(CORNER(theElement,j)),1);
689 }
690 #ifdef ModelP
691 assert(i==0 || EGHOST(theElement));
692 #endif
693 }
694
695 for (i=0; i<=TOPLEVEL(theMG); i++)
696 for (theElement=PFIRSTELEMENT(GRID_ON_LEVEL(theMG,i)); theElement!=NULL; theElement=SUCCE(theElement))
697 if (EFATHER(theElement)!=NULL && THEFLAG(theElement)==0)
698 ID(theElement) = eid++;
699 if (nboe!=NULL) *nboe = n_boe;
700 if (nioe!=NULL) *nioe = n_ioe;
701
702 /* renumber vertices */
703 vid=n_iov=n_bov=0;
704 for (i=0; i<=TOPLEVEL(theMG); i++)
705 for (theNode=PFIRSTNODE(GRID_ON_LEVEL(theMG,i)); theNode!=NULL; theNode=SUCCN(theNode))
706 if (THEFLAG(MYVERTEX(theNode))==0 && USED(MYVERTEX(theNode)) && OBJT(MYVERTEX(theNode))==BVOBJ)
707 {
708 ID(MYVERTEX(theNode)) = vid++;
709 n_bov++;
710 SETTHEFLAG(MYVERTEX(theNode),1);
711 }
712 for (i=0; i<=TOPLEVEL(theMG); i++)
713 for (theNode=PFIRSTNODE(GRID_ON_LEVEL(theMG,i)); theNode!=NULL; theNode=SUCCN(theNode))
714 if (THEFLAG(MYVERTEX(theNode))==0 && USED(MYVERTEX(theNode)) && OBJT(MYVERTEX(theNode))==IVOBJ)
715 {
716 ID(MYVERTEX(theNode)) = vid++;
717 n_iov++;
718 SETTHEFLAG(MYVERTEX(theNode),1);
719 }
720 if (vid_n!=NULL)
721 {
722 *vid_n = (NODE**)GetTmpMem(MGHEAP(theMG),(n_iov+n_bov)*sizeof(NODE*),MarkKey);
723 for (i=0; i<n_iov+n_bov; i++) (*vid_n)[i] = NULL;
724 for (i=0; i<=TOPLEVEL(theMG); i++)
725 for (theNode=PFIRSTNODE(GRID_ON_LEVEL(theMG,i)); theNode!=NULL; theNode=SUCCN(theNode))
726 if (USED(theNode))
727 {
728 assert(ID(MYVERTEX(theNode))<n_iov+n_bov);
729 if ((*vid_n)[ID(MYVERTEX(theNode))]!=NULL) continue;
730 (*vid_n)[ID(MYVERTEX(theNode))] = theNode;
731 }
732 IFDEBUG(gm,4)
733 for (i=0; i<n_iov+n_bov; i++)
734 assert((*vid_n)[i] != NULL);
735 ENDDEBUG
736 }
737 for (i=0; i<=TOPLEVEL(theMG); i++) /* not neccessary for i/o */
738 for (theNode=PFIRSTNODE(GRID_ON_LEVEL(theMG,i)); theNode!=NULL; theNode=SUCCN(theNode))
739 if (THEFLAG(MYVERTEX(theNode))==0 && !USED(MYVERTEX(theNode)))
740 {
741 ID(MYVERTEX(theNode)) = vid++;
742 SETTHEFLAG(MYVERTEX(theNode),1);
743 }
744 if (nbov!=NULL) *nbov = n_bov;
745 if (niov!=NULL) *niov = n_iov;
746
747 /* renumber nodes */
748 if (RenumberNodes(theMG,&lfoid,&lnon)) return (1);
749 if (foid!=NULL) *foid = lfoid;
750 if (non!=NULL) *non = lnon;
751
752 return (0);
753 }
754
755 #if !EXTRACT_RULES
Write_RefRules(MULTIGRID * theMG,INT * RefRuleOffset,INT MarkKey)756 static INT Write_RefRules (MULTIGRID *theMG, INT *RefRuleOffset, INT MarkKey)
757 {
758 MGIO_RR_GENERAL rr_general;
759 INT i,j,k,t,nRules;
760 HEAP *theHeap;
761 MGIO_RR_RULE *Refrule;
762 REFRULE * ug_refrule;
763 struct mgio_sondata *sonData;
764
765 /* init */
766 if (theMG==NULL) REP_ERR_RETURN(1);
767 theHeap = MGHEAP(theMG);
768
769 /* write refrules general */
770 nRules = 0; RefRuleOffset[0] = 0;
771 for (i=0; i<TAGS; i++)
772 {
773 nRules += MaxRules[i];
774 if (i>0) RefRuleOffset[i] = RefRuleOffset[i-1] + MaxRules[i-1];
775 rr_general.RefRuleOffset[i] = RefRuleOffset[i];
776 }
777 rr_general.nRules = nRules;
778 if (Write_RR_General(&rr_general)) REP_ERR_RETURN(1);
779
780 /* write refrules */
781 rr_rules = Refrule = (MGIO_RR_RULE *)GetTmpMem(theHeap,nRules*sizeof(MGIO_RR_RULE),MarkKey);
782 if (Refrule==NULL) {UserWriteF("ERROR: cannot allocate %d bytes for refrules\n",(int)nRules*sizeof(MGIO_RR_RULE)); REP_ERR_RETURN(1);}
783 for (t=0; t<TAGS; t++)
784 {
785 ug_refrule = RefRules[t];
786 for (i=0; i<MaxRules[t]; i++)
787 {
788 Refrule->rclass = ug_refrule->rclass;
789 Refrule->nsons = ug_refrule->nsons;
790 for (j=0; j<MGIO_MAX_NEW_CORNERS; j++)
791 Refrule->pattern[j] = ug_refrule->pattern[j];
792 for (j=0; j<MGIO_MAX_NEW_CORNERS; j++)
793 {
794 Refrule->sonandnode[j][0] = ug_refrule->sonandnode[j][0];
795 Refrule->sonandnode[j][1] = ug_refrule->sonandnode[j][1];
796 }
797 for (j=0; j<Refrule->nsons; j++)
798 {
799 sonData = &(Refrule->sons[j]);
800 sonData->tag = ug_refrule->sons[j].tag;
801 for (k=0; k<MGIO_MAX_CORNERS_OF_ELEM; k++)
802 sonData->corners[k] = ug_refrule->sons[j].corners[k];
803 for (k=0; k<MGIO_MAX_SIDES_OF_ELEM; k++)
804 sonData->nb[k] = ug_refrule->sons[j].nb[k];
805 sonData->path = ug_refrule->sons[j].path;
806 }
807 Refrule++;
808 ug_refrule++;
809 }
810 }
811 if (Write_RR_Rules(nRules,rr_rules)) REP_ERR_RETURN(1);
812
813 return (0);
814 }
815 #endif
816
SetRefinement(GRID * theGrid,ELEMENT * theElement,NODE ** NodeContext,ELEMENT * SonList[MAX_SONS],INT nmax,MGIO_REFINEMENT * refinement,INT * RefRuleOffset)817 static INT SetRefinement (GRID *theGrid, ELEMENT *theElement,
818 NODE **NodeContext, ELEMENT *SonList[MAX_SONS],
819 INT nmax, MGIO_REFINEMENT *refinement,
820 INT *RefRuleOffset)
821 {
822 MGIO_RR_RULE *theRule;
823 INT i,j,n,sonRefined,sonex,nex;
824 /* ELEMENT *SonSonList[MAX_SONS]; */
825
826 refinement->refclass = REFINECLASS(theElement);
827 refinement->refrule = REFINE(theElement) + RefRuleOffset[TAG(theElement)];
828 theRule = rr_rules + rr_rule_offsets[TAG(theElement)] + REFINE(theElement);
829
830 if (MGIO_PARFILE)
831 {
832 nex=0;
833 for (i=0; i<nmax; i++)
834 if (SonList[i] != NULL)
835 for (j=0; j<CORNERS_OF_TAG(theRule->sons[i].tag); j++)
836 nex |= (1<<theRule->sons[i].corners[j]);
837 }
838 else
839 nex = ~0;
840
841 /* store new corner ids */
842 n=0;
843 for (i=0; i<CORNERS_OF_ELEM(theElement)+EDGES_OF_ELEM(theElement)+SIDES_OF_ELEM(theElement); i++)
844 if (NodeContext[i]!=NULL && ((nex>>i)&0x1))
845 refinement->newcornerid[n++] = ID(NodeContext[i]);
846 i = CORNERS_OF_ELEM(theElement)+CENTER_NODE_INDEX(theElement);
847 if (NodeContext[i]!=NULL && ((nex>>i)&0x1))
848 refinement->newcornerid[n++] = ID(NodeContext[i]);
849 refinement->nnewcorners = n;
850
851 /* sons are refined ? */
852 sonRefined=sonex=refinement->nbid_ex=0;
853 for (i=0,n=0; i<nmax; i++)
854 {
855 if (SonList[i]==NULL) continue;
856
857 /*
858 GetAllSons(SonList[i],SonSonList);
859
860 refined = 0;
861 for (j=0; SonSonList[j]!=NULL; j++)
862 if (!EORPHAN(SonSonList[j])) refined = 1;
863 if (REFINE(SonList[i])!=NO_REFINEMENT && refined) sonRefined |= (1<<i);
864 */
865
866 if (REFINE(SonList[i])!=NO_REFINEMENT) sonRefined |= (1<<i);
867 if (MGIO_PARFILE)
868 {
869 sonex |= (1<<i);
870 if (WriteElementParInfo(theGrid,SonList[i],&refinement->pinfo[i]))
871 REP_ERR_RETURN(1);
872 for (j=0; j<SIDES_OF_ELEM(SonList[i]); j++)
873 if (NBELEM(SonList[i],j)!=NULL && EORPHAN(NBELEM(SonList[i],j)))
874 {
875 refinement->nbid_ex |= (1<<i);
876 refinement->nbid[i][j] = ID(NBELEM(SonList[i],j));
877 }
878 else
879 refinement->nbid[i][j] = -1;
880 }
881 }
882 refinement->sonref = sonRefined;
883 if (MGIO_PARFILE) refinement->sonex = sonex;
884
885 #ifdef __MGIO_PE_INFO__
886 if (npe_info > 1)
887 {
888 int i;
889
890 if (nmax == theRule->nsons)
891 refinement->pe_info = 0;
892 else
893 refinement->pe_info = 1;
894
895 for (i=0; i<theRule->nsons; i++)
896 {
897 if (i<nmax && SonList[i]!=NULL)
898 {
899 refinement->pes[i] = PARTITION(SonList[i]);
900 if (!EMASTER(SonList[i])) refinement->pe_info = 1;
901 }
902 else
903 {
904 refinement->pe_info = 1;
905 refinement->pes[i] = -1;
906 ASSERT(MGIO_PARFILE);
907 }
908 }
909 }
910 else
911 {
912 refinement->pe_info = 0;
913 }
914 #endif
915
916 /* not movable at the moment */
917 refinement->nmoved = 0;
918
919 /* write orphan-sonsonson-node info */
920 if (MGIO_PARFILE)
921 {
922 refinement->orphanid_ex = n = 0;
923 for (i=0; i<CORNERS_OF_ELEM(theElement)+EDGES_OF_ELEM(theElement)+SIDES_OF_ELEM(theElement); i++)
924 {
925 if (NodeContext[i]!=NULL && ((nex>>i)&0x1))
926 {
927 if (ID(MYVERTEX(NodeContext[i]))<nov
928 && LEVEL(NodeContext[i])<LEVEL(vid_n[ID(MYVERTEX(NodeContext[i]))]))
929 {
930 /* TODO remove the next assert; Christian Wrobel 980127
931 it seems that the assert fails for valid load balancings
932 consider the following situation:
933
934 level 3: |ooo* * is orphan node
935 level 2: &------| & is corner node and i points to it; thus i<=CORNERS_OF_ELEM(theElement) is valid
936 level 1: |-------------|
937 level 0: |---------------------| this is also an orphan element but this doesn't influence the current case
938 */
939 /*assert(i>=CORNERS_OF_ELEM(theElement)); */
940
941 refinement->orphanid[n] = ID(vid_n[ID(MYVERTEX(NodeContext[i]))]);
942 refinement->orphanid_ex = 1;
943 }
944 else
945 {
946 refinement->orphanid[n] = -1;
947 }
948 n++;
949 }
950 }
951 i = CORNERS_OF_ELEM(theElement)+CENTER_NODE_INDEX(theElement);
952 if (NodeContext[i]!=NULL && ((nex>>i)&0x1))
953 {
954 if (ID(MYVERTEX(NodeContext[i]))<nov
955 && LEVEL(NodeContext[i])<LEVEL(vid_n[ID(MYVERTEX(NodeContext[i]))]))
956 {
957 refinement->orphanid[n++] = ID(vid_n[ID(MYVERTEX(NodeContext[i]))]);
958 refinement->orphanid_ex = 1;
959 }
960 else
961 {
962 refinement->orphanid[n++] = -1;
963 }
964 }
965 }
966
967 #if (MGIO_DEBUG>0)
968 /* write mykey */
969 refinement->mykey = KeyForObject((KEY_OBJECT *)theElement);
970
971 /* write myfatherkey */
972 if ( THEFLAG(theElement) )
973 refinement->myfatherkey = -1; /* orphan element has no father */
974 else
975 refinement->myfatherkey = KeyForObject((KEY_OBJECT *)EFATHER(theElement));
976
977 /* write nbkey[] */
978 for (j=0; j<SIDES_OF_ELEM(theElement); j++)
979 refinement->nbkey[j] = KeyForObject((KEY_OBJECT *)NBELEM(theElement,j));
980
981 n=0;
982
983 refinement->mycorners = CORNERS_OF_ELEM(theElement);
984
985 /* write mycornerkey[], mycornerfatherkey[] and mycornersonkey[] */
986 for (i=0; i<CORNERS_OF_ELEM(theElement); i++)
987 {
988 refinement->mycornerkey[n] = KeyForObject((KEY_OBJECT *)CORNER(theElement,i));
989 refinement->mycornerfatherkey[n] = KeyForObject((KEY_OBJECT *)NFATHER(CORNER(theElement,i)));
990 refinement->mycornersonkey[n] = KeyForObject((KEY_OBJECT *)SONNODE(CORNER(theElement,i)));
991 n++;
992 }
993
994 if (MGIO_PARFILE)
995 {
996 for (i=0; i<nmax; i++)
997 {
998 if ((sonex>>i)&1)
999 { /* son exists */
1000 refinement->sonskey[i] = KeyForObject((KEY_OBJECT *)SonList[i]);
1001 for (j=0; j<SIDES_OF_ELEM(SonList[i]); j++)
1002 refinement->sonsnbkey[i][j] = KeyForObject((KEY_OBJECT *)NBELEM(SonList[i],j));
1003 }
1004 else
1005 refinement->sonskey[i] = -1;
1006 }
1007 }
1008 #endif
1009
1010 return (0);
1011 }
1012
RemoveOrphanSons(ELEMENT ** SonList,INT * nmax)1013 static INT RemoveOrphanSons (ELEMENT **SonList, INT *nmax)
1014 {
1015 INT i,max;
1016
1017 if (nmax != NULL)
1018 {
1019 for (max=0,i=0; i<*nmax; i++)
1020 {
1021 if (SonList[i]!=NULL)
1022 {
1023 if (THEFLAG(SonList[i]))
1024 SonList[i] = NULL;
1025 else
1026 max = i+1;
1027 }
1028 }
1029 *nmax = max;
1030 }
1031 else
1032 {
1033 for (i=0; SonList[i]!=NULL; i++)
1034 if (SonList[i]!=NULL && THEFLAG(SonList[i])) SonList[i] = NULL;
1035 for (max=0,i=0; i<MAX_SONS; i++)
1036 if (SonList[i]!=NULL)
1037 {
1038 if (max<i) SonList[max] = SonList[i];
1039 max++;
1040 }
1041 }
1042
1043 return(0);
1044 }
1045
SetHierRefinement(GRID * theGrid,ELEMENT * theElement,MGIO_REFINEMENT * refinement,INT * RefRuleOffset)1046 static INT SetHierRefinement (GRID *theGrid, ELEMENT *theElement, MGIO_REFINEMENT *refinement, INT *RefRuleOffset)
1047 {
1048 INT i,nmax;
1049 ELEMENT *theSon;
1050 ELEMENT *SonList[MAX_SONS];
1051 NODE *NodeContext[MAX_NEW_CORNERS_DIM+MAX_CORNERS_OF_ELEM];
1052 MGIO_RR_RULE *theRule;
1053
1054 /*PRINTDEBUG(gm,0,("SetHierRefinement(): level=%d elem=" EID_FMTX "\n",LEVEL(theGrid),EID_PRTX(theElement)));*/
1055
1056 /* sequential part */
1057 if (REFINE(theElement)==NO_REFINEMENT) return (0);
1058 if (GetNodeContext(theElement,NodeContext)) REP_ERR_RETURN(1);
1059 theRule = rr_rules + rr_rule_offsets[TAG(theElement)] + REFINE(theElement);
1060 if (GetOrderedSons(theElement,theRule,NodeContext,SonList,&nmax)) REP_ERR_RETURN(1);
1061 if (RemoveOrphanSons(SonList,&nmax)) REP_ERR_RETURN(1);
1062 if (SetRefinement (theGrid,theElement,NodeContext,SonList,nmax,refinement,RefRuleOffset)) REP_ERR_RETURN(1);
1063 if (Write_Refinement (refinement,rr_rules)) REP_ERR_RETURN(1);
1064
1065 /* recursive call */
1066 for (i=0; i<nmax; i++)
1067 {
1068 theSon = SonList[i];
1069 if (theSon==NULL) continue;
1070 if (REFINE(theSon)!=NO_REFINEMENT)
1071 if (SetHierRefinement(theGrid,theSon,refinement,RefRuleOffset))
1072 REP_ERR_RETURN(1);
1073 }
1074
1075 return (0);
1076 }
1077
nRefinements(ELEMENT * theElement,INT * n)1078 static INT nRefinements (ELEMENT *theElement, INT *n)
1079 {
1080 INT i,nmax;
1081 ELEMENT *theSon;
1082 ELEMENT *SonList[MAX_SONS];
1083 NODE *NodeContext[MAX_NEW_CORNERS_DIM+MAX_CORNERS_OF_ELEM];
1084 MGIO_RR_RULE *theRule;
1085
1086 if (REFINE(theElement)==NO_REFINEMENT) return (0);
1087 if (GetNodeContext(theElement,NodeContext)) REP_ERR_RETURN(1);
1088 theRule = rr_rules + rr_rule_offsets[TAG(theElement)] + REFINE(theElement);
1089 if (GetOrderedSons(theElement,theRule,NodeContext,SonList,&nmax)) REP_ERR_RETURN(1);
1090 if (RemoveOrphanSons(SonList,&nmax)) REP_ERR_RETURN(1);
1091 (*n)++;
1092 for (i=0; i<nmax; i++)
1093 {
1094 theSon = SonList[i];
1095 if (theSon==NULL) continue;
1096 if (REFINE(theSon)!=NO_REFINEMENT)
1097 if (nRefinements(theSon,n)) REP_ERR_RETURN(1);
1098 }
1099
1100 return (0);
1101 }
1102
WriteCG_Vertices(MULTIGRID * theMG,INT MarkKey)1103 static INT WriteCG_Vertices (MULTIGRID *theMG, INT MarkKey)
1104 {
1105 INT i,j,n;
1106 MGIO_CG_POINT *cg_point,*cgp;
1107 HEAP *theHeap;
1108 VERTEX *theVertex;
1109
1110 theHeap = MGHEAP(theMG);
1111 n = nov*(MGIO_CG_POINT_SIZE);
1112 cg_point = (MGIO_CG_POINT *)GetTmpMem(theHeap,n,MarkKey);
1113 if (cg_point==NULL) {UserWriteF("ERROR: cannot allocate %d bytes for cg_points\n",(int)n); REP_ERR_RETURN(1);}
1114 for (i=0; i<=TOPLEVEL(theMG); i++)
1115 for (theVertex=PFIRSTVERTEX(GRID_ON_LEVEL(theMG,i)); theVertex!=NULL; theVertex=SUCCV(theVertex))
1116 if (ID(theVertex)<nov)
1117 {
1118 assert(USED(theVertex));
1119 cgp = MGIO_CG_POINT_PS(cg_point,ID(theVertex));
1120 for (j=0; j<MGIO_DIM; j++)
1121 cgp->position[j] = CVECT(theVertex)[j];
1122 if (MGIO_PARFILE)
1123 {
1124 cgp->level = LEVEL(theVertex);
1125 cgp->prio = 0;
1126 }
1127 }
1128
1129 if (Write_CG_Points((int)nov,cg_point)) REP_ERR_RETURN(1);
1130
1131 return (0);
1132 }
1133
1134 #ifdef ModelP
WriteElementParInfo(GRID * theGrid,ELEMENT * theElement,MGIO_PARINFO * pinfo)1135 static INT WriteElementParInfo (GRID *theGrid,
1136 ELEMENT *theElement, MGIO_PARINFO *pinfo)
1137 {
1138 INT i,j,k,s,n_max;
1139 int *pl;
1140 NODE *theNode;
1141 VERTEX *theVertex;
1142 EDGE *theEdge;
1143 #ifdef ModelP
1144 auto& dddContext = theGrid->dddContext();
1145 #endif
1146
1147 memset(pinfo,0,sizeof(MGIO_PARINFO));
1148
1149 n_max = PROCLISTSIZE-(ActProcListPos-ProcList);
1150
1151 s=0;
1152 pinfo->prio_elem = EPRIO(theElement);
1153 pinfo->ncopies_elem = ENCOPIES(dddContext, theElement);
1154 if (n_max<pinfo->ncopies_elem)
1155 {
1156 PrintErrorMessage('E',"WriteElementParInfo","increase PROCLISTSIZE in gm/ugio.c\n");
1157 REP_ERR_RETURN(1);
1158 }
1159 if (pinfo->ncopies_elem>0)
1160 {
1161 pl = EPROCLIST(dddContext, theElement);
1162 for (i=0,j=2; i<pinfo->ncopies_elem; i++,j+=2)
1163 ActProcListPos[s++] = pl[j];
1164 }
1165 pinfo->e_ident = EGID(theElement);
1166 for (k=0; k<CORNERS_OF_ELEM(theElement); k++)
1167 {
1168 theNode = CORNER(theElement,k);
1169 pinfo->prio_node[k] = PRIO(theNode);
1170 pinfo->ncopies_node[k] = NCOPIES(dddContext, theNode);
1171 if (n_max<pinfo->ncopies_node[k]+s)
1172 {
1173 PrintErrorMessage('E',"WriteElementParInfo","increase PROCLISTSIZE in gm/ugio.c\n");
1174 REP_ERR_RETURN(1);
1175 }
1176 if (pinfo->ncopies_node[k]>0)
1177 {
1178 pl = PROCLIST(dddContext, theNode);
1179 for (i=0,j=2; i<pinfo->ncopies_node[k]; i++,j+=2)
1180 ActProcListPos[s++] = pl[j];
1181 }
1182 pinfo->n_ident[k] = GID(theNode);
1183 }
1184 for (k=0; k<CORNERS_OF_ELEM(theElement); k++)
1185 {
1186 theVertex = MYVERTEX(CORNER(theElement,k));
1187 pinfo->prio_vertex[k] = VXPRIO(theVertex);
1188 pinfo->ncopies_vertex[k] = VXNCOPIES(dddContext, theVertex);
1189 if (n_max<pinfo->ncopies_vertex[k]+s)
1190 {
1191 PrintErrorMessage('E',"WriteElementParInfo","increase PROCLISTSIZE in gm/ugio.c\n");
1192 REP_ERR_RETURN(1);
1193 }
1194 if (pinfo->ncopies_vertex[k]>0)
1195 {
1196 pl = VXPROCLIST(dddContext, theVertex);
1197 for (i=0,j=2; i<pinfo->ncopies_vertex[k]; i++,j+=2)
1198 ActProcListPos[s++] = pl[j];
1199 }
1200 pinfo->v_ident[k] = VXGID(theVertex);
1201 }
1202
1203 #ifdef __TWODIM__
1204 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,EDGEVEC))
1205 {
1206 VECTOR *v;
1207 for (k=0; k<EDGES_OF_ELEM(theElement); k++) {
1208 theEdge=GetEdge(CORNER(theElement,CORNER_OF_EDGE(theElement,k,0)),
1209 CORNER(theElement,CORNER_OF_EDGE(theElement,k,1)));
1210 v = EDVECTOR(theEdge);
1211 pinfo->prio_edge[k] = PRIO(v);
1212 pinfo->ncopies_edge[k] = NCOPIES(dddContext, v);
1213 if (n_max<pinfo->ncopies_edge[k]+s)
1214 {
1215 PrintErrorMessage('E',"WriteElementParInfo","increase PROCLISTSIZE in gm/ugio.c\n");
1216 REP_ERR_RETURN(1);
1217 }
1218 pinfo->ed_ident[k] = GID(v);
1219 if (pinfo->ncopies_edge[k]>0) {
1220 pl = PROCLIST(dddContext, v);
1221 for (i=0,j=2; i<pinfo->ncopies_edge[k]; i++,j+=2)
1222 ActProcListPos[s++] = pl[j];
1223 }
1224 }
1225 }
1226 #endif
1227
1228 #ifdef __THREEDIM__
1229 for (k=0; k<EDGES_OF_ELEM(theElement); k++)
1230 {
1231 theEdge = GetEdge(CORNER(theElement,CORNER_OF_EDGE(theElement,k,0)),
1232 CORNER(theElement,CORNER_OF_EDGE(theElement,k,1)));
1233 pinfo->prio_edge[k] = PRIO(theEdge);
1234 pinfo->ncopies_edge[k] = NCOPIES(dddContext, theEdge);
1235 if (n_max<pinfo->ncopies_edge[k]+s)
1236 {
1237 PrintErrorMessage('E',"WriteElementParInfo","increase PROCLISTSIZE in gm/ugio.c\n");
1238 REP_ERR_RETURN(1);
1239 }
1240 if (pinfo->ncopies_edge[k]>0) {
1241 pl = PROCLIST(dddContext, theEdge);
1242 for (i=0,j=2; i<pinfo->ncopies_edge[k]; i++,j+=2)
1243 ActProcListPos[s++] = pl[j];
1244 }
1245 pinfo->ed_ident[k] = GID(theEdge);
1246 }
1247 #endif
1248
1249 if (s>0) pinfo->proclist = ActProcListPos;
1250 else pinfo->proclist = NULL;
1251 ActProcListPos += s;
1252
1253 return (0);
1254 }
1255 #else
WriteElementParInfo(GRID * theGrid,ELEMENT * theElement,MGIO_PARINFO * pinfo)1256 static INT WriteElementParInfo (GRID *theGrid, ELEMENT *theElement, MGIO_PARINFO *pinfo)
1257 {
1258 REP_ERR_RETURN(1);
1259 }
1260 #endif
1261
SaveMultiGrid_SPF(MULTIGRID * theMG,const char * name,const char * type,const char * comment,INT autosave,INT rename)1262 static INT SaveMultiGrid_SPF (MULTIGRID *theMG, const char *name, const char *type, const char *comment, INT autosave, INT rename)
1263 {
1264 GRID *theGrid;
1265 NODE *theNode;
1266 EDGE *theEdge;
1267 ELEMENT *theElement;
1268 HEAP *theHeap;
1269 MGIO_MG_GENERAL mg_general;
1270 BVP *theBVP;
1271 BVP_DESC theBVPDesc;
1272 MGIO_GE_GENERAL ge_general;
1273 MGIO_GE_ELEMENT ge_element[TAGS];
1274 MGIO_CG_GENERAL cg_general;
1275 MGIO_CG_ELEMENT *cg_element,*cge;
1276 MGIO_REFINEMENT *refinement;
1277 MGIO_BD_GENERAL bd_general;
1278 MGIO_PARINFO cg_pinfo;
1279 INT i,j,k,niv,nbv,nie,nbe,n,nref,hr_max,mode,level,id,foid,non,tl,saved;
1280 int *vidlist;
1281 char *p,*f,*s,*l;
1282 BNDP **BndPList;
1283 char filename[NAMESIZE];
1284 char buf[64],itype[10];
1285 int lastnumber;
1286 INT MarkKey;
1287 #ifdef ModelP
1288 int error;
1289 #ifdef STAT_OUT
1290 DOUBLE ugio_begin,ugio_end;
1291 #endif
1292 #endif
1293
1294
1295 #ifdef STAT_OUT
1296 ugio_begin = CURRENT_TIME;
1297 #endif
1298
1299 /* check */
1300 if (theMG==NULL) REP_ERR_RETURN(1);
1301 if (!MG_COARSE_FIXED(theMG))
1302 {
1303 PrintErrorMessage('E',"SaveMultiGrid_SPF","coarse grid not fixed\n");
1304 REP_ERR_RETURN(1);
1305 }
1306
1307 const auto& procs = theMG->ppifContext().procs();
1308
1309 theHeap = MGHEAP(theMG);
1310 MarkTmpMem(theHeap,&MarkKey);
1311
1312 /* something to do ? */
1313 saved = MG_SAVED(theMG);
1314 #ifdef ModelP
1315 saved = UG_GlobalMinINT(theMG->ppifContext(), saved);
1316 proc_list_size = PROCLISTSIZE_VALUE;
1317 #endif
1318 if (saved)
1319 {
1320 UserWriteF("WARNING: multigrid already saved as %s\n",MG_FILENAME(theMG));
1321 }
1322 /* open file */
1323 nparfiles = procs;
1324 if (autosave)
1325 {
1326 if (name==NULL)
1327 {
1328 if (type!=NULL) REP_ERR_RETURN(1);
1329 strcpy(filename,MG_FILENAME(theMG));
1330 f = strtok(filename,"."); if (f==NULL) REP_ERR_RETURN(1);
1331 s = strtok(NULL,"."); if (s==NULL) REP_ERR_RETURN(1);
1332 l = strtok(NULL,"."); if (l==NULL) REP_ERR_RETURN(1);
1333 l = strtok(NULL,"."); if (l==NULL) REP_ERR_RETURN(1);
1334 l = strtok(NULL,"."); if (l==NULL) REP_ERR_RETURN(1);
1335 strtok(l,"/");
1336 if (sscanf(s,"%d",&lastnumber)!=1) REP_ERR_RETURN(1);
1337 if (lastnumber<0) REP_ERR_RETURN(1);
1338 lastnumber++;
1339 strcpy(itype,l);
1340 sprintf(buf,".%06d",lastnumber);
1341 strcat(filename,buf);
1342 }
1343 else
1344 {
1345 if (type==NULL) REP_ERR_RETURN(1);
1346 strcpy(itype,type);
1347 if (name==NULL) REP_ERR_RETURN(1);
1348 strcpy(filename,name);
1349 strcat(filename,".000000");
1350 }
1351 }
1352 else
1353 {
1354 if (type==NULL) REP_ERR_RETURN(1);
1355 strcpy(itype,type);
1356 if (name==NULL) REP_ERR_RETURN(1);
1357 strcpy(filename,name);
1358 }
1359 if (strcmp(itype,"asc")==0) mode = BIO_ASCII;
1360 else if (strcmp(itype,"bin")==0) mode = BIO_BIN;
1361 else REP_ERR_RETURN(1);
1362 sprintf(buf,".ug.mg.");
1363 strcat(filename,buf);
1364 strcat(filename,itype);
1365 #ifdef ModelP
1366 error = 0;
1367 #ifndef LOCAL_FILE_SYSTEM
1368 if (theMG->ppifContext().isMaster())
1369 #endif
1370 if (MGIO_PARFILE)
1371 if (MGIO_dircreate(filename,(int)rename))
1372 error = -1;
1373 #ifdef LOCAL_FILE_SYSTEM
1374 error = UG_GlobalMinINT(theMG->ppifContext(), error);
1375 #endif
1376 Broadcast(theMG->ppifContext(), &error,sizeof(int));
1377 if (error == -1)
1378 {
1379 UserWriteF("SaveMultiGrid_SPF(): error during file/directory creation\n");
1380 REP_ERR_RETURN(1);
1381 }
1382 if (MGIO_PARFILE)
1383 {
1384 sprintf(buf,"/mg.%04d",(int)theMG->ppifContext().me());
1385 strcat(filename,buf);
1386 }
1387 #endif
1388 if (Write_OpenMGFile (filename,(int)rename)) REP_ERR_RETURN(1);
1389
1390 /* write general information */
1391 theBVP = MG_BVP(theMG);
1392 if (BVP_SetBVPDesc(theBVP,&theBVPDesc)) REP_ERR_RETURN(1);
1393 mg_general.mode = mode;
1394 mg_general.dim = DIM;
1395 Broadcast(theMG->ppifContext(), &MG_MAGIC_COOKIE(theMG),sizeof(INT));
1396 mg_general.magic_cookie = MG_MAGIC_COOKIE(theMG);
1397 mg_general.heapsize = MGHEAP(theMG)->size/KBYTE;
1398 mg_general.nLevel = TOPLEVEL(theMG) + 1;
1399 mg_general.nNode = mg_general.nPoint = mg_general.nElement = 0;
1400 for (i=0; i<=TOPLEVEL(theMG); i++)
1401 {
1402 theGrid = GRID_ON_LEVEL(theMG,i);
1403 mg_general.nNode += NN(theGrid);
1404 mg_general.nPoint += NV(theGrid);
1405 mg_general.nElement += NT(theGrid);
1406 }
1407 strcpy(mg_general.version,MGIO_VERSION);
1408 p = GetStringVar (":IDENTIFICATION");
1409 if (p!=NULL) strcpy(mg_general.ident,p);
1410 else strcpy(mg_general.ident,"---");
1411 strcpy(mg_general.DomainName,BVPD_NAME(&theBVPDesc));
1412 strcpy(mg_general.MultiGridName,MGNAME(theMG));
1413 std::string formatName = "DuneFormat" + std::to_string( DIM ) + "d";
1414 strcpy(mg_general.Formatname,formatName.c_str());
1415 mg_general.VectorTypes = 0;
1416
1417 /* parallel part */
1418 mg_general.nparfiles = nparfiles;
1419 mg_general.me = theMG->ppifContext().me();
1420 #ifdef __MGIO_PE_INFO__
1421 mg_general.npe_info = nparfiles;
1422 npe_info = mg_general.npe_info;
1423 #endif
1424
1425
1426 if (Write_MG_General(&mg_general)) REP_ERR_RETURN(1);
1427
1428 /* write information about general elements */
1429 ge_general.nGenElement = TAGS;
1430 if (Write_GE_General(&ge_general)) REP_ERR_RETURN(1);
1431 for (i=0; i<TAGS; i++)
1432 {
1433 ge_element[i].tag = i;
1434 if (element_descriptors[i]!=NULL)
1435 {
1436 ge_element[i].nCorner = element_descriptors[i]->corners_of_elem;
1437 ge_element[i].nEdge = element_descriptors[i]->edges_of_elem;
1438 ge_element[i].nSide = element_descriptors[i]->sides_of_elem;
1439 for (j=0; j<MGIO_MAX_EDGES_OF_ELEM; j++)
1440 {
1441 ge_element[i].CornerOfEdge[j][0] = element_descriptors[i]->corner_of_edge[j][0];
1442 ge_element[i].CornerOfEdge[j][1] = element_descriptors[i]->corner_of_edge[j][1];
1443 }
1444 for (j=0; j<MGIO_MAX_SIDES_OF_ELEM; j++)
1445 for (k=0; k<MGIO_MAX_CORNERS_OF_SIDE; k++)
1446 ge_element[i].CornerOfSide[j][k] = element_descriptors[i]->corner_of_side[j][k];
1447 }
1448 else
1449 {
1450 ge_element[i].nCorner = 0;
1451 ge_element[i].nEdge = 0;
1452 ge_element[i].nSide = 0;
1453 }
1454 }
1455 if (Write_GE_Elements(TAGS,ge_element)) REP_ERR_RETURN(1);
1456
1457 /* write information about refrules used */
1458 /* TODO (HRR 971209): drop old call of Write_RefRules */
1459 #if EXTRACT_RULES
1460 if (NEW_Write_RefRules(theMG,rr_rule_offsets,MarkKey,&rr_rules)) REP_ERR_RETURN(1);
1461 #else
1462 if (Write_RefRules(theMG,rr_rule_offsets,MarkKey)) REP_ERR_RETURN(1);
1463 #endif
1464
1465 i = OrphanCons(theMG);
1466 if (i)
1467 {
1468 PRINTDEBUG(gm,1,("OrphanCons() returned %d errors\n",i));
1469 fflush(stdout);
1470 assert(0);
1471 }
1472
1473 /* renumber objects */
1474 if (MGIO_PARFILE)
1475 {
1476 if (RenumberMultiGrid (theMG,&nbe,&nie,&nbv,&niv,&vid_n,&foid,&non,MarkKey)) {UserWriteF("ERROR: cannot renumber multigrid\n"); REP_ERR_RETURN(1);}
1477 }
1478 else
1479 {
1480 if (RenumberMultiGrid (theMG,&nbe,&nie,&nbv,&niv,NULL,&foid,&non,MarkKey)) {UserWriteF("ERROR: cannot renumber multigrid\n"); REP_ERR_RETURN(1);}
1481 }
1482 nov = nbv+niv;
1483
1484 /* write general information about coarse grid */
1485 theGrid = GRID_ON_LEVEL(theMG,0);
1486 cg_general.nPoint = nbv+niv;
1487 cg_general.nBndPoint = nbv;
1488 cg_general.nInnerPoint = niv;
1489 cg_general.nElement = nbe+nie;
1490 cg_general.nBndElement = nbe;
1491 cg_general.nInnerElement = nie;
1492 if (Write_CG_General(&cg_general)) REP_ERR_RETURN(1);
1493
1494 #ifdef ModelP
1495 /* this proc has no elements */
1496 if (cg_general.nElement == 0)
1497 {
1498 /* release tmp mem */
1499 ReleaseTmpMem(theHeap,MarkKey);
1500
1501 /* close file */
1502 if (CloseMGFile ()) REP_ERR_RETURN(1);
1503
1504 /* saved */
1505 MG_SAVED(theMG) = 1;
1506 strcpy(MG_FILENAME(theMG),filename);
1507
1508 return(0);
1509 }
1510 #endif
1511
1512 /* write coarse grid points */
1513 if (WriteCG_Vertices(theMG,MarkKey)) REP_ERR_RETURN(1);
1514
1515 /* write mapping: node-id --> vertex-id for orphan-nodes, if(MGIO_PARFILE) */
1516 if (MGIO_PARFILE)
1517 {
1518 vidlist = (int*)GetTmpMem(theHeap,(2+non)*sizeof(int),MarkKey);
1519 vidlist[0] = non;
1520 vidlist[1] = foid;
1521 for (i=0; i<=TOPLEVEL(theMG); i++)
1522 for (theNode=PFIRSTNODE(GRID_ON_LEVEL(theMG,i)); theNode!=NULL; theNode=SUCCN(theNode))
1523 if (USED(theNode))
1524 vidlist[ID(theNode)+2-foid] = ID(MYVERTEX(theNode));
1525 if (Bio_Write_mint((int)(2+non),vidlist)) REP_ERR_RETURN(1);
1526 }
1527
1528 /* write orphan elements */
1529 n = cg_general.nElement; hr_max=0;
1530 cg_element = (MGIO_CG_ELEMENT*)GetTmpMem(theHeap,n*MGIO_CG_ELEMENT_SIZE,MarkKey);
1531 if (cg_element==NULL) {UserWriteF("ERROR: cannot allocate %d bytes for cg_elements\n",(int)n*MGIO_CG_ELEMENT_SIZE); REP_ERR_RETURN(1);}
1532 if (MGIO_PARFILE)
1533 {
1534 ActProcListPos = ProcList = (unsigned short*)GetTmpMem(theHeap,PROCLISTSIZE*sizeof(unsigned short),MarkKey);
1535 if (ProcList==NULL) {UserWriteF("ERROR: cannot allocate %d bytes for ProcList\n",(int)PROCLISTSIZE*sizeof(unsigned short)); REP_ERR_RETURN(1);}
1536 }
1537 else
1538 ActProcListPos = ProcList = NULL;
1539 id=i=0;
1540 for (level=0; level<=TOPLEVEL(theMG); level++)
1541 for (theElement = PFIRSTELEMENT(GRID_ON_LEVEL(theMG,level)); theElement!=NULL; theElement=SUCCE(theElement))
1542 {
1543 /* entry in cg_element-list */
1544 cge = MGIO_CG_ELEMENT_PS(cg_element,i);
1545
1546 /* only orphan elements */
1547 if (!EORPHAN(theElement)) continue;
1548 assert(id==ID(theElement));
1549 assert(id<n);
1550 PRINTDEBUG(gm,1,("write orphan element " EID_FMTX "\n", EID_PRTX(theElement)));
1551
1552 if (MGIO_PARFILE)
1553 {
1554 /* parallel part */
1555 cge->level = LEVEL(theElement);
1556 }
1557
1558 /* --- */
1559 cge->ge = TAG(theElement);
1560 nref=0;
1561 if (nRefinements (theElement,&nref)) REP_ERR_RETURN(1);
1562 hr_max = MAX(hr_max,nref);
1563 cge->nref = nref;
1564 for (j=0; j<CORNERS_OF_ELEM(theElement); j++)
1565 cge->cornerid[j] = ID(CORNER(theElement,j));
1566 for (j=0; j<SIDES_OF_ELEM(theElement); j++)
1567 if (NBELEM(theElement,j)!=NULL)
1568 cge->nbid[j] = ID(NBELEM(theElement,j));
1569 else
1570 cge->nbid[j] = -1;
1571 cge->se_on_bnd = 0;
1572 if (OBJT(theElement)==BEOBJ)
1573 for (j=0; j<SIDES_OF_ELEM(theElement); j++)
1574 if (SIDE_ON_BND(theElement,j))
1575 cge->se_on_bnd |= (1<<j);
1576 for (j=0; j<EDGES_OF_ELEM(theElement); j++)
1577 {
1578 theEdge=GetEdge(CORNER_OF_EDGE_PTR(theElement,j,0),CORNER_OF_EDGE_PTR(theElement,j,1));
1579 assert(theEdge!=NULL);
1580 if (EDSUBDOM(theEdge)==0)
1581 cge->se_on_bnd |= (1<<(SIDES_OF_ELEM(theElement)+j));
1582 }
1583 cge->subdomain = SUBDOMAIN(theElement);
1584 #ifdef __MGIO_PE_INFO__
1585 cge->pe = PARTITION(theElement);
1586 #endif
1587
1588
1589 #if (MGIO_DEBUG>0)
1590 /* write debug extension */
1591 cge->mykey = KeyForObject((KEY_OBJECT *)theElement);
1592 for (j=0; j<CORNERS_OF_ELEM(theElement); j++)
1593 cge->nodekey[j] = KeyForObject((KEY_OBJECT *)CORNER(theElement,j));
1594 for (j=0; j<SIDES_OF_ELEM(theElement); j++)
1595 cge->neighborkey[j] = KeyForObject((KEY_OBJECT *)NBELEM(theElement,j));
1596 #endif
1597
1598 /* increment counters */
1599 i++; id++; /* id this is the id of theElement, according to RenumberElements (ugm.c) */
1600 }
1601 if (Write_CG_Elements((int)i,cg_element)) REP_ERR_RETURN(1);
1602
1603 /* write general bnd information */
1604 if (Bio_Jump_From ()) REP_ERR_RETURN(1);
1605 bd_general.nBndP = nbv;
1606 if (Write_BD_General (&bd_general)) REP_ERR_RETURN(1);
1607
1608 /* write bnd information */
1609 if (nbv > 0) {
1610 BndPList = (BNDP**)GetTmpMem(theHeap,nbv*sizeof(BNDP*),MarkKey);
1611 if (BndPList==NULL) {
1612 UserWriteF("ERROR: cannot allocate %d bytes for BndPList\n",
1613 (int)nbv*sizeof(BNDP*));
1614 REP_ERR_RETURN(1);
1615 }
1616 if (procs>1) i=TOPLEVEL(theMG);
1617 else i=0;
1618 for (level=0; level<=i; level++)
1619 for (theNode=PFIRSTNODE(GRID_ON_LEVEL(theMG,level)); theNode!=NULL; theNode=SUCCN(theNode))
1620 if (USED(MYVERTEX(theNode)) && OBJT(MYVERTEX(theNode))==BVOBJ) {
1621 /* see ugm.c: RenumberVertices */
1622 id = ID(MYVERTEX(theNode));
1623 if (id<0 || id>=nbv) REP_ERR_RETURN(1);
1624 BndPList[id] = V_BNDP(MYVERTEX(theNode));
1625 SETUSED(MYVERTEX(theNode),0);
1626 }
1627 if (Write_PBndDesc (nbv,BndPList)) REP_ERR_RETURN(1);
1628 }
1629 if (Bio_Jump_To ()) REP_ERR_RETURN(1);
1630
1631 /* write parinfo of coarse-grid */
1632 if (MGIO_PARFILE)
1633 {
1634 cg_pinfo.proclist = ProcList;
1635 for (level=0; level<=TOPLEVEL(theMG); level++)
1636 for (theElement = PFIRSTELEMENT(GRID_ON_LEVEL(theMG,level)); theElement!=NULL; theElement=SUCCE(theElement))
1637 if (EORPHAN(theElement))
1638 {
1639 if (WriteElementParInfo(theGrid, theElement,&cg_pinfo)) REP_ERR_RETURN(1);
1640 if (Write_pinfo (TAG(theElement),&cg_pinfo)) REP_ERR_RETURN(1);
1641 }
1642 }
1643
1644 /* save refinement */
1645 refinement = (MGIO_REFINEMENT *)GetTmpMem(theHeap,MGIO_REFINEMENT_SIZE,MarkKey); /* size according to procs==1 or procs>1 (see mgio.h) */
1646 if (refinement==NULL) {UserWriteF("ERROR: cannot allocate %ld bytes for refinement\n",(long)hr_max*MGIO_REFINEMENT_SIZE); REP_ERR_RETURN(1);}
1647 if (procs>1) tl=TOPLEVEL(theMG);
1648 else tl=0;
1649
1650 id=0;
1651 for (level=0; level<=tl; level++) {
1652 theGrid = GRID_ON_LEVEL(theMG,level);
1653 for (theElement = PFIRSTELEMENT(theGrid);
1654 theElement!=NULL; theElement=SUCCE(theElement)) {
1655 if (!EORPHAN(theElement)) continue;
1656 assert(id==ID(theElement));
1657 if (REFINE(theElement)!=NO_REFINEMENT)
1658 /* write refinement only if it is neccessary */
1659 if(SetHierRefinement(theGrid,theElement,refinement,rr_rule_offsets))
1660 REP_ERR_RETURN(1);
1661 id++;
1662 }
1663 }
1664
1665 /* release tmp mem */
1666 ReleaseTmpMem(theHeap,MarkKey);
1667
1668 /* close file */
1669 if (CloseMGFile ()) REP_ERR_RETURN(1);
1670
1671 /* saved */
1672 MG_SAVED(theMG) = 1;
1673 strcpy(MG_FILENAME(theMG),filename);
1674
1675 #if EXTRACT_RULES
1676 ResetRefineTagsBeyondRuleManager(theMG);
1677 #endif
1678
1679 #ifdef STAT_OUT
1680 ugio_end = CURRENT_TIME;
1681
1682 UserWriteF("UGIO: t_iosave=%.2lf\n", ugio_end-ugio_begin);
1683 #endif
1684
1685 return (0);
1686 }
1687
1688 /****************************************************************************/
1689 /** \brief Save complete multigrid structure in a text file
1690
1691 \param theMG - pointer to multigrid
1692 \param name - name of the text file
1693 \param comment - to be included at beginning of file
1694
1695 This function saves the grid on level 0 to a text file.
1696 The text file can be used in a script to load the grid.
1697
1698 \return <ul>
1699 <li> 0 if ok </li>
1700 <li> >0 if error occured </li>
1701 </ul>
1702 */
1703 /****************************************************************************/
1704
SaveMultiGrid(MULTIGRID * theMG,const char * name,const char * type,const char * comment,INT autosave,INT rename)1705 INT NS_DIM_PREFIX SaveMultiGrid (MULTIGRID *theMG, const char *name, const char *type, const char *comment, INT autosave, INT rename)
1706 {
1707 /* check name convention */
1708 if (name==NULL || strcmp(name+strlen(name)-4,".scr")!=0)
1709 {
1710 if (SaveMultiGrid_SPF (theMG,name,type,comment,autosave,rename)) REP_ERR_RETURN(1);
1711 }
1712 else
1713 {
1714 if (SaveMultiGrid_SCR (theMG,name,comment)) REP_ERR_RETURN(1);
1715 }
1716 return (0);
1717 }
1718
Evaluate_pinfo(GRID * theGrid,ELEMENT * theElement,MGIO_PARINFO * pinfo)1719 static INT Evaluate_pinfo (GRID *theGrid, ELEMENT *theElement, MGIO_PARINFO *pinfo)
1720 {
1721 INT i,j,s,prio,where,oldwhere;
1722 INT evec,nvec,edvec,svec;
1723 GRID *vgrid;
1724 ELEMENT *theFather,*After,*Next,*Succe;
1725 NODE *theNode;
1726 VERTEX *theVertex;
1727 VECTOR *theVector;
1728 EDGE *theEdge;
1729
1730 #ifdef ModelP
1731 auto& dddContext = theGrid->dddContext();
1732 #endif
1733
1734 evec = VEC_DEF_IN_OBJ_OF_MG(MYMG(theGrid),ELEMVEC);
1735 nvec = VEC_DEF_IN_OBJ_OF_MG(MYMG(theGrid),NODEVEC);
1736 edvec = VEC_DEF_IN_OBJ_OF_MG(MYMG(theGrid),EDGEVEC);
1737 #ifdef __THREEDIM__
1738 svec = VEC_DEF_IN_OBJ_OF_MG(MYMG(theGrid),SIDEVEC);
1739 #else
1740 svec = 0;
1741 #endif
1742 /* this function does not support side vectors */
1743 /* proclist and identificator need to be stored for each side vector */
1744 if (svec) assert(0);
1745
1746 theFather = EFATHER(theElement);
1747 s = 0;
1748 if ((prio = pinfo->prio_elem) != PrioMaster)
1749 {
1750 oldwhere = PRIO2INDEX(EPRIO(theElement));
1751 Succe = SUCCE(theElement);
1752 GRID_UNLINK_ELEMENT(theGrid,theElement);
1753 SETEPRIO(dddContext, theElement,prio);
1754 if (theFather != NULL)
1755 {
1756 if (theElement == SON(theFather,oldwhere))
1757 {
1758 Next = NULL;
1759 if (Succe != NULL)
1760 if (EFATHER(Succe)==theFather && PRIO2INDEX(EPRIO(Succe))==oldwhere)
1761 Next = Succe;
1762 SET_SON(theFather,oldwhere,Next);
1763 }
1764 where = PRIO2INDEX(prio);
1765 After = SON(theFather,where);
1766 if (After == NULL) SET_SON(theFather,where,theElement);
1767 GRID_LINKX_ELEMENT(theGrid,theElement,prio,After);
1768 }
1769 else
1770 GRID_LINK_ELEMENT(theGrid,theElement,prio);
1771 if (evec)
1772 {
1773 theVector = EVECTOR(theElement);
1774 GRID_UNLINK_VECTOR(theGrid,theVector);
1775 SETPRIO(dddContext, EVECTOR(theElement),prio);
1776 GRID_LINK_VECTOR(theGrid,theVector,prio);
1777 }
1778 }
1779 for (i=0; i<pinfo->ncopies_elem; i++)
1780 {
1781 DDD_IdentifyNumber(dddContext, PARHDRE(theElement),pinfo->proclist[s],pinfo->e_ident);
1782 if (evec)
1783 DDD_IdentifyNumber(dddContext, PARHDR(EVECTOR(theElement)),pinfo->proclist[s],pinfo->e_ident);
1784 s++;
1785 }
1786
1787 for (j=0; j<CORNERS_OF_ELEM(theElement); j++)
1788 {
1789 theNode = CORNER(theElement,j);
1790 if (USED(theNode) == 0)
1791 {
1792 if ((prio = pinfo->prio_node[j]) != PrioMaster)
1793 {
1794 GRID_UNLINK_NODE(theGrid,theNode);
1795 SETPRIO(dddContext, theNode,prio);
1796 GRID_LINK_NODE(theGrid,theNode,prio);
1797 if (nvec)
1798 {
1799 theVector = NVECTOR(theNode);
1800 GRID_UNLINK_VECTOR(theGrid,theVector);
1801 SETPRIO(dddContext, NVECTOR(theNode),prio);
1802 GRID_LINK_VECTOR(theGrid,theVector,prio);
1803 }
1804 }
1805 PRINTDEBUG(gm,1,("Evaluate-pinfo():nid=%d prio=%d\n",ID(theNode),prio);fflush(stdout));
1806 for (i=0; i<pinfo->ncopies_node[j]; i++)
1807 {
1808 DDD_IdentifyNumber(dddContext, PARHDR(theNode),pinfo->proclist[s],pinfo->n_ident[j]);
1809 if (nvec)
1810 DDD_IdentifyNumber(dddContext, PARHDR(NVECTOR(theNode)),pinfo->proclist[s],pinfo->n_ident[j]);
1811 s++;
1812 }
1813 SETUSED(theNode,1);
1814 }
1815 else
1816 s += pinfo->ncopies_node[j];
1817 }
1818 for (j=0; j<CORNERS_OF_ELEM(theElement); j++)
1819 {
1820 theVertex = MYVERTEX(CORNER(theElement,j));
1821 if (USED(theVertex) == 0)
1822 {
1823 vgrid = GRID_ON_LEVEL(MYMG(theGrid),LEVEL(theVertex));
1824 if ((prio = pinfo->prio_vertex[j]) != PrioMaster)
1825 {
1826 GRID_UNLINK_VERTEX(vgrid,theVertex);
1827 SETVXPRIO(dddContext, theVertex,prio);
1828 GRID_LINK_VERTEX(vgrid,theVertex,prio);
1829 }
1830 PRINTDEBUG(gm,1,("Evaluate-pinfo():vid=%d prio=%d\n",ID(theVertex),prio);fflush(stdout));
1831 for (i=0; i<pinfo->ncopies_vertex[j]; i++)
1832 {
1833 DDD_IdentifyNumber(dddContext, PARHDRV(theVertex),pinfo->proclist[s],pinfo->v_ident[j]);
1834 s++;
1835 }
1836 SETUSED(theVertex,1);
1837 }
1838 else
1839 s += pinfo->ncopies_vertex[j];
1840 }
1841
1842 #ifdef __TWODIM__
1843 if (edvec) {
1844 for (j=0; j<EDGES_OF_ELEM(theElement); j++) {
1845 theEdge=GetEdge(CORNER(theElement,CORNER_OF_EDGE(theElement,j,0)),
1846 CORNER(theElement,CORNER_OF_EDGE(theElement,j,1)));
1847 if (USED(theEdge) == 0) {
1848 theVector = EDVECTOR(theEdge);
1849 if ((prio = pinfo->prio_edge[j]) != PrioMaster) {
1850 GRID_UNLINK_VECTOR(theGrid,theVector);
1851 SETPRIO(dddContext, theVector,prio);
1852 GRID_LINK_VECTOR(theGrid,theVector,prio);
1853 }
1854 for (i=0; i<pinfo->ncopies_edge[j]; i++) {
1855 DDD_IdentifyNumber(dddContext, PARHDR(theVector),
1856 pinfo->proclist[s],pinfo->ed_ident[j]);
1857 s++;
1858 }
1859 SETUSED(theEdge,1);
1860 }
1861 else
1862 s += pinfo->ncopies_edge[j];
1863 }
1864 }
1865 #endif
1866
1867 #if (MGIO_DIM==3)
1868 for (j=0; j<EDGES_OF_ELEM(theElement); j++)
1869 {
1870 theEdge = GetEdge(CORNER_OF_EDGE_PTR(theElement,j,0),
1871 CORNER_OF_EDGE_PTR(theElement,j,1));
1872 if (USED(theEdge) == 0)
1873 {
1874 if ((prio = pinfo->prio_edge[j]) != PrioMaster)
1875 {
1876 SETPRIO(dddContext, theEdge,prio);
1877 if (edvec)
1878 {
1879 theVector = EDVECTOR(theEdge);
1880 GRID_UNLINK_VECTOR(theGrid,theVector);
1881 SETPRIO(dddContext, EDVECTOR(theEdge),prio);
1882 GRID_LINK_VECTOR(theGrid,theVector,prio);
1883 }
1884 }
1885 for (i=0; i<pinfo->ncopies_edge[j]; i++)
1886 {
1887 DDD_IdentifyNumber(dddContext, PARHDR(theEdge),pinfo->proclist[s],pinfo->ed_ident[j]);
1888 if (edvec)
1889 DDD_IdentifyNumber(dddContext, PARHDR(EDVECTOR(theEdge)),pinfo->proclist[s],pinfo->ed_ident[j]);
1890 s++;
1891 }
1892 SETUSED(theEdge,1);
1893 }
1894 else
1895 s += pinfo->ncopies_edge[j];
1896 }
1897 #endif
1898
1899 return(0);
1900 }
1901
1902 #ifdef ModelP
Gather_RefineInfo(DDD::DDDContext &,DDD_OBJ obj,void * data)1903 static int Gather_RefineInfo (DDD::DDDContext&, DDD_OBJ obj, void *data)
1904 {
1905 ELEMENT *theElement = (ELEMENT *)obj;
1906
1907 ((int *)data)[0] = REFINECLASS(theElement);
1908 ((int *)data)[1] = REFINE(theElement);
1909 ((int *)data)[2] = MARKCLASS(theElement);
1910 ((int *)data)[3] = MARK(theElement);
1911
1912 return(GM_OK);
1913 }
1914
Scatter_RefineInfo(DDD::DDDContext &,DDD_OBJ obj,void * data)1915 static int Scatter_RefineInfo (DDD::DDDContext&, DDD_OBJ obj, void *data)
1916 {
1917 ELEMENT *theElement = (ELEMENT *)obj;
1918
1919 SETREFINECLASS(theElement,((int *)data)[0]);
1920 SETREFINE(theElement,((int *)data)[1]);
1921 SETMARKCLASS(theElement,((int *)data)[2]);
1922 SETMARK(theElement,((int *)data)[3]);
1923
1924 return(GM_OK);
1925 }
1926
SpreadRefineInfo(GRID * theGrid)1927 static INT SpreadRefineInfo(GRID *theGrid)
1928 {
1929 auto& context = theGrid->dddContext();
1930 const auto& dddctrl = ddd_ctrl(context);
1931
1932 DDD_IFAOneway(context,
1933 dddctrl.ElementIF,GRID_ATTR(theGrid),IF_FORWARD,4*sizeof(INT),
1934 Gather_RefineInfo,Scatter_RefineInfo);
1935 return(GM_OK);
1936 }
1937
Gather_NodeType(DDD::DDDContext &,DDD_OBJ obj,void * data)1938 static int Gather_NodeType (DDD::DDDContext&, DDD_OBJ obj, void *data)
1939 {
1940 NODE *theNode = (NODE *)obj;
1941
1942 ((int *)data)[0] = NTYPE(theNode);
1943
1944 return(GM_OK);
1945 }
1946
Scatter_NodeType(DDD::DDDContext &,DDD_OBJ obj,void * data)1947 static int Scatter_NodeType (DDD::DDDContext&, DDD_OBJ obj, void *data)
1948 {
1949 NODE *theNode = (NODE *)obj;
1950
1951 SETNTYPE(theNode,((int *)data)[0]);
1952
1953 return(GM_OK);
1954 }
1955
SpreadGridNodeTypes(GRID * theGrid)1956 static INT SpreadGridNodeTypes(GRID *theGrid)
1957 {
1958 auto& context = theGrid->dddContext();
1959 const auto& dddctrl = ddd_ctrl(context);
1960
1961 DDD_IFAOneway(context,
1962 dddctrl.NodeIF,GRID_ATTR(theGrid),IF_FORWARD,sizeof(INT),
1963 Gather_NodeType,Scatter_NodeType);
1964 return(GM_OK);
1965 }
1966 #endif
1967
IO_GridCons(MULTIGRID * theMG)1968 static INT IO_GridCons(MULTIGRID *theMG)
1969 {
1970 INT i;
1971 int *proclist;
1972 GRID *theGrid;
1973 ELEMENT *theElement;
1974 VECTOR *theVector;
1975 #ifdef ModelP
1976 auto& dddContext = theMG->dddContext();
1977 #endif
1978 [[maybe_unused]] const auto& me = theMG->ppifContext().me();
1979
1980 for (i=TOPLEVEL(theMG); i>=0; i--) /* propagate information top-down */
1981 {
1982 theGrid = GRID_ON_LEVEL(theMG,i);
1983 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL; theElement=SUCCE(theElement))
1984 {
1985 proclist = EPROCLIST(dddContext, theElement);
1986 while (proclist[0] != -1)
1987 {
1988 if (EMASTERPRIO(proclist[1])) PARTITION(theElement) = proclist[0];
1989 proclist += 2;
1990 }
1991 ASSERT((PARTITION(theElement)==me && EMASTER(theElement))
1992 || (PARTITION(theElement)!=me && !EMASTER(theElement)));
1993 }
1994 for (theVector=PFIRSTVECTOR(theGrid); theVector!=NULL; theVector=SUCCVC(theVector))
1995 if (!MASTER(theVector))
1996 DisposeConnectionFromVector(theGrid,theVector);
1997
1998 #ifdef ModelP
1999 /* spread element refine info */
2000 if (SpreadRefineInfo(theGrid) != GM_OK) RETURN(GM_FATAL);
2001
2002 /* spread nodetypes from master to its copies */
2003 if (SpreadGridNodeTypes(theGrid) != GM_OK) RETURN(GM_FATAL);
2004
2005 /* repair parallel information */
2006 #ifndef OPTIMIZED_IO
2007 ConstructConsistentGrid(theGrid);
2008 #endif
2009 #endif
2010 }
2011 #ifdef ModelP
2012 #ifdef OPTIMIZED_IO
2013 ConstructConsistentMultiGrid(theMG);
2014 #endif
2015 #endif
2016
2017 return(GM_OK);
2018 }
2019
2020
2021 #if (MGIO_DEBUG>0)
CheckNeighborElement(ELEMENT * theElement,INT nb_nr,INT expected_nb_key,const char * messagePrefix)2022 static void CheckNeighborElement( ELEMENT *theElement, INT nb_nr, INT expected_nb_key, const char *messagePrefix )
2023 /* Checks wheather the relation between an element and a given neighbor is ok.
2024 If an error is detected, the program stopps with 'assert'.
2025 */
2026 {
2027 INT k, key;
2028 ELEMENT *nb_back, *nb_element;
2029
2030 nb_element = NBELEM(theElement,nb_nr);
2031
2032 if (nb_element!=NULL)
2033 {
2034 key = KeyForObject((KEY_OBJECT *)nb_element);
2035 if (key!=expected_nb_key)
2036 {
2037 if ( EGHOST(theElement) && EGHOST(nb_element) )
2038 { /* if both elements are ghosts, then the neighbor-pointers are
2039 optionally; but if there is one pointer then also the
2040 other must exist.
2041 if at save time none of this pointers have existed,
2042 at load time they may exist and it is ok.
2043 */
2044
2045 if ( expected_nb_key == -1 )
2046 { /* at save time there were no neighbor pointers
2047 ( i.e. expected_nb_key == -1 ), but now.
2048 ensure that the pointers are symmetric.
2049 */
2050
2051 /* search for the neighbor pointer from nb_element
2052 back to theElement */
2053 nb_back = NULL;
2054 for (k=0; k<SIDES_OF_ELEM(nb_element); k++)
2055 if (NBELEM(nb_element,k)==theElement)
2056 {
2057 nb_back = theElement;
2058 break;
2059 }
2060
2061 if ( nb_back != theElement )
2062 {
2063 printf(" %s neighbor relation unsymmetric, element: " EID_FMTX "has neighbor[%d]: " EID_FMTX " but not vice versa (this pointer has not been in the coarse grid at save time, but it's ok)\n",
2064 messagePrefix,EID_PRTX(theElement),nb_nr,EID_PRTX(nb_element));
2065 assert(0);
2066 }
2067 else
2068 { /* else the neighbor pointer relation is symmetric and hence ok */
2069 return;
2070 }
2071 }
2072 else /* it is an error because there are 2 valid (!=-1) but different keys */
2073 {
2074 printf(" %s neighbor relation with different keys for ghosts, element: " EID_FMTX "has neighbor[%d]: " EID_FMTX "\n",
2075 messagePrefix,EID_PRTX(theElement),nb_nr,EID_PRTX(nb_element));
2076 assert(0);
2077 }
2078 }
2079 else
2080 {
2081 /* else it is an error because at least one of the 2 considered elements
2082 is a master and masters must always have pointers to its
2083 existing neighbors. Thus especially it can't be, that at
2084 save time expected_nb_key==-1 has held and a key difference
2085 is in this case always an error.
2086 */
2087 printf(" %s neighbor relation with different keys for master/ghost, element: " EID_FMTX "has neighbor[%d]: " EID_FMTX ":expected key %d does not match\n",
2088 messagePrefix,EID_PRTX(theElement),nb_nr,EID_PRTX(nb_element),expected_nb_key);
2089 assert(0);
2090 }
2091 }
2092 }
2093 /* you cannot expect that the neighbors are already existing
2094 else
2095 {
2096 if ( expected_nb_key!=-1 && must_exist )
2097 {
2098 printf(" %s neighbor relation, element: " EID_FMTX ", neighbor %d doesn't exist but key %d expected\n",messagePrefix,EID_PRTX(theElement),nb_nr,ref->mykey);
2099 assert(0);
2100 }
2101 }
2102 */
2103 }
2104
CheckLocalElementKeys(ELEMENT * theElement,MGIO_REFINEMENT * ref,INT must_exist)2105 static INT CheckLocalElementKeys (ELEMENT *theElement, MGIO_REFINEMENT *ref, INT must_exist)
2106 /* check as many consistencies as possible. Various keys (for neibhbours,
2107 sons, ...) are stored in the debug mode and are now compared to the
2108 existing objects in the memory.
2109 if must_exist == true all sons must exist. */
2110 {
2111 INT i, j, k, key, nmax;
2112 NODE *corner_node;
2113 ELEMENT *SonList[MAX_SONS], *nb_element, *nb_back;
2114
2115 /* check consistency of element key */
2116 key = KeyForObject((KEY_OBJECT *)theElement);
2117 if (key!=ref->mykey)
2118 {
2119 printf(" IO_Loc: element key, element: " EID_FMTX ": exp.key %d does not match\n",EID_PRTX(theElement),ref->mykey);
2120 assert(0);
2121 }
2122
2123 /* check consistency of father element key;
2124 the check is ok also for orphan elements (they have fatherkey -1) */
2125 key = KeyForObject((KEY_OBJECT *)EFATHER(theElement));
2126 if (key!=ref->myfatherkey)
2127 {
2128 printf(" IO_Loc: father element key, element: " EID_FMTX ", father element: " EID_FMTX ": exp.key %d does not match\n",EID_PRTX(theElement),EID_PRTX(EFATHER(theElement)),ref->mykey);
2129 assert(0);
2130 }
2131
2132 /* check corners */
2133 if (CORNERS_OF_ELEM(theElement)!=ref->mycorners)
2134 {
2135 printf(" IO_Loc: number of corners, element: " EID_FMTX ": exp.key %d does not matchas %d corners but expected %dh\n",EID_PRTX(theElement),CORNERS_OF_ELEM(theElement),ref->mycorners);
2136 assert(0);
2137 }
2138 for (j=0; j<ref->mycorners; j++)
2139 {
2140 corner_node = CORNER(theElement,j);
2141 if (corner_node!=NULL)
2142 {
2143 if (ref->mycornerkey[j] != KeyForObject((KEY_OBJECT *)corner_node))
2144 {
2145 printf(" IO_Loc: corner relation, element: " EID_FMTX ", corner %d: " ID_FMTX " :exp.key %d does not match\n",
2146 EID_PRTX(theElement),j,ID_PRTX(corner_node),ref->mycornerkey[j]);
2147 assert(0);
2148 }
2149
2150 /* check father corner */
2151 if (NFATHER(corner_node)!=NULL)
2152 {
2153 if ( ref->mycornerfatherkey[j] != KeyForObject((KEY_OBJECT *)NFATHER(corner_node)) )
2154 {
2155 if ( ref->mycornerfatherkey[j] != -1 )
2156 {
2157 printf(" IO_Loc: father corner relation, element: " EID_FMTX ", corner[%d]: " ID_FMTX " fathercorner: " ID_FMTX ":exp.key %d does not match\n",
2158 EID_PRTX(theElement),j,ID_PRTX(corner_node),ID_PRTX(NFATHER(corner_node)),ref->mycornerfatherkey[j]);
2159 assert(0);
2160 }
2161 else
2162 {
2163 IFDEBUG(gm,1)
2164 printf(" IO_Loc NOTE: father corner relation, element: " EID_FMTX ", corner[%d]: " ID_FMTX " fathercorner: " ID_FMTX " wasn't here at save time\n",
2165 EID_PRTX(theElement),j,ID_PRTX(corner_node),ID_PRTX(NFATHER(corner_node)));
2166 ENDDEBUG
2167 }
2168 }
2169 }
2170 else
2171 {
2172 if ( ref->mycornerfatherkey[j]!=-1 )
2173 {
2174 if ( !EORPHAN(theElement) )
2175 {
2176 printf(" IO_Loc: father corner relation, element: " EID_FMTX " corner %d doesn't exist but key %d expected\n",EID_PRTX(theElement),j,ref->mycornerfatherkey[j]);
2177 assert(0);
2178 }
2179 /* else for orphan element: at save time there was a father-corner. Such father-corner pointers are
2180 reconstructed not from the orphan element itself, but from the father-corner which is
2181 located in a neighbor element-tree; and since it is possible, that this neighbor element-tree
2182 is not already constructed and the father-corner pointer is still missing. It's OK.
2183 */
2184 }
2185 }
2186
2187 /* check son corner */
2188 if (SONNODE(corner_node)!=NULL)
2189 {
2190 if ( ref->mycornersonkey[j] != KeyForObject((KEY_OBJECT *)SONNODE(corner_node)) )
2191 {
2192 if ( ref->mycornersonkey[j] != -1 )
2193 {
2194 printf(" IO_Loc: son corner relation, element: " EID_FMTX ", corner[%d]: " ID_FMTX " soncorner: " ID_FMTX ":exp.key %d does not match\n",
2195 EID_PRTX(theElement),j,ID_PRTX(corner_node),ID_PRTX(SONNODE(corner_node)),ref->mycornersonkey[j]);
2196 assert(0);
2197 }
2198 else
2199 {
2200 IFDEBUG(gm,1)
2201 printf(" IO_Loc NOTE: son corner relation, element: " EID_FMTX ", corner[%d]: " ID_FMTX " soncorner: " ID_FMTX " wasn't here at save time\n",
2202 EID_PRTX(theElement),j,ID_PRTX(corner_node),ID_PRTX(SONNODE(corner_node)));
2203 ENDDEBUG
2204 }
2205
2206 }
2207 }
2208 /* you cannot expect that the sonnodes are already existing
2209 else
2210 {
2211 if ( ref->mycornersonkey[j]!=-1 && must_exist)
2212 {
2213 printf(" IO_Loc: son corner relation, element: " EID_FMTX ", corner[%d] " ID_FMTX ": son corner doesn't exist but key %d expected\n",EID_PRTX(theElement),j,ID_PRTX(corner_node),ref->mycornersonkey[j]);
2214 assert(0);
2215 }
2216 }
2217 */
2218 }
2219 else
2220 {
2221 if ( ref->mycornerkey[j]!=-1 )
2222 {
2223 printf(" IO_Loc: corner relation, element: " EID_FMTX ", corner %d doesn't exist but key %d expected\n",EID_PRTX(theElement),j,ref->mycornerkey[j]);
2224 assert(0);
2225 }
2226 }
2227
2228 }
2229
2230 /* check nbkey[] */
2231 for (j=0; j<SIDES_OF_ELEM(theElement); j++)
2232 CheckNeighborElement( theElement, j, ref->nbkey[j], "IO_Loc:" );
2233
2234 if (MGIO_PARFILE)
2235 {
2236 NODE *NodeContext[MAX_NEW_CORNERS_DIM+MAX_CORNERS_OF_ELEM];
2237 MGIO_RR_RULE *theRule;
2238
2239 if (GetNodeContext(theElement,NodeContext)) REP_ERR_RETURN(1);
2240 theRule = rr_rules + rr_rule_offsets[TAG(theElement)] + REFINE(theElement);
2241 if (GetOrderedSons(theElement,theRule,NodeContext,SonList,&nmax)) REP_ERR_RETURN(1);
2242
2243 /* check the sons */
2244 for (i=0; i<nmax; i++)
2245 {
2246 if ( SonList[i]==NULL )
2247 {
2248 if (!must_exist)
2249 continue;
2250
2251 if (((1<<i) & ref->sonex)!=0)
2252 {
2253 printf(" IO_Loc: sonlist relation, element: " EID_FMTX ", son %d doesn't exist but sonex 0x%x set\n",EID_PRTX(theElement),i,ref->sonex);
2254 assert(0);
2255 }
2256 if ( ref->sonskey[i] != -1 )
2257 {
2258 printf(" IO_Loc: sonlist relation, element: " EID_FMTX ", son %d doesn't exist but key %d expected\n",EID_PRTX(theElement),i,ref->sonskey[i]);
2259 assert(0);
2260 }
2261 continue;
2262 }
2263
2264 if (((1<<i) & ref->sonex)==0)
2265 {
2266 /* son exists, but is not set in sonex */
2267 /* TODO: check whether SonList[i] is really a (orphan) coarse grid element */
2268 continue;
2269 }
2270
2271
2272 /* check consistency of son element key */
2273 key = KeyForObject((KEY_OBJECT *)SonList[i]);
2274 if (key!=ref->sonskey[i])
2275 {
2276 printf(" IO_Loc: son element key, element: " EID_FMTX ", son[%d] " EID_FMTX ": exp.key %d does not match\n",EID_PRTX(theElement),i,EID_PRTX(SonList[i]),ref->sonskey[i]);
2277 assert(0);
2278 }
2279
2280 for (j=0; j<SIDES_OF_ELEM(SonList[i]); j++)
2281 CheckNeighborElement( SonList[i], j, ref->sonsnbkey[i][j], "IO_Loc: son" );
2282 }
2283 }
2284
2285 return (0);
2286 }
2287 #endif
2288
2289
InsertLocalTree(GRID * theGrid,ELEMENT * theElement,MGIO_REFINEMENT * ref,int * RefRuleOffset)2290 static INT InsertLocalTree (GRID *theGrid, ELEMENT *theElement, MGIO_REFINEMENT *ref, int *RefRuleOffset)
2291 {
2292 INT i,j,k,r_index,nedge,type,sonRefined,n0,n1,Sons_of_Side,SonSides[MAX_SONS],offset;
2293 ELEMENT *theSonElem[MAX_SONS];
2294 ELEMENT *theSonList[MAX_SONS];
2295 VERTEX *theVertex;
2296 NODE *NodeList[MAX_NEW_CORNERS_DIM+MAX_CORNERS_OF_ELEM];
2297 NODE *SonNodeList[MAX_CORNERS_OF_ELEM];
2298 GRID *upGrid;
2299 EDGE *theEdge;
2300 MGIO_RR_RULE *theRule;
2301 struct mgio_sondata *SonData;
2302 INT nbside,nex;
2303 #if (MGIO_DEBUG>0)
2304 MGIO_REFINEMENT ref_copy;
2305 #endif
2306
2307 /* read refinement */
2308 if (Read_Refinement(ref,rr_rules)) REP_ERR_RETURN(1);
2309 #if (MGIO_DEBUG>0)
2310 ref_copy = *ref; /* copy the read ref. into local buffer for subsequent checks */
2311 #endif
2312
2313 /*PRINTDEBUG(gm,0,("InsertLocalTree(): level=%d elem=" EID_FMTX " REFINECLASS %d\n",LEVEL(theGrid),EID_PRTX(theElement),ref->refclass));*/
2314
2315 /* init */
2316 if (ref->refrule==-1) REP_ERR_RETURN(1);
2317 SETREFINE(theElement,ref->refrule-RefRuleOffset[TAG(theElement)]);
2318 SETREFINECLASS(theElement,ref->refclass);
2319 SETMARK(theElement,ref->refrule-RefRuleOffset[TAG(theElement)]);
2320 SETMARKCLASS(theElement,ref->refclass);
2321 theRule = rr_rules+ref->refrule;
2322 upGrid = UPGRID(theGrid);
2323
2324 /* insert nodes */
2325 if (MGIO_PARFILE)
2326 {
2327 nex=0;
2328 for (i=0; i<theRule->nsons; i++)
2329 if ((ref->sonex>>i)&1)
2330 for (j=0; j<CORNERS_OF_TAG(theRule->sons[i].tag); j++)
2331 nex |= (1<<theRule->sons[i].corners[j]);
2332 }
2333 r_index = 0;
2334 for (i=0; i<CORNERS_OF_ELEM(theElement); i++)
2335 {
2336 if (MGIO_PARFILE && !((nex>>i)&1))
2337 {
2338 NodeList[i] = NULL;
2339 continue;
2340 }
2341 NodeList[i] = SONNODE(CORNER(theElement,i));
2342 if (NodeList[i]==NULL)
2343 {
2344 if (ref->newcornerid[r_index]-foid>=0 && ref->newcornerid[r_index]-foid<non)
2345 {
2346 NodeList[i] = nid_n[ref->newcornerid[r_index]-foid];
2347 assert(NodeList[i]!=NULL);
2348 assert(ID(NodeList[i]) == ref->newcornerid[r_index]);
2349 SONNODE(CORNER(theElement,i)) = NodeList[i];
2350 SETNFATHER(NodeList[i],(GEOM_OBJECT *)CORNER(theElement,i));
2351 }
2352 else
2353 {
2354 NodeList[i] = CreateSonNode(upGrid,CORNER(theElement,i));
2355 if (NodeList[i]==NULL) REP_ERR_RETURN(1);
2356 ID(NodeList[i]) = ref->newcornerid[r_index];
2357 }
2358 }
2359 else
2360 {
2361 ASSERT(ID(NodeList[i]) == ref->newcornerid[r_index]);
2362 }
2363 SETNTYPE(NodeList[i],CORNER_NODE);
2364 r_index++;
2365 }
2366 offset = i;
2367
2368 nedge = EDGES_OF_ELEM(theElement);
2369 for (; i<nedge+offset; i++)
2370 {
2371 theVertex = NULL;
2372 if (MGIO_PARFILE)
2373 {
2374 if (!((nex>>i)&1))
2375 {
2376 NodeList[i] = NULL;
2377 continue;
2378 }
2379 else if (ref->orphanid_ex==1)
2380 {
2381 if (ref->orphanid[r_index]>=0)
2382 {
2383 assert(ref->orphanid[r_index]>=foid && ref->orphanid[r_index]<non+foid);
2384 theVertex = MYVERTEX(nid_n[ref->orphanid[r_index]-foid]);
2385 }
2386 }
2387 }
2388 if (theRule->pattern[i-offset]!=1)
2389 {
2390 NodeList[i] = NULL;
2391 continue;
2392 }
2393 n0 = CORNER_OF_EDGE(theElement,i-offset,0);
2394 n1 = CORNER_OF_EDGE(theElement,i-offset,1);
2395 theEdge = GetEdge(CORNER(theElement,n0),CORNER(theElement,n1));
2396 if (theEdge==NULL) REP_ERR_RETURN(1);
2397 NodeList[i] = MIDNODE(theEdge);
2398 if (NodeList[i]==NULL)
2399 {
2400 if (ref->newcornerid[r_index]-foid>=0 && ref->newcornerid[r_index]-foid<non)
2401 {
2402 NodeList[i] = nid_n[ref->newcornerid[r_index]-foid];
2403 assert(NodeList[i]!=NULL);
2404 assert(ID(NodeList[i]) == ref->newcornerid[r_index]);
2405 MIDNODE(theEdge) = NodeList[i];
2406 SETNFATHER(NodeList[i],(GEOM_OBJECT *)theEdge);
2407 }
2408 else
2409 {
2410 NodeList[i] = CreateMidNode(upGrid,theElement,theVertex,i-offset);
2411 if (NodeList[i]==NULL) REP_ERR_RETURN(1);
2412 ID(NodeList[i]) = ref->newcornerid[r_index];
2413 }
2414 }
2415 else
2416 {
2417 ASSERT(ID(NodeList[i]) == ref->newcornerid[r_index]);
2418 }
2419 SETNTYPE(NodeList[i],MID_NODE);
2420 r_index++;
2421 }
2422 offset = i;
2423
2424 #ifdef __THREEDIM__
2425 for (; i<SIDES_OF_ELEM(theElement)+offset; i++)
2426 {
2427 if (theRule->pattern[i-CORNERS_OF_ELEM(theElement)]!=1)
2428 {
2429 NodeList[i] = NULL;
2430 continue;
2431 }
2432 theVertex = NULL;
2433 if (MGIO_PARFILE)
2434 {
2435 if (!((nex>>i)&1))
2436 {
2437 NodeList[i] = NULL;
2438 continue;
2439 }
2440 else if (ref->orphanid_ex==1)
2441 {
2442 if (ref->orphanid[r_index]>=0)
2443 {
2444 assert(ref->orphanid[r_index]>=foid && ref->orphanid[r_index]<non+foid);
2445 theVertex = MYVERTEX(nid_n[ref->orphanid[r_index]-foid]);
2446 }
2447 }
2448 }
2449 if (ref->newcornerid[r_index]-foid>=0 && ref->newcornerid[r_index]-foid<non)
2450 {
2451 NodeList[i] = nid_n[ref->newcornerid[r_index]-foid];
2452 assert(NodeList[i]!=NULL);
2453 assert(ID(NodeList[i]) == ref->newcornerid[r_index]);
2454 }
2455 else
2456 NodeList[i] = GetSideNode(theElement,i-offset);
2457 if (NodeList[i]==NULL)
2458 {
2459 NodeList[i] = CreateSideNode(upGrid,theElement,theVertex,i-offset);
2460 if (NodeList[i]==NULL) REP_ERR_RETURN(1);
2461 ID(NodeList[i]) = ref->newcornerid[r_index];
2462 }
2463 else
2464 SETONNBSIDE(MYVERTEX(NodeList[i]),i-offset);
2465 SETNTYPE(NodeList[i],SIDE_NODE);
2466 r_index++;
2467 }
2468 offset = i;
2469 #endif
2470
2471 i = CORNERS_OF_ELEM(theElement)+CENTER_NODE_INDEX(theElement); /* using consistency within rm.c (id of CenterNode in SonData corresponds to CORNERS_OF_ELEM(theElement)+CENTER_NODE_INDEX(theElement)) */
2472 if (theRule->pattern[CENTER_NODE_INDEX(theElement)]==1 && (!MGIO_PARFILE || ((nex>>(i))&1)))
2473 {
2474 theVertex = NULL;
2475 if (MGIO_PARFILE && ref->orphanid_ex==1 && ref->orphanid[r_index]>=0)
2476 {
2477 assert(ref->orphanid[r_index]>=foid && ref->orphanid[r_index]<non+foid);
2478 theVertex = MYVERTEX(nid_n[ref->orphanid[r_index]-foid]);
2479 }
2480 if (ref->newcornerid[r_index]-foid>=0 && ref->newcornerid[r_index]-foid<non)
2481 {
2482 NodeList[i] = nid_n[ref->newcornerid[r_index]-foid];
2483 assert(NodeList[i]!=NULL);
2484 assert(ID(NodeList[i]) == ref->newcornerid[r_index]);
2485 }
2486 else {
2487 if (theVertex != NULL) /* NEU Ch. Wieners */
2488 /* if (!VXGHOST(theVertex)) */ /* NEU Ch. Wieners */
2489 VFATHER(theVertex) = theElement; /* NEU Ch. Wieners */
2490 NodeList[i] = CreateCenterNode(upGrid,theElement,theVertex);
2491 }
2492 if (NodeList[i]==NULL) REP_ERR_RETURN(1);
2493 SETNTYPE(NodeList[i],CENTER_NODE);
2494 ID(NodeList[i]) = ref->newcornerid[r_index++];
2495 }
2496 else
2497 NodeList[i] = NULL;
2498
2499 /* insert elements */
2500 for (i=0; i<theRule->nsons; i++)
2501 {
2502 if (MGIO_PARFILE && !((ref->sonex>>i)&1))
2503 {
2504 theSonList[i] = NULL;
2505 continue;
2506 }
2507 SonData = theRule->sons+i;
2508 type = IEOBJ;
2509 if (OBJT(theElement)==BEOBJ)
2510 for (j=0; j<SIDES_OF_TAG(SonData->tag); j++)
2511 if (SonData->nb[j]>=MGIO_FATHER_SIDE_OFFSET && SIDE_ON_BND(theElement,SonData->nb[j]-MGIO_FATHER_SIDE_OFFSET))
2512 {
2513 type = BEOBJ;
2514 break;
2515 }
2516 for (j=0; j<CORNERS_OF_TAG(SonData->tag); j++)
2517 SonNodeList[j] = NodeList[SonData->corners[j]];
2518 theSonList[i] = CreateElement(upGrid,SonData->tag,type,SonNodeList,theElement,1);
2519 if (theSonList[i]==NULL) REP_ERR_RETURN(1);
2520 SETECLASS(theSonList[i],ref->refclass);
2521 SETSUBDOMAIN(theSonList[i],SUBDOMAIN(theElement));
2522 }
2523
2524 /* identify objects */
2525 if (MGIO_PARFILE)
2526 for (i=0; i<theRule->nsons; i++)
2527 {
2528 if (theSonList[i] != NULL)
2529 Evaluate_pinfo(upGrid,theSonList[i],ref->pinfo+i);
2530 }
2531
2532 /* set neighbor relation between sons */
2533 for (i=0; i<theRule->nsons; i++)
2534 {
2535 if (MGIO_PARFILE && !((ref->sonex>>i)&1)) continue;
2536 SonData = theRule->sons+i;
2537 for (j=0; j<SIDES_OF_TAG(SonData->tag); j++)
2538 if (SonData->nb[j]<MGIO_FATHER_SIDE_OFFSET && (!MGIO_PARFILE || ((ref->sonex>>(SonData->nb[j]))&1)))
2539 SET_NBELEM(theSonList[i],j,theSonList[SonData->nb[j]]);
2540 else
2541 SET_NBELEM(theSonList[i],j,NULL);
2542 }
2543
2544 /* connect to neighbors */
2545 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
2546 {
2547 for (k=0,j=0; j<theRule->nsons; j++)
2548 {
2549 if (!MGIO_PARFILE || ((ref->sonex>>j)&1))
2550 theSonElem[k++] = theSonList[j];
2551 }
2552 theSonElem[k] = NULL;
2553 if (Get_Sons_of_ElementSide (theElement,i,&Sons_of_Side,theSonElem,SonSides,0,1)) REP_ERR_RETURN(1);
2554
2555 if (Connect_Sons_of_ElementSide(UPGRID(theGrid),theElement,i,Sons_of_Side,theSonElem,SonSides,1)) REP_ERR_RETURN(1);
2556 }
2557
2558 /* connect to orphan-neighbors */
2559 if (MGIO_PARFILE)
2560 for (i=0; i<theRule->nsons; i++)
2561 if (((ref->sonex>>i)&1) && ((ref->nbid_ex>>i)&1))
2562 {
2563 SonData = theRule->sons+i;
2564 for (j=0; j<SIDES_OF_TAG(SonData->tag); j++)
2565 if (ref->nbid[i][j]>=0)
2566 {
2567 SET_NBELEM(theSonList[i],j,eid_e[ref->nbid[i][j]]);
2568 GetNbSideByNodes(eid_e[ref->nbid[i][j]],&nbside,theSonList[i],j);
2569 SET_NBELEM(eid_e[ref->nbid[i][j]],nbside,theSonList[i]);
2570 }
2571 }
2572
2573 #if (MGIO_DEBUG>0)
2574 if (CheckLocalElementKeys(theElement, &ref_copy, false)==0)
2575 {
2576 PRINTDEBUG(gm,4,("InsertLocalTree: preliminary CheckLocalElementKeys ok.\n"));
2577 }
2578 else
2579 {
2580 PrintErrorMessage('E',"InsertLocalTree","preliminary CheckLocalElementKeys failed\n");
2581 }
2582 #endif
2583
2584 /* jump to the sons ? */
2585 sonRefined = ref->sonref;
2586
2587 /* call recursively */
2588 for (i=0; i<theRule->nsons; i++)
2589 {
2590 if (sonRefined & (1<<i))
2591 {
2592 if (InsertLocalTree (upGrid,theSonList[i],ref,RefRuleOffset)) REP_ERR_RETURN(1);
2593 }
2594 else if (theSonList[i] != NULL)
2595 {
2596 SETREFINE(theSonList[i],NO_REFINEMENT);
2597 }
2598
2599 #if (MGIO_DEBUG>0)
2600 if (CheckLocalElementKeys(theElement, &ref_copy, true)==0)
2601 {
2602 PRINTDEBUG(gm,4,("InsertLocalTree: after son %d intermediate CheckLocalElementKeys ok.\n",i));
2603 }
2604 else
2605 {
2606 PrintErrorMessage('E',"InsertLocalTree","intermediate CheckLocalElementKeys failed\n");
2607 }
2608 #endif
2609 }
2610
2611 #if (MGIO_DEBUG>0)
2612 if (CheckLocalElementKeys(theElement, &ref_copy, true)==0)
2613 {
2614 PRINTDEBUG(gm,4,("InsertLocalTree: final CheckLocalElementKeys ok.\n"));
2615 }
2616 else
2617 {
2618 PrintErrorMessage('E',"InsertLocalTree","final CheckLocalElementKeys failed\n");
2619 }
2620 #endif
2621
2622 return (0);
2623 }
2624
2625 #ifdef ModelP
Gather_EClasses(DDD::DDDContext &,DDD_OBJ obj,void * data)2626 static int Gather_EClasses (DDD::DDDContext&, DDD_OBJ obj, void *data)
2627 {
2628 ELEMENT *p;
2629 int *d;
2630
2631 p = (ELEMENT *)obj;
2632 d = (int *)data;
2633 *d = ECLASS(p);
2634
2635 return (0);
2636 }
2637
Scatter_EClasses(DDD::DDDContext &,DDD_OBJ obj,void * data)2638 static int Scatter_EClasses(DDD::DDDContext&, DDD_OBJ obj, void *data)
2639 {
2640 ELEMENT *p;
2641 int *d;
2642
2643 p = (ELEMENT *)obj;
2644 d = (int *)data;
2645 SETECLASS(p,*d);
2646
2647 return (0);
2648 }
2649
CommunicateEClasses(MULTIGRID * theMG)2650 void CommunicateEClasses (MULTIGRID *theMG)
2651 {
2652 auto& context = theMG->dddContext();
2653 const auto& dddctrl = ddd_ctrl(context);
2654
2655 DDD_IFOneway(context,
2656 dddctrl.ElementVHIF,IF_FORWARD,sizeof(int),
2657 Gather_EClasses, Scatter_EClasses);
2658 return;
2659 }
2660 #endif
2661
2662 #if (MGIO_DEBUG>0)
CheckCGKeys(INT ne,ELEMENT ** eid_e,MGIO_CG_ELEMENT * cg_elem)2663 static INT CheckCGKeys (INT ne, ELEMENT** eid_e, MGIO_CG_ELEMENT *cg_elem)
2664 {
2665 INT i,j,k,key;
2666 MGIO_CG_ELEMENT *cge;
2667 ELEMENT *theElement, *nb_element, *nb_back;
2668
2669 for (i=0; i<ne; i++)
2670 {
2671 cge = MGIO_CG_ELEMENT_PS(cg_elem,i);
2672 theElement = eid_e[i];
2673
2674 /* check consistency of keys */
2675 key = KeyForObject((KEY_OBJECT *)theElement);
2676 if (key!=-1 && key!=cge->mykey)
2677 {
2678 printf(" IO_CG: element key, element: " EID_FMTX ": exp.key %d does not match\n",EID_PRTX(theElement),cge->mykey);
2679 assert(0);
2680 }
2681
2682 key = KeyForObject((KEY_OBJECT *)EFATHER(theElement));
2683 if (key!=-1)
2684 {
2685 printf(" IO_CG: father relation, element: " EID_FMTX ", father: " EID_FMTX " ; but no father expected because element is orphan\n",
2686 EID_PRTX(theElement),EID_PRTX(EFATHER(theElement)));
2687
2688 assert(0);
2689 }
2690
2691 for (j=0; j<CORNERS_OF_ELEM(theElement); j++)
2692 {
2693 key = KeyForObject((KEY_OBJECT *)CORNER(theElement,j));
2694 if (key!=-1 && key!=cge->nodekey[j])
2695 {
2696 printf(" IO_CG: corner relation, element: " EID_FMTX ", corner[%d]: " ID_FMTX " :exp.key %d does not match\n",
2697 EID_PRTX(theElement),j,ID_PRTX(CORNER(theElement,j)),cge->nodekey[j]);
2698 assert(0);
2699 }
2700 }
2701
2702 for (j=0; j<SIDES_OF_ELEM(theElement); j++)
2703 CheckNeighborElement( theElement, j, cge->neighborkey[j], "IO_CG:" );
2704 }
2705
2706 return (0);
2707 }
2708 #else
CheckCGKeys(INT ne,ELEMENT ** eid_e,MGIO_CG_ELEMENT * cg_elem)2709 static INT CheckCGKeys (INT ne, ELEMENT** eid_e, MGIO_CG_ELEMENT *cg_elem)
2710 {
2711 return (0);
2712 }
2713 #endif
2714
2715 /****************************************************************************/
2716 /** \brief Load complete multigrid structure from a text file
2717
2718 \param MultigridName - Name of the new 'MULTIGRID' structure in memory.
2719 \param FileName - Name of the file to be read.
2720 \param BVPName - `Name` of the BVP used for the 'MULTIGRID'.
2721 \param format - `Name` of the 'FORMAT' to be used for the 'MULTIGRID'.
2722 \param heapSize - Size of the heap in bytes that will be allocated for the 'MULTIGRID'.
2723
2724 This function can read grid files produced with the 'SaveMultiGrid' function.
2725
2726 \return
2727 NULL if an error occured else a pointer to new 'MULTIGRID'
2728 */
2729 /****************************************************************************/
2730
2731
LoadMultiGrid(const char * MultigridName,const char * name,const char * type,const char * BVPName,const char * format,unsigned long heapSize,INT force,INT optimizedIE,INT autosave,std::shared_ptr<PPIF::PPIFContext> ppifContext)2732 MULTIGRID * NS_DIM_PREFIX LoadMultiGrid (const char *MultigridName,
2733 const char *name,
2734 const char *type,
2735 const char *BVPName,
2736 const char *format,
2737 unsigned long heapSize,
2738 INT force,
2739 INT optimizedIE,
2740 INT autosave,
2741 std::shared_ptr<PPIF::PPIFContext> ppifContext)
2742 /* Documentation of the intended program flow resp. communication requirements.
2743 Functions introducing a global communication (all processors without any exception)
2744 are:
2745 UG_GlobalMin* and DDD_XferBegin/End
2746 the corresponding topics in the following lists are preceeded by "#"
2747 The sequence of main topics is:
2748 read coarse grid data
2749 CreateMultiGrid
2750 DisposeGrid(grid_on_level_0) (caution: for gridlevels>0 DisposeTopLevel is called which involves UG_GlobalMinINT)
2751 create all grid levels
2752 insert coarse mesh
2753
2754 #CreateAlgebra (incl. DDD_IF* and
2755 SetSurfaceClasses incl. UG_GlobalMinINT)
2756 PrepareAlgebraModification
2757
2758 DDD_IdentifyBegin
2759 read parinfo of coarse-grid
2760 read refinements and build local tree (fine grid)
2761 DDD_IdentifyEnd
2762
2763 repair inconsistencies:
2764 CommunicateEClasses (incl. DDD_IF*)
2765 #IO_GridCons (incl. ConstructConsistentGrid
2766 incl. DDD_IF* and DDD_XferBegin/End )
2767 Propagate*VectorClasses (incl. DDD_IF*)
2768 PrepareAlgebraModification
2769 #SetSurfaceClasses (incl. UG_GlobalMinINT)
2770
2771 if you must/will leave this sequence premature ENSURE THE COMMUNICATION PATTERN!!!
2772 i.e. do all global communications as dummies and
2773 all element orientated communication (like DDD_IF*) if there exist elements
2774 */
2775 {
2776 MULTIGRID *theMG;
2777 GRID *theGrid;
2778 ELEMENT *theElement,*ENext;
2779 NODE *theNode;
2780 EDGE *theEdge;
2781 VECTOR *theVector;
2782 HEAP *theHeap;
2783 MGIO_MG_GENERAL mg_general;
2784 MGIO_GE_GENERAL ge_general;
2785 MGIO_GE_ELEMENT ge_element[TAGS];
2786 MGIO_RR_GENERAL rr_general;
2787 MGIO_CG_GENERAL cg_general;
2788 MGIO_CG_POINT *cg_point,*cgp;
2789 MGIO_CG_ELEMENT *cg_element,*cge;
2790 MGIO_BD_GENERAL bd_general;
2791 MGIO_PARINFO cg_pinfo;
2792 MGIO_REFINEMENT *refinement;
2793 BNDP **BndPList;
2794 DOUBLE *Positions;
2795 BVP *theBVP;
2796 BVP_DESC theBVPDesc;
2797 MESH theMesh;
2798 char FormatName[NAMESIZE], BndValName[NAMESIZE], MGName[NAMESIZE], filename[NAMESIZE];
2799 INT i,j,*Element_corner_uniq_subdom, *Ecusdp[2],**Enusdp[2],**Ecidusdp[2],
2800 **Element_corner_ids_uniq_subdom,*Element_corner_ids,max,**Element_nb_uniq_subdom,
2801 *Element_nb_ids,level;
2802 INT *Element_SideOnBnd_uniq_subdom,*ESoBusdp[2];
2803 char buf[64],itype[10];
2804 int *vidlist;
2805 INT MarkKey;
2806 #ifdef __THREEDIM__
2807 ELEMENT *theNeighbor;
2808 INT k;
2809 #endif
2810 #ifdef STAT_OUT
2811 DOUBLE ugio_begin, ugio_end;
2812
2813 ugio_begin = CURRENT_TIME;
2814 #endif
2815
2816 if (not ppifContext)
2817 ppifContext = std::make_shared<PPIF::PPIFContext>();
2818
2819 const auto& me = ppifContext->me();
2820 const auto& procs = ppifContext->procs();
2821
2822 if (autosave)
2823 {
2824 if (name==NULL)
2825 {
2826 return (NULL);
2827 }
2828 else
2829 {
2830 if (type==NULL) return (NULL);
2831 strcpy(itype,type);
2832 if (name==NULL) return (NULL);
2833 strcpy(filename,name);
2834 strcat(filename,".000000");
2835 }
2836 }
2837 else
2838 {
2839 if (type==NULL) return (NULL);
2840 strcpy(itype,type);
2841 if (name==NULL) return (NULL);
2842 strcpy(filename,name);
2843 }
2844 sprintf(buf,".ug.mg.");
2845 strcat(filename,buf);
2846 strcat(filename,itype);
2847
2848 #ifdef ModelP
2849 proc_list_size = PROCLISTSIZE_VALUE;
2850 if (ppifContext->isMaster())
2851 {
2852 #endif
2853 nparfiles = 1;
2854 if (MGIO_filetype(filename) == FT_DIR)
2855 {
2856 sprintf(buf,"/mg.%04d",(int)me);
2857 strcat(filename,buf);
2858 if (Read_OpenMGFile (filename)) {nparfiles = -1;}
2859 else
2860 if (Read_MG_General(&mg_general))
2861 {
2862 CloseMGFile ();
2863 nparfiles = -1;
2864 #ifdef ModelP
2865 assert(0);
2866 #endif
2867 }
2868
2869 nparfiles = mg_general.nparfiles;
2870 if (mg_general.nparfiles>procs) {UserWrite("ERROR: too many processors needed\n");CloseMGFile (); nparfiles = -1;}
2871 assert(mg_general.me == me);
2872
2873 }
2874 else if(MGIO_filetype(filename) == FT_FILE)
2875 {
2876 if (Read_OpenMGFile (filename)) {nparfiles = -1;}
2877 else
2878 if (Read_MG_General(&mg_general))
2879 {
2880 CloseMGFile ();
2881 nparfiles = -1;
2882 #ifdef ModelP
2883 assert(0);
2884 #endif
2885 }
2886 }
2887 else
2888 nparfiles = -1;
2889 #ifdef ModelP
2890 Broadcast(*ppifContext, &nparfiles,sizeof(int));
2891 }
2892 else
2893 {
2894 sprintf(buf,"/mg.%04d",(int)me); /* Also me>=nparfiles needs its filename to be stored in the multigrid */
2895 strcat(filename,buf);
2896 Broadcast(*ppifContext, &nparfiles,sizeof(int));
2897 if (me < nparfiles)
2898 {
2899 if (Read_OpenMGFile (filename)) {nparfiles = -1;}
2900 else
2901 if (Read_MG_General(&mg_general)) {CloseMGFile (); nparfiles = -1;}
2902
2903 }
2904 }
2905 nparfiles = UG_GlobalMinINT(*ppifContext, nparfiles);
2906 #endif
2907 if (nparfiles == -1)
2908 {
2909 UserWriteF("ERROR in LoadMultiGrid: could not open file %s.\n", filename );
2910 return(NULL);
2911 }
2912
2913 if (procs>nparfiles)
2914 {
2915 Broadcast(*ppifContext, &mg_general,sizeof(MGIO_MG_GENERAL));
2916 if (me < nparfiles)
2917 mg_general.me = me;
2918 }
2919 if (mg_general.dim!=DIM) {
2920 UserWrite("ERROR: wrong dimension\n");
2921 CloseMGFile ();
2922 return (NULL);
2923 }
2924 if (strcmp(mg_general.version,MGIO_VERSION)!=0 && force==0) {
2925 UserWrite("ERROR: wrong version\n");
2926 CloseMGFile ();
2927 return (NULL);
2928 }
2929 /* BVP and format */
2930 if (BVPName==NULL) strcpy(BndValName,mg_general.DomainName);
2931 else strcpy(BndValName,BVPName);
2932 if (MultigridName==NULL) strcpy(MGName,mg_general.MultiGridName);
2933 else strcpy(MGName,MultigridName);
2934 if (format==NULL) strcpy(FormatName,mg_general.Formatname);
2935 else strcpy(FormatName,format);
2936 if (heapSize==0) heapSize = mg_general.heapsize * KBYTE;
2937
2938 /* create a virginenal multigrid on the BVP */
2939 theMG = CreateMultiGrid(MGName,BndValName,FormatName,true,false, ppifContext);
2940 if (theMG==NULL) {
2941 UserWrite("ERROR(ugio): cannot create multigrid\n");
2942 CloseMGFile ();
2943 return (NULL);
2944 }
2945 MG_MAGIC_COOKIE(theMG) = mg_general.magic_cookie;
2946
2947 if (me >= nparfiles)
2948 {
2949 /* in this case no CloseMGFile() may be used because no Read_OpenMGFile() was executed */
2950
2951 if (DisposeGrid(GRID_ON_LEVEL(theMG,0))) {DisposeMultiGrid(theMG); return (NULL);}
2952
2953 for (i=0; i<mg_general.nLevel; i++)
2954 if (CreateNewLevel(theMG)==NULL) {DisposeMultiGrid(theMG); return (NULL);}
2955
2956 /* no coarse mesh */
2957
2958 if (CreateAlgebra (theMG)) {DisposeMultiGrid(theMG); return (NULL);}
2959 if (DisposeBottomHeapTmpMemory(theMG)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
2960 if (PrepareAlgebraModification(theMG)) {DisposeMultiGrid(theMG); return (NULL);}
2961
2962 #ifdef ModelP
2963 if (DisposeBottomHeapTmpMemory(theMG)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
2964
2965 DDD_IdentifyBegin(theMG->dddContext());
2966 /* no elements to insert */
2967 if (MGCreateConnection(theMG)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
2968 DDD_IdentifyEnd(theMG->dddContext());
2969
2970 if (MGIO_PARFILE)
2971 {
2972 /* replaced by IO_GridCons because of uniformity; Christian Wrobel 980311. TODO remove:
2973 for (i=0; i<mg_general.nLevel; i++)
2974 ConstructConsistentGrid(GRID_ON_LEVEL(theMG,i));
2975 */
2976 if (IO_GridCons(theMG)) {DisposeMultiGrid(theMG); return (NULL);}
2977 }
2978 #endif
2979
2980 if (SetSurfaceClasses (theMG)) {DisposeMultiGrid(theMG); return (NULL);}
2981
2982 /* saved */
2983 MG_SAVED(theMG) = 1;
2984 strcpy(MG_FILENAME(theMG),filename);
2985
2986 return (theMG);
2987
2988 } /* else if (me < nparfiles) */
2989 theHeap = MGHEAP(theMG);
2990 MarkKey = MG_MARK_KEY(theMG);
2991 if (DisposeGrid(GRID_ON_LEVEL(theMG,0))) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
2992 if (CreateNewLevel(theMG)==NULL) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
2993 theHeap = MGHEAP(theMG);
2994 theBVP = MG_BVP(theMG);
2995 if (theBVP==NULL) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
2996 if (BVP_SetBVPDesc(theBVP,&theBVPDesc)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
2997
2998 /* read general element information */
2999 if (Read_GE_General(&ge_general)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3000 if (Read_GE_Elements(TAGS,ge_element)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3001
3002 /* read refrule information */
3003 if (Read_RR_General(&rr_general)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3004 rr_rules = (MGIO_RR_RULE *)GetTmpMem(theHeap,rr_general.nRules*sizeof(MGIO_RR_RULE),MarkKey);
3005 if (rr_rules==NULL) {UserWriteF("ERROR: cannot allocate %d bytes for rr_rules\n",(int)rr_general.nRules*sizeof(MGIO_RR_RULE)); CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3006 if (Read_RR_Rules(rr_general.nRules,rr_rules)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3007
3008 /* read general information about coarse grid */
3009 if (Read_CG_General(&cg_general)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3010
3011 #ifdef ModelP
3012 /* this proc has no elements */
3013 if (cg_general.nElement == 0)
3014 { /* no elements for me; may skip element orientated communication routines */
3015
3016 ReleaseTmpMem(theHeap,MarkKey);
3017
3018 /* no coarse mesh */
3019
3020 /* replaced by IO_GridCons because of uniformity; Christian Wrobel 980311. TODO remove:
3021 ConstructConsistentGrid(GRID_ON_LEVEL(theMG,0));
3022 */
3023 /* create levels */
3024 for (i=1; i<mg_general.nLevel; i++)
3025 {
3026 if (CreateNewLevel(theMG)==NULL) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3027 /* replaced by IO_GridCons because of uniformity; Christian Wrobel 980311. TODO remove:
3028 ConstructConsistentGrid(GRID_ON_LEVEL(theMG,i));
3029 */
3030 }
3031
3032 if (CreateAlgebra (theMG)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3033 if (DisposeBottomHeapTmpMemory(theMG)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3034 if (PrepareAlgebraModification(theMG)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3035
3036 #ifdef ModelP
3037 if (DisposeBottomHeapTmpMemory(theMG)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3038 #endif
3039
3040 DDD_IdentifyBegin(theMG->dddContext());
3041 /* no elements to insert */
3042 if (MGCreateConnection(theMG)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3043 DDD_IdentifyEnd(theMG->dddContext());
3044
3045 if (MGIO_PARFILE)
3046 if (IO_GridCons(theMG)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3047
3048
3049 if (SetSurfaceClasses (theMG)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3050
3051 if (CloseMGFile ()) {DisposeMultiGrid(theMG); return (NULL);}
3052
3053 /* saved */
3054 MG_SAVED(theMG) = 1;
3055 strcpy(MG_FILENAME(theMG),filename);
3056
3057 return (theMG);
3058 }
3059 #endif
3060
3061 /* read coarse grid points and elements */
3062 cg_point = (MGIO_CG_POINT *)GetTmpMem(theHeap,cg_general.nPoint*sizeof(MGIO_CG_POINT),MarkKey);
3063 if (cg_point==NULL) {UserWriteF("ERROR: cannot allocate %d bytes for cg_points\n",(int)cg_general.nPoint*sizeof(MGIO_CG_POINT)); CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3064 if (Read_CG_Points(cg_general.nPoint,cg_point)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3065
3066 if (MGIO_PARFILE)
3067 {
3068 /* read mapping: node-id --> vertex-id for orphan-nodes */
3069 if (Bio_Read_mint(1,&non)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3070 if (Bio_Read_mint(1,&foid)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3071 vidlist = (int*)GetTmpMem(theHeap,non*sizeof(int),MarkKey);
3072 if (Bio_Read_mint(non,vidlist)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3073 #ifdef Debug
3074 for (i=0; i<non; i++) PRINTDEBUG(gm,1,("LoadMG(): vidList[%d]=%d\n",i,vidlist[i]));
3075 #endif
3076 }
3077 else
3078 {
3079 foid = 0;
3080 non = cg_general.nPoint;
3081 }
3082
3083 cg_element = (MGIO_CG_ELEMENT *)GetTmpMem(theHeap,cg_general.nElement*sizeof(MGIO_CG_ELEMENT),MarkKey);
3084 if (cg_element==NULL) {UserWriteF("ERROR: cannot allocate %d bytes for cg_elements\n",(int)cg_general.nElement*sizeof(MGIO_CG_ELEMENT)); CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3085 if (Read_CG_Elements(cg_general.nElement,cg_element)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3086
3087 /* read general bnd information */
3088 if (Bio_Jump (0)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3089 if (Read_BD_General (&bd_general)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3090
3091 /* read bnd points */
3092 if (bd_general.nBndP > 0) {
3093 BndPList = (BNDP**)GetTmpMem(theHeap,bd_general.nBndP*sizeof(BNDP*),MarkKey);
3094 if (BndPList==NULL) {
3095 UserWriteF("ERROR: cannot allocate %d bytes for BndPList\n",
3096 (int)bd_general.nBndP*sizeof(BNDP*));
3097 CloseMGFile ();
3098 DisposeMultiGrid(theMG);
3099 return (NULL);
3100 }
3101 if (Read_PBndDesc (theBVP,theHeap,bd_general.nBndP,BndPList)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3102 }
3103
3104 /* create and insert coarse mesh: prepare */
3105 theMesh.nBndP = cg_general.nBndPoint;
3106 theMesh.theBndPs = BndPList;
3107 theMesh.nInnP = cg_general.nInnerPoint;
3108 if (cg_general.nInnerPoint>0)
3109 {
3110 theMesh.Position = (DOUBLE**)GetTmpMem(theHeap,cg_general.nInnerPoint*sizeof(DOUBLE*),MarkKey);
3111 if (theMesh.Position==NULL) {UserWriteF("ERROR: cannot allocate %d bytes for theMesh.Position\n",(int)cg_general.nInnerPoint*sizeof(DOUBLE*)); CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3112 Positions = (DOUBLE*)GetTmpMem(theHeap,MGIO_DIM*cg_general.nInnerPoint*sizeof(DOUBLE),MarkKey);
3113 if (Positions==NULL) {UserWriteF("ERROR: cannot allocate %d bytes for Positions\n",(int)MGIO_DIM*cg_general.nInnerPoint*sizeof(DOUBLE)); CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3114 }
3115 if (MGIO_PARFILE)
3116 {
3117 theMesh.VertexLevel = (char*)GetTmpMem(theHeap,(cg_general.nBndPoint+cg_general.nInnerPoint)*sizeof(char),MarkKey);
3118 if (theMesh.VertexLevel==NULL) {UserWriteF("ERROR: cannot allocate %d bytes for VertexLevel\n",(int)(cg_general.nBndPoint+cg_general.nInnerPoint)*sizeof(char)); CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3119 }
3120 else
3121 theMesh.VertexLevel = NULL;
3122 theMesh.nSubDomains = theBVPDesc.nSubDomains;
3123 theMesh.nElements = (INT*)GetTmpMem(theHeap,(theMesh.nSubDomains+1)*sizeof(INT),MarkKey);
3124 if (theMesh.nElements==NULL) {UserWriteF("ERROR: cannot allocate %d bytes for theMesh.nElements\n",(int)(theMesh.nSubDomains+1)*sizeof(INT)); CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3125 theMesh.ElementLevel = (char**)GetTmpMem(theHeap,(theMesh.nSubDomains+1)*sizeof(char*),MarkKey);
3126 if (theMesh.ElementLevel==NULL) {UserWriteF("ERROR: cannot allocate %d bytes for theMesh.ElementLevel\n",(int)(theMesh.nSubDomains+1)*sizeof(char*)); CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3127 for (i=0; i<=theMesh.nSubDomains; i++)
3128 {
3129 theMesh.nElements[i] = 0;
3130 theMesh.ElementLevel[i] = NULL;
3131 }
3132 theMesh.nElements[1] = cg_general.nElement;
3133 theMesh.ElementLevel[1] = (char*)GetTmpMem(theHeap,cg_general.nElement*sizeof(char),MarkKey);
3134 if (theMesh.ElementLevel[1]==NULL) {UserWriteF("ERROR: cannot allocate %d bytes for theMesh.ElementLevel[1]\n",(int)cg_general.nElement*sizeof(char)); CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3135 for (i=0; i<cg_general.nInnerPoint; i++)
3136 theMesh.Position[i] = Positions+MGIO_DIM*i;
3137 for (i=0; i<cg_general.nInnerPoint; i++)
3138 {
3139 cgp = MGIO_CG_POINT_PS(cg_point,cg_general.nBndPoint+i);
3140 for (j=0; j<MGIO_DIM; j++)
3141 theMesh.Position[i][j] = cgp->position[j];
3142 }
3143 if (MGIO_PARFILE)
3144 for (i=0; i<cg_general.nBndPoint+cg_general.nInnerPoint; i++)
3145 {
3146 cgp = MGIO_CG_POINT_PS(cg_point,i);
3147 theMesh.VertexLevel[i] = cgp->level;
3148 }
3149
3150 /* nb of corners of elements */
3151 Element_corner_uniq_subdom = (INT*)GetTmpMem(theHeap,cg_general.nElement*sizeof(INT),MarkKey);
3152 if (Element_corner_uniq_subdom==NULL)
3153 {
3154 UserWriteF("ERROR: cannot allocate %d bytes for Element_corner_uniq_subdom\n",(int)cg_general.nElement*sizeof(INT));
3155 CloseMGFile ();
3156 DisposeMultiGrid(theMG);
3157 return (NULL);
3158 }
3159 Element_SideOnBnd_uniq_subdom = (INT*)GetTmpMem(theHeap,cg_general.nElement*sizeof(INT),MarkKey);
3160 if (Element_SideOnBnd_uniq_subdom==NULL)
3161 {
3162 UserWriteF("ERROR: cannot allocate %d bytes for Element_SideOnBnd_uniq_subdom\n",(int)cg_general.nElement*sizeof(INT));
3163 CloseMGFile ();
3164 DisposeMultiGrid(theMG);
3165 return (NULL);
3166 }
3167 max = 0;
3168 for (i=0; i<cg_general.nElement; i++)
3169 {
3170 cge = MGIO_CG_ELEMENT_PS(cg_element,i);
3171 Element_corner_uniq_subdom[i] = ge_element[cge->ge].nCorner;
3172 max = MAX(max,ge_element[cge->ge].nCorner);
3173 max = MAX(max,ge_element[cge->ge].nSide);
3174 if (MGIO_PARFILE) theMesh.ElementLevel[1][i] = cge->level;
3175 else theMesh.ElementLevel[1][i] = 0;
3176 Element_SideOnBnd_uniq_subdom[i] = cge->se_on_bnd;
3177 }
3178 Ecusdp[1] = Element_corner_uniq_subdom;
3179 theMesh.Element_corners = Ecusdp;
3180 ESoBusdp[1] = Element_SideOnBnd_uniq_subdom;
3181 theMesh.ElemSideOnBnd = ESoBusdp;
3182
3183 /* corners ids of elements */
3184 Element_corner_ids_uniq_subdom = (INT**)GetTmpMem(theHeap,cg_general.nElement*sizeof(INT*),MarkKey);
3185 if (Element_corner_ids_uniq_subdom==NULL) {UserWriteF("ERROR: cannot allocate %d bytes for Element_corner_ids_uniq_subdom\n",(int)cg_general.nElement*sizeof(INT*)); CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3186 Element_corner_ids = (INT*)GetTmpMem(theHeap,max*cg_general.nElement*sizeof(INT),MarkKey);
3187 if (Element_corner_ids==NULL) {UserWriteF("ERROR: cannot allocate %d bytes for Element_corner_ids\n",(int)max*cg_general.nElement*sizeof(INT)); CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3188 for (i=0; i<cg_general.nElement; i++)
3189 Element_corner_ids_uniq_subdom[i] = Element_corner_ids+max*i;
3190 for (i=0; i<cg_general.nElement; i++)
3191 {
3192 cge = MGIO_CG_ELEMENT_PS(cg_element,i);
3193 if (MGIO_PARFILE)
3194 for (j=0; j<Element_corner_uniq_subdom[i]; j++)
3195 {
3196 Element_corner_ids_uniq_subdom[i][j] = vidlist[cge->cornerid[j]-foid];
3197 PRINTDEBUG(gm,1,("LoadMultiGrid(): cg_elem=%d cg_nid[%d]=%d foid=%d\n",i,j,cge->cornerid[j],foid));
3198 }
3199 else
3200 for (j=0; j<Element_corner_uniq_subdom[i]; j++)
3201 Element_corner_ids_uniq_subdom[i][j] = cge->cornerid[j];
3202 }
3203 Ecidusdp[1] = Element_corner_ids_uniq_subdom;
3204 theMesh.Element_corner_ids = Ecidusdp;
3205
3206 /* nb of elements */
3207 Element_nb_uniq_subdom = (INT**)GetTmpMem(theHeap,cg_general.nElement*sizeof(INT*),MarkKey);
3208 if (Element_nb_uniq_subdom==NULL) {UserWriteF("ERROR: cannot allocate %d bytes for Element_nb_uniq_subdom\n",(int)cg_general.nElement*sizeof(INT*)); CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3209 Element_nb_ids = (INT*)GetTmpMem(theHeap,max*cg_general.nElement*sizeof(INT),MarkKey);
3210 if (Element_nb_ids==NULL) {UserWriteF("ERROR: cannot allocate %d bytes for Element_nb_ids\n",(int)max*cg_general.nElement*sizeof(INT)); CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3211 for (i=0; i<cg_general.nElement; i++)
3212 Element_nb_uniq_subdom[i] = Element_nb_ids+max*i;
3213 for (i=0; i<cg_general.nElement; i++)
3214 {
3215 cge = MGIO_CG_ELEMENT_PS(cg_element,i);
3216 for (j=0; j<ge_element[cge->ge].nSide; j++)
3217 Element_nb_uniq_subdom[i][j] = cge->nbid[j];
3218 }
3219 Enusdp[1] = Element_nb_uniq_subdom;
3220 theMesh.nbElements = Enusdp;
3221
3222 /* create levels */
3223 for (i=1; i<mg_general.nLevel; i++)
3224 if (CreateNewLevel(theMG)==NULL) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3225
3226 /* insert coarse mesh */
3227 if (InsertMesh(theMG,&theMesh)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3228 for (i=0; i<=TOPLEVEL(theMG); i++)
3229 for (theElement = PFIRSTELEMENT(GRID_ON_LEVEL(theMG,i)); theElement!=NULL; theElement=SUCCE(theElement))
3230 {
3231 cge = MGIO_CG_ELEMENT_PS(cg_element,ID(theElement));
3232 SETREFINE(theElement,0);
3233 SETREFINECLASS(theElement,NO_CLASS);
3234 SETMARK(theElement,0);
3235 SETMARKCLASS(theElement,NO_CLASS);
3236 SETSUBDOMAIN(theElement,cge->subdomain);
3237 for (j=0; j<CORNERS_OF_ELEM(theElement); j++)
3238 {
3239 theNode = CORNER(theElement,j);
3240 if (OBJT(MYVERTEX(theNode))==BVOBJ) SETNSUBDOM(theNode,0);
3241 else SETNSUBDOM(theNode,cge->subdomain);
3242 ID(CORNER(theElement,j)) = cge->cornerid[j];
3243 }
3244 for (j=0; j<EDGES_OF_ELEM(theElement); j++)
3245 {
3246 theEdge=GetEdge(CORNER_OF_EDGE_PTR(theElement,j,0),CORNER_OF_EDGE_PTR(theElement,j,1));
3247 ASSERT(theEdge!=NULL);
3248 if ((cge->se_on_bnd & (1<<(SIDES_OF_ELEM(theElement)+j))))
3249 SETEDSUBDOM(theEdge,0);
3250 else
3251 SETEDSUBDOM(theEdge,cge->subdomain);
3252 }
3253 }
3254
3255 /* now CreateAlgebra is necessary to have the coarse grid nodevectors for DDD identification in Evaluate_pinfo */
3256 if (CreateAlgebra (theMG)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3257 if (DisposeBottomHeapTmpMemory(theMG)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3258 if (PrepareAlgebraModification(theMG)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3259 #ifdef ModelP
3260 if (DisposeBottomHeapTmpMemory(theMG)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3261 #endif
3262
3263
3264 i = MG_ELEMUSED | MG_NODEUSED | MG_EDGEUSED | MG_VERTEXUSED | MG_VECTORUSED;
3265 ClearMultiGridUsedFlags(theMG,0,TOPLEVEL(theMG),i);
3266
3267 /* open identification context */
3268 DDD_IdentifyBegin(theMG->dddContext());
3269
3270 /* read parinfo of coarse-grid */
3271 if (MGIO_PARFILE)
3272 {
3273 ActProcListPos = ProcList = (unsigned short*)GetTmpMem(theHeap,PROCLISTSIZE*sizeof(unsigned short),MarkKey);
3274 if (ProcList==NULL) {UserWriteF("ERROR: cannot allocate %d bytes for ProcList\n",(int)PROCLISTSIZE*sizeof(unsigned short)); CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3275
3276 cg_pinfo.proclist = ProcList;
3277 for (level=0; level<=TOPLEVEL(theMG); level++)
3278 {
3279 theGrid = GRID_ON_LEVEL(theMG,level);
3280 for (theElement = PFIRSTELEMENT(theGrid); theElement!=NULL; theElement=ENext)
3281 {
3282 /* store succe because Evaluate_pinfo() relinks the element list */
3283 ENext = SUCCE(theElement);
3284 cge = MGIO_CG_ELEMENT_PS(cg_element,ID(theElement));
3285 if (Read_pinfo (TAG(theElement),&cg_pinfo)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3286 if (Evaluate_pinfo(theGrid,theElement,&cg_pinfo)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3287 PRINTDEBUG(gm,1,("read orphan element " EID_FMTX "\n", EID_PRTX(theElement)));
3288 }
3289 }
3290 }
3291
3292 /* are we ready ? */
3293 if (mg_general.nLevel==1)
3294 { /* if we are here now, all other processors (wich haven't returned yet) are here TOO
3295 because mg_general.nLevel is a global quantity. */
3296
3297 /* no fine grid elements */
3298 if (MGCreateConnection(theMG)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3299
3300 /* close identification context */
3301 DDD_IdentifyEnd(theMG->dddContext());
3302
3303 /* repair inconsistencies */
3304 if (MGIO_PARFILE)
3305 if (IO_GridCons(theMG)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3306
3307 for (theVector=PFIRSTVECTOR(GRID_ON_LEVEL(theMG,0));
3308 theVector!=NULL; theVector=SUCCVC(theVector)) {
3309 SETVCLASS(theVector,3);
3310 SETVNCLASS(theVector,0);
3311 SETNEW_DEFECT(theVector,1);
3312 SETFINE_GRID_DOF(theVector,1);
3313 }
3314
3315 if (SetSurfaceClasses (theMG)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3316
3317 ReleaseTmpMem(theHeap,MarkKey);
3318 if (CloseMGFile ()) {DisposeMultiGrid(theMG); return (NULL);}
3319
3320 /* saved */
3321 MG_SAVED(theMG) = 1;
3322 strcpy(MG_FILENAME(theMG),filename);
3323
3324 return (theMG);
3325 }
3326
3327 /* list: node-id --> node */
3328 nid_n = (NODE**)GetTmpMem(theHeap,non*sizeof(NODE*),MarkKey);
3329 for (i=0; i<=TOPLEVEL(theMG); i++)
3330 for (theNode=PFIRSTNODE(GRID_ON_LEVEL(theMG,i)); theNode!=NULL; theNode=SUCCN(theNode))
3331 {
3332 assert(ID(theNode)-foid<non);
3333 nid_n[ID(theNode)-foid] = theNode;
3334 }
3335
3336 /* list: elem-id --> elem */
3337 eid_e = (ELEMENT**)GetTmpMem(theHeap,cg_general.nElement*sizeof(ELEMENT*),MarkKey);
3338 for (i=0; i<=TOPLEVEL(theMG); i++)
3339 for (theElement=PFIRSTELEMENT(GRID_ON_LEVEL(theMG,i)); theElement!=NULL; theElement=SUCCE(theElement))
3340 {
3341 assert(ID(theElement)>=0);
3342 assert(ID(theElement)<cg_general.nElement);
3343 eid_e[ID(theElement)] = theElement;
3344 }
3345
3346 if (CheckCGKeys(cg_general.nElement, eid_e, cg_element)==0)
3347 {
3348 PRINTDEBUG(gm,4,("CheckCGKeys ok.\n"));
3349 }
3350 else
3351 {
3352 PrintErrorMessage('E',"LoadMultiGrid","CheckCGKeys failed\n");
3353 }
3354
3355 /* read hierarchical elements */
3356 refinement = (MGIO_REFINEMENT*)malloc(MGIO_REFINEMENT_SIZE);
3357 if (refinement==NULL) {UserWriteF("ERROR: cannot allocate %d bytes for refinement\n",(int)MGIO_REFINEMENT_SIZE); CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3358 if (MGIO_PARFILE)
3359 {
3360 ProcList = (unsigned short*)malloc(PROCLISTSIZE*sizeof(unsigned short));
3361 if (ProcList==NULL) {UserWriteF("ERROR: cannot allocate %d bytes for ProcList\n",(int)PROCLISTSIZE*sizeof(unsigned short)); CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3362 for (i=0; i<MAX_SONS; i++) refinement->pinfo[i].proclist = ProcList+i*ELEMPROCLISTSIZE;
3363 }
3364 for (j=0; j<cg_general.nElement; j++)
3365 {
3366 theElement = eid_e[j];
3367 theGrid = GRID_ON_LEVEL(theMG,LEVEL(theElement));
3368 cge = MGIO_CG_ELEMENT_PS(cg_element,j);
3369 if (cge->nref==0)
3370 {
3371 SETREFINE(theElement,NO_REFINEMENT);
3372 SETREFINECLASS(theElement,NO_CLASS);
3373 SETMARK(theElement,NO_REFINEMENT);
3374 SETMARKCLASS(theElement,NO_CLASS);
3375 if (LEVEL(theElement)==0) SETECLASS(theElement,RED_CLASS);
3376 #ifdef ModelP
3377 else assert(EGHOST(theElement)); /* masters elements must have a father or be on level 0*/
3378 #else
3379 else assert(0); /* all orphans must be on level 0 in sequential mode */
3380 #endif
3381 continue;
3382 }
3383 if (InsertLocalTree(theGrid,theElement,refinement,rr_general.RefRuleOffset)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3384 }
3385
3386 /* close identification context */
3387 #ifdef OPTIMIZED_IO
3388 DDD_SetOption(theMG->dddContext(), OPT_IF_CREATE_EXPLICIT,OPT_ON);
3389 #endif
3390 DDD_IdentifyEnd(theMG->dddContext());
3391 #ifdef OPTIMIZED_IO
3392 DDD_SetOption(theMG->dddContext(), OPT_IF_CREATE_EXPLICIT,OPT_OFF);
3393 #endif
3394
3395 /* repair inconsistencies */
3396 if (MGIO_PARFILE)
3397 {
3398 #ifdef ModelP
3399 CommunicateEClasses(theMG);
3400 #endif
3401 if (IO_GridCons(theMG)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3402 }
3403
3404 /* postprocess */
3405 for (i=0; i<TOPLEVEL(theMG); i++)
3406 for (theElement = FIRSTELEMENT(GRID_ON_LEVEL(theMG,i)); theElement!=NULL; theElement=SUCCE(theElement))
3407 {
3408 SETMARK(theElement,0);
3409 SETMARKCLASS(theElement,NO_CLASS);
3410 SETEBUILDCON(theElement,1);
3411 }
3412 for (i=0; i<=TOPLEVEL(theMG); i++)
3413 {
3414 theGrid = GRID_ON_LEVEL(theMG,i);
3415
3416 #ifdef __THREEDIM__
3417 if (VEC_DEF_IN_OBJ_OF_MG(MYMG(theGrid),SIDEVEC))
3418 for (theElement = FIRSTELEMENT(theGrid); theElement!=NULL; theElement=SUCCE(theElement))
3419 for (j=0; j<SIDES_OF_ELEM(theElement); j++)
3420 {
3421 theNeighbor = NBELEM(theElement,j);
3422 if (theNeighbor==NULL) continue;
3423 for (k=0; k<SIDES_OF_ELEM(theNeighbor); k++)
3424 if (NBELEM(theNeighbor,k)==theElement)
3425 break;
3426 if (k==SIDES_OF_ELEM(theNeighbor)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3427 if (DisposeDoubledSideVector (theGrid,theElement,j,theNeighbor,k)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3428 }
3429 #endif
3430 }
3431 if (MGCreateConnection(theMG)) {CloseMGFile (); DisposeMultiGrid(theMG); return (NULL);}
3432 theGrid = GRID_ON_LEVEL(theMG,0);
3433 ClearNextNodeClasses(theGrid);
3434 for (theElement=FIRSTELEMENT(theGrid);
3435 theElement!=NULL; theElement=SUCCE(theElement))
3436 if (REFINECLASS(theElement)>=GREEN_CLASS)
3437 SeedNextNodeClasses(theElement);
3438 PropagateNextNodeClasses(theGrid);
3439 for (i=1; i<TOPLEVEL(theMG); i++) {
3440 theGrid = GRID_ON_LEVEL(theMG,i);
3441 ClearNodeClasses(theGrid);
3442 ClearNextNodeClasses(theGrid);
3443 for (theElement=FIRSTELEMENT(theGrid);
3444 theElement!=NULL; theElement=SUCCE(theElement)) {
3445 if (ECLASS(theElement)>=GREEN_CLASS)
3446 SeedNodeClasses(theElement);
3447 if (REFINECLASS(theElement)>=GREEN_CLASS)
3448 SeedNextNodeClasses(theElement);
3449 }
3450 PropagateNodeClasses(theGrid);
3451 PropagateNextNodeClasses(theGrid);
3452 }
3453 theGrid = GRID_ON_LEVEL(theMG,TOPLEVEL(theMG));
3454 ClearNodeClasses(theGrid);
3455 ClearNextNodeClasses(theGrid);
3456 for (theElement=FIRSTELEMENT(theGrid);
3457 theElement!=NULL; theElement=SUCCE(theElement))
3458 if (ECLASS(theElement)>=GREEN_CLASS)
3459 SeedNodeClasses(theElement);
3460 PropagateNodeClasses(theGrid);
3461 if (PrepareAlgebraModification(theMG)) {DisposeMultiGrid(theMG); return (NULL);}
3462
3463 /* set DOFs on vectors */
3464 if (SetSurfaceClasses (theMG)) {DisposeMultiGrid(theMG); return (NULL);}
3465
3466 /* close file */
3467 ReleaseTmpMem(theHeap,MarkKey);
3468 if (CloseMGFile ()) {DisposeMultiGrid(theMG); return (NULL);}
3469
3470
3471 /* saved */
3472 MG_SAVED(theMG) = 1;
3473 strcpy(MG_FILENAME(theMG),filename);
3474
3475 #if EXTRACT_RULES
3476 ResetRefineTagsBeyondRuleManager(theMG);
3477 #endif
3478
3479 return (theMG);
3480 }
3481
InitUgio()3482 INT NS_DIM_PREFIX InitUgio ()
3483 {
3484 /* read gridpaths from defaults file (iff) */
3485 gridpaths_set = false;
3486 if (ReadSearchingPaths("defaults","gridpaths")==0)
3487 gridpaths_set = true;
3488
3489 if (MGIO_Init ()) return(1);
3490
3491 return (0);
3492 }
3493