1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 /*! \file algebra.c
4 * \ingroup gm
5 */
6
7 /** \addtogroup gm
8 *
9 * @{
10 */
11
12 /****************************************************************************/
13 /* */
14 /* File: algebra.c */
15 /* */
16 /* Purpose: management for algebraic structures */
17 /* */
18 /* Author: Klaus Johannsen */
19 /* Interdisziplinaeres Zentrum fuer Wissenschaftliches Rechnen */
20 /* Universitaet Heidelberg */
21 /* Im Neuenheimer Feld 294 */
22 /* 6900 Heidelberg */
23 /* email: ug@ica3.uni-stuttgart.de */
24 /* */
25 /* blockvector data structure: */
26 /* Christian Wrobel */
27 /* Institut fuer Computeranwendungen III */
28 /* Universitaet Stuttgart */
29 /* Pfaffenwaldring 27 */
30 /* 70569 Stuttgart */
31 /* email: ug@ica3.uni-stuttgart.de */
32 /* */
33 /* History: 1.12.93 begin, ug 3d */
34 /* 26.10.94 begin combination 2D/3D version */
35 /* 28.09.95 blockvector implemented (Christian Wrobel) */
36 /* */
37 /* Remarks: */
38 /* */
39 /****************************************************************************/
40
41
42 /****************************************************************************/
43 /* */
44 /* defines to exclude functions */
45 /* */
46 /****************************************************************************/
47
48
49 /****************************************************************************/
50 /* */
51 /* include files */
52 /* system include files */
53 /* application include files */
54 /* */
55 /****************************************************************************/
56
57 #include <config.h>
58 #include <cassert>
59 #include <cmath>
60 #include <cstdio>
61 #include <cstdlib>
62 #include <cstring>
63
64 #include <dune/uggrid/low/architecture.h>
65 #include <dune/uggrid/low/debug.h>
66 #include <dune/uggrid/low/fifo.h>
67 #include <dune/uggrid/low/heaps.h>
68 #include <dune/uggrid/low/misc.h>
69 #include <dune/uggrid/low/namespace.h>
70 #include <dune/uggrid/low/ugenv.h>
71 #include <dune/uggrid/low/ugtypes.h>
72
73 #include <dune/uggrid/ugdevices.h>
74
75 #include "algebra.h"
76 #include "cw.h"
77 #include "dlmgr.h"
78 #include "gm.h"
79 #include "elements.h"
80 #include "refine.h"
81 #include "ugm.h"
82 #include "evm.h"
83 #include "dlmgr.h"
84 #include "mgheapmgr.h"
85
86 #ifdef ModelP
87 #include <dune/uggrid/parallel/dddif/parallel.h>
88 #endif
89
90
91 USING_UG_NAMESPACE
92 USING_UGDIM_NAMESPACE
93 using namespace PPIF;
94
95
96 /****************************************************************************/
97 /* */
98 /* defines in the following order */
99 /* */
100 /* compile time constants defining static data size (i.e. arrays) */
101 /* other constants */
102 /* macros */
103 /* */
104 /****************************************************************************/
105
106 /** \brief Resolution for LexAlgDep */
107 #define ORDERRES 1e-3
108
109 /** \brief For GetDomainPart indicating an element is meant rather than an element side */
110 #define NOSIDE -1
111
112 /** \brief Constants for the direction of domain halfening */
113 #define BV_VERTICAL 0
114 #define BV_HORIZONTAL 1
115
116 /****************************************************************************/
117 /* */
118 /* data structures used in this source file (exported data structures are */
119 /* in the corresponding include file!) */
120 /* */
121 /****************************************************************************/
122
123 /****************************************************************************/
124 /* */
125 /* definition of exported global variables */
126 /* */
127 /****************************************************************************/
128
129 const char* NS_DIM_PREFIX ObjTypeName[MAXVOBJECTS];
130
131 /****************************************************************************/
132 /* */
133 /* definition of variables global to this source file only (static!) */
134 /* */
135 /****************************************************************************/
136
137 static INT theAlgDepDirID; /* env type for Alg Dep dir */
138 static INT theAlgDepVarID; /* env type for ALG_DEP vars */
139
140 static INT theFindCutDirID; /* env type for FindCut dir */
141 static INT theFindCutVarID; /* env type for FIND_CUT vars */
142
143 #ifdef __TWODIM__
144 static MULTIGRID *GBNV_mg; /* multigrid */
145 static INT GBNV_MarkKey; /* key for Mark/Release */
146 #endif
147 static INT GBNV_n; /* list items */
148 static INT GBNV_curr; /* curr pos */
149 static VECTOR **GBNV_list=NULL; /* list pointer */
150
151 /****************************************************************************/
152 /* */
153 /* definition of variables global to this source file only (static!) */
154 /* */
155 /****************************************************************************/
156
157 /* for LexOrderVectorsInGrid */
158 static DOUBLE InvMeshSize;
159
160 /****************************************************************************/
161 /** \brief Compute part information of geometrical object
162 *
163 * @param s2p - table translating subdomain to domain part
164 * @param obj - geometric object (node, element or edge)
165 * @param side - if element side is meant for obj==element side has to be >=0, negative else
166
167 Compute the part info for a geometrical object including element sides.
168
169 * @return <ul>
170 * <li> part if ok </li>
171 * <li> -n else </li>
172 </ul>
173 */
174 /****************************************************************************/
175
GetDomainPart(const INT s2p[],const GEOM_OBJECT * obj,INT side)176 INT NS_DIM_PREFIX GetDomainPart (const INT s2p[], const GEOM_OBJECT *obj, INT side)
177 {
178 NODE *nd,*n0,*n1;
179 EDGE *ed;
180 ELEMENT *elem;
181 VERTEX *v0,*v1;
182 BNDS *bs;
183 INT part=-1,subdom,move,left,right;
184
185 switch (OBJT(obj))
186 {
187 case NDOBJ :
188 nd = (NODE*)obj;
189 v0 = MYVERTEX(nd);
190 if (OBJT(v0)==IVOBJ) {
191 subdom = NSUBDOM(nd);
192 ASSERT(subdom>0);
193 part = s2p[subdom];
194 }
195 else
196 {
197 /* get part info from domain module */
198 if (BNDP_BndPDesc(V_BNDP(v0),&move,&part))
199 REP_ERR_RETURN(-2);
200 ASSERT(NSUBDOM(nd) == 0);
201 }
202 break;
203
204 case IEOBJ :
205 case BEOBJ :
206 elem = (ELEMENT*)obj;
207 if (side==NOSIDE)
208 {
209 /* get info for element */
210
211 subdom = SUBDOMAIN(elem);
212 ASSERT(subdom>0);
213 part = s2p[subdom];
214 }
215 else
216 {
217 /* get info for element side */
218
219 ASSERT(side<SIDES_OF_ELEM(elem));
220 ASSERT(side>=0);
221
222 if ((OBJT(elem)==BEOBJ) && ((bs = ELEM_BNDS(elem,side))!=NULL))
223 {
224 /* this is a boundary side: ask domain module */
225 if (BNDS_BndSDesc(ELEM_BNDS(elem,side),&left,&right,&part))
226 REP_ERR_RETURN(-3);
227 }
228 else
229 {
230 /* there still is the possibility that the side is a boundary side
231 because VECTORs are created while CreateElement but boundary sides later
232 by CreateSonElementSide.
233 The vector eventually will be reallocated by ReinspectSonSideVector later */
234 subdom = SUBDOMAIN(elem);
235
236 // The following assertion is out-commented, for the following reason:
237 // It fails in test-ug, which is likely to indicate a bug that I don't understand.
238 // However, Dune UGGrid does not use the 'part' information anyway,
239 // and therefore I am not motivated to actually go and find the bug.
240 // ASSERT(subdom>0);
241 part = s2p[subdom];
242 }
243 }
244 break;
245
246 case EDOBJ :
247 ed = (EDGE*)obj;
248 n0 = NBNODE(LINK0(ed));
249 n1 = NBNODE(LINK1(ed));
250 v0 = MYVERTEX(n0);
251 v1 = MYVERTEX(n1);
252 if ((OBJT(v0)==BVOBJ) && (OBJT(v1)==BVOBJ))
253 if (BNDP_BndEDesc(V_BNDP(v0),V_BNDP(v1),&part) == 0)
254 return(part);
255 subdom = EDSUBDOM(ed);
256 if (subdom > 0)
257 return(s2p[subdom]);
258 subdom = NSUBDOM(n0);
259 if (subdom > 0)
260 return(s2p[subdom]);
261 subdom = NSUBDOM(n1);
262 if (subdom > 0)
263 return(s2p[subdom]);
264 REP_ERR_RETURN(-4);
265
266 default : REP_ERR_RETURN (-5);
267 }
268 return (part);
269 }
270
271 #ifdef ModelP
GetVectorSize(GRID * theGrid,INT VectorObjType,GEOM_OBJECT * object)272 INT NS_DIM_PREFIX GetVectorSize (GRID *theGrid, INT VectorObjType, GEOM_OBJECT *object)
273 {
274 MULTIGRID *mg;
275 INT part,vtype;
276
277 mg = MYMG(theGrid);
278 part = GetDomainPart(BVPD_S2P_PTR(MG_BVPD(mg)),object,NOSIDE);
279 if (part < 0)
280 REP_ERR_RETURN(-1);
281 vtype = FMT_PO2T(MGFORMAT(mg),part,VectorObjType);
282
283 return(FMT_S_VEC_TP(MGFORMAT(mg),vtype));
284 }
285 #endif
286
287 /****************************************************************************/
288 /** \brief Return pointer to a new vector structure
289 *
290 * @param theGrid - grid where vector should be inserted
291 * @param DomPart - part of the domain where vector is created
292 * @param ObjType - one of the types defined in gm
293 * @param object - associate vector with this object
294 * @param vHandle - handle of new vector, i.e. a pointer to a pointer where
295 a pointer to the new vector is placed.
296
297 This function returns a pointer to a new vector structure.
298 The vector type is determined by DomPart and ObjType
299 First the free list is checked for a free entry, if none
300 is available, a new structure is allocated from the heap.
301
302 * @return <ul>
303 * <li> 0 if ok </li>
304 * <li> 1 if error occured. </li>
305 </ul>
306 */
307 /****************************************************************************/
308
CreateVectorInPart(GRID * theGrid,INT DomPart,INT VectorObjType,GEOM_OBJECT * object,VECTOR ** vHandle)309 static INT CreateVectorInPart (GRID *theGrid, INT DomPart, INT VectorObjType,
310 GEOM_OBJECT *object, VECTOR **vHandle)
311 {
312 MULTIGRID *theMG;
313 VECTOR *pv;
314 INT ds, Size, vtype;
315
316 *vHandle = NULL;
317
318 theMG = MYMG(theGrid);
319 vtype = FMT_PO2T(MGFORMAT(theMG),DomPart,VectorObjType);
320 ds = FMT_S_VEC_TP(MGFORMAT(theMG),vtype);
321 if (ds == 0)
322 return (0); /* HRR: this is ok now, no XXXXVEC in part of the domain */
323
324 Size = sizeof(VECTOR)-sizeof(DOUBLE)+ds;
325 pv = (VECTOR *)GetMemoryForObject(theMG,Size,VEOBJ);
326 if (pv==NULL)
327 REP_ERR_RETURN(1);
328
329 /* initialize data */
330 SETOBJT(pv,VEOBJ);
331 SETVTYPE(pv,vtype);
332 SETVPART(pv,DomPart);
333 ds = VPART(pv);
334 if (ds!=DomPart)
335 return (1);
336 SETVDATATYPE(pv,BITWISE_TYPE(vtype));
337 SETVOTYPE(pv,VectorObjType);
338 SETVCLASS(pv,3);
339 SETVNCLASS(pv,0);
340 SETVBUILDCON(pv,1);
341 SETVNEW(pv,1);
342 /* SETPRIO(dddContext, pv,PrioMaster); */
343
344 #ifndef ModelP
345 // Dune uses the id field for face indices in sequential grids
346 pv->id = (theGrid->mg->vectorIdCounter)++;
347 #endif
348
349 #ifdef ModelP
350 DDD_AttrSet(PARHDR(pv),GRID_ATTR(theGrid));
351 #endif
352
353 VOBJECT(pv) = object;
354 VINDEX(pv) = (long)NVEC(theGrid);
355 SUCCVC(pv) = FIRSTVECTOR(theGrid);
356 VSTART(pv) = NULL;
357
358 GRID_LINK_VECTOR(theGrid,pv,PrioMaster);
359
360 *vHandle = pv;
361
362 PRINTDEBUG(gm,1,("%s-vector created (%d): p=%d, t=%d\n",ObjTypeName[VOTYPE(pv)],ID(VOBJECT(pv)),VPART(pv),VTYPE(pv)));
363
364 return (0);
365 }
366
367 /****************************************************************************/
368 /** \brief Return pointer to a new vector structure
369 *
370 * @param theGrid - grid where vector should be inserted
371 * @param VectorObjType - one of the types defined in gm
372 * @param object - associate vector with this object
373 * @param vHandle - handle of new vector, i.e. a pointer to a pointer where
374 a pointer to the new vector is placed.
375
376 This function returns a pointer to a new vector structure. First the free list
377 is checked for a free entry, if none is available, a new structure is allocated
378 from the heap.
379
380 * @return <ul>
381 * <li> 0 if ok </li>
382 * <li> 1 if error occured. </li>
383 </ul>
384 */
385 /****************************************************************************/
386
CreateVector(GRID * theGrid,INT VectorObjType,GEOM_OBJECT * object,VECTOR ** vHandle)387 INT NS_DIM_PREFIX CreateVector (GRID *theGrid, INT VectorObjType, GEOM_OBJECT *object, VECTOR **vHandle)
388 {
389 MULTIGRID *mg;
390 INT part;
391
392 *vHandle = NULL;
393 mg = MYMG(theGrid);
394
395 part = GetDomainPart(BVPD_S2P_PTR(MG_BVPD(mg)),object,NOSIDE);
396
397 if (part < 0)
398 REP_ERR_RETURN(1);
399 if (CreateVectorInPart(theGrid,part,VectorObjType,object,vHandle)) {
400 REP_ERR_RETURN(1);
401 }
402 else
403 return (0);
404 }
405
CreateSideVector(GRID * theGrid,INT side,GEOM_OBJECT * object,VECTOR ** vHandle)406 INT NS_DIM_PREFIX CreateSideVector (GRID *theGrid, INT side, GEOM_OBJECT *object, VECTOR **vHandle)
407 {
408 MULTIGRID *mg;
409 INT part;
410
411 *vHandle = NULL;
412 mg = MYMG(theGrid);
413 part = GetDomainPart(BVPD_S2P_PTR(MG_BVPD(mg)),object,side);
414 if (part<0)
415 REP_ERR_RETURN(1);
416
417 if (CreateVectorInPart(theGrid,part,SIDEVEC,object,vHandle))
418 REP_ERR_RETURN(1);
419
420 SETVECTORSIDE(*vHandle,side);
421 SETVCOUNT(*vHandle,1);
422
423 return (0);
424 }
425
426 /****************************************************************************/
427 /** \brief Return pointer to a new connection structure
428 *
429 * @param theGrid - grid where matrix should be inserted
430 * @param from - source vector
431 * @param to - destination vector
432
433 This function allocates a new CONNECTION and inserts the two
434 MATRIX structures in the lists of from and to vector.
435 Since the operation is symmetric, the order of from and to
436 is not important.
437
438 * @return <ul>
439 * <li> NULL if error occured. </li>
440 * <li> else a pointer to the new CONNECTION is returned. </li>
441 * </ul>
442 */
443 /****************************************************************************/
444
CreateConnection(GRID * theGrid,VECTOR * from,VECTOR * to)445 CONNECTION * NS_DIM_PREFIX CreateConnection (GRID *theGrid, VECTOR *from, VECTOR *to)
446 {
447 MULTIGRID *theMG;
448 CONNECTION *pc;
449 MATRIX *pm;
450 INT RootType, DestType, MType, ds, Diag, Size;
451
452 /* set Diag, RootType and DestType */
453 Diag = ((from == to) ? 1 : 0);
454 RootType = VTYPE(from);
455 DestType = VTYPE(to);
456 if (Diag)
457 MType=DIAGMATRIXTYPE(RootType);
458 else
459 MType=MATRIXTYPE(RootType,DestType);
460
461 /* check expected size */
462 theMG = MYMG(theGrid);
463 ds = FMT_S_MAT_TP(MGFORMAT(theMG),MType);
464 if (ds == 0)
465 return (NULL);
466 Size = sizeof(MATRIX)-sizeof(DOUBLE)+ds;
467 if (MSIZEMAX<Size) return (NULL);
468
469 /* is there already the desired connection ? */
470 pc = GetConnection(from,to);
471 if (pc != NULL)
472 {
473 SETCEXTRA(pc,0);
474 return (pc);
475 }
476
477 if (Diag)
478 pc = (CONNECTION*)GetMemoryForObject(theMG,Size,MAOBJ);
479 else
480 pc = (CONNECTION*)GetMemoryForObject(theMG,2*Size,COOBJ);
481 if (pc==NULL) return (NULL);
482
483 /* initialize data */
484 pm = CMATRIX0(pc);
485 SETOBJT(pm,MAOBJ);
486 SETMROOTTYPE(pm,RootType);
487 SETMDESTTYPE(pm,DestType);
488 SETMDIAG(pm,Diag);
489 SETMOFFSET(pm,0);
490 SETMSIZE(pm,Size);
491 SETMNEW(pm,1);
492 SETCEXTRA(pc,0);
493 MDEST(pm) = to;
494 if (!Diag)
495 {
496 pm = CMATRIX1(pc);
497 CTRL(pm) = 0;
498 SETOBJT(pm,MAOBJ);
499 SETMROOTTYPE(pm,DestType);
500 SETMDESTTYPE(pm,RootType);
501 SETMDIAG(pm,Diag);
502 SETMOFFSET(pm,1);
503 SETMSIZE(pm,Size);
504 SETMNEW(pm,1);
505 MDEST(pm) = from;
506 }
507
508 /* set sizes */
509 if (!Diag)
510 {
511 Size = (char*)pm - (char*)pc;
512 SETMSIZE(pc,Size);
513 SETMSIZE(pm,Size);
514 }
515
516 /* put in matrix list */
517 if (Diag)
518 {
519 /* insert at first place in the list (only one matrix) */
520 MNEXT(CMATRIX0(pc)) = VSTART(from);
521 VSTART(from) = CMATRIX0(pc);
522 }
523 else
524 {
525 /* insert at second place in the list (both matrices) */
526 pm = VSTART(from);
527 if (pm == NULL)
528 {
529 MNEXT(CMATRIX0(pc)) = NULL;
530 VSTART(from) = CMATRIX0(pc);
531 }
532 else
533 {
534 MNEXT(CMATRIX0(pc)) = MNEXT(pm);
535 MNEXT(pm) = CMATRIX0(pc);
536 }
537
538 pm = VSTART(to);
539 if (pm == NULL)
540 {
541 MNEXT(CMATRIX1(pc)) = NULL;
542 VSTART(to) = CMATRIX1(pc);
543 }
544 else
545 {
546 MNEXT(CMATRIX1(pc)) = MNEXT(pm);
547 MNEXT(pm) = CMATRIX1(pc);
548 }
549 }
550
551 /* counters */
552 theGrid->nCon++;
553
554 return(pc);
555 }
556
557 /****************************************************************************/
558 /** \brief Return pointer to a new matrix structure with extra flag set.
559 *
560 * @param theGrid - grid level where connection will be inserted.
561 * @param from,to - Pointers to vectors where connection is inserted.
562
563 This function returns a pointer to a new CONNECTION
564 structure with extra flag set. This e.g. for a direct solver
565 or ILU with fill in. The new connections can be distinguished
566 from the connections necessary for the stiffness matrix.
567
568 * @return <ul>
569 * <li> NULL if error occured. </li>
570 * <li> else pointer to new CONNECTION </li>
571 </ul>
572 */
573 /****************************************************************************/
574
CreateExtraConnection(GRID * theGrid,VECTOR * from,VECTOR * to)575 CONNECTION *NS_DIM_PREFIX CreateExtraConnection (GRID *theGrid, VECTOR *from, VECTOR *to)
576 {
577 CONNECTION *pc;
578
579 pc = CreateConnection(theGrid,from,to);
580 if (pc==NULL) return(NULL);
581 SETCEXTRA(pc,1);
582 return(pc);
583 }
584
CreateElementList(GRID * theGrid,NODE * theNode,ELEMENT * theElement)585 INT NS_DIM_PREFIX CreateElementList (GRID *theGrid, NODE *theNode, ELEMENT *theElement)
586 {
587 ELEMENTLIST *pel;
588
589 for (pel=NODE_ELEMENT_LIST(theNode); pel!=NULL; pel=NEXT(pel))
590 if(pel->el==theElement)
591 return(0);
592
593 pel = (ELEMENTLIST *)GetMemoryForObject(theGrid->mg,
594 sizeof(ELEMENTLIST),MAOBJ);
595 if (pel == NULL)
596 return(1);
597
598 pel->next = NODE_ELEMENT_LIST(theNode);
599 pel->el = theElement;
600 NDATA(theNode) = (void *) pel;
601
602 return(0);
603 }
604
605 /****************************************************************************/
606 /** \brief Remove vector from the data structure
607 *
608 * @param theGrid - grid level where theVector is in.
609 * @param theVector - VECTOR to be disposed.
610
611 This function removes vector from the data structure and places
612 it in the free list.
613
614 * @return <ul>
615 * <li> 0 if ok </li>
616 * <li> 1 if error occured. </li>
617 </ul>
618 */
619 /****************************************************************************/
620
DisposeVector(GRID * theGrid,VECTOR * theVector)621 INT NS_DIM_PREFIX DisposeVector (GRID *theGrid, VECTOR *theVector)
622 {
623 MATRIX *theMatrix, *next;
624 INT Size;
625
626 if (theVector == NULL)
627 return(0);
628
629 /* remove all connections concerning the vector */
630 for (theMatrix=VSTART(theVector); theMatrix!=NULL; theMatrix=next)
631 {
632 next = MNEXT(theMatrix);
633 if (DisposeConnection(theGrid,MMYCON(theMatrix)))
634 RETURN (1);
635 }
636
637 /* now remove vector from vector list */
638 GRID_UNLINK_VECTOR(theGrid,theVector);
639
640 /* reset count flags */
641 SETVCOUNT(theVector,0);
642
643
644 /* delete the vector itself */
645 Size = sizeof(VECTOR)-sizeof(DOUBLE)
646 + FMT_S_VEC_TP(MGFORMAT(MYMG(theGrid)),VTYPE(theVector));
647 if (PutFreeObject(theGrid->mg,theVector,Size,VEOBJ))
648 RETURN(1);
649
650 return(0);
651 }
652
653 /****************************************************************************/
654 /* \brief Change vector allocated with wrong part
655
656 * @param g - grid level where the element is in.
657 * @param elem - element of side vector
658 * @param side - element side
659 * @param vHandle - handle to side vector (inialized with old, may be changed)
660
661 This changes a side vector which was allocated with the wrong part (maybe unknown for side vectors
662 when creating an element).
663
664 * @return <ul>
665 * <li> GM_OK if ok
666 * <li> GM_ERROR if error occured.
667 */
668 /****************************************************************************/
669
ReinspectSonSideVector(GRID * g,ELEMENT * elem,INT side,VECTOR ** vHandle)670 INT NS_DIM_PREFIX ReinspectSonSideVector (GRID *g, ELEMENT *elem, INT side, VECTOR **vHandle)
671 {
672 MULTIGRID *mg;
673 VECTOR *vold,*vnew;
674 INT partnew,partold,vtnew,vtold,dsnew,dsold;
675
676 mg = MYMG(g);
677 const FORMAT* fmt = mg->theFormat.get();
678
679 vold = *vHandle;
680
681 /* check whether part has actually changed */
682 partold = (vold!=NULL) ? VPART(vold) : BVPD_S2P(MG_BVPD(mg),SUBDOMAIN(elem));
683 partnew = GetDomainPart(BVPD_S2P_PTR(MG_BVPD(mg)),(GEOM_OBJECT*)elem,side);
684 if (partnew<0)
685 REP_ERR_RETURN(GM_ERROR);
686 if (partnew==partold)
687 return (GM_OK);
688
689 /* check whether vtype has actually changed */
690 vtold = (vold!=NULL) ? VTYPE(vold) : FMT_PO2T(fmt,partold,SIDEVEC);
691 vtnew = FMT_PO2T(fmt,partnew,SIDEVEC);
692 if (vtnew==vtold)
693 {
694 if (vold!=NULL)
695 {
696 /* just change part */
697 SETVPART(vold,partnew);
698 }
699 PRINTDEBUG(gm,1,("SIDEVEC (%d,%d): part\n",ID(elem),side));
700 return (GM_OK);
701 }
702
703 /* check whether size has actually changed */
704 dsold = FMT_S_VEC_TP(fmt,vtold);
705 dsnew = FMT_S_VEC_TP(fmt,vtnew);
706 if (dsold==dsnew)
707 {
708 if (vold!=NULL)
709 {
710 /* just change part and type */
711 SETVTYPE(vold,vtnew);
712 SETVPART(vold,partnew);
713
714 DisposeConnectionFromVector(g,vold);
715
716 SETVBUILDCON(vold,1);
717 }
718 PRINTDEBUG(gm,1,("SIDEVEC (%d,%d): part and type\n",ID(elem),side));
719 return (GM_OK);
720 }
721
722 PRINTDEBUG(gm,1,("SIDEVEC (%d,%d): part, type and size\n",ID(elem),side));
723
724 /* create new vector */
725 if (CreateVectorInPart(g,partnew,SIDEVEC,(GEOM_OBJECT*)elem,&vnew))
726 REP_ERR_RETURN(GM_ERROR);
727
728 /* now we can dispose the old vector */
729 if (DisposeVector(g,vold))
730 REP_ERR_RETURN(GM_ERROR);
731
732 *vHandle = vnew;
733
734 return (GM_OK);
735 }
736
737 /****************************************************************************/
738 /** \brief Remove connection from the data structure
739 *
740 * @param theGrid - the grid to remove from
741 * @param theConnection - connection to dispose
742
743 This function removes a connection from the data structure. The connection
744 is removed from the list of the two vectors and is placed in the
745 free list.
746
747 * @return <ul>
748 * <li> 0 if ok </li>
749 * <li> 1 if error occured. </li>
750 </ul>
751 */
752 /****************************************************************************/
753
DisposeConnection(GRID * theGrid,CONNECTION * theConnection)754 INT NS_DIM_PREFIX DisposeConnection (GRID *theGrid, CONNECTION *theConnection)
755 {
756 VECTOR *from, *to;
757 MATRIX *Matrix, *ReverseMatrix, *SearchMatrix;
758
759 /* remove matrix(s) from their list(s) */
760 Matrix = CMATRIX0(theConnection);
761 to = MDEST(Matrix);
762 if (MDIAG(Matrix))
763 {
764 from = to;
765 VSTART(to) = MNEXT(Matrix);
766 }
767 else
768 {
769 ReverseMatrix = CMATRIX1(theConnection);
770 from = MDEST(ReverseMatrix);
771 if (VSTART(from) == Matrix)
772 VSTART(from) = MNEXT(Matrix);
773 else
774 for (SearchMatrix=VSTART(from); SearchMatrix!=NULL; SearchMatrix=MNEXT(SearchMatrix))
775 if (MNEXT(SearchMatrix) == Matrix)
776 MNEXT(SearchMatrix) = MNEXT(Matrix);
777 if (VSTART(to) == ReverseMatrix)
778 VSTART(to) = MNEXT(ReverseMatrix);
779 else
780 for (SearchMatrix=VSTART(to); SearchMatrix!=NULL; SearchMatrix=MNEXT(SearchMatrix))
781 if (MNEXT(SearchMatrix) == ReverseMatrix)
782 MNEXT(SearchMatrix) = MNEXT(ReverseMatrix);
783 }
784
785 /* free connection object */
786 if (MDIAG(Matrix))
787 PutFreeObject(MYMG(theGrid),Matrix,UG_MSIZE(Matrix),MAOBJ);
788 else
789 PutFreeObject(MYMG(theGrid),Matrix,2*UG_MSIZE(Matrix),COOBJ);
790
791 theGrid->nCon--;
792
793 /* return ok */
794 return(0);
795 }
796
797 /****************************************************************************/
798 /** \brief Dispose one of two vectors associated to a face
799 *
800 * @param theGrid - pointer to grid
801 * @param Elem0,side0 - first element and side
802 * @param Elem1,side1 - second element and side
803
804 After grid refinement it may happen that two side vector are associated to a face
805 in the grid. The elements on both sides of the face each have their own side
806 vector. This is an inconsistent state -- there should be only one vector per face.
807 This routine deletes one of the two side vectors and makes the data structure
808 consistent again. This is easier than only creating one unique side vector
809 to begin with.
810
811 * @return <ul>
812 * <li> 0 if ok </li>
813 * <li> 1 if error occured. </li>
814 </ul>
815 */
816 /****************************************************************************/
817
818 #ifdef __THREEDIM__
DisposeDoubledSideVector(GRID * theGrid,ELEMENT * Elem0,INT Side0,ELEMENT * Elem1,INT Side1)819 INT NS_DIM_PREFIX DisposeDoubledSideVector (GRID *theGrid, ELEMENT *Elem0, INT Side0, ELEMENT *Elem1, INT Side1)
820 {
821 VECTOR *Vector0, *Vector1;
822
823 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,SIDEVEC))
824 {
825 assert(NBELEM(Elem0,Side0)==Elem1 && NBELEM(Elem1,Side1)==Elem0);
826 Vector0 = SVECTOR(Elem0,Side0);
827 Vector1 = SVECTOR(Elem1,Side1);
828 if (Vector0 == Vector1)
829 return (0);
830 if ((Vector0==NULL) || (Vector1==NULL))
831 /* this is the case at boundaries between different parts
832 the part not using vectors in SIDEs will not need a pointer
833 to the side vector */
834 return (0);
835 assert(VCOUNT(Vector0)==1 && VCOUNT(Vector1)==1);
836 assert(VSTART(Vector0)==NULL || VSTART(Vector1)==NULL);
837 if (VSTART(Vector0)==NULL)
838 {
839 SET_SVECTOR(Elem0,Side0,Vector1);
840 SETVCOUNT(Vector1,2);
841 if (DisposeVector (theGrid,Vector0))
842 RETURN (1);
843 }
844 else
845 {
846 SET_SVECTOR(Elem1,Side1,Vector0);
847 SETVCOUNT(Vector0,2);
848 if (DisposeVector (theGrid,Vector1))
849 RETURN (1);
850 }
851 return (0);
852
853 }
854
855 RETURN (1);
856 }
857 #endif
858
859 /****************************************************************************/
860 /** \brief Remove all connections associated with a vector
861 *
862 * @param theGrid - grid level where vector belongs to
863 * @param theVector - vector where connections are disposed from
864
865 This function removes all connections from a vector.
866
867 * @return <ul>
868 * <li> 0 if ok </li>
869 * <li> 1 if error occured. </li>
870 </ul>
871 */
872 /****************************************************************************/
873
DisposeConnectionFromVector(GRID * theGrid,VECTOR * theVector)874 INT NS_DIM_PREFIX DisposeConnectionFromVector (GRID *theGrid, VECTOR *theVector)
875 {
876 while(VSTART(theVector) != NULL)
877 if (DisposeConnection (theGrid,MMYCON(VSTART(theVector))))
878 return (1);
879
880 return (0);
881 }
882
883
884 /****************************************************************************/
885 /** \brief Removes all connections from all vectors associated with an element
886 *
887 * @param theGrid - grid level where element is on
888 * @param theElement - element from which to dispose connections
889
890 This function removes all connections from all vectors
891 associated with an element.
892
893 * @return <ul>
894 * <li> 0 if ok </li>
895 * <li> 1 if error occured. </li>
896 </ul>
897 */
898 /****************************************************************************/
899
DisposeConnectionFromElement(GRID * theGrid,ELEMENT * theElement)900 INT NS_DIM_PREFIX DisposeConnectionFromElement (GRID *theGrid, ELEMENT *theElement)
901 {
902 INT i;
903 VECTOR *vList[20];
904 INT cnt;
905
906 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,ELEMVEC))
907 {
908 GetVectorsOfElement(theElement,&cnt,vList);
909 for (i=0; i<cnt; i++)
910 {
911 if (DisposeConnectionFromVector(theGrid,vList[i])) RETURN(GM_ERROR);
912 SETVBUILDCON(vList[i],1);
913 }
914 }
915 #ifdef __THREEDIM__
916 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,SIDEVEC))
917 {
918 GetVectorsOfSides(theElement,&cnt,vList);
919 for (i=0; i<cnt; i++)
920 {
921 if (DisposeConnectionFromVector(theGrid,vList[i])) RETURN(GM_ERROR);
922 SETVBUILDCON(vList[i],1);
923 }
924 }
925 #endif
926 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,EDGEVEC))
927 {
928 GetVectorsOfEdges(theElement,&cnt,vList);
929 for (i=0; i<cnt; i++)
930 {
931 if (DisposeConnectionFromVector(theGrid,vList[i])) RETURN(GM_ERROR);
932 SETVBUILDCON(vList[i],1);
933 }
934 }
935 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,NODEVEC))
936 {
937 GetVectorsOfNodes(theElement,&cnt,vList);
938 for (i=0; i<cnt; i++)
939 {
940 if (DisposeConnectionFromVector(theGrid,vList[i])) RETURN(GM_ERROR);
941 SETVBUILDCON(vList[i],1);
942 }
943 }
944
945 return(GM_OK);
946 }
947
948 /****************************************************************************/
949 /** \brief Remove matrices
950 *
951 * @param theGrid - the grid to remove from
952 * @param theElement - that element
953 * @param Depth - that many slices around the element
954
955 This function removes connections concerning an element from the data structure
956 and stores flags saying: "connection has to be rebuild",
957 it does this in a neighborhood of the elem of depth Depth, where depth
958 is the distance in the element-neighborship-graph (see also FORMAT).
959
960 * @return <ul>
961 * <li> 0 if ok </li>
962 * <li> 1 if error occured. </li>
963 </ul>
964 */
965 /****************************************************************************/
966
DisposeConnectionFromElementInNeighborhood(GRID * theGrid,ELEMENT * theElement,INT Depth)967 static INT DisposeConnectionFromElementInNeighborhood (GRID *theGrid, ELEMENT *theElement, INT Depth)
968 {
969 INT i;
970
971 if (Depth < 0) RETURN (GM_ERROR);
972
973 if (theElement==NULL) RETURN (GM_OK);
974
975 /* create connection at that depth */
976 if (DisposeConnectionFromElement(theGrid,theElement))
977 RETURN (GM_ERROR);
978 SETEBUILDCON(theElement,1);
979
980 /* dispose connection in neighborhood */
981 if (Depth > 0)
982 {
983 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
984 if (DisposeConnectionFromElementInNeighborhood(theGrid,NBELEM(theElement,i),Depth-1))
985 RETURN (GM_ERROR);
986 }
987
988 RETURN (GM_OK);
989 }
990
DisposeConnectionsInNeighborhood(GRID * theGrid,ELEMENT * theElement)991 INT NS_DIM_PREFIX DisposeConnectionsInNeighborhood (GRID *theGrid, ELEMENT *theElement)
992 {
993 INT Depth;
994 Depth = (INT)(floor(0.5*(double)FMT_CONN_DEPTH_MAX(MGFORMAT(MYMG(theGrid)))));
995 return(DisposeConnectionFromElementInNeighborhood(theGrid,theElement,Depth));
996 }
997
DisposeConnectionsFromMultiGrid(MULTIGRID * theMG)998 INT NS_DIM_PREFIX DisposeConnectionsFromMultiGrid (MULTIGRID *theMG)
999 {
1000 INT i;
1001
1002 for (i=0; i<=TOPLEVEL(theMG); i++)
1003 {
1004 GRID *theGrid = GRID_ON_LEVEL(theMG,i);
1005 ELEMENT *theElement;
1006
1007 theGrid = GRID_ON_LEVEL(theMG,i);
1008 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL;
1009 theElement=SUCCE(theElement))
1010 if (DisposeConnectionsInNeighborhood(theGrid,theElement))
1011 REP_ERR_RETURN(1);
1012 }
1013
1014 return(0);
1015 }
1016
DisposeElementFromElementList(GRID * theGrid,NODE * theNode,ELEMENT * theElement)1017 INT NS_DIM_PREFIX DisposeElementFromElementList (GRID *theGrid, NODE *theNode,
1018 ELEMENT *theElement)
1019 {
1020 ELEMENTLIST *pel,*next;
1021
1022 pel = NODE_ELEMENT_LIST(theNode);
1023 if (pel == NULL) return(0);
1024 if (pel->el == theElement) {
1025 NDATA(theNode) = (void *) pel->next;
1026 return(PutFreeObject(theGrid->mg,pel,sizeof(ELEMENTLIST),MAOBJ));
1027 }
1028 next = pel->next;
1029 while (next != NULL) {
1030 if (next->el == theElement) {
1031 pel->next = next->next;
1032 return(PutFreeObject(theGrid->mg,next,sizeof(ELEMENTLIST),MAOBJ));
1033 }
1034 pel = next;
1035 next = pel->next;
1036 }
1037
1038 return(0);
1039 }
1040
DisposeElementList(GRID * theGrid,NODE * theNode)1041 INT NS_DIM_PREFIX DisposeElementList (GRID *theGrid, NODE *theNode)
1042 {
1043 ELEMENTLIST *pel,*next;
1044
1045 pel = NODE_ELEMENT_LIST(theNode);
1046 while (pel != NULL) {
1047 next = pel->next;
1048 if (PutFreeObject(theGrid->mg,pel,sizeof(ELEMENTLIST),MAOBJ))
1049 return(1);
1050 pel = next;
1051 }
1052 NDATA(theNode) = NULL;
1053
1054 return(0);
1055 }
1056
1057 /****************************************************************************/
1058 /** \brief Return pointer to matrix if it exists
1059
1060 * @param FromVector - starting vector of the Matrix
1061 * @param ToVector - destination vector of the Matrix
1062
1063 This function returns pointer to matrix if it exists. The function
1064 runs through the single linked list, since the list is
1065 assumed to be small (sparse matrix!) the cost is assumed to be negligible.
1066
1067 * @return <ul>
1068 * <li> pointer to Matrix, </li>
1069 * <li> NULL if Matrix does not exist. </li>
1070 </ul>
1071 */
1072 /****************************************************************************/
1073
GetMatrix(const VECTOR * FromVector,const VECTOR * ToVector)1074 MATRIX * NS_DIM_PREFIX GetMatrix (const VECTOR *FromVector, const VECTOR *ToVector)
1075 {
1076 MATRIX *theMatrix;
1077
1078 for (theMatrix=VSTART(FromVector); theMatrix!=NULL; theMatrix = MNEXT(theMatrix))
1079 if (MDEST(theMatrix)==ToVector)
1080 return (theMatrix);
1081
1082 /* return not found */
1083 return (NULL);
1084 }
1085
1086 /****************************************************************************/
1087 /** \brief Return pointer to connection if it exists
1088
1089 * @param FromVector - starting vector of the con
1090 * @param ToVector - destination vector of the con
1091
1092 This function returns pointer to connection if it exists.
1093
1094 * @return <ul>
1095 * <li> pointer to </li>
1096 * <li> NULL if connection does not exist. </li>
1097 </ul>
1098 */
1099 /****************************************************************************/
1100
GetConnection(const VECTOR * FromVector,const VECTOR * ToVector)1101 CONNECTION * NS_DIM_PREFIX GetConnection (const VECTOR *FromVector, const VECTOR *ToVector)
1102 {
1103 MATRIX *Matrix;
1104
1105 Matrix = GetMatrix(FromVector,ToVector);
1106 if (Matrix != NULL)
1107 return (MMYCON(Matrix));
1108
1109 /* return not found */
1110 return (NULL);
1111 }
1112
1113 /****************************************************************************/
1114 /** \brief Get a pointer list to all element data
1115
1116 * @param theElement - that element
1117 * @param cnt - how many vectors
1118 * @param vList - array to store vector list
1119
1120 This function returns a pointer to the VECTOR associated with the
1121 element (if the element is allowed to have one). cnt will either
1122 be 0 or 1.
1123
1124 * @return <ul>
1125 * <li> GM_OK if ok </li>
1126 * <li> GM_ERROR if error occured. </li>
1127 </ul>
1128 */
1129 /****************************************************************************/
1130
GetVectorsOfElement(const ELEMENT * theElement,INT * cnt,VECTOR ** vList)1131 INT NS_DIM_PREFIX GetVectorsOfElement (const ELEMENT *theElement, INT *cnt, VECTOR **vList)
1132 {
1133 *cnt = 0;
1134 if (EVECTOR(theElement) != NULL)
1135 vList[(*cnt)++] = EVECTOR(theElement);
1136
1137 IFDEBUG(gm,0)
1138 if (*cnt)
1139 {
1140 assert(vList[0] != NULL);
1141 assert(VOTYPE(vList[0]) == ELEMVEC);
1142 }
1143 ENDDEBUG
1144
1145 return(GM_OK);
1146 }
1147
1148 /****************************************************************************/
1149 /** \brief Get a pointer list to all side data
1150
1151 * @param theElement - that element
1152 * @param cnt - how many vectors
1153 * @param vList - array to store vector list
1154
1155 This function gets a pointer array to all VECTORs in sides of the given element.
1156
1157 * @return <ul>
1158 * <li> GM_OK if ok </li>
1159 * <li> GM_ERROR if error occured. </li>
1160 </ul>
1161 */
1162 /****************************************************************************/
1163
1164 #ifdef __THREEDIM__
GetVectorsOfSides(const ELEMENT * theElement,INT * cnt,VECTOR ** vList)1165 INT NS_DIM_PREFIX GetVectorsOfSides (const ELEMENT *theElement, INT *cnt, VECTOR **vList)
1166 {
1167 INT i;
1168
1169 *cnt = 0;
1170
1171 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
1172 if (SVECTOR(theElement,i) != NULL)
1173 vList[(*cnt)++] = SVECTOR(theElement,i);
1174
1175 IFDEBUG(gm,0)
1176 for (i=0; i<(*cnt); i++)
1177 {
1178 assert(vList[i] != NULL);
1179 assert(VOTYPE(vList[i]) == SIDEVEC);
1180 }
1181 ENDDEBUG
1182
1183 return(GM_OK);
1184 }
1185 #endif
1186
1187 /****************************************************************************/
1188 /** \brief Get a pointer list to all edge data
1189
1190 * @param theElement - that element
1191 * @param cnt - how many vectors
1192 * @param vList - array to store vector list
1193
1194 This function gets a pointer array to all VECTORs in edges of the element.
1195
1196 * @return <ul>
1197 * <li> GM_OK if ok </li>
1198 * <li> GM_ERROR if error occured. </li>
1199 </ul>
1200 */
1201 /****************************************************************************/
1202
GetVectorsOfEdges(const ELEMENT * theElement,INT * cnt,VECTOR ** vList)1203 INT NS_DIM_PREFIX GetVectorsOfEdges (const ELEMENT *theElement, INT *cnt, VECTOR **vList)
1204 {
1205 EDGE *theEdge;
1206 INT i;
1207
1208 *cnt = 0;
1209 for (i=0; i<EDGES_OF_ELEM(theElement); i++)
1210 if ((theEdge =
1211 GetEdge(CORNER(theElement,CORNER_OF_EDGE(theElement,i,0)),
1212 CORNER(theElement,CORNER_OF_EDGE(theElement,i,1)))) != NULL)
1213 if (EDVECTOR(theEdge) != NULL)
1214 vList[(*cnt)++] = EDVECTOR(theEdge);
1215
1216 IFDEBUG(gm,0)
1217 for (i=0; i<(*cnt); i++)
1218 {
1219 assert(vList[i] != NULL);
1220 assert(VOTYPE(vList[i]) == EDGEVEC);
1221 }
1222 ENDDEBUG
1223
1224 return(GM_OK);
1225 }
1226
1227 /****************************************************************************/
1228 /** \brief Get a pointer list to all node data
1229
1230 * @param theElement - that element
1231 * @param cnt - how many vectors
1232 * @param vList - array to store vector list
1233
1234 This function gets a pointer array to all VECTORs in nodes of the element.
1235
1236 * @return <ul>
1237 * <li> GM_OK if ok </li>
1238 * <li> GM_ERROR if error occured. </li>
1239 </ul>
1240 */
1241 /****************************************************************************/
1242
GetVectorsOfNodes(const ELEMENT * theElement,INT * cnt,VECTOR ** vList)1243 INT NS_DIM_PREFIX GetVectorsOfNodes (const ELEMENT *theElement, INT *cnt, VECTOR **vList)
1244 {
1245 INT i;
1246
1247 *cnt = 0;
1248 for (i=0; i<CORNERS_OF_ELEM(theElement); i++)
1249 if (NVECTOR(CORNER(theElement,i)) != NULL)
1250 vList[(*cnt)++] = NVECTOR(CORNER(theElement,i));
1251
1252 IFDEBUG(gm,0)
1253 for (i=0; i<(*cnt); i++)
1254 {
1255 assert(vList[i] != NULL);
1256 assert(VOTYPE(vList[i]) == NODEVEC);
1257 }
1258 ENDDEBUG
1259
1260 return (GM_OK);
1261 }
1262
1263 /****************************************************************************/
1264 /** \brief Get a pointer list to all vector data of specified object type
1265
1266 * @param theElement - that element
1267 * @param type - fill array with vectors of this type
1268 * @param cnt - how many vectors
1269 * @param vList - array to store vector list
1270
1271 This function returns a pointer array to all VECTORs of the element of the specified type.
1272
1273 * @return <ul>
1274 * <li> GM_OK if ok </li>
1275 * <li> GM_ERROR if error occured. </li>
1276 </ul>
1277 */
1278 /****************************************************************************/
1279
GetVectorsOfOType(const ELEMENT * theElement,INT type,INT * cnt,VECTOR ** vList)1280 INT NS_DIM_PREFIX GetVectorsOfOType (const ELEMENT *theElement, INT type, INT *cnt, VECTOR **vList)
1281 {
1282 switch (type)
1283 {
1284 case NODEVEC :
1285 return (GetVectorsOfNodes(theElement,cnt,vList));
1286 case ELEMVEC :
1287 return (GetVectorsOfElement(theElement,cnt,vList));
1288 case EDGEVEC :
1289 return (GetVectorsOfEdges(theElement,cnt,vList));
1290 #ifdef __THREEDIM__
1291 case SIDEVEC :
1292 return (GetVectorsOfSides(theElement,cnt,vList));
1293 #endif
1294 }
1295 RETURN (GM_ERROR);
1296 }
1297
1298 /****************************************************************************/
1299 /** \brief Remove vectors with non-matching type from list
1300
1301 * @param dt - data type including all vtypes needed
1302 * @param vec - vector list
1303 * @param cnt - number of vectors return in the list
1304
1305 This function removes vectors with non-matching types from the list.
1306
1307 * @return
1308 * Number of remaining vectors in the list
1309 */
1310 /****************************************************************************/
1311
DataTypeFilterVList(INT dt,VECTOR ** vec,INT * cnt)1312 INT NS_DIM_PREFIX DataTypeFilterVList (INT dt, VECTOR **vec, INT *cnt)
1313 {
1314 INT i,n;
1315
1316 n = *cnt;
1317 *cnt = 0;
1318 for (i=0; i<n; i++)
1319 if (VDATATYPE(vec[i]) & dt)
1320 vec[(*cnt)++] = vec[i];
1321
1322 return (*cnt);
1323 }
1324
1325 /****************************************************************************/
1326 /** \brief Get vector list including all vectors of specified vtypes
1327
1328 * @param theElement - pointer to an element
1329 * @param dt - data type including all vtypes needed
1330 * @param obj - flags for objects types needed
1331 * @param cnt - number of vectors return in the list
1332 * @param vec - vector list
1333
1334 This function gets a list of vectors of the specified vtypes corresponding to an element.
1335
1336 * @return <ul>
1337 * <li> GM_OK if ok </li>
1338 * <li> GM_ERROR if error occured </li>
1339 </ul>
1340 */
1341 /****************************************************************************/
1342
GetVectorsOfDataTypesInObjects(const ELEMENT * theElement,INT dt,INT obj,INT * cnt,VECTOR * VecList[])1343 INT NS_DIM_PREFIX GetVectorsOfDataTypesInObjects (const ELEMENT *theElement, INT dt, INT obj, INT *cnt, VECTOR *VecList[])
1344 {
1345 INT i,n;
1346
1347 *cnt = n = 0;
1348
1349 if (obj & BITWISE_TYPE(NODEVEC))
1350 {
1351 if (GetVectorsOfNodes(theElement,&i,VecList) != GM_OK)
1352 return(GM_ERROR);
1353 n += i;
1354 }
1355
1356 if (obj & BITWISE_TYPE(EDGEVEC))
1357 {
1358 if (GetVectorsOfEdges(theElement,&i,VecList+n) != GM_OK)
1359 return(GM_ERROR);
1360 n += i;
1361 }
1362
1363 if (obj & BITWISE_TYPE(ELEMVEC))
1364 {
1365 if (GetVectorsOfElement(theElement,&i,VecList+n) != GM_OK)
1366 return(GM_ERROR);
1367 n += i;
1368 }
1369
1370 #ifdef __THREEDIM__
1371 if (obj & BITWISE_TYPE(SIDEVEC))
1372 {
1373 if (GetVectorsOfSides(theElement,&i,VecList+n) != GM_OK)
1374 return(GM_ERROR);
1375 n += i;
1376 }
1377 #endif
1378
1379 *cnt = n;
1380
1381 /* remove vectors of types not needed */
1382 DataTypeFilterVList(dt,VecList,cnt);
1383
1384 return (GM_OK);
1385 }
1386
1387
1388
PrintVectorTriple(int i)1389 static void PrintVectorTriple (int i)
1390 {
1391 VECTOR *vec0 = GBNV_list[i];
1392 VECTOR *vec1 = GBNV_list[i+1];
1393 VECTOR *vec2 = GBNV_list[i+2];
1394 [[maybe_unused]] VERTEX *vtx0 = MYVERTEX((NODE*)VOBJECT(vec0));
1395 [[maybe_unused]] VERTEX *vtx1 = MYVERTEX((NODE*)VOBJECT(vec1));
1396 [[maybe_unused]] VERTEX *vtx2 = MYVERTEX((NODE*)VOBJECT(vec2));
1397 PrintDebug("0: VTYPE=%d XC=%.5g YC=%.5g\n",VTYPE(vec0),XC(vtx0),YC(vtx0));
1398 PrintDebug("1: VTYPE=%d XC=%.5g YC=%.5g\n",VTYPE(vec1),XC(vtx1),YC(vtx1));
1399 PrintDebug("2: VTYPE=%d XC=%.5g YC=%.5g\n",VTYPE(vec2),XC(vtx2),YC(vtx2));
1400 }
1401
1402 /****************************************************************************/
1403 /** \brief Prepare GetBoundaryNeighbourVectors
1404
1405 * @param theGrid - grid level
1406 * @param MaxListLen - max size of neighbourhood and therefore max list lenght of the
1407 VecList array of GetBoundaryNeighbourVectors
1408
1409 This function stores lists of boundary vector neighbourhoods in temp storage
1410 of the multigrid. The neighborhoods of the boundary vectors can be received
1411 via a call of GetBoundaryNeighbourVectors one by one.
1412
1413 * @return <ul>
1414 * <li> 0 if ok </li>
1415 * <li> >0 else </li>
1416 </ul>
1417
1418 * \sa
1419 GetBoundaryNeighbourVectors, ResetGetBoundaryNeighbourVectors, FinishBoundaryNeighbourVectors
1420 */
1421 /****************************************************************************/
1422
PrepareGetBoundaryNeighbourVectors(GRID * theGrid,INT * MaxListLen)1423 INT NS_DIM_PREFIX PrepareGetBoundaryNeighbourVectors (GRID *theGrid, INT *MaxListLen)
1424 {
1425 #ifdef __TWODIM__
1426 ELEMENT *elem;
1427 VECTOR *vec,*v0,*v1;
1428 INT i;
1429
1430 if (GBNV_list!=NULL)
1431 /* last process not finished properly by call of GetBoundaryNeighbourVectors */
1432 REP_ERR_RETURN(1);
1433
1434 /* count boundary node vectors */
1435 GBNV_n = 0;
1436 for (vec=FIRSTVECTOR(theGrid); vec!=NULL; vec=SUCCVC(vec))
1437 if (VOTYPE(vec)==NODEVEC)
1438 if (OBJT(MYVERTEX((NODE*)VOBJECT(vec)))==BVOBJ)
1439 GBNV_n++;
1440
1441 /* allocate list storage: 3 pointers each */
1442 GBNV_mg = MYMG(theGrid);
1443 MarkTmpMem(MGHEAP(GBNV_mg),&GBNV_MarkKey);
1444 GBNV_list = (VECTOR **) GetTmpMem(MGHEAP(GBNV_mg),3*GBNV_n*sizeof(VECTOR *),GBNV_MarkKey);
1445 if (GBNV_list==NULL)
1446 REP_ERR_RETURN(1);
1447
1448 /* store offset in vector index field */
1449 i = 0;
1450 for (vec=FIRSTVECTOR(theGrid); vec!=NULL; vec=SUCCVC(vec))
1451 if (VOTYPE(vec)==NODEVEC)
1452 if (OBJT(MYVERTEX((NODE*)VOBJECT(vec)))==BVOBJ)
1453 {
1454 VINDEX(vec) = i;
1455 GBNV_list[i] = vec;
1456 i += 3;
1457 }
1458
1459 /* loop elements and fill neighbours */
1460 /* TODO: parallel also orphan(?) elements to be complete */
1461 for (elem=FIRSTELEMENT(theGrid); elem!=NULL; elem=SUCCE(elem))
1462 if (OBJT(elem)==BEOBJ)
1463 for (i=0; i<SIDES_OF_ELEM(elem); i++)
1464 if (ELEM_BNDS(elem,i)!=NULL)
1465 {
1466 /* 2D: two corners */
1467 v0 = NVECTOR(CORNER(elem,CORNER_OF_SIDE(elem,i,0)));
1468 v1 = NVECTOR(CORNER(elem,CORNER_OF_SIDE(elem,i,1)));
1469
1470 ASSERT(GBNV_list[VINDEX(v0)]==v0);
1471 ASSERT(GBNV_list[VINDEX(v1)]==v1);
1472 GBNV_list[VINDEX(v0)+2] = v1;
1473 GBNV_list[VINDEX(v1)+1] = v0;
1474 }
1475 GBNV_curr = 0;
1476
1477 /* this is simple in 2D: center, pred and succ in positive sense */
1478 *MaxListLen = 3;
1479
1480 IFDEBUG(gm,2)
1481 PrintDebug("PrepareGetBoundaryNeighbourVectors:\n");
1482 for (i=0; i<GBNV_n; i++)
1483 PrintVectorTriple(3*i);
1484 ENDDEBUG
1485
1486 return (0);
1487
1488 #endif /* __TWODIM__ */
1489
1490 #ifdef __THREEDIM__
1491 /* 3D requires somewhat more work! */
1492
1493 /* but it should be possible to create an oriented list
1494 of neighbours for each boundary vector */
1495
1496 REP_ERR_RETURN (1);
1497 #endif /* __THREEDIM__ */
1498 }
1499
1500 /****************************************************************************/
1501 /** \brief Reset current neighbourhood to begin of list
1502
1503 This function resets the pointer to the current neighbourhood to the beginning of
1504 the list. GetBoundaryNeighbourVectors will start again with the first one.
1505
1506 * @return <ul>
1507 * <li> 0 if ok </li>
1508 * <li> >0 else </li>
1509 </ul>
1510
1511 * \sa
1512 GetBoundaryNeighbourVectors, PrepareGetBoundaryNeighbourVectors, FinishBoundaryNeighbourVectors
1513 */
1514 /****************************************************************************/
1515
ResetGetBoundaryNeighbourVectors(void)1516 INT NS_DIM_PREFIX ResetGetBoundaryNeighbourVectors (void)
1517 {
1518 if (GBNV_list==NULL)
1519 REP_ERR_RETURN(1);
1520
1521 GBNV_curr = 0;
1522 return (0);
1523 }
1524
1525 /****************************************************************************/
1526 /** \brief Get a neighbourhood of boundary vectors
1527
1528 * @param dt - datatypes of center vectors (bitwise)
1529 * @param obj - object types of center vectors (bitwise)
1530 * @param cnt - vector list length
1531 * @param VecList - vector list
1532 * @param end - if YES the currently returned list was the last one
1533
1534 This function returns a neighbourhood of boundary vectors, center first and the
1535 remainder oriented in positiv sense. The boundary vector heighbourhood lists
1536 are created by a call of PrepareGetBoundaryNeighbourVectors. Use
1537 ResetGetBoundaryNeighbourVectors to begin again with the first
1538 neighbourhood and finish processing the boundary vectors with a call of
1539 FinishBoundaryNeighbourVectors which releases the temporary storage
1540 occupied by PrepareGetBoundaryNeighbourVectors..
1541
1542 * @return <ul>
1543 * <li> 0 if ok </li>
1544 * <li> >0 else </li>
1545 </ul>
1546 * \sa
1547 PrepareGetBoundaryNeighbourVectors, ResetGetBoundaryNeighbourVectors, FinishBoundaryNeighbourVectors
1548 */
1549 /****************************************************************************/
1550
GetBoundaryNeighbourVectors(INT dt,INT obj,INT * cnt,VECTOR * VecList[])1551 INT NS_DIM_PREFIX GetBoundaryNeighbourVectors (INT dt, INT obj, INT *cnt, VECTOR *VecList[])
1552 {
1553 VECTOR *vec;
1554
1555 *cnt = 0;
1556
1557 if (GBNV_list==NULL)
1558 REP_ERR_RETURN(1);
1559
1560 /* find next center vec matching data type */
1561 for (vec=GBNV_list[GBNV_curr]; GBNV_curr<3*GBNV_n; GBNV_curr+=3, vec=GBNV_list[GBNV_curr])
1562 if (BITWISE_TYPE(VTYPE(vec)) & dt)
1563 break;
1564 if (GBNV_curr>=3*GBNV_n)
1565 /* no (more) vector found */
1566 return (0);
1567
1568 if (VOTYPE(vec)!=NODEVEC)
1569 REP_ERR_RETURN(1);
1570
1571 IFDEBUG(gm,2)
1572 PrintDebug("GetBoundaryNeighbourVectors:\n");
1573 PrintVectorTriple(GBNV_curr);
1574 ENDDEBUG
1575
1576 /* vector, pre and succ in positive sense */
1577 VecList[(*cnt)++] = GBNV_list[GBNV_curr];
1578 VecList[(*cnt)++] = GBNV_list[GBNV_curr+1];
1579 VecList[(*cnt)++] = GBNV_list[GBNV_curr+2];
1580
1581 /* move on to next position */
1582 GBNV_curr += 3;
1583
1584 return (0);
1585 }
1586
1587 /****************************************************************************/
1588 /** \brief Finish processing of boundary vectors
1589
1590 This function releases the temporary memory allocated by PrepareGetBoundaryNeighbourVectors
1591 from the multigrid heap.
1592
1593 * @return <ul>
1594 * <li> 0 if ok </li>
1595 * <li> >0 else </li>
1596 </ul>
1597
1598 * \sa
1599 PrepareGetBoundaryNeighbourVectors, ResetGetBoundaryNeighbourVectors, GetBoundaryNeighbourVectors
1600 */
1601 /****************************************************************************/
1602
FinishBoundaryNeighbourVectors()1603 INT NS_DIM_PREFIX FinishBoundaryNeighbourVectors ()
1604 {
1605 #ifdef __TWODIM
1606 if (ReleaseTmpMem(MGHEAP(GBNV_mg)),GBNV_MarkKey)
1607 REP_ERR_RETURN(1);
1608 #endif
1609
1610 GBNV_list = NULL;
1611
1612 return (0);
1613 }
1614
1615 /****************************************************************************/
1616 /** \brief Get vector list
1617
1618 * @param theGrid - pointer to a grid
1619 * @param theElement - pointer to an element
1620 * @param vec - vector list
1621
1622 This function gets a list of vectors corresponding to an element.
1623
1624 * @return <ul>
1625 * <li> number of components </li>
1626 * <li> -1 if error occured </li>
1627 </ul>
1628 */
1629 /****************************************************************************/
1630
GetAllVectorsOfElement(GRID * theGrid,ELEMENT * theElement,VECTOR ** vec)1631 INT NS_DIM_PREFIX GetAllVectorsOfElement (GRID *theGrid, ELEMENT *theElement, VECTOR **vec)
1632 {
1633 INT i;
1634 INT cnt;
1635
1636 cnt = 0;
1637 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,NODEVEC))
1638 {
1639 if (GetVectorsOfNodes(theElement,&i,vec) == GM_ERROR)
1640 RETURN(-1);
1641 cnt += i;
1642 }
1643 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,EDGEVEC))
1644 {
1645 if (GetVectorsOfEdges(theElement,&i,vec+cnt) == GM_ERROR)
1646 RETURN(-1);
1647 cnt += i;
1648 }
1649 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,ELEMVEC))
1650 {
1651 if (GetVectorsOfElement(theElement,&i,vec+cnt) == GM_ERROR)
1652 RETURN(-1);
1653 cnt += i;
1654 }
1655 #ifdef __THREEDIM__
1656 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,SIDEVEC))
1657 {
1658 if (GetVectorsOfSides(theElement,&i,vec+cnt) == GM_ERROR)
1659 RETURN(-1);
1660 cnt += i;
1661 }
1662 #endif
1663
1664 return (cnt);
1665 }
1666
1667 /****************************************************************************/
1668 /** \brief Remove all extra connections from the grid
1669
1670 * @param theGrid - grid to remove from
1671
1672 This function removes all extra connections from the grid, i.e. those
1673 that have been allocated with the CreateExtraConnection function.
1674
1675 * @return <ul>
1676 * <li> 0 if ok
1677 * <li> 1 if error occured
1678 </ul>
1679 */
1680 /****************************************************************************/
1681
DisposeExtraConnections(GRID * theGrid)1682 INT NS_DIM_PREFIX DisposeExtraConnections (GRID *theGrid)
1683 {
1684 VECTOR *theVector;
1685 MATRIX *theMatrix, *nextMatrix;
1686 CONNECTION *theCon;
1687
1688 for (theVector=FIRSTVECTOR(theGrid); theVector!=NULL; theVector=SUCCVC(theVector))
1689 {
1690 theMatrix = VSTART(theVector);
1691 while (theMatrix!=NULL)
1692 {
1693 nextMatrix = MNEXT(theMatrix);
1694 theCon = MMYCON(theMatrix);
1695 if (CEXTRA(theCon)) DisposeConnection(theGrid,theCon);
1696 theMatrix = nextMatrix;
1697 }
1698 }
1699 return(GM_OK);
1700 }
1701
DisposeConnectionsInGrid(GRID * theGrid)1702 INT NS_DIM_PREFIX DisposeConnectionsInGrid (GRID *theGrid)
1703 {
1704 VECTOR *theVector;
1705 MATRIX *theMatrix, *nextMatrix;
1706 CONNECTION *theCon;
1707
1708 for (theVector=PFIRSTVECTOR(theGrid); theVector!=NULL; theVector=SUCCVC(theVector))
1709 {
1710 theMatrix = VSTART(theVector);
1711 while (theMatrix!=NULL)
1712 {
1713 nextMatrix = MNEXT(theMatrix);
1714 theCon = MMYCON(theMatrix);
1715 DisposeConnection(theGrid,theCon);
1716 theMatrix = nextMatrix;
1717 }
1718 }
1719 return(GM_OK);
1720 }
1721
1722 /****************************************************************************/
1723 /** \brief Create connections of two elements
1724
1725 * @param theGrid - pointer to grid
1726 * @param Elem0,Elem1 - elements to be connected
1727 * @param ActDepth - distance of the two elements in the element neighborship graph
1728 * @param ConDepth - Array containing the connection depth desired
1729 * @param MatSize - Array containing the connection size. This is constructed from the information in the FORMAT.
1730
1731 This function creates connections between all the VECTORs associated
1732 with two given elements according to the specifications in the ConDepth array.
1733
1734 * @return <ul>
1735 * <li> 0 if ok
1736 * <li> 1 if error occured.
1737 </ul>
1738 */
1739 /****************************************************************************/
1740
ElementElementCreateConnection(GRID * theGrid,ELEMENT * Elem0,ELEMENT * Elem1,INT ActDepth,INT * ConDepth,INT * MatSize)1741 static INT ElementElementCreateConnection (GRID *theGrid, ELEMENT *Elem0, ELEMENT *Elem1, INT ActDepth, INT *ConDepth, INT *MatSize)
1742 {
1743 INT cnt0,cnt1,i,j,itype,jtype,mtype,size;
1744 VECTOR *vec0[MAX_SIDES_OF_ELEM+MAX_EDGES_OF_ELEM+MAX_CORNERS_OF_ELEM+1];
1745 VECTOR *vec1[MAX_SIDES_OF_ELEM+MAX_EDGES_OF_ELEM+MAX_CORNERS_OF_ELEM+1];
1746
1747 cnt0 = GetAllVectorsOfElement(theGrid,Elem0,vec0);
1748 if (Elem0 == Elem1)
1749 {
1750 for (i=0; i<cnt0; i++)
1751 {
1752 itype=VTYPE(vec0[i]);
1753 for (j=i; j<cnt0; j++)
1754 {
1755 if (i==j)
1756 {
1757 mtype = DIAGMATRIXTYPE(itype);
1758 size = MatSize[mtype];
1759 }
1760 else
1761 {
1762 jtype = VTYPE(vec0[j]);
1763 size = MatSize[MATRIXTYPE(jtype,itype)];
1764 mtype = MATRIXTYPE(itype,jtype);
1765 if (MatSize[mtype]>size) size=MatSize[mtype];
1766 }
1767 if (size > 0)
1768 if (ActDepth <= ConDepth[mtype])
1769 if (CreateConnection(theGrid,vec0[i],vec0[j])==NULL)
1770 RETURN(GM_ERROR);
1771 }
1772 }
1773 return (0);
1774 }
1775
1776 cnt1 = GetAllVectorsOfElement(theGrid,Elem1,vec1);
1777 for (i=0; i<cnt0; i++)
1778 {
1779 itype=VTYPE(vec0[i]);
1780 for (j=0; j<cnt1; j++)
1781 {
1782 if (vec0[i]==vec1[j])
1783 {
1784 mtype = DIAGMATRIXTYPE(itype);
1785 size = MatSize[mtype];
1786 }
1787 else
1788 {
1789 jtype = VTYPE(vec1[j]);
1790 size = MatSize[MATRIXTYPE(jtype,itype)];
1791 mtype = MATRIXTYPE(itype,jtype);
1792 if (MatSize[mtype]>size) size=MatSize[mtype];
1793 }
1794 if (size > 0)
1795 if (ActDepth <= ConDepth[mtype])
1796 if (CreateConnection(theGrid,vec0[i],vec1[j])==NULL)
1797 RETURN(GM_ERROR);
1798 }
1799 }
1800
1801 return (0);
1802 }
1803
1804 /****************************************************************************/
1805 /** \brief Get pointers to elements having a common side
1806
1807 * @param theVector - given vector associated with a side of an element in 3D
1808 * @param Elements - array to be filled with two pointers of elements
1809 * @param Sides - array to be filled with side number within respective element
1810
1811 Given a VECTOR associated with the side of an element, this function
1812 retrieves pointers to at most two elements having this side in common.
1813 If the side is part of the exterior boundary of the domain, then
1814 Elements[1] will be the NULL pointer.
1815
1816 * @return <ul>
1817 * <li> 0 if ok
1818 * <li> 1 if error occured.
1819 </ul>
1820 */
1821 /****************************************************************************/
1822
1823 #ifdef __THREEDIM__
GetElementInfoFromSideVector(const VECTOR * theVector,ELEMENT ** Elements,INT * Sides)1824 INT NS_DIM_PREFIX GetElementInfoFromSideVector (const VECTOR *theVector, ELEMENT **Elements, INT *Sides)
1825 {
1826 INT i;
1827 ELEMENT *theNeighbor;
1828
1829 if (VOTYPE(theVector) != SIDEVEC)
1830 RETURN (1);
1831 Elements[0] = (ELEMENT *)VOBJECT(theVector);
1832 Sides[0] = VECTORSIDE(theVector);
1833
1834 /* find neighbor */
1835 Elements[1] = theNeighbor = NBELEM(Elements[0],Sides[0]);
1836 if (theNeighbor == NULL) return (0);
1837
1838 /* search the side */
1839 for (i=0; i<SIDES_OF_ELEM(theNeighbor); i++)
1840 if (NBELEM(theNeighbor,i) == Elements[0])
1841 break;
1842
1843 /* found ? */
1844 if (i<SIDES_OF_ELEM(theNeighbor))
1845 {
1846 Sides[1] = i;
1847 return (0);
1848 }
1849 RETURN (1);
1850 }
1851 #endif
1852
1853 /****************************************************************************/
1854 /** \brief Reset all USED flags in neighborhood of an element
1855
1856 * @param theElement - given element
1857 * @param ActDepth - recursion depth
1858 * @param MaxDepth - end of recursion
1859
1860 This function calls itself recursively and resets all USED flags in the
1861 neighborhood of depth MaxDepth of theElement. For the first call
1862 ActDepth should be zero.
1863
1864 * @return <ul>
1865 * <li> 0 if ok
1866 * <li> 1 if error occured.
1867 </ul>
1868 */
1869 /****************************************************************************/
1870
ResetUsedFlagInNeighborhood(ELEMENT * theElement,INT ActDepth,INT MaxDepth)1871 static INT ResetUsedFlagInNeighborhood (ELEMENT *theElement, INT ActDepth, INT MaxDepth)
1872 {
1873 int i;
1874
1875 /* is anything to do ? */
1876 if (theElement==NULL) return (0);
1877
1878 /* action */
1879 if (ActDepth>=0) SETUSED(theElement,0);
1880
1881 /* call all neighbors recursively */
1882 if (ActDepth<MaxDepth)
1883 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
1884 if (ResetUsedFlagInNeighborhood(NBELEM(theElement,i),ActDepth+1,MaxDepth)) RETURN (1);
1885
1886 return (0);
1887 }
1888
ConnectWithNeighborhood(ELEMENT * theElement,GRID * theGrid,ELEMENT * centerElement,INT * ConDepth,INT * MatSize,INT ActDepth,INT MaxDepth)1889 static INT ConnectWithNeighborhood (ELEMENT *theElement, GRID *theGrid, ELEMENT *centerElement, INT *ConDepth, INT *MatSize, INT ActDepth, INT MaxDepth)
1890 {
1891 int i;
1892
1893 /* is anything to do ? */
1894 if (theElement==NULL) return (0);
1895
1896 /* action */
1897 if (ActDepth>=0)
1898 if (ElementElementCreateConnection(theGrid,centerElement,theElement,
1899 ActDepth,ConDepth,MatSize))
1900 RETURN (1);
1901
1902 /* call all neighbors recursively */
1903 if (ActDepth<MaxDepth)
1904 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
1905 if (ConnectWithNeighborhood(NBELEM(theElement,i),theGrid,
1906 centerElement,ConDepth,MatSize,
1907 ActDepth+1,MaxDepth)) RETURN (1);
1908
1909 return (0);
1910 }
1911
1912 /****************************************************************************/
1913 /** \brief Create connection of an element
1914
1915 * @param theGrid - pointer to grid
1916 * @param theElement - pointer to an element
1917
1918 This function creates connection for all VECTORs of an element
1919 with the depth specified in the FORMAT structure.
1920
1921 * @return <ul>
1922 * <li> 0 if ok
1923 * <li> 1 if error occured.
1924 </ul>
1925 */
1926 /****************************************************************************/
1927
CreateConnectionsInNeighborhood(GRID * theGrid,ELEMENT * theElement)1928 INT NS_DIM_PREFIX CreateConnectionsInNeighborhood (GRID *theGrid, ELEMENT *theElement)
1929 {
1930 INT MaxDepth;
1931 INT *ConDepth;
1932 INT *MatSize;
1933
1934 /* set pointers */
1935 FORMAT* theFormat = theGrid->mg->theFormat.get();
1936 MaxDepth = FMT_CONN_DEPTH_MAX(theFormat);
1937 ConDepth = FMT_CONN_DEPTH_PTR(theFormat);
1938 MatSize = FMT_S_MATPTR(theFormat);
1939
1940 /* reset used flags in neighborhood */
1941 if (ResetUsedFlagInNeighborhood(theElement,0,MaxDepth))
1942 RETURN (1);
1943
1944 /* create connection in neighborhood */
1945 if (ConnectWithNeighborhood(theElement,theGrid,theElement,ConDepth,
1946 MatSize,0,MaxDepth))
1947 RETURN (1);
1948
1949 return (0);
1950 }
1951
1952 /****************************************************************************/
1953 /** \brief Create connection of an inserted element
1954
1955 This function creates connection of an inserted element ,
1956 i.e. an element is inserted interactively by the user.
1957
1958 * @return <ul>
1959 * <li> 0 if ok
1960 * <li> 1 if error occured.
1961 </ul>
1962 */
1963 /****************************************************************************/
1964
ConnectInsertedWithNeighborhood(ELEMENT * theElement,GRID * theGrid,INT ActDepth,INT MaxDepth)1965 static INT ConnectInsertedWithNeighborhood (ELEMENT *theElement, GRID *theGrid, INT ActDepth, INT MaxDepth)
1966 {
1967 int i;
1968
1969 /* is anything to do ? */
1970 if (theElement==NULL) return (0);
1971
1972 /* action */
1973 if (ActDepth>=0)
1974 if (CreateConnectionsInNeighborhood(theGrid,theElement))
1975 RETURN (1);
1976
1977 /* call all neighbors recursively */
1978 if (ActDepth<MaxDepth)
1979 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
1980 if (ConnectInsertedWithNeighborhood(NBELEM(theElement,i),theGrid,ActDepth+1,MaxDepth)) RETURN (1);
1981
1982 return (0);
1983 }
1984
1985 /****************************************************************************/
1986 /** \brief Create connection of an inserted element
1987
1988 * @param theGrid - grid level
1989 * @param theElement - pointer to an element
1990
1991 This function creates connections of an inserted element ,
1992 i.e. when an element is inserted interactively by the user. This function
1993 is inefficient and only intended for that purpose.
1994
1995 * @return <ul>
1996 * <li> 0 if ok
1997 * <li> 1 if error occured.
1998 </ul>
1999 */
2000 /****************************************************************************/
2001
InsertedElementCreateConnection(GRID * theGrid,ELEMENT * theElement)2002 INT NS_DIM_PREFIX InsertedElementCreateConnection (GRID *theGrid, ELEMENT *theElement)
2003 {
2004 INT MaxDepth;
2005
2006 if (!MG_COARSE_FIXED(MYMG(theGrid)))
2007 RETURN (1);
2008
2009 /* set pointers */
2010 FORMAT* theFormat = theGrid->mg->theFormat.get();
2011 MaxDepth = (INT)(floor(0.5*(double)FMT_CONN_DEPTH_MAX(theFormat)));
2012
2013 /* reset used flags in neighborhood */
2014 if (ResetUsedFlagInNeighborhood(theElement,0,MaxDepth))
2015 RETURN (1);
2016
2017 /* call 'CreateConnectionsInNeighborhood'
2018 in a neighborhood of theElement */
2019 if (ConnectInsertedWithNeighborhood (theElement,theGrid,0,MaxDepth))
2020 RETURN (1);
2021
2022 return (0);
2023 }
2024
2025 /****************************************************************************/
2026 /** \brief Create all connections needed on a grid level
2027
2028 * @param theGrid - pointer to grid
2029
2030 This function creates all connections needed on the grid.
2031
2032 * @return <ul>
2033 * <li> 0 if ok
2034 * <li> 1 if error occured.
2035 </ul>
2036 */
2037 /****************************************************************************/
2038
GridCreateConnection(GRID * theGrid)2039 INT NS_DIM_PREFIX GridCreateConnection (GRID *theGrid)
2040 {
2041 ELEMENT *theElement;
2042 VECTOR *vList[20];
2043 INT i,cnt;
2044
2045 #if defined(ModelP) and defined(__THREEDIM__)
2046 auto& dddContext = theGrid->dddContext();
2047 #endif
2048
2049 if (!MG_COARSE_FIXED(MYMG(theGrid)))
2050 RETURN (1);
2051
2052 /* lets see if there's something to do */
2053 if (theGrid == NULL)
2054 return (0);
2055
2056 #if defined(ModelP) and defined(__THREEDIM__)
2057 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,EDGEVEC))
2058 DDD_XferBegin(dddContext);
2059 #endif
2060
2061 /* set EBUILDCON-flags also in elements accessing a vector with VBUILDCON true */
2062 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL; theElement=SUCCE(theElement))
2063 {
2064 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,EDGEVEC))
2065 {
2066 EDGE *ed;
2067
2068 for (i=0; i<EDGES_OF_ELEM(theElement); i++) {
2069 ed = GetEdge(CORNER(theElement,
2070 CORNER_OF_EDGE(theElement,i,0)),
2071 CORNER(theElement,
2072 CORNER_OF_EDGE(theElement,i,1)));
2073 ASSERT (ed != NULL);
2074 if (EDVECTOR(ed) == NULL) {
2075 CreateVector(theGrid,EDGEVEC,(GEOM_OBJECT *)ed,vList);
2076 EDVECTOR(ed) = vList[0];
2077 #if defined(ModelP) and defined(__THREEDIM__)
2078 SETPRIO(dddContext, EDVECTOR(ed),PRIO(ed));
2079 #endif
2080 }
2081 }
2082 }
2083
2084 #ifdef ModelP
2085 if (!EMASTER(theElement)) continue;
2086 #endif
2087
2088 /* see if it is set */
2089 if (EBUILDCON(theElement)) continue;
2090
2091 /* check flags in vectors */
2092 #ifdef __THREEDIM__
2093 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,SIDEVEC))
2094 {
2095 GetVectorsOfSides(theElement,&cnt,vList);
2096 for (i=0; i<cnt; i++)
2097 if (VBUILDCON(vList[i])) {SETEBUILDCON(theElement,1); break;}
2098 }
2099 if (EBUILDCON(theElement)) continue;
2100 #endif
2101 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,EDGEVEC))
2102 {
2103 GetVectorsOfEdges(theElement,&cnt,vList);
2104 for (i=0; i<cnt; i++)
2105 if (VBUILDCON(vList[i])) {SETEBUILDCON(theElement,1); break;}
2106 }
2107 if (EBUILDCON(theElement)) continue;
2108 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,NODEVEC))
2109 {
2110 GetVectorsOfNodes(theElement,&cnt,vList);
2111 for (i=0; i<cnt; i++)
2112 if (VBUILDCON(vList[i])) {SETEBUILDCON(theElement,1); break;}
2113 }
2114 }
2115
2116 #if defined(ModelP) and defined(__THREEDIM__)
2117 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,EDGEVEC))
2118 DDD_XferEnd(theGrid->dddContext());
2119 #endif
2120
2121 /* run over all elements with EBUILDCON true and build connections */
2122 for (theElement=FIRSTELEMENT(theGrid); theElement!=NULL;
2123 theElement=SUCCE(theElement))
2124 if (EBUILDCON(theElement)) /* this is the trigger ! */
2125 if (CreateConnectionsInNeighborhood(theGrid,theElement))
2126 RETURN (1);
2127
2128 return(GM_OK);
2129 }
2130
2131
2132 #ifdef ModelP
Gather_VectorVNew(DDD::DDDContext &,DDD_OBJ obj,void * data)2133 static int Gather_VectorVNew (DDD::DDDContext&, DDD_OBJ obj, void *data)
2134 {
2135 VECTOR *theVector = (VECTOR *)obj;
2136
2137 ((INT *)data)[0] = VNEW(theVector);
2138
2139 return(0);
2140 }
2141
Scatter_VectorVNew(DDD::DDDContext &,DDD_OBJ obj,void * data)2142 static int Scatter_VectorVNew (DDD::DDDContext&, DDD_OBJ obj, void *data)
2143 {
2144 VECTOR *theVector = (VECTOR *)obj;
2145
2146 SETVNEW(theVector,MAX(VNEW(theVector),((INT *)data)[0]));
2147
2148 return(0);
2149 }
2150
Scatter_GhostVectorVNew(DDD::DDDContext &,DDD_OBJ obj,void * data)2151 static int Scatter_GhostVectorVNew (DDD::DDDContext&, DDD_OBJ obj, void *data)
2152 {
2153 VECTOR *theVector = (VECTOR *)obj;
2154
2155 SETVNEW(theVector,((INT *)data)[0]);
2156
2157 return(0);
2158 }
2159 #endif
2160
SetSurfaceClasses(MULTIGRID * theMG)2161 INT NS_DIM_PREFIX SetSurfaceClasses (MULTIGRID *theMG)
2162 {
2163 GRID *theGrid;
2164 ELEMENT *theElement; VECTOR *v;
2165 INT level,fullrefine;
2166
2167 level = TOPLEVEL(theMG);
2168 if (level > 0) {
2169 theGrid = GRID_ON_LEVEL(theMG,TOPLEVEL(theMG));
2170 ClearVectorClasses(theGrid);
2171 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL;
2172 theElement=SUCCE(theElement))
2173 {
2174 if (MinNodeClass(theElement)==3)
2175 SeedVectorClasses(theGrid,theElement);
2176 }
2177 PropagateVectorClasses(theGrid);
2178 theGrid = GRID_ON_LEVEL(theMG,0);
2179 ClearNextVectorClasses(theGrid);
2180 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL;
2181 theElement=SUCCE(theElement))
2182 {
2183 if (MinNextNodeClass(theElement)==3)
2184 SeedNextVectorClasses(theGrid,theElement);
2185 }
2186 PropagateNextVectorClasses(theGrid);
2187 }
2188 for (level--; level>0; level--)
2189 {
2190 theGrid = GRID_ON_LEVEL(theMG,level);
2191 ClearVectorClasses(theGrid);
2192 ClearNextVectorClasses(theGrid);
2193 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL;
2194 theElement=SUCCE(theElement)) {
2195 if (MinNodeClass(theElement)==3)
2196 SeedVectorClasses(theGrid,theElement);
2197 if (MinNextNodeClass(theElement)==3)
2198 SeedNextVectorClasses(theGrid,theElement);
2199 }
2200 PropagateVectorClasses(theGrid);
2201 PropagateNextVectorClasses(theGrid);
2202 }
2203
2204 fullrefine = TOPLEVEL(theMG);
2205 for (level=TOPLEVEL(theMG); level>=0; level--)
2206 {
2207 theGrid = GRID_ON_LEVEL(theMG,level);
2208 for (v=PFIRSTVECTOR(theGrid); v!= NULL; v=SUCCVC(v)) {
2209 SETNEW_DEFECT(v,(VCLASS(v)>=2));
2210 SETFINE_GRID_DOF(v,((VCLASS(v)>=2)&&(VNCLASS(v)<=1)));
2211 if (FINE_GRID_DOF(v))
2212 fullrefine = level;
2213 }
2214 }
2215 #ifdef ModelP
2216 fullrefine = UG_GlobalMinINT(theMG->ppifContext(), fullrefine);
2217 #endif
2218
2219 FULLREFINELEVEL(theMG) = fullrefine;
2220
2221 return(0);
2222 }
2223
2224 /****************************************************************************/
2225 /** \brief Creates the algebra for a grid
2226
2227 * @param theGrid - pointer to grid
2228
2229 This function allocates VECTORs in all geometrical objects of the grid
2230 (where described by format) and then calls GridCreateConnection.
2231
2232 * @return <ul>
2233 * <li> GM_OK if ok
2234 </ul>
2235 */
2236 /****************************************************************************/
2237
CreateAlgebra(MULTIGRID * theMG)2238 INT NS_DIM_PREFIX CreateAlgebra (MULTIGRID *theMG)
2239 {
2240 GRID *g;
2241 FORMAT *fmt;
2242 VECTOR *vec;
2243 #ifdef __THREEDIM__
2244 VECTOR *nbvec;
2245 ELEMENT *nbelem;
2246 INT j,n;
2247 #endif
2248 ELEMENT *elem;
2249 NODE *nd;
2250 LINK *li;
2251 EDGE *ed;
2252 INT side,i;
2253
2254
2255 if (MG_COARSE_FIXED(theMG) == false) {
2256 for (i=0; i<=TOPLEVEL(theMG); i++) {
2257 g = GRID_ON_LEVEL(theMG,i);
2258
2259 if (NVEC(g)>0)
2260 continue; /* skip this level */
2261
2262 fmt = g->mg->theFormat.get();
2263
2264 /* loop nodes and edges */
2265 for (nd=PFIRSTNODE(g); nd!=NULL; nd=SUCCN(nd)) {
2266 /* node vector */
2267 if (FMT_USES_OBJ(fmt,NODEVEC))
2268 {
2269 ASSERT(NVECTOR(nd)==NULL);
2270 if (CreateVector (g,NODEVEC,(GEOM_OBJECT *)nd,&vec))
2271 REP_ERR_RETURN (GM_ERROR);
2272 NVECTOR(nd) = vec;
2273 }
2274
2275 /* edge vectors */
2276 if (FMT_USES_OBJ(fmt,EDGEVEC))
2277 for (li=START(nd); li!=NULL; li=NEXT(li)) {
2278 ed = MYEDGE(li);
2279 if (li==LINK0(ed)) /* to avoid double access of edges*/
2280 {
2281 ASSERT(EDVECTOR(ed)==NULL);
2282 if (CreateVector(g,EDGEVEC,(GEOM_OBJECT *)ed,&vec))
2283 REP_ERR_RETURN (GM_ERROR);
2284 EDVECTOR(ed) = vec;
2285 }
2286 }
2287 }
2288
2289 /* loop elements and element sides */
2290 for (elem=PFIRSTELEMENT(g); elem!=NULL; elem=SUCCE(elem)) {
2291 /* to tell GridCreateConnection to build connections */
2292 if (EMASTER(elem)) SETEBUILDCON(elem,1);
2293
2294 /* element vector */
2295 if (FMT_USES_OBJ(fmt,ELEMVEC)) {
2296 ASSERT(EVECTOR(elem)==NULL);
2297 if (CreateVector (g,ELEMVEC,(GEOM_OBJECT *)elem,&vec))
2298 REP_ERR_RETURN (GM_ERROR);
2299 SET_EVECTOR(elem,vec);
2300 }
2301
2302 /* side vectors */
2303 if (FMT_USES_OBJ(fmt,SIDEVEC))
2304 for (side=0; side<SIDES_OF_ELEM(elem); side++)
2305 if (SVECTOR(elem,side)==NULL) {
2306 if (CreateSideVector (g,side,
2307 (GEOM_OBJECT *)elem,&vec))
2308 REP_ERR_RETURN (GM_ERROR);
2309 SET_SVECTOR(elem,side,vec);
2310 }
2311 }
2312 }
2313 #ifdef __THREEDIM__
2314 /* dispose doubled side vectors */
2315 if (FMT_USES_OBJ(fmt,SIDEVEC))
2316 for (elem=PFIRSTELEMENT(g); elem!=NULL; elem=SUCCE(elem))
2317 {
2318 for(side=0; side<SIDES_OF_ELEM(elem); side++)
2319 {
2320 if(OBJT(elem)==BEOBJ)
2321 {
2322 if(INNER_SIDE(elem,side))
2323 {
2324 nbelem = NBELEM(elem,side);
2325 ASSERT(nbelem!=NULL);
2326 vec=SVECTOR(elem,side);
2327 n=0;
2328 for(j=0; j<SIDES_OF_ELEM(nbelem); j++)
2329 {
2330 nbvec=SVECTOR(nbelem,j);
2331
2332 /* doubled sidevectors ? */
2333 if(NBELEM(nbelem,j)==elem)
2334 {
2335 if(vec!=nbvec)
2336 {
2337 if(DisposeVector(g,nbvec))
2338 REP_ERR_RETURN(GM_ERROR);
2339 SET_SVECTOR(nbelem,j,vec);
2340 SETVCOUNT(vec,2); /* PB, 25 Sep 2005: changed from 1 to 2 */
2341 }
2342 }
2343 }
2344 n++;
2345 ASSERT(n==1);
2346 }
2347 }else
2348 {
2349 nbelem = NBELEM(elem,side);
2350 ASSERT(nbelem!=NULL);
2351 vec=SVECTOR(elem,side);
2352 n=0;
2353 for(j=0; j<SIDES_OF_ELEM(nbelem); j++)
2354 {
2355 nbvec=SVECTOR(nbelem,j);
2356 ASSERT(nbvec!=NULL);
2357
2358 /* doubled sidevectors ? */
2359 if(NBELEM(nbelem,j)==elem)
2360 {
2361 if(vec!=nbvec)
2362 {
2363 if(DisposeVector(g,nbvec))
2364 REP_ERR_RETURN(GM_ERROR);
2365 SET_SVECTOR(nbelem,j,vec);
2366 SETVCOUNT(vec,2); /* PB, 25 Sep 2005: changed from 1 to 2 */
2367 }
2368 }
2369 }
2370 n++;
2371 ASSERT(n==1);
2372 }
2373 }
2374 }
2375 #endif
2376 MG_COARSE_FIXED(theMG) = true;
2377
2378 /* now connections */
2379 if (MGCreateConnection(theMG))
2380 REP_ERR_RETURN (1);
2381 }
2382
2383 /* now we should be safe to clear the InsertElement face map */
2384 theMG->facemap.clear();
2385
2386 #ifdef ModelP
2387 #ifndef __EXCHANGE_CONNECTIONS__
2388 MGCreateConnection(theMG);
2389 #endif
2390 /* update VNEW-flags */
2391 auto& context = theMG->dddContext();
2392 const auto& dddctrl = ddd_ctrl(context);
2393 DDD_IFExchange(context,
2394 dddctrl.BorderVectorSymmIF,sizeof(INT),
2395 Gather_VectorVNew,Scatter_VectorVNew);
2396 DDD_IFOneway(context,
2397 dddctrl.VectorIF,IF_FORWARD,sizeof(INT),
2398 Gather_VectorVNew,Scatter_GhostVectorVNew);
2399 #else
2400 MGCreateConnection(theMG);
2401 #endif
2402
2403 SetSurfaceClasses(theMG);
2404
2405 return(GM_OK);
2406 }
2407
2408 /****************************************************************************/
2409 /** \brief Create all connections in multigrid
2410
2411 * @param theMG - pointer to mulrigrid
2412
2413 This function creates all connections in multigrid.
2414
2415 * @return <ul>
2416 * <li> 0 if ok
2417 * <li> 1 if error occured.
2418 </ul>
2419 */
2420 /****************************************************************************/
2421
MGCreateConnection(MULTIGRID * theMG)2422 INT NS_DIM_PREFIX MGCreateConnection (MULTIGRID *theMG)
2423 {
2424 INT i;
2425 GRID *theGrid;
2426 ELEMENT *theElement;
2427
2428 if (!MG_COARSE_FIXED(theMG))
2429 RETURN (1);
2430
2431 for (i=0; i<=theMG->topLevel; i++)
2432 {
2433 theGrid = GRID_ON_LEVEL(theMG,i);
2434 for (theElement=FIRSTELEMENT(theGrid); theElement!=NULL; theElement=SUCCE(theElement))
2435 SETEBUILDCON(theElement,1);
2436 if (GridCreateConnection(theGrid)) RETURN (1);
2437 }
2438
2439 return (0);
2440 }
2441
2442
PrepareAlgebraModification(MULTIGRID * theMG)2443 INT NS_DIM_PREFIX PrepareAlgebraModification (MULTIGRID *theMG)
2444 {
2445 int j,k;
2446 ELEMENT *theElement;
2447 VECTOR *theVector;
2448 MATRIX *theMatrix;
2449
2450 j = theMG->topLevel;
2451 for (k=0; k<=j; k++)
2452 {
2453 for (theElement=PFIRSTELEMENT(GRID_ON_LEVEL(theMG,k)); theElement!=NULL; theElement=SUCCE(theElement))
2454 {
2455 SETUSED(theElement,0);
2456 SETEBUILDCON(theElement,0);
2457 }
2458 for (theVector=PFIRSTVECTOR(GRID_ON_LEVEL(theMG,k)); theVector!= NULL; theVector=SUCCVC(theVector))
2459 SETVBUILDCON(theVector,0);
2460 for (theVector=PFIRSTVECTOR(GRID_ON_LEVEL(theMG,k)); theVector!= NULL; theVector=SUCCVC(theVector))
2461 {
2462 SETVNEW(theVector,0);
2463 for (theMatrix=VSTART(theVector); theMatrix!=NULL; theMatrix = MNEXT(theMatrix))
2464 SETMNEW(theMatrix,0);
2465 }
2466 }
2467
2468 return (0);
2469 }
2470
2471 /****************************************************************************/
2472 /** \brief Check connection of two elements
2473
2474 * @param theGrid - grid level to check
2475 * @param Elem0,Elem1 - elements to be checked.
2476 * @param ActDepth - recursion depth
2477 * @param ConDepth - connection depth as provided by format.
2478 * @param MatSize - matrix sizes as provided by format
2479
2480 This function recursively checks connection of two elements.
2481
2482 * @return <ul>
2483 * <li> 0 if ok
2484 * <li> != 0 if errors occured.
2485 </ul>
2486 */
2487 /****************************************************************************/
2488
ElementElementCheck(GRID * theGrid,ELEMENT * Elem0,ELEMENT * Elem1,INT ActDepth,INT * ConDepth,INT * MatSize)2489 static INT ElementElementCheck (GRID *theGrid, ELEMENT *Elem0, ELEMENT *Elem1, INT ActDepth, INT *ConDepth, INT *MatSize)
2490 {
2491 INT cnt0,cnt1,i,j,itype,jtype,mtype,size;
2492 VECTOR *vec0[MAX_SIDES_OF_ELEM+MAX_EDGES_OF_ELEM+MAX_CORNERS_OF_ELEM+1];
2493 VECTOR *vec1[MAX_SIDES_OF_ELEM+MAX_EDGES_OF_ELEM+MAX_CORNERS_OF_ELEM+1];
2494 CONNECTION *theCon;
2495 char msg[128];
2496 INT errors = 0;
2497
2498 sprintf(msg, " ERROR: missing connection between elem0=" EID_FMTX " elem1=" EID_FMTX,
2499 EID_PRTX(Elem0),EID_PRTX(Elem1));
2500
2501 cnt0 = GetAllVectorsOfElement(theGrid,Elem0,vec0);
2502 if (Elem0 == Elem1)
2503 {
2504 for (i=0; i<cnt0; i++)
2505 {
2506 itype=VTYPE(vec0[i]);
2507 for (j=0; j<cnt0; j++)
2508 {
2509 if (i==j)
2510 {
2511 mtype = DIAGMATRIXTYPE(itype);
2512 size = MatSize[mtype];
2513 }
2514 else
2515 {
2516 jtype = VTYPE(vec0[j]);
2517 size = MatSize[MATRIXTYPE(jtype,itype)];
2518 mtype = MATRIXTYPE(itype,jtype);
2519 if (MatSize[mtype]>size) size=MatSize[mtype];
2520 }
2521 if (size > 0)
2522 if (ActDepth <= ConDepth[mtype])
2523 {
2524 /* check connection in both directions */
2525 theCon = GetConnection(vec0[i],vec0[j]);
2526 if (theCon==NULL)
2527 {
2528 errors++;
2529 UserWriteF("%s vec0[%d]=" VINDEX_FMTX
2530 " to vec0[%d]=" VINDEX_FMTX "\n",
2531 msg,i,VINDEX_PRTX(vec0[i]),
2532 j,VINDEX_PRTX(vec0[j]));
2533 }
2534 else
2535 {
2536 theCon = GetConnection(vec0[j],vec0[i]);
2537 if (theCon==NULL)
2538 {
2539 errors++;
2540 UserWriteF("%s vec0[%d]=" VINDEX_FMTX
2541 " to vec0[%d]=" VINDEX_FMTX "\n",
2542 msg,j,VINDEX_PRTX(vec0[j]),
2543 i,VINDEX_PRTX(vec0[i]));
2544 }
2545 else
2546 {
2547 SETCUSED(theCon,1);
2548 }
2549 }
2550 }
2551 }
2552 }
2553 return (errors);
2554 }
2555
2556 cnt1 = GetAllVectorsOfElement(theGrid,Elem1,vec1);
2557 for (i=0; i<cnt0; i++)
2558 {
2559 itype=VTYPE(vec0[i]);
2560 for (j=0; j<cnt1; j++)
2561 {
2562 if (i==j)
2563 {
2564 mtype = DIAGMATRIXTYPE(itype);
2565 size = MatSize[mtype];
2566 }
2567 else
2568 {
2569 jtype = VTYPE(vec1[j]);
2570 size = MatSize[MATRIXTYPE(jtype,itype)];
2571 mtype = MATRIXTYPE(itype,jtype);
2572 if (MatSize[mtype]>size) size=MatSize[mtype];
2573 }
2574 if (size > 0)
2575 if (ActDepth <= ConDepth[mtype])
2576 {
2577 /* check connection in both directions */
2578 theCon = GetConnection(vec0[i],vec1[j]);
2579 if (theCon==NULL)
2580 {
2581 errors++;
2582 UserWriteF("%s vec0[%d]=" VINDEX_FMTX
2583 " to vec1[%d]=" VINDEX_FMTX "\n",
2584 msg,i,VINDEX_PRTX(vec0[i]),
2585 j,VINDEX_PRTX(vec1[j]));
2586 }
2587 else
2588 {
2589 theCon = GetConnection(vec1[j],vec0[i]);
2590 if (theCon == NULL)
2591 {
2592 errors++;
2593 UserWriteF("%s vec1[%d]=" VINDEX_FMTX
2594 " to vec0[%d]=%x/" VINDEX_FMTX "\n",
2595 msg,j,VINDEX_PRTX(vec1[j]),
2596 i,VINDEX_PRTX(vec0[i]));
2597 }
2598 else
2599 SETCUSED(theCon,1);
2600 }
2601 }
2602 }
2603 }
2604
2605 return (errors);
2606 }
2607
CheckNeighborhood(GRID * theGrid,ELEMENT * theElement,ELEMENT * centerElement,INT * ConDepth,INT ActDepth,INT MaxDepth,INT * MatSize)2608 static INT CheckNeighborhood (GRID *theGrid, ELEMENT *theElement, ELEMENT *centerElement, INT *ConDepth, INT ActDepth, INT MaxDepth, INT *MatSize)
2609 {
2610 INT i,errors = 0;
2611
2612 /* is anything to do ? */
2613 if (theElement==NULL) return (0);
2614
2615 /* action */
2616 if (ActDepth>=0)
2617 if ((errors+=ElementElementCheck(theGrid,centerElement,theElement,ActDepth,ConDepth,MatSize))!=0)
2618 return (errors);
2619
2620 /* call all neighbors recursively */
2621 if (ActDepth<MaxDepth)
2622 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
2623 if ((errors+=CheckNeighborhood(theGrid,NBELEM(theElement,i),centerElement,ConDepth,ActDepth+1,MaxDepth,MatSize))!=0)
2624 return(errors);
2625
2626 return (0);
2627 }
2628
2629 /****************************************************************************/
2630 /** \brief Check connection of the element
2631
2632 * @param theGrid - pointer to grid
2633 * @param theElement - pointer to element
2634
2635 This function checks all connections of the given element.
2636
2637 * @return <ul>
2638 * <li> number of errors
2639 * <li> 0 if ok
2640 * <li> != 0 if errors occured in that connection.
2641 </ul>
2642 */
2643 /****************************************************************************/
2644
ElementCheckConnection(GRID * theGrid,ELEMENT * theElement)2645 INT NS_DIM_PREFIX ElementCheckConnection (GRID *theGrid, ELEMENT *theElement)
2646 {
2647 INT MaxDepth;
2648 INT *ConDepth;
2649 INT *MatSize;
2650
2651 /* set pointers */
2652 FORMAT* theFormat = theGrid->mg->theFormat.get();
2653 MaxDepth = FMT_CONN_DEPTH_MAX(theFormat);
2654 ConDepth = FMT_CONN_DEPTH_PTR(theFormat);
2655 MatSize = FMT_S_MATPTR(theFormat);
2656
2657 /* call elements recursivly */
2658 return(CheckNeighborhood(theGrid,theElement,theElement,ConDepth,0,MaxDepth,MatSize));
2659 }
2660
2661 /****************************************************************************/
2662 /** \brief Check if connections are correctly allocated
2663
2664 * @param theGrid - grid level to check
2665
2666 This function checks if connections are correctly allocated.
2667
2668 * @return <ul>
2669 * <li> GM_OK if ok
2670 * <li> GM_ERROR if error occured.
2671 </ul>
2672 */
2673 /****************************************************************************/
2674
CheckConnections(GRID * theGrid)2675 static INT CheckConnections (GRID *theGrid)
2676 {
2677 ELEMENT *theElement;
2678 INT errors;
2679 INT error;
2680
2681 errors = 0;
2682
2683 for (theElement=FIRSTELEMENT(theGrid);
2684 theElement!=NULL;
2685 theElement=SUCCE(theElement))
2686 {
2687 if ((error=ElementCheckConnection(theGrid,theElement))!=0)
2688 {
2689 UserWriteF("element=" EID_FMTX " has bad connections\n",
2690 EID_PRTX(theElement));
2691 errors+=error;
2692 }
2693 }
2694
2695 return(errors);
2696 }
2697
2698
2699 /****************************************************************************/
2700 /** \brief Checks validity of geom_object and its vector
2701
2702 * @param fmt - FORMAT of associated multigrid
2703 * @param s2p - subdomain to part table
2704 * @param theObject - the object which points to theVector
2705 * @param ObjectString - for error message
2706 * @param theVector - the vector of theObject
2707 * @param VectorObjType - NODEVEC,...
2708 * @param side - element side for SIDEVEC, NOSIDE else
2709
2710 This function checks the consistency between an geom_object and
2711 its vector.
2712
2713 * @return <ul>
2714 * <li> GM_OK if ok
2715 * <li> GM_ERROR if error occured.
2716 </ul>
2717 */
2718 /****************************************************************************/
2719
CheckVector(const FORMAT * fmt,const INT s2p[],GEOM_OBJECT * theObject,const char * ObjectString,VECTOR * theVector,INT VectorObjType,INT side)2720 static INT CheckVector (const FORMAT *fmt, const INT s2p[], GEOM_OBJECT *theObject, const char *ObjectString,
2721 VECTOR *theVector, INT VectorObjType, INT side)
2722 {
2723 GEOM_OBJECT *VecObject;
2724 MATRIX *mat;
2725 INT errors = 0,vtype,DomPart,ds;
2726
2727 if (theVector == NULL)
2728 {
2729 /* check if size is really 0 */
2730 DomPart = GetDomainPart(s2p,theObject,side);
2731 vtype = FMT_PO2T(fmt,DomPart,VectorObjType);
2732 ds = FMT_S_VEC_TP(fmt,vtype);
2733 if (ds>0)
2734 {
2735 errors++;
2736 UserWriteF("%s ID=%ld has NO VECTOR", ObjectString,
2737 ID(theObject));
2738 #ifdef ModelP
2739 #ifdef __THREEDIM__
2740 if (VectorObjType == EDGEVEC)
2741 UserWriteF(" prio=%d",PRIO((EDGE*)theObject));
2742 #endif
2743 #endif
2744 UserWrite("\n");
2745 }
2746 }
2747 else
2748 {
2749 ds = FMT_S_VEC_TP(fmt,VTYPE(theVector));
2750 if (ds==0)
2751 {
2752 errors++;
2753 UserWriteF("%s ID=%ld exists but should not\n", ObjectString,
2754 ID(theObject));
2755 }
2756 SETVCUSED(theVector,1);
2757
2758 VecObject = VOBJECT(theVector);
2759 if (VecObject == NULL)
2760 {
2761 errors++;
2762 UserWriteF("vector=" VINDEX_FMTX " %s GID=" GID_FMT " has NO BACKPTR\n",
2763 VINDEX_PRTX(theVector), ObjectString,
2764 (OBJT(theObject)==BEOBJ || OBJT(theObject)==IEOBJ) ?
2765 EGID(&(theObject->el)) : (OBJT(theObject)==NDOBJ) ?
2766 GID(&(theObject->nd)) :
2767 GID(&(theObject->ed))
2768 );
2769 }
2770 else
2771 {
2772 if (VOTYPE(theVector) != VectorObjType)
2773 {
2774 errors++;
2775 UserWriteF("%s vector=" VINDEX_FMTX " has incompatible type=%d, "
2776 "should be type=%s\n",
2777 ObjectString, VINDEX_PRTX(theVector), VTYPE(theVector),
2778 ObjTypeName[VectorObjType]);
2779 }
2780
2781 if (VecObject != theObject)
2782 {
2783 if (OBJT(VecObject) != OBJT(theObject))
2784 {
2785 int error = 1;
2786
2787 /* both objects may be elements */
2788 if ((OBJT(VecObject)==BEOBJ || OBJT(VecObject)==IEOBJ) &&
2789 (OBJT(theObject)==BEOBJ || OBJT(theObject)==IEOBJ) )
2790 {
2791 ELEMENT *theElement = (ELEMENT *)theObject;
2792 ELEMENT *vecElement = (ELEMENT *)VecObject;
2793 int i;
2794
2795 #ifdef ModelP
2796 if (EMASTER(theElement) ||
2797 EMASTER(vecElement) )
2798 {
2799 #endif
2800 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
2801 if (NBELEM(theElement,i) == vecElement)
2802 {
2803 /* they are neighbors -> ok */
2804 error = 0;
2805 break;
2806 }
2807 #ifdef ModelP
2808 }
2809 #endif
2810 if (error)
2811 {
2812 UserWriteF("vector=" VINDEX_FMTX " has type %s, but points "
2813 "to wrong vecobj=" EID_FMTX " NO NB of obj=" EID_FMTX "\n",
2814 VINDEX_PRTX(theVector),ObjectString,
2815 EID_PRTX(vecElement),EID_PRTX(theElement));
2816 }
2817 }
2818 else
2819 {
2820 errors++;
2821 UserWriteF("vector=" VINDEX_FMTX " has type %s, but points "
2822 "to wrong obj=%d type OBJT=%d\n",
2823 VINDEX_PRTX(theVector),ObjectString,ID(VecObject),
2824 OBJT(VecObject));
2825 }
2826 }
2827 else
2828 {
2829 #ifdef __THREEDIM__
2830 if (VectorObjType == SIDEVEC)
2831 {
2832 /* TODO: check side vector */
2833 }
2834 else
2835 #endif
2836 {
2837 errors++;
2838 UserWriteF("%s vector=" VINDEX_FMTX " is referenced by "
2839 "obj0=%x, but points to wrong obj1=%x\n",
2840 ObjectString, VINDEX_PRTX(theVector),
2841 theObject, VecObject);
2842 #ifdef ModelP
2843 if (strcmp(ObjectString,"EDGE")==0)
2844 UserWriteF("obj0: n0=%d n1=%d obj1: "
2845 "n0=%d n1=%d\n",
2846 GID(NBNODE(LINK0(&(theObject->ed)))),
2847 GID(NBNODE(LINK1(&(theObject->ed)))),
2848 GID(NBNODE(LINK0(&(VecObject->ed)))),
2849 GID(NBNODE(LINK1(&(VecObject->ed)))) );
2850 #endif
2851 }
2852 }
2853 }
2854 }
2855 /* check connectivity of matrices */
2856 for (mat=VSTART(theVector); mat!=NULL; mat=MNEXT(mat))
2857 if (MDEST(mat)==NULL)
2858 {
2859 errors++;
2860 UserWriteF("%s vector=" VINDEX_FMTX ": matrix dest==NULL\n",
2861 ObjectString, VINDEX_PRTX(theVector));
2862 }
2863 else if (MDEST(MADJ(mat))!=theVector)
2864 {
2865 errors++;
2866 UserWriteF("%s vector=" VINDEX_FMTX ": adj matrix dest does not coincide"
2867 " with vector conn=%x mat=%x mdest=%x\n",
2868 ObjectString, VINDEX_PRTX(theVector),MMYCON(mat),MDEST(mat),MDEST(MADJ(mat)));
2869 }
2870
2871 }
2872
2873 return(errors);
2874 }
2875
2876 /****************************************************************************/
2877 /** \brief Check the algebraic part of the data structure
2878
2879 * @param theGrid - grid level to check
2880
2881 This function checks the consistency of the algebraic data structures
2882 including the interconnections between the geometric part.
2883 This function assumes a correct geometric data structure.
2884
2885 * @return <ul>
2886 * <li> GM_OK if ok </li>
2887 * <li> GM_ERROR if error occured. </li>
2888 </ul>
2889 */
2890 /****************************************************************************/
2891
CheckAlgebra(GRID * theGrid)2892 INT NS_DIM_PREFIX CheckAlgebra (GRID *theGrid)
2893 {
2894 ELEMENT *theElement;
2895 NODE *theNode;
2896 VECTOR *theVector;
2897 EDGE *theEdge;
2898 LINK *theLink;
2899 INT errors;
2900 INT *s2p;
2901 # ifdef __THREEDIM__
2902 INT i;
2903 # endif
2904
2905 errors = 0;
2906
2907 if ((GLEVEL(theGrid)==0) && !MG_COARSE_FIXED(MYMG(theGrid)))
2908 {
2909 if ((NVEC(theGrid)>0) || (NC(theGrid)>0))
2910 {
2911 errors++;
2912 UserWriteF("coarse grid not fixed but vectors allocated\n");
2913 }
2914 return(errors);
2915 }
2916
2917 s2p = BVPD_S2P_PTR(MG_BVPD(MYMG(theGrid)));
2918 FORMAT* fmt = theGrid->mg->theFormat.get();
2919
2920 /* reset USED flag */
2921 for (theVector=PFIRSTVECTOR(theGrid); theVector!=NULL;
2922 theVector=SUCCVC(theVector))
2923 {
2924 SETVCUSED(theVector,0);
2925 }
2926
2927 /* check pointers to element, side, edge vector */
2928 for (theElement=PFIRSTELEMENT(theGrid); theElement!=NULL;
2929 theElement=SUCCE(theElement))
2930 {
2931
2932 /* check element vectors */
2933 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,ELEMVEC))
2934 {
2935 theVector = EVECTOR(theElement);
2936 errors += CheckVector(fmt,s2p,(GEOM_OBJECT *) theElement, "ELEMENT",
2937 theVector, ELEMVEC,NOSIDE);
2938 }
2939
2940 #ifdef __THREEDIM__
2941 /* check side vectors */
2942 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,SIDEVEC))
2943 {
2944 for (i=0; i<SIDES_OF_ELEM(theElement); i++)
2945 {
2946 theVector = SVECTOR(theElement,i);
2947 errors += CheckVector(fmt,s2p,(GEOM_OBJECT *) theElement, "ELEMSIDE",
2948 theVector, SIDEVEC,i);
2949 }
2950 }
2951 #endif
2952 }
2953
2954 for (theNode=PFIRSTNODE(theGrid); theNode!=NULL; theNode=SUCCN(theNode))
2955 {
2956
2957 /* check node vectors */
2958 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,NODEVEC))
2959 {
2960 theVector = NVECTOR(theNode);
2961 errors += CheckVector(fmt,s2p,(GEOM_OBJECT *) theNode, "NODE", theVector,
2962 NODEVEC,NOSIDE);
2963 }
2964
2965 /* check edge vectors */
2966 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,EDGEVEC))
2967 {
2968 for (theLink=START(theNode); theLink!=NULL; theLink=NEXT(theLink))
2969 {
2970 theEdge = GetEdge(theNode,NBNODE(theLink));
2971 if (theEdge != NULL) {
2972 theVector = EDVECTOR(theEdge);
2973 errors += CheckVector(fmt,s2p,(GEOM_OBJECT *) theEdge, "EDGE",
2974 theVector, EDGEVEC,NOSIDE);
2975 }
2976 }
2977 }
2978 }
2979
2980 /* check USED flag */
2981 for (theVector=PFIRSTVECTOR(theGrid); theVector!=NULL;
2982 theVector=SUCCVC(theVector))
2983 {
2984 if (VCUSED(theVector) != 1)
2985 {
2986 errors++;
2987 UserWriteF("vector" VINDEX_FMTX " NOT referenced by an geom_object: "
2988 "vtype=%d, objptr=%x",
2989 VINDEX_PRTX(theVector), VTYPE(theVector), VOBJECT(theVector));
2990 if (VOBJECT(theVector) != NULL)
2991 UserWriteF(" objtype=%d\n",OBJT(VOBJECT(theVector)));
2992 else
2993 UserWrite("\n");
2994 }
2995 else
2996 SETVCUSED(theVector,0);
2997 }
2998
2999 /* check validity of all defined connections */
3000 errors += CheckConnections(theGrid);
3001
3002 /* reset flags in connections */
3003 for (theVector=PFIRSTVECTOR(theGrid); theVector!=NULL;
3004 theVector=SUCCVC(theVector))
3005 {
3006 MATRIX *theMatrix;
3007
3008 for (theMatrix=VSTART(theVector); theMatrix!=NULL;
3009 theMatrix = MNEXT(theMatrix)) SETCUSED(MMYCON(theMatrix),0);
3010 }
3011
3012 /* set flags in connections */
3013 #if defined __OVERLAP2__
3014 for (theVector=PFIRSTVECTOR(theGrid); theVector!=NULL; theVector=SUCCVC(theVector))
3015 #else
3016 for (theVector=FIRSTVECTOR(theGrid); theVector!=NULL; theVector=SUCCVC(theVector))
3017 #endif
3018 {
3019 MATRIX *theMatrix;
3020
3021 for (theMatrix=VSTART(theVector); theMatrix!=NULL;
3022 theMatrix = MNEXT(theMatrix)) SETMUSED(MADJ(theMatrix),1);
3023 }
3024
3025 /* check matrices and connections */
3026 for (theVector=PFIRSTVECTOR(theGrid);
3027 theVector!=NULL;
3028 theVector=SUCCVC(theVector))
3029 {
3030 MATRIX *theMatrix;
3031 #ifdef ModelP
3032 INT prio = PRIO(theVector);
3033 #endif
3034
3035 for (theMatrix=VSTART(theVector);
3036 theMatrix!=NULL;
3037 theMatrix = MNEXT(theMatrix))
3038 {
3039 MATRIX *Adj = MADJ(theMatrix);
3040
3041 if (MDEST(theMatrix) == NULL)
3042 {
3043 errors++;
3044 UserWriteF("ERROR: matrix %x has no dest, start vec="
3045 VINDEX_FMTX "\n",
3046 theMatrix,VINDEX_PRTX(theVector));
3047 }
3048 if (MDEST(Adj) != theVector)
3049 {
3050 errors++;
3051 UserWriteF("ERROR: dest=%x of adj matrix "
3052 " unequal vec=" VINDEX_FMTX "\n",
3053 MDEST(Adj),VINDEX_PRTX(theVector));
3054 }
3055
3056 #if defined ModelP && ! defined __OVERLAP2__
3057 if (prio != PrioHGhost)
3058 #endif
3059 if (MUSED(theMatrix)!=1 && !CEXTRA(MMYCON(theMatrix)))
3060 {
3061 errors++;
3062 UserWriteF("ERROR: connection dead vec=" VINDEX_FMTX
3063 " vector=" VINDEX_FMTX " con=%x mat=%x matadj=%x level(vec)=%d is_extra_con %d\n",
3064 VINDEX_PRTX(theVector),VINDEX_PRTX(MDEST(theMatrix)),
3065 MMYCON(theMatrix),MDEST(theMatrix),MDEST(MADJ(theMatrix)),
3066 GLEVEL(theGrid),CEXTRA(MMYCON(theMatrix)));
3067 }
3068
3069 #if defined ModelP && ! defined __OVERLAP2__
3070 if (GHOSTPRIO(prio) && !CEXTRA(MMYCON(theMatrix)))
3071 {
3072 errors++;
3073 UserWriteF("ERROR: ghost vector has matrix vec="
3074 VINDEX_FMTX " con=%x mat=%x\n",
3075 VINDEX_PRTX(theVector),MMYCON(theMatrix),theMatrix);
3076 }
3077 #endif
3078 }
3079 }
3080
3081 return(errors);
3082 }
3083
3084 /****************************************************************************/
3085 /** \brief Decide whether a vector corresponds to an element or not
3086
3087 * @param theElement - pointer to element
3088 * @param theVector - pointer to a vector
3089
3090 This function decides whether a given vector belongs to the given element, or
3091 one of its sides, edges or nodes.
3092
3093 * @return <ul>
3094 * <li> 0 if does not correspond </li>
3095 * <li> 1 if does correspond. </li>
3096 </ul>
3097 */
3098 /****************************************************************************/
3099
VectorInElement(ELEMENT * theElement,VECTOR * theVector)3100 INT NS_DIM_PREFIX VectorInElement (ELEMENT *theElement, VECTOR *theVector)
3101 {
3102 INT i;
3103 VECTOR *vList[20];
3104 INT cnt;
3105
3106 if (VOTYPE(theVector) == ELEMVEC)
3107 {
3108 GetVectorsOfElement(theElement,&cnt,vList);
3109 for (i=0; i<cnt; i++)
3110 if (vList[i]==theVector) RETURN(1);
3111 }
3112 #ifdef __THREEDIM__
3113 if (VOTYPE(theVector) == SIDEVEC)
3114 {
3115 GetVectorsOfSides(theElement,&cnt,vList);
3116 for (i=0; i<cnt; i++)
3117 if (vList[i]==theVector) RETURN(1);
3118 }
3119 #endif
3120 if (VOTYPE(theVector) == EDGEVEC)
3121 {
3122 GetVectorsOfEdges(theElement,&cnt,vList);
3123 for (i=0; i<cnt; i++)
3124 if (vList[i]==theVector) RETURN(1);
3125 }
3126 if (VOTYPE(theVector) == NODEVEC)
3127 {
3128 GetVectorsOfNodes(theElement,&cnt,vList);
3129 for (i=0; i<cnt; i++)
3130 if (vList[i]==theVector) RETURN(1);
3131 }
3132
3133 return (0);
3134 }
3135
3136 /****************************************************************************/
3137 /** \brief Calc coordinate position of vector
3138
3139 * @param theVector - a given vector
3140 * @param position - array to be filled
3141
3142 This function calcs physical position of vector. For edge vectors the
3143 midpoint is returned, and for sides and elements the center of mass.
3144
3145 * @return <ul>
3146 * <li> 0 if ok </li>
3147 * <li> 1 if error occured. </li>
3148 </ul>
3149 */
3150 /****************************************************************************/
3151
VectorPosition(const VECTOR * theVector,DOUBLE * position)3152 INT NS_DIM_PREFIX VectorPosition (const VECTOR *theVector, DOUBLE *position)
3153 {
3154 INT i;
3155 EDGE *theEdge;
3156 #ifdef __THREEDIM__
3157 ELEMENT *theElement;
3158 INT theSide,j;
3159 #endif
3160
3161 ASSERT(theVector != NULL);
3162
3163 #if defined __OVERLAP2__
3164 if( VOBJECT(theVector) == NULL )
3165 {
3166 for (i=0; i<DIM; i++)
3167 position[i] = -MAX_D;
3168 return (0);
3169 }
3170 #endif
3171
3172 switch (VOTYPE(theVector))
3173 {
3174 case NODEVEC :
3175 for (i=0; i<DIM; i++)
3176 position[i] = CVECT(MYVERTEX((NODE*)VOBJECT(theVector)))[i];
3177 return (0);
3178
3179 case EDGEVEC :
3180 theEdge = (EDGE*)VOBJECT(theVector);
3181 for (i=0; i<DIM; i++)
3182 position[i] = 0.5*(CVECT(MYVERTEX(NBNODE(LINK0(theEdge))))[i] +
3183 CVECT(MYVERTEX(NBNODE(LINK1(theEdge))))[i] );
3184 return (0);
3185 #ifdef __THREEDIM__
3186 case SIDEVEC :
3187 theElement = (ELEMENT *)VOBJECT(theVector);
3188 theSide = VECTORSIDE(theVector);
3189 for (i=0; i<DIM; i++)
3190 {
3191 position[i] = 0.0;
3192 for(j=0; j<CORNERS_OF_SIDE(theElement,theSide); j++)
3193 position[i] += CVECT(MYVERTEX(CORNER(theElement,CORNER_OF_SIDE(theElement,theSide,j))))[i];
3194 position[i] /= CORNERS_OF_SIDE(theElement,theSide);
3195 }
3196 return (0);
3197 #endif
3198 case ELEMVEC :
3199 /* calculate center of mass */
3200 CalculateCenterOfMass( (ELEMENT *) VOBJECT(theVector), position );
3201 return (0);
3202
3203 default : PrintErrorMessage('E',"VectorPosition","unrecognized object type for vector");
3204 assert(0);
3205 }
3206
3207 RETURN (GM_ERROR);
3208 }
3209
3210
3211 /****************************************************************************/
3212 /** \brief Initialize vector classes
3213
3214 * @param theGrid - given grid
3215 * @param theElement - given element
3216
3217 Initialize vector class in all vectors associated with given element with 3.
3218
3219 * @return <ul>
3220 * <li> 0 if ok </li>
3221 * <li> 1 if error occured. </li>
3222 </ul>
3223 */
3224 /****************************************************************************/
3225
SeedVectorClasses(GRID * theGrid,ELEMENT * theElement)3226 INT NS_DIM_PREFIX SeedVectorClasses (GRID *theGrid, ELEMENT *theElement)
3227 {
3228 INT i;
3229 VECTOR *vList[20];
3230 INT cnt;
3231
3232 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,ELEMVEC))
3233 {
3234 GetVectorsOfElement(theElement,&cnt,vList);
3235 for (i=0; i<cnt; i++) SETVCLASS(vList[i],3);
3236 }
3237 #ifdef __THREEDIM__
3238 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,SIDEVEC))
3239 {
3240 GetVectorsOfSides(theElement,&cnt,vList);
3241 for (i=0; i<cnt; i++) SETVCLASS(vList[i],3);
3242 }
3243 #endif
3244 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,EDGEVEC))
3245 {
3246 GetVectorsOfEdges(theElement,&cnt,vList);
3247 for (i=0; i<cnt; i++) SETVCLASS(vList[i],3);
3248 }
3249 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,NODEVEC))
3250 {
3251 GetVectorsOfNodes(theElement,&cnt,vList);
3252 for (i=0; i<cnt; i++) SETVCLASS(vList[i],3);
3253 }
3254 return (0);
3255 }
3256
3257 /****************************************************************************/
3258 /** \brief Reset vector classes
3259
3260 * @param theGrid - pointer to grid
3261
3262 Reset all vector classes in all vectors of given grid to 0.
3263
3264 * @return <ul>
3265 * <li> 0 if ok </li>
3266 * <li> 1 if error occured. </li>
3267 </ul>
3268 */
3269 /****************************************************************************/
3270
ClearVectorClasses(GRID * theGrid)3271 INT NS_DIM_PREFIX ClearVectorClasses (GRID *theGrid)
3272 {
3273 VECTOR *theVector;
3274
3275 /* reset class of each vector to 0 */
3276 for (theVector=PFIRSTVECTOR(theGrid); theVector!=NULL; theVector=SUCCVC(theVector))
3277 SETVCLASS(theVector,0);
3278
3279 return(0);
3280 }
3281
3282
3283 #ifdef ModelP
Gather_VectorVClass(DDD::DDDContext &,DDD_OBJ obj,void * data)3284 static int Gather_VectorVClass (DDD::DDDContext&, DDD_OBJ obj, void *data)
3285 {
3286 VECTOR *theVector = (VECTOR *)obj;
3287
3288 PRINTDEBUG(gm,1,("Gather_VectorVClass(): v=" VINDEX_FMTX " vclass=%d\n",
3289 VINDEX_PRTX(theVector),VCLASS(theVector)))
3290
3291 ((INT *)data)[0] = VCLASS(theVector);
3292
3293 return(0);
3294 }
3295
Scatter_VectorVClass(DDD::DDDContext &,DDD_OBJ obj,void * data)3296 static int Scatter_VectorVClass (DDD::DDDContext&, DDD_OBJ obj, void *data)
3297 {
3298 VECTOR *theVector = (VECTOR *)obj;
3299
3300 SETVCLASS(theVector,MAX(VCLASS(theVector),((INT *)data)[0]));
3301
3302 PRINTDEBUG(gm,2,("Scatter_VectorVClass(): v=" VINDEX_FMTX " vclass=%d\n",
3303 VINDEX_PRTX(theVector),VCLASS(theVector)))
3304
3305 return(0);
3306 }
3307
Scatter_GhostVectorVClass(DDD::DDDContext &,DDD_OBJ obj,void * data)3308 static int Scatter_GhostVectorVClass (DDD::DDDContext&, DDD_OBJ obj, void *data)
3309 {
3310 VECTOR *theVector = (VECTOR *)obj;
3311
3312 SETVCLASS(theVector,((INT *)data)[0]);
3313
3314 PRINTDEBUG(gm,1,("Scatter_GhostVectorVClass(): v=" VINDEX_FMTX " vclass=%d\n",
3315 VINDEX_PRTX(theVector),VCLASS(theVector)))
3316
3317 return(0);
3318 }
3319 #endif
3320
PropagateVectorClass(GRID * theGrid,INT vclass)3321 static INT PropagateVectorClass (GRID *theGrid, INT vclass)
3322 {
3323 VECTOR *theVector;
3324 MATRIX *theMatrix;
3325
3326 /* set vector classes in the algebraic neighborhood to vclass-1 */
3327 /* use matrices to determine next vectors!!!!! */
3328 for (theVector=FIRSTVECTOR(theGrid); theVector!=NULL;
3329 theVector=SUCCVC(theVector))
3330 if ((VCLASS(theVector)==vclass)&&(VSTART(theVector)!=NULL))
3331 for (theMatrix=MNEXT(VSTART(theVector)); theMatrix!=NULL;
3332 theMatrix=MNEXT(theMatrix))
3333 if ((VCLASS(MDEST(theMatrix))<vclass)
3334 &&(CEXTRA(MMYCON(theMatrix))!=1))
3335 SETVCLASS(MDEST(theMatrix),vclass-1);
3336
3337 /* only for this values valid */
3338 ASSERT(vclass==3 || vclass==2);
3339
3340 return(0);
3341 }
3342
3343
3344 /****************************************************************************/
3345 /** \brief Compute vector classes after initialization
3346
3347 * @param theGrid - pointer to grid
3348
3349 After vector classes have been reset and initialized, this function
3350 now computes the class 2 and class 1 vectors.
3351
3352 * @return <ul>
3353 * <li> 0 if ok </li>
3354 * <li> 1 if error occured </li>
3355 </ul>
3356 */
3357 /****************************************************************************/
PropagateVectorClasses(GRID * theGrid)3358 INT NS_DIM_PREFIX PropagateVectorClasses (GRID *theGrid)
3359 {
3360 #ifdef ModelP
3361 auto& context = theGrid->dddContext();
3362 const auto& dddctrl = ddd_ctrl(context);
3363
3364 PRINTDEBUG(gm,1,("\nPropagateVectorClasses():"
3365 " 1. communication on level %d\n",GLEVEL(theGrid)))
3366 /* exchange VCLASS of vectors */
3367 DDD_IFAExchange(context,
3368 dddctrl.BorderVectorSymmIF,GRID_ATTR(theGrid),sizeof(INT),
3369 Gather_VectorVClass, Scatter_VectorVClass);
3370 #endif
3371
3372 /* set vector classes in the algebraic neighborhood to 2 */
3373 if (PropagateVectorClass(theGrid,3)) REP_ERR_RETURN(1);
3374
3375 #ifdef ModelP
3376 PRINTDEBUG(gm,1,("\nPropagateVectorClasses(): 2. communication\n"))
3377 /* exchange VCLASS of vectors */
3378 DDD_IFAExchange(context,
3379 dddctrl.BorderVectorSymmIF,GRID_ATTR(theGrid),sizeof(INT),
3380 Gather_VectorVClass, Scatter_VectorVClass);
3381 #endif
3382
3383 /* set vector classes in the algebraic neighborhood to 1 */
3384 if (PropagateVectorClass(theGrid,2)) REP_ERR_RETURN(1);
3385
3386 #ifdef ModelP
3387 PRINTDEBUG(gm,1,("\nPropagateVectorClasses(): 3. communication\n"))
3388 /* exchange VCLASS of vectors */
3389 DDD_IFAExchange(context,
3390 dddctrl.BorderVectorSymmIF,GRID_ATTR(theGrid),sizeof(INT),
3391 Gather_VectorVClass, Scatter_VectorVClass);
3392 #endif
3393
3394 #ifdef ModelP
3395 /* send VCLASS to ghosts */
3396 DDD_IFAOneway(context,
3397 dddctrl.VectorIF,GRID_ATTR(theGrid),IF_FORWARD,sizeof(INT),
3398 Gather_VectorVClass, Scatter_GhostVectorVClass);
3399 #endif
3400
3401 return(0);
3402 }
3403
3404 /****************************************************************************/
3405 /** \brief Reset class of the vectors on the next level
3406
3407 * @param theGrid - pointer to grid
3408
3409 This function clears VNCLASS flag in all vectors. This is the first step to
3410 compute the class of the dofs on the *NEXT* level, which
3411 is also the basis for determining copies.
3412
3413 * @return <ul>
3414 * <li> 0 if ok </li>
3415 * <li> 1 if error occured. </li>
3416 </ul>
3417 */
3418 /****************************************************************************/
3419
ClearNextVectorClasses(GRID * theGrid)3420 INT NS_DIM_PREFIX ClearNextVectorClasses (GRID *theGrid)
3421 {
3422 VECTOR *theVector;
3423
3424 /* reset class of each vector to 0 */
3425 for (theVector=PFIRSTVECTOR(theGrid); theVector!=NULL; theVector=SUCCVC(theVector))
3426 SETVNCLASS(theVector,0);
3427
3428 /* now the refinement algorithm will initialize the class 3 vectors */
3429 /* on the *NEXT* level. */
3430 return(0);
3431 }
3432
3433 /****************************************************************************/
3434 /** \brief Set VNCLASS in all vectors associated with element
3435
3436 * @param theGrid - given grid
3437 * @param theElement - pointer to element
3438
3439 Set VNCLASS in all vectors associated with the element to 3.
3440
3441 * @return <ul>
3442 * <li> 0 if ok </li>
3443 * <li> 1 if error occured. </li>
3444 </ul>
3445 */
3446 /****************************************************************************/
3447
SeedNextVectorClasses(GRID * theGrid,ELEMENT * theElement)3448 INT NS_DIM_PREFIX SeedNextVectorClasses (GRID *theGrid, ELEMENT *theElement)
3449 {
3450 INT i;
3451 VECTOR *vList[20];
3452 INT cnt;
3453
3454 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,ELEMVEC))
3455 {
3456 GetVectorsOfElement(theElement,&cnt,vList);
3457 for (i=0; i<cnt; i++) SETVNCLASS(vList[i],3);
3458 }
3459 #ifdef __THREEDIM__
3460 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,SIDEVEC))
3461 {
3462 GetVectorsOfSides(theElement,&cnt,vList);
3463 for (i=0; i<cnt; i++) SETVNCLASS(vList[i],3);
3464 }
3465 #endif
3466 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,EDGEVEC))
3467 {
3468 GetVectorsOfEdges(theElement,&cnt,vList);
3469 for (i=0; i<cnt; i++) SETVNCLASS(vList[i],3);
3470 }
3471 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,NODEVEC))
3472 {
3473 GetVectorsOfNodes(theElement,&cnt,vList);
3474 for (i=0; i<cnt; i++) SETVNCLASS(vList[i],3);
3475 }
3476 return (0);
3477 }
3478
3479
3480
3481
3482 #ifdef ModelP
Gather_VectorVNClass(DDD::DDDContext &,DDD_OBJ obj,void * data)3483 static int Gather_VectorVNClass (DDD::DDDContext&, DDD_OBJ obj, void *data)
3484 {
3485 VECTOR *theVector = (VECTOR *)obj;
3486
3487 PRINTDEBUG(gm,2,("Gather_VectorVNClass(): v=" VINDEX_FMTX " vnclass=%d\n",
3488 VINDEX_PRTX(theVector),VNCLASS(theVector)))
3489
3490 ((INT *)data)[0] = VNCLASS(theVector);
3491
3492 return(GM_OK);
3493 }
3494
Scatter_VectorVNClass(DDD::DDDContext &,DDD_OBJ obj,void * data)3495 static int Scatter_VectorVNClass (DDD::DDDContext&, DDD_OBJ obj, void *data)
3496 {
3497 VECTOR *theVector = (VECTOR *)obj;
3498
3499 SETVNCLASS(theVector,MAX(VNCLASS(theVector),((INT *)data)[0]));
3500
3501 PRINTDEBUG(gm,2,("Scatter_VectorVNClass(): v=" VINDEX_FMTX " vnclass=%d\n",
3502 VINDEX_PRTX(theVector),VNCLASS(theVector)))
3503
3504 return(GM_OK);
3505 }
3506
Scatter_GhostVectorVNClass(DDD::DDDContext &,DDD_OBJ obj,void * data)3507 static int Scatter_GhostVectorVNClass (DDD::DDDContext&, DDD_OBJ obj, void *data)
3508 {
3509 VECTOR *theVector = (VECTOR *)obj;
3510
3511 SETVNCLASS(theVector,((INT *)data)[0]);
3512
3513 PRINTDEBUG(gm,2,("Scatter_GhostVectorVNClass(): v=" VINDEX_FMTX " vnclass=%d\n",
3514 VINDEX_PRTX(theVector),VNCLASS(theVector)))
3515
3516 return(GM_OK);
3517 }
3518 #endif
3519
PropagateNextVectorClass(GRID * theGrid,INT vnclass)3520 static INT PropagateNextVectorClass (GRID *theGrid, INT vnclass)
3521 {
3522 VECTOR *theVector;
3523 MATRIX *theMatrix;
3524
3525 /* set vector classes in the algebraic neighborhood to vnclass-1 */
3526 /* use matrices to determine next vectors!!!!! */
3527 for (theVector=FIRSTVECTOR(theGrid); theVector!=NULL;
3528 theVector=SUCCVC(theVector))
3529 if ((VNCLASS(theVector)==vnclass)&&(VSTART(theVector)!=NULL))
3530 for (theMatrix=MNEXT(VSTART(theVector)); theMatrix!=NULL;
3531 theMatrix=MNEXT(theMatrix))
3532 if ((VNCLASS(MDEST(theMatrix))<vnclass)
3533 &&(CEXTRA(MMYCON(theMatrix))!=1))
3534 SETVNCLASS(MDEST(theMatrix),vnclass-1);
3535
3536 /* only for this values valid */
3537 ASSERT(vnclass==3 || vnclass==2);
3538
3539 return(0);
3540 }
3541
3542 /****************************************************************************/
3543 /** \brief Compute VNCLASS in all vectors of a grid level
3544
3545 * @param theGrid - pointer to grid
3546
3547 Computes values of VNCLASS field in all vectors after seed.
3548
3549 * @return <ul>
3550 * <li> 0 if ok </li>
3551 * <li> 1 if error occured </li>
3552 </ul>
3553 */
3554 /****************************************************************************/
PropagateNextVectorClasses(GRID * theGrid)3555 INT NS_DIM_PREFIX PropagateNextVectorClasses (GRID *theGrid)
3556 {
3557 #ifdef ModelP
3558 auto& context = theGrid->dddContext();
3559 const auto& dddctrl = ddd_ctrl(context);
3560
3561 PRINTDEBUG(gm,1,("\nPropagateNextVectorClasses(): 1. communication\n"))
3562 /* exchange VNCLASS of vectors */
3563 DDD_IFAExchange(context,
3564 dddctrl.BorderVectorSymmIF,GRID_ATTR(theGrid),sizeof(INT),
3565 Gather_VectorVNClass, Scatter_VectorVNClass);
3566 #endif
3567
3568 if (PropagateNextVectorClass(theGrid,3)) REP_ERR_RETURN(1);
3569
3570 #ifdef ModelP
3571 PRINTDEBUG(gm,1,("\nPropagateNextVectorClasses(): 2. communication\n"))
3572 /* exchange VNCLASS of vectors */
3573 DDD_IFAExchange(context,
3574 dddctrl.BorderVectorSymmIF,GRID_ATTR(theGrid),sizeof(INT),
3575 Gather_VectorVNClass, Scatter_VectorVNClass);
3576 #endif
3577
3578 if (PropagateNextVectorClass(theGrid,2)) REP_ERR_RETURN(1);
3579
3580 #ifdef ModelP
3581 PRINTDEBUG(gm,1,("\nPropagateNextVectorClasses(): 3. communication\n"))
3582 /* exchange VNCLASS of vectors */
3583 DDD_IFAExchange(context,
3584 dddctrl.BorderVectorSymmIF,GRID_ATTR(theGrid),sizeof(INT),
3585 Gather_VectorVNClass, Scatter_VectorVNClass);
3586
3587 /* send VCLASS to ghosts */
3588 DDD_IFAOneway(context,
3589 dddctrl.VectorIF,GRID_ATTR(theGrid),IF_FORWARD,sizeof(INT),
3590 Gather_VectorVNClass, Scatter_GhostVectorVNClass);
3591 #endif
3592
3593 return(0);
3594 }
3595
3596 /****************************************************************************/
3597 /** \brief Returns highest vector class of a dof on next level
3598
3599 * @param theGrid - pointer to a grid
3600 * @param theElement - pointer to a element
3601
3602 This function returns highest VNCLASS of a vector associated with the
3603 element.
3604
3605 * @return <ul>
3606 * <li> 0 if ok </li>
3607 * <li> 1 if error occured. </li>
3608 </ul>
3609 */
3610 /****************************************************************************/
3611
MaxNextVectorClass(GRID * theGrid,ELEMENT * theElement)3612 INT NS_DIM_PREFIX MaxNextVectorClass (GRID *theGrid, ELEMENT *theElement)
3613 {
3614 INT i,m;
3615 VECTOR *vList[20];
3616 INT cnt;
3617
3618 m = 0;
3619 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,ELEMVEC))
3620 {
3621 GetVectorsOfElement(theElement,&cnt,vList);
3622 for (i=0; i<cnt; i++) m = MAX(m,VNCLASS(vList[i]));
3623 }
3624 #ifdef __THREEDIM__
3625 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,SIDEVEC))
3626 {
3627 GetVectorsOfSides(theElement,&cnt,vList);
3628 for (i=0; i<cnt; i++) m = MAX(m,VNCLASS(vList[i]));
3629 }
3630 #endif
3631 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,EDGEVEC))
3632 {
3633 GetVectorsOfEdges(theElement,&cnt,vList);
3634 for (i=0; i<cnt; i++) m = MAX(m,VNCLASS(vList[i]));
3635 }
3636 if (VEC_DEF_IN_OBJ_OF_GRID(theGrid,NODEVEC))
3637 {
3638 GetVectorsOfNodes(theElement,&cnt,vList);
3639 for (i=0; i<cnt; i++) m = MAX(m,VNCLASS(vList[i]));
3640 }
3641 return (m);
3642 }
3643
3644 /****************************************************************************/
3645 /** \brief Create a new find cut procedure in environment
3646
3647 * @param name - name
3648 * @param FindCutProc - the find cut procedure
3649
3650 This function creates a new find cut dependency in environement.
3651
3652 */
3653 /****************************************************************************/
3654
CreateFindCutProc(const char * name,FindCutProcPtr FindCutProc)3655 FIND_CUT * NS_DIM_PREFIX CreateFindCutProc (const char *name, FindCutProcPtr FindCutProc)
3656 {
3657 FIND_CUT *newFindCut;
3658
3659 if (ChangeEnvDir("/FindCut")==NULL)
3660 {
3661 UserWrite("cannot change to dir '/FindCut'\n");
3662 return(NULL);
3663 }
3664 newFindCut = (FIND_CUT *) MakeEnvItem (name,theFindCutVarID,sizeof(FIND_CUT));
3665 if (newFindCut==NULL) return(NULL);
3666
3667 newFindCut->FindCutProc = FindCutProc;
3668
3669 return (newFindCut);
3670 }
3671
3672 /****************************************************************************/
3673 /** \brief Create a new algebraic dependency in environment
3674
3675 * @param name - name
3676 * @param DependencyProc - the dependency procedure
3677
3678 This function creates a new algebraic dependency in environement.
3679
3680 */
3681 /****************************************************************************/
3682
CreateAlgebraicDependency(const char * name,DependencyProcPtr DependencyProc)3683 ALG_DEP * NS_DIM_PREFIX CreateAlgebraicDependency (const char *name, DependencyProcPtr DependencyProc)
3684 {
3685 ALG_DEP *newAlgDep;
3686
3687 if (ChangeEnvDir("/Alg Dep")==NULL)
3688 {
3689 UserWrite("cannot change to dir '/Alg Dep'\n");
3690 return(NULL);
3691 }
3692 newAlgDep = (ALG_DEP *) MakeEnvItem (name,theAlgDepVarID,sizeof(ALG_DEP));
3693 if (newAlgDep==NULL) return(NULL);
3694
3695 newAlgDep->DependencyProc = DependencyProc;
3696
3697 return (newAlgDep);
3698 }
3699
3700 /****************************************************************************/
3701 /** \brief Put `every` vector !VCUSED into the new order list
3702
3703 * @param theGrid - pointer to grid
3704 . CutVectors -
3705 . nb -
3706
3707 This function can be used as a find-cut-procedure for the streamwise ordering
3708 algorithm. But it just puts `all` the remaining vectors into the new order
3709 list instead of removing only as much as necessary to get rid of cyclic
3710 dependencies.
3711
3712 It is used as default when no cut procedure is specified.
3713
3714 * @return <ul>
3715 * <li> 0 if ok
3716 * <li> 1 if error occured.
3717 </ul>
3718 */
3719 /****************************************************************************/
3720
FeedbackVertexVectors(GRID * theGrid,VECTOR * LastVector,INT * nb)3721 static VECTOR *FeedbackVertexVectors (GRID *theGrid, VECTOR *LastVector, INT *nb)
3722 {
3723 VECTOR *theVector;
3724
3725 /* push all remaining vectors */
3726 *nb = 0;
3727 for (theVector=FIRSTVECTOR(theGrid); theVector!=NULL; theVector=SUCCVC(theVector))
3728 if (!VCUSED(theVector))
3729 {
3730 (*nb)++;
3731 PREDVC(LastVector) = theVector;
3732 LastVector = theVector;
3733 SETVCUSED(theVector,1);
3734 }
3735 return (LastVector);
3736 }
3737
3738
3739 /****************************************************************************/
3740 /** \brief Dependency function for lexicographic ordering
3741
3742 * @param theGrid - pointer to grid
3743 * @param data - option string from orderv command
3744
3745 This function defines a dependency function for lexicographic ordering.
3746
3747 * @return <ul>
3748 * <li> 0 if ok
3749 * <li> 1 if error occured.
3750 </ul>
3751 */
3752 /****************************************************************************/
3753
LexAlgDep(GRID * theGrid,const char * data)3754 static INT LexAlgDep (GRID *theGrid, const char *data)
3755 {
3756 MULTIGRID *theMG;
3757 VECTOR *theVector,*NBVector;
3758 MATRIX *theMatrix;
3759 DOUBLE_VECTOR pos,nbpos;
3760 DOUBLE diff[DIM];
3761 INT i,order,res;
3762 INT Sign[DIM],Order[DIM],xused,yused,zused,error;
3763 char ord[3];
3764
3765 /* read ordering directions */
3766 #ifdef __TWODIM__
3767 res = sscanf(data,expandfmt("%2[rlud]"),ord);
3768 #else
3769 res = sscanf(data,expandfmt("%3[rlbfud]"),ord);
3770 #endif
3771 if (res!=1)
3772 {
3773 PrintErrorMessage('E',"LexAlgDep","could not read order type");
3774 return(1);
3775 }
3776 if (strlen(ord)!=DIM)
3777 {
3778 #ifdef __TWODIM__
3779 PrintErrorMessage('E',"LexAlgDep","specify 2 chars out of 'rlud'");
3780 #else
3781 PrintErrorMessage('E',"LexAlgDep","specify 3 chars out of 'rlbfud'");
3782 #endif
3783 return(1);
3784 }
3785 error = xused = yused = zused = false;
3786 for (i=0; i<DIM; i++)
3787 switch (ord[i])
3788 {
3789 case 'r' :
3790 if (xused) error = true;
3791 xused = true;
3792 Order[i] = _X_; Sign[i] = 1; break;
3793 case 'l' :
3794 if (xused) error = true;
3795 xused = true;
3796 Order[i] = _X_; Sign[i] = -1; break;
3797
3798 #ifdef __TWODIM__
3799 case 'u' :
3800 if (yused) error = true;
3801 yused = true;
3802 Order[i] = _Y_; Sign[i] = 1; break;
3803 case 'd' :
3804 if (yused) error = true;
3805 yused = true;
3806 Order[i] = _Y_; Sign[i] = -1; break;
3807 #else
3808 case 'b' :
3809 if (yused) error = true;
3810 yused = true;
3811 Order[i] = _Y_; Sign[i] = 1; break;
3812 case 'f' :
3813 if (yused) error = true;
3814 yused = true;
3815 Order[i] = _Y_; Sign[i] = -1; break;
3816
3817 case 'u' :
3818 if (zused) error = true;
3819 zused = true;
3820 Order[i] = _Z_; Sign[i] = 1; break;
3821 case 'd' :
3822 if (zused) error = true;
3823 zused = true;
3824 Order[i] = _Z_; Sign[i] = -1; break;
3825 #endif
3826 }
3827 if (error)
3828 {
3829 PrintErrorMessage('E',"LexAlgDep","bad combination of 'rludr' or 'rlbfud' resp.");
3830 return(1);
3831 }
3832
3833 theMG = MYMG(theGrid);
3834
3835 /* find an approximate measure for the mesh size */
3836 // The following method wants the domain radius, which has been removed.
3837 // Dune has been setting this radius to 1.0 for years now, so I don't think it matters.
3838 DOUBLE BVPD_RADIUS = 1.0;
3839 InvMeshSize = POW2(GLEVEL(theGrid)) * pow(NN(GRID_ON_LEVEL(theMG,0)),1.0/DIM) / BVPD_RADIUS;
3840
3841 for (theVector=FIRSTVECTOR(theGrid); theVector!=NULL; theVector=SUCCVC(theVector))
3842 {
3843 VectorPosition(theVector,pos);
3844
3845 for (theMatrix=MNEXT(VSTART(theVector)); theMatrix!=NULL; theMatrix=MNEXT(theMatrix))
3846 {
3847 NBVector = MDEST(theMatrix);
3848
3849 SETMUP(theMatrix,0);
3850 SETMDOWN(theMatrix,0);
3851
3852 VectorPosition(NBVector,nbpos);
3853
3854 V_DIM_SUBTRACT(nbpos,pos,diff);
3855 V_DIM_SCALE(InvMeshSize,diff);
3856
3857 if (fabs(diff[Order[DIM-1]])<ORDERRES)
3858 {
3859 #ifdef __THREEDIM__
3860 if (fabs(diff[Order[DIM-2]])<ORDERRES)
3861 {
3862 if (diff[Order[DIM-3]]>0.0) order = -Sign[DIM-3];
3863 else order = Sign[DIM-3];
3864 }
3865 else
3866 #endif
3867 if (diff[Order[DIM-2]]>0.0) order = -Sign[DIM-2];
3868 else order = Sign[DIM-2];
3869 }
3870 else
3871 {
3872 if (diff[Order[DIM-1]]>0.0) order = -Sign[DIM-1];
3873 else order = Sign[DIM-1];
3874 }
3875
3876 if (order==1) SETMUP(theMatrix,1);
3877 else SETMDOWN(theMatrix,1);
3878 }
3879 }
3880
3881 return (0);
3882 }
3883
3884 /****************************************************************************/
3885 /** \brief Dependency function for lexicographic ordering
3886
3887 * @param theGrid - pointer to grid
3888 * @param data - option string from orderv command
3889
3890 This function defines a dependency function for lexicographic ordering.
3891
3892 * @return <ul>
3893 * <li> 0 if ok
3894 * <li> 1 if error occured.
3895 </ul>
3896 */
3897 /****************************************************************************/
3898
StrongLexAlgDep(GRID * theGrid,const char * data)3899 static INT StrongLexAlgDep (GRID *theGrid, const char *data)
3900 {
3901 MULTIGRID *theMG;
3902 VECTOR *theVector,*NBVector;
3903 MATRIX *theMatrix;
3904 DOUBLE_VECTOR pos,nbpos;
3905 DOUBLE diff[DIM];
3906 INT i,order,res;
3907 INT Sign[DIM],Order[DIM],xused,yused,zused,error;
3908 char ord[3];
3909
3910 /* read ordering directions */
3911 #ifdef __TWODIM__
3912 res = sscanf(data,expandfmt("%2[rlud]"),ord);
3913 #else
3914 res = sscanf(data,expandfmt("%3[rlbfud]"),ord);
3915 #endif
3916 if (res!=1)
3917 {
3918 PrintErrorMessage('E',"LexAlgDep","could not read order type");
3919 return(1);
3920 }
3921 if (strlen(ord)!=DIM)
3922 {
3923 #ifdef __TWODIM__
3924 PrintErrorMessage('E',"LexAlgDep","specify 2 chars out of 'rlud'");
3925 #else
3926 PrintErrorMessage('E',"LexAlgDep","specify 3 chars out of 'rlbfud'");
3927 #endif
3928 return(1);
3929 }
3930 error = xused = yused = zused = false;
3931 for (i=0; i<DIM; i++)
3932 switch (ord[i])
3933 {
3934 case 'r' :
3935 if (xused) error = true;
3936 xused = true;
3937 Order[i] = _X_; Sign[i] = 1; break;
3938 case 'l' :
3939 if (xused) error = true;
3940 xused = true;
3941 Order[i] = _X_; Sign[i] = -1; break;
3942
3943 #ifdef __TWODIM__
3944 case 'u' :
3945 if (yused) error = true;
3946 yused = true;
3947 Order[i] = _Y_; Sign[i] = 1; break;
3948 case 'd' :
3949 if (yused) error = true;
3950 yused = true;
3951 Order[i] = _Y_; Sign[i] = -1; break;
3952 #else
3953 case 'b' :
3954 if (yused) error = true;
3955 yused = true;
3956 Order[i] = _Y_; Sign[i] = 1; break;
3957 case 'f' :
3958 if (yused) error = true;
3959 yused = true;
3960 Order[i] = _Y_; Sign[i] = -1; break;
3961
3962 case 'u' :
3963 if (zused) error = true;
3964 zused = true;
3965 Order[i] = _Z_; Sign[i] = 1; break;
3966 case 'd' :
3967 if (zused) error = true;
3968 zused = true;
3969 Order[i] = _Z_; Sign[i] = -1; break;
3970 #endif
3971 }
3972 if (error)
3973 {
3974 PrintErrorMessage('E',"LexAlgDep","bad combination of 'rludr' or 'rlbfud' resp.");
3975 return(1);
3976 }
3977
3978 theMG = MYMG(theGrid);
3979
3980 /* find an approximate measure for the mesh size */
3981 InvMeshSize = POW2(GLEVEL(theGrid)) * pow(NN(GRID_ON_LEVEL(theMG,0)),1.0/DIM);
3982
3983 for (theVector=FIRSTVECTOR(theGrid); theVector!=NULL; theVector=SUCCVC(theVector))
3984 {
3985 VectorPosition(theVector,pos);
3986
3987 for (theMatrix=MNEXT(VSTART(theVector)); theMatrix!=NULL; theMatrix=MNEXT(theMatrix))
3988 {
3989 NBVector = MDEST(theMatrix);
3990
3991 SETMUP(theMatrix,0);
3992 SETMDOWN(theMatrix,0);
3993 SETMUSED(theMatrix,0);
3994
3995 VectorPosition(NBVector,nbpos);
3996
3997 V_DIM_SUBTRACT(nbpos,pos,diff);
3998 V_DIM_SCALE(InvMeshSize,diff);
3999
4000 if (fabs(diff[Order[DIM-1]])<ORDERRES)
4001 {
4002 #ifdef __THREEDIM__
4003 if (fabs(diff[Order[DIM-2]])<ORDERRES)
4004 {
4005 if (diff[Order[DIM-3]]>0.0) order = -Sign[DIM-3];
4006 else order = Sign[DIM-3];
4007 }
4008 else
4009 #endif
4010 if (diff[Order[DIM-2]]>0.0) order = -Sign[DIM-2];
4011 else order = Sign[DIM-2];
4012 SETMUSED(theMatrix,1);
4013 }
4014 else
4015 {
4016 if (diff[Order[DIM-1]]>0.0) order = -Sign[DIM-1];
4017 else order = Sign[DIM-1];
4018 }
4019 switch (order)
4020 {
4021 case 1 : SETMUP(theMatrix,1); break;
4022 case -1 : SETMDOWN(theMatrix,1); break;
4023 case 0 : SETMUP(theMatrix,1); SETMDOWN(theMatrix,1); break;
4024 }
4025 }
4026 }
4027 for (theVector=FIRSTVECTOR(theGrid); theVector!=NULL; theVector=SUCCVC(theVector))
4028 {
4029 SETVCUSED(theVector,0);
4030 SETVCFLAG(theVector,0);
4031
4032 for (theMatrix=MNEXT(VSTART(theVector)); theMatrix!=NULL; theMatrix=MNEXT(theMatrix))
4033 if (MUP(theMatrix) && !MUSED(theMatrix))
4034 break;
4035 if (theMatrix==NULL)
4036 SETVCUSED(theVector,1);
4037 }
4038 for (theVector=FIRSTVECTOR(theGrid); theVector!=NULL; theVector=SUCCVC(theVector))
4039 for (theMatrix=MNEXT(VSTART(theVector)); theMatrix!=NULL; theMatrix=MNEXT(theMatrix))
4040 if (MUSED(theMatrix) && MUSED(MADJ(theMatrix)))
4041 {
4042 SETMUP(theMatrix,1);
4043 SETMDOWN(theMatrix,1);
4044 }
4045
4046 return (0);
4047 }
4048
4049 /****************************************************************************/
4050 /** \brief Init algebra
4051 *
4052 * This function inits algebra.
4053 *
4054 * @return <ul>
4055 * <li> 0 if ok </li>
4056 * <li> 1 if error occured. </li>
4057 * </ul>
4058 */
4059 /****************************************************************************/
4060
InitAlgebra(void)4061 INT NS_DIM_PREFIX InitAlgebra (void)
4062 {
4063 INT i;
4064
4065 /* install the /Alg Dep directory */
4066 if (ChangeEnvDir("/")==NULL)
4067 {
4068 PrintErrorMessage('F',"InitAlgebra","could not changedir to root");
4069 return(__LINE__);
4070 }
4071 theAlgDepDirID = GetNewEnvDirID();
4072 if (MakeEnvItem("Alg Dep",theAlgDepDirID,sizeof(ENVDIR))==NULL)
4073 {
4074 PrintErrorMessage('F',"InitAlgebra","could not install '/Alg Dep' dir");
4075 return(__LINE__);
4076 }
4077 theAlgDepVarID = GetNewEnvVarID();
4078
4079 /* install the /FindCut directory */
4080 if (ChangeEnvDir("/")==NULL)
4081 {
4082 PrintErrorMessage('F',"InitAlgebra","could not changedir to root");
4083 return(__LINE__);
4084 }
4085 theFindCutDirID = GetNewEnvDirID();
4086 if (MakeEnvItem("FindCut",theFindCutDirID,sizeof(ENVDIR))==NULL)
4087 {
4088 PrintErrorMessage('F',"InitAlgebra","could not install '/FindCut' dir");
4089 return(__LINE__);
4090 }
4091 theFindCutVarID = GetNewEnvVarID();
4092
4093 /* init standard algebraic dependencies */
4094 if (CreateAlgebraicDependency ("lex",LexAlgDep)==NULL) return(__LINE__);
4095 if (CreateAlgebraicDependency ("stronglex",StrongLexAlgDep)==NULL) return(__LINE__);
4096
4097 /* init default find cut proc */
4098 if (CreateFindCutProc ("lex",FeedbackVertexVectors)==NULL) return(__LINE__);
4099
4100 for (i=0; i<MAXVOBJECTS; i++)
4101 switch (i)
4102 {
4103 case NODEVEC : ObjTypeName[i] = "nd"; break;
4104 case EDGEVEC : ObjTypeName[i] = "ed"; break;
4105 case ELEMVEC : ObjTypeName[i] = "el"; break;
4106 case SIDEVEC : ObjTypeName[i] = "si"; break;
4107 default : ObjTypeName[i] = "";
4108 }
4109
4110 return (0);
4111 }
4112
4113
4114 /*@}*/
4115