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: ugm.c */
6 /* */
7 /* Purpose: unstructured grid manager */
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: 09.03.92 begin, ug version 2.0 */
17 /* Aug 28 1996, ug version 3.4 */
18 /* */
19 /* Remarks: */
20 /* */
21 /****************************************************************************/
22
23 /****************************************************************************/
24 /* */
25 /* defines to exclude functions */
26 /* */
27 /****************************************************************************/
28
29 /****************************************************************************/
30 /* */
31 /* include files */
32 /* system include files */
33 /* application include files */
34 /* */
35 /****************************************************************************/
36
37 #include <config.h>
38 #include <cstdlib>
39 #include <cstddef>
40 #include <cstdio>
41 #include <cstring>
42 #include <cmath>
43 #include <cassert>
44 #include <algorithm>
45
46 #include <errno.h>
47 #include <vector>
48
49 #include <dune/uggrid/low/architecture.h>
50 #include <dune/uggrid/low/heaps.h>
51 #include <dune/uggrid/low/debug.h>
52 #include <dune/uggrid/low/fifo.h>
53 #include <dune/uggrid/low/misc.h>
54 #include <dune/uggrid/low/ugenv.h>
55 #include <dune/uggrid/low/ugstruct.h>
56 #include <dune/uggrid/low/ugtypes.h>
57
58 #include <dune/uggrid/ugdevices.h>
59
60 #include "cw.h"
61 #include "evm.h"
62 #include "gm.h"
63 #include "rm.h"
64 #include "dlmgr.h"
65 #include "algebra.h"
66 #include "ugm.h"
67 #include "elements.h"
68 #include "shapes.h"
69 #include "refine.h"
70 #include <dune/uggrid/domain/domain.h>
71 #include "pargm.h"
72 #include "mgheapmgr.h"
73
74 #ifdef ModelP
75 #include <dune/uggrid/parallel/dddif/identify.h>
76 #include <dune/uggrid/parallel/ppif/ppif.h>
77 #endif
78
79 #include <dune/uggrid/parallel/ppif/ppifcontext.hh>
80
81 USING_UG_NAMESPACE
82 USING_UGDIM_NAMESPACE
83 using namespace PPIF;
84
85 /****************************************************************************/
86 /* */
87 /* defines in the following order */
88 /* */
89 /* compile time constants defining static data size (i.e. arrays) */
90 /* other constants */
91 /* macros */
92 /* */
93 /****************************************************************************/
94
95 #define RESOLUTION 20 /* resolution for creating boundary midnode */
96 #define SMALL1 0.001
97
98 #define LINKTABLESIZE 32 /* max number of inks per node for ordering */
99
100 /** \brief macro for controlling debugging output by conditions on objects */
101 #define UGM_CDBG(x,y)
102
103 /****************************************************************************/
104 /* */
105 /* data structures used in this source file (exported data structures are */
106 /* in the corresponding include file!) */
107 /* */
108 /****************************************************************************/
109
110 /****************************************************************************/
111 /* */
112 /* definition of exported global variables */
113 /* */
114 /****************************************************************************/
115
116 #if defined ModelP && defined __OVERLAP2__
117 INT ce_NO_DELETE_OVERLAP2 = -1;
118 #endif
119
120
121 /****************************************************************************/
122 /* */
123 /* definition of variables global to this source file only (static!) */
124 /* */
125 /****************************************************************************/
126
127 /** \brief General purpose text buffer */
128 static char buffer[4*256];
129
130 static INT theMGDirID; /* env var ID for the multigrids */
131 static INT theMGRootDirID; /* env dir ID for the multigrids */
132
133 static UINT UsedOBJT; /* for the dynamic OBJECT management */
134
135 /****************************************************************************/
136 /* */
137 /* forward declarations of functions used before they are defined */
138 /* */
139 /****************************************************************************/
140
141 static INT DisposeVertex (GRID *theGrid, VERTEX *theVertex);
142 static INT DisposeEdge (GRID *theGrid, EDGE *theEdge);
143
144
145 /****************************************************************************/
146 /** \brief Get an object type id not occupied in theMG
147 *
148 * This function gets an object type id not occupied in theMG.
149 *
150 * @return <ul>
151 * <li> id of object type if ok </li>
152 * <li> -1 when error occurred </li>
153 * </ul>
154 */
155 /****************************************************************************/
156
GetFreeOBJT()157 GM_OBJECTS NS_DIM_PREFIX GetFreeOBJT ()
158 {
159 INT i;
160
161 /* skip predefined object types, they cannot be re-allocated */
162 for (i=NPREDEFOBJ; i<MAXOBJECTS; i++)
163 if (!READ_FLAG(UsedOBJT,1<<i))
164 break;
165
166 if (i<MAXOBJECTS)
167 {
168 SET_FLAG(UsedOBJT,1<<i);
169 return (GM_OBJECTS)i;
170 }
171 else
172 return NOOBJ;
173 }
174
175 /****************************************************************************/
176 /** \brief Release an object type id not needed anymore
177 *
178 * @param type - object type
179 *
180 * This function releases an object type id not needed anymore.
181 *
182 * @return <ul>
183 * <li> GM_OK if ok </li>
184 * <li> GM_ERROR when error occured. </li>
185 * </ul>
186 */
187 /****************************************************************************/
188
ReleaseOBJT(GM_OBJECTS type)189 INT NS_DIM_PREFIX ReleaseOBJT (GM_OBJECTS type)
190 {
191 if (type>=MAXOBJECTS)
192 RETURN (GM_ERROR);
193
194 /* we cannot release predefined object types! */
195 if (type<NPREDEFOBJ)
196 RETURN (GM_ERROR);
197
198 CLEAR_FLAG(UsedOBJT,1<<type);
199
200 return (GM_OK);
201 }
202
203 /****************************************************************************/
204 /** \brief Get an object from free list if possible
205 \fn GetMemoryForObject
206
207 * @param theMG - pointer to multigrid
208 * @param size - size of the object
209 * @param type - type of the requested object
210
211 This function gets an object of type `type` from free list if possible,
212 otherwise it allocates memory from the multigrid heap using 'GetMem'.
213
214 @return <ul>
215 <li> pointer to an object of the requested type </li>
216 <li> NULL if object of requested type is not available </li>
217 * </ul>
218 */
219 /****************************************************************************/
220
221 #ifdef ModelP
ConstructDDDObject(DDD::DDDContext & context,void * obj,INT size,INT type)222 static void ConstructDDDObject (DDD::DDDContext& context, void *obj, INT size, INT type)
223 {
224 if (obj!=NULL && type!=NOOBJ)
225 {
226 memset(obj,0,size);
227 /* link this object to DDD management */
228 if (HAS_DDDHDR(context, type))
229 {
230 DDD_TYPE dddtype = DDDTYPE(context, type);
231 DDD_HDR dddhdr = (DDD_HDR)(((char *)obj) + DDD_InfoHdrOffset(context, dddtype));
232 DDD_HdrConstructor(context, dddhdr, dddtype, PrioMaster, 0);
233 }
234 }
235 return;
236 }
237 #endif
238
GetMemoryForObject(MULTIGRID * theMG,INT size,INT type)239 void * NS_DIM_PREFIX GetMemoryForObject (MULTIGRID *theMG, INT size, INT type)
240 {
241 void * obj = GetMem(MGHEAP(theMG),size);
242 if (obj != NULL)
243 memset(obj,0,size);
244
245 #ifdef ModelP
246 if (type!=MAOBJ && type!=COOBJ)
247 ConstructDDDObject(theMG->dddContext(), obj,size,type);
248 #endif
249
250 return obj;
251 }
252
253 /****************************************************************************/
254 /** \brief Put an object in the free list
255 \fn PutFreeObject
256
257 * @param theMG - pointer to multigrid
258 * @param object - object to insert in free list
259 * @param size - size of the object
260 * @param type - type of the requested object
261
262 This function puts an object in the free list.
263
264 @return <ul>
265 <li> 0 if ok </li>
266 <li> 1 when error occured. </li>
267 * </ul> */
268 /****************************************************************************/
269
270 #ifdef ModelP
DestructDDDObject(DDD::DDDContext & context,void * object,INT type)271 static void DestructDDDObject(DDD::DDDContext& context, void *object, INT type)
272 {
273 if (type!=NOOBJ)
274 {
275 /* unlink object from DDD management */
276 if (HAS_DDDHDR(context, type))
277 {
278 DDD_HDR dddhdr = (DDD_HDR)
279 (((char *)object)+DDD_InfoHdrOffset(context, DDDTYPE(context, type)));
280 DDD_HdrDestructor(context, dddhdr);
281 }
282 }
283 return;
284 }
285 #endif
286
PutFreeObject(MULTIGRID * theMG,void * object,INT size,GM_OBJECTS type)287 INT NS_DIM_PREFIX PutFreeObject (MULTIGRID *theMG, void *object, INT size, GM_OBJECTS type)
288 {
289 #ifdef ModelP
290 if (type!=MAOBJ && type!=COOBJ)
291 DestructDDDObject(theMG->dddContext(), object,type);
292 #endif
293
294 DisposeMem(MGHEAP(theMG), object);
295 return 0;
296 }
297
298 /****************************************************************************/
299 /** \brief Return pointer to a new boundary vertex structure
300 *
301 * @param theGrid grid where vertex should be inserted
302 *
303 * This function creates and initializes a new boundary vertex structure
304 * and returns a pointer to it.
305 *
306 * @return <ul>
307 * <li> pointer to requested object </li>
308 * <li> NULL if out of memory </li>
309 * </ul>
310 */
311 /****************************************************************************/
312
CreateBoundaryVertex(GRID * theGrid)313 static VERTEX *CreateBoundaryVertex (GRID *theGrid)
314 {
315 VERTEX *pv;
316 INT i;
317
318 pv = (VERTEX*)GetMemoryForObject(MYMG(theGrid),sizeof(struct bvertex),BVOBJ);
319 if (pv==NULL) return(NULL);
320 VDATA(pv) = NULL;
321
322 /* initialize data */
323 CTRL(pv) = 0;
324 SETOBJT(pv,BVOBJ);
325 SETNOOFNODE(pv,0);
326 SETLEVEL(pv,theGrid->level);
327 ID(pv) = (theGrid->mg->vertIdCounter)++;
328 VFATHER(pv) = NULL;
329 #ifdef TOPNODE
330 TOPNODE(pv) = NULL;
331 #endif
332 for (i=0; i<DIM; i++) LCVECT(pv)[i] = 0.0;
333 SETONEDGE(pv,0);
334 SETMOVE(pv,DIM_OF_BND);
335 #ifdef ModelP
336 DDD_AttrSet(PARHDRV(pv),GRID_ATTR(theGrid));
337 /* SETVXPRIO(pv,PrioMaster); */
338 #endif
339
340 /* insert in vertex list */
341 GRID_LINK_VERTEX(theGrid,pv,PrioMaster);
342
343 return(pv);
344 }
345
346 /****************************************************************************/
347 /** \brief Return pointer to a new inner vertex structure
348 *
349 * @param theGrid grid where vertex should be inserted
350 *
351 * This function creates and initializes a new inner vertex structure
352 * and returns a pointer to it.
353 *
354 * @return <ul>
355 <li> pointer to requested object </li>
356 <li> NULL if out of memory </li>
357 </ul> */
358 /****************************************************************************/
359
CreateInnerVertex(GRID * theGrid)360 static VERTEX *CreateInnerVertex (GRID *theGrid)
361 {
362 VERTEX *pv;
363 INT i;
364
365 pv = (VERTEX*)GetMemoryForObject(MYMG(theGrid),sizeof(struct ivertex),IVOBJ);
366 if (pv==NULL) return(NULL);
367 VDATA(pv) = NULL;
368
369 /* initialize data */
370 CTRL(pv) = 0;
371 SETOBJT(pv,IVOBJ);
372 SETNOOFNODE(pv,0);
373 SETLEVEL(pv,theGrid->level);
374 ID(pv) = (theGrid->mg->vertIdCounter)++;
375 VFATHER(pv) = NULL;
376 #ifdef TOPNODE
377 TOPNODE(pv) = NULL;
378 #endif
379 SETMOVE(pv,DIM);
380 #ifdef ModelP
381 DDD_AttrSet(PARHDRV(pv),GRID_ATTR(theGrid));
382 /* SETVXPRIO(pv,PrioMaster); */
383 #endif
384 for (i=0; i<DIM; i++) LCVECT(pv)[i] = 0.0;
385
386 /* insert in vertex list */
387 GRID_LINK_VERTEX(theGrid,pv,PrioMaster);
388
389 return(pv);
390 }
391
392 /****************************************************************************/
393 /** \brief Return pointer to a new node structure
394
395 * @param theGrid - grid where vertex should be inserted
396 * @param vertex - vertex of the node
397 * @param FatherNode - father node (may be NULL)
398 * @param NodeType - node type (CORNER_NODE..., cf. gm.h)
399 * @param with_vector
400
401 This function creates and initializes a new node structure
402 and returns a pointer to it.
403
404 @return <ul>
405 <li> pointer to requested object </li>
406 <li> NULL if out of memory </li>
407 </ul> */
408 /****************************************************************************/
409
CreateNode(GRID * theGrid,VERTEX * vertex,GEOM_OBJECT * Father,INT NodeType,INT with_vector)410 static NODE *CreateNode (GRID *theGrid, VERTEX *vertex,
411 GEOM_OBJECT *Father, INT NodeType, INT with_vector)
412 {
413 NODE *pn;
414 VECTOR *pv;
415 INT size;
416
417 size = sizeof(NODE);
418 if (!VEC_DEF_IN_OBJ_OF_GRID(theGrid,NODEVEC)) size -= sizeof(VECTOR *);
419
420 pn = (NODE *)GetMemoryForObject(MYMG(theGrid),size,NDOBJ);
421 if (pn==NULL) return(NULL);
422
423 /* initialize data */
424 SETOBJT(pn,NDOBJ);
425 SETLEVEL(pn,theGrid->level);
426 #ifdef ModelP
427 DDD_AttrSet(PARHDR(pn),GRID_ATTR(theGrid));
428 /* SETPRIO(pn,PrioMaster); */
429 pn->message_buffer_ = nullptr;
430 pn->message_buffer_size_ = 0;
431 #endif
432 ID(pn) = (theGrid->mg->nodeIdCounter)++;
433 START(pn) = NULL;
434 SONNODE(pn) = NULL;
435 MYVERTEX(pn) = vertex;
436 if (NOOFNODE(vertex)<NOOFNODEMAX)
437 INCNOOFNODE(vertex);
438 else
439 ASSERT(0);
440 /* priliminary */
441 if (Father != NULL)
442 if ((OBJT(Father) == IEOBJ) || (OBJT(Father) == BEOBJ))
443 Father = NULL;
444 SETNFATHER(pn,Father);
445 SETNTYPE(pn,NodeType);
446 SETNCLASS(pn,3),
447 SETNNCLASS(pn,0);
448 if (OBJT(vertex) == BVOBJ)
449 SETNSUBDOM(pn,0);
450 else if (VFATHER(vertex) != NULL)
451 SETNSUBDOM(pn,SUBDOMAIN(VFATHER(vertex)));
452 else if (Father != NULL) {
453 if (OBJT(Father) == NDOBJ)
454 SETNSUBDOM(pn,NSUBDOM((NODE *)Father));
455 else if (OBJT(Father) == EDOBJ)
456 SETNSUBDOM(pn,EDSUBDOM((EDGE *)Father));
457 }
458 else
459 SETNSUBDOM(pn,0);
460
461 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,NODEVEC))
462 {
463 if (with_vector)
464 {
465 if (CreateVector (theGrid,NODEVEC,(GEOM_OBJECT *)pn,&pv))
466 {
467 DisposeNode (theGrid,pn);
468 return (NULL);
469 }
470 NVECTOR(pn) = pv;
471 }
472 else
473 NVECTOR(pn) = NULL;
474 }
475
476 theGrid->status |= 1; /* recalculate stiffness matrix */
477
478 /* insert in vertex list */
479 GRID_LINK_NODE(theGrid,pn,PrioMaster);
480
481 return(pn);
482 }
483
484 /****************************************************************************/
485 /** \brief Return pointer to a new node structure on an edge
486
487 * @param theGrid - grid where vertex should be inserted
488 * @param FatherNode - node father
489
490 This function creates and initializes a new node structure
491 at the midpoint of an element edge and returns a pointer to it.
492
493 @return <ul>
494 <li> pointer to requested object </li>
495 <li> NULL if out of memory </li>
496 </ul> */
497 /****************************************************************************/
498
CreateSonNode(GRID * theGrid,NODE * FatherNode)499 NODE *NS_DIM_PREFIX CreateSonNode (GRID *theGrid, NODE *FatherNode)
500 {
501 NODE *pn;
502 VERTEX *theVertex;
503
504 theVertex = MYVERTEX(FatherNode);
505
506 pn = CreateNode(theGrid,theVertex,(GEOM_OBJECT *)FatherNode,CORNER_NODE,1);
507 if (pn == NULL)
508 return(NULL);
509 SONNODE(FatherNode) = pn;
510
511 #ifdef TOPNODE
512 TOPNODE(theVertex) = pn;
513 #endif
514
515 return(pn);
516 }
517
518 /****************************************************************************/
519 /** \brief Return pointer to a new node structure on an edge
520
521 * @param theGrid - grid where node should be inserted
522 * @param theElement - pointer to an element
523 * @param theVertex - pointer to vertex if already existing
524 * @param edge - id of an element edge
525
526 This function creates and initializes a new node structure
527 at the midpoint of an element edge and returns a pointer to it.
528
529 @return <ul>
530 <li> pointer to requested object </li>
531 <li> NULL if out of memory </li>
532 </ul> */
533 /****************************************************************************/
534
CreateMidNode(GRID * theGrid,ELEMENT * theElement,VERTEX * theVertex,INT edge)535 NODE *NS_DIM_PREFIX CreateMidNode (GRID *theGrid, ELEMENT *theElement, VERTEX *theVertex, INT edge)
536 {
537 BNDP *bndp;
538 DOUBLE *local,*x[MAX_CORNERS_OF_ELEM];
539 DOUBLE_VECTOR bnd_global,global;
540 DOUBLE diff;
541 INT n,move,part;
542
543 const INT co0 = CORNER_OF_EDGE(theElement,edge,0);
544 const INT co1 = CORNER_OF_EDGE(theElement,edge,1);
545 VERTEX* v0 = MYVERTEX(CORNER(theElement,co0));
546 VERTEX* v1 = MYVERTEX(CORNER(theElement,co1));
547 V_DIM_LINCOMB(0.5, CVECT(v0), 0.5, CVECT(v1), global);
548
549 /* set MIDNODE pointer */
550 EDGE* theEdge = GetEdge(CORNER(theElement,co0),CORNER(theElement,co1));
551 ASSERT(theEdge!=NULL);
552
553 /* allocate vertex */
554 const bool vertex_null = (theVertex==NULL);
555 if (theVertex==NULL)
556 {
557 if ((OBJT(v0) == BVOBJ) && (OBJT(v1) == BVOBJ))
558 #ifdef __TWODIM__
559 if (OBJT(theElement) == BEOBJ)
560 if (SIDE_ON_BND(theElement,edge))
561 #endif
562 #ifdef __THREEDIM__
563 if (EDSUBDOM(theEdge) == 0)
564 #endif
565 {
566 bndp = BNDP_CreateBndP(MGHEAP(MYMG(theGrid)),V_BNDP(v0),V_BNDP(v1),0.5);
567 if (bndp != NULL)
568 {
569 theVertex = CreateBoundaryVertex(theGrid);
570 if (theVertex == NULL)
571 return(NULL);
572 if (BNDP_Global(bndp,bnd_global))
573 return(NULL);
574 if (BNDP_BndPDesc(bndp,&move,&part))
575 return(NULL);
576 SETMOVE(theVertex,move);
577 V_BNDP(theVertex) = bndp;
578 V_DIM_COPY(bnd_global,CVECT(theVertex));
579 local = LCVECT(theVertex);
580 V_DIM_EUKLIDNORM_OF_DIFF(bnd_global,global,diff);
581 if (diff > MAX_PAR_DIST)
582 {
583 SETMOVED(theVertex,1);
584 CORNER_COORDINATES(theElement,n,x);
585 UG_GlobalToLocal(n,(const DOUBLE **)x,bnd_global,local);
586 }
587 else
588 V_DIM_LINCOMB(0.5, LOCAL_COORD_OF_ELEM(theElement,co0),
589 0.5, LOCAL_COORD_OF_ELEM(theElement,co1),local);
590 PRINTDEBUG(gm,1,("local = %f %f %f\n",local[0],local[1],local[2]));
591 }
592 }
593 if (theVertex == NULL)
594 {
595 /* we need an inner vertex */
596 theVertex = CreateInnerVertex(theGrid);
597 if (theVertex==NULL) return(NULL);
598 V_DIM_COPY(global,CVECT(theVertex));
599 V_DIM_LINCOMB(0.5, LOCAL_COORD_OF_ELEM(theElement,co0),
600 0.5, LOCAL_COORD_OF_ELEM(theElement,co1),
601 LCVECT(theVertex));
602 }
603 VFATHER(theVertex) = theElement;
604 SETONEDGE(theVertex,edge);
605 }
606 /* allocate node */
607 NODE* theNode = CreateNode(theGrid,theVertex,(GEOM_OBJECT *)theEdge,MID_NODE,1);
608 if (theNode==NULL && vertex_null)
609 {
610 DisposeVertex(theGrid,theVertex);
611 return(NULL);
612 }
613
614 MIDNODE(theEdge) = theNode;
615 #ifdef TOPNODE
616 if (TOPNODE(theVertex)==NULL || LEVEL(TOPNODE(theVertex))<LEVEL(theNode))
617 TOPNODE(theVertex) = theNode;
618 #endif
619
620 if (OBJT(theVertex) == BVOBJ)
621 PRINTDEBUG(dom,1,(" MidPoint %d %f %f %f\n",ID(theNode),
622 bnd_global[0],
623 bnd_global[1],
624 bnd_global[2]));
625
626 PRINTDEBUG(dddif,1,(PFMT " CreateMidNode(): n=" ID_FMTX
627 " NTYPE=%d OBJT=%d father " ID_FMTX " \n",
628 theGrid->ppifContext().me(),ID_PRTX(theNode),NTYPE(theNode),
629 OBJT(NFATHER(theNode)),NFATHER(theNode)));
630
631 return(theNode);
632 }
633
634
GetMidNode(const ELEMENT * theElement,INT edge)635 NODE * NS_DIM_PREFIX GetMidNode (const ELEMENT *theElement, INT edge)
636 {
637 EDGE *theEdge;
638 NODE *theNode;
639 VERTEX *theVertex;
640
641 theEdge = GetEdge(CORNER(theElement,CORNER_OF_EDGE(theElement,edge,0)),
642 CORNER(theElement,CORNER_OF_EDGE(theElement,edge,1)));
643 if (theEdge == NULL) return(NULL);
644 theNode = MIDNODE(theEdge);
645 if (theNode == NULL) return(NULL);
646
647 /** \todo This is a bad place for the following code (s.l. 981015) */
648 theVertex = MYVERTEX(theNode);
649 if (theVertex!=NULL && VFATHER(theVertex) == NULL) {
650 /** \todo Strange that this cast has to be here. O.S. 060902 */
651 VFATHER(theVertex) = (ELEMENT*)theElement;
652 SETONEDGE(theVertex,edge);
653 V_DIM_LINCOMB(0.5,
654 LOCAL_COORD_OF_ELEM(theElement,
655 CORNER_OF_EDGE(theElement,edge,0)),
656 0.5,
657 LOCAL_COORD_OF_ELEM(theElement,
658 CORNER_OF_EDGE(theElement,edge,1)),
659 LCVECT(theVertex));
660 }
661 return(theNode);
662 }
663
664
665 /****************************************************************************/
666 /** \brief ???
667 */
668 /****************************************************************************/
669
SideOfNbElement(const ELEMENT * theElement,INT side)670 static INT SideOfNbElement(const ELEMENT *theElement, INT side)
671 {
672 ELEMENT *nb;
673 NODE *nd[MAX_CORNERS_OF_SIDE];
674 INT i,j,m,n,num;
675
676 nb = NBELEM(theElement,side);
677 if (nb == NULL) return(MAX_SIDES_OF_ELEM);
678
679 for (j=0; j<SIDES_OF_ELEM(nb); j++)
680 if (NBELEM(nb,j) == theElement) return(j);
681
682 n = CORNERS_OF_SIDE(theElement,side);
683 for (i=0; i<n; i++)
684 nd[i] = CORNER(theElement,CORNER_OF_SIDE(theElement,side,i));
685
686 for (j=0; j<SIDES_OF_ELEM(nb); j++) {
687 num = 0;
688 for (i=0; i<n; i++)
689 for (m=0; m<CORNERS_OF_SIDE(nb,j); m++)
690 if (nd[i] == CORNER(nb,CORNER_OF_SIDE(nb,j,m))) num++;
691 if (num == n) return(j);
692 }
693
694 return(MAX_SIDES_OF_ELEM);
695 }
696
697 #ifdef __THREEDIM__
698
699
700 /****************************************************************************/
701 /** \brief Return pointer to a new node structure on a side (3d)
702
703 * @param theGrid - grid where vertex should be inserted
704 * @param theElement - pointer to an element
705 * @param theVertex - pointer vertex
706 * @param side - id of an element side
707
708 This function creates and initializes a new node structure
709 at the midpoint of an element side and returns a pointer to it.
710
711 @return <ul>
712 <li> pointer to requested object </li>
713 <li> NULL if out of memory </li>
714 </ul> */
715 /****************************************************************************/
716
CreateSideNode(GRID * theGrid,ELEMENT * theElement,VERTEX * theVertex,INT side)717 NODE *NS_DIM_PREFIX CreateSideNode (GRID *theGrid, ELEMENT *theElement, VERTEX *theVertex, INT side)
718 {
719 DOUBLE_VECTOR bnd_global,global,local,bnd_local;
720 DOUBLE *x[MAX_CORNERS_OF_ELEM];
721 NODE *theNode;
722 BNDP *bndp;
723 BNDS *bnds;
724 DOUBLE fac, diff;
725 INT n,j,k,move,part,vertex_null;
726
727 n = CORNERS_OF_SIDE(theElement,side);
728 fac = 1.0 / n;
729 V_DIM_CLEAR(local);
730 V_DIM_CLEAR(global);
731 for (j=0; j<n; j++)
732 {
733 k = CORNER_OF_SIDE(theElement,side,j);
734 V_DIM_LINCOMB(1.0,local,1.0,
735 LOCAL_COORD_OF_ELEM(theElement,k),local);
736 V_DIM_LINCOMB(1.0,global,1.0,
737 CVECT(MYVERTEX(CORNER(theElement,k))),global);
738 }
739 V_DIM_SCALE(fac,local);
740 V_DIM_SCALE(fac,global);
741
742 /* check if boundary vertex */
743 vertex_null = (theVertex==NULL);
744 if (theVertex==NULL)
745 {
746 if (OBJT(theElement) == BEOBJ) {
747 bnds = ELEM_BNDS(theElement,side);
748 if (bnds != NULL) {
749 if (n == 3)
750 bnd_local[0] = bnd_local[1] = 0.33333333333333;
751 else if (n == 4)
752 bnd_local[0] = bnd_local[1] = 0.5;
753 bndp = BNDS_CreateBndP(MGHEAP(MYMG(theGrid)),bnds,bnd_local);
754 if (bndp != NULL) {
755 theVertex = CreateBoundaryVertex(theGrid);
756 if (theVertex == NULL)
757 return(NULL);
758 if (BNDP_BndPDesc(bndp,&move,&part))
759 return(NULL);
760 SETMOVE(theVertex,move);
761 if (BNDP_Global(bndp,bnd_global))
762 return(NULL);
763 V_BNDP(theVertex) = bndp;
764 V_DIM_COPY(bnd_global,CVECT(theVertex));
765 V_DIM_EUKLIDNORM_OF_DIFF(bnd_global,global,diff);
766 if (diff > MAX_PAR_DIST) {
767 SETMOVED(theVertex,1);
768 CORNER_COORDINATES(theElement,k,x);
769 UG_GlobalToLocal(k,(const DOUBLE **)x,bnd_global,local);
770 PRINTDEBUG(gm,1,("local = %f %f %f\n",local[0],local[1],local[2]));
771 }
772 }
773 }
774 }
775
776 if (theVertex == NULL)
777 {
778 theVertex = CreateInnerVertex(theGrid);
779 if (theVertex == NULL) return(NULL);
780 V_DIM_COPY(global,CVECT(theVertex));
781 }
782 VFATHER(theVertex) = theElement;
783 SETONSIDE(theVertex,side);
784 SETONNBSIDE(theVertex,SideOfNbElement(theElement,side));
785 V_DIM_COPY(local,LCVECT(theVertex));
786 }
787 /* create node */
788 theNode = CreateNode(theGrid,theVertex,
789 (GEOM_OBJECT *)theElement,SIDE_NODE,1);
790 if (theNode==NULL && vertex_null)
791 {
792 DisposeVertex(theGrid,theVertex);
793 return(NULL);
794 }
795 #ifdef TOPNODE
796 if (TOPNODE(theVertex) == NULL || LEVEL(TOPNODE(theVertex))<LEVEL(theNode))
797 TOPNODE(theVertex) = theNode;
798 #endif
799 theGrid->status |= 1;
800
801 return(theNode);
802 }
803
GetSideNodeX(const ELEMENT * theElement,INT side,INT n,NODE ** MidNodes)804 static NODE *GetSideNodeX (const ELEMENT *theElement, INT side, INT n,
805 NODE **MidNodes)
806 {
807 ELEMENT *theFather;
808 NODE *theNode;
809 VERTEX *theVertex;
810 LINK *theLink0,*theLink1,*theLink2,*theLink3;
811 DOUBLE fac,*local;
812 INT i;
813
814 if (n == 4) {
815 for (theLink0=START(MidNodes[0]); theLink0!=NULL;
816 theLink0=NEXT(theLink0)) {
817 theNode = NBNODE(theLink0);
818 if (NTYPE(theNode) != SIDE_NODE)
819 continue;
820 for (theLink1=START(MidNodes[1]); theLink1!=NULL;
821 theLink1=NEXT(theLink1)) {
822 if (theNode != NBNODE(theLink1))
823 continue;
824 for (theLink2=START(MidNodes[2]); theLink2!=NULL;
825 theLink2=NEXT(theLink2)) {
826 if (theNode != NBNODE(theLink2))
827 continue;
828 for (theLink3=START(MidNodes[3]); theLink3!=NULL;
829 theLink3=NEXT(theLink3)) {
830 if (theNode != NBNODE(theLink3))
831 continue;
832 theVertex = MYVERTEX(theNode);
833 theFather = VFATHER(theVertex);
834 if (theFather == theElement) {
835 #ifndef ModelP
836 /* HEAPFAULT in theFather possible,
837 if in a previous call of DisposeElement
838 some son is not reached by GetAllSons */
839 assert(ONSIDE(theVertex) == side);
840 #endif
841 SETONSIDE(theVertex,side);
842 return(theNode);
843 }
844 else if (theFather == NBELEM(theElement,side)) {
845 SETONNBSIDE(theVertex,side);
846 return(theNode);
847 }
848 else if (theFather == NULL) {
849 /** \todo Strange that this cast has to be here */
850 VFATHER(theVertex) = (ELEMENT*)theElement;
851 SETONSIDE(theVertex,side);
852 SETONNBSIDE(theVertex,
853 SideOfNbElement(theElement,side));
854 fac = 1.0 / n;
855 local = LCVECT(theVertex);
856 V_DIM_CLEAR(local);
857 for (i=0; i<n; i++) {
858 V_DIM_LINCOMB(1.0,local,fac,
859 LOCAL_COORD_OF_ELEM(theElement,
860 CORNER_OF_SIDE(theElement,side,i)),
861 local);
862 }
863 return(theNode);
864 }
865 #ifndef ModelP
866 /* HEAPFAULT in theFather possible,
867 if in a previous call of DisposeElement
868 some son is not reached by GetAllSons */
869 else
870 assert(0);
871 #endif
872 return(theNode);
873 }
874 }
875 }
876 }
877 }
878 else if (n == 3) {
879 for (theLink0=START(MidNodes[0]); theLink0!=NULL;
880 theLink0=NEXT(theLink0)) {
881 theNode = NBNODE(theLink0);
882 if (NTYPE(theNode) != SIDE_NODE)
883 continue;
884 for (theLink1=START(MidNodes[1]); theLink1!=NULL;
885 theLink1=NEXT(theLink1)) {
886 if (theNode != NBNODE(theLink1))
887 continue;
888 for (theLink2=START(MidNodes[2]); theLink2!=NULL;
889 theLink2=NEXT(theLink2)) {
890 if (theNode != NBNODE(theLink2))
891 continue;
892 theVertex = MYVERTEX(theNode);
893 theFather = VFATHER(theVertex);
894 if (theFather == theElement) {
895 if (ONSIDE(theVertex) == side)
896 return(theNode);
897 #ifdef ModelP
898 SETONSIDE(theVertex,side);
899 return(theNode);
900 #endif
901 }
902 else if (theFather == NBELEM(theElement,side))
903 {
904 INT nbside = SideOfNbElement(theElement,side);
905 if (nbside==ONSIDE(theVertex))
906 {
907 SETONNBSIDE(theVertex,side);
908 return(theNode);
909 }
910 #ifdef ModelP
911 /** \todo Strange that this cast has to be here */
912 VFATHER(theVertex) = (ELEMENT*)theElement;
913 SETONSIDE(theVertex,side);
914 SETONNBSIDE(theVertex,nbside);
915 return(theNode);
916 #endif
917 }
918 else if (theFather == NULL) {
919 /** \todo Strange that this cast has to be here */
920 VFATHER(theVertex) = (ELEMENT*)theElement;
921 SETONSIDE(theVertex,side);
922 SETONNBSIDE(theVertex,
923 SideOfNbElement(theElement,side));
924 fac = 1.0 / n;
925 local = LCVECT(theVertex);
926 V_DIM_CLEAR(local);
927 for (i=0; i<n; i++) {
928 V_DIM_LINCOMB(1.0,local,fac,
929 LOCAL_COORD_OF_ELEM(theElement,
930 CORNER_OF_SIDE(theElement,side,i)),
931 local);
932 }
933 return(theNode);
934 }
935 #ifdef ModelP
936 else {
937 return(theNode);
938 }
939 #endif
940 }
941 }
942 }
943 }
944 #ifdef ModelP
945 else if (n == 2) {
946 for (theLink0=START(MidNodes[0]); theLink0!=NULL;
947 theLink0=NEXT(theLink0)) {
948 theNode = NBNODE(theLink0);
949 if (NTYPE(theNode) != SIDE_NODE)
950 continue;
951 for (theLink1=START(MidNodes[1]); theLink1!=NULL;
952 theLink1=NEXT(theLink1)) {
953 if (theNode != NBNODE(theLink1))
954 continue;
955 theVertex = MYVERTEX(theNode);
956 theFather = VFATHER(theVertex);
957 if (theFather == theElement) {
958 if (ONSIDE(theVertex) == side)
959 return(theNode);
960 SETONSIDE(theVertex,side);
961 return(theNode);
962 }
963 else if (theFather == NBELEM(theElement,side)) {
964 SETONNBSIDE(theVertex,side);
965 return(theNode);
966 }
967 return(theNode);
968 }
969 }
970 }
971 #endif
972
973 return(NULL);
974 }
975
GetSideNode(const ELEMENT * theElement,INT side)976 NODE * NS_DIM_PREFIX GetSideNode (const ELEMENT *theElement, INT side)
977 {
978 NODE *theNode;
979 NODE *MidNodes[MAX_EDGES_OF_SIDE];
980 INT i,n;
981 #ifdef ModelP
982 INT k;
983 #endif
984
985 n = 0;
986 for (i=0; i<EDGES_OF_SIDE(theElement,side); i++) {
987 theNode = GetMidNode(theElement,EDGE_OF_SIDE(theElement,side,i));
988 if (theNode != NULL)
989 MidNodes[n++] = theNode;
990 #ifndef ModelP
991 else return(NULL);
992 #endif
993 }
994 PRINTDEBUG(gm,2,("GetSideNode(): elem=" EID_FMTX
995 " side=%d nb. of midnodes=%d\n",
996 EID_PRTX(theElement),side,n));
997 #ifdef ModelP
998 if (TAG(theElement)==PYRAMID && side!=0) return(NULL);
999 #endif
1000 theNode = GetSideNodeX(theElement,side,n,MidNodes);
1001 #ifdef ModelP
1002 if (theNode != NULL)
1003 return(theNode);
1004 if (n < 3)
1005 return(NULL);
1006 for (i=0; i<n; i++)
1007 {
1008 NODE *MidNodes1[MAX_EDGES_OF_SIDE-1];
1009 INT j,m;
1010
1011 m = 0;
1012 for (j=0; j<n; j++) {
1013 if (i == j) continue;
1014 MidNodes1[m++] = MidNodes[j];
1015 }
1016 theNode = GetSideNodeX(theElement,side,n-1,MidNodes1);
1017 if (theNode != NULL)
1018 return(theNode);
1019 }
1020 if (n < 4)
1021 return(NULL);
1022 for (i=1; i<n; i++)
1023 for (k=0; k<i; k++)
1024 {
1025 NODE *MidNodes1[MAX_EDGES_OF_SIDE-2];
1026 INT j,m;
1027
1028 m = 0;
1029 for (j=0; j<n; j++) {
1030 if (i == j) continue;
1031 if (k == j) continue;
1032 MidNodes1[m++] = MidNodes[j];
1033 }
1034 theNode = GetSideNodeX(theElement,side,n-2,MidNodes1);
1035 if (theNode != NULL)
1036 return(theNode);
1037 }
1038 #endif
1039
1040 return(theNode);
1041 }
1042
CountSideNodes(ELEMENT * e)1043 static int CountSideNodes (ELEMENT *e)
1044 {
1045 int i,side;
1046 NODE *n;
1047
1048 side = 0;
1049 for (i=0; i<CORNERS_OF_ELEM(e); i++)
1050 {
1051 n = CORNER(e,i);
1052 if (SIDETYPE(n)) side++;
1053 }
1054 return(side);
1055 }
1056
GetSideIDFromScratchSpecialRule17Pyr(ELEMENT * theElement,NODE * theNode)1057 int GetSideIDFromScratchSpecialRule17Pyr (ELEMENT *theElement, NODE *theNode)
1058 {
1059 int i,k,l,nodes;
1060 #ifdef Debug
1061 INT cnodes,snodes;
1062 #endif
1063 ELEMENT *f = EFATHER(theElement);
1064 NODE *fnode,*enode;
1065 int side = SIDES_OF_ELEM(f);
1066
1067 #ifdef Debug
1068 assert(TAG(theElement)==PYRAMID);
1069 snodes = cnodes = 0;
1070 for (l=0; l<CORNERS_OF_ELEM(theElement); l++)
1071 {
1072 enode = CORNER(theElement,l);
1073 if (CORNERTYPE(enode)) cnodes++;
1074 if (SIDETYPE(enode)) snodes++;
1075 }
1076 assert(snodes == 1);
1077 assert(cnodes == 4);
1078 #endif
1079
1080 for (i=0; i<SIDES_OF_ELEM(f); i++)
1081 {
1082 nodes = 0;
1083 for (k=0; k<CORNERS_OF_SIDE(f,i); k++)
1084 {
1085 fnode = CORNER(f,CORNER_OF_SIDE(f,i,k));
1086 for (l=0; l<CORNERS_OF_ELEM(theElement); l++)
1087 {
1088 enode = CORNER(theElement,l);
1089 if (enode == SONNODE(fnode)) nodes++;
1090 }
1091 }
1092 assert(nodes==0 || nodes==2 || nodes==4);
1093 #ifdef Debug
1094 if (nodes == 0) side = i;
1095 #else
1096 if (nodes == 0) return(i);
1097 #endif
1098 }
1099
1100 assert(side<SIDES_OF_ELEM(f));
1101 return(side);
1102 }
1103
1104
GetSideIDFromScratchSpecialRule22Tet(ELEMENT * theElement,NODE * theNode)1105 int GetSideIDFromScratchSpecialRule22Tet (ELEMENT *theElement, NODE *theNode)
1106 {
1107 int i,k,l,nodes,midnodes;
1108 #ifdef Debug
1109 INT cnodes, snodes, mnodes;
1110 #endif
1111 ELEMENT *f = EFATHER(theElement);
1112 NODE *fnode,*enode;
1113 EDGE *edge;
1114 int side = SIDES_OF_ELEM(f);
1115
1116 #ifdef Debug
1117 assert(TAG(theElement)==TETRAHEDRON);
1118 snodes = cnodes = mnodes = 0;
1119 for (l=0; l<CORNERS_OF_ELEM(theElement); l++)
1120 {
1121 enode = CORNER(theElement,l);
1122 if (CORNERTYPE(enode)) cnodes++;
1123 if (MIDTYPE(enode)) mnodes++;
1124 if (SIDETYPE(enode)) snodes++;
1125 }
1126 assert(cnodes == 2);
1127 assert(mnodes == 1);
1128 assert(snodes == 1);
1129 #endif
1130
1131 for (i=0; i<SIDES_OF_ELEM(f); i++)
1132 {
1133 nodes = 0;
1134 midnodes = 0;
1135 for (k=0; k<CORNERS_OF_SIDE(f,i); k++)
1136 {
1137 fnode = CORNER(f,CORNER_OF_SIDE(f,i,k));
1138
1139 edge = GetEdge(CORNER_OF_SIDE_PTR(f,i,k),
1140 CORNER_OF_SIDE_PTR(f,i,(k+1)%CORNERS_OF_SIDE(f,i)));
1141 assert(edge != NULL);
1142
1143 for (l=0; l<CORNERS_OF_ELEM(theElement); l++)
1144 {
1145 enode = CORNER(theElement,l);
1146 if (enode == SONNODE(fnode)) nodes++;
1147 if (enode == MIDNODE(edge)) midnodes++;
1148 }
1149 }
1150 assert(nodes==0 || nodes==1 || nodes==2 || nodes==4);
1151 #ifdef Debug
1152 if (nodes==0 && midnodes==1) side = i;
1153 #else
1154 if (nodes==0 && midnodes==1) return(i);
1155 #endif
1156 }
1157
1158 assert(side<SIDES_OF_ELEM(f));
1159 return(side);
1160 }
1161
1162
GetSideIDFromScratchSpecialRule(ELEMENT * theElement,NODE * theNode)1163 INT GetSideIDFromScratchSpecialRule (ELEMENT *theElement, NODE *theNode)
1164 {
1165 int j,l;
1166 #ifndef NDEBUG
1167 ELEMENT *f = EFATHER(theElement);
1168 #endif
1169
1170 assert(TAG(f)==HEXAHEDRON);
1171 assert(ECLASS(theElement)==GREEN_CLASS);
1172 assert(NSONS(f)==9 || NSONS(f)==11 || EHGHOST(theElement));
1173
1174 if (TAG(theElement)==PYRAMID)
1175 {
1176 return(GetSideIDFromScratchSpecialRule17Pyr(theElement,theNode));
1177 }
1178
1179 assert(TAG(theElement)==TETRAHEDRON);
1180 /* centroid tetrahedron of special rule 22 */
1181 if (CountSideNodes(theElement) == 2)
1182 {
1183 /* if side not found search over neighbor */
1184 for (j=0; j<SIDES_OF_ELEM(theElement); j++)
1185 {
1186 ELEMENT *nb = NBELEM(theElement,j);
1187
1188 if (nb == NULL) continue;
1189
1190 for (l=0; l<CORNERS_OF_ELEM(nb); l++)
1191 if (theNode == CORNER(nb,l))
1192 return(GetSideIDFromScratch(nb,theNode));
1193 }
1194 }
1195
1196 assert(CountSideNodes(theElement)==1);
1197
1198 return(GetSideIDFromScratchSpecialRule22Tet(theElement,theNode));
1199 }
1200
GetSideIDFromScratch(ELEMENT * theElement,NODE * theNode)1201 INT NS_DIM_PREFIX GetSideIDFromScratch (ELEMENT *theElement, NODE *theNode)
1202 {
1203 ELEMENT *theFather;
1204 NODE *nd[MAX_EDGES_OF_ELEM];
1205 EDGE *edge;
1206 INT i,j,k,l,cnt;
1207
1208 ASSERT(NTYPE(theNode) == SIDE_NODE);
1209
1210 theFather = EFATHER(theElement);
1211
1212 /* determine midnodes of father */
1213 for (i=0; i<EDGES_OF_ELEM(theFather); i++)
1214 {
1215 edge = GetEdge(CORNER_OF_EDGE_PTR(theFather,i,0),
1216 CORNER_OF_EDGE_PTR(theFather,i,1));
1217 nd[i] = MIDNODE(edge);
1218 }
1219
1220 for (j=0; j<SIDES_OF_ELEM(theElement); j++)
1221 {
1222 if (3 == CORNERS_OF_SIDE(theElement,j)) continue;
1223
1224 for (l=0; l<CORNERS_OF_SIDE(theElement,j); l++)
1225 if (theNode == CORNER(theElement,CORNER_OF_SIDE(theElement,j,l)))
1226 break;
1227 if (l == CORNERS_OF_SIDE(theElement,j)) continue;
1228
1229 for (i=0; i<SIDES_OF_ELEM(theFather); i++)
1230 {
1231 #ifdef DUNE_UGGRID_DUNE_UGGRID_TET_RULESET
1232 if (3 == CORNERS_OF_SIDE(theFather,i)) continue;
1233 #endif
1234
1235 cnt = 0;
1236 for (k=0; k<EDGES_OF_SIDE(theFather,i); k++)
1237 for (l=0; l<CORNERS_OF_SIDE(theElement,j); l++)
1238 {
1239 if (nd[EDGE_OF_SIDE(theFather,i,k)] ==
1240 CORNER(theElement,CORNER_OF_SIDE(theElement,j,l)))
1241 cnt++;
1242 if (cnt == 2)
1243 return(i);
1244 }
1245 }
1246 }
1247
1248
1249 /* if side not found search over neighbor */
1250 for (j=0; j<SIDES_OF_ELEM(theElement); j++)
1251 {
1252 ELEMENT *nb = NBELEM(theElement,j);
1253
1254 if (3 == CORNERS_OF_SIDE(theElement,j))
1255 continue;
1256
1257 if (nb == NULL) continue;
1258
1259 for (l=0; l<CORNERS_OF_ELEM(nb); l++)
1260 if (theNode == CORNER(nb,l))
1261 return(GetSideIDFromScratch(nb,theNode));
1262 }
1263
1264
1265 for (j=0; j<SIDES_OF_ELEM(theElement); j++)
1266 {
1267 if (4 != CORNERS_OF_SIDE(theElement,j)) continue;
1268 for (l=0; l<4; l++)
1269 if (theNode == CORNER(theElement,CORNER_OF_SIDE(theElement,j,l)))
1270 break;
1271 if (l < 4)
1272 {
1273 INT l1 = (l+1) % 4;
1274
1275 for (i=0; i<SIDES_OF_ELEM(theFather); i++) {
1276 if (3 == CORNERS_OF_SIDE(theFather,i)) continue;
1277 for (k=0; k<EDGES_OF_SIDE(theFather,i); k++) {
1278 if (nd[EDGE_OF_SIDE(theFather,i,k)] ==
1279 CORNER(theElement,CORNER_OF_SIDE(theElement,j,l1)))
1280 return(i);
1281 if (nd[EDGE_OF_SIDE(theFather,i,k)] ==
1282 CORNER(theElement,CORNER_OF_SIDE(theElement,j,l1)))
1283 return(i);
1284 }
1285 }
1286 }
1287 }
1288
1289 return(GetSideIDFromScratchSpecialRule(theElement,theNode));
1290
1291 return(SIDES_OF_ELEM(theFather));
1292 }
1293
GetSideIDFromScratchOld(ELEMENT * theElement,NODE * theNode)1294 INT GetSideIDFromScratchOld (ELEMENT *theElement, NODE *theNode)
1295 {
1296 ELEMENT *theFather;
1297 NODE *nd[MAX_EDGES_OF_ELEM];
1298 EDGE *edge;
1299 INT i,j,k,l,cnt;
1300
1301 ASSERT(NTYPE(theNode) == SIDE_NODE);
1302
1303 theFather = EFATHER(theElement);
1304
1305 /* determine midnodes of father */
1306 for (i=0; i<EDGES_OF_ELEM(theFather); i++)
1307 {
1308 edge = GetEdge(CORNER_OF_EDGE_PTR(theFather,i,0),
1309 CORNER_OF_EDGE_PTR(theFather,i,1));
1310 nd[i] = MIDNODE(edge);
1311 }
1312
1313 for (j=0; j<SIDES_OF_ELEM(theElement); j++)
1314 {
1315 if (3 == CORNERS_OF_SIDE(theElement,j)) continue;
1316
1317 for (l=0; l<CORNERS_OF_SIDE(theElement,j); l++)
1318 if (theNode == CORNER(theElement,CORNER_OF_SIDE(theElement,j,l)))
1319 break;
1320 if (l == CORNERS_OF_SIDE(theElement,j)) continue;
1321
1322 for (i=0; i<SIDES_OF_ELEM(theFather); i++)
1323 {
1324 if (3 == CORNERS_OF_SIDE(theFather,i)) continue;
1325
1326 cnt = 0;
1327 for (k=0; k<EDGES_OF_SIDE(theFather,i); k++)
1328 for (l=0; l<CORNERS_OF_SIDE(theElement,j); l++)
1329 {
1330 if (nd[EDGE_OF_SIDE(theFather,i,k)] ==
1331 CORNER(theElement,CORNER_OF_SIDE(theElement,j,l)))
1332 cnt++;
1333 if (cnt == 2)
1334 return(i);
1335 }
1336 }
1337 }
1338
1339
1340 /* if side not found search over neighbor */
1341 for (j=0; j<SIDES_OF_ELEM(theElement); j++)
1342 {
1343 ELEMENT *nb = NBELEM(theElement,j);
1344
1345 if (3 == CORNERS_OF_SIDE(theElement,j))
1346 {
1347
1348 /* treatment of special green rule 17 and 22 */
1349 if (((TAG(theElement)==PYRAMID && NSONS(theFather)==9) ||
1350 (TAG(theElement)==TETRAHEDRON && NSONS(theFather)==11))
1351 && 2==CountSideNodes(theElement) &&
1352 TAG(theFather)==HEXAHEDRON &&
1353 ECLASS(theElement)==GREEN_CLASS)
1354 /* not continue */;
1355 else
1356 continue;
1357 }
1358
1359 if (nb == NULL) continue;
1360
1361 for (l=0; l<CORNERS_OF_ELEM(nb); l++)
1362 if (theNode == CORNER(nb,l))
1363 return(GetSideIDFromScratch(nb,theNode));
1364 }
1365
1366
1367 for (j=0; j<SIDES_OF_ELEM(theElement); j++)
1368 {
1369 if (4 != CORNERS_OF_SIDE(theElement,j)) continue;
1370 for (l=0; l<4; l++)
1371 if (theNode == CORNER(theElement,CORNER_OF_SIDE(theElement,j,l)))
1372 break;
1373 if (l < 4)
1374 {
1375 INT l1 = (l+1) % 4;
1376
1377 for (i=0; i<SIDES_OF_ELEM(theFather); i++) {
1378 if (3 == CORNERS_OF_SIDE(theFather,i)) continue;
1379 for (k=0; k<EDGES_OF_SIDE(theFather,i); k++) {
1380 if (nd[EDGE_OF_SIDE(theFather,i,k)] ==
1381 CORNER(theElement,CORNER_OF_SIDE(theElement,j,l1)))
1382 return(i);
1383 if (nd[EDGE_OF_SIDE(theFather,i,k)] ==
1384 CORNER(theElement,CORNER_OF_SIDE(theElement,j,l1)))
1385 return(i);
1386 }
1387 }
1388 }
1389 }
1390
1391 /* treatment of special green rule 17 and 22 */
1392 for (j=0; j<SIDES_OF_ELEM(theElement); j++)
1393 {
1394 for (l=0; l<CORNERS_OF_SIDE(theElement,j); l++)
1395 if (theNode == CORNER(theElement,CORNER_OF_SIDE(theElement,j,l)))
1396 break;
1397 if (l == CORNERS_OF_SIDE(theElement,j)) continue;
1398
1399 for (i=0; i<SIDES_OF_ELEM(theFather); i++)
1400 {
1401 if (3 == CORNERS_OF_SIDE(theFather,i)) continue;
1402
1403 cnt = 0;
1404 for (k=0; k<EDGES_OF_SIDE(theFather,i); k++)
1405 for (l=0; l<CORNERS_OF_SIDE(theElement,j); l++)
1406 {
1407 if (nd[EDGE_OF_SIDE(theFather,i,k)] ==
1408 CORNER(theElement,CORNER_OF_SIDE(theElement,j,l)))
1409 cnt++;
1410 if (cnt==1 && ECLASS(theElement)==GREEN_CLASS &&
1411 TAG(theElement)==TETRAHEDRON &&
1412 TAG(theFather)==HEXAHEDRON &&
1413 (NSONS(theFather)==9 || NSONS(theFather)==11))
1414 {
1415 return(i);
1416 }
1417
1418 }
1419 }
1420 }
1421
1422 UserWriteF("GetSideIDFromScratch(): e=" EID_FMTX " f=" EID_FMTX "\n",
1423 EID_PRTX(theElement),EID_PRTX(theFather));
1424 return(0);
1425 return(SIDES_OF_ELEM(theFather));
1426 }
1427
1428 #endif /* __THREEDIM__ */
1429
1430
1431 /****************************************************************************/
1432 /** \brief Get the center node of an element of next finer level
1433 *
1434 * This function gets the center node of an element of next finer level
1435
1436 @return <ul>
1437 <li> pointer to center node </li>
1438 <li> NULL no node found </li>
1439 </ul> */
1440 /****************************************************************************/
1441
GetCenterNode(const ELEMENT * theElement)1442 NODE * NS_DIM_PREFIX GetCenterNode (const ELEMENT *theElement)
1443 {
1444 INT i,j;
1445 NODE *theNode;
1446 ELEMENT *SonList[MAX_SONS],*theSon;
1447
1448 #ifdef __CENTERNODE__
1449 return(CENTERNODE(theElement));
1450 #endif
1451
1452 theNode = NULL;
1453 if (GetAllSons(theElement,SonList) != GM_OK) assert(0);
1454
1455 for (i=0; SonList[i]!=NULL; i++)
1456 {
1457 theSon = SonList[i];
1458 for (j=0; j<CORNERS_OF_ELEM(theSon); j++)
1459 {
1460 theNode = CORNER(theSon,j);
1461 if (NTYPE(theNode) == CENTER_NODE)
1462 {
1463 if (EMASTER(theElement))
1464 assert(VFATHER(MYVERTEX(theNode)) == theElement);
1465 return (theNode);
1466 }
1467 }
1468 }
1469 return (NULL);
1470 }
1471
1472 /****************************************************************************/
1473 /** \brief Allocate a new node on a side of an element
1474 *
1475 * Includes vertex
1476 * best fit boundary coordinates and local coordinates.
1477 *
1478 * @return <ul>
1479 * <li> pointer to new node </li>
1480 * <li> NULL: could not allocate </li>
1481 </ul> */
1482 /****************************************************************************/
1483 /* #define MOVE_MIDNODE */
CreateCenterNode(GRID * theGrid,ELEMENT * theElement,VERTEX * theVertex)1484 NODE * NS_DIM_PREFIX CreateCenterNode (GRID *theGrid, ELEMENT *theElement, VERTEX *theVertex)
1485 {
1486 DOUBLE *global,*local;
1487 DOUBLE_VECTOR diff;
1488 INT n,j,moved,vertex_null;
1489 VERTEX *VertexOnEdge[MAX_EDGES_OF_ELEM];
1490 NODE *theNode;
1491 EDGE *theEdge;
1492 DOUBLE fac;
1493 DOUBLE *x[MAX_CORNERS_OF_ELEM];
1494 #ifdef MOVE_MIDNODE
1495 #ifndef ModelP
1496 DOUBLE len_opp,len_bnd;
1497 #endif
1498 #endif
1499
1500 /* check if moved side nodes exist */
1501 CORNER_COORDINATES(theElement,n,x);
1502 moved = 0;
1503 vertex_null = (theVertex==NULL);
1504 if (theVertex==NULL && OBJT(theElement) == BEOBJ) {
1505 for (j=0; j<EDGES_OF_ELEM(theElement); j++) {
1506 theEdge=GetEdge(CORNER(theElement,CORNER_OF_EDGE(theElement,j,0)),
1507 CORNER(theElement,CORNER_OF_EDGE(theElement,j,1)));
1508 ASSERT(theEdge != NULL);
1509 theNode = MIDNODE(theEdge);
1510 if (theNode == NULL)
1511 VertexOnEdge[j] = NULL;
1512 else {
1513 VertexOnEdge[j] = MYVERTEX(theNode);
1514 moved += MOVED(VertexOnEdge[j]);
1515 }
1516 }
1517 #ifdef MOVE_MIDNODE
1518 #ifndef ModelP
1519 if (moved == 1) {
1520 for (j=0; j<EDGES_OF_ELEM(theElement); j++)
1521 if (VertexOnEdge[j] != NULL)
1522 if (MOVED(VertexOnEdge[j])) break;
1523 theVertex = VertexOnEdge[OPPOSITE_EDGE(theElement,j)];
1524 if (theVertex != NULL) {
1525 V_DIM_LINCOMB(0.5,CVECT(MYVERTEX(CORNER(theElement,CORNER_OF_EDGE(theElement,j,0)))),
1526 0.5,CVECT(MYVERTEX(CORNER(theElement,CORNER_OF_EDGE(theElement,j,1)))),
1527 diff);
1528 V_DIM_LINCOMB(1.0,CVECT(VertexOnEdge[j]),-1.0,diff,diff);
1529 /* scale diff according to length of edges */
1530 V_DIM_EUKLIDNORM_OF_DIFF(
1531 CVECT(MYVERTEX(CORNER(theElement,CORNER_OF_EDGE(theElement,j,0)))),
1532 CVECT(MYVERTEX(CORNER(theElement,CORNER_OF_EDGE(theElement,j,1)))),len_bnd);
1533 V_DIM_EUKLIDNORM_OF_DIFF(
1534 CVECT(MYVERTEX(CORNER(theElement,
1535 CORNER_OF_EDGE(theElement,OPPOSITE_EDGE(theElement,j),0)))),
1536 CVECT(MYVERTEX(CORNER(theElement,
1537 CORNER_OF_EDGE(theElement,OPPOSITE_EDGE(theElement,j),1)))),
1538 len_opp);
1539 V_DIM_SCALE(len_opp/len_bnd,diff);
1540 global = CVECT(theVertex);
1541 PRINTDEBUG(gm,1,("CreateCenterNode: global_orig = %f %f %f\n",global[0],global[1],global[2]));
1542 PRINTDEBUG(gm,1,("CreateCenterNode: diff = %f %f %f\n",diff[0],diff[1],diff[2]));
1543 V_DIM_LINCOMB(1.0,global,0.5,diff,global);
1544 PRINTDEBUG(gm,1,("CreateCenterNode: global_mod = %f %f %f\n",global[0],global[1],global[2]));
1545 SETMOVED(VertexOnEdge[OPPOSITE_EDGE(theElement,j)],1);
1546 UG_GlobalToLocal(n,(const DOUBLE **)x,global,LCVECT(theVertex));
1547 PRINTDEBUG(gm,1,("CreateCenterNode: local = %f %f %f\n",LCVECT(theVertex)[0],LCVECT(theVertex)[1],LCVECT(theVertex)[2]));
1548 SETONEDGE(theVertex,OPPOSITE_EDGE(theElement,j));
1549 VFATHER(theVertex) = theElement;
1550 }
1551 }
1552 #endif
1553 #endif
1554 }
1555
1556 if (vertex_null)
1557 {
1558 theVertex = CreateInnerVertex(theGrid);
1559 if (theVertex==NULL)
1560 return(NULL);
1561 VFATHER(theVertex) = theElement;
1562 }
1563
1564 theNode = CreateNode(theGrid,theVertex,(GEOM_OBJECT *)theElement,CENTER_NODE,1);
1565 if (theNode==NULL && vertex_null)
1566 {
1567 DisposeVertex(theGrid,theVertex);
1568 return(NULL);
1569 }
1570
1571 #ifdef TOPNODE
1572 if (TOPNODE(theVertex) == NULL || LEVEL(TOPNODE(theVertex))<LEVEL(theNode))
1573 TOPNODE(theVertex) = theNode;
1574 #endif
1575 theGrid->status |= 1;
1576
1577 if (!vertex_null) return(theNode);
1578
1579 global = CVECT(theVertex);
1580 local = LCVECT(theVertex);
1581 V_DIM_CLEAR(local);
1582 fac = 1.0 / n;
1583 for (j=0; j<n; j++)
1584 V_DIM_LINCOMB(1.0,local,
1585 fac,LOCAL_COORD_OF_ELEM(theElement,j),local);
1586 LOCAL_TO_GLOBAL(n,x,local,global);
1587 if (moved) {
1588 V_DIM_CLEAR(diff);
1589 for (j=0; j<EDGES_OF_ELEM(theElement); j++)
1590 if (VertexOnEdge[j] != NULL) {
1591 V_DIM_COPY(CVECT(VertexOnEdge[j]),diff);
1592 V_DIM_LINCOMB(1.0,diff,-0.5,CVECT(MYVERTEX(CORNER(theElement,CORNER_OF_EDGE(theElement,j,0)))),diff);
1593 V_DIM_LINCOMB(1.0,diff,-0.5,CVECT(MYVERTEX(CORNER(theElement,CORNER_OF_EDGE(theElement,j,1)))),diff);
1594 V_DIM_LINCOMB(0.5,diff,1.0,global,global);
1595 }
1596 UG_GlobalToLocal(n,(const DOUBLE **)x,global,local);
1597 LOCAL_TO_GLOBAL(n,x,local,diff);
1598 SETMOVED(theVertex,1);
1599 }
1600 return(theNode);
1601 }
1602
1603
1604 /****************************************************************************/
1605 /** \brief Get all nodes related to this element on next level
1606
1607 * @param theElement - element for context
1608 * @param theElementContext - node context of this element
1609
1610 This function returns the nodes related to the element on the next
1611 finer level. The ordering is according to the reference numbering.
1612
1613 @return <ul>
1614 <li> GM_OK if ok </li>
1615 <li> != GM_OK if not ok </li>
1616 </ul> */
1617 /****************************************************************************/
1618
GetNodeContext(const ELEMENT * theElement,NODE ** theElementContext)1619 INT NS_DIM_PREFIX GetNodeContext (const ELEMENT *theElement, NODE **theElementContext)
1620 {
1621 NODE *theNode, **MidNodes, **CenterNode;
1622 EDGE *theEdge;
1623 INT i,Corner0, Corner1;
1624 #ifdef __THREEDIM__
1625 NODE **SideNodes;
1626 #endif
1627
1628 /* reset context */
1629 for(i=0; i<MAX_CORNERS_OF_ELEM+MAX_NEW_CORNERS_DIM; i++)
1630 theElementContext[i] = NULL;
1631
1632 /* is element to refine */
1633 if (!IS_REFINED(theElement)) return(GM_OK);
1634
1635 /* get corner nodes */
1636 for (i=0; i<CORNERS_OF_ELEM(theElement); i++)
1637 {
1638 theNode = CORNER(theElement,i);
1639 theElementContext[i] = SONNODE(theNode);
1640 }
1641
1642 /* check for midpoint nodes */
1643 MidNodes = theElementContext+CORNERS_OF_ELEM(theElement);
1644 for (i=0; i<EDGES_OF_ELEM(theElement); i++)
1645 {
1646 Corner0 = CORNER_OF_EDGE(theElement,i,0);
1647 Corner1 = CORNER_OF_EDGE(theElement,i,1);
1648
1649 theEdge = GetEdge(CORNER(theElement,Corner0),
1650 CORNER(theElement,Corner1));
1651 ASSERT(theEdge != NULL);
1652
1653 MidNodes[i] = MIDNODE(theEdge);
1654 }
1655
1656 #ifdef __THREEDIM__
1657 SideNodes = theElementContext+CORNERS_OF_ELEM(theElement)+
1658 EDGES_OF_ELEM(theElement);
1659 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
1660 {
1661 #ifdef DUNE_UGGRID_DUNE_UGGRID_TET_RULESET
1662 /* no side nodes for triangular sides yet */
1663 if (CORNERS_OF_SIDE(theElement,i) == 3) continue;
1664 #endif
1665 /* check for side node */
1666 SideNodes[i] = GetSideNode(theElement,i);
1667 }
1668 #endif
1669
1670 /* check for center node */
1671 CenterNode = MidNodes+CENTER_NODE_INDEX(theElement);
1672 CenterNode[0] = GetCenterNode(theElement);
1673
1674 return(GM_OK);
1675 }
1676
1677 /****************************************************************************/
1678 /** \brief Return matching side of the neighboring element
1679
1680 * @param theNeighbor - element to test for matching side
1681 * @param nbside - the matching side
1682 * @param theElement - element with side to match
1683 * @param side - side of element to match
1684
1685 This function computes the matching side of neighboring element.
1686
1687 */
1688 /****************************************************************************/
1689
GetNbSideByNodes(ELEMENT * theNeighbor,INT * nbside,ELEMENT * theElement,INT side)1690 void NS_DIM_PREFIX GetNbSideByNodes (ELEMENT *theNeighbor, INT *nbside, ELEMENT *theElement, INT side)
1691 {
1692 INT i,k,l,ec,nc;
1693
1694 ec = CORNERS_OF_SIDE(theElement,side);
1695
1696 for (i=0; i<SIDES_OF_ELEM(theNeighbor); i++)
1697 {
1698 nc = CORNERS_OF_SIDE(theNeighbor,i);
1699 if (ec != nc) continue;
1700
1701 for (k=0; k<nc; k++)
1702 if (CORNER_OF_SIDE_PTR(theElement,side,0) ==
1703 CORNER_OF_SIDE_PTR(theNeighbor,i,k))
1704 break;
1705 if (k == nc) continue;
1706
1707 for (l=1; l<ec; l++)
1708 {
1709 if (CORNER_OF_SIDE_PTR(theElement,side,l) !=
1710 CORNER_OF_SIDE_PTR(theNeighbor,i,(nc+k-l)%nc))
1711 break;
1712 }
1713 if (l == ec)
1714 {
1715 *nbside = i;
1716 return;
1717 }
1718 }
1719
1720 /* no side of the neighbor matches */
1721 assert(0);
1722 }
1723
1724 /****************************************************************************/
1725 /** \brief Return pointer to son edge if it exists
1726 *
1727 * @param theEdge edge for which son is searched
1728 *
1729 * This function returns the pointer to the son edge if it exists.
1730 *
1731 * @return <ul>
1732 <li> pointer to specified object </li>
1733 <li> NULL if not found </li>
1734 </ul> */
1735 /****************************************************************************/
1736
GetSonEdge(const EDGE * theEdge)1737 EDGE * NS_DIM_PREFIX GetSonEdge (const EDGE *theEdge)
1738 {
1739 EDGE *SonEdge=NULL;
1740 NODE *Node0,*Node1,*SonNode0,*SonNode1;
1741
1742 Node0 = NBNODE(LINK0(theEdge));
1743 Node1 = NBNODE(LINK1(theEdge));
1744
1745 SonNode0 = SONNODE(Node0);
1746 SonNode1 = SONNODE(Node1);
1747
1748 if (SonNode0!=NULL && SonNode1!=NULL)
1749 SonEdge = GetEdge(SonNode0,SonNode1);
1750
1751 return(SonEdge);
1752 }
1753
1754
1755 /****************************************************************************/
1756 /** \brief Return pointer to son edges if it exists
1757
1758 * @param theEdge - edge for which son is searched
1759 * @param SonEdges - array of pointers will be filled with son edges
1760
1761 This function returns the pointer to the son edges if existing.
1762
1763 @return <ul>
1764 <li> number of found edges (0,1,2) </li>
1765 </ul> */
1766 /****************************************************************************/
1767
GetSonEdges(const EDGE * theEdge,EDGE * SonEdges[MAX_SON_EDGES])1768 INT NS_DIM_PREFIX GetSonEdges (const EDGE *theEdge, EDGE *SonEdges[MAX_SON_EDGES])
1769 {
1770 INT nedges;
1771 NODE *Node0,*Node1,*SonNode0,*SonNode1,*MidNode;
1772
1773 nedges = 0;
1774 SonEdges[0] = NULL;
1775 SonEdges[1] = NULL;
1776
1777 Node0 = NBNODE(LINK0(theEdge));
1778 Node1 = NBNODE(LINK1(theEdge));
1779
1780 if (GID(Node0)<GID(Node1))
1781 {
1782 SonNode0 = SONNODE(Node0);
1783 SonNode1 = SONNODE(Node1);
1784 }
1785 else
1786 {
1787 SonNode0 = SONNODE(Node1);
1788 SonNode1 = SONNODE(Node0);
1789 }
1790 MidNode = MIDNODE(theEdge);
1791
1792 /* parallel note: */
1793 /* since existance of MidNode decides whether for one SonEdge or */
1794 /* two half SonEdges is searched, the data structure must be */
1795 /* consistent in a way that if the MidNode exists also the */
1796 /* MIDNODE pointer is set to MidNode. (s.l. 980227) */
1797 if (MidNode == NULL)
1798 {
1799 if (SonNode0!=NULL && SonNode1!=NULL)
1800 SonEdges[0] = GetEdge(SonNode0,SonNode1);
1801 }
1802 else
1803 {
1804 if (SonNode0!=NULL)
1805 SonEdges[0] = GetEdge(SonNode0,MidNode);
1806
1807 if (SonNode1!=NULL)
1808 SonEdges[1] = GetEdge(MidNode,SonNode1);
1809 }
1810
1811 if (SonEdges[0] != NULL) nedges++;
1812 if (SonEdges[1] != NULL) nedges++;
1813
1814 return(nedges);
1815 }
1816
1817 /****************************************************************************/
1818 /** \brief Return pointer to father edge if it exists
1819
1820 * @param theEdge - edge for which father is searched
1821
1822 This function returns the pointer to the father edge if it exists.
1823
1824 @return <ul>
1825 <li> pointer to specified object </li>
1826 <li> NULL if not found </li>
1827 </ul> */
1828 /****************************************************************************/
1829
GetFatherEdge(const EDGE * theEdge)1830 EDGE * NS_DIM_PREFIX GetFatherEdge (const EDGE *theEdge)
1831 {
1832 NODE *theNode0 = NBNODE(LINK0(theEdge));
1833 NODE *theNode1 = NBNODE(LINK1(theEdge));
1834 EDGE *FatherEdge = NULL;
1835
1836 /* one node is center node -> no father edge */
1837 if (CENTERTYPE(theNode0) || CENTERTYPE(theNode1)) return(NULL);
1838
1839 #ifdef __THREEDIM__
1840 /* one node is side node -> no father edge */
1841 if (SIDETYPE(theNode0) || SIDETYPE(theNode1)) return(NULL);
1842 #endif
1843
1844 /* both nodes are mid nodes -> no father edge */
1845 if (MIDTYPE(theNode0) && MIDTYPE(theNode1)) return(NULL);
1846
1847 /* one node is mid node -> no father edge */
1848 if (MIDTYPE(theNode0) || MIDTYPE(theNode1))
1849 {
1850 NODE *FatherNode0,*FatherNode1, *theNode;
1851
1852 if (MIDTYPE(theNode1))
1853 {
1854 theNode = theNode0; theNode0 = theNode1; theNode1 = theNode;
1855 }
1856 FatherEdge = (EDGE *) NFATHER(theNode0);
1857 if (FatherEdge == NULL) return(NULL);
1858
1859 FatherNode0 = NBNODE(LINK0(FatherEdge));
1860 FatherNode1 = NBNODE(LINK1(FatherEdge));
1861 if (SONNODE(FatherNode0)==theNode1 || SONNODE(FatherNode1)==theNode1)
1862 return(FatherEdge);
1863 else
1864 return(NULL);
1865 }
1866
1867 /* both nodes are corner nodes -> try to get the edge */
1868 if (CORNERTYPE(theNode0) && CORNERTYPE(theNode1))
1869 {
1870 if (NFATHER(theNode0)!=NULL && NFATHER(theNode1)!=NULL)
1871 return(GetEdge(NFATHER(theNode0),NFATHER(theNode1)));
1872 else
1873 return(NULL);
1874 }
1875
1876 /* No father available */
1877 return NULL;
1878 }
1879
1880 #ifdef __THREEDIM__
1881
1882 /****************************************************************************/
1883 /** \brief Return pointer to father edge if it exists
1884
1885 * @param SideNodes - nodes of the side
1886 * @param ncorners - number of sidenodes
1887 * @param Nodes - corners of edge for which father is searched
1888 * @param theEdge - edge for which father is searched
1889
1890 This function returns the pointer to the father edge if it exists.
1891
1892 @return <ul>
1893 <li> pointer to specified object </li>
1894 <li> NULL if not found </li>
1895 </ul> */
1896 /****************************************************************************/
1897
FatherEdge(NODE ** SideNodes,INT ncorners,NODE ** Nodes,EDGE * theEdge)1898 EDGE * NS_DIM_PREFIX FatherEdge (NODE **SideNodes, INT ncorners, NODE **Nodes, EDGE *theEdge)
1899 {
1900 INT pos0,pos1;
1901 EDGE *fatherEdge = NULL;
1902
1903 ASSERT(Nodes[0]!=NULL);
1904 ASSERT(Nodes[1]!=NULL);
1905
1906 /* one node is side node -> no father edge */
1907 if (NTYPE(Nodes[0])==SIDE_NODE || NTYPE(Nodes[1])==SIDE_NODE) return(NULL);
1908
1909 /* both nodes are side nodes -> no father edge */
1910 if (NTYPE(Nodes[0])==MID_NODE && NTYPE(Nodes[1])==MID_NODE) return(NULL);
1911
1912 for (pos0=0; pos0<MAX_SIDE_NODES; pos0++) {
1913 if (SideNodes[pos0] == Nodes[0])
1914 break;
1915 }
1916 ASSERT(pos0<MAX_SIDE_NODES);
1917
1918 for (pos1=0; pos1<MAX_SIDE_NODES; pos1++)
1919 if (SideNodes[pos1] == Nodes[1])
1920 break;
1921 ASSERT(pos1<MAX_SIDE_NODES);
1922
1923 switch (NTYPE(Nodes[0]))
1924 {
1925 case (CORNER_NODE) :
1926
1927 ASSERT(pos0<ncorners);
1928 if ( ((pos0+1)%ncorners == pos1) ||
1929 (pos0+ncorners == pos1) )
1930 {
1931 ASSERT(OBJT(NFATHER(SideNodes[(pos0+1)%ncorners])) == NDOBJ);
1932 fatherEdge = GetEdge((NODE *)NFATHER(Nodes[0]),
1933 (NODE *)NFATHER(SideNodes[(pos0+1)%ncorners]));
1934 ASSERT(fatherEdge!=NULL);
1935 }
1936
1937 if ( ((pos0-1+ncorners)%ncorners == pos1) ||
1938 ((pos0-1+ncorners)%ncorners+ncorners == pos1) )
1939 {
1940 ASSERT(OBJT(NFATHER(SideNodes[(pos0-1+ncorners)%ncorners])) == NDOBJ);
1941 fatherEdge = GetEdge((NODE *)NFATHER(Nodes[0]),
1942 (NODE *)NFATHER(SideNodes[(pos0-1+ncorners)%ncorners]));
1943 ASSERT(fatherEdge!=NULL);
1944 }
1945
1946 break;
1947
1948 case (MID_NODE) :
1949
1950 ASSERT(pos0>=ncorners);
1951 ASSERT(pos0<2*ncorners);
1952
1953 if ((pos0+1)%ncorners == pos1)
1954 {
1955 ASSERT(OBJT(NFATHER(SideNodes[pos0%ncorners])) == NDOBJ);
1956 fatherEdge = GetEdge((NODE *)NFATHER(SideNodes[pos0%ncorners]),
1957 (NODE *)NFATHER(Nodes[1]));
1958 ASSERT(fatherEdge!=NULL);
1959 }
1960
1961 if (pos0%ncorners == pos1)
1962 {
1963 ASSERT(OBJT(NFATHER(SideNodes[(pos0+1)%ncorners])) == NDOBJ);
1964 fatherEdge = GetEdge((NODE *)NFATHER(SideNodes[(pos0+1)%ncorners]),
1965 (NODE *)NFATHER(Nodes[1]));
1966 ASSERT(fatherEdge!=NULL);
1967 }
1968
1969 break;
1970
1971 case (SIDE_NODE) :
1972
1973 /* this edge has no father edge */
1974 fatherEdge = NULL;
1975 break;
1976
1977 default :
1978 assert(0);
1979 break;
1980 }
1981
1982 IFDEBUG(dddif,1)
1983 INT i;
1984 EDGE* edge0, *edge1;
1985
1986 edge0 = edge1 = NULL;
1987
1988 /* test whether theEdge lies above fatherEdge */
1989 if (fatherEdge!=NULL)
1990 {
1991 if (MIDNODE(fatherEdge)!=NULL)
1992 {
1993 edge0 = GetEdge(MIDNODE(fatherEdge),SONNODE(NBNODE(LINK0(fatherEdge))));
1994 edge1 = GetEdge(MIDNODE(fatherEdge),SONNODE(NBNODE(LINK1(fatherEdge))));
1995 }
1996 else
1997 edge0 = GetEdge(SONNODE(NBNODE(LINK0(fatherEdge))),
1998 SONNODE(NBNODE(LINK1(fatherEdge))));
1999
2000 IFDEBUG(dddif,5)
2001 UserWriteF("fatherEdge=%x theEdge=%x edge0=%x edge1=%x\n",
2002 fatherEdge,theEdge,edge0,edge1);
2003 UserWriteF("Nodes[0]=%d Nodes[1]=%d\n",ID(Nodes[0]),ID(Nodes[1]));
2004
2005 UserWriteF("SideNodes\n");
2006 for (i=0; i<MAX_SIDE_NODES; i++) UserWriteF(" %5d",i);
2007 UserWriteF("\n");
2008 for (i=0; i<MAX_SIDE_NODES; i++)
2009 if (SideNodes[i]!=NULL) UserWriteF(" %5d",ID(SideNodes[i]));
2010 UserWriteF("\n");
2011 ENDDEBUG
2012
2013 assert(edge0==theEdge || edge1==theEdge);
2014 }
2015 ENDDEBUG
2016
2017 return(fatherEdge);
2018 }
2019 #endif
2020
2021 /****************************************************************************/
2022 /** \brief Return pointer to edge if it exists
2023
2024 * @param from - starting node of edge
2025 * @param to - end node of edge
2026
2027 This function returns the pointer to the specified edge if it exists.
2028
2029 @return <ul>
2030 <li> pointer to specified object </li>
2031 <li> NULL if not found </li>
2032 </ul> */
2033 /****************************************************************************/
2034
GetEdge(const NODE * from,const NODE * to)2035 EDGE * NS_DIM_PREFIX GetEdge (const NODE *from, const NODE *to)
2036 {
2037 LINK *pl;
2038
2039 /* run through neighbor list */
2040 for (pl=START(from); pl!=NULL; pl = NEXT(pl))
2041 if (NBNODE(pl)==to)
2042 return(MYEDGE(pl));
2043
2044 /* return not found */
2045 return(NULL);
2046 }
2047
2048 /****************************************************************************/
2049 /** \brief Return pointer to a new edge structure
2050
2051 * @param theGrid - grid where vertex should be inserted
2052 * @param theElement - pointer to element
2053 * @param edge - number of edge
2054 * @param with_vector - also create vector for edge (true/false)
2055
2056 This function returns a pointer to a new edge structure.
2057
2058 @return <ul>
2059 <li> pointer to requested object </li>
2060 <li> NULL if out of memory </li>
2061 </ul> */
2062 /****************************************************************************/
2063
2064 #ifndef ModelP
2065 static
2066 #endif
2067 EDGE *
2068 #ifdef ModelP
2069 NS_DIM_PREFIX
2070 #endif
CreateEdge(GRID * theGrid,ELEMENT * theElement,INT edge,bool with_vector)2071 CreateEdge (GRID *theGrid, ELEMENT *theElement, INT edge, bool with_vector)
2072 {
2073 ELEMENT *theFather;
2074 EDGE *pe,*father_edge;
2075 NODE *from,*to,*n1,*n2;
2076 LINK *link0,*link1;
2077 VECTOR *pv;
2078 #ifdef __THREEDIM__
2079 VERTEX *theVertex;
2080 NODE *nbn1,*nbn2,*nbn3,*nbn4;
2081 INT sc,found,side,k,j;
2082 #endif
2083
2084 from = CORNER(theElement,CORNER_OF_EDGE(theElement,edge,0));
2085 to = CORNER(theElement,CORNER_OF_EDGE(theElement,edge,1));
2086
2087 /* check if edge exists already */
2088 if( (pe = GetEdge(from, to)) != NULL ) {
2089 if (NO_OF_ELEM(pe)<NO_OF_ELEM_MAX-1)
2090 INC_NO_OF_ELEM(pe);
2091 else
2092 ASSERT(0);
2093
2094 return(pe);
2095 }
2096
2097 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,EDGEVEC))
2098 pe = (EDGE*)GetMemoryForObject(theGrid->mg,sizeof(EDGE),EDOBJ);
2099 else
2100 pe = (EDGE*)GetMemoryForObject(theGrid->mg,sizeof(EDGE)-sizeof(VECTOR*),EDOBJ);
2101 if (pe==NULL) return(NULL);
2102
2103 /* initialize data */
2104 link0 = LINK0(pe);
2105 link1 = LINK1(pe);
2106 SETOBJT(pe,EDOBJ);
2107 SETLOFFSET(link0,0);
2108 #ifdef _DEBUG_CW_
2109 SETOBJT(link1,LIOBJ);
2110 #endif
2111 SETLOFFSET(link1,1);
2112
2113 pe->id = (theGrid->mg->edgeIdCounter)++;
2114
2115 SETLEVEL(pe,GLEVEL(theGrid));
2116 #ifdef ModelP
2117 DDD_AttrSet(PARHDR(pe), GRID_ATTR(theGrid));
2118 /* SETPRIO(pe,PrioMaster); */
2119 #endif
2120 #ifdef IDENT_ONLY_NEW
2121 if (GET_IDENT_MODE() == IDENT_ON)
2122 SETNEW_EDIDENT(pe,1);
2123 #endif
2124
2125 UGM_CDBG(pe,
2126 UserWriteF(PFMT "create edge=" EDID_FMTX " from=" ID_FMTX "tf=%d to=" ID_FMTX "tt=%d"
2127 "elem=" EID_FMTX "edge=%d\n",
2128 theGrid->ppifContext().me(),EDID_PRTX(pe),ID_PRTX(from),NTYPE(from),ID_PRTX(to),NTYPE(to),
2129 EID_PRTX(theElement),edge);
2130 if (0)
2131 UserWriteF(PFMT "nfatherf=" ID_FMTX "nfathert=" ID_FMTX " fatheredge=" EDID_FMTX "\n",
2132 theGrid->ppifContext().me(),ID_PRTX((NODE*)NFATHER(from)),ID_PRTX((NODE*)NFATHER(to)),
2133 EDID_PRTX(GetEdge((NODE*)NFATHER(from),(NODE*)NFATHER(to))));)
2134
2135 NBNODE(link0) = to;
2136 NBNODE(link1) = from;
2137 SET_NO_OF_ELEM(pe,1);
2138 SETEDGENEW(pe,1);
2139
2140 /* set edge-subdomain from topological information with respect to father-element */
2141 SETEDSUBDOM(pe,SUBDOMAIN(theElement));
2142 theFather = EFATHER(theElement);
2143 if (theFather!=NULL)
2144 {
2145 SETEDSUBDOM(pe,SUBDOMAIN(theFather));
2146 if (NTYPE(from)<NTYPE(to))
2147 {
2148 n1 = from;
2149 n2 = to;
2150 }
2151 else
2152 {
2153 n1 = to;
2154 n2 = from;
2155 }
2156 switch(NTYPE(n1)|(NTYPE(n2)<<4))
2157 {
2158 #ifdef __TWODIM__
2159 case (CORNER_NODE | (CORNER_NODE<<4)) :
2160 father_edge = GetEdge(NFATHER(n1),NFATHER(n2));
2161 if (father_edge!=NULL) SETEDSUBDOM(pe,EDSUBDOM(father_edge));
2162 break;
2163 case (CORNER_NODE | (MID_NODE<<4)) :
2164 father_edge = NFATHEREDGE(n2);
2165 #ifdef ModelP
2166 if (father_edge==NULL)
2167 {
2168 /* TODO: check this after priority set:
2169 assert( GHOST(n1) || GHOST(n2) ); */
2170 break;
2171 }
2172 #endif
2173 assert(father_edge!=NULL);
2174 if (NBNODE(LINK0(father_edge))==NFATHER(n1) || NBNODE(LINK1(father_edge))==NFATHER(n1)) SETEDSUBDOM(pe,EDSUBDOM(father_edge));
2175 break;
2176 #endif
2177 #ifdef __THREEDIM__
2178 case (CORNER_NODE | (CORNER_NODE<<4)) :
2179 father_edge = GetEdge(NFATHER(n1),NFATHER(n2));
2180 if (father_edge!=NULL) SETEDSUBDOM(pe,EDSUBDOM(father_edge));
2181 else
2182 {
2183 /* do fathers of n1, n2 lies on a side (of the father) which has BNDS? */
2184 for (j=0; j<SIDES_OF_ELEM(theFather); j++)
2185 {
2186 found=0;
2187 for (k=0; k<CORNERS_OF_SIDE(theFather,j); k++)
2188 {
2189 sc = CORNER_OF_SIDE(theFather,j,k);
2190 if (CORNER(theFather,sc)==NFATHER(n1) || CORNER(theFather,sc)==NFATHER(n2)) found++;
2191 }
2192 if (found==2 && (OBJT(theFather)==BEOBJ) && SIDE_ON_BND(theFather,j))
2193 {
2194 SETEDSUBDOM(pe,0);
2195 break;
2196 }
2197 }
2198 }
2199 break;
2200
2201 case (CORNER_NODE | (MID_NODE<<4)) :
2202 father_edge = NFATHEREDGE(n2);
2203 assert(father_edge!=NULL);
2204 nbn1 = NBNODE(LINK0(father_edge));
2205 nbn2 = NBNODE(LINK1(father_edge));
2206 if (nbn1==NFATHER(n1) || nbn2==NFATHER(n1)) SETEDSUBDOM(pe,EDSUBDOM(father_edge));
2207 else
2208 {
2209 /* do all nodes n1, nbn1, nbn2 ly on the same side of father? */
2210 side=-1;
2211 for (j=0; j<SIDES_OF_ELEM(theFather); j++)
2212 {
2213 found=0;
2214 for (k=0; k<CORNERS_OF_SIDE(theFather,j); k++)
2215 {
2216 sc = CORNER_OF_SIDE(theFather,j,k);
2217 if (CORNER(theFather,sc)==NFATHER(n1) || CORNER(theFather,sc)==nbn1 || CORNER(theFather,sc)==nbn2) found++;
2218 }
2219 if (found==3)
2220 {
2221 side = j;
2222 break;
2223 }
2224 }
2225 if (side>=0 && (OBJT(theFather)==BEOBJ) && SIDE_ON_BND(theFather,side)) SETEDSUBDOM(pe,0);
2226 }
2227 break;
2228
2229 case (MID_NODE | (MID_NODE<<4)) :
2230 father_edge = NFATHEREDGE(n1);
2231 assert(father_edge!=NULL);
2232 nbn1 = NBNODE(LINK0(father_edge));
2233 nbn2 = NBNODE(LINK1(father_edge));
2234 father_edge = NFATHEREDGE(n2);
2235 assert(father_edge!=NULL);
2236 nbn3 = NBNODE(LINK0(father_edge));
2237 nbn4 = NBNODE(LINK1(father_edge));
2238
2239 /* do all nodes nbn1, nbn2, nbn3, nbn4 ly on the same side of father? */
2240 side=-1;
2241 for (j=0; j<SIDES_OF_ELEM(theFather); j++)
2242 {
2243 found=0;
2244 for (k=0; k<CORNERS_OF_SIDE(theFather,j); k++)
2245 {
2246 sc = CORNER_OF_SIDE(theFather,j,k);
2247 if (CORNER(theFather,sc)==nbn1) found++;
2248 if (CORNER(theFather,sc)==nbn2) found++;
2249 if (CORNER(theFather,sc)==nbn3) found++;
2250 if (CORNER(theFather,sc)==nbn4) found++;
2251 }
2252 if (found==4)
2253 {
2254 side = j;
2255 break;
2256 }
2257 }
2258 if (side>=0 && (OBJT(theFather)==BEOBJ) && SIDE_ON_BND(theFather,side)) SETEDSUBDOM(pe,0);
2259 break;
2260
2261 case (CORNER_NODE | (SIDE_NODE<<4)) :
2262 theVertex = MYVERTEX(n2);
2263 if (VFATHER(theVertex) == theFather)
2264 side = ONSIDE(theVertex);
2265 else
2266 side = ONNBSIDE(theVertex);
2267 if ((OBJT(theFather)==BEOBJ) && SIDE_ON_BND(theFather,side))
2268 for (k=0; k<CORNERS_OF_SIDE(theFather,side); k++)
2269 if (CORNER(theFather,CORNER_OF_SIDE(theFather,side,k))==NFATHER(n1))
2270 {
2271 SETEDSUBDOM(pe,0);
2272 break;
2273 }
2274 break;
2275
2276 case (MID_NODE | (SIDE_NODE<<4)) :
2277 theVertex = MYVERTEX(n2);
2278 if (VFATHER(theVertex) == theFather)
2279 side = ONSIDE(theVertex);
2280 else
2281 side = ONNBSIDE(theVertex);
2282 if ((OBJT(theFather)==BEOBJ) && SIDE_ON_BND(theFather,side))
2283 {
2284 found=0;
2285 father_edge = NFATHEREDGE(n1);
2286 assert(father_edge!=NULL);
2287 nbn1 = NBNODE(LINK0(father_edge));
2288 nbn2 = NBNODE(LINK1(father_edge));
2289 for (k=0; k<CORNERS_OF_SIDE(theFather,side); k++)
2290 {
2291 if (CORNER(theFather,CORNER_OF_SIDE(theFather,side,k))==nbn1 || CORNER(theFather,CORNER_OF_SIDE(theFather,side,k))==nbn2)
2292 found++;
2293 }
2294 if (found==2) SETEDSUBDOM(pe,0);
2295 }
2296 break;
2297 #endif
2298 } /* end switch */
2299 } /* end (theFather!=NULL) */
2300
2301 /* create vector if */
2302 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,EDGEVEC))
2303 {
2304 if (with_vector) {
2305 if (CreateVector (theGrid,EDGEVEC,(GEOM_OBJECT *)pe,&pv)) {
2306 DisposeEdge (theGrid,pe);
2307 return (NULL);
2308 }
2309 EDVECTOR(pe) = pv;
2310 }
2311 else
2312 EDVECTOR(pe) = NULL;
2313 }
2314
2315 /* put in neighbor lists */
2316 NEXT(link0) = START(from);
2317 START(from) = link0;
2318 NEXT(link1) = START(to);
2319 START(to) = link1;
2320
2321 /* counters */
2322 NE(theGrid)++;
2323
2324 /* return ok */
2325 return(pe);
2326 }
2327
2328 /****************************************************************************/
2329 /** \brief Return pointer to link if it exists
2330
2331 * @param from - starting node of link
2332 * @param to - end node of link
2333
2334 This function returns the pointer to the specified link if it exists.
2335
2336 @return <ul>
2337 <li> pointer to specified link </li>
2338 <li> NULL if not found. </li>
2339 </ul> */
2340 /****************************************************************************/
2341
GetLink(const NODE * from,const NODE * to)2342 LINK *GetLink (const NODE *from, const NODE *to)
2343 {
2344 LINK *pl;
2345
2346 /* run through neighbor list */
2347 for (pl=START(from); pl!=NULL; pl = NEXT(pl))
2348 if (NBNODE(pl)==to)
2349 return(pl);
2350
2351 /* return not found */
2352 return(NULL);
2353 }
2354
2355 /****************************************************************************/
2356 /** \brief Return a pointer to a new element structure
2357
2358 * @param theGrid - grid structure to extend
2359 * @param tag - the element type
2360 * @param objtype - inner element (IEOBJ) or boundary element (BEOBJ)
2361 * @param nodes - list of corner nodes in reference numbering
2362 * @param Father - pointer to father element (NULL on base level)
2363 * @param with_vector -
2364
2365 This function creates and initializes a new element and returns a pointer to it.
2366
2367 @return <ul>
2368 <li> pointer to requested object </li>
2369 <li> NULL if out of memory </li>
2370 </ul> */
2371 /****************************************************************************/
2372
CreateElement(GRID * theGrid,INT tag,INT objtype,NODE ** nodes,ELEMENT * Father,bool with_vector)2373 ELEMENT * NS_DIM_PREFIX CreateElement (GRID *theGrid, INT tag, INT objtype, NODE **nodes,
2374 ELEMENT *Father, bool with_vector)
2375 {
2376 ELEMENT *pe;
2377 INT i,s_id;
2378 VECTOR *pv;
2379
2380 if (objtype == IEOBJ)
2381 pe = (ELEMENT*)GetMemoryForObject(MYMG(theGrid),INNER_SIZE_TAG(tag),
2382 MAPPED_INNER_OBJT_TAG(tag));
2383 else if (objtype == BEOBJ)
2384 pe = (ELEMENT*)GetMemoryForObject(MYMG(theGrid),BND_SIZE_TAG(tag),
2385 MAPPED_BND_OBJT_TAG(tag));
2386 else
2387 std::abort();
2388
2389 if (pe==NULL) return(NULL);
2390
2391 /* initialize data */
2392 SETNEWEL(pe,1);
2393 SETOBJT(pe,objtype);
2394 SETTAG(pe,tag);
2395 SETLEVEL(pe,theGrid->level);
2396 #ifdef ModelP
2397 DDD_AttrSet(PARHDRE(pe),GRID_ATTR(theGrid));
2398 /* SETEPRIO(theGrid->dddContext(), pe,PrioMaster); */
2399 PARTITION(pe) = theGrid->ppifContext().me();
2400 #endif
2401 SETEBUILDCON(pe,1);
2402 ID(pe) = (theGrid->mg->elemIdCounter)++;
2403
2404 /* subdomain id */
2405 s_id = (Father != NULL) ? SUBDOMAIN(Father) : 0;
2406 SETSUBDOMAIN(pe,s_id);
2407
2408 #ifdef __CENTERNODE__
2409 SET_CENTERNODE(pe,NULL);
2410 #endif
2411
2412 SET_EFATHER(pe,Father);
2413
2414 /* set corner nodes */
2415 for (i=0; i<CORNERS_OF_ELEM(pe); i++)
2416 SET_CORNER(pe,i,nodes[i]);
2417
2418 /* create edges */
2419 for (i=0; i<EDGES_OF_ELEM(pe); i++)
2420 if (CreateEdge (theGrid,pe,i,with_vector) == NULL) {
2421 DisposeElement(theGrid,pe,true);
2422 return(NULL);
2423 }
2424
2425 UGM_CDBG(pe,
2426 UserWriteF(PFMT "create elem=" EID_FMTX,
2427 theGrid->ppifContext().me(),EID_PRTX(pe));
2428 for (i=0; i<CORNERS_OF_ELEM(pe); i++)
2429 UserWriteF(" n%d=" ID_FMTX, i,ID_PRTX(CORNER(pe,i)));
2430 UserWriteF("\n");
2431 for (i=0; i<EDGES_OF_ELEM(pe); i++)
2432 {
2433 EDGE *theEdge;
2434
2435 theEdge = GetEdge(CORNER_OF_EDGE_PTR(pe,i,0),
2436 CORNER_OF_EDGE_PTR(pe,i,1));
2437 UserWriteF(" e%d=" EDID_FMTX, i,EDID_PRTX(theEdge));
2438 }
2439 UserWriteF("\n");)
2440
2441
2442 /* create element vector if */
2443 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,ELEMVEC))
2444 {
2445 if (with_vector)
2446 {
2447 if (CreateVector (theGrid,ELEMVEC,(GEOM_OBJECT *)pe,&pv))
2448 {
2449 DisposeElement(theGrid,pe,true);
2450 return (NULL);
2451 }
2452 SET_EVECTOR(pe,pv);
2453 }
2454 else
2455 SET_EVECTOR(pe,NULL);
2456 }
2457
2458 /* create side vectors if */
2459 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,SIDEVEC))
2460 {
2461 for (i=0; i<SIDES_OF_ELEM(pe); i++)
2462 if (with_vector)
2463 {
2464 if (CreateSideVector (theGrid,i,(GEOM_OBJECT *)pe,&pv))
2465 {
2466 DisposeElement(theGrid,pe,true);
2467 return (NULL);
2468 }
2469 SET_SVECTOR(pe,i,pv);
2470 }
2471 else
2472 SET_SVECTOR(pe,i,NULL);
2473 }
2474
2475 /* insert in element list */
2476 GRID_LINK_ELEMENT(theGrid,pe,PrioMaster);
2477
2478 if (theGrid->level>0)
2479 {
2480 INT where = PRIO2INDEX(PrioMaster);
2481
2482 #ifndef ModelP
2483 ASSERT(Father != NULL);
2484 #endif
2485 if (Father != NULL)
2486 {
2487 if (SON(Father,where) == NULL)
2488 SET_SON(Father,where,pe);
2489 SETNSONS(Father,NSONS(Father)+1);
2490 }
2491 }
2492
2493 /* return ok */
2494 return(pe);
2495 }
2496
2497 /****************************************************************************/
2498 /** \brief Creates the element sides of son elements
2499
2500 * @param theGrid - grid for which to create
2501 * @param theElement - pointer to a boundary element
2502 * @param side - side id of a side of the element
2503 * @param theSon - pointer to a son element
2504 * @param son_side - side id of a side of the son
2505
2506 This function creates and initializes an element side of a son element.
2507 Here also the side vector (iff at all) is inspected in 'ReinspectSonSideVector'.
2508 The latter function eventually reallocates the vector if its size has changed and
2509 sets the VBUILDCON flag in the vector. The connections of the old vector are
2510 thereby disposed. The refine-module which is calling
2511 'CreateSonElementSide' will finally call 'GridCreateConnection' to reinstall
2512 the connections of the side-vector.
2513
2514 @return <ul>
2515 <li> GM_OK if ok </li>
2516 <li> GM_ERROR when error occured. </li>
2517 </ul> */
2518 /****************************************************************************/
2519
CreateSonElementSide(GRID * theGrid,ELEMENT * theElement,INT side,ELEMENT * theSon,INT son_side)2520 INT NS_DIM_PREFIX CreateSonElementSide (GRID *theGrid, ELEMENT *theElement, INT side,
2521 ELEMENT *theSon, INT son_side)
2522 {
2523 INT n,i;
2524 BNDS *bnds;
2525 BNDP *bndp[MAX_CORNERS_OF_ELEM];
2526 VECTOR *vec;
2527 [[maybe_unused]] EDGE *theEdge;
2528
2529 ASSERT (OBJT(theElement) == BEOBJ);
2530
2531 ASSERT (ELEM_BNDS(theElement,side) != NULL);
2532
2533 /* check if Edges of theElement, which are on the side 'side' have EDSUBDOM 0 */
2534 n = CORNERS_OF_SIDE(theElement,side);
2535 for (i=0; i<n; i++)
2536 {
2537 theEdge = GetEdge(CORNER(theElement,CORNER_OF_SIDE(theElement,side,i)),CORNER(theElement,CORNER_OF_SIDE(theElement,side,(i+1)%n)));
2538 assert(EDSUBDOM(theEdge)==0);
2539 }
2540
2541 n = CORNERS_OF_SIDE(theSon,son_side);
2542 for (i=0; i<n; i++)
2543 {
2544 /* check if vertices of Son lie on boundary */
2545 if (OBJT(MYVERTEX(CORNER(theSon,CORNER_OF_SIDE(theSon,son_side,i))))!=BVOBJ)
2546 {
2547 NODE *theNode;
2548 EDGE *theFatherEdge;
2549 INT t1,t2;
2550
2551 theNode = CORNER(theSon,CORNER_OF_SIDE(theSon,son_side,i));
2552 printf("ID=%d\n",(int)ID(theNode));
2553 switch (NTYPE(theNode))
2554 {
2555 case CORNER_NODE :
2556 printf("NTYPE = CORNER_NODE");
2557 break;
2558 case MID_NODE :
2559 printf(PFMT "el " EID_FMTX " son " EID_FMTX " vertex " VID_FMTX "\n",theGrid->ppifContext().me(),EID_PRTX(theElement),EID_PRTX(theSon),VID_PRTX(MYVERTEX(CORNER(theSon,CORNER_OF_SIDE(theSon,son_side,i)))));
2560 printf(PFMT "NTYPE = MID_NODE\n",theGrid->ppifContext().me());
2561 theFatherEdge = NFATHEREDGE(theNode);
2562 printf(PFMT "EDSUBDOM = %d\n",theGrid->ppifContext().me(),(int)EDSUBDOM(theFatherEdge));
2563 t1 = (OBJT(MYVERTEX(NBNODE(LINK0(theFatherEdge))))==BVOBJ);
2564 t2 = (OBJT(MYVERTEX(NBNODE(LINK1(theFatherEdge))))==BVOBJ);
2565 printf(PFMT "BVOBJ(theFatherEdge): %d %d\n",theGrid->ppifContext().me(),(int)t1,(int)t2);
2566 break;
2567 case SIDE_NODE :
2568 printf("NTYPE = SIDE_NODE");
2569 break;
2570 case CENTER_NODE :
2571 printf("NTYPE = CENTER_NODE");
2572 break;
2573 }
2574 ASSERT(0);
2575 }
2576 bndp[i] = V_BNDP(MYVERTEX(CORNER(theSon,CORNER_OF_SIDE(theSon,son_side,i))));
2577 }
2578 bnds = BNDP_CreateBndS(MGHEAP(MYMG(theGrid)),bndp,n);
2579 if (bnds == NULL)
2580 RETURN(GM_ERROR);
2581 SET_BNDS(theSon,son_side,bnds);
2582
2583 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,SIDEVEC))
2584 {
2585 vec = SVECTOR(theSon,son_side);
2586 ReinspectSonSideVector(theGrid,theSon,son_side,&vec);
2587 SET_SVECTOR(theSon,son_side,vec);
2588 }
2589 #ifdef __TWODIM__
2590 theEdge = GetEdge(CORNER(theSon,CORNER_OF_EDGE(theSon,son_side,0)),
2591 CORNER(theSon,CORNER_OF_EDGE(theSon,son_side,1)));
2592 ASSERT(theEdge != NULL);
2593 SETEDSUBDOM(theEdge,0);
2594 #endif
2595
2596 #ifdef __THREEDIM__
2597 /** \todo is this necessary?
2598 for (i=0; i<EDGES_OF_SIDE(theSon,son_side); i++) {
2599 int k = EDGE_OF_SIDE(theSon,son_side,i);
2600 theEdge = GetEdge(CORNER(theSon,CORNER_OF_EDGE(theSon,k,0)),
2601 CORNER(theSon,CORNER_OF_EDGE(theSon,k,1)));
2602 ASSERT(theEdge != NULL);
2603 SETEDSUBDOM(theEdge,0);
2604 } */
2605 #endif
2606
2607 return(GM_OK);
2608 }
2609
2610 /****************************************************************************/
2611 /** \brief Return pointer to new grid structure
2612
2613 * @param theMG - multigrid structure
2614
2615 This function creates and initialized a new grid structure for top level + 1
2616 and returns a pointer to it.
2617
2618 @return <ul>
2619 <li> pointer to requested object </li>
2620 <li> NULL if out of memory </li>
2621 </ul> */
2622 /****************************************************************************/
2623
CreateNewLevel(MULTIGRID * theMG)2624 GRID * NS_DIM_PREFIX CreateNewLevel (MULTIGRID *theMG)
2625 {
2626 GRID *theGrid;
2627
2628 if (TOPLEVEL(theMG)+1>=MAXLEVEL) return(NULL);
2629 INT l = TOPLEVEL(theMG)+1;
2630
2631 /* allocate grid object */
2632 theGrid = (GRID*)GetMemoryForObject(theMG,sizeof(GRID),GROBJ);
2633 if (theGrid==NULL) return(NULL);
2634
2635 /* fill in data */
2636 CTRL(theGrid) = 0;
2637 SETOBJT(theGrid,GROBJ);
2638 GLEVEL(theGrid) = l;
2639 GATTR(theGrid) = GRID_ATTR(theGrid);
2640 NE(theGrid) = 0;
2641 NC(theGrid) = 0;
2642 /* other counters are init in INIT fcts below */
2643
2644 GSTATUS(theGrid,0);
2645 GRID_INIT_ELEMENT_LIST(theGrid);
2646 GRID_INIT_NODE_LIST(theGrid);
2647 GRID_INIT_VERTEX_LIST(theGrid);
2648 GRID_INIT_VECTOR_LIST(theGrid);
2649 if (l>0)
2650 {
2651 DOWNGRID(theGrid) = GRID_ON_LEVEL(theMG,l-1);
2652 UPGRID(GRID_ON_LEVEL(theMG,l-1)) = theGrid;
2653 UPGRID(theGrid) = NULL;
2654 }
2655 else if (l==0)
2656 {
2657 DOWNGRID(theGrid) = NULL;
2658 UPGRID(theGrid) = NULL;
2659 }
2660 else
2661 {
2662 UPGRID(theGrid) = GRID_ON_LEVEL(theMG,l+1);
2663 DOWNGRID(theGrid) = NULL;
2664 DOWNGRID(GRID_ON_LEVEL(theMG,l+1)) = theGrid;
2665 }
2666 MYMG(theGrid) = theMG;
2667 GRID_ON_LEVEL(theMG,l) = theGrid;
2668 TOPLEVEL(theMG) = l;
2669 CURRENTLEVEL(theMG) = l;
2670
2671 return(theGrid);
2672 }
2673
2674
2675 /****************************************************************************/
2676 /** \brief Create a multigrid environment item
2677
2678 * @param name - name of the multigrid
2679
2680 This function creates a multigrid environment directory.
2681
2682 @return <ul>
2683 <li> pointer to new MULTIGRID </li>
2684 <li> NULL if error occured </li>
2685 </ul> */
2686 /****************************************************************************/
2687
MakeMGItem(const char * name,std::shared_ptr<PPIF::PPIFContext> ppifContext)2688 MULTIGRID * NS_DIM_PREFIX MakeMGItem (const char *name, std::shared_ptr<PPIF::PPIFContext> ppifContext)
2689 {
2690 MULTIGRID *theMG;
2691
2692 if (ChangeEnvDir("/Multigrids") == NULL) return (NULL);
2693 if (strlen(name)>=NAMESIZE || strlen(name)<=1) return (NULL);
2694 theMG = (MULTIGRID *) MakeEnvItem(name,theMGDirID,sizeof(MULTIGRID));
2695 if (theMG == NULL) return(NULL);
2696
2697 new(theMG) multigrid;
2698
2699 #if ModelP
2700 theMG->ppifContext_ = ppifContext;
2701 theMG->dddContext_ = std::make_shared<DDD::DDDContext>(
2702 theMG->ppifContext_,
2703 std::make_shared<DDD_CTRL>()
2704 );
2705
2706 InitDDD(theMG->dddContext());
2707
2708 globalDDDContext(theMG->dddContext_);
2709 #else
2710 theMG->ppifContext_ = std::make_shared<PPIF::PPIFContext>();
2711 #endif
2712
2713 return (theMG);
2714 }
2715
2716 /****************************************************************************/
2717 /** \todo Please doc me!
2718
2719 * @param theMG
2720 * @param FromLevel
2721 * @param ToLevel
2722 * @param mask
2723
2724 DESCRIPTION:
2725
2726 @return
2727 INT
2728 */
2729 /****************************************************************************/
2730
ClearMultiGridUsedFlags(MULTIGRID * theMG,INT FromLevel,INT ToLevel,INT mask)2731 INT NS_DIM_PREFIX ClearMultiGridUsedFlags (MULTIGRID *theMG, INT FromLevel, INT ToLevel, INT mask)
2732 {
2733 int i,level,elem,node,edge,vertex,vector,matrix;
2734 GRID *theGrid;
2735 ELEMENT *theElement;
2736 NODE *theNode;
2737 EDGE *theEdge;
2738 VECTOR *theVector;
2739 MATRIX *theMatrix;
2740
2741 elem = mask & MG_ELEMUSED;
2742 node = mask & MG_NODEUSED;
2743 edge = mask & MG_EDGEUSED;
2744 vertex = mask & MG_VERTEXUSED;
2745 vector = mask & MG_VECTORUSED;
2746 matrix = mask & MG_MATRIXUSED;
2747
2748 for (level=FromLevel; level<=ToLevel; level++)
2749 {
2750 theGrid = GRID_ON_LEVEL(theMG,level);
2751 if (elem || edge)
2752 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL; theElement=SUCCE(theElement))
2753 {
2754 if (elem) SETUSED(theElement,0);
2755 if (edge)
2756 {
2757 for (i=0; i<EDGES_OF_ELEM(theElement); i++)
2758 {
2759 theEdge = GetEdge(CORNER_OF_EDGE_PTR(theElement,i,0),
2760 CORNER_OF_EDGE_PTR(theElement,i,1));
2761 SETUSED(theEdge,0);
2762 }
2763 }
2764 }
2765 if (node || vertex)
2766 {
2767 for (theNode=PFIRSTNODE(theGrid); theNode!=NULL; theNode=SUCCN(theNode))
2768 {
2769 if (node) SETUSED(theNode,0);
2770 if (vertex) SETUSED(MYVERTEX(theNode),0);
2771 }
2772 }
2773 if (vector || matrix)
2774 {
2775 for (theVector=PFIRSTVECTOR(theGrid); theVector!=NULL; theVector=SUCCVC(theVector))
2776 {
2777 if (vector) SETUSED(theVector,0);
2778 if (matrix)
2779 for (theMatrix=VSTART(theVector); theMatrix!=NULL; theMatrix=MNEXT(theMatrix))
2780 SETUSED(theMatrix,0);
2781 }
2782 }
2783
2784 }
2785
2786 return(GM_OK);
2787 }
2788
2789
2790 /****************************************************************************/
2791 /** \brief Find the multigrid environment item with name
2792
2793 * @param name - name of the multigrid to find
2794
2795 This function find the multigrid environment item with `name` and
2796 returns a pointer to the multigrid structure.
2797
2798 @return <ul>
2799 <li> pointer to MULTIGRID </li>
2800 <li> NULL if not found. </li>
2801 </ul> */
2802 /****************************************************************************/
2803
GetMultigrid(const char * name)2804 MULTIGRID * NS_DIM_PREFIX GetMultigrid (const char *name)
2805 {
2806 return ((MULTIGRID *) SearchEnv(name,"/Multigrids",
2807 theMGDirID,theMGRootDirID));
2808 }
2809
2810 /****************************************************************************/
2811 /** \brief Return a pointer to the first multigrid
2812
2813 This function returns a pointer to the first multigrid in the /Multigrids
2814 directory.
2815
2816 @return <ul>
2817 <li> pointer to MULTIGRID </li>
2818 <li> NULL if not found. </li>
2819 </ul> */
2820 /****************************************************************************/
2821
GetFirstMultigrid()2822 MULTIGRID * NS_DIM_PREFIX GetFirstMultigrid ()
2823 {
2824 ENVDIR *theMGRootDir;
2825 MULTIGRID *theMG;
2826
2827 theMGRootDir = ChangeEnvDir("/Multigrids");
2828
2829 assert (theMGRootDir!=NULL);
2830
2831 theMG = (MULTIGRID *) ENVDIR_DOWN(theMGRootDir);
2832
2833 if (theMG != NULL)
2834 if (InitElementTypes(theMG) != GM_OK) {
2835 PrintErrorMessage('E',"GetFirstMultigrid",
2836 "error in InitElementTypes");
2837 return(NULL);
2838 }
2839
2840 return (theMG);
2841 }
2842
2843 /****************************************************************************/
2844 /** \brief Return a pointer to the next multigrid
2845
2846 * @param theMG - multigrid structure
2847
2848 This function returns a pointer to the next multigrid in the /Multigrids
2849 directory.
2850
2851 @return <ul>
2852 <li> pointer to MULTIGRID </li>
2853 <li> NULL if not found. </li>
2854 </ul> */
2855 /****************************************************************************/
2856
GetNextMultigrid(const MULTIGRID * theMG)2857 MULTIGRID * NS_DIM_PREFIX GetNextMultigrid (const MULTIGRID *theMG)
2858 {
2859 MULTIGRID *MG;
2860
2861 MG = (MULTIGRID *) NEXT_ENVITEM(theMG);
2862
2863 if (MG != NULL)
2864 if (InitElementTypes(MG)!=GM_OK) {
2865 PrintErrorMessage('E',"GetNextMultigrid",
2866 "error in InitElementTypes");
2867 return(NULL);
2868 }
2869
2870 return (MG);
2871 }
2872
2873 /****************************************************************************/
2874 /** \brief Return a pointer to new multigrid structure
2875
2876 * @param MultigridName - name of multigrid
2877 * @param domain - name of domain description from environment
2878 * @param problem - name of problem description from environment
2879 * @param format - name of format description from environment
2880 * @param optimizedIE - allocate NodeElementList
2881
2882 This function creates and initializes a new multigrid structure including
2883 allocation of heap, combining the domain and the boundary conditions
2884 and creation of the fixed corners of the domain.
2885
2886 @return <ul>
2887 <li> pointer to new object </li>
2888 <li> NULL if an error occured. </li>
2889 </ul> */
2890 /****************************************************************************/
2891
CreateMultiGrid(char * MultigridName,char * BndValProblem,const char * format,INT optimizedIE,INT insertMesh,std::shared_ptr<PPIF::PPIFContext> ppifContext)2892 MULTIGRID * NS_DIM_PREFIX CreateMultiGrid (char *MultigridName, char *BndValProblem,
2893 const char *format, INT optimizedIE, INT insertMesh,
2894 std::shared_ptr<PPIF::PPIFContext> ppifContext)
2895 {
2896 HEAP *theHeap;
2897 MULTIGRID *theMG;
2898 INT i;
2899 BVP *theBVP;
2900 BVP_DESC *theBVPDesc;
2901 MESH mesh;
2902 INT MarkKey;
2903
2904 if (not ppifContext)
2905 ppifContext = std::make_shared<PPIF::PPIFContext>();
2906
2907 std::unique_ptr<FORMAT> theFormat = CreateFormat();
2908 if (theFormat==NULL)
2909 {
2910 PrintErrorMessage('E',"CreateMultiGrid","format not found");
2911 return(NULL);
2912 }
2913
2914
2915 /* allocate multigrid envitem */
2916 theMG = MakeMGItem(MultigridName, ppifContext);
2917 if (theMG==NULL) return(NULL);
2918 MGFORMAT(theMG) = std::move(theFormat);
2919 if (InitElementTypes(theMG)!=GM_OK)
2920 {
2921 PrintErrorMessage('E',"CreateMultiGrid","error in InitElementTypes");
2922 return(NULL);
2923 }
2924
2925 /* allocate the heap */
2926 /* When using the system heap: allocate just enough memory for the actual bookkeeping data structure */
2927 theHeap = NewHeap(SIMPLE_HEAP, sizeof(HEAP), malloc(sizeof(HEAP)));
2928 if (theHeap==NULL)
2929 {
2930 UserWriteF("CreateMultiGrid: cannot allocate %ld bytes\n", sizeof(HEAP));
2931 PrintErrorMessage('E', "CreateMultiGrid","Cannot allocate heap!");
2932
2933 DisposeMultiGrid(theMG);
2934 return(NULL);
2935 }
2936
2937 /* mark temp memory here, release it after coarse grid construction in FixCoarseGrid */
2938 MarkTmpMem(theHeap,&MarkKey);
2939 MG_MARK_KEY(theMG) = MarkKey;
2940
2941 if (insertMesh)
2942 theBVP = BVP_Init(BndValProblem,theHeap,&mesh,MarkKey);
2943 else
2944 theBVP = BVP_Init(BndValProblem,theHeap,NULL,MarkKey);
2945 if (theBVP==NULL)
2946 {
2947 PrintErrorMessage('E',"CreateMultiGrid","BVP not found");
2948 return(NULL);
2949 }
2950 if (BVP_SetBVPDesc(theBVP,&theMG->theBVPD))
2951 {
2952 PrintErrorMessage('E',"CreateMultiGrid","BVP not evaluated");
2953 return(NULL);
2954 }
2955 theBVPDesc = MG_BVPD(theMG);
2956
2957 /* 1: general user data space */
2958 // As we are using this version only with DUNE, we will never have UG user data
2959 /* 2: user heap */
2960 // As we are using this version only with DUNE, we will never need the user heap
2961
2962 /* fill multigrid structure */
2963 theMG->status = 0;
2964 MG_COARSE_FIXED(theMG) = 0;
2965 theMG->vertIdCounter = 0;
2966 theMG->nodeIdCounter = 0;
2967 theMG->elemIdCounter = 0;
2968 theMG->edgeIdCounter = 0;
2969 #ifndef ModelP
2970 theMG->vectorIdCounter = 0;
2971 #endif
2972 theMG->topLevel = -1;
2973 MG_BVP(theMG) = theBVP;
2974 MG_NPROPERTY(theMG) = BVPD_NSUBDOM(theBVPDesc);
2975 RESETMGSTATUS(theMG);
2976
2977 theMG->theHeap = theHeap;
2978 for (i=0; i<MAXLEVEL; i++)
2979 GRID_ON_LEVEL(theMG,i) = NULL;
2980
2981 /* allocate level 0 grid */
2982 if (CreateNewLevel(theMG)==NULL)
2983 {
2984 DisposeMultiGrid(theMG);
2985 return(NULL);
2986 }
2987
2988 /* allocate predefined mesh, e. g. corner vertices pointers */
2989 if (insertMesh)
2990 {
2991 #ifdef ModelP
2992 if (theMG->ppifContext().isMaster())
2993 {
2994 #endif
2995 if (InsertMesh(theMG,&mesh))
2996 {
2997 DisposeMultiGrid(theMG);
2998 return(NULL);
2999 }
3000 #ifdef ModelP
3001 }
3002 #endif
3003
3004 ASSERT(mesh.mesh_status!=MESHSTAT_NOTINIT);
3005 if (mesh.mesh_status==MESHSTAT_MESH)
3006 if (FixCoarseGrid(theMG))
3007 {
3008 DisposeMultiGrid(theMG);
3009 return(NULL);
3010 }
3011 }
3012 /* return ok */
3013 return(theMG);
3014 }
3015
3016 /****************************************************************************/
3017 /** \brief Remove edge from the data structure
3018
3019 * @param theGrid - grid to remove from
3020 * @param theEdge - edge to remove
3021
3022 This function remove an edge from the data structure including its
3023 vector (if one) and inserts them into the free list.
3024
3025 @return <ul>
3026 <li> 0 if ok </li>
3027 <li> 1 if an error occured. </li>
3028 </ul> */
3029 /****************************************************************************/
3030
DisposeEdge(GRID * theGrid,EDGE * theEdge)3031 static INT DisposeEdge (GRID *theGrid, EDGE *theEdge)
3032 {
3033 LINK *link0,*link1,*pl;
3034 NODE *from,*to;
3035 INT found;
3036
3037 /* reconstruct data */
3038 link0 = LINK0(theEdge);
3039 link1 = LINK1(theEdge);
3040 from = NBNODE(link1);
3041 to = NBNODE(link0);
3042 found = 0;
3043
3044 /* delete link0 in from vertex */
3045 if (START(from)==link0)
3046 {
3047 START(from) = NEXT(link0);
3048 found++;
3049 }
3050 else
3051 {
3052 for (pl=START(from); pl!=NULL; pl = NEXT(pl))
3053 {
3054 if (NEXT(pl)==link0)
3055 {
3056 NEXT(pl) = NEXT(link0);
3057 found++;
3058 break;
3059 }
3060 }
3061 }
3062
3063 /* delete link1 in to vertex */
3064 if (START(to)==link1)
3065 {
3066 START(to) = NEXT(link1);
3067 found++;
3068 }
3069 else
3070 {
3071 for (pl=START(to); pl!=NULL; pl = NEXT(pl))
3072 {
3073 if (NEXT(pl)==link1)
3074 {
3075 NEXT(pl) = NEXT(link1);
3076 found++;
3077 break;
3078 }
3079 }
3080 }
3081
3082 /* reset pointer of midnode to edge */
3083 if (MIDNODE(theEdge) != NULL)
3084 SETNFATHER(MIDNODE(theEdge),NULL);
3085
3086 /* dispose vector and its matrices from edge-vector if */
3087 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,EDGEVEC))
3088 {
3089 if (DisposeVector (theGrid,EDVECTOR(theEdge)))
3090 RETURN(1);
3091 PutFreeObject(theGrid->mg,theEdge,sizeof(EDGE),EDOBJ);
3092 }
3093 else
3094 PutFreeObject(theGrid->mg,theEdge,sizeof(EDGE)-sizeof(VECTOR*),EDOBJ);
3095
3096 /* check error condition */
3097 if (found!=2) RETURN(1);
3098
3099 /* return ok */
3100 NE(theGrid)--;
3101 return(0);
3102 }
3103
3104 /****************************************************************************/
3105 /** \brief Remove node including its edges from the data structure
3106
3107 * @param theGrid - grid to remove from
3108 * @param theNode - node to remove
3109
3110 This function removes node including its edges and vector (if one)
3111 from the data structure and inserts all objects into the free list.
3112
3113 @return <ul>
3114 <li> 0 if ok </li>
3115 <li> 1 when error occured. </li>
3116 </ul> */
3117 /****************************************************************************/
3118
DisposeNode(GRID * theGrid,NODE * theNode)3119 INT NS_DIM_PREFIX DisposeNode (GRID *theGrid, NODE *theNode)
3120 {
3121 VERTEX *theVertex;
3122 GEOM_OBJECT *father;
3123 INT size;
3124
3125 /* call DisposeElement first! */
3126 assert(START(theNode) == NULL);
3127 #ifdef ModelP
3128 if (SONNODE(theNode) != NULL)
3129 {
3130 SETNFATHER(SONNODE(theNode),NULL);
3131 }
3132 #else
3133 assert(SONNODE(theNode) == NULL);
3134 #endif
3135
3136 /* remove node from node list */
3137 GRID_UNLINK_NODE(theGrid,theNode);
3138
3139 theVertex = MYVERTEX(theNode);
3140 father = (GEOM_OBJECT *)NFATHER(theNode);
3141 if (father != NULL)
3142 {
3143 switch (NTYPE(theNode))
3144 {
3145
3146 case (CORNER_NODE) :
3147 ASSERT(OBJT(father) == NDOBJ);
3148 SONNODE((NODE *)father) = NULL;
3149 #ifdef TOPNODE
3150 if (theVertex != NULL)
3151 TOPNODE(theVertex) = (NODE *)father;
3152 #endif
3153 break;
3154
3155 case (MID_NODE) :
3156 ASSERT(OBJT(father) == EDOBJ);
3157 MIDNODE((EDGE *)father) = NULL;
3158 break;
3159
3160 #ifdef __CENTERNODE__
3161 case (CENTER_NODE) :
3162 ASSERT(OBJT(father)==IEOBJ || OBJT(father)==BEOBJ);
3163 SET_CENTERNODE((ELEMENT *)father,NULL);
3164 break;
3165 #endif
3166
3167 default :
3168 ASSERT(0);
3169 break;
3170 }
3171 }
3172
3173 /** \todo delete old vertex handling */
3174 if (0)
3175 if (theVertex != NULL)
3176 {
3177 #ifdef ModelP
3178 /* vertices have to be linked and unlinked */
3179 /* relative to the level they are created for */
3180 INT levelofvertex = LEVEL(theVertex);
3181 MULTIGRID *MG = MYMG(theGrid);
3182 GRID *GridOfVertex = GRID_ON_LEVEL(MG,levelofvertex);
3183
3184 if (SONNODE(theNode) == NULL)
3185 DisposeVertex(GridOfVertex,theVertex);
3186 #else
3187 DisposeVertex(theGrid,theVertex);
3188 #endif
3189 }
3190
3191 if (NOOFNODE(theVertex)<1)
3192 RETURN(GM_ERROR);
3193 if (NOOFNODE(theVertex)==1)
3194 DisposeVertex(theGrid,theVertex);
3195 else
3196 DECNOOFNODE(theVertex);
3197
3198 #ifdef ModelP
3199 /* free message buffer */
3200 theNode->message_buffer_free();
3201 #endif
3202
3203 /* dispose vector and its matrices from node-vector */
3204 size = sizeof(NODE);
3205 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,NODEVEC))
3206 {
3207 if (DisposeVector (theGrid,NVECTOR(theNode)))
3208 RETURN(1);
3209 }
3210 else
3211 size -= sizeof(VECTOR *);
3212 PutFreeObject(theGrid->mg,theNode,size,NDOBJ);
3213
3214 /* return ok */
3215 return(0);
3216 }
3217
3218 /****************************************************************************/
3219 /** \brief Remove vertex from the data structure
3220
3221 * @param theGrid - grid to remove from
3222 * @param theVertex - vertex to remove
3223
3224 This function removes a vertex from the data structure
3225 and puts it into the free list.
3226
3227 @return <ul>
3228 <li> 0 if ok </li>
3229 <li> 1 no valid object number </li>
3230 </ul> */
3231 /****************************************************************************/
3232
DisposeVertex(GRID * theGrid,VERTEX * theVertex)3233 static INT DisposeVertex (GRID *theGrid, VERTEX *theVertex)
3234 {
3235 // The following call to HEAPFAULT triggers a failing assertion in some
3236 // distributed settings. I don't know whether this is the sign of a hidden bug
3237 // somewhere, or whether HEAPFAULT is an obsolete left-over of UG3's
3238 // hand-written manual memory heap.
3239 //HEAPFAULT(theVertex);
3240
3241 PRINTDEBUG(gm,1,(PFMT "DisposeVertex(): Gridlevel=%d theVertex=" VID_FMTX "\n",
3242 theGrid->ppifContext().me(),GLEVEL(theGrid),VID_PRTX(theVertex)));
3243
3244 theGrid = GRID_ON_LEVEL(MYMG(theGrid),LEVEL(theVertex));
3245
3246 /* remove vertex from vertex list */
3247 GRID_UNLINK_VERTEX(theGrid,theVertex);
3248
3249 if( OBJT(theVertex) == BVOBJ )
3250 {
3251 BNDP_Dispose(MGHEAP(MYMG(theGrid)),V_BNDP(theVertex));
3252 PutFreeObject(MYMG(theGrid),theVertex,sizeof(struct bvertex),BVOBJ);
3253 }
3254 else
3255 PutFreeObject(MYMG(theGrid),theVertex,sizeof(struct ivertex),IVOBJ);
3256
3257 return(0);
3258 }
3259
3260 /****************************************************************************/
3261 /** \brief Remove element from the data structure
3262
3263 * @param theGrid - grid to remove from
3264 * @param theElement - element to remove
3265 * @param dispose_connections - also dispose connections (true/false)
3266
3267 This function removes an element from the data structure and inserts it
3268 into the free list. This includes all elementsides, sidevectors and the
3269 elementvector if they exist.
3270
3271 @return <ul>
3272 <li> 0 if ok </li>
3273 <li> 1 no valid object number. </li>
3274 </ul> */
3275 /****************************************************************************/
3276
DisposeElement(GRID * theGrid,ELEMENT * theElement,INT dispose_connections)3277 INT NS_DIM_PREFIX DisposeElement (GRID *theGrid, ELEMENT *theElement, INT dispose_connections)
3278 {
3279 INT i,j,tag;
3280 NODE *theNode;
3281 VERTEX *theVertex;
3282 EDGE *theEdge;
3283 BNDS *bnds;
3284 ELEMENT *theFather;
3285 ELEMENT *succe = SUCCE(theElement);
3286 #ifdef __THREEDIM__
3287 VECTOR *theVector;
3288 DOUBLE *local,fac;
3289 INT k,m,o,l;
3290 #endif
3291
3292 GRID_UNLINK_ELEMENT(theGrid,theElement);
3293
3294 #ifdef __CENTERNODE__
3295 {
3296 theNode = CENTERNODE(theElement);
3297
3298 if (theNode != NULL) SETNFATHER(theNode,NULL);
3299 }
3300 #endif
3301
3302 theFather = EFATHER(theElement);
3303
3304 if (LEVEL(theElement)>0)
3305 {
3306 #ifndef ModelP
3307 ASSERT(theFather != NULL);
3308 #endif
3309
3310 /* check intergrid pointer from father */
3311 if (theFather != NULL)
3312 {
3313 #ifdef ModelP
3314 int index = PRIO2INDEX(EPRIO(theElement));
3315 #else
3316 int index = 0;
3317 #endif
3318 ELEMENT *Next = NULL;
3319 ASSERT(index!=-1 && index<2);
3320
3321 if (SON(theFather,index) == theElement)
3322 {
3323 if (succe != NULL)
3324 {
3325 if (EFATHER(succe)==theFather)
3326 #ifdef ModelP
3327 if (PRIO2INDEX(EPRIO(succe)) == PRIO2INDEX(EPRIO(theElement)))
3328 #endif
3329 {
3330 Next = succe;
3331 }
3332 }
3333 SET_SON(theFather,index,Next);
3334 }
3335
3336 SETNSONS(theFather,NSONS(theFather)-1);
3337
3338 PRINTDEBUG(gm,2,(PFMT "DisposeElement(): elem=" EID_FMTX
3339 " father=" EID_FMTX " son0=%x son1=%x\n",
3340 theGrid->ppifContext().me(),EID_PRTX(theElement),EID_PRTX(theFather),
3341 SON(theFather,0),SON(theFather,1)));
3342 }
3343 }
3344
3345 #ifdef ModelP
3346 /* reset father pointers of sons */
3347 /** \todo possibly some son cannot be reached by GetAllSons, */
3348 /* because their father has not been on this proc and */
3349 /* they lost their father pointers */
3350 if (NSONS(theElement)>0)
3351 {
3352 INT i,j;
3353 ELEMENT *SonList[MAX_SONS];
3354
3355 if (GetAllSons(theElement,SonList)) RETURN(GM_FATAL);
3356
3357 i = 0;
3358 while (SonList[i] != NULL)
3359 {
3360 PRINTDEBUG(gm,2,(PFMT "DisposeElement(): elem=" EID_FMTX
3361 " deleting fatherpointer of son=" EID_FMTX "\n",
3362 theGrid->ppifContext().me(),EID_PRTX(theElement),EID_PRTX(SonList[i])));
3363 SET_EFATHER(SonList[i],NULL);
3364
3365 /* reset VFATHER of centernode vertex */
3366 for (j=0; j<CORNERS_OF_ELEM(SonList[i]); j++)
3367 {
3368 theNode = CORNER(SonList[i],j);
3369 #ifndef __CENTERNODE__
3370 if (CENTERTYPE(theNode) && NFATHER(theNode)!=NULL)
3371 SETNFATHER(theNode,NULL);
3372 #endif
3373 theVertex = MYVERTEX(theNode);
3374 if (VFATHER(theVertex) != NULL && VFATHER(theVertex) == theElement)
3375 VFATHER(theVertex) = NULL;
3376 }
3377 i++;
3378 }
3379 }
3380 #endif
3381
3382 /* remove element sides if it's a boundary element */
3383 if (OBJT(theElement)==BEOBJ)
3384 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
3385 {
3386 bnds = ELEM_BNDS(theElement,i);
3387 if (bnds != NULL)
3388 BNDS_Dispose(MGHEAP(MYMG(theGrid)),bnds);
3389 }
3390
3391 #ifdef __THREEDIM__
3392 /* reset VFATHER of sidenodes */
3393 for (j=0; j<SIDES_OF_ELEM(theElement); j++) {
3394 theNode = GetSideNode(theElement,j);
3395 if (theNode == NULL) continue;
3396 theVertex = MYVERTEX(theNode);
3397 if (VFATHER(MYVERTEX(theNode)) == theElement) {
3398 ELEMENT *theNb = NBELEM(theElement,j);
3399
3400 VFATHER(theVertex) = theNb;
3401 if (theNb != NULL) {
3402 /* calculate new local coords */
3403 k = ONNBSIDE(theVertex);
3404 SETONSIDE(theVertex,k);
3405 m = CORNERS_OF_SIDE(theNb,k);
3406 local = LCVECT(theVertex);
3407 fac = 1.0 / m;
3408 V_DIM_CLEAR(local);
3409 for (o=0; o<m; o++) {
3410 l = CORNER_OF_SIDE(theNb,k,o);
3411 V_DIM_LINCOMB(1.0,local,1.0,
3412 LOCAL_COORD_OF_ELEM(theNb,l),local);
3413 }
3414 V_DIM_SCALE(fac,local);
3415 }
3416 }
3417 SETONNBSIDE(theVertex,MAX_SIDES_OF_ELEM);
3418 }
3419 #endif
3420
3421 for (j=0; j<EDGES_OF_ELEM(theElement); j++)
3422 {
3423 theEdge=GetEdge(CORNER(theElement,CORNER_OF_EDGE(theElement,j,0)),
3424 CORNER(theElement,CORNER_OF_EDGE(theElement,j,1)));
3425 ASSERT(theEdge!=NULL);
3426
3427 if (NO_OF_ELEM(theEdge)<1)
3428 RETURN(GM_ERROR);
3429
3430 /* edit VFATHER of midnode */
3431 if (MIDNODE(theEdge) != NULL)
3432 {
3433 theVertex = MYVERTEX(MIDNODE(theEdge));
3434 if (VFATHER(theVertex) == theElement) {
3435 #ifdef __TWODIM__
3436 theFather = NBELEM(theElement,j);
3437 VFATHER(theVertex) = theFather;
3438 if (theFather != NULL)
3439 {
3440 /* calculate new local coords */
3441 int co0,co1;
3442
3443 /* reconstruct local coordinates of vertex */
3444 co0 = CORNER_OF_EDGE(theFather,j,0);
3445 co1 = CORNER_OF_EDGE(theFather,j,1);
3446
3447 /* local coordinates have to be local towards pe */
3448 V_DIM_LINCOMB(0.5, LOCAL_COORD_OF_ELEM(theFather,co0),
3449 0.5, LOCAL_COORD_OF_ELEM(theFather,co1),
3450 LCVECT(theVertex));
3451 SETONEDGE(theVertex,j);
3452 }
3453 #endif
3454 #ifdef __THREEDIM__
3455 VFATHER(theVertex) = NULL;
3456 #endif
3457 }
3458 }
3459
3460 if (NO_OF_ELEM(theEdge)==1)
3461 DisposeEdge(theGrid,theEdge);
3462 else
3463 DEC_NO_OF_ELEM(theEdge);
3464 }
3465
3466 /* dispose matrices from element-vector */
3467 if (dispose_connections)
3468 if (DisposeConnectionFromElement(theGrid,theElement))
3469 RETURN(1);
3470
3471 for (j=0; j<CORNERS_OF_ELEM(theElement); j++)
3472 {
3473 theNode = CORNER(theElement,j);
3474
3475 #ifdef __OVERLAP2__
3476 if( ce_NO_DELETE_OVERLAP2 != -1 && NO_DELETE_OVERLAP2(theNode) )
3477 {
3478 continue;
3479 }
3480 #endif
3481
3482 if (START(theNode) == NULL)
3483 {
3484 if (NTYPE(theNode)==MID_NODE)
3485 {
3486 if (NFATHER(theNode)!=NULL)
3487 {
3488 MIDNODE((EDGE *)NFATHER(theNode)) = NULL;
3489 }
3490 #ifndef ModelP
3491 /* HEAPFAULT in theFather possible, if in a previous call
3492 some son is not reached by GetAllSons */
3493 else
3494 {
3495 theVertex = MYVERTEX(theNode);
3496 theFather = VFATHER(theVertex);
3497 if (theFather != NULL)
3498 {
3499 INT edge = ONEDGE(theVertex);
3500 theEdge = GetEdge(CORNER(theFather,
3501 CORNER_OF_EDGE(theFather,edge,0)),
3502 CORNER(theFather,
3503 CORNER_OF_EDGE(theFather,edge,1)));
3504 ASSERT(theEdge!=NULL);
3505 MIDNODE(theEdge) = NULL;
3506 }
3507 }
3508 #endif
3509 }
3510 DisposeNode(theGrid,theNode);
3511 }
3512 }
3513
3514 /* reset neighbor pointers referencing element and dispose vectors in sides if */
3515 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
3516 {
3517 ELEMENT *theNeighbor = NBELEM(theElement,i);
3518
3519 #ifdef __THREEDIM__
3520 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,SIDEVEC))
3521 {
3522 theVector = SVECTOR(theElement,i);
3523 if (theVector!=NULL)
3524 {
3525 assert (VCOUNT(theVector) != 0);
3526 assert (VCOUNT(theVector) != 3);
3527 if (VCOUNT(theVector) == 1)
3528 {
3529 if (DisposeVector (theGrid,theVector))
3530 RETURN (1);
3531 }
3532 else
3533 {
3534 if (!FindNeighborElement (theElement,i,&theNeighbor,&j))
3535 RETURN (1);
3536 VOBJECT(theVector) = (GEOM_OBJECT *)theNeighbor;
3537 SETVECTORSIDE(theVector,j);
3538 SETVCOUNT(SVECTOR(theElement,i),1);
3539 }
3540 }
3541 }
3542 #endif
3543 if (theNeighbor!=NULL)
3544 {
3545 for (j=0; j<SIDES_OF_ELEM(theNeighbor); j++)
3546 if (NBELEM(theNeighbor,j)==theElement)
3547 {
3548 SET_NBELEM(theNeighbor,j,NULL);
3549 break;
3550 }
3551 #ifdef ModelP
3552 ASSERT(j<SIDES_OF_ELEM(theNeighbor) || EGHOST(theElement));
3553 #else
3554 ASSERT(j<SIDES_OF_ELEM(theNeighbor));
3555 #endif
3556 }
3557 }
3558
3559 /* dispose vector in center of element */
3560 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,ELEMVEC))
3561 if (DisposeVector (theGrid,EVECTOR(theElement)))
3562 RETURN(1);
3563
3564 #ifdef ModelP
3565 /* free message buffer */
3566 theElement->message_buffer_free();
3567 #endif
3568
3569 /* dispose element */
3570 /* give it a new tag ! (I know this is somewhat ugly) */
3571 tag = TAG(theElement);
3572 if (OBJT(theElement)==BEOBJ)
3573 {
3574 SETOBJT(theElement,MAPPED_BND_OBJT_TAG(tag));
3575 PutFreeObject(theGrid->mg,theElement,
3576 BND_SIZE_TAG(tag),MAPPED_BND_OBJT_TAG(tag));
3577 }
3578 else
3579 {
3580 SETOBJT(theElement,MAPPED_INNER_OBJT_TAG(tag));
3581 PutFreeObject(theGrid->mg,theElement,INNER_SIZE_TAG(tag),
3582 MAPPED_INNER_OBJT_TAG(tag));
3583 }
3584
3585 return(0);
3586 }
3587
3588 #ifndef ModelP
3589 #define DO_NOT_DISPOSE return (2)
3590 #else
3591 #define DO_NOT_DISPOSE dispose=0
3592 #endif
3593
3594 /****************************************************************************/
3595 /** \brief Construct coarse grid from surface
3596
3597 * @param theMG - multigrid to collapse
3598
3599 This function constructs coarse grid from surface. ATTENTION: Use refine $g
3600 to cover always the whole domain with the grid on each level.
3601
3602 @return <ul>
3603 <li> 0 if ok </li>
3604 <li> 1 no valid object number </li>
3605 <li> 2 grid structure not empty or level 0 </li>
3606 </ul> */
3607 /****************************************************************************/
3608
Collapse(MULTIGRID * theMG)3609 INT NS_DIM_PREFIX Collapse (MULTIGRID *theMG)
3610 {
3611 GRID *theGrid;
3612 ELEMENT *theElement;
3613 NODE *theNode;
3614 EDGE *theEdge;
3615 VERTEX *theVertex;
3616 #ifdef ModelP
3617 VECTOR *vec;
3618 #endif
3619 INT tl = TOPLEVEL(theMG);
3620 INT l,i;
3621
3622 if (MG_COARSE_FIXED(theMG))
3623 if (DisposeBottomHeapTmpMemory(theMG))
3624 REP_ERR_RETURN(1);
3625
3626 #ifdef ModelP
3627 DDD_XferBegin(theMG->dddContext());
3628 #ifdef DDDOBJMGR
3629 DDD_ObjMgrBegin();
3630 #endif
3631 #endif
3632 for (l=tl-1; l>=0; l--) {
3633 theGrid = GRID_ON_LEVEL(theMG,l);
3634 for (theNode=PFIRSTNODE(theGrid); theNode != NULL;
3635 theNode = SUCCN(theNode)) {
3636 SONNODE(theNode) = NULL;
3637 SETNFATHER(theNode,NULL);
3638 }
3639 for (theElement=PFIRSTELEMENT(theGrid); theElement != NULL;
3640 theElement = SUCCE(theElement)) {
3641 SETNSONS(theElement,0);
3642 SET_SON(theElement,0,NULL);
3643 #ifdef ModelP
3644 SET_SON(theElement,1,NULL);
3645 #endif
3646 for (i=0; i<EDGES_OF_ELEM(theElement); i++) {
3647 theEdge = GetEdge(CORNER(theElement,
3648 CORNER_OF_EDGE(theElement,i,0)),
3649 CORNER(theElement,
3650 CORNER_OF_EDGE(theElement,i,1)));
3651 MIDNODE(theEdge) = NULL;
3652 }
3653 }
3654 while (PFIRSTELEMENT(theGrid)!=NULL)
3655 if (DisposeElement(theGrid,PFIRSTELEMENT(theGrid),1))
3656 return(1);
3657 while (PFIRSTNODE(theGrid)!=NULL)
3658 {
3659 if (DisposeNode(theGrid,PFIRSTNODE(theGrid)))
3660 return(1);
3661 }
3662 while (PFIRSTVERTEX(theGrid)!=NULL) {
3663 theVertex = PFIRSTVERTEX(theGrid);
3664 GRID_UNLINK_VERTEX(theGrid,theVertex);
3665 GRID_LINK_VERTEX(GRID_ON_LEVEL(theMG,tl),
3666 theVertex,VXPRIO(theVertex));
3667 }
3668 GRID_ON_LEVEL(theMG,l) = NULL;
3669 }
3670
3671 #ifdef ModelP
3672 #ifdef DDDOBJMGR
3673 DDD_ObjMgrEnd();
3674 #endif
3675 DDD_XferEnd(theMG->dddContext());
3676 #endif
3677
3678 /* move top level grid to bottom (level 0) */
3679 theGrid = GRID_ON_LEVEL(theMG,tl);
3680 theGrid->finer = NULL;
3681 theGrid->coarser = NULL;
3682 theGrid->level = 0;
3683 GATTR(theGrid) = GRID_ATTR(theGrid);
3684 GRID_ON_LEVEL(theMG,tl) = NULL;
3685 GRID_ON_LEVEL(theMG,0) = theGrid;
3686 theMG->topLevel = 0;
3687 theMG->fullrefineLevel = 0;
3688 theMG->currentLevel = 0;
3689
3690 for (theNode=PFIRSTNODE(theGrid); theNode != NULL;
3691 theNode = SUCCN(theNode)) {
3692 SETNFATHER(theNode,NULL);
3693 SETNTYPE(theNode,LEVEL_0_NODE);
3694 SETNCLASS(theNode,3);
3695 SETNNCLASS(theNode,0);
3696 SETLEVEL(theNode,0);
3697 VFATHER(MYVERTEX(theNode)) = NULL;
3698 #ifdef ModelP
3699 DDD_AttrSet(PARHDR(theNode),GRID_ATTR(theGrid));
3700 #endif
3701 }
3702 for (theElement=PFIRSTELEMENT(theGrid); theElement != NULL;
3703 theElement = SUCCE(theElement)) {
3704 SETECLASS(theElement,RED_CLASS);
3705 SET_EFATHER(theElement,NULL);
3706 SETLEVEL(theElement,0);
3707 #ifdef ModelP
3708 DDD_AttrSet(PARHDRE(theElement),GRID_ATTR(theGrid));
3709 #endif
3710 for (i=0; i<EDGES_OF_ELEM(theElement); i++) {
3711 theEdge = GetEdge(CORNER(theElement,
3712 CORNER_OF_EDGE(theElement,i,0)),
3713 CORNER(theElement,
3714 CORNER_OF_EDGE(theElement,i,1)));
3715 SETLEVEL(theEdge,0);
3716 #if (defined ModelP) && (defined __THREEDIM__)
3717 DDD_AttrSet(PARHDR(theEdge),GRID_ATTR(theGrid));
3718 #endif
3719 }
3720 }
3721 for (theVertex=PFIRSTVERTEX(theGrid); theVertex != NULL;
3722 theVertex = SUCCV(theVertex)) {
3723 SETLEVEL(theVertex,0);
3724 #ifdef ModelP
3725 DDD_AttrSet(PARHDRV(theVertex),GRID_ATTR(theGrid));
3726 #endif
3727 ASSERT(NOOFNODE(theVertex)==1);
3728 }
3729
3730 #ifdef ModelP
3731 for (vec=PFIRSTVECTOR(theGrid); vec != NULL; vec = SUCCVC(vec))
3732 DDD_AttrSet(PARHDR(vec),GRID_ATTR(theGrid));
3733
3734 /* rebiuld all DDD interfaces due to removed objects and changed attributes */
3735 DDD_IFRefreshAll(theGrid->dddContext());
3736 #endif
3737
3738 if (MG_COARSE_FIXED(theMG))
3739 if (CreateAlgebra(theMG))
3740 REP_ERR_RETURN(1);
3741
3742 return(0);
3743 }
3744
3745 /****************************************************************************/
3746 /** \brief Remove top level grid from multigrid structure
3747
3748 * @param theMG - multigrid to remove from
3749
3750 This function removes the top level grid from multigrid structure.
3751
3752 @return <ul>
3753 <li> 0 if ok </li>
3754 <li> 1 no valid object number </li>
3755 <li> 2 grid structure not empty or level 0 </li>
3756 </ul> */
3757 /****************************************************************************/
3758
DisposeTopLevel(MULTIGRID * theMG)3759 INT NS_DIM_PREFIX DisposeTopLevel (MULTIGRID *theMG)
3760 {
3761 int l;
3762 GRID *theGrid;
3763 #ifdef ModelP
3764 int dispose = 1;
3765 #endif
3766
3767 /* level 0 can not be deleted */
3768 l = theMG->topLevel;
3769 if (l<=0) DO_NOT_DISPOSE;
3770 theGrid = GRID_ON_LEVEL(theMG,l);
3771
3772 /* is level empty */
3773 if (PFIRSTELEMENT(theGrid)!=NULL) DO_NOT_DISPOSE;
3774 if (PFIRSTVERTEX(theGrid)!=NULL) DO_NOT_DISPOSE;
3775 if (PFIRSTNODE(theGrid)!=NULL) DO_NOT_DISPOSE;
3776
3777 #ifdef ModelP
3778 dispose = UG_GlobalMinINT(theMG->ppifContext(), dispose);
3779 if (!dispose) return(2);
3780 #endif
3781
3782 /* remove from grids array */
3783 GRID_ON_LEVEL(theMG,l) = NULL;
3784 GRID_ON_LEVEL(theMG,l-1)->finer = NULL;
3785 (theMG->topLevel)--;
3786 if (theMG->currentLevel>theMG->topLevel)
3787 theMG->currentLevel = theMG->topLevel;
3788
3789 PutFreeObject(theMG,theGrid,sizeof(GRID),GROBJ);
3790
3791 return(0);
3792 }
3793
3794 /****************************************************************************/
3795 /** \brief Dispose top level grid
3796
3797 * @param theGrid - grid to be removed
3798
3799 This function removes the top level grid from multigrid structure.
3800
3801 @return <ul>
3802 <li> 0 if ok </li>
3803 <li> 1 no valid object number </li>
3804 <li> 2 grid structure not empty or level 0 </li>
3805 </ul> */
3806 /****************************************************************************/
3807
DisposeGrid(GRID * theGrid)3808 INT NS_DIM_PREFIX DisposeGrid (GRID *theGrid)
3809 {
3810 MULTIGRID *theMG;
3811
3812 if (theGrid == NULL)
3813 return(0);
3814
3815 theMG = MYMG(theGrid);
3816
3817 if (GLEVEL(theGrid)<0)
3818 return (1);
3819
3820 if (theGrid->finer != NULL)
3821 return(1);
3822
3823 /* clear level */
3824 while (PFIRSTELEMENT(theGrid)!=NULL)
3825 if (DisposeElement(theGrid,PFIRSTELEMENT(theGrid),1))
3826 return(2);
3827
3828 while (PFIRSTNODE(theGrid)!=NULL)
3829 if (DisposeNode(theGrid,PFIRSTNODE(theGrid)))
3830 return(2);
3831
3832 while (PFIRSTVERTEX(theGrid)!=NULL)
3833 if (DisposeVertex(theGrid,PFIRSTVERTEX(theGrid)))
3834 return(4);
3835
3836 /* level 0 can not be deleted */
3837 if (GLEVEL(theGrid) > 0)
3838 return(DisposeTopLevel(theMG));
3839
3840 /* remove from grids array */
3841 GRID_ON_LEVEL(theMG,0) = NULL;
3842 theMG->currentLevel = theMG->topLevel = -1;
3843 theMG->nodeIdCounter = 0;
3844 theMG->vertIdCounter = 0;
3845 theMG->elemIdCounter = 0;
3846
3847 PutFreeObject(theMG,theGrid,sizeof(GRID),GROBJ);
3848
3849 return(0);
3850 }
3851
3852
3853 /****************************************************************************/
3854 /** \brief Release memory for the whole multigrid structure
3855
3856 * @param theMG - multigrid to remove
3857
3858 This function releases the memory for the whole multigrid structure.
3859
3860 @return <ul>
3861 <li> GM_OK if ok </li>
3862 <li> GM_ERROR when error occured. </li>
3863 </ul> */
3864 /****************************************************************************/
3865
DisposeMultiGrid(MULTIGRID * theMG)3866 INT NS_DIM_PREFIX DisposeMultiGrid (MULTIGRID *theMG)
3867 {
3868 INT level;
3869
3870 if (DisposeBottomHeapTmpMemory(theMG)) REP_ERR_RETURN(1);
3871
3872 #ifdef ModelP
3873 /* tell DDD that we will 'inconsistently' delete objects.
3874 this is a dangerous mode as it switches DDD warnings off. */
3875 DDD_SetOption(theMG->dddContext(), OPT_WARNING_DESTRUCT_HDR, OPT_OFF);
3876 #endif
3877
3878 for (level = TOPLEVEL(theMG); level >= 0; level --)
3879 if (DisposeGrid(GRID_ON_LEVEL(theMG,level)))
3880 RETURN(1);
3881
3882 #ifdef ModelP
3883 /* stop dangerous mode. from now on DDD will issue warnings again. */
3884 DDD_SetOption(theMG->dddContext(), OPT_WARNING_DESTRUCT_HDR, OPT_ON);
3885
3886 /* rebuild DDD-interfaces because distributed vectors have been
3887 deleted without communication */
3888 DDD_IFRefreshAll(theMG->dddContext());
3889 #endif
3890
3891 /** \todo Normally the MG-heap should be cleaned-up before freeing.
3892 DDD depends on storage in the heap, even if no DDD objects
3893 are allocated!! (due to free-lists, DDD type definitions
3894 etc.) therefore, repeated new/close commands are inhibited
3895 explicitly in dune/uggrid/parallel/dddif/initddd.c(InitCurrMG()). */
3896 DisposeHeap(MGHEAP(theMG));
3897
3898 /* dispose BVP */
3899 if (MG_BVP(theMG)!=NULL)
3900 if (BVP_Dispose(MG_BVP(theMG))) return (GM_ERROR);
3901
3902 /* first unlock the mg */
3903 ((ENVITEM*) theMG)->v.locked = false;
3904
3905 #ifdef ModelP
3906 ExitDDD(theMG->dddContext());
3907 globalDDDContext(nullptr);
3908 #endif
3909 theMG->~multigrid();
3910
3911 /* delete mg */
3912 if (ChangeEnvDir("/Multigrids")==NULL) RETURN (GM_ERROR);
3913 if (RemoveEnvDir ((ENVITEM *)theMG)) RETURN (GM_ERROR);
3914
3915 return(GM_OK);
3916 }
3917
3918 /****************************************************************************/
3919 /** \brief Determine neighbor and side of neighbor that goes back to element
3920 *
3921 * @param theElement - considered element
3922 * @param Side - side of that element
3923 * @param theNeighbor - handle to neighbor
3924 * @param NeighborSide - number of side of neighbor that goes back to elem
3925 *
3926 This function determines the neighbor and side of the neighbor that goes back to elem.
3927
3928 @return <ul>
3929 <li> 0 if ok </li>
3930 <li> 1 when error occured. </li>
3931 </ul> */
3932 /****************************************************************************/
3933
FindNeighborElement(const ELEMENT * theElement,INT Side,ELEMENT ** theNeighbor,INT * NeighborSide)3934 INT NS_DIM_PREFIX FindNeighborElement (const ELEMENT *theElement, INT Side, ELEMENT **theNeighbor, INT *NeighborSide)
3935 {
3936 INT i;
3937
3938 /* find neighbor */
3939 *theNeighbor = NBELEM(theElement,Side);
3940 if (*theNeighbor == NULL) return (0);
3941
3942 /* search the side */
3943 for (i=0; i<SIDES_OF_ELEM(*theNeighbor); i++)
3944 if (NBELEM(*theNeighbor,i) == theElement)
3945 break;
3946
3947 /* found ? */
3948 if (i<SIDES_OF_ELEM(*theNeighbor))
3949 {
3950 *NeighborSide = i;
3951 return (1);
3952 }
3953 return (0);
3954 }
3955
3956
3957 /****************************************************************************/
3958 /** \brief Insert an inner node
3959 *
3960 * @param theGrid grid structure
3961 * @param pos array containing position
3962 *
3963 * This function inserts a inner node into level 0.
3964 *
3965 * @return <ul>
3966 * <li> pointer to new node if ok </li>
3967 * <li> NULL when error occured </li>
3968 * </ul>
3969 */
3970 /****************************************************************************/
3971
InsertInnerNode(GRID * theGrid,const DOUBLE * pos)3972 NODE * NS_DIM_PREFIX InsertInnerNode (GRID *theGrid, const DOUBLE *pos)
3973 {
3974 VERTEX *theVertex;
3975 NODE *theNode;
3976 INT i;
3977
3978 /* create objects */
3979 theVertex = CreateInnerVertex(theGrid);
3980 if (theVertex==NULL)
3981 {
3982 PrintErrorMessage('E',"InsertInnerNode","cannot create vertex");
3983 return(NULL);
3984 }
3985 theNode = CreateNode(theGrid,theVertex,NULL,LEVEL_0_NODE,0);
3986 if (theNode==NULL)
3987 {
3988 DisposeVertex(theGrid,theVertex);
3989 PrintErrorMessage('E',"InsertInnerNode","cannot create node");
3990 return(NULL);
3991 }
3992
3993 /* fill data */
3994 for (i=0; i<DIM; i++) CVECT(theVertex)[i] = pos[i];
3995 SETMOVE(theVertex,DIM);
3996
3997 return(theNode);
3998 }
3999
4000 /****************************************************************************/
4001 /** \brief Insert a boundary node
4002
4003 * @param theGrid - grid structure
4004 * @param bndp - boundary point descriptor
4005
4006 This function inserts a boundary node into level 0.
4007
4008 @return <ul>
4009 <li> GM_OK if ok </li>
4010 <li> GM_ERROR when error occured. </li>
4011 </ul> */
4012 /****************************************************************************/
4013
InsertBoundaryNode(GRID * theGrid,BNDP * bndp)4014 NODE * NS_DIM_PREFIX InsertBoundaryNode (GRID *theGrid, BNDP *bndp)
4015 {
4016 NODE *theNode;
4017 VERTEX *theVertex;
4018 INT move,part;
4019
4020 /* create objects */
4021 theVertex = CreateBoundaryVertex(theGrid);
4022 if (theVertex==NULL)
4023 {
4024 BNDP_Dispose(MGHEAP(MYMG(theGrid)),bndp);
4025 PrintErrorMessage('E',"InsertBoundaryNode","cannot create vertex");
4026 REP_ERR_RETURN(NULL);
4027 }
4028 if (BNDP_Global(bndp,CVECT(theVertex)))
4029 {
4030 DisposeVertex(theGrid,theVertex);
4031 return(NULL);
4032 }
4033
4034 if (BNDP_BndPDesc(bndp,&move,&part))
4035 {
4036 DisposeVertex(theGrid,theVertex);
4037 return(NULL);
4038 }
4039 SETMOVE(theVertex,move);
4040 V_BNDP(theVertex) = bndp;
4041
4042 theNode = CreateNode(theGrid,theVertex,NULL,LEVEL_0_NODE,0);
4043 if (theNode==NULL)
4044 {
4045 DisposeVertex(theGrid,theVertex);
4046 PrintErrorMessage('E',"InsertBoundaryNode","cannot create node");
4047 REP_ERR_RETURN(NULL);
4048 }
4049 #ifdef TOPNODE
4050 TOPNODE(theVertex) = theNode;
4051 #endif
4052
4053 PRINTDEBUG(dom,1,(" ipn %ld nd %x bndp %x \n",
4054 ID(theNode),theNode,V_BNDP(theVertex)));
4055
4056 SetStringValue(":bndp0",XC(theVertex));
4057 SetStringValue(":bndp1",YC(theVertex));
4058 #ifdef __THREEDIM__
4059 SetStringValue(":bndp2",ZC(theVertex));
4060 #endif
4061
4062 return(theNode);
4063 }
4064
4065 /****************************************************************************/
4066 /** \brief Delete a node
4067
4068 * @param theGrid - grid structure
4069 * @param theNode - node to delete
4070
4071 This function deletes a node from level 0.
4072
4073 @return <ul>
4074 <li> GM_OK if ok </li>
4075 <li> GM_ERROR when error occured. </li>
4076 </ul> */
4077 /****************************************************************************/
4078
DeleteNode(GRID * theGrid,NODE * theNode)4079 INT NS_DIM_PREFIX DeleteNode (GRID *theGrid, NODE *theNode)
4080 {
4081 VERTEX *theVertex;
4082 ELEMENT *theElement;
4083 INT i;
4084
4085 if (theNode==NULL)
4086 {
4087 PrintErrorMessage('E',"DeleteNode","node not found");
4088 RETURN(GM_ERROR);
4089 }
4090
4091 /* check corner */
4092 theVertex = MYVERTEX(theNode);
4093 if (MOVE(theVertex)==0)
4094 {
4095 PrintErrorMessage('E',"DeleteNode","corners cannot be deleted");
4096 RETURN(GM_ERROR);
4097 }
4098
4099 /* check if some element needs that node */
4100 for (theElement=FIRSTELEMENT(theGrid); theElement!=NULL; theElement=SUCCE(theElement))
4101 for (i=0; i<CORNERS_OF_ELEM(theElement); i++)
4102 if (CORNER(theElement,i)==theNode)
4103 {
4104 PrintErrorMessage('E',"DeleteNode","there is an element needing that node");
4105 RETURN(GM_ERROR);
4106 }
4107
4108 /* now allowed to delete */
4109 DisposeNode(theGrid,theNode);
4110
4111 return(GM_OK);
4112 }
4113
4114 #ifdef __TWODIM__
4115
4116
4117 /****************************************************************************/
4118 /** \todo Please doc me!
4119
4120 CheckOrientation -
4121
4122 SYNOPSIS:
4123 INT CheckOrientation (INT n, VERTEX **vertices);
4124
4125
4126 * @param n
4127 * @param vertices
4128
4129 DESCRIPTION:
4130
4131 @return
4132 INT
4133 */
4134 /****************************************************************************/
4135
CheckOrientation(INT n,VERTEX ** vertices)4136 INT NS_DIM_PREFIX CheckOrientation (INT n, VERTEX **vertices)
4137 {
4138 int i;
4139 DOUBLE x1,x2,y1,y2;
4140
4141 for (i=0; i<n; i++)
4142 {
4143 x1 = XC(vertices[(i+1)%n])-XC(vertices[i]);
4144 x2 = XC(vertices[(i+n-1)%n])-XC(vertices[i]);
4145 y1 = YC(vertices[(i+1)%n])-YC(vertices[i]);
4146 y2 = YC(vertices[(i+n-1)%n])-YC(vertices[i]);
4147 if (vp(x1,y1,x2,y2)<SMALL_C)
4148 {
4149 return(0);
4150 }
4151 }
4152 return(1);
4153 }
4154
4155 #define SWAP_IJ(a,i,j,t) {t = a[i]; a[i] = a[j]; a[j] = t;}
4156 #endif
4157
4158 #ifdef __THREEDIM__
4159
4160
4161 /****************************************************************************/
4162 /** \todo Please doc me!
4163 CheckOrientation -
4164
4165 SYNOPSIS:
4166 INT CheckOrientation (INT n, VERTEX **vertices);
4167
4168
4169 * @param n
4170 * @param vertices
4171
4172 DESCRIPTION:
4173
4174 @return
4175 INT
4176 */
4177 /****************************************************************************/
4178
CheckOrientation(INT n,VERTEX ** vertices)4179 INT NS_DIM_PREFIX CheckOrientation (INT n, VERTEX **vertices)
4180 {
4181 DOUBLE_VECTOR diff[3],rot;
4182 DOUBLE det;
4183 INT i;
4184
4185 /* TODO: this case */
4186 if (n == 8 || n==6 || n==5)
4187 return(1);
4188
4189 for (i=1; i<n; i++)
4190 V3_SUBTRACT(CVECT(vertices[i]),CVECT(vertices[0]),diff[i-1]);
4191 V3_VECTOR_PRODUCT(diff[0],diff[1],rot);
4192 V3_SCALAR_PRODUCT(rot,diff[2],det);
4193
4194 if (det < 0.0)
4195 return(0);
4196
4197 return(1);
4198 }
4199 #endif
4200
4201
4202 /****************************************************************************/
4203 /** \todo Please doc me!
4204 CheckOrientationInGrid -
4205
4206 SYNOPSIS:
4207 INT CheckOrientationInGrid (GRID *theGrid);
4208
4209
4210 * @param theGrid
4211
4212 DESCRIPTION:
4213
4214 @return
4215 INT
4216 */
4217 /****************************************************************************/
4218
CheckOrientationInGrid(GRID * theGrid)4219 INT NS_DIM_PREFIX CheckOrientationInGrid (GRID *theGrid)
4220 {
4221 ELEMENT *theElement;
4222 NODE *theNode;
4223 VERTEX *vertices[MAX_CORNERS_OF_ELEM];
4224 INT i;
4225
4226 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL; theElement=SUCCE(theElement))
4227 {
4228 for (i=0; i<CORNERS_OF_ELEM(theElement); i++)
4229 {
4230 theNode = CORNER(theElement,i);
4231 if (theNode==NULL) return (1);
4232 vertices[i] = MYVERTEX(theNode);
4233 if (vertices[i]==NULL) return (1);
4234 }
4235 if (!CheckOrientation (CORNERS_OF_ELEM(theElement),vertices)) return (1);
4236 }
4237
4238 return (0);
4239 }
4240
4241
4242 /****************************************************************************/
4243 /** \todo Please doc me!
4244
4245 * @param[in] n Number of element corners
4246 * @param[in] Node
4247 * @param[in] theMG
4248 * @param[out] NbrS
4249 * @param[out] Nbr
4250
4251 @return 0 if all went well, 1 if an error occurred
4252 */
4253 /****************************************************************************/
4254
NeighborSearch_O_n(INT n,ELEMENT * theElement,NODE ** Node,MULTIGRID * theMG,INT * NbrS,ELEMENT ** Nbr)4255 static INT NeighborSearch_O_n(INT n, ELEMENT *theElement, NODE **Node, MULTIGRID *theMG, INT *NbrS, ELEMENT **Nbr)
4256 {
4257
4258 #if 0
4259 // this is (as a reference) the O(n^2) version
4260 INT i,jj,k,m,num;
4261 NODE *sideNode[MAX_CORNERS_OF_SIDE];
4262 ELEMENT *theOther;
4263 NODE *NeighborNode;
4264 MULTIGRID *theMG = MYMG(theGrid);
4265
4266 /*O(n*n)InsertElement ...*/
4267 /* for all sides of the element to be created */
4268 for (i=0; i<SIDES_OF_REF(n); i++)
4269 {
4270 for(jj=0; jj<CORNERS_OF_SIDE_REF(n,i); jj++ )
4271 sideNode[jj] = Node[CORNER_OF_SIDE_REF(n,i,jj)];
4272
4273 /* for all neighbouring elements already inserted */
4274 for (theOther=FIRSTELEMENT(theGrid); theOther!=NULL;
4275 theOther=SUCCE(theOther))
4276 {
4277 if (theOther == theElement)
4278 continue;
4279 /* for all sides of the neighbour element */
4280 for (jj=0; jj<SIDES_OF_ELEM(theOther); jj++)
4281 {
4282 num = 0;
4283 /* for all corners of the side of the neighbour */
4284 for (m=0; m<CORNERS_OF_SIDE(theOther,jj); m++)
4285 {
4286 NeighborNode = CORNER(theOther,
4287 CORNER_OF_SIDE(theOther,jj,m));
4288 /* for all corners of the side of the
4289 element to be created */
4290 for (k=0; k<CORNERS_OF_SIDE_REF(n,i); k++)
4291 if(NeighborNode==sideNode[k])
4292 {
4293 num++;
4294 break;
4295 }
4296 }
4297 if(num==CORNERS_OF_SIDE_REF(n,i))
4298 {
4299 if (NBELEM(theOther,jj)!=NULL)
4300 {
4301 PrintErrorMessage('E',"InsertElement -> NeighborSearch_O_nn",
4302 "neighbor relation inconsistent");
4303 return(1);
4304 }
4305 Nbr[i] = theOther;
4306 NbrS[i] = jj;
4307 }
4308 }
4309 }
4310 }
4311 /* ... O(n*n)InsertElement */
4312 #else
4313 /*O(n*n)InsertElement ...*/
4314 /* for all sides of the element to be created */
4315 MULTIGRID::FaceNodes faceNodes;
4316 for (int i=0; i<SIDES_OF_REF(n); i++)
4317 {
4318 int j = 0;
4319 for(j=0; j<CORNERS_OF_SIDE_REF(n,i); j++)
4320 faceNodes[j] = Node[CORNER_OF_SIDE_REF(n,i,j)];
4321 for(; j<MAX_CORNERS_OF_SIDE; j++)
4322 faceNodes[j] = 0;
4323 std::sort(faceNodes.begin(), faceNodes.begin()+CORNERS_OF_SIDE_REF(n,i));
4324
4325 // try to write myself...
4326 auto result = theMG->facemap.emplace(faceNodes,std::make_pair(theElement,i));
4327 // if this failed (i.e. result.second == false) an entry already exists
4328 if (! result.second)
4329 {
4330 // update neighbor list the other entry, stored in result.first->second
4331 auto & data = result.first->second;
4332 ELEMENT* theOther = data.first;
4333 int idx = data.second;
4334 Nbr[i] = theOther;
4335 NbrS[i] = idx;
4336 theMG->facemap.erase(faceNodes);
4337 }
4338
4339 }
4340
4341 #endif
4342
4343 return(0);
4344 /*... O(n)InsertElement ...*/
4345 } /*of static INT NeighborSearch_O_n()*/
4346
4347
4348
4349 /****************************************************************************/
4350 /** \todo Please doc me!
4351 Neighbor_Direct_Insert -
4352
4353 SYNOPSIS:
4354 static INT Neighbor_Direct_Insert(INT n, ELEMENT **ElemList, INT *NbgSdList, INT* NbrS, ELEMENT **Nbr);
4355
4356
4357 * @param n
4358 * @param ElemList
4359 * @param NbgSdList
4360 * @param NbrS
4361 * @param Nbr
4362
4363 DESCRIPTION:
4364
4365 @return
4366 INT
4367 */
4368 /****************************************************************************/
4369
Neighbor_Direct_Insert(INT n,ELEMENT ** ElemList,INT * NbgSdList,INT * NbrS,ELEMENT ** Nbr)4370 static INT Neighbor_Direct_Insert(INT n, ELEMENT **ElemList, INT *NbgSdList, INT* NbrS, ELEMENT **Nbr)
4371 {
4372 INT i;
4373
4374 for (i=0; i<SIDES_OF_REF(n); i++)
4375 Nbr[i] = ElemList[i];
4376 if (NbgSdList != NULL)
4377 for (i=0; i<SIDES_OF_REF(n); i++)
4378 NbrS[i] = NbgSdList[i];
4379
4380 return(0);
4381 }
4382
4383
4384 /****************************************************************************/
4385 /** \brief Insert an element
4386
4387 * @param theGrid - grid structure
4388 * @param[in] n Number of vertices of the element to be inserted
4389 * @param Node
4390 * @param ElemList
4391 * @param NbgSdList
4392 * @param bnds_flag
4393
4394 This function inserts an element
4395
4396 \return Pointer to the newly created element, NULL if an error occurred
4397
4398 */
4399 /****************************************************************************/
4400
InsertElement(GRID * theGrid,INT n,NODE ** Node,ELEMENT ** ElemList,INT * NbgSdList,INT * bnds_flag)4401 ELEMENT * NS_DIM_PREFIX InsertElement (GRID *theGrid, INT n, NODE **Node, ELEMENT **ElemList, INT *NbgSdList, INT *bnds_flag)
4402 {
4403 MULTIGRID *theMG;
4404 INT i,j,k,m,rv,tag,ElementType;
4405 INT NeighborSide[MAX_SIDES_OF_ELEM];
4406 [[maybe_unused]] NODE *sideNode[MAX_CORNERS_OF_SIDE];
4407 VERTEX *Vertex[MAX_CORNERS_OF_ELEM],*sideVertex[MAX_CORNERS_OF_SIDE];
4408 ELEMENT *theElement,*Neighbor[MAX_SIDES_OF_ELEM];
4409 BNDS *bnds[MAX_SIDES_OF_ELEM];
4410 BNDP *bndp[MAX_CORNERS_OF_ELEM];
4411 #ifdef __TWODIM__
4412 VERTEX *theVertex;
4413 NODE *theNode;
4414 #endif
4415
4416 theMG = MYMG(theGrid);
4417
4418 // nodes are already inserted, so we know how many there are...
4419 if (theMG->facemap.bucket_count() <= 1)
4420 {
4421 // try to allocate the right size a-priori to avoid rehashing
4422 theMG->facemap.rehash(theMG->nodeIdCounter);
4423 // theMG->facemap.max_load_factor(1000);
4424 }
4425
4426 /* check parameters */
4427 #ifdef __TWODIM__
4428 switch (n)
4429 {
4430 case 3 :
4431 tag = TRIANGLE;
4432 break;
4433 case 4 :
4434 tag = QUADRILATERAL;
4435 break;
4436 default :
4437 PrintErrorMessage('E',"InsertElement","only triangles and quadrilaterals allowed in 2D");
4438 return(NULL);
4439 }
4440 #endif
4441
4442 #ifdef __THREEDIM__
4443 switch (n)
4444 {
4445 case 4 :
4446 tag = TETRAHEDRON;
4447 break;
4448 case 5 :
4449 tag = PYRAMID;
4450 break;
4451 case 6 :
4452 tag = PRISM;
4453 break;
4454 case 8 :
4455 tag = HEXAHEDRON;
4456 break;
4457 default :
4458 PrintErrorMessage('E',"InsertElement","only tetrahedra, prisms, pyramids, and hexahedra are allowed in the 3D coarse grid");
4459 return(NULL);
4460 }
4461 #endif
4462
4463 /* init vertices */
4464 for (i=0; i<n; i++)
4465 {
4466 PRINTDEBUG(gm,1,("InsertElement(): node[%d]=" ID_FMTX "vertex[%d]=" VID_FMTX "\n",
4467 i,ID_PRTX(Node[i]),i,VID_PRTX(MYVERTEX(Node[i]))))
4468 Vertex[i] = MYVERTEX(Node[i]);
4469 }
4470
4471 #ifdef __TWODIM__
4472 /* find orientation */
4473 if (!CheckOrientation(n,Vertex))
4474 {
4475 /* flip order */
4476 SWAP_IJ(Node, 0,n/2,theNode);
4477 SWAP_IJ(Vertex,0,n/2,theVertex);
4478
4479 if (!CheckOrientation(n,Vertex))
4480 {
4481 /* this was the only possibility for a triangle: so is a nonconvex quadrilateral */
4482 /* interchange first two nodes and try again */
4483 SWAP_IJ(Node, 0,1,theNode);
4484 SWAP_IJ(Vertex,0,1,theVertex);
4485 if (!CheckOrientation(n,Vertex))
4486 {
4487 /* flip order */
4488 SWAP_IJ(Node, 0,n/2,theNode);
4489 SWAP_IJ(Vertex,0,n/2,theVertex);
4490 if (!CheckOrientation(n,Vertex))
4491 {
4492 /* flip order back */
4493 SWAP_IJ(Node, 0,n/2,theNode);
4494 SWAP_IJ(Vertex,0,n/2,theVertex);
4495 /* interchange second two nodes and try again */
4496 SWAP_IJ(Node, 1,2,theNode);
4497 SWAP_IJ(Vertex,1,2,theVertex);
4498 if (!CheckOrientation(n,Vertex))
4499 {
4500 /* flip order */
4501 SWAP_IJ(Node, 0,n/2,theNode);
4502 SWAP_IJ(Vertex,0,n/2,theVertex);
4503 if (!CheckOrientation(n,Vertex))
4504 {
4505 PrintErrorMessage('E',"InsertElement",
4506 "cannot find orientation");
4507 return(NULL);
4508 }
4509 }
4510 }
4511 }
4512 }
4513 }
4514 #endif
4515
4516 #ifdef __THREEDIM__
4517 if (!CheckOrientation (n,Vertex))
4518 {
4519 sideNode[0] = Node[0];
4520 sideVertex[0] = Vertex[0];
4521 Node[0] = Node[1];
4522 Vertex[0] = Vertex[1];
4523 Node[1] = sideNode[0];
4524 Vertex[1] = sideVertex[0];
4525 }
4526 #endif
4527
4528 /* init pointers */
4529 for (i=0; i<SIDES_OF_REF(n); i++)
4530 {
4531 Neighbor[i] = NULL;
4532 bnds[i] = NULL;
4533 }
4534
4535 /* compute side information (theSeg[i]==NULL) means inner side */
4536 ElementType = IEOBJ;
4537 for (i=0; i<SIDES_OF_REF(n); i++)
4538 {
4539 m = CORNERS_OF_SIDE_REF(n,i);
4540 for(j=0; j<m; j++ )
4541 {
4542 k = CORNER_OF_SIDE_REF(n,i,j);
4543 sideNode[j] = Node[k];
4544 sideVertex[j] = Vertex[k];
4545 }
4546 bool found = false;
4547 for(j=0; j<m; j++ )
4548 {
4549 if( OBJT(sideVertex[j]) == IVOBJ ) found = true;
4550 }
4551 if( found ) continue;
4552
4553 /* all vertices of side[i] are on the boundary now */
4554
4555 /* We now assume, that: */
4556 /* if bnds_flag!=NULL && bnds_flag[i]!=0 there has to be a bnds */
4557 /* so, if not -->error */
4558 /* or: if bnds_flag==NULL, the domain decides whether there */
4559 /* should be a bnds or not (never an error) */
4560
4561 for (j=0; j<m; j++)
4562 bndp[j] = V_BNDP(sideVertex[j]);
4563
4564 if (bnds_flag==NULL)
4565 {
4566 bnds[i] = BNDP_CreateBndS(MGHEAP(theMG),bndp,m);
4567 if (bnds[i] != NULL) ElementType = BEOBJ;
4568 }
4569 else if (bnds_flag[i]!=0)
4570 {
4571 bnds[i] = BNDP_CreateBndS(MGHEAP(theMG),bndp,m);
4572 assert(bnds[i]!=NULL);
4573 ElementType = BEOBJ;
4574 }
4575 }
4576
4577 /* create the new Element */
4578 theElement = CreateElement(theGrid,tag,ElementType,Node,NULL,0);
4579 if (theElement==NULL)
4580 {
4581 PrintErrorMessage('E',"InsertElement","cannot allocate element");
4582 return(NULL);
4583 }
4584
4585 if (ElemList == NULL)
4586 {
4587 /* using the fast O(n) algorithm */
4588 NeighborSearch_O_n(n, theElement, Node, theMG, NeighborSide, Neighbor);
4589 }
4590 else
4591 {
4592 /* use given neighboring elements */
4593 if ( (rv = Neighbor_Direct_Insert(n, ElemList, NbgSdList, NeighborSide, Neighbor)) == 1)
4594 {
4595 PrintErrorMessage('E',"InsertElement"," ERROR by calling Neighbor_Direct_Insert()");
4596 return(NULL);
4597 }
4598 }
4599
4600 /* create element sides if necessary */
4601 if (OBJT(theElement)==BEOBJ)
4602 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
4603 SET_BNDS(theElement,i,bnds[i]);
4604
4605 /* fill element data */
4606 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
4607 {
4608 SET_NBELEM(theElement,i,Neighbor[i]);
4609 if (Neighbor[i]!=NULL)
4610 {
4611 if (NbgSdList == NULL)
4612 NeighborSide[i] = SideOfNbElement(theElement,i);
4613 if (NeighborSide[i] >= MAX_SIDES_OF_ELEM) {
4614 PrintErrorMessage('E',"InsertElement",
4615 "neighbor relation inconsistent");
4616 return(NULL);
4617 }
4618 SET_NBELEM(Neighbor[i],NeighborSide[i],theElement);
4619 #ifdef __THREEDIM__
4620 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,SIDEVEC))
4621 if (DisposeDoubledSideVector(theGrid,Neighbor[i],
4622 NeighborSide[i],theElement,i))
4623 return(NULL);
4624 #endif
4625 }
4626 }
4627
4628 SET_EFATHER(theElement,NULL);
4629 SETECLASS(theElement,RED_CLASS);
4630
4631 return(theElement);
4632 }
4633
4634 /****************************************************************************/
4635 /** \brief Delete an element
4636
4637 * @param theMG - multigrid structure
4638 * @param theElement - element to delete
4639
4640 This function deletes an element from level 0.
4641
4642 @return <ul>
4643 <li> GM_OK if ok </li>
4644 <li> GM_ERROR when error occured. </li>
4645 </ul> */
4646 /****************************************************************************/
4647
DeleteElement(MULTIGRID * theMG,ELEMENT * theElement)4648 INT NS_DIM_PREFIX DeleteElement (MULTIGRID *theMG, ELEMENT *theElement) /* 3D VERSION */
4649 {
4650 GRID *theGrid;
4651 ELEMENT *theNeighbor;
4652 INT i,j,found;
4653
4654 /* check level */
4655 if ((CURRENTLEVEL(theMG)!=0)||(TOPLEVEL(theMG)!=0))
4656 {
4657 PrintErrorMessage('E',"DeleteElement",
4658 "only a multigrid with exactly one level can be edited");
4659 RETURN(GM_ERROR);
4660 }
4661 theGrid = GRID_ON_LEVEL(theMG,0);
4662
4663 /* delete pointers in neighbors */
4664 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
4665 {
4666 theNeighbor = NBELEM(theElement,i);
4667 if (theNeighbor!=NULL)
4668 {
4669 found = 0;
4670 for (j=0; j<SIDES_OF_ELEM(theNeighbor); j++)
4671 if (NBELEM(theNeighbor,j)==theElement)
4672 {
4673 found++;
4674 SET_NBELEM(theNeighbor,j,NULL);
4675 }
4676 if (found!=1) RETURN(GM_ERROR);
4677 }
4678 }
4679
4680 /* delete element now */
4681 DisposeElement(theGrid,theElement,true);
4682
4683 return(GM_OK);
4684 }
4685
4686
4687 /****************************************************************************/
4688 /** \brief Insert a mesh described by the domain
4689
4690 * @param theMG - multigrid structure
4691 * @param theMesh - mesh structure
4692
4693 This function inserts all nodes and elements given by the mesh.
4694
4695 @return <ul>
4696 <li> GM_OK if ok </li>
4697 <li> GM_ERROR when error occured. </li>
4698 </ul> */
4699 /****************************************************************************/
4700
InsertMesh(MULTIGRID * theMG,MESH * theMesh)4701 INT NS_DIM_PREFIX InsertMesh (MULTIGRID *theMG, MESH *theMesh)
4702 {
4703 GRID *theGrid;
4704 ELEMENT *theElement;
4705 NODE **NList,*Nodes[MAX_CORNERS_OF_ELEM],*ListNode;
4706 VERTEX **VList;
4707 INT i,k,n,nv,j,maxlevel,l,move,part;
4708 INT ElemSideOnBnd[MAX_SIDES_OF_ELEM];
4709 INT MarkKey = MG_MARK_KEY(theMG);
4710
4711 if (theMesh == NULL) return(GM_OK);
4712 if (theMesh->nElements == NULL)
4713 {
4714 assert(theMesh->VertexLevel==NULL);
4715 theGrid = GRID_ON_LEVEL(theMG,0);
4716 for (i=0; i<theMesh->nBndP; i++)
4717 if (InsertBoundaryNode(theGrid,theMesh->theBndPs[i]) == NULL)
4718 REP_ERR_RETURN(GM_ERROR);
4719
4720 for (i=0; i<theMesh->nInnP; i++)
4721 if (InsertInnerNode(theGrid,theMesh->Position[i]) == NULL)
4722 REP_ERR_RETURN(GM_ERROR);
4723 return(GM_OK);
4724 }
4725
4726 /* prepare */
4727 nv = theMesh->nBndP + theMesh->nInnP;
4728 VList = (VERTEX **) GetTmpMem(MGHEAP(theMG),nv*sizeof(VERTEX *),MarkKey);
4729 if (VList == NULL) return(GM_ERROR);
4730 NList = (NODE **) GetTmpMem(MGHEAP(theMG),nv*sizeof(NODE *),MarkKey);
4731 if (NList == NULL) return(GM_ERROR);
4732 for (j=0; j<nv; j++) NList[j] = NULL;
4733
4734 maxlevel = 0;
4735 if (theMesh->VertexLevel!=NULL)
4736 {
4737 for (i=0; i<theMesh->nBndP; i++)
4738 {
4739 theGrid = GRID_ON_LEVEL(theMG,theMesh->VertexLevel[i]);
4740 VList[i] = CreateBoundaryVertex(theGrid);
4741 assert(VList[i]!=NULL);
4742 if (BNDP_Global(theMesh->theBndPs[i],CVECT(VList[i]))) assert(0);
4743 if (BNDP_BndPDesc(theMesh->theBndPs[i],&move,&part))
4744 return(GM_OK);
4745 SETMOVE(VList[i],move);
4746 V_BNDP(VList[i]) = theMesh->theBndPs[i];
4747 maxlevel = MAX(maxlevel,theMesh->VertexLevel[i]);
4748 }
4749 for (i=theMesh->nBndP; i<nv; i++)
4750 {
4751 theGrid = GRID_ON_LEVEL(theMG,theMesh->VertexLevel[i]);
4752 VList[i] = CreateInnerVertex(theGrid);
4753 V_DIM_COPY(theMesh->Position[i-theMesh->nBndP],CVECT(VList[i]));
4754 maxlevel = MAX(maxlevel,theMesh->VertexLevel[i]);
4755 }
4756 }
4757 else
4758 {
4759 theGrid = GRID_ON_LEVEL(theMG,0);
4760 for (i=0; i<theMesh->nBndP; i++)
4761 {
4762 VList[i] = CreateBoundaryVertex(theGrid);
4763 assert(VList[i]!=NULL);
4764 if (BNDP_Global(theMesh->theBndPs[i],CVECT(VList[i]))) assert(0);
4765 if (BNDP_BndPDesc(theMesh->theBndPs[i],&move,&part))
4766 return(GM_OK);
4767 SETMOVE(VList[i],move);
4768 V_BNDP(VList[i]) = theMesh->theBndPs[i];
4769 }
4770 for (i=theMesh->nBndP; i<nv; i++)
4771 {
4772 VList[i] = CreateInnerVertex(theGrid);
4773 V_DIM_COPY(theMesh->Position[i-theMesh->nBndP],CVECT(VList[i]));
4774 }
4775 }
4776 if (theMesh->nElements == NULL)
4777 return(GM_OK);
4778 for (j=1; j<=theMesh->nSubDomains; j++)
4779 for (k=0; k<theMesh->nElements[j]; k++)
4780 {
4781 if (theMesh->ElementLevel!=NULL) i = theMesh->ElementLevel[j][k];
4782 else i=0;
4783 theGrid = GRID_ON_LEVEL(theMG,i);
4784 n = theMesh->Element_corners[j][k];
4785 for (l=0; l<n; l++)
4786 {
4787 ListNode = NList[theMesh->Element_corner_ids[j][k][l]];
4788 if (ListNode==NULL || LEVEL(ListNode)<i)
4789 {
4790 Nodes[l] = CreateNode(theGrid,VList[theMesh->Element_corner_ids[j][k][l]],NULL,LEVEL_0_NODE,0);
4791 if (Nodes[l]==NULL) assert(0);
4792 NList[theMesh->Element_corner_ids[j][k][l]] = Nodes[l];
4793 if (ListNode==NULL || LEVEL(ListNode)<i-1)
4794 {
4795 SETNFATHER(Nodes[l],NULL);
4796 }
4797 else
4798 {
4799 SETNFATHER(Nodes[l],(GEOM_OBJECT *)ListNode);
4800 SONNODE(ListNode) = Nodes[l];
4801 }
4802 }
4803 else
4804 {
4805 Nodes[l] = ListNode;
4806 }
4807 }
4808 if (theMesh->ElemSideOnBnd==NULL)
4809 theElement = InsertElement (theGrid,n,Nodes,NULL,NULL,NULL);
4810 else
4811 {
4812 for (l=0; l<SIDES_OF_REF(n); l++) ElemSideOnBnd[l] = (theMesh->ElemSideOnBnd[j][k]&(1<<l));
4813 theElement = InsertElement (theGrid,n,Nodes,NULL,NULL,ElemSideOnBnd);
4814 }
4815 SETSUBDOMAIN(theElement,j);
4816 }
4817
4818 return(GM_OK);
4819 }
4820
4821 /****************************************************************************/
4822 /** \brief Determine whether point is contained in element
4823
4824 * @param x - coordinates of given point
4825 * @param theElement - element to scan
4826
4827 This function determines whether a given point specified by coordinates `x`
4828 is contained in an element.
4829
4830 @return <ul>
4831 <li> false: point is not in the element</li>
4832 <li> true: point is in the element</li>
4833 </ul> */
4834 /****************************************************************************/
4835
4836 #ifdef __TWODIM__
PointInElement(const DOUBLE * x,const ELEMENT * theElement)4837 bool NS_DIM_PREFIX PointInElement (const DOUBLE *x, const ELEMENT *theElement) /* 2D version */
4838 {
4839 COORD_POINT point[MAX_CORNERS_OF_ELEM],thePoint;
4840 int n,i;
4841
4842 /* check element */
4843 if (theElement==NULL)
4844 return false;
4845
4846 /* load geometrical data of the corners */
4847 n = CORNERS_OF_ELEM(theElement);
4848 for (i=0; i<n; i++)
4849 {
4850 point[i].x = XC(MYVERTEX(CORNER(theElement,i)));
4851 point[i].y = YC(MYVERTEX(CORNER(theElement,i)));
4852 }
4853 thePoint.x = (DOUBLE)x[0];
4854 thePoint.y = (DOUBLE)x[1];
4855
4856 return PointInPolygon(point,n,thePoint);
4857 }
4858 #endif
4859
4860 #ifdef __THREEDIM__
PointInElement(const DOUBLE * global,const ELEMENT * theElement)4861 bool NS_DIM_PREFIX PointInElement (const DOUBLE *global, const ELEMENT *theElement)
4862 {
4863 DOUBLE *x[MAX_CORNERS_OF_ELEM];
4864 DOUBLE_VECTOR a,b,rot;
4865 DOUBLE det;
4866 INT n,i;
4867
4868 /* check element */
4869 if (theElement==NULL)
4870 return false;
4871
4872 CORNER_COORDINATES(theElement,n,x);
4873
4874 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
4875 {
4876 V3_SUBTRACT(x[CORNER_OF_SIDE(theElement,i,1)],
4877 x[CORNER_OF_SIDE(theElement,i,0)],a);
4878 V3_SUBTRACT(x[CORNER_OF_SIDE(theElement,i,2)],
4879 x[CORNER_OF_SIDE(theElement,i,0)],b);
4880 V3_VECTOR_PRODUCT(a,b,rot);
4881 V3_SUBTRACT(global,x[CORNER_OF_SIDE(theElement,i,0)],b);
4882 V3_SCALAR_PRODUCT(rot,b,det);
4883 if (det > SMALL_C)
4884 return false;
4885 }
4886
4887 return true;
4888 }
4889
4890 #endif
4891
4892
4893 /****************************************************************************/
4894 /** \brief Determine whether point is on an element side
4895
4896 * @param x - coordinates of given point
4897 * @param theElement - element to scan
4898 * @param side - the element side
4899
4900 This function determines whether a given point specified by coordinates `x`
4901 is contained in an element side.
4902
4903 Beware: The function only tests if the Point is in the plane spawned by the element side.
4904 The point could be outside the element side area.
4905
4906 @return <ul>
4907 <li> 0 not on side </li>
4908 <li> 1 x is on side </li>
4909 </ul> */
4910 /****************************************************************************/
4911
4912 #ifdef __TWODIM__
PointOnSide(const DOUBLE * global,const ELEMENT * theElement,INT side)4913 INT NS_DIM_PREFIX PointOnSide(const DOUBLE *global, const ELEMENT *theElement, INT side)
4914 {
4915 INT n;
4916 DOUBLE *x[MAX_CORNERS_OF_ELEM];
4917 DOUBLE M[DIM+DIM];
4918 DOUBLE *a, *b;
4919 DOUBLE det;
4920
4921 a = &M[0];
4922 b = &M[DIM];
4923
4924 CORNER_COORDINATES(theElement,n,x);
4925
4926 V2_SUBTRACT(x[CORNER_OF_SIDE(theElement,side,1)], x[CORNER_OF_SIDE(theElement,side,0)], a);
4927 V2_SUBTRACT(global, x[CORNER_OF_SIDE(theElement,side,0)], b);
4928 det = M2_DET(M);
4929 if (fabs(det) < SMALL_C)
4930 return 1;
4931
4932 return 0;
4933 }
4934 #else
PointOnSide(const DOUBLE * global,const ELEMENT * theElement,INT side)4935 INT NS_DIM_PREFIX PointOnSide(const DOUBLE *global, const ELEMENT *theElement, INT side)
4936 {
4937 INT n;
4938 DOUBLE *x[MAX_CORNERS_OF_ELEM];
4939 DOUBLE M[DIM*DIM];
4940 DOUBLE *a, *b, *c;
4941 DOUBLE det;
4942
4943 a = &M[0];
4944 b = &M[DIM];
4945 c = &M[2*DIM];
4946
4947 CORNER_COORDINATES(theElement,n,x);
4948
4949 V3_SUBTRACT(x[CORNER_OF_SIDE(theElement,side,1)], x[CORNER_OF_SIDE(theElement,side,0)], a);
4950 V3_SUBTRACT(x[CORNER_OF_SIDE(theElement,side,2)], x[CORNER_OF_SIDE(theElement,side,0)], b);
4951 V3_SUBTRACT(global, x[CORNER_OF_SIDE(theElement,side,0)], c);
4952 det = M3_DET(M);
4953 if (fabs(det) < SMALL_C)
4954 return 1;
4955
4956 return 0;
4957 }
4958 #endif
4959
4960 /****************************************************************************/
4961 /** \brief Determine distance of a point to an element side
4962
4963 * @param x - coordinates of given point
4964 * @param theElement - element to scan
4965 * @param side - the element side
4966
4967 This function determines the distance of a given point specified by coordinates `x`
4968 from an element side.
4969
4970 Beware: The function only tests if the Point is in the plane spawned by the element side.
4971 The point could be outside the element side area.
4972
4973 @return <ul>
4974 <li> 0 not on side </li>
4975 <li> 1 x is on side </li>
4976 </ul> */
4977 /****************************************************************************/
4978
4979 #ifdef __TWODIM__
DistanceFromSide(const DOUBLE * global,const ELEMENT * theElement,INT side)4980 DOUBLE NS_DIM_PREFIX DistanceFromSide(const DOUBLE *global, const ELEMENT *theElement, INT side)
4981 {
4982 INT n;
4983 DOUBLE *x[MAX_CORNERS_OF_ELEM];
4984 DOUBLE M[DIM+DIM];
4985 DOUBLE *a, *b;
4986 DOUBLE det;
4987
4988 a = &M[0];
4989 b = &M[DIM];
4990
4991 CORNER_COORDINATES(theElement,n,x);
4992
4993 V2_SUBTRACT(x[CORNER_OF_SIDE(theElement,side,1)], x[CORNER_OF_SIDE(theElement,side,0)], a);
4994 V2_SUBTRACT(global, x[CORNER_OF_SIDE(theElement,side,0)], b);
4995 det = M2_DET(M);
4996
4997 return det;
4998 }
4999 #else
DistanceFromSide(const DOUBLE * global,const ELEMENT * theElement,INT side)5000 DOUBLE NS_DIM_PREFIX DistanceFromSide(const DOUBLE *global, const ELEMENT *theElement, INT side)
5001 {
5002 INT n;
5003 DOUBLE *x[MAX_CORNERS_OF_ELEM];
5004 DOUBLE M[DIM*DIM];
5005 DOUBLE *a, *b, *c;
5006 DOUBLE det;
5007
5008 a = &M[0];
5009 b = &M[DIM];
5010 c = &M[2*DIM];
5011
5012 CORNER_COORDINATES(theElement,n,x);
5013
5014 V3_SUBTRACT(x[CORNER_OF_SIDE(theElement,side,1)], x[CORNER_OF_SIDE(theElement,side,0)], a);
5015 V3_SUBTRACT(x[CORNER_OF_SIDE(theElement,side,2)], x[CORNER_OF_SIDE(theElement,side,0)], b);
5016 V3_SUBTRACT(global, x[CORNER_OF_SIDE(theElement,side,0)], c);
5017 det = M3_DET(M);
5018
5019 return det;
5020 }
5021 #endif
5022
5023 /****************************************************************************/
5024 /** \brief Find element containing position
5025
5026 * @param theGrid - grid level to search
5027 * @param pos - given position
5028
5029 This function finds the first element containing the position `pos`.
5030
5031 @return <ul>
5032 <li> pointer to ELEMENT </li>
5033 <li> NULL if not found. </li>
5034 </ul> */
5035 /****************************************************************************/
5036
FindElementOnSurface(MULTIGRID * theMG,DOUBLE * global)5037 ELEMENT * NS_DIM_PREFIX FindElementOnSurface (MULTIGRID *theMG, DOUBLE *global)
5038 {
5039 ELEMENT *t;
5040 INT k;
5041
5042 for (k=0; k<=TOPLEVEL(theMG); k++)
5043 for (t=FIRSTELEMENT(GRID_ON_LEVEL(theMG,k)); t!=NULL; t=SUCCE(t))
5044 if (EstimateHere(t))
5045 if (PointInElement(global,t)) return(t);
5046
5047 return(NULL);
5048 }
5049
5050 /****************************************************************************/
5051 /** \todo Please doc me!
5052 InnerBoundary -
5053
5054 SYNOPSIS:
5055 INT InnerBoundary (ELEMENT *t, INT side);
5056
5057
5058 * @param t
5059 * @param side
5060
5061 DESCRIPTION:
5062
5063 @return
5064 INT
5065 */
5066 /****************************************************************************/
5067
InnerBoundary(ELEMENT * t,INT side)5068 INT NS_DIM_PREFIX InnerBoundary (ELEMENT *t, INT side)
5069 {
5070 INT left,right,part;
5071
5072 ASSERT(OBJT(t) == BEOBJ);
5073 ASSERT(SIDE_ON_BND(t,side));
5074
5075 BNDS_BndSDesc(ELEM_BNDS(t,side),&left,&right,&part);
5076
5077 return((left != 0) && (right != 0));
5078 }
5079
5080
5081 /****************************************************************************/
5082 /** \brief Calculate the center of mass for an element
5083 *
5084 * @param theElement the element
5085 * @param center_of_mass center of mass as the result
5086 *
5087 * This function calculates the center of mass for an arbitrary element.
5088 * DOUBLE_VECTOR is an array for a 2D resp. 3D coordinate.
5089 *
5090 * \sa DOUBLE_VECTOR, ELEMENT
5091 */
5092 /****************************************************************************/
5093
CalculateCenterOfMass(ELEMENT * theElement,DOUBLE_VECTOR center_of_mass)5094 void NS_DIM_PREFIX CalculateCenterOfMass(ELEMENT *theElement, DOUBLE_VECTOR center_of_mass)
5095 {
5096 DOUBLE *corner;
5097 INT i, nr_corners;
5098
5099 nr_corners = CORNERS_OF_ELEM(theElement);
5100 V_DIM_CLEAR(center_of_mass);
5101
5102 for (i=0; i<nr_corners; i++)
5103 {
5104 corner = CVECT(MYVERTEX(CORNER(theElement,i)));
5105 V_DIM_ADD(center_of_mass,corner,center_of_mass);
5106 }
5107
5108 V_DIM_SCALE(1.0/nr_corners,center_of_mass);
5109 }
5110
5111 /****************************************************************************/
5112 /** \brief Calculate an (hopefully) unique key for the geometric object
5113
5114 * @param obj - geometric object which from the key is needed (can be one of VERTEX, ELEMENT, NODE or VECTOR)
5115
5116 This function calculates an (hopefully) unique key for VERTEX,
5117 ELEMENT, NODE, EDGE and VECTOR typed objects.
5118
5119 The heuristic is: calculate a 2D/3D position for the geometric object and
5120 transform this position to a single number by a weighted summation of the
5121 leading digits of the 2 resp. 3 coordinates and taking from this again
5122 the sigificant digits and adding the level number.
5123
5124 APPLICATION:
5125 Use always an explicit cast to avoid compiler warnings, e.g.
5126 NODE *theNode;
5127 KeyForObject((KEY_OBJECT *)theNode);
5128
5129 * \sa VERTEX, ELEMENT, NODE, EDGE, VECTOR
5130
5131 @return
5132 * the resulting key
5133 */
5134 /****************************************************************************/
5135
KeyForObject(KEY_OBJECT * obj)5136 INT NS_DIM_PREFIX KeyForObject( KEY_OBJECT *obj )
5137 {
5138 int dummy,i; /* dummy variable */
5139 DOUBLE_VECTOR coord;
5140
5141 if (obj==NULL) return (-1);
5142 switch( OBJT(obj) )
5143 {
5144 /* vertex */
5145 case BVOBJ :
5146 case IVOBJ : /* both together cover all vertex types */
5147 return LEVEL(obj)+COORDINATE_TO_KEY(CVECT((VERTEX*)obj),&dummy);
5148
5149 /* element */
5150 case BEOBJ :
5151 case IEOBJ : for (i=0; i<CORNERS_OF_ELEM((ELEMENT*)obj); i++)
5152 {
5153 if(CORNER((ELEMENT*)obj,i)==NULL)
5154 return (-1);
5155 if(MYVERTEX(CORNER((ELEMENT*)obj,i))==NULL)
5156 return (-1);
5157 }
5158 /* both together cover all element types */
5159 CalculateCenterOfMass( (ELEMENT*)obj, coord );
5160 return LEVEL(obj)+COORDINATE_TO_KEY(coord,&dummy);
5161
5162 /* node */
5163 case NDOBJ : if( MYVERTEX((NODE*)obj) == NULL )
5164 return (-1);
5165 return LEVEL(obj)+COORDINATE_TO_KEY(CVECT(MYVERTEX((NODE*)obj)),&dummy);
5166
5167 /* vector */
5168 case VEOBJ : if( VOBJECT((VECTOR*)obj) == NULL )
5169 return (-1);
5170 VectorPosition( (VECTOR*)obj, coord );
5171 return LEVEL(obj)+COORDINATE_TO_KEY(coord,&dummy);
5172
5173 /* edge */
5174 case EDOBJ : if( NBNODE(LINK0((EDGE*)obj)) == NULL )
5175 return (-1);
5176 if( MYVERTEX(NBNODE(LINK0((EDGE*)obj))) == NULL )
5177 return (-1);
5178 if( NBNODE(LINK1((EDGE*)obj)) == NULL )
5179 return (-1);
5180 if( MYVERTEX(NBNODE(LINK1((EDGE*)obj))) == NULL )
5181 return (-1);
5182 V_DIM_CLEAR(coord);
5183 /* sum of the coordinates of the 2 edge corners */
5184 V_DIM_ADD(coord,CVECT(MYVERTEX(NBNODE(LINK0((EDGE*)obj)))),coord);
5185 V_DIM_ADD(coord,CVECT(MYVERTEX(NBNODE(LINK1((EDGE*)obj)))),coord);
5186 /* the midpoint of the line is half of the sum */
5187 V_DIM_SCALE(0.5,coord);
5188 /* return the key of the midpoint as the key for the edge */
5189 return LEVEL(obj)+COORDINATE_TO_KEY(coord,&dummy);
5190
5191 default : sprintf( buffer, "unrecognized object type %d", OBJT(obj) );
5192 PrintErrorMessage('E',"KeyForObject",buffer);
5193 return(0);
5194 assert(0);
5195 }
5196 return (GM_ERROR);
5197 }
5198
ListMultiGridHeader(const INT longformat)5199 void NS_DIM_PREFIX ListMultiGridHeader (const INT longformat)
5200 {
5201 if (longformat)
5202 sprintf(buffer," %-20.20s %-20.20s %-20.20s %10.10s %10.10s\n","mg name","domain name","problem name","heap size","heap used");
5203 else
5204 sprintf(buffer," %-20.20s\n","mg name");
5205 }
5206
5207 /****************************************************************************/
5208 /** \brief List general information about multigrid structure
5209
5210 * @param theMG - structure to list
5211 * @param isCurrent - is `theMG` current multigrid
5212 * @param longformat - print all information or only name of `theMG`
5213
5214 This function lists general information about a multigrid structure.
5215
5216 */
5217 /****************************************************************************/
5218
ListMultiGrid(const MULTIGRID * theMG,const INT isCurrent,const INT longformat)5219 void NS_DIM_PREFIX ListMultiGrid (const MULTIGRID *theMG, const INT isCurrent, const INT longformat)
5220 {
5221 char c;
5222 const BVP_DESC *theBVPDesc;
5223
5224 /* get BVP description */
5225 theBVPDesc = MG_BVPD(theMG);
5226
5227 c = isCurrent ? '*' : ' ';
5228
5229 if (longformat)
5230 UserWriteF(" %c %-20.20s %-20.20s\n",c,ENVITEM_NAME(theMG),
5231 BVPD_NAME(theBVPDesc));
5232 else
5233 UserWriteF(" %c %-20.20s\n",c,ENVITEM_NAME(theMG));
5234
5235 return;
5236 }
5237
5238 /****************************************************************************/
5239 /** \brief List information about refinement type distribution
5240
5241 * @param theMG - structure to list
5242 * @param gridflag -
5243 * @param greenflag
5244 * @param lbflag
5245 * @param verbose
5246
5247 This function lists information about multigrids element types.
5248
5249 * \todo Please return value!
5250 */
5251 /****************************************************************************/
5252
MultiGridStatus(const MULTIGRID * theMG,INT gridflag,INT greenflag,INT lbflag,INT verbose)5253 INT NS_DIM_PREFIX MultiGridStatus (const MULTIGRID *theMG, INT gridflag, INT greenflag, INT lbflag, INT verbose)
5254 {
5255 INT i,j,sons,maxsons;
5256 INT red, green, yellow;
5257 INT mg_red,mg_green,mg_yellow;
5258 INT mg_greenrulesons[MAXLEVEL+1][MAX_SONS+1],mg_greenrules[MAXLEVEL+1];
5259 INT markcount[MAXLEVEL+1],closuresides[MAXLEVEL+1];
5260 FLOAT sum,sum_div_red,redplusgreen_div_red;
5261 FLOAT mg_sum,mg_sum_div_red,mg_redplusgreen_div_red;
5262 ELEMENT *theElement;
5263 GRID *theGrid;
5264 #ifdef ModelP
5265 INT MarkKey;
5266 INT total_elements,sum_elements;
5267 INT master_elements,hghost_elements,vghost_elements,vhghost_elements;
5268 #endif
5269
5270 const auto& ppifContext = theMG->ppifContext();
5271 const int me = ppifContext.me();
5272 const int procs = ppifContext.procs();
5273
5274 mg_red = mg_green = mg_yellow = mg_sum = 0;
5275 mg_sum_div_red = mg_redplusgreen_div_red = 0.0;
5276
5277 for (i=0; i<MAXLEVEL+1; i++)
5278 {
5279 /* for greenflag */
5280 mg_greenrules[i] = 0;
5281 for (j=0; j<MAX_SONS+1; j++) mg_greenrulesons[i][j] = 0;
5282 maxsons = 0;
5283
5284 /* for gridflag */
5285 markcount[i] = 0;
5286 closuresides[i] = 0;
5287 }
5288
5289 #ifdef ModelP
5290 MarkTmpMem(MGHEAP(theMG),&MarkKey);
5291 std::vector<int> infobuffer((procs+1)*(MAXLEVEL+1)*ELEMENT_PRIOS, 0);
5292 std::vector<int*> lbinfo(procs+1);
5293 for (i=0; i<procs+1; i++)
5294 lbinfo[i] = &infobuffer[i*(MAXLEVEL+1)*ELEMENT_PRIOS];
5295
5296 total_elements = sum_elements = 0;
5297 master_elements = hghost_elements = vghost_elements = vhghost_elements = 0;
5298 #endif
5299
5300 if (verbose && gridflag)
5301 {
5302 UserWriteF("\nMULTIGRID STATISTICS:\n");
5303 UserWriteF("LEVEL RED GREEN YELLOW SUM SUM/RED (RED+GREEN)/RED\n");
5304 }
5305
5306 /* compute multi grid infos */
5307 for (i=0; i<=TOPLEVEL(theMG); i++)
5308 {
5309 theGrid = GRID_ON_LEVEL(theMG,i);
5310 red = green = yellow = 0;
5311 sum = sum_div_red = redplusgreen_div_red = 0.0;
5312
5313 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL;
5314 theElement=SUCCE(theElement))
5315 {
5316 SETUSED(theElement,0);
5317 /* count eclasses */
5318 switch (ECLASS(theElement))
5319 {
5320 case RED_CLASS : red++; break;
5321 case GREEN_CLASS : green++; break;
5322 case YELLOW_CLASS : yellow++; break;
5323 default : assert(0);
5324 }
5325 /* count marks and closuresides */
5326 if (EstimateHere(theElement))
5327 {
5328 ELEMENT *MarkElement = ELEMENT_TO_MARK(theElement);
5329 INT marktype = GetRefinementMarkType(theElement);
5330
5331 if (marktype==1 &&
5332 USED(MarkElement)==0)
5333 {
5334 markcount[LEVEL(MarkElement)]++;
5335 markcount[MAXLEVEL]++;
5336 for (j=0; j<SIDES_OF_ELEM(MarkElement); j++)
5337 {
5338 ELEMENT *NbElement = NBELEM(MarkElement,j);
5339 if (NbElement!=NULL && MARKCLASS(NbElement)==RED_CLASS)
5340 {
5341 closuresides[LEVEL(MarkElement)]++;
5342 closuresides[MAXLEVEL]++;
5343 }
5344 }
5345 SETUSED(MarkElement,1);
5346 }
5347 }
5348 /* green refinement statistics */
5349 switch (REFINECLASS(theElement))
5350 {
5351 case GREEN_CLASS :
5352 sons = NSONS(theElement);
5353 mg_greenrulesons[i][sons]++;
5354 mg_greenrulesons[i][MAX_SONS]+=sons;
5355 mg_greenrules[i]++;
5356 mg_greenrulesons[MAXLEVEL][sons]++;
5357 mg_greenrulesons[MAXLEVEL][MAX_SONS]+=sons;
5358 mg_greenrules[MAXLEVEL]++;
5359 if (maxsons < sons) maxsons = sons;
5360 break;
5361 default :
5362 break;
5363 }
5364 #ifdef ModelP
5365 /* count master, hghost, vghost and vhghost elements */
5366 switch (EPRIO(theElement))
5367 {
5368 case PrioMaster :
5369 lbinfo[me][ELEMENT_PRIOS*i]++;
5370 lbinfo[me][ELEMENT_PRIOS*MAXLEVEL]++;
5371 break;
5372 case PrioHGhost :
5373 lbinfo[me][ELEMENT_PRIOS*i+1]++;
5374 lbinfo[me][ELEMENT_PRIOS*MAXLEVEL+1]++;
5375 break;
5376 case PrioVGhost :
5377 lbinfo[me][ELEMENT_PRIOS*i+2]++;
5378 lbinfo[me][ELEMENT_PRIOS*MAXLEVEL+2]++;
5379 break;
5380 case PrioVHGhost :
5381 lbinfo[me][ELEMENT_PRIOS*i+3]++;
5382 lbinfo[me][ELEMENT_PRIOS*MAXLEVEL+3]++;
5383 break;
5384 default :
5385 printf( PFMT "MultiGridStatus: wrong element prio %d\n",me,EPRIO(theElement));
5386 assert(0);
5387 break;
5388 }
5389 #endif
5390 }
5391 sum = red + green + yellow;
5392 if (red > 0)
5393 {
5394 sum_div_red = sum / red;
5395 redplusgreen_div_red = ((float)(red+green)) / red;
5396 }
5397 else
5398 {
5399 sum_div_red = 0.0;
5400 redplusgreen_div_red = 0.0;
5401 }
5402
5403 if (verbose && gridflag)
5404 UserWriteF(" %2d %9d %9d %9d %9.0f %2.3f %2.3f\n",
5405 i,red,green,yellow,sum,sum_div_red,redplusgreen_div_red);
5406
5407 mg_red += red;
5408 mg_green += green;
5409 mg_yellow += yellow;
5410 mg_sum += sum;
5411 }
5412 if (mg_red > 0)
5413 {
5414 mg_sum_div_red = mg_sum / mg_red;
5415 mg_redplusgreen_div_red = ((float)(mg_red + mg_green)) / mg_red;
5416 }
5417 else
5418 {
5419 mg_sum_div_red = 0.0;
5420 mg_redplusgreen_div_red = 0.0;
5421 }
5422
5423 if (verbose && gridflag)
5424 UserWriteF(" ALL %9d %9d %9d %9.0f %2.3f %2.3f\n",
5425 mg_red,mg_green,mg_yellow,mg_sum,mg_sum_div_red,mg_redplusgreen_div_red);
5426
5427 /* set heap info in refine info */
5428 if (gridflag)
5429 {
5430 float New;
5431 float newpergreen;
5432 float predmax;
5433
5434 SETMARKCOUNT(REFINEINFO(theMG),markcount[MAXLEVEL]);
5435
5436 New = markcount[MAXLEVEL]*(2<<(DIM-1))*mg_sum_div_red;
5437 SETPREDNEW0(REFINEINFO(theMG),New);
5438
5439 if (mg_greenrules[MAXLEVEL] > 0)
5440 newpergreen = ((float)mg_greenrulesons[MAXLEVEL][MAX_SONS])/mg_greenrules[MAXLEVEL];
5441 else
5442 newpergreen = 0;
5443 New = markcount[MAXLEVEL]*(2<<(DIM-1))+newpergreen*closuresides[MAXLEVEL];
5444 SETPREDNEW1(REFINEINFO(theMG),New);
5445
5446 SETREAL(REFINEINFO(theMG),mg_sum);
5447 SETPREDMAX(REFINEINFO(theMG),predmax);
5448 }
5449
5450 /* list heap info */
5451 if (verbose && gridflag)
5452 {
5453 UserWriteF(" EST %2d ELEMS=%9.0f MARKCOUNT=%9.0f PRED_NEW0=%9.0f PRED_NEW1=%9.0f PRED_MAX=%9.0f\n",
5454 REFINESTEP(REFINEINFO(theMG)),REAL(REFINEINFO(theMG)),MARKCOUNT(REFINEINFO(theMG)),
5455 PREDNEW0(REFINEINFO(theMG)),PREDNEW1(REFINEINFO(theMG)),PREDMAX(REFINEINFO(theMG)));
5456 UserWriteF(" EST TRACE step=%d\n",refine_info.step);
5457 for (i=0; i<refine_info.step; i++)
5458 UserWriteF(" EST %2d ELEMS=%9.0f MARKS=%9.0f REAL=%9.0f PRED0=%9.0f PRED1=%9.0f PRED_MAX=%9.0f\n",
5459 i,refine_info.real[i],refine_info.markcount[i],
5460 ((i<refine_info.step) ? refine_info.real[i+1]-refine_info.real[i] : 0),
5461 refine_info.predicted_new[i][0],
5462 refine_info.predicted_new[i][1],refine_info.predicted_max[i]);
5463 }
5464
5465 /* compute and list green rule info */
5466 if (verbose && greenflag)
5467 {
5468 UserWriteF("\nGREEN RULE STATISTICS:\n");
5469 UserWriteF(" LEVEL GREENSONS RULES GREENSONS/RUL");
5470 for (j=0; j<8 && j<maxsons; j++) UserWriteF(" %1d/%2d/...",j,j+8);
5471 UserWriteF("\n");
5472
5473 for (i=0; i<=TOPLEVEL(theMG); i++)
5474 {
5475 UserWriteF(" %2d %9d %9d %2.3f",i,mg_greenrulesons[i][MAX_SONS],
5476 mg_greenrules[i],
5477 (mg_greenrules[i]!=0) ? ((float)mg_greenrulesons[i][MAX_SONS])/mg_greenrules[i] : 0);
5478 for (j=0; j<maxsons; j++)
5479 {
5480 UserWriteF(" %9d",mg_greenrulesons[i][j]);
5481 if ((j+1)%8 == 0) UserWriteF("\n%41s"," ");
5482 }
5483 UserWriteF("\n");
5484 }
5485 UserWriteF(" ALL %9d %9d %2.3f",mg_greenrulesons[MAXLEVEL][MAX_SONS],
5486 mg_greenrules[MAXLEVEL],
5487 (mg_greenrules[MAXLEVEL]) ? ((float)mg_greenrulesons[MAXLEVEL][MAX_SONS])/
5488 mg_greenrules[MAXLEVEL] : 0.0);
5489 for (j=0; j<maxsons; j++)
5490 {
5491 UserWriteF(" %9d",mg_greenrulesons[MAXLEVEL][j]);
5492 if ((j+1)%8 == 0) UserWriteF("\n%41s"," ");
5493 }
5494 UserWriteF("\n");
5495 }
5496
5497 #ifdef ModelP
5498 /* compute and list loadbalancing info */
5499 if (verbose && lbflag)
5500 {
5501 UserWriteF("\nLB INFO:\n");
5502 /* now collect lb info on master */
5503 if (ppifContext.isMaster())
5504 {
5505 std::vector<VChannelPtr> mych(procs, nullptr);
5506
5507 for (i=1; i<procs; i++)
5508 {
5509 mych[i] =ConnSync(ppifContext, i,3917);
5510 RecvSync(ppifContext, mych[i],(void *)lbinfo[i],(MAXLEVEL+1)*ELEMENT_PRIOS*sizeof(INT));
5511 }
5512 Synchronize(ppifContext);
5513 for (i=1; i<procs; i++)
5514 {
5515 DiscSync(ppifContext, mych[i]);
5516 }
5517 }
5518 else
5519 {
5520 VChannelPtr mych;
5521
5522 mych = ConnSync(ppifContext, ppifContext.master(), 3917);
5523 SendSync(ppifContext, mych,(void *)lbinfo[me],(MAXLEVEL+1)*ELEMENT_PRIOS*sizeof(INT));
5524 Synchronize(ppifContext);
5525 DiscSync(ppifContext, mych);
5526 ReleaseTmpMem(MGHEAP(theMG),MarkKey);
5527 return(GM_OK);
5528 }
5529
5530 /* sum levels over procs */
5531 for (i=0; i<procs; i++)
5532 {
5533 for (j=0; j<TOPLEVEL(theMG)+1; j++)
5534 {
5535 lbinfo[procs][ELEMENT_PRIOS*j] += lbinfo[i][ELEMENT_PRIOS*j];
5536 lbinfo[procs][ELEMENT_PRIOS*j+1] += lbinfo[i][ELEMENT_PRIOS*j+1];
5537 lbinfo[procs][ELEMENT_PRIOS*j+2] += lbinfo[i][ELEMENT_PRIOS*j+2];
5538 lbinfo[procs][ELEMENT_PRIOS*j+3] += lbinfo[i][ELEMENT_PRIOS*j+3];
5539 }
5540 }
5541
5542 /* only master */
5543 if (lbflag >= 3)
5544 {
5545 UserWriteF(" LEVEL");
5546 for (i=0; i<ELEMENT_PRIOS*(TOPLEVEL(theMG)+1); i++)
5547 {
5548 UserWriteF(" %9d",i/ELEMENT_PRIOS);
5549 }
5550 UserWrite("\n");
5551 UserWriteF("PROC ");
5552 for (i=0; i<ELEMENT_PRIOS*(TOPLEVEL(theMG)+1); i++)
5553 {
5554 UserWriteF(" %9s",(i%ELEMENT_PRIOS==0) ? "MASTER" : (i%ELEMENT_PRIOS==1) ? "HGHOST" :
5555 (i%ELEMENT_PRIOS==2) ? "VGHOST" : "VHGHOST");
5556 }
5557 UserWrite("\n");
5558 for (i=0; i<procs; i++)
5559 {
5560 UserWriteF("%4d ",i);
5561 for (j=0; j<ELEMENT_PRIOS*(TOPLEVEL(theMG)+1); j++)
5562 {
5563 UserWriteF(" %9d",lbinfo[i][j]);
5564 }
5565 UserWrite("\n");
5566 }
5567 UserWriteF("\n");
5568 }
5569
5570 if (lbflag >= 2)
5571 {
5572 float memeff;
5573
5574 UserWriteF("%5s %9s %9s %9s %9s %9s %6s\n",
5575 "LEVEL","SUM","MASTER","HGHOST","VGHOST","VHGHOST","MEMEFF");
5576 for (i=0; i<=TOPLEVEL(theMG); i++)
5577 {
5578 sum_elements = lbinfo[procs][ELEMENT_PRIOS*i]+lbinfo[procs][ELEMENT_PRIOS*i+1]+
5579 lbinfo[procs][ELEMENT_PRIOS*i+2]+lbinfo[procs][ELEMENT_PRIOS*i+3];
5580 if (sum_elements > 0)
5581 memeff = ((float)lbinfo[procs][ELEMENT_PRIOS*i])/sum_elements*100;
5582 else
5583 memeff = 0.0;
5584 UserWriteF("%4d %9d %9d %9d %9d %9d %3.2f\n",i,sum_elements,
5585 lbinfo[procs][ELEMENT_PRIOS*i],lbinfo[procs][ELEMENT_PRIOS*i+1],
5586 lbinfo[procs][ELEMENT_PRIOS*i+2],lbinfo[procs][ELEMENT_PRIOS*i+3],memeff);
5587 }
5588 UserWrite("\n");
5589
5590 UserWriteF("%4s %9s %9s %9s %9s %9s %6s\n",
5591 "PROC","SUM","MASTER","HGHOST","VGHOST","VHGHOST","MEMEFF");
5592 for (i=0; i<procs; i++)
5593 {
5594 sum_elements = lbinfo[i][ELEMENT_PRIOS*MAXLEVEL]+lbinfo[i][ELEMENT_PRIOS*MAXLEVEL+1]+
5595 lbinfo[i][ELEMENT_PRIOS*MAXLEVEL+2]+lbinfo[i][ELEMENT_PRIOS*MAXLEVEL+3];
5596 if (sum_elements > 0)
5597 memeff = ((float)lbinfo[i][ELEMENT_PRIOS*MAXLEVEL])/sum_elements*100;
5598 else
5599 memeff = 0.0;
5600 UserWriteF("%4d %9d %9d %9d %9d %9d %3.2f\n",i,sum_elements,
5601 lbinfo[i][ELEMENT_PRIOS*MAXLEVEL],lbinfo[i][ELEMENT_PRIOS*MAXLEVEL+1],
5602 lbinfo[i][ELEMENT_PRIOS*MAXLEVEL+2],lbinfo[i][ELEMENT_PRIOS*MAXLEVEL+3],memeff);
5603 }
5604 UserWrite("\n");
5605 }
5606
5607 if (lbflag >= 1)
5608 {
5609 float memeff;
5610
5611 for (i=0; i<procs; i++)
5612 {
5613 master_elements += lbinfo[i][ELEMENT_PRIOS*MAXLEVEL];
5614 hghost_elements += lbinfo[i][ELEMENT_PRIOS*MAXLEVEL+1];
5615 vghost_elements += lbinfo[i][ELEMENT_PRIOS*MAXLEVEL+2];
5616 vhghost_elements += lbinfo[i][ELEMENT_PRIOS*MAXLEVEL+3];
5617 }
5618 total_elements = master_elements + hghost_elements + vghost_elements;
5619 if (total_elements > 0)
5620 memeff = ((float)master_elements)/total_elements*100;
5621 else
5622 memeff = 0.0;
5623 UserWriteF("%9s %9s %9s %9s %9s %6s\n","TOTAL","MASTER","HGHOST","VGHOST","VHGHOST","MEMEFF");
5624 UserWriteF("%9d %9d %9d %9d %9d %3.2f\n",total_elements,master_elements,hghost_elements,
5625 vghost_elements,vhghost_elements,memeff);
5626 }
5627
5628 }
5629 ReleaseTmpMem(MGHEAP(theMG),MarkKey);
5630 #endif
5631
5632 return (GM_OK);
5633 }
5634
5635 /****************************************************************************/
5636 /** \brief List general information about grids of multigrid
5637
5638 * @param theMG - multigrid structure
5639
5640 This function lists general information about the grids of a multigrid.
5641
5642 */
5643 /****************************************************************************/
5644
ListGrids(const MULTIGRID * theMG)5645 void NS_DIM_PREFIX ListGrids (const MULTIGRID *theMG)
5646 {
5647 GRID *theGrid;
5648 ELEMENT *theElement,*NBElem;
5649 VERTEX *myVertex,*nbVertex,*v0,*v1;
5650 NODE *theNode,*n0,*n1;
5651 EDGE *theEdge;
5652 LINK *theLink;
5653 VECTOR *vec;
5654 MATRIX *mat;
5655 char c;
5656 DOUBLE hmin,hmax,h;
5657 INT l,cl,minl,i,soe,eos,coe,side,e;
5658 INT nn,ne,nt,ns,nvec,nc;
5659
5660 cl = CURRENTLEVEL(theMG);
5661
5662 UserWriteF("grids of '%s':\n",ENVITEM_NAME(theMG));
5663
5664 UserWrite("level maxlevel #vert #node #edge #elem #side #vect #conn");
5665 UserWrite(" minedge maxedge\n");
5666
5667 for (l=0; l<=TOPLEVEL(theMG); l++)
5668 {
5669 theGrid = GRID_ON_LEVEL(theMG,l);
5670
5671 c = (l==cl) ? '*' : ' ';
5672
5673 /* calculate minimal and maximal edge */
5674 hmin = MAX_C;
5675 hmax = 0.0;
5676 for (theNode=FIRSTNODE(theGrid); theNode!=NULL; theNode=SUCCN(theNode))
5677 {
5678 myVertex = MYVERTEX(theNode);
5679 for (theLink=START(theNode); theLink!=NULL; theLink=NEXT(theLink))
5680 {
5681 nbVertex = MYVERTEX(NBNODE(theLink));
5682 V_DIM_EUKLIDNORM_OF_DIFF(CVECT(myVertex),CVECT(nbVertex),h);
5683 hmin = MIN(hmin,h);
5684 hmax = MAX(hmax,h);
5685 }
5686 }
5687 ns = 0;
5688 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL;
5689 theElement=SUCCE(theElement))
5690 if (OBJT(theElement) == BEOBJ)
5691 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
5692 if (SIDE_ON_BND(theElement,i))
5693 ns++;
5694
5695 UserWriteF("%c %3d %8d %8ld %8ld %8ld %8ld %8ld %8ld %8ld %9.3e %9.3e\n",c,l,(int)TOPLEVEL(theMG),
5696 (long)NV(theGrid),(long)NN(theGrid),(long)NE(theGrid),(long)NT(theGrid),
5697 (long)ns,(long)NVEC(theGrid),(long)NC(theGrid),(float)hmin,(float)hmax);
5698 }
5699
5700 /* surface grid up to current level */
5701 minl = cl;
5702 hmin = MAX_C;
5703 hmax = 0.0;
5704 nn = ne = nt = ns = nvec = nc = 0;
5705 for (l=0; l<=cl; l++)
5706 {
5707 theGrid = GRID_ON_LEVEL(theMG,l);
5708
5709 /* reset USED flags in all objects to be counted */
5710 for (theNode=FIRSTNODE(theGrid); theNode!=NULL; theNode=SUCCN(theNode))
5711 {
5712 SETUSED(theNode,0);
5713 for (theLink=START(theNode); theLink!=NULL; theLink=NEXT(theLink))
5714 SETUSED(MYEDGE(theLink),0);
5715 }
5716 for (vec=FIRSTVECTOR(theGrid); vec!=NULL; vec=SUCCVC(vec))
5717 for (mat=VSTART(vec); mat!=NULL; mat=MNEXT(mat))
5718 SETCUSED(MMYCON(mat),0);
5719
5720 /* count vectors and connections */
5721 for (vec=FIRSTVECTOR(theGrid); vec!=NULL; vec=SUCCVC(vec))
5722 if ((l==cl) || (VNCLASS(vec)<1))
5723 {
5724 nvec++;
5725 for (mat=VSTART(vec); mat!=NULL; mat=MNEXT(mat))
5726 {
5727 if (MUSED(mat)) continue;
5728 SETCUSED(MMYCON(mat),1);
5729
5730 if ((l==cl) || (VNCLASS(MDEST(mat))<1))
5731 nc++;
5732 }
5733 }
5734
5735 /* count other objects */
5736 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL; theElement=SUCCE(theElement))
5737 if ((NSONS(theElement)==0) || (l==cl))
5738 {
5739 nt++;
5740
5741 minl = MIN(minl,l);
5742
5743 coe = CORNERS_OF_ELEM(theElement);
5744 for (i=0; i<coe; i++)
5745 {
5746 theNode = CORNER(theElement,i);
5747 if (USED(theNode)) continue;
5748 SETUSED(theNode,1);
5749
5750 if ((SONNODE(theNode)==NULL) || (l==cl))
5751 nn++;
5752 }
5753
5754 soe = SIDES_OF_ELEM(theElement);
5755 for (side=0; side<soe; side++)
5756 {
5757 if (OBJT(theElement)==BEOBJ)
5758 if (ELEM_BNDS(theElement,side)!=NULL) ns++;
5759
5760 /* check neighbour element */
5761 if (l<cl)
5762 if ((NBElem=NBELEM(theElement,side))!=NULL)
5763 if (NSONS(NBElem)>0)
5764 continue; /* objects of this side will be counted by the neighbour */
5765
5766 eos = EDGES_OF_SIDE(theElement,side);
5767 for (i=0; i<eos; i++)
5768 {
5769 e = EDGE_OF_SIDE(theElement,side,i);
5770 n0 = CORNER(theElement,CORNER_OF_EDGE(theElement,e,0));
5771 v0 = MYVERTEX(n0);
5772 n1 = CORNER(theElement,CORNER_OF_EDGE(theElement,e,1));
5773 v1 = MYVERTEX(n1);
5774
5775 if ((theEdge=GetEdge(n0,n1))==NULL) continue;
5776 if (USED(theEdge)) continue;
5777 SETUSED(theEdge,1);
5778
5779 /* any sons ? */
5780 if (SONNODE(n0)!=NULL && SONNODE(n1)!=NULL)
5781 if (GetEdge(SONNODE(n0),SONNODE(n1))!=NULL) continue;
5782 if (MIDNODE(theEdge) != NULL)
5783 {
5784 if (SONNODE(n0)!=NULL)
5785 if (GetEdge(MIDNODE(theEdge),SONNODE(n0))!=NULL) continue;
5786 if (SONNODE(n1)!=NULL)
5787 if (GetEdge(MIDNODE(theEdge),SONNODE(n1))!=NULL) continue;
5788 }
5789 ne++;
5790
5791 V_DIM_EUKLIDNORM_OF_DIFF(CVECT(v0),CVECT(v1),h);
5792 hmin = MIN(hmin,h);
5793 hmax = MAX(hmax,h);
5794 }
5795 }
5796 }
5797 }
5798
5799 UserWrite("\nsurface grid up to current level:\n");
5800 UserWriteF("%c %3d %8d %8s %8ld %8ld %8ld %8ld %8ld %8ld %9.3e %9.3e\n",' ',minl,(int)cl,
5801 "---",(long)nn,(long)ne,(long)nt,
5802 (long)ns,(long)nvec,(long)nc,(float)hmin,(float)hmax);
5803
5804 #ifdef ModelP
5805 /* surface grid up to current level */
5806 minl = cl;
5807 hmin = MAX_C;
5808 hmax = 0.0;
5809 nn = ne = nt = ns = nvec = nc = 0;
5810 for (l=0; l<=cl; l++)
5811 {
5812 theGrid = GRID_ON_LEVEL(theMG,l);
5813
5814 /* reset USED flags in all objects to be counted */
5815 for (theNode=FIRSTNODE(theGrid); theNode!=NULL; theNode=SUCCN(theNode))
5816 {
5817 SETUSED(theNode,0);
5818 for (theLink=START(theNode); theLink!=NULL; theLink=NEXT(theLink))
5819 SETUSED(MYEDGE(theLink),0);
5820 }
5821 /* count vectors and connections */
5822 for (vec=FIRSTVECTOR(theGrid); vec!=NULL; vec=SUCCVC(vec))
5823 if ((l==cl) || (VNCLASS(vec)<1))
5824 if (PRIO(vec) == PrioMaster) nvec++;
5825
5826 /* count other objects */
5827 for (theElement=FIRSTELEMENT(theGrid);
5828 theElement!=NULL; theElement=SUCCE(theElement))
5829 if (EstimateHere(theElement))
5830 {
5831 nt++;
5832
5833 minl = MIN(minl,l);
5834
5835 coe = CORNERS_OF_ELEM(theElement);
5836 for (i=0; i<coe; i++)
5837 {
5838 theNode = CORNER(theElement,i);
5839 if (USED(theNode)) continue;
5840 SETUSED(theNode,1);
5841
5842 if ((SONNODE(theNode)==NULL) || (l==cl))
5843 if (PRIO(theNode) == PrioMaster) nn++;
5844 }
5845
5846 soe = SIDES_OF_ELEM(theElement);
5847 for (side=0; side<soe; side++)
5848 {
5849 if (OBJT(theElement)==BEOBJ)
5850 if (ELEM_BNDS(theElement,side)!=NULL) ns++;
5851
5852 /* check neighbour element */
5853 if (l<cl)
5854 if ((NBElem=NBELEM(theElement,side))!=NULL)
5855 if (NSONS(NBElem)>0)
5856 continue; /* objects of this side will be counted by the neighbour */
5857
5858 eos = EDGES_OF_SIDE(theElement,side);
5859 for (i=0; i<eos; i++)
5860 {
5861 e = EDGE_OF_SIDE(theElement,side,i);
5862 n0 = CORNER(theElement,CORNER_OF_EDGE(theElement,e,0));
5863 v0 = MYVERTEX(n0);
5864 n1 = CORNER(theElement,CORNER_OF_EDGE(theElement,e,1));
5865 v1 = MYVERTEX(n1);
5866 V_DIM_EUKLIDNORM_OF_DIFF(CVECT(v0),CVECT(v1),h);
5867 hmin = MIN(hmin,h);
5868 hmax = MAX(hmax,h);
5869 }
5870 }
5871 }
5872 }
5873 nn = UG_GlobalSumINT(theMG->ppifContext(), nn);
5874 ne = UG_GlobalSumINT(theMG->ppifContext(), ne);
5875 nt = UG_GlobalSumINT(theMG->ppifContext(), nt);
5876 ns = UG_GlobalSumINT(theMG->ppifContext(), ns);
5877 nvec = UG_GlobalSumINT(theMG->ppifContext(), nvec);
5878 nc = UG_GlobalSumINT(theMG->ppifContext(), nc);
5879 hmin = UG_GlobalMinDOUBLE(theMG->ppifContext(), hmin);
5880 hmax = UG_GlobalMaxDOUBLE(theMG->ppifContext(), hmax);
5881 UserWrite("\nsurface of all processors up to current level:\n");
5882 UserWriteF("%c %3d %8d %8s %8ld %8s %8ld %8ld %8ld %8s %9.3e %9.3e\n",
5883 ' ',minl,(int)cl,
5884 "---",(long)nn," ",(long)nt,
5885 (long)ns,(long)nvec," ",(float)hmin,(float)hmax);
5886 #endif
5887 }
5888
5889 /****************************************************************************/
5890 /** \brief List information about node in multigrid
5891
5892 * @param theMG - structure containing the node
5893 * @param theNode - node to list
5894 * @param dataopt - list user data if true
5895 * @param bopt - list boundary info if true
5896 * @param nbopt - list info about neighbors if true
5897 * @param vopt - list more information
5898
5899 This function lists information about a node in a multigrid.
5900
5901 */
5902 /****************************************************************************/
5903
ListNode(const MULTIGRID * theMG,const NODE * theNode,INT dataopt,INT bopt,INT nbopt,INT vopt)5904 void NS_DIM_PREFIX ListNode (const MULTIGRID *theMG, const NODE *theNode, INT dataopt, INT bopt, INT nbopt, INT vopt)
5905 {
5906 VERTEX *theVertex;
5907 LINK *theLink;
5908 INT i,part;
5909
5910 theVertex = MYVERTEX(theNode);
5911
5912 /******************************/
5913 /* print standard information */
5914 /******************************/
5915 /* line 1 */ UserWriteF("NODEID=" ID_FFMTE " CTRL=%8lx VEID="
5916 VID_FMTX " LEVEL=%2d",
5917 ID_PRTE(theNode),(long)CTRL(theNode),
5918 VID_PRTX(theVertex),LEVEL(theNode));
5919
5920 /* print coordinates of that node */
5921 for(i=0; i<DIM; i++)
5922 {
5923 UserWriteF(" x%1d=%11.4E",i, (float)(CVECT(theVertex)[i]) );
5924 }
5925 UserWrite("\n");
5926
5927 if (vopt) /* verbose: print all information */
5928 {
5929 /* print nfather information */
5930 if (NFATHER(theNode)!=NULL)
5931 {
5932 switch (NTYPE(theNode))
5933 {
5934 case (CORNER_NODE) :
5935 UserWriteF(" NFATHER(Node)=" ID_FMTX "\n",
5936 ID_PRTX((NODE *)NFATHER(theNode)));
5937 break;
5938 case (MID_NODE) :
5939 UserWriteF(" NFATHER(Edge)=" EDID_FMTX "\n",
5940 EDID_PRTX((EDGE *)NFATHER(theNode)));
5941 break;
5942 default :
5943 break;
5944 }
5945 }
5946 /* print nfather information */
5947 if (SONNODE(theNode)!=NULL)
5948 {
5949 UserWriteF(" SONNODE=" ID_FMTX "\n",ID_PRTX(SONNODE(theNode)));
5950 }
5951
5952 /* line 3 */ /* print vertex father information */
5953 if (VFATHER(theVertex)!=NULL)
5954 {
5955 UserWriteF(" VERTEXFATHER=" EID_FMTX " ",
5956 EID_PRTX(VFATHER(theVertex)));
5957 for(i=0; i<DIM; i++)
5958 {
5959 UserWriteF("XI[%d]=%11.4E ",i, (float)(LCVECT(theVertex)[i]) );
5960 }
5961 }
5962
5963 UserWriteF(" key=%d\n", KeyForObject((KEY_OBJECT *)theNode) );
5964
5965 if (NVECTOR(theNode) != NULL)
5966 UserWriteF(" vec=" VINDEX_FMTX "\n",
5967 VINDEX_PRTX(NVECTOR(theNode)));
5968
5969 UserWriteF(" classes: NCLASS = %d NNCLASS = %d\n",NCLASS(theNode),NNCLASS(theNode));
5970 }
5971
5972 /******************************/
5973 /* print boundary information */
5974 /******************************/
5975 if (bopt)
5976 {
5977 if (OBJT(theVertex) == BVOBJ)
5978 {
5979 if (BNDP_BndPDesc(V_BNDP(theVertex),&i,&part))
5980 UserWrite("Error in boundary point\n");
5981 else
5982 UserWriteF("boundary point: move %d moved %d\n",i,
5983 MOVED(theVertex));
5984 }
5985 }
5986
5987 if (nbopt)
5988 {
5989 for (theLink=START(theNode); theLink!=NULL; theLink=NEXT(theLink))
5990 {
5991 #if defined __THREEDIM__ && defined ModelP
5992 UserWriteF(" EDGE=%x/%08x ",MYEDGE(theLink),
5993 DDD_InfoGlobalId(PARHDR(MYEDGE(theLink))));
5994 #else
5995 UserWrite(" ");
5996 #endif
5997 UserWriteF("NB=" ID_FMTX " CTRL=%8lx NO_OF_ELEM=%3d",
5998 ID_PRTX(NBNODE(theLink)),(long)CTRL(theLink),NO_OF_ELEM(MYEDGE(theLink)));
5999 if (MIDNODE(MYEDGE(theLink))!=NULL)
6000 UserWriteF(" MIDNODE=" ID_FMTX, ID_PRTX(MIDNODE(MYEDGE(theLink))));
6001 theVertex=MYVERTEX(NBNODE(theLink));
6002 for(i=0; i<DIM; i++)
6003 {
6004 UserWriteF(" x%1d=%11.4E",i, (float)(CVECT(theVertex)[i]) );
6005 }
6006 UserWrite("\n");
6007 }
6008 }
6009 return;
6010 }
6011
6012
6013 /****************************************************************************/
6014 /** \brief
6015 ListElement - List information about element
6016
6017 * @param theMG - structure to list
6018 * @param theElement - element to list
6019 * @param dataopt - list user data if true
6020 * @param vopt - list more information
6021
6022 This function lists information about an element
6023
6024 */
6025 /****************************************************************************/
6026
ListElement(const MULTIGRID * theMG,const ELEMENT * theElement,INT dataopt,INT bopt,INT nbopt,INT vopt)6027 void NS_DIM_PREFIX ListElement (const MULTIGRID *theMG, const ELEMENT *theElement, INT dataopt, INT bopt, INT nbopt, INT vopt)
6028 {
6029 char etype[10];
6030 char ekind[8];
6031 int i,j;
6032 ELEMENT *SonList[MAX_SONS];
6033
6034 if (DIM==2)
6035 switch (TAG(theElement))
6036 {
6037 case TRIANGLE : strcpy(etype,"TRI"); break;
6038 case QUADRILATERAL : strcpy(etype,"QUA"); break;
6039 default : strcpy(etype,"???"); break;
6040 }
6041 else
6042 switch (TAG(theElement))
6043 {
6044 case TETRAHEDRON : strcpy(etype,"TET"); break;
6045 case PYRAMID : strcpy(etype,"PYR"); break;
6046 case PRISM : strcpy(etype,"PRI"); break;
6047 case HEXAHEDRON : strcpy(etype,"HEX"); break;
6048 default : strcpy(etype,"???"); break;
6049 }
6050 switch (ECLASS(theElement))
6051 {
6052 case YELLOW_CLASS : strcpy(ekind,"YELLOW "); break;
6053 case GREEN_CLASS : strcpy(ekind,"GREEN "); break;
6054 case RED_CLASS : strcpy(ekind,"RED "); break;
6055 default : strcpy(ekind,"??? "); break;
6056 }
6057 UserWriteF("ELEMID=" EID_FFMTE " %5s %5s CTRL=%8lx CTRL2=%8lx REFINE=%2d MARK=%2d LEVEL=%2d",
6058 EID_PRTE(theElement),ekind,etype,
6059 (long)CTRL(theElement),(long)FLAG(theElement),REFINE(theElement),MARK(theElement),LEVEL(theElement));
6060 if (COARSEN(theElement)) UserWrite(" COARSEN");
6061 UserWrite("\n");
6062
6063 if (vopt)
6064 {
6065 UserWriteF("subdomain=%d \n", SUBDOMAIN(theElement));
6066 for (i=0; i<CORNERS_OF_ELEM(theElement); i++)
6067 {
6068 UserWriteF(" N%d=" ID_FMTX,i,ID_PRTX(CORNER(theElement,i)));
6069 }
6070 UserWriteF("\n");
6071 if (EFATHER(theElement))
6072 UserWriteF(" FA=" EID_FMTX ,EID_PRTX(EFATHER(theElement)));
6073 else
6074 UserWriteF(" FA=NULL");
6075
6076 UserWriteF(" NSONS=%d\n",NSONS(theElement));
6077
6078 if (GetAllSons(theElement,SonList)!=0) return;
6079 for (i=0; SonList[i] != NULL; i++)
6080 {
6081 UserWriteF(" S%d=" EID_FMTX ,i,EID_PRTX(SonList[i]));
6082 if ((i+1)%4 == 0) UserWrite("\n");
6083 }
6084
6085 }
6086 if (nbopt)
6087 {
6088 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
6089 if (NBELEM(theElement,i)!=NULL)
6090 {
6091 UserWriteF(" NB%d=" EID_FMTX ,
6092 i,EID_PRTX(NBELEM(theElement,i)));
6093 }
6094 UserWrite("\n");
6095 }
6096 if (bopt)
6097 {
6098 UserWrite(" ");
6099 if (OBJT(theElement)==BEOBJ)
6100 {
6101 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
6102 {
6103 for(j=0; j<CORNERS_OF_SIDE(theElement,i); j++)
6104 {
6105 #if defined(ModelP) && defined(__THREEDIM__)
6106 UserWriteF(" NODE[ID=%ld]: ",
6107 (long)(ID(CORNER(theElement,
6108 CORNER_OF_SIDE(theElement,i,j)))));
6109 #endif
6110 UserWrite("\n");
6111 }
6112 }
6113 }
6114 UserWrite("\n");
6115 }
6116
6117 return;
6118 }
6119
6120
6121 /****************************************************************************/
6122 /** \brief List information about vector
6123
6124 * @param theMG - multigrid structure to list
6125 * @param theVector - vector to list
6126 * @param matrixopt - list line of matrix corresponding to theVector
6127 * @param dataopt - list user data if true
6128 * @param modifiers - flags modifying output style and verbose level
6129
6130 This function lists information about a vector.
6131
6132 */
6133 /****************************************************************************/
6134
ListVector(const MULTIGRID * theMG,const VECTOR * theVector,INT matrixopt,INT dataopt,INT modifiers)6135 void NS_DIM_PREFIX ListVector (const MULTIGRID *theMG, const VECTOR *theVector, INT matrixopt, INT dataopt, INT modifiers)
6136 {
6137 NODE *theNode;
6138 EDGE *theEdge;
6139 ELEMENT *theElement;
6140 MATRIX *theMatrix;
6141 DOUBLE_VECTOR pos;
6142
6143 FORMAT* theFormat = theMG->theFormat.get();
6144
6145 /* print index and type of vector */
6146 UserWriteF("IND=" VINDEX_FFMTE " VTYPE=%d(%c) ",
6147 VINDEX_PRTE(theVector),
6148 VTYPE(theVector),
6149 FMT_T2N(theFormat,VTYPE(theVector)));
6150
6151 if (READ_FLAG(modifiers,LV_POS))
6152 {
6153 if (VectorPosition(theVector,pos))
6154 return;
6155 #ifdef __TWODIM__
6156 UserWriteF("POS=(%10.2e,%10.2e)",pos[_X_],pos[_Y_]);
6157 #endif
6158 #ifdef __THREEDIM__
6159 UserWriteF("POS=(%10.2e,%10.2e,%10.2e)",pos[_X_],pos[_Y_],pos[_Z_]);
6160 #endif
6161 }
6162
6163 /* print object type of vector */
6164 if (READ_FLAG(modifiers,LV_VO_INFO))
6165 switch (VOTYPE(theVector))
6166 {
6167 case NODEVEC :
6168 {
6169 theNode = (NODE*)VOBJECT(theVector);
6170 #if defined __OVERLAP2__
6171 if ( theNode == NULL )
6172 UserWriteF("NODE-V NULL ");
6173 else
6174 #endif
6175 UserWriteF("NODE-V nodeID=" ID_FMTX
6176 " ",
6177 ID_PRTX(theNode));
6178 }
6179 break;
6180 case EDGEVEC :
6181 {
6182 theEdge = (EDGE*)VOBJECT(theVector);
6183 UserWriteF("EDGE-V fromID=" ID_FFMT
6184 " to__ID=%7ld ",
6185 ID_PRT(NBNODE(LINK0(theEdge))),
6186 ID(NBNODE(LINK1(theEdge))));
6187 }
6188 break;
6189 #ifdef __THREEDIM__
6190 case SIDEVEC :
6191 {
6192 theElement = (ELEMENT*)VOBJECT(theVector);
6193 UserWriteF("SIDE-V elemID=" EID_FFMT
6194 " ",
6195 EID_PRT(theElement));
6196 }
6197 break;
6198 #endif
6199 case ELEMVEC :
6200 {
6201 theElement = (ELEMENT*)VOBJECT(theVector);
6202 UserWriteF("ELEM-V elemID=" EID_FFMT
6203 " ",
6204 EID_PRT(theElement));
6205 }
6206 break;
6207
6208 default : PrintErrorMessage( 'E', "ListVector", "unrecognized VECTOR type" );
6209 assert(0);
6210 }
6211
6212 UserWriteF("VCLASS=%1d VNCLASS=%1d",VCLASS(theVector),VNCLASS(theVector));
6213 UserWriteF(" key=%d\n", KeyForObject((KEY_OBJECT *)theVector) );
6214
6215 /* print matrix list if */
6216 if (matrixopt > 0)
6217 for (theMatrix = VSTART(theVector); theMatrix!=NULL; theMatrix=MNEXT(theMatrix))
6218 {
6219 UserWrite(" DEST(MATRIX): ");
6220 ListVector(theMG,MDEST(theMatrix),0,0,modifiers);
6221 }
6222 }
6223
6224 /****************************************************************************/
6225 /** \brief Returns highest Node class of a dof on next level
6226
6227 * @param theElement - pointer to a element
6228
6229 This function returns highest 'NCLASS' of a Node associated with the
6230 element.
6231
6232 @return <ul>
6233 <li> 0 if ok </li>
6234 <li> 1 if error occured. </li>
6235 </ul> */
6236 /****************************************************************************/
6237
MaxNodeClass(const ELEMENT * theElement)6238 static INT MaxNodeClass (const ELEMENT *theElement)
6239 {
6240 INT m = 0;
6241 INT i;
6242
6243 for (i=0; i<CORNERS_OF_ELEM(theElement); i++) {
6244 INT c = NCLASS(CORNER(theElement,i));
6245
6246 m = MAX(m,c);
6247 }
6248
6249 return (m);
6250 }
6251
6252 /****************************************************************************/
6253 /** \brief Returns highest Node class of a dof on next level
6254
6255 * @param theElement - pointer to a element
6256
6257 This function returns highest 'NNCLASS' of a Node associated with the
6258 element.
6259
6260 @return <ul>
6261 <li> 0 if ok </li>
6262 <li> 1 if error occured. </li>
6263 </ul> */
6264 /****************************************************************************/
6265
MaxNextNodeClass(const ELEMENT * theElement)6266 INT NS_DIM_PREFIX MaxNextNodeClass (const ELEMENT *theElement)
6267 {
6268 INT m = 0;
6269 INT i;
6270
6271 for (i=0; i<CORNERS_OF_ELEM(theElement); i++) {
6272 INT c = NNCLASS(CORNER(theElement,i));
6273
6274 m = MAX(m,c);
6275 }
6276
6277 return (m);
6278 }
6279
6280 /****************************************************************************/
6281 /** \brief Returns minimal Node class of a dof on next level
6282
6283 * @param theElement - pointer to a element
6284
6285 This function returns highest 'NNCLASS' of a Node associated with the
6286 element.
6287
6288 @return <ul>
6289 <li> 0 if ok </li>
6290 <li> 1 if error occured. </li>
6291 </ul> */
6292 /****************************************************************************/
6293
MinNodeClass(const ELEMENT * theElement)6294 INT NS_DIM_PREFIX MinNodeClass (const ELEMENT *theElement)
6295 {
6296 INT m = 3;
6297 INT i;
6298
6299 for (i=0; i<CORNERS_OF_ELEM(theElement); i++) {
6300 INT c = NCLASS(CORNER(theElement,i));
6301
6302 m = MIN(m,c);
6303 }
6304
6305 return (m);
6306 }
6307
6308 /****************************************************************************/
6309 /** \brief Returns minimal Node class of a dof on next level
6310
6311 * @param theElement - pointer to a element
6312
6313 This function returns highest 'NNCLASS' of a Node associated with the
6314 element.
6315
6316 @return <ul>
6317 <li> 0 if ok </li>
6318 <li> 1 if error occured. </li>
6319 </ul> */
6320 /****************************************************************************/
6321
MinNextNodeClass(const ELEMENT * theElement)6322 INT NS_DIM_PREFIX MinNextNodeClass (const ELEMENT *theElement)
6323 {
6324 INT m = 3;
6325 INT i;
6326
6327 for (i=0; i<CORNERS_OF_ELEM(theElement); i++) {
6328 INT c = NNCLASS(CORNER(theElement,i));
6329
6330 m = MIN(m,c);
6331 }
6332
6333 return (m);
6334 }
6335
6336 /****************************************************************************/
6337 /** \brief Initialize node classes
6338
6339 * @param theGrid - given grid
6340 * @param theElement - given element
6341
6342 Initialize Node class in all nodes associated with given element with 3.
6343
6344 @return <ul>
6345 <li> 0 if ok </li>
6346 <li> 1 if error occured. </li>
6347 </ul> */
6348 /****************************************************************************/
6349
SeedNodeClasses(ELEMENT * theElement)6350 INT NS_DIM_PREFIX SeedNodeClasses (ELEMENT *theElement)
6351 {
6352 INT i;
6353
6354 for (i=0; i<CORNERS_OF_ELEM(theElement); i++)
6355 SETNCLASS(CORNER(theElement,i),3);
6356
6357 return (0);
6358 }
6359
6360 /****************************************************************************/
6361 /** \brief Reset node classes
6362
6363 * @param theGrid - pointer to grid
6364
6365 Reset all node classes in all nodes of given grid to 0.
6366
6367 @return <ul>
6368 <li> 0 if ok </li>
6369 <li> 1 if error occured. </li>
6370 </ul> */
6371 /****************************************************************************/
6372
ClearNodeClasses(GRID * theGrid)6373 INT NS_DIM_PREFIX ClearNodeClasses (GRID *theGrid)
6374 {
6375 NODE *theNode;
6376
6377 /* reset class of each Node to 0 */
6378 for (theNode=PFIRSTNODE(theGrid); theNode!=NULL; theNode=SUCCN(theNode))
6379 SETNCLASS(theNode,0);
6380
6381 return(0);
6382 }
6383
6384 #ifdef ModelP
Gather_NodeClass(DDD::DDDContext &,DDD_OBJ obj,void * data)6385 static int Gather_NodeClass (DDD::DDDContext&, DDD_OBJ obj, void *data)
6386 {
6387 NODE *theNode = (NODE *)obj;
6388
6389 ((INT *)data)[0] = NCLASS(theNode);
6390
6391 return(0);
6392 }
6393
Scatter_NodeClass(DDD::DDDContext &,DDD_OBJ obj,void * data)6394 static int Scatter_NodeClass (DDD::DDDContext&, DDD_OBJ obj, void *data)
6395 {
6396 NODE *theNode = (NODE *)obj;
6397
6398 SETNCLASS(theNode,MAX(NCLASS(theNode),((INT *)data)[0]));
6399
6400 return(0);
6401 }
6402
Scatter_GhostNodeClass(DDD::DDDContext &,DDD_OBJ obj,void * data)6403 static int Scatter_GhostNodeClass (DDD::DDDContext&, DDD_OBJ obj, void *data)
6404 {
6405 NODE *theNode = (NODE *)obj;
6406
6407 SETNCLASS(theNode,((INT *)data)[0]);
6408
6409 return(0);
6410 }
6411 #endif
6412
6413 /****************************************************************************/
6414 /** \brief Compute Node classes after initialization
6415
6416 * @param theGrid - pointer to grid
6417
6418 After Node classes have been reset and initialized, this function
6419 now computes the class 2 and class 1 Nodes.
6420
6421 @return <ul>
6422 <li> 0 if ok </li>
6423 <li> 1 if error occured </li>
6424 </ul> */
6425 /****************************************************************************/
PropagateNodeClass(GRID * theGrid,INT nclass)6426 static INT PropagateNodeClass (GRID *theGrid, INT nclass)
6427 {
6428 ELEMENT *theElement;
6429 INT i;
6430
6431 for (theElement=FIRSTELEMENT(theGrid);
6432 theElement!= NULL; theElement = SUCCE(theElement))
6433 if (MaxNodeClass(theElement) == nclass)
6434 for (i=0; i<CORNERS_OF_ELEM(theElement); i++)
6435 {
6436 NODE *theNode = CORNER(theElement,i);
6437
6438 if (NCLASS(theNode) < nclass)
6439 SETNCLASS(theNode,nclass-1);
6440 }
6441
6442 /* only for this values valid */
6443 ASSERT(nclass==3 || nclass==2);
6444
6445 return(0);
6446 }
6447
PropagateNodeClasses(GRID * theGrid)6448 INT NS_DIM_PREFIX PropagateNodeClasses (GRID *theGrid)
6449 {
6450 #ifdef ModelP
6451 auto& context = theGrid->dddContext();
6452 const auto& dddctrl = ddd_ctrl(context);
6453
6454 PRINTDEBUG(gm,1,("\n" PFMT "PropagateNodeClasses():"
6455 " 1. communication on level %d\n",theGrid->ppifContext().me(),GLEVEL(theGrid)))
6456 /* exchange NCLASS of Nodes */
6457 DDD_IFAExchange(context,
6458 dddctrl.BorderNodeSymmIF,GRID_ATTR(theGrid),sizeof(INT),
6459 Gather_NodeClass, Scatter_NodeClass);
6460 #endif
6461
6462 /* set Node classes in the algebraic neighborhood to 2 */
6463 if (PropagateNodeClass(theGrid,3)) REP_ERR_RETURN(1);
6464
6465 #ifdef ModelP
6466 PRINTDEBUG(gm,1,("\n" PFMT "PropagateNodeClasses(): 2. communication\n",
6467 theGrid->ppifContext().me()))
6468 /* exchange NCLASS of Nodes */
6469 DDD_IFAExchange(context,
6470 dddctrl.BorderNodeSymmIF,GRID_ATTR(theGrid),sizeof(INT),
6471 Gather_NodeClass, Scatter_NodeClass);
6472 #endif
6473
6474 /* set Node classes in the algebraic neighborhood to 1 */
6475 if (PropagateNodeClass(theGrid,2)) REP_ERR_RETURN(1);
6476
6477 #ifdef ModelP
6478 PRINTDEBUG(gm,1,("\n" PFMT "PropagateNodeClasses(): 3. communication\n",
6479 theGrid->ppifContext().me()))
6480 /* exchange NCLASS of Nodes */
6481 DDD_IFAExchange(context,
6482 dddctrl.BorderNodeSymmIF,GRID_ATTR(theGrid),sizeof(INT),
6483 Gather_NodeClass, Scatter_NodeClass);
6484
6485 /* send NCLASS to ghosts */
6486 DDD_IFAOneway(context,
6487 dddctrl.NodeIF,GRID_ATTR(theGrid),IF_FORWARD,sizeof(INT),
6488 Gather_NodeClass, Scatter_GhostNodeClass);
6489 #endif
6490
6491 return(0);
6492 }
6493
6494 /****************************************************************************/
6495 /** \brief Reset class of the Nodes on the next level
6496
6497 * @param theGrid - pointer to grid
6498
6499 This function clears NNCLASS flag in all Nodes. This is the first step to
6500 compute the class of the dofs on the *NEXT* level, which
6501 is also the basis for determining copies.
6502
6503 @return <ul>
6504 <li> 0 if ok </li>
6505 <li> 1 if error occured. </li>
6506 </ul> */
6507 /****************************************************************************/
6508
ClearNextNodeClasses(GRID * theGrid)6509 INT NS_DIM_PREFIX ClearNextNodeClasses (GRID *theGrid)
6510 {
6511 NODE *theNode;
6512
6513 /* reset class of each Node to 0 */
6514 for (theNode=PFIRSTNODE(theGrid); theNode!=NULL; theNode=SUCCN(theNode))
6515 SETNNCLASS(theNode,0);
6516
6517 /* now the refinement algorithm will initialize the class 3 Nodes */
6518 /* on the *NEXT* level. */
6519 return(0);
6520 }
6521
6522 /****************************************************************************/
6523 /** \brief Set 'NNCLASS' in all Nodes associated with element
6524
6525 * @param theElement - pointer to element
6526
6527 Set 'NNCLASS' in all nodes associated with the element to 3.
6528
6529 @return <ul>
6530 <li> 0 if ok </li>
6531 <li> 1 if error occured. </li>
6532 </ul> */
6533 /****************************************************************************/
6534
SeedNextNodeClasses(ELEMENT * theElement)6535 INT NS_DIM_PREFIX SeedNextNodeClasses (ELEMENT *theElement)
6536 {
6537 INT i;
6538
6539 for (i=0; i<CORNERS_OF_ELEM(theElement); i++)
6540 SETNNCLASS(CORNER(theElement,i),3);
6541
6542 return (0);
6543 }
6544
6545 /****************************************************************************/
6546 /** \brief
6547 PropagateNextNodeClasses - Compute 'NNCLASS' in all Nodes of a grid level
6548
6549 * @param theGrid - pointer to grid
6550
6551 Computes values of 'NNCLASS' field in all nodes after seed.
6552
6553 @return <ul>
6554 <li> 0 if ok </li>
6555 <li> 1 if error occured </li>
6556 </ul> */
6557 /****************************************************************************/
6558
6559 #ifdef ModelP
Gather_NextNodeClass(DDD::DDDContext &,DDD_OBJ obj,void * data)6560 static int Gather_NextNodeClass (DDD::DDDContext&, DDD_OBJ obj, void *data)
6561 {
6562 NODE *theNode = (NODE *)obj;
6563
6564 ((INT *)data)[0] = NNCLASS(theNode);
6565
6566 return(GM_OK);
6567 }
6568
Scatter_NextNodeClass(DDD::DDDContext &,DDD_OBJ obj,void * data)6569 static int Scatter_NextNodeClass (DDD::DDDContext&, DDD_OBJ obj, void *data)
6570 {
6571 NODE *theNode = (NODE *)obj;
6572
6573 SETNNCLASS(theNode,MAX(NNCLASS(theNode),((INT *)data)[0]));
6574
6575 return(GM_OK);
6576 }
6577
Scatter_GhostNextNodeClass(DDD::DDDContext &,DDD_OBJ obj,void * data)6578 static int Scatter_GhostNextNodeClass (DDD::DDDContext&, DDD_OBJ obj, void *data)
6579 {
6580 NODE *theNode = (NODE *)obj;
6581
6582 SETNNCLASS(theNode,((INT *)data)[0]);
6583
6584 return(GM_OK);
6585 }
6586 #endif
6587
PropagateNextNodeClass(GRID * theGrid,INT nnclass)6588 static INT PropagateNextNodeClass (GRID *theGrid, INT nnclass)
6589 {
6590 ELEMENT *theElement;
6591 INT i;
6592
6593 for (theElement=FIRSTELEMENT(theGrid);
6594 theElement!= NULL; theElement = SUCCE(theElement))
6595 if (MaxNextNodeClass(theElement) == nnclass)
6596 for (i=0; i<CORNERS_OF_ELEM(theElement); i++)
6597 {
6598 NODE *theNode = CORNER(theElement,i);
6599
6600 if (NNCLASS(theNode) < nnclass)
6601 SETNNCLASS(theNode,nnclass-1);
6602 }
6603
6604 /* only for this values valid */
6605 ASSERT(nnclass==3 || nnclass==2);
6606
6607 return(0);
6608 }
6609
PropagateNextNodeClasses(GRID * theGrid)6610 INT NS_DIM_PREFIX PropagateNextNodeClasses (GRID *theGrid)
6611 {
6612 #ifdef ModelP
6613 auto& context = theGrid->dddContext();
6614 const auto& dddctrl = ddd_ctrl(context);
6615
6616 PRINTDEBUG(gm,1,("\n" PFMT "PropagateNextNodeClasses(): 1. communication\n",theGrid->ppifContext().me()))
6617 /* exchange NNCLASS of Nodes */
6618 DDD_IFAExchange(context,
6619 dddctrl.BorderNodeSymmIF,GRID_ATTR(theGrid),sizeof(INT),
6620 Gather_NextNodeClass, Scatter_NextNodeClass);
6621 #endif
6622
6623 if (PropagateNextNodeClass(theGrid,3)) REP_ERR_RETURN(1);
6624
6625 #ifdef ModelP
6626 PRINTDEBUG(gm,1,("\n" PFMT "PropagateNextNodeClasses(): 2. communication\n",theGrid->ppifContext().me()))
6627 /* exchange NNCLASS of Nodes */
6628 DDD_IFAExchange(context,
6629 dddctrl.BorderNodeSymmIF,GRID_ATTR(theGrid),sizeof(INT),
6630 Gather_NextNodeClass, Scatter_NextNodeClass);
6631 #endif
6632
6633 if (PropagateNextNodeClass(theGrid,2)) REP_ERR_RETURN(1);
6634
6635 #ifdef ModelP
6636 PRINTDEBUG(gm,1,("\n" PFMT "PropagateNextNodeClasses(): 3. communication\n",theGrid->ppifContext().me()))
6637 /* exchange NNCLASS of Nodes */
6638 DDD_IFAExchange(context,
6639 dddctrl.BorderNodeSymmIF,GRID_ATTR(theGrid),sizeof(INT),
6640 Gather_NextNodeClass, Scatter_NextNodeClass);
6641
6642 /* send NNCLASSn to ghosts */
6643 DDD_IFAOneway(context,
6644 dddctrl.NodeIF,GRID_ATTR(theGrid),IF_FORWARD,sizeof(INT),
6645 Gather_NextNodeClass, Scatter_GhostNextNodeClass);
6646 #endif
6647
6648 return(0);
6649 }
6650
6651 /****************************************************************************/
6652 /** \brief Set subdomain id on level 0 edges
6653
6654 * @param id - the id of the block to be allocated
6655
6656 This function sets the subdomain id taken from the elements for level 0 edges.
6657
6658 @return <ul>
6659 <li> GM_OK if ok </li>
6660 <li> GM_ERROR if error occured </li>
6661 </ul> */
6662 /****************************************************************************/
6663
SetEdgeAndNodeSubdomainFromElements(GRID * theGrid)6664 static INT SetEdgeAndNodeSubdomainFromElements (GRID *theGrid)
6665 {
6666 ELEMENT *theElement;
6667 NODE *n0,*n1;
6668 EDGE *ed;
6669 INT s_id,s,i,k;
6670
6671 /* first set subdomain id for all edges */
6672 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL;
6673 theElement=SUCCE(theElement))
6674 {
6675 /* all edges of the element acquire the subdomain id of the element */
6676 s_id = SUBDOMAIN(theElement);
6677 for (k=0; k<EDGES_OF_ELEM(theElement); k++)
6678 {
6679 n0 = CORNER(theElement,CORNER_OF_EDGE(theElement,k,0));
6680 n1 = CORNER(theElement,CORNER_OF_EDGE(theElement,k,1));
6681 ed = GetEdge(n0,n1);
6682 ASSERT(ed!=NULL);
6683 SETEDSUBDOM(ed,s_id);
6684 }
6685 for (i=0; i<CORNERS_OF_ELEM(theElement); i++)
6686 SETNSUBDOM(CORNER(theElement,i),s_id);
6687 }
6688
6689 /* now change subdomain id for boundary edges and nodes to 0 */
6690 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL;
6691 theElement=SUCCE(theElement))
6692 if (OBJT(theElement)==BEOBJ)
6693 for (s=0; s<SIDES_OF_ELEM(theElement); s++)
6694 {
6695 if (ELEM_BNDS(theElement,s)==NULL)
6696 continue;
6697
6698 for (i=0; i<EDGES_OF_SIDE(theElement,s); i++)
6699 {
6700 k = EDGE_OF_SIDE(theElement,s,i);
6701 n0 = CORNER(theElement,CORNER_OF_EDGE(theElement,k,0));
6702 n1 = CORNER(theElement,CORNER_OF_EDGE(theElement,k,1));
6703 SETNSUBDOM(n0,0);
6704 ASSERT(OBJT(MYVERTEX(n0)) == BVOBJ);
6705 SETNSUBDOM(n1,0);
6706 ASSERT(OBJT(MYVERTEX(n1)) == BVOBJ);
6707 ed = GetEdge(n0,n1);
6708 ASSERT(ed!=NULL);
6709 SETEDSUBDOM(ed,0);
6710 }
6711 }
6712 IFDEBUG(gm,1)
6713 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL;
6714 theElement=SUCCE(theElement))
6715 {
6716 PRINTDEBUG(gm,1,("el(%d)-sd=%d\n",ID(theElement),
6717 SUBDOMAIN(theElement)));
6718 for (k=0; k<EDGES_OF_ELEM(theElement); k++)
6719 {
6720 n0 = CORNER(theElement,CORNER_OF_EDGE(theElement,k,0));
6721 n1 = CORNER(theElement,CORNER_OF_EDGE(theElement,k,1));
6722 ed = GetEdge(n0,n1);
6723 PRINTDEBUG(gm,1,(" ed(%d,%d)-sd=%d nsub %d %d\n",ID(n0),ID(n1),
6724 EDSUBDOM(ed),NSUBDOM(n0),NSUBDOM(n1)));
6725 }
6726 }
6727 ENDDEBUG
6728
6729 return (GM_OK);
6730 }
6731
6732 /****************************************************************************/
6733 /** \brief
6734 RemoveSpuriousBoundarySides - remove boundary side of element and neighbour
6735
6736 * @param elem - element
6737 * @param side - boundary side to remove
6738
6739 This function removes the boundary side of element and neighbour.
6740
6741 @return <ul>
6742 <li> GM_OK if ok </li>
6743 <li> GM_ERROR if error occured </li>
6744 </ul> */
6745 /****************************************************************************/
6746
RemoveSpuriousBoundarySides(HEAP * heap,ELEMENT * elem,INT side)6747 static INT RemoveSpuriousBoundarySides (HEAP *heap, ELEMENT *elem, INT side)
6748 {
6749 ELEMENT *nb=NBELEM(elem,side);
6750 BNDS *nbbside,*bside=ELEM_BNDS(elem,side);
6751 INT nbside;
6752
6753 ASSERT(bside!=NULL);
6754 ASSERT(OBJT(elem)==BEOBJ);
6755 ASSERT(nb!=NULL);
6756 ASSERT(OBJT(nb)==BEOBJ);
6757
6758 /* search nbside */
6759 for (nbside=0; nbside<SIDES_OF_ELEM(nb); nbside++)
6760 if (NBELEM(nb,nbside)==elem)
6761 break;
6762 ASSERT(nbside<SIDES_OF_ELEM(nb));
6763 nbbside = ELEM_BNDS(nb,nbside);
6764 ASSERT(nbbside!=NULL);
6765
6766 PRINTDEBUG(gm,1,("spurious bsides between elem %ld and elem %ld removed",(long)ID(elem),(long)ID(nb)));
6767
6768 if (BNDS_Dispose(heap,bside))
6769 REP_ERR_RETURN(1);
6770 SET_BNDS(elem,side,NULL);
6771
6772 if (BNDS_Dispose(heap,nbbside))
6773 REP_ERR_RETURN(2);
6774 SET_BNDS(nb,nbside,NULL);
6775
6776 return (0);
6777 }
6778
6779 /****************************************************************************/
6780 /** \brief Replace boundary by inner element
6781 *
6782 * @param grid grid where elem resides
6783 * @param elemH handle to element to be replaced: will point to new element
6784 *
6785 * This function replaces a boundary by an inner element.
6786 *
6787 * @return <ul>
6788 * <li> GM_OK if ok </li>
6789 * <li> GM_ERROR if error occured </li>
6790 **</ul>
6791 */
6792 /****************************************************************************/
6793
BElem2IElem(GRID * grid,ELEMENT ** elemH)6794 static INT BElem2IElem (GRID *grid, ELEMENT **elemH)
6795 {
6796 ELEMENT *nb[MAX_SIDES_OF_ELEM],*elem=*elemH,*ielem;
6797 NODE *nodes[MAX_CORNERS_OF_ELEM];
6798 INT i,j,nbside[MAX_SIDES_OF_ELEM],s_id;
6799
6800 ASSERT(GLEVEL(grid)==0);
6801
6802 /* save context */
6803 for (i=0; i<CORNERS_OF_ELEM(elem); i++)
6804 nodes[i] = CORNER(elem,i);
6805
6806 for (i=0; i<SIDES_OF_ELEM(elem); i++)
6807 {
6808 nb[i] = NBELEM(elem,i);
6809 for (j=0; j<SIDES_OF_ELEM(nb[i]); j++)
6810 if (NBELEM(nb[i],j)==elem)
6811 break;
6812 ASSERT(j<SIDES_OF_ELEM(nb[i]));
6813 nbside[i] = j;
6814 }
6815
6816 s_id = SUBDOMAIN(elem);
6817
6818 /* create/dispose */
6819 ielem = CreateElement(grid,TAG(elem),IEOBJ,nodes,EFATHER(elem),NO);
6820 if (ielem==NULL)
6821 REP_ERR_RETURN(1);
6822
6823 if (DisposeElement(grid,elem,NO))
6824 REP_ERR_RETURN(1);
6825
6826 *elemH = ielem;
6827
6828 /* set context */
6829 for (i=0; i<SIDES_OF_ELEM(ielem); i++)
6830 {
6831 SET_NBELEM(ielem,i,nb[i]);
6832
6833 SET_NBELEM(nb[i],nbside[i],ielem);
6834 }
6835 SETSUBDOMAIN(ielem,s_id);
6836 SETECLASS(ielem,RED_CLASS);
6837
6838 return (0);
6839 }
6840
6841 /****************************************************************************/
6842 /** \brief
6843 FinishGrid - remove erroneously introduced bsides and propagate sub domain IDs
6844
6845 * @param mg - multigrid
6846
6847 This function removes erroneously introduced bsides and propagates sub domain IDs.
6848
6849 @return <ul>
6850 <li> GM_OK if ok </li>
6851 <li> GM_ERROR if error occured </li>
6852 </ul> */
6853 /****************************************************************************/
6854
FinishGrid(MULTIGRID * mg)6855 static INT FinishGrid (MULTIGRID *mg)
6856 {
6857 GRID *grid;
6858 ELEMENT *elem,*nb,*succ;
6859 HEAP *heap=MGHEAP(mg);
6860 FIFO unused,shell;
6861 INT MarkKey = MG_MARK_KEY(mg);
6862 INT i,side,id,nbid,part,nsd,found,s_id;
6863 INT *sd_table;
6864 void *buffer;
6865
6866 /* prepare */
6867 if (TOPLEVEL(mg)<0)
6868 REP_ERR_RETURN (GM_ERROR);
6869 grid = GRID_ON_LEVEL(mg,0);
6870 if (!NT(grid))
6871 return (GM_OK);
6872
6873 for (elem=PFIRSTELEMENT(grid); elem!=NULL; elem=SUCCE(elem))
6874 {
6875 SETUSED(elem,false);
6876 SETTHEFLAG(elem,false);
6877 }
6878
6879 /* table for subdomain ids */
6880 nsd = 1 + BVPD_NSUBDOM(MG_BVPD(mg));
6881 sd_table = (INT*)GetTmpMem(heap,nsd*sizeof(INT),MarkKey);
6882 if (sd_table==NULL)
6883 REP_ERR_RETURN (GM_ERROR);
6884
6885 /* init two fifos */
6886 buffer=(void *)GetTmpMem(heap,sizeof(ELEMENT*)*NT(grid),MarkKey);
6887 if (buffer==NULL)
6888 REP_ERR_RETURN (GM_ERROR);
6889 fifo_init(&unused,buffer,sizeof(ELEMENT*)*NT(grid));
6890 buffer=(void *)GetTmpMem(heap,sizeof(ELEMENT*)*NT(grid),MarkKey);
6891 if (buffer==NULL)
6892 REP_ERR_RETURN (GM_ERROR);
6893 fifo_init(&shell,buffer,sizeof(ELEMENT*)*NT(grid));
6894
6895 /* outermost loop handles nonconnected domains */
6896 while (true)
6897 {
6898 for (elem=PFIRSTELEMENT(grid); elem!=NULL; elem=SUCCE(elem))
6899 if (!USED(elem))
6900 break;
6901 if (elem!=NULL)
6902 fifo_in(&unused,elem);
6903 else
6904 break;
6905
6906 while (!fifo_empty(&unused))
6907 {
6908 /* grab next !USED element */
6909 do
6910 elem = (ELEMENT*) fifo_out(&unused);
6911 while (USED(elem) && !fifo_empty(&unused));
6912 if (USED(elem))
6913 /* we are done */
6914 break;
6915
6916 /* shell algo (using FLAG): neighbours, but not across bside */
6917 fifo_clear(&shell);
6918 fifo_in(&shell,elem);
6919 SETTHEFLAG(elem,true);
6920 for (i=0; i<=nsd; i++) sd_table[i] = 0;
6921 found = false;
6922 while (!fifo_empty(&shell))
6923 {
6924 elem = (ELEMENT*) fifo_out(&shell);
6925
6926 if (OBJT(elem)==BEOBJ)
6927 for (side=0; side<SIDES_OF_ELEM(elem); side++)
6928 if (SIDE_ON_BND(elem,side))
6929 {
6930 if (BNDS_BndSDesc(ELEM_BNDS(elem,side),&id,&nbid,&part))
6931 REP_ERR_RETURN (GM_ERROR);
6932
6933 if ((nb=NBELEM(elem,side))==NULL)
6934 {
6935 /* this bside must be ok (outer boundary) */
6936 /* TODO (HRR 971012): parallel? */
6937 ASSERT(nbid==0);
6938 s_id = id;
6939 found = true;
6940 break;
6941 }
6942 else
6943 if (USED(nb))
6944 {
6945 /* he must know! */
6946 if (nbid==SUBDOMAIN(nb))
6947 s_id = id;
6948 else if (id==SUBDOMAIN(nb))
6949 s_id = nbid;
6950 else
6951 ASSERT(false);
6952 }
6953
6954 /* handle outer boundary cases */
6955 if (id==0)
6956 {
6957 ASSERT(nbid>0);
6958 s_id = nbid;
6959 found = true;
6960 break;
6961 }
6962 if (nbid==0)
6963 {
6964 ASSERT(id>0);
6965 s_id = id;
6966 found = true;
6967 break;
6968 }
6969
6970 ++sd_table[id];
6971 if (sd_table[id]>1)
6972 {
6973 s_id = id;
6974 found = true;
6975 break;
6976 }
6977 }
6978 if (found)
6979 break;
6980
6981 /* push neighbours not across boundary */
6982 if (OBJT(elem)==BEOBJ)
6983 {
6984 for (side=0; side<SIDES_OF_ELEM(elem); side++)
6985 if (!SIDE_ON_BND(elem,side))
6986 if ((nb=NBELEM(elem,side))!=NULL)
6987 if (!USED(nb) && !THEFLAG(nb))
6988 {
6989 fifo_in(&shell,nb);
6990 SETTHEFLAG(nb,true);
6991 }
6992 }
6993 else
6994 {
6995 for (side=0; side<SIDES_OF_ELEM(elem); side++)
6996 if ((nb=NBELEM(elem,side))!=NULL)
6997 if (!USED(nb) && !THEFLAG(nb))
6998 {
6999 fifo_in(&shell,nb);
7000 SETTHEFLAG(nb,true);
7001 }
7002 }
7003 }
7004
7005 /* count occurences of subdom ids (max 2 different) */
7006 for (found=0, i=0; i<=nsd; i++)
7007 if (sd_table[i])
7008 found++;
7009 if (found>2)
7010 /* FATAL: algorithm relies on assumptions obviously not fulfilled! */
7011 ASSERT(false);
7012
7013 /* again shell algo starting from last element */
7014 /* set USED, propagate subdomain ids and remove spurious bsides (has to have NB!) */
7015 /* use PRINTDEBUG */
7016 fifo_clear(&shell);
7017 fifo_in(&shell,elem);
7018 SETUSED(elem,true);
7019 SETSUBDOMAIN(elem,s_id);
7020 while (!fifo_empty(&shell))
7021 {
7022 elem = (ELEMENT*) fifo_out(&shell);
7023
7024 if (OBJT(elem)==BEOBJ)
7025 {
7026 for (side=0; side<SIDES_OF_ELEM(elem); side++)
7027 if (SIDE_ON_BND(elem,side))
7028 {
7029 if ((nb=NBELEM(elem,side))==NULL)
7030 continue;
7031 if (!USED(nb))
7032 /* push unused neighbour across boundary to unused fifo */
7033 fifo_in(&unused,nb);
7034
7035 if (BNDS_BndSDesc(ELEM_BNDS(elem,side),&id,&nbid,&part))
7036 REP_ERR_RETURN (GM_ERROR);
7037
7038 if (id!=s_id || nbid==0)
7039 {
7040 /* remove spurious bside of both elements */
7041 if (RemoveSpuriousBoundarySides(heap,elem,side))
7042 REP_ERR_RETURN(1);
7043 }
7044 }
7045 }
7046
7047 /* push neighbours not across boundary */
7048 if (OBJT(elem)==BEOBJ)
7049 {
7050 for (side=0; side<SIDES_OF_ELEM(elem); side++)
7051 if (!SIDE_ON_BND(elem,side))
7052 {
7053 if ((nb=NBELEM(elem,side))!=NULL)
7054 {
7055 if (!USED(nb))
7056 {
7057 fifo_in(&shell,nb);
7058 SETUSED(nb,true);
7059 SETSUBDOMAIN(nb,s_id);
7060 }
7061 }
7062 else
7063 /* TODO (HRR 971012): ModelP: no error if EGHOST? */
7064 /* grid not closed */
7065 REP_ERR_RETURN(1);
7066 }
7067 }
7068 else
7069 {
7070 for (side=0; side<SIDES_OF_ELEM(elem); side++)
7071 if ((nb=NBELEM(elem,side))!=NULL)
7072 {
7073 if (!USED(nb))
7074 {
7075 fifo_in(&shell,nb);
7076 SETUSED(nb,true);
7077 SETSUBDOMAIN(nb,s_id);
7078 }
7079 }
7080 else
7081 /* TODO (HRR 971012): ModelP: no error if EGHOST? */
7082 /* grid not closed */
7083 REP_ERR_RETURN(1);
7084 }
7085 }
7086 }
7087 }
7088
7089 for (elem=PFIRSTELEMENT(grid); elem!=NULL; elem=succ)
7090 {
7091 succ = SUCCE(elem);
7092
7093 if (OBJT(elem)!=BEOBJ) continue;
7094
7095 /* check whether element still has bsides */
7096 for (side=0; side<SIDES_OF_ELEM(elem); side++)
7097 if (ELEM_BNDS(elem,side)!=NULL)
7098 break;
7099 if (side>=SIDES_OF_ELEM(elem))
7100 if (BElem2IElem(grid,&elem))
7101 REP_ERR_RETURN(1);
7102 }
7103
7104 if (SetEdgeAndNodeSubdomainFromElements(grid))
7105 REP_ERR_RETURN (GM_ERROR);
7106
7107 return (GM_OK);
7108 }
7109
7110 /****************************************************************************/
7111 /** \brief
7112 SetSubdomainIDfromBndInfo - set subdomain id on level 0 elements and edges
7113
7114 * @param id - the id of the block to be allocated
7115
7116 This function sets the subdomain for level 0 elements and edges.
7117
7118 @return <ul>
7119 <li> GM_OK if ok </li>
7120 <li> GM_ERROR if error occured </li>
7121 </ul> */
7122 /****************************************************************************/
7123
SetSubdomainIDfromBndInfo(MULTIGRID * theMG)7124 INT NS_DIM_PREFIX SetSubdomainIDfromBndInfo (MULTIGRID *theMG)
7125 {
7126 HEAP *theHeap;
7127 GRID *theGrid;
7128 ELEMENT *theElement, *theNeighbor;
7129 NODE *theNode;
7130 void *buffer;
7131 INT i,n,id,nbid,part,j;
7132 FIFO myfifo;
7133 INT MarkKey = MG_MARK_KEY(theMG);
7134
7135 /* prepare */
7136 if (TOPLEVEL(theMG)<0) REP_ERR_RETURN (GM_ERROR);
7137 theGrid = GRID_ON_LEVEL(theMG,0);
7138 n = NT(theGrid); if (n==0) return(0);
7139
7140 /* allocate fifo and init */
7141 theHeap = MYMG(theGrid)->theHeap;
7142 buffer=(void *)GetTmpMem(theHeap,sizeof(ELEMENT*)*n,MarkKey);
7143 fifo_init(&myfifo,buffer,sizeof(ELEMENT*)*n);
7144 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL;
7145 theElement=SUCCE(theElement))
7146 SETUSED(theElement,0);
7147
7148 /* insert all boundary elements */
7149 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL;
7150 theElement=SUCCE(theElement))
7151 if (OBJT(theElement)==BEOBJ && !USED(theElement))
7152 {
7153 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
7154 if (ELEM_BNDS(theElement,i)!=NULL)
7155 break;
7156 assert(i<SIDES_OF_ELEM(theElement));
7157
7158 /* set id from BNDS */
7159 if (BNDS_BndSDesc(ELEM_BNDS(theElement,i),&id,&nbid,&part))
7160 REP_ERR_RETURN (GM_ERROR);
7161 assert(id>0);
7162 SETSUBDOMAIN(theElement,id);
7163 SETUSED(theElement,1);
7164 fifo_in(&myfifo,(void *)theElement);
7165 PRINTDEBUG(gm,1,("elem %3d sid %d\n",ID(theElement),SUBDOMAIN(theElement)));
7166 for (i=0; i<CORNERS_OF_ELEM(theElement); i++)
7167 {
7168 theNode = CORNER(theElement,i);
7169 if (OBJT(MYVERTEX(theNode))==IVOBJ)
7170 SETNSUBDOM(theNode,id);
7171 }
7172 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
7173 {
7174 if (NBELEM(theElement,i)==NULL || SIDE_ON_BND(theElement,i)) continue;
7175 theNeighbor = NBELEM(theElement,i);
7176 if (USED(theNeighbor))
7177 assert(SUBDOMAIN(theElement)==SUBDOMAIN(theNeighbor));
7178 }
7179 }
7180
7181 /* set subdomain id for all elements */
7182 while(!fifo_empty(&myfifo))
7183 {
7184 theElement = (ELEMENT*)fifo_out(&myfifo);
7185 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
7186 {
7187 if (NBELEM(theElement,i)==NULL) continue;
7188 theNeighbor = NBELEM(theElement,i);
7189 if (USED(theNeighbor))
7190 {
7191 if (INNER_SIDE(theElement,i))
7192 assert(SUBDOMAIN(theElement)==SUBDOMAIN(theNeighbor));
7193 continue;
7194 }
7195 SETSUBDOMAIN(theNeighbor,SUBDOMAIN(theElement));
7196 SETUSED(theNeighbor,1);
7197 for (j=0; j<CORNERS_OF_ELEM(theElement); j++)
7198 {
7199 theNode = CORNER(theElement,j);
7200 if (OBJT(MYVERTEX(theNode))==IVOBJ)
7201 SETNSUBDOM(theNode,SUBDOMAIN(theElement));
7202 }
7203 fifo_in(&myfifo,(void *)theNeighbor);
7204 }
7205 }
7206
7207 IFDEBUG(gm,1)
7208 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL; theElement=SUCCE(theElement))
7209 assert(USED(theElement));
7210 ENDDEBUG
7211
7212 if (SetEdgeAndNodeSubdomainFromElements(theGrid))
7213 REP_ERR_RETURN (GM_ERROR);
7214
7215 return (GM_OK);
7216 }
7217
7218 /****************************************************************************/
7219 /** \brief Do all that is necessary to complete the coarse grid
7220
7221 * @param id - the id of the block to be allocated
7222
7223 This function does all that is necessary to complete the coarse grid.
7224 Finally the MG_COARSE_FIXED flag is set.
7225
7226 @return <ul>
7227 <li> GM_OK if ok </li>
7228 <li> GM_ERROR if error occured </li>
7229 * </ul>
7230 */
7231 /****************************************************************************/
7232
FixCoarseGrid(MULTIGRID * theMG)7233 INT NS_DIM_PREFIX FixCoarseGrid (MULTIGRID *theMG)
7234 {
7235 if (MG_COARSE_FIXED(theMG)) return (GM_OK);
7236
7237 /** \todo (HRR 971031): check that before check-in!
7238 if (FinishGrid(theMG))
7239 REP_ERR_RETURN (GM_ERROR);*/
7240
7241 /** \todo (HRR 971031): remove if above works */
7242 if (SetSubdomainIDfromBndInfo(theMG))
7243 REP_ERR_RETURN (GM_ERROR);
7244
7245 /* set this flag here because it is checked by CreateAlgebra */
7246 if (CreateAlgebra(theMG) != GM_OK)
7247 REP_ERR_RETURN (GM_ERROR);
7248
7249 /* here all temp memory since CreateMultiGrid is released */
7250 ReleaseTmpMem(MGHEAP(theMG),MG_MARK_KEY(theMG));
7251 MG_MARK_KEY(theMG) = 0;
7252
7253 return (GM_OK);
7254 }
7255
7256 /****************************************************************************/
7257 /** \brief Init what is necessary
7258 *
7259 * This function initializes the grid manager.
7260 *
7261 * @return <ul>
7262 * <li> GM_OK if ok </li>
7263 * <li> > 0 line in which error occured </li>
7264 * </ul>
7265 */
7266 /****************************************************************************/
7267
InitUGManager()7268 INT NS_DIM_PREFIX InitUGManager ()
7269 {
7270 INT i;
7271
7272 /* install the /Multigrids directory */
7273 if (ChangeEnvDir("/")==NULL)
7274 {
7275 PrintErrorMessage('F',"InitUGManager","could not changedir to root");
7276 return(__LINE__);
7277 }
7278 theMGRootDirID = GetNewEnvDirID();
7279 if (MakeEnvItem("Multigrids",theMGRootDirID,sizeof(ENVDIR))==NULL)
7280 {
7281 PrintErrorMessage('F',"InitUGManager","could not install /Multigrids dir");
7282 return(__LINE__);
7283 }
7284 theMGDirID = GetNewEnvDirID();
7285
7286 /* init the OBJT management */
7287 UsedOBJT = 0;
7288 for (i=0; i<NPREDEFOBJ; i++)
7289 SET_FLAG(UsedOBJT,1<<i);
7290
7291 return (GM_OK);
7292 }
7293
ExitUGManager()7294 INT NS_DIM_PREFIX ExitUGManager ()
7295 {
7296 return 0;
7297 }
7298
7299 /* nur temporaer zum Debuggen drin (Christian Wrobel): */
7300 /* TODO: entfernen nach Debuggphase */
PrintElementInfo(ELEMENT * theElement,INT full)7301 char *PrintElementInfo (ELEMENT *theElement,INT full)
7302 {
7303 static char out[2000];
7304 char tmp[200];
7305 char etype[10];
7306 char ekind[8];
7307 int i,j;
7308 ELEMENT *SonList[MAX_SONS];
7309
7310 if (theElement==NULL)
7311 {
7312 printf( "PrintElementInfo: element == NULL\n");
7313 return (NULL);
7314 }
7315
7316 if (DIM==2)
7317 switch (TAG(theElement))
7318 {
7319 case TRIANGLE : strcpy(etype,"TRI"); break;
7320 case QUADRILATERAL : strcpy(etype,"QUA"); break;
7321 default : strcpy(etype,"???"); break;
7322 }
7323 else
7324 switch (TAG(theElement))
7325 {
7326 case TETRAHEDRON : strcpy(etype,"TET"); break;
7327 case PYRAMID : strcpy(etype,"PYR"); break;
7328 case PRISM : strcpy(etype,"PRI"); break;
7329 case HEXAHEDRON : strcpy(etype,"HEX"); break;
7330 default : strcpy(etype,"???"); break;
7331 }
7332 switch (ECLASS(theElement))
7333 {
7334 case YELLOW_CLASS : strcpy(ekind,"YELLOW "); break;
7335 case GREEN_CLASS : strcpy(ekind,"GREEN "); break;
7336 case RED_CLASS : strcpy(ekind,"RED "); break;
7337 default : strcpy(ekind,"??? "); break;
7338 }
7339 if(full)
7340 sprintf(out,"ELEMID=" EID_FFMTE " %5s %5s CTRL=%8lx CTRL2=%8lx REFINE=%2d MARK=%2d LEVEL=%2d",
7341 EID_PRTE(theElement),ekind,etype,
7342 (long)CTRL(theElement),(long)FLAG(theElement),REFINE(theElement),MARK(theElement),LEVEL(theElement));
7343 else
7344 sprintf(out,"ELEMID=" EID_FFMTE, EID_PRTE(theElement));
7345
7346 if (COARSEN(theElement)) strcat(out," COARSEN");
7347 strcat(out,"\n");
7348 for (i=0; i<CORNERS_OF_ELEM(theElement); i++)
7349 {
7350 #ifdef __TWODIM__
7351 sprintf(tmp," N%d=" ID_FMTX " x=%g y=%g\n",
7352 i,
7353 ID_PRTX(CORNER(theElement,i)),
7354 CVECT(MYVERTEX(CORNER(theElement,i)))[0],
7355 CVECT(MYVERTEX(CORNER(theElement,i)))[1]
7356 );
7357 #endif
7358 #ifdef __THREEDIM__
7359 sprintf(tmp," N%d=" ID_FMTX " x=%g y=%g z=%g\n",
7360 i,
7361 ID_PRTX(CORNER(theElement,i)),
7362 CVECT(MYVERTEX(CORNER(theElement,i)))[0],
7363 CVECT(MYVERTEX(CORNER(theElement,i)))[1],
7364 CVECT(MYVERTEX(CORNER(theElement,i)))[2]
7365 );
7366 #endif
7367 strcat( out, tmp );
7368 }
7369
7370 if (EFATHER(theElement))
7371 {
7372 sprintf(tmp," FA=" EID_FMTX "\n" ,EID_PRTX(EFATHER(theElement)));
7373 strcat( out, tmp );
7374 }
7375 else
7376 strcat( out," FA=NULL\n");
7377
7378 if( full)
7379 {
7380 UserWriteF(" NSONS=%d\n",NSONS(theElement));
7381 if (GetAllSons(theElement,SonList)==0)
7382 {
7383 for (i=0; SonList[i] != NULL; i++)
7384 {
7385 sprintf(tmp," SON%d " EID_FMTX "\n" ,i,EID_PRTX(SonList[i]));
7386 strcat( out, tmp );
7387
7388 for (j=0; j<CORNERS_OF_ELEM(SonList[i]); j++)
7389 {
7390 #ifdef __TWODIM__
7391 sprintf(tmp," N%d= " ID_FMTX " x=%g y=%g\n",
7392 j,
7393 ID_PRTX(CORNER(SonList[i],j)),
7394 CVECT(MYVERTEX(CORNER(SonList[i],j)))[0],
7395 CVECT(MYVERTEX(CORNER(SonList[i],j)))[1]
7396 );
7397 #endif
7398 #ifdef __THREEDIM__
7399 sprintf(tmp," N%d= " ID_FMTX " x=%g y=%g z=%g\n",
7400 j,
7401 ID_PRTX(CORNER(SonList[i],j)),
7402 CVECT(MYVERTEX(CORNER(SonList[i],j)))[0],
7403 CVECT(MYVERTEX(CORNER(SonList[i],j)))[1],
7404 CVECT(MYVERTEX(CORNER(SonList[i],j)))[2]
7405 );
7406 #endif
7407 strcat( out, tmp );
7408 }
7409 }
7410 }
7411 }
7412 sprintf(tmp," key=%d\n", KeyForObject((KEY_OBJECT *)theElement) );
7413 strcat( out, tmp );
7414
7415 if(full)
7416 {
7417 if (OBJT(theElement)==BEOBJ)
7418 strcat( out," boundary element\n" );
7419 else
7420 strcat( out," no boundary element\n" );
7421
7422 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
7423 {
7424 for(j=0; j<CORNERS_OF_SIDE(theElement,i); j++)
7425 {
7426 #ifdef __TWODIM__
7427 sprintf(tmp," NODE[ID=%ld]: x=%g y=%g",
7428 (long)(ID(CORNER(theElement,CORNER_OF_SIDE(theElement,i,j)))),
7429 CVECT(MYVERTEX(CORNER(theElement,CORNER_OF_SIDE(theElement,i,j))))[0],
7430 CVECT(MYVERTEX(CORNER(theElement,CORNER_OF_SIDE(theElement,i,j))))[1]);
7431 #endif
7432 #ifdef __THREEDIM__
7433 sprintf(tmp," NODE[ID=%ld]: x=%g y=%g z=%g",
7434 (long)(ID(CORNER(theElement,CORNER_OF_SIDE(theElement,i,j)))),
7435 CVECT(MYVERTEX(CORNER(theElement,CORNER_OF_SIDE(theElement,i,j))))[0],
7436 CVECT(MYVERTEX(CORNER(theElement,CORNER_OF_SIDE(theElement,i,j))))[1],
7437 CVECT(MYVERTEX(CORNER(theElement,CORNER_OF_SIDE(theElement,i,j))))[2]);
7438 #endif
7439 strcat( out, tmp );
7440 }
7441 strcat( out,"\n");
7442 }
7443 }
7444 #ifdef ModelP
7445 /*UserWriteF(PFMT"%s", me,out );*/
7446 printf("%s",out);
7447 #else
7448 UserWrite(out);
7449 #endif
7450
7451 return out;
7452 }
7453