1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 /*! \file std_domain.c
4 * \ingroup std
5 */
6
7 /** \addtogroup std
8 *
9 * @{
10 */
11
12 /****************************************************************************/
13 /* */
14 /* File: std_domain.c */
15 /* */
16 /* Purpose: standard ug domain description */
17 /* */
18 /* Author: Klaus Johannsen / Christian Wieners */
19 /* Institut fuer Computeranwendungen III */
20 /* Universitaet Stuttgart */
21 /* Pfaffenwaldring 27 */
22 /* 70550 Stuttgart */
23 /* email: ug@ica3.uni-stuttgart.de */
24 /* */
25 /* History: Feb 18 1996 begin, ug version 3.1 */
26 /* Sep 12 1996 ug version 3.4 */
27 /* Apr 9 1998 first step to Marc */
28 /* */
29 /* Remarks: */
30 /* */
31 /****************************************************************************/
32
33 /****************************************************************************/
34 /* */
35 /* include files */
36 /* system include files */
37 /* application include files */
38 /* */
39 /****************************************************************************/
40
41 #include <config.h>
42
43 /* standard C library */
44 #include <cstdlib>
45 #include <cstddef>
46 #include <cstdio>
47 #include <cstring>
48 #include <cassert>
49 #include <cmath>
50
51 /* standard C++ library */
52 /* set needed in BVP_Init */
53 #include <set>
54
55 /* low modules */
56 #include <dune/uggrid/low/architecture.h>
57 #include <dune/uggrid/low/bio.h>
58 #include <dune/uggrid/low/debug.h>
59 #include <dune/uggrid/low/fileopen.h>
60 #include <dune/uggrid/low/heaps.h>
61 #include <dune/uggrid/low/misc.h>
62 #include <dune/uggrid/low/namespace.h>
63 #include <dune/uggrid/low/scan.h>
64 #include <dune/uggrid/low/ugenv.h>
65 #include <dune/uggrid/low/ugtypes.h>
66
67 #include <dune/uggrid/gm/evm.h>
68
69 /* dev modules */
70 #include <dune/uggrid/ugdevices.h>
71
72 /* domain module */
73 #include "domain.h"
74 #include "std_internal.h"
75
76 USING_UGDIM_NAMESPACE
77 USING_UG_NAMESPACE
78
79 #ifdef ModelP
80 using namespace PPIF;
81 #endif
82
83 /****************************************************************************/
84 /* */
85 /* defines in the following order */
86 /* */
87 /* compile time constants defining static data size (i.e. arrays) */
88 /* other constants */
89 /* macros */
90 /* */
91 /****************************************************************************/
92
93 #define SMALL_DIFF SMALL_C*100
94 #define RESOLUTION 100
95
96 #define DEFAULTDOMMEMORY 50000
97
98 #define OPTIONLEN 32
99
100 #define V2_LINCOMB(a,A,b,B,C) {(C)[0] = (a)*(A)[0] + (b)*(B)[0];\
101 (C)[1] = (a)*(A)[1] + (b)*(B)[1];}
102
103 #define V2_EUKLIDNORM_OF_DIFF(A,B,b) (b) = sqrt((double)(((A)[0]-(B)[0])*((A)[0]-(B)[0])+((A)[1]-(B)[1])*((A)[1]-(B)[1])));
104
105 #define V3_EUKLIDNORM_OF_DIFF(A,B,b) (b) = (sqrt((double)(((A)[0]-(B)[0])*((A)[0]-(B)[0])+((A)[1]-(B)[1])*((A)[1]-(B)[1])+((A)[2]-(B)[2])*((A)[2]-(B)[2]))));
106
107 #ifdef __TWODIM__
108 #define V_DIM_EUKLIDNORM_OF_DIFF(A,B,b) V2_EUKLIDNORM_OF_DIFF(A,B,b)
109 #else
110 #define V_DIM_EUKLIDNORM_OF_DIFF(A,B,b) V3_EUKLIDNORM_OF_DIFF(A,B,b)
111 #endif
112
113 /****************************************************************************/
114 /* */
115 /* data structures used in this source file (exported data structures are */
116 /* in the corresponding include file!) */
117 /* */
118 /****************************************************************************/
119
120 /****************************************************************************/
121 /* */
122 /* definition of exported global variables */
123 /* */
124 /****************************************************************************/
125
126 /****************************************************************************/
127 /* */
128 /* definition of variables global to this source file only (static!) */
129 /* */
130 /****************************************************************************/
131
132 static INT theProblemDirID; /*!< env type for Problem dir */
133 static INT theBdryCondVarID; /*!< env type for Problem vars */
134
135 static INT theDomainDirID; /*!< env type for Domain dir */
136 static INT theBdrySegVarID; /*!< env type for bdry segment vars */
137 static INT theLinSegVarID; /*!< env type for linear segment vars */
138
139 static INT theBVPDirID; /*!< env type for BVP dir */
140
141 static STD_BVP *currBVP;
142
143 /****************************************************************************/
144 /* */
145 /* forward declarations of functions used before they are defined */
146 /* */
147 /****************************************************************************/
148
149 static INT BndPointGlobal (const BNDP * aBndP, DOUBLE * global);
150 static INT PatchGlobal (const PATCH * p, DOUBLE * lambda, DOUBLE * global);
151
152 /****************************************************************************/
153 /** \brief Allocate a new BNDCOND structure
154 *
155 * @param theBCond - the boundary condition
156 *
157 * This function gets next BNDCOND structure
158 *
159 * @return <ul>
160 * <li> pointer to BOUNDARY_CONDITION </li>
161 * <li> NULL if not found or error. </li>
162 * </ul>
163 */
164 /****************************************************************************/
165 static BOUNDARY_CONDITION *
GetNextBoundaryCondition(BOUNDARY_CONDITION * theBCond)166 GetNextBoundaryCondition (BOUNDARY_CONDITION * theBCond)
167 {
168 ENVITEM *theItem;
169
170 theItem = (ENVITEM *) theBCond;
171
172 do
173 theItem = NEXT_ENVITEM (theItem);
174 while ((theItem != NULL) && (ENVITEM_TYPE (theItem) != theBdryCondVarID));
175
176 return ((BOUNDARY_CONDITION *) theItem);
177 }
178
179 /****************************************************************************/
180 /** \brief Get first BNDCOND structure of `theProblem`
181 *
182 * @param theProblem - pointer to PROBLEM
183 *
184 * This function gets the first BNDCOND structure of a problem
185 *
186 * @return <ul>
187 * <li> pointer to BOUNDARY_CONDITION </li>
188 * <li> NULL if not found or error. </li>
189 * </ul>
190 */
191 /****************************************************************************/
192 static BOUNDARY_CONDITION *
GetFirstBoundaryCondition(PROBLEM * theProblem)193 GetFirstBoundaryCondition (PROBLEM * theProblem)
194 {
195 ENVITEM *theItem;
196
197 theItem = ENVITEM_DOWN (theProblem);
198
199 if (ENVITEM_TYPE (theItem) == theBdryCondVarID)
200 return ((BOUNDARY_CONDITION *) theItem);
201 else
202 return (GetNextBoundaryCondition ((BOUNDARY_CONDITION *) theItem));
203 }
204
205 /****************************************************************************/
206 /** \brief Create a new DOMAIN data structure with part description
207 *
208 * @param name - name of the domain
209 * @param MidPoint - coordinates of some inner point
210 * @param radius - radius of a circle, containing the domain
211 * @param segments - number of the boundary segments
212 * @param corners - number of corners
213 * @param Convex - 0, if convex, 1 - else
214 * @param nParts - number of parts in the domain
215 * @param dpi - description of the parts for lines, segments, points
216 *
217 * This function allocates and initializes a new DOMAIN data structure in the
218 * /domains directory in the environment.
219 * Additinally domain parts will defined.
220 *
221 * @return <ul>
222 * <li> pointer to a DOMAIN </li>
223 * <li> NULL if out of memory. </li>
224 * </ul>
225 */
226 /****************************************************************************/
227
228 void * NS_DIM_PREFIX
CreateDomainWithParts(const char * name,INT segments,INT corners,INT nParts,const DOMAIN_PART_INFO * dpi)229 CreateDomainWithParts (const char *name,
230 INT segments, INT corners, INT nParts,
231 const DOMAIN_PART_INFO * dpi)
232 {
233 DOMAIN *newDomain;
234
235 /* change to /domains directory */
236 if (ChangeEnvDir ("/Domains") == NULL)
237 return (NULL);
238
239 /* allocate new domain structure */
240 newDomain = (DOMAIN *) MakeEnvItem (name, theDomainDirID, sizeof (DOMAIN));
241 if (newDomain == NULL)
242 return (NULL);
243
244 /* fill in data */
245 DOMAIN_NSEGMENT (newDomain) = segments;
246 DOMAIN_NCORNER (newDomain) = corners;
247 DOMAIN_NPARTS (newDomain) = nParts;
248 DOMAIN_PARTINFO (newDomain) = dpi;
249
250 if (ChangeEnvDir (name) == NULL)
251 return (NULL);
252 UserWrite ("domain ");
253 UserWrite (name);
254 UserWrite (" installed\n");
255
256 return (newDomain);
257 }
258
259 /****************************************************************************/
260 /** \brief Create a new DOMAIN data structure
261 *
262 * @param name - name of the domain
263 * @param MidPoint - coordinates of some inner point
264 * @param radius - radius of a circle, containing the domain
265 * @param segments - number of the boundary segments
266 * @param corners - number of corners
267 * @param Convex - 0, if convex, 1 - else
268 *
269 * This function allocates and initializes a new DOMAIN data structure in the
270 * /domains directory in the environment.
271 *
272 * @return <ul>
273 * <li> pointer to a DOMAIN </li>
274 * <li> NULL if out of memory. </li>
275 * </ul>
276 */
277 /****************************************************************************/
278
279 void *NS_DIM_PREFIX
CreateDomain(const char * name,INT segments,INT corners)280 CreateDomain (const char *name, INT segments,
281 INT corners)
282 {
283 return (CreateDomainWithParts
284 (name, segments, corners, 1, NULL));
285 }
286
287 /****************************************************************************/
288 /** \brief Get a pointer to a domain structure
289 *
290 * @param name - name of the domain
291 *
292 * This function searches the environment for a domain with the name `name`
293 * and return a pointer to the domain structure.
294 *
295 * @return <ul>
296 * <li> pointer to a DOMAIN </li>
297 * <li> NULL if not found or error. </li>
298 * </ul>
299 */
300 /****************************************************************************/
301
302 DOMAIN *NS_DIM_PREFIX
GetDomain(const char * name)303 GetDomain (const char *name)
304 {
305 return ((DOMAIN *)
306 SearchEnv (name, "/Domains", theDomainDirID, theDomainDirID));
307 }
308
309 /****************************************************************************/
310 /** \brief Remove a domain structure
311 *
312 * @param name - name of the domain
313 *
314 * This function searches the environment for a domain with the name `name`
315 * and removes it.
316 *
317 */
318 /****************************************************************************/
319
320 void NS_DIM_PREFIX
RemoveDomain(const char * name)321 RemoveDomain (const char *name)
322 {
323 ENVITEM* d = SearchEnv (name, "/Domains", theDomainDirID, theDomainDirID);
324 if (d!=0) {
325 d->v.locked = false;
326 RemoveEnvDir(d);
327 }
328 }
329
330 /****************************************************************************/
331 /** \brief Create a new BOUNDARY_SEGMENT
332 *
333 * @param name - name of the boundary segment
334 * @param left - id of left subdomain
335 * @param right - id of right subdomain
336 * @param id - id of this boundary segment
337 * @param type - type of the boundary segment
338 * @param point - the endpoints of the boundary segment
339 * @param alpha - list where the parameter interval begins
340 * @param beta - list where the parameter interval ends
341 * @param BndSegFunc - function mapping parameters
342 * @param data - user defined space
343 *
344 * This function allocates and initializes a new BOUNDARY_SEGMENT for the previously
345 * allocated DOMAIN structure.
346 *
347 * @return <ul>
348 * <li> Pointer to a BOUNDARY_SEGMENT </li>
349 * <li> NULL if out of memory. </li>
350 * </ul>
351 */
352 /****************************************************************************/
353
354 void *NS_DIM_PREFIX
CreateBoundarySegment(const char * name,INT left,INT right,INT id,enum BoundaryType type,const INT * point,const DOUBLE * alpha,const DOUBLE * beta,BndSegFuncPtr BndSegFunc,void * data)355 CreateBoundarySegment (const char *name,
356 INT left, INT right, INT id, enum BoundaryType type,
357 const INT * point, const DOUBLE * alpha, const DOUBLE * beta,
358 BndSegFuncPtr BndSegFunc, void *data)
359 {
360 BOUNDARY_SEGMENT *newSegment;
361 INT i;
362
363 /* allocate the boundary segment */
364 newSegment =
365 (BOUNDARY_SEGMENT *) MakeEnvItem (name, theBdrySegVarID,
366 sizeof (BOUNDARY_SEGMENT));
367 if (newSegment == NULL)
368 return (NULL);
369
370 /* fill in data */
371 newSegment->left = left;
372 newSegment->right = right;
373 newSegment->id = id;
374 newSegment->segType = type;
375 for (i = 0; i < CORNERS_OF_BND_SEG; i++)
376 newSegment->points[i] = point[i];
377 for (i = 0; i < DIM_OF_BND; i++)
378 {
379 newSegment->alpha[i] = alpha[i];
380 newSegment->beta[i] = beta[i];
381 }
382 newSegment->BndSegFunc = BndSegFunc;
383 newSegment->data = data;
384
385 return (newSegment);
386 }
387
388 /****************************************************************************/
389 /** \brief Get next boundary segment of a domain
390 *
391 * @param theBSeg - pointer to a boundary segment
392 *
393 * This function gets the next boundary segment of a domain.
394 *
395 * @return <ul>
396 * <li> pointer to next BOUNDARY_SEGMENT </li>
397 * <li> NULL if no more found. </li>
398 * </ul>
399 */
400 /****************************************************************************/
401
402 static BOUNDARY_SEGMENT *
GetNextBoundarySegment(BOUNDARY_SEGMENT * theBSeg)403 GetNextBoundarySegment (BOUNDARY_SEGMENT * theBSeg)
404 {
405 ENVITEM *theItem;
406
407 theItem = (ENVITEM *) theBSeg;
408
409 do
410 theItem = NEXT_ENVITEM (theItem);
411 while ((theItem != NULL) && (ENVITEM_TYPE (theItem) != theBdrySegVarID));
412
413 return ((BOUNDARY_SEGMENT *) theItem);
414 }
415
416 /****************************************************************************/
417 /** \brief Get first boundary segment of a domain
418 *
419 * @param theDomain - pointer to the domain structure
420 *
421 * This function returns the first boundary segment of a domain.
422 *
423 * @return <ul>
424 * <li> pointer to a BOUNDARY_SEGMENT </li>
425 * <li> NULL if not found or error. </li>
426 * </ul>
427 */
428 /****************************************************************************/
429
430 static BOUNDARY_SEGMENT *
GetFirstBoundarySegment(DOMAIN * theDomain)431 GetFirstBoundarySegment (DOMAIN * theDomain)
432 {
433 ENVITEM *theItem;
434
435 theItem = ENVITEM_DOWN (theDomain);
436
437 if (ENVITEM_TYPE (theItem) == theBdrySegVarID)
438 return ((BOUNDARY_SEGMENT *) theItem);
439 else
440 return (GetNextBoundarySegment ((BOUNDARY_SEGMENT *) theItem));
441 }
442
443 /****************************************************************************/
444 /** \brief Create a new LINEAR_SEGMENT
445 *
446 * @param name - name of the boundary segment
447 * @param left - id of left subdomain
448 * @param right - id of right subdomain
449 * @param id - id of this boundary segment
450 * @param n - number of corners
451 * @param point - the endpoints of the boundary segment
452 * @param x - coordinates
453 *
454 * This function allocates and initializes a new LINEAR_SEGMENT
455 * for the previously allocated DOMAIN structure.
456 *
457 * @return <ul>
458 * <li> Pointer to a LINEAR_SEGMENT </li>
459 * <li> NULL if out of memory. </li>
460 * </ul>
461 */
462 /****************************************************************************/
463
464 void *NS_DIM_PREFIX
CreateLinearSegment(const char * name,INT left,INT right,INT id,INT n,const INT * point,DOUBLE x[CORNERS_OF_BND_SEG][DIM])465 CreateLinearSegment (const char *name,
466 INT left, INT right, INT id,
467 INT n, const INT * point,
468 DOUBLE x[CORNERS_OF_BND_SEG][DIM])
469 {
470 LINEAR_SEGMENT *newSegment;
471 INT i, k;
472
473 if (n > CORNERS_OF_BND_SEG)
474 return (NULL);
475
476 /* allocate the boundary segment */
477 newSegment = (LINEAR_SEGMENT *)
478 MakeEnvItem (name, theLinSegVarID, sizeof (LINEAR_SEGMENT));
479
480 if (newSegment == NULL)
481 return (NULL);
482
483 /* fill in data */
484 newSegment->left = left;
485 newSegment->right = right;
486 newSegment->id = id;
487 newSegment->n = n;
488 for (i = 0; i < n; i++)
489 {
490 newSegment->points[i] = point[i];
491 for (k = 0; k < DIM; k++)
492 newSegment->x[i][k] = x[i][k];
493 }
494
495 return (newSegment);
496 }
497
GetBoundarySegmentId(BNDS * boundarySegment)498 UINT NS_DIM_PREFIX GetBoundarySegmentId(BNDS* boundarySegment)
499 {
500 BND_PS *ps;
501 PATCH *patch;
502
503 ps = (BND_PS *) boundarySegment;
504 patch = currBVP->patches[ps->patch_id];
505 if (patch == NULL) {
506 PrintErrorMessageF ('E', "GetBoundarySegmentId", "invalid argument");
507 return 0;
508 }
509
510 /* The ids in the patch data structure are consecutive but they
511 start at currBVP->sideoffset instead of zero. */
512 return PATCH_ID(patch) - currBVP->sideoffset;
513 }
514
515
516 static LINEAR_SEGMENT *
GetNextLinearSegment(LINEAR_SEGMENT * theBSeg)517 GetNextLinearSegment (LINEAR_SEGMENT * theBSeg)
518 {
519 ENVITEM *theItem;
520
521 theItem = (ENVITEM *) theBSeg;
522
523 do
524 theItem = NEXT_ENVITEM (theItem);
525 while ((theItem != NULL) && (ENVITEM_TYPE (theItem) != theLinSegVarID));
526
527 return ((LINEAR_SEGMENT *) theItem);
528 }
529
530 static LINEAR_SEGMENT *
GetFirstLinearSegment(DOMAIN * theDomain)531 GetFirstLinearSegment (DOMAIN * theDomain)
532 {
533 ENVITEM *theItem;
534
535 theItem = ENVITEM_DOWN (theDomain);
536
537 if (ENVITEM_TYPE (theItem) == theLinSegVarID)
538 return ((LINEAR_SEGMENT *) theItem);
539 else
540 return (GetNextLinearSegment ((LINEAR_SEGMENT *) theItem));
541 }
542
543 /* configuring a domain */
STD_BVP_Configure(INT argc,char ** argv)544 static INT STD_BVP_Configure(INT argc, char **argv)
545 {
546 STD_BVP *theBVP;
547 DOMAIN *theDomain;
548 char BVPName[NAMESIZE];
549 char DomainName[NAMESIZE];
550 INT i;
551
552 /* get BVP name */
553 if (sscanf(argv[0], expandfmt(" configure %" NAMELENSTR "[ -~]"), BVPName) != 1
554 || strlen(BVPName) == 0)
555 return 1;
556
557 theBVP = (STD_BVP *) BVP_GetByName(BVPName);
558 if (theBVP == nullptr)
559 return 1;
560
561 for (i=0; i<argc; i++)
562 if (argv[i][0] == 'd' && argv[i][1] == ' ')
563 if (sscanf(argv[i], expandfmt("d %" NAMELENSTR "[ -~]"), DomainName) != 1
564 || strlen(DomainName) == 0)
565 continue;
566
567 theDomain = GetDomain(DomainName);
568 if (theDomain == nullptr)
569 return 1;
570
571 theBVP->Domain = theDomain;
572
573 return 0;
574 }
575
576 BVP *NS_DIM_PREFIX
CreateBoundaryValueProblem(const char * BVPName,BndCondProcPtr theBndCond,int numOfCoeffFct,CoeffProcPtr coeffs[],int numOfUserFct,UserProcPtr userfct[])577 CreateBoundaryValueProblem (const char *BVPName, BndCondProcPtr theBndCond,
578 int numOfCoeffFct, CoeffProcPtr coeffs[],
579 int numOfUserFct, UserProcPtr userfct[])
580 {
581 STD_BVP *theBVP;
582 INT i, n;
583
584 /* change to /BVP directory */
585 if (ChangeEnvDir ("/BVP") == NULL)
586 return (NULL);
587
588 /* allocate new domain structure */
589 n = (numOfCoeffFct + numOfUserFct - 1) * sizeof (void *);
590 theBVP =
591 (STD_BVP *) MakeEnvItem (BVPName, theBVPDirID, sizeof (STD_BVP) + n);
592 if (theBVP == NULL)
593 return (NULL);
594 if (ChangeEnvDir (BVPName) == NULL)
595 return (NULL);
596
597 theBVP->numOfCoeffFct = numOfCoeffFct;
598 theBVP->numOfUserFct = numOfUserFct;
599 for (i = 0; i < numOfCoeffFct; i++)
600 theBVP->CU_ProcPtr[i] = (void *) (coeffs[i]);
601 for (i = 0; i < numOfUserFct; i++)
602 theBVP->CU_ProcPtr[i + numOfCoeffFct] = (void *) (userfct[i]);
603
604 STD_BVP_S2P_PTR (theBVP) = NULL;
605
606 theBVP->Domain = NULL;
607 theBVP->Problem = NULL;
608 theBVP->ConfigProc = STD_BVP_Configure;
609 theBVP->GeneralBndCond = theBndCond;
610
611 UserWriteF ("BVP %s installed.\n", BVPName);
612
613 return ((BVP *) theBVP);
614 }
615
616 static INT
GetNumberOfPatches(PATCH * p)617 GetNumberOfPatches (PATCH * p)
618 {
619 switch (PATCH_TYPE (p))
620 {
621 case PARAMETRIC_PATCH_TYPE :
622 case LINEAR_PATCH_TYPE :
623 return (1);
624 case POINT_PATCH_TYPE :
625 return (POINT_PATCH_N (p));
626 #ifdef __THREEDIM__
627 case LINE_PATCH_TYPE :
628 return (LINE_PATCH_N (p));
629 #endif
630 }
631
632 return (-1);
633 }
634
635 static INT
GetPatchId(PATCH * p,INT i)636 GetPatchId (PATCH * p, INT i)
637 {
638 switch (PATCH_TYPE (p))
639 {
640 case LINEAR_PATCH_TYPE :
641 case PARAMETRIC_PATCH_TYPE :
642 return (PATCH_ID (p));
643 case POINT_PATCH_TYPE :
644 return (POINT_PATCH_PID (p, i));
645 #ifdef __THREEDIM__
646 case LINE_PATCH_TYPE :
647 return (LINE_PATCH_PID (p, i));
648 #endif
649 }
650
651 assert(0);
652 return (-1);
653 }
654
655 static BNDP *
CreateBndPOnPoint(HEAP * Heap,PATCH * p)656 CreateBndPOnPoint (HEAP * Heap, PATCH * p)
657 {
658 BND_PS *ps;
659 PATCH *pp;
660 INT j, l, m;
661
662 if (PATCH_TYPE (p) != POINT_PATCH_TYPE)
663 return (NULL);
664
665 PRINTDEBUG (dom, 1, (" p %d\n", PATCH_ID (p)));
666 for (l = 0; l < GetNumberOfPatches (p); l++)
667 PRINTDEBUG (dom, 1, (" bp pid %d\n", GetPatchId (p, l)));
668
669 m = POINT_PATCH_N (p);
670 ps = (BND_PS *) GetFreelistMemory (Heap, sizeof (BND_PS)
671 + (m - 1) * sizeof (COORD_BND_VECTOR));
672 if (ps == NULL)
673 REP_ERR_RETURN (NULL);
674 ps->n = m;
675 ps->patch_id = PATCH_ID (p);
676
677 for (j = 0; j < m; j++)
678 {
679 pp = currBVP->patches[POINT_PATCH_PID (p, j)];
680 if (PATCH_TYPE (pp) == PARAMETRIC_PATCH_TYPE)
681 {
682 PRINTDEBUG (dom, 1, ("cp r %f %f %f %f\n",
683 PARAM_PATCH_RANGE (pp)[0][0],
684 PARAM_PATCH_RANGE (pp)[0][1],
685 PARAM_PATCH_RANGE (pp)[1][0],
686 PARAM_PATCH_RANGE (pp)[1][1]));
687 switch (POINT_PATCH_CID (p, j))
688 {
689 #ifdef __TWODIM__
690 case 0 :
691 ps->local[j][0] = PARAM_PATCH_RANGE (pp)[0][0];
692 break;
693 case 1 :
694 ps->local[j][0] = PARAM_PATCH_RANGE (pp)[1][0];
695 break;
696 #endif
697 #ifdef __THREEDIM__
698 case 0 :
699 ps->local[j][0] = PARAM_PATCH_RANGE (pp)[0][0];
700 ps->local[j][1] = PARAM_PATCH_RANGE (pp)[0][1];
701 break;
702 case 1 :
703 ps->local[j][0] = PARAM_PATCH_RANGE (pp)[1][0];
704 ps->local[j][1] = PARAM_PATCH_RANGE (pp)[0][1];
705 break;
706 case 2 :
707 ps->local[j][0] = PARAM_PATCH_RANGE (pp)[1][0];
708 ps->local[j][1] = PARAM_PATCH_RANGE (pp)[1][1];
709 break;
710 case 3 :
711 ps->local[j][0] = PARAM_PATCH_RANGE (pp)[0][0];
712 ps->local[j][1] = PARAM_PATCH_RANGE (pp)[1][1];
713 break;
714 #endif
715 }
716 PRINTDEBUG (dom, 1, ("mesh loc j %d pid %d cid %d loc %f %f\n",
717 j,
718 POINT_PATCH_PID (p, j),
719 POINT_PATCH_CID (p, j),
720 ps->local[j][0], ps->local[j][1]));
721 }
722 else if (PATCH_TYPE (pp) == LINEAR_PATCH_TYPE)
723 {
724 switch (POINT_PATCH_CID (p, j))
725 {
726 #ifdef __TWODIM__
727 case 0 :
728 ps->local[j][0] = 0.0;
729 break;
730 case 1 :
731 ps->local[j][0] = 1.0;
732 break;
733 #endif
734 #ifdef __THREEDIM__
735 case 0 :
736 ps->local[j][0] = 0.0;
737 ps->local[j][1] = 0.0;
738 break;
739 case 1 :
740 ps->local[j][0] = 1.0;
741 ps->local[j][1] = 0.0;
742 break;
743 case 2 :
744 /* This coordinate depends on whether we're looking at a triangle or a quadrilateral */
745 ps->local[j][0] = (LINEAR_PATCH_N(pp)==3) ? 0.0 : 1.0;
746 ps->local[j][1] = 1.0;
747 break;
748 case 3 :
749 ps->local[j][0] = 0.0;
750 ps->local[j][1] = 1.0;
751 break;
752 #endif
753 }
754 }
755 }
756 if (!PATCH_IS_FIXED (p))
757 {
758 /* store global coordinates */
759 BND_DATA (ps) = GetFreelistMemory (Heap, DIM * sizeof (DOUBLE));
760 if (BND_DATA (ps) == NULL)
761 REP_ERR_RETURN (NULL);
762
763 if (BndPointGlobal ((BNDP *) ps, (DOUBLE *) BND_DATA (ps)))
764 REP_ERR_RETURN (NULL);
765 }
766 return ((BNDP *) ps);
767 }
768
769 static INT
CreateCornerPoints(HEAP * Heap,STD_BVP * theBVP,BNDP ** bndp)770 CreateCornerPoints (HEAP * Heap, STD_BVP * theBVP, BNDP ** bndp)
771 {
772 int i;
773
774 for (i = 0; i < theBVP->ncorners; i++)
775 {
776 bndp[i] = CreateBndPOnPoint (Heap, theBVP->patches[i]);
777 if (bndp[i] == NULL)
778 REP_ERR_RETURN (1);
779 }
780
781 for (i = 0; i < theBVP->ncorners; i++)
782 PRINTDEBUG (dom, 1, (" id %d\n", PATCH_ID (theBVP->patches[i])));
783
784 return (0);
785 }
786
787 #ifdef __THREEDIM__
788 /** \brief Add a line between two vertices to the boundary data structure
789
790 \param i, j The line goes from vertex i to vertex j
791 \param corners Array with pointers to all vertices
792 */
793 static void
CreateLine(INT i,INT j,HEAP * Heap,PATCH * thePatch,PATCH ** corners,PATCH ** lines,PATCH ** sides,INT * nlines,INT * err)794 CreateLine(INT i, INT j, HEAP *Heap, PATCH *thePatch, PATCH **corners, PATCH **lines, PATCH **sides,
795 INT *nlines, INT *err)
796 {
797 INT k, n, m;
798 INT freePatches;
799
800 k = 0;
801 for (n = 0; n < POINT_PATCH_N (corners[i]); n++)
802 for (m = 0; m < POINT_PATCH_N (corners[j]); m++)
803 if (POINT_PATCH_PID (corners[i], n) ==
804 POINT_PATCH_PID (corners[j], m))
805 k++;
806 /* points share one patch only and lie on opposite corners of this patch */
807 if (k < 2)
808 return;
809
810 thePatch =
811 (PATCH *) GetFreelistMemory (Heap, sizeof (LINE_PATCH)
812 + (k -
813 1) * sizeof (struct line_on_patch));
814 if (thePatch == NULL)
815 return;
816 PATCH_TYPE (thePatch) = LINE_PATCH_TYPE;
817 PATCH_ID (thePatch) = *nlines;
818 LINE_PATCH_C0 (thePatch) = i;
819 LINE_PATCH_C1 (thePatch) = j;
820 k = 0;
821 freePatches = 0;
822 for (n = 0; n < POINT_PATCH_N (corners[i]); n++)
823 for (m = 0; m < POINT_PATCH_N (corners[j]); m++)
824 if (POINT_PATCH_PID (corners[i], n) ==
825 POINT_PATCH_PID (corners[j], m))
826 {
827 LINE_PATCH_PID (thePatch, k) =
828 POINT_PATCH_PID (corners[i], n);
829 LINE_PATCH_CID0 (thePatch, k) =
830 POINT_PATCH_CID (corners[i], n);
831 LINE_PATCH_CID1 (thePatch, k) =
832 POINT_PATCH_CID (corners[j], m);
833 if (PATCH_IS_FREE (sides[LINE_PATCH_PID (thePatch, k)]))
834 freePatches++;
835 k++;
836 }
837 LINE_PATCH_N (thePatch) = k;
838
839 for (n = 0; n < LINE_PATCH_N (thePatch); n++)
840 PRINTDEBUG (dom, 1, (" pid %d cid %d %d",
841 LINE_PATCH_PID (thePatch, n),
842 LINE_PATCH_CID0 (thePatch, n),
843 LINE_PATCH_CID1 (thePatch, n)));
844 PRINTDEBUG (dom, 1, ("\n"));
845
846 IFDEBUG (dom, 10) if (k == 2)
847 {
848 INT o0, o1, s0, s1;
849
850 s0 = LINE_PATCH_PID (thePatch, 0);
851 s1 = LINE_PATCH_PID (thePatch, 1);
852 o0 =
853 (LINE_PATCH_CID0 (thePatch, 0) ==
854 ((LINE_PATCH_CID1 (thePatch, 0) + 1) % (2 * DIM_OF_BND)));
855 o1 =
856 (LINE_PATCH_CID0 (thePatch, 1) ==
857 ((LINE_PATCH_CID1 (thePatch, 1) + 1) % (2 * DIM_OF_BND)));
858 if (o0 != o1)
859 {
860 if ((PARAM_PATCH_LEFT (sides[s0]) !=
861 PARAM_PATCH_LEFT (sides[s1]))
862 || ((PARAM_PATCH_RIGHT (sides[s0]) !=
863 PARAM_PATCH_RIGHT (sides[s1]))))
864 {
865 PRINTDEBUG (dom, 0, ("patch %d and patch %d:"
866 "orientation not maches\n", s0, s1));
867 (*err)++;
868 }
869 }
870 else
871 {
872 if ((PARAM_PATCH_LEFT (sides[s0]) !=
873 PARAM_PATCH_RIGHT (sides[s1]))
874 || ((PARAM_PATCH_RIGHT (sides[s0]) !=
875 PARAM_PATCH_LEFT (sides[s1]))))
876 {
877 PRINTDEBUG (dom, 0, ("patch %d and patch %d:"
878 "orientation not maches\n", s0, s1));
879 (*err)++;
880 }
881 }
882 }
883 ENDDEBUG if (freePatches == k)
884 PATCH_STATE (thePatch) = PATCH_FREE;
885 else if (freePatches == 0)
886 PATCH_STATE (thePatch) = PATCH_FIXED;
887 else
888 PATCH_STATE (thePatch) = PATCH_BND_OF_FREE;
889 lines[(*nlines)++] = thePatch;
890 PRINTDEBUG (dom, 1, ("lines id %d type %d n %d\n",
891 PATCH_ID (thePatch), PATCH_TYPE (thePatch),
892 LINE_PATCH_N (thePatch)));
893
894
895 }
896 #endif
897
898 BVP *NS_DIM_PREFIX
BVP_Init(const char * name,HEAP * Heap,MESH * Mesh,INT MarkKey)899 BVP_Init (const char *name, HEAP * Heap, MESH * Mesh, INT MarkKey)
900 {
901 STD_BVP *theBVP;
902 DOMAIN *theDomain;
903 PROBLEM *theProblem;
904 BOUNDARY_SEGMENT *theSegment;
905 LINEAR_SEGMENT *theLinSegment;
906 BOUNDARY_CONDITION *theBndCond;
907 PATCH **corners, **sides, *thePatch;
908 unsigned short* segmentsPerPoint, *freeSegmentsPerPoint, *cornerCounters;
909 INT i, j, n, m, maxSubDomains, ncorners, nlines, nsides;
910 # ifdef __THREEDIM__
911 PATCH **lines;
912 INT err;
913 INT nn;
914 # endif
915
916 theBVP = (STD_BVP *) BVP_GetByName (name);
917 if (theBVP == NULL)
918 return (NULL);
919 currBVP = theBVP;
920
921 theDomain = theBVP->Domain;
922 if (theDomain == NULL)
923 return (NULL);
924 theProblem = theBVP->Problem;
925
926 /* fill in data of domain */
927 ncorners = theDomain->numOfCorners;
928 nsides = theDomain->numOfSegments;
929
930 /* create parameter patches */
931 maxSubDomains = 0;
932 sides = (PATCH **) GetTmpMem (Heap, nsides * sizeof (PATCH *), MarkKey);
933 if (sides == NULL)
934 return (NULL);
935
936 for (i = 0; i < nsides; i++)
937 sides[i] = NULL;
938 theBVP->nsides = nsides;
939 for (theSegment = GetFirstBoundarySegment (theDomain); theSegment != NULL;
940 theSegment = GetNextBoundarySegment (theSegment))
941 {
942 if ((theSegment->id < 0) || (theSegment->id >= nsides))
943 return (NULL);
944 thePatch = (PATCH *) GetFreelistMemory (Heap, sizeof (PARAMETER_PATCH));
945 if (thePatch == NULL)
946 return (NULL);
947 PATCH_TYPE (thePatch) = PARAMETRIC_PATCH_TYPE;
948 PATCH_ID (thePatch) = theSegment->id;
949 if (theSegment->segType == FREE)
950 PATCH_STATE (thePatch) = PATCH_FREE;
951 else
952 PATCH_STATE (thePatch) = PATCH_FIXED;
953 PARAM_PATCH_LEFT (thePatch) = theSegment->left;
954 PARAM_PATCH_RIGHT (thePatch) = theSegment->right;
955 PARAM_PATCH_BC (thePatch) = NULL;
956 PARAM_PATCH_BCD (thePatch) = NULL;
957 for (i = 0; i < 2 * DIM_OF_BND; i++)
958 PARAM_PATCH_POINTS (thePatch, i) = theSegment->points[i];
959 for (i = 0; i < DIM_OF_BND; i++)
960 {
961 PARAM_PATCH_RANGE (thePatch)[0][i] = theSegment->alpha[i];
962 PARAM_PATCH_RANGE (thePatch)[1][i] = theSegment->beta[i];
963 }
964 PARAM_PATCH_BS (thePatch) = theSegment->BndSegFunc;
965 PARAM_PATCH_BSD (thePatch) = theSegment->data;
966 maxSubDomains = MAX (maxSubDomains, theSegment->left);
967 maxSubDomains = MAX (maxSubDomains, theSegment->right);
968 sides[theSegment->id] = thePatch;
969 PRINTDEBUG (dom, 1, ("sides id %d type %d left %d right %d\n",
970 PATCH_ID (thePatch), PATCH_TYPE (thePatch),
971 PARAM_PATCH_LEFT (thePatch),
972 PARAM_PATCH_RIGHT (thePatch)));
973 for (i = 0; i < 2 * DIM_OF_BND; i++)
974 PRINTDEBUG (dom, 1,
975 (" corners %d", PARAM_PATCH_POINTS (thePatch, i)));
976 PRINTDEBUG (dom, 1, ("\n"));
977 }
978 for (theLinSegment = GetFirstLinearSegment (theDomain);
979 theLinSegment != NULL;
980 theLinSegment = GetNextLinearSegment (theLinSegment))
981 {
982 if ((theLinSegment->id < 0) || (theLinSegment->id >= nsides))
983 return (NULL);
984 thePatch = (PATCH *) GetFreelistMemory (Heap, sizeof (LINEAR_PATCH));
985 if (thePatch == NULL)
986 return (NULL);
987 PATCH_TYPE (thePatch) = LINEAR_PATCH_TYPE;
988 PATCH_ID (thePatch) = theLinSegment->id;
989 LINEAR_PATCH_LEFT (thePatch) = theLinSegment->left;
990 LINEAR_PATCH_RIGHT (thePatch) = theLinSegment->right;
991 LINEAR_PATCH_N (thePatch) = theLinSegment->n;
992 for (j = 0; j < theLinSegment->n; j++)
993 {
994 LINEAR_PATCH_POINTS (thePatch, j) = theLinSegment->points[j];
995 for (i = 0; i < DIM; i++)
996 LINEAR_PATCH_POS (thePatch, j)[i] = theLinSegment->x[j][i];
997 }
998 maxSubDomains = MAX (maxSubDomains, theLinSegment->left);
999 maxSubDomains = MAX (maxSubDomains, theLinSegment->right);
1000 sides[theLinSegment->id] = thePatch;
1001 PRINTDEBUG (dom, 1, ("sides id %d type %d left %d right %d\n",
1002 PATCH_ID (thePatch), PATCH_TYPE (thePatch),
1003 LINEAR_PATCH_LEFT (thePatch),
1004 LINEAR_PATCH_RIGHT (thePatch)));
1005 /** \todo why this here??? (CVS-merge mess-up?) */
1006 if (theProblem != NULL)
1007 {
1008 UserWrite ("Use CreateBoundaryValueProblem!");
1009 return (NULL);
1010 }
1011 }
1012 theBVP->numOfSubdomains = maxSubDomains;
1013 PRINTDEBUG (dom, 1, (" bvp nsubcf %x\n", theBVP->numOfSubdomains));
1014 for (i = 0; i < nsides; i++)
1015 if (sides[i] == NULL)
1016 return (NULL);
1017
1018 if (theProblem != NULL)
1019 for (theBndCond = GetFirstBoundaryCondition (theProblem);
1020 theBndCond != NULL;
1021 theBndCond = GetNextBoundaryCondition (theBndCond))
1022 {
1023 i = theBndCond->id;
1024 if ((i < 0) || (i >= nsides))
1025 return (NULL);
1026 thePatch = sides[i];
1027 PARAM_PATCH_BC (thePatch) = theBndCond->BndCond;
1028 PARAM_PATCH_BCD (thePatch) = theBndCond->data;
1029 }
1030
1031 /* create point patches */
1032 corners = (PATCH **) GetTmpMem (Heap, ncorners * sizeof (PATCH *), MarkKey);
1033 if (corners == NULL)
1034 return (NULL);
1035 theBVP->ncorners = ncorners;
1036
1037 /* precompute the number of segments at each point patch */
1038 segmentsPerPoint = (unsigned short*)calloc(ncorners,sizeof(unsigned short));
1039 freeSegmentsPerPoint = (unsigned short*)calloc(ncorners,sizeof(unsigned short));
1040
1041 for (j = 0; j < nsides; j++) {
1042
1043 if (PATCH_TYPE (sides[j]) == LINEAR_PATCH_TYPE) {
1044 for (n = 0; n < LINEAR_PATCH_N (sides[j]); n++)
1045 segmentsPerPoint[LINEAR_PATCH_POINTS (sides[j], n)]++;
1046
1047 if (PATCH_IS_FREE(sides[j]))
1048 for (n = 0; n < LINEAR_PATCH_N (sides[j]); n++)
1049 freeSegmentsPerPoint[LINEAR_PATCH_POINTS (sides[j], n)]++;
1050
1051 } else if (PATCH_TYPE (sides[j]) == PARAMETRIC_PATCH_TYPE) {
1052 for (n = 0; n < 2 * DIM_OF_BND; n++)
1053 if (PARAM_PATCH_POINTS (sides[j], n) >= 0) /* Entry may be -1 for triangular faces */
1054 segmentsPerPoint[PARAM_PATCH_POINTS (sides[j], n)]++;
1055
1056 if (PATCH_IS_FREE(sides[j]))
1057 for (n = 0; n < 2 * DIM_OF_BND; n++)
1058 if (PARAM_PATCH_POINTS (sides[j], n) >= 0) /* Entry may be -1 for triangular faces */
1059 freeSegmentsPerPoint[PARAM_PATCH_POINTS (sides[j], n)]++;
1060
1061 }
1062
1063 }
1064
1065 /* Allocate point patches */
1066 for (i = 0; i < ncorners; i++) {
1067
1068 m = segmentsPerPoint[i];
1069
1070 thePatch =
1071 (PATCH *) GetFreelistMemory (Heap, sizeof (POINT_PATCH)
1072 + (m-1) * sizeof (struct point_on_patch));
1073 if (thePatch == NULL)
1074 return (NULL);
1075
1076 PATCH_TYPE (thePatch) = POINT_PATCH_TYPE;
1077 PATCH_ID (thePatch) = i;
1078
1079 POINT_PATCH_N (thePatch) = m;
1080 corners[i] = thePatch;
1081
1082 }
1083
1084 cornerCounters = (unsigned short*)calloc(ncorners,sizeof(unsigned short));
1085
1086 for (j = 0; j < nsides; j++)
1087 {
1088 if (PATCH_TYPE (sides[j]) == LINEAR_PATCH_TYPE)
1089 {
1090 for (n = 0; n < LINEAR_PATCH_N (sides[j]); n++)
1091 {
1092 i = LINEAR_PATCH_POINTS (sides[j], n);
1093 POINT_PATCH_PID (corners[i], cornerCounters[i]) = j;
1094 POINT_PATCH_CID (corners[i], cornerCounters[i]++) = n;
1095 }
1096 }
1097 else if (PATCH_TYPE (sides[j]) == PARAMETRIC_PATCH_TYPE)
1098 {
1099 for (n = 0; n < 2 * DIM_OF_BND; n++)
1100 {
1101 i = PARAM_PATCH_POINTS (sides[j], n);
1102 if (i>=0 && i<ncorners) {
1103 POINT_PATCH_PID (corners[i], cornerCounters[i]) = j;
1104 POINT_PATCH_CID (corners[i], cornerCounters[i]++) = n;
1105 }
1106 }
1107 }
1108 }
1109
1110 for (i = 0; i < ncorners; i++) {
1111 /* Shouldn't the 'thePatch' be replaced 'corners[i]' in this loop? OS */
1112 if (freeSegmentsPerPoint[i] == cornerCounters[i])
1113 PATCH_STATE (thePatch) = PATCH_FREE;
1114 else if (freeSegmentsPerPoint[i] == 0)
1115 PATCH_STATE (thePatch) = PATCH_FIXED;
1116 else
1117 PATCH_STATE (thePatch) = PATCH_BND_OF_FREE;
1118
1119 PRINTDEBUG (dom, 1, ("corners id %d type %d n %d\n",
1120 PATCH_ID (thePatch), PATCH_TYPE (thePatch),
1121 POINT_PATCH_N (thePatch)));
1122 for (j = 0; j < POINT_PATCH_N (thePatch); j++)
1123 PRINTDEBUG (dom, 1, (" pid %d cid %d\n",
1124 POINT_PATCH_PID (thePatch, j),
1125 POINT_PATCH_CID (thePatch, j)));
1126 }
1127
1128 /* Return memory that will not be used anymore */
1129 free(segmentsPerPoint);
1130 free(freeSegmentsPerPoint);
1131 free(cornerCounters);
1132
1133 /* create line patches */
1134 nlines = 0;
1135 #ifdef __THREEDIM__
1136 /* The maximum number of boundary lines is equal to the number of
1137 boundary faces (nsides) times the max number of lines per face (=4)
1138 divided by two. */
1139 lines =
1140 (PATCH **) GetTmpMem (Heap, nsides * 2 * sizeof (PATCH *),
1141 MarkKey);
1142 if (lines == NULL)
1143 return (NULL);
1144 err = 0;
1145
1146 /* We create the set of all boundary lines by looping over the sides
1147 and for each side loop over the edges of this side. That way, we
1148 meet each boundary line twice. We cannot assume
1149 that the sides are properly oriented. Thus we introduce a c++-std::set
1150 which saves the index pair of the nodes defining the line.
1151 When a line is visited the second time, the pair is expunged from the
1152 set in order to save memory.
1153 */
1154 typedef std::set<std::pair<long,long> > set;
1155 set bnd_edges;
1156
1157 for (int s=0; s<nsides; s++) {
1158
1159 if (PATCH_TYPE (sides[s]) == LINEAR_PATCH_TYPE) {
1160
1161 for (nn = 0; nn < LINEAR_PATCH_N (sides[s]); nn++)
1162 {
1163 i = LINEAR_PATCH_POINTS (sides[s], nn);
1164 j = LINEAR_PATCH_POINTS (sides[s], (nn+1)%LINEAR_PATCH_N(sides[s]));
1165
1166 long max = std::max(i,j);
1167 long min = i + j - max;
1168 std::pair<long, long> z(min,max);
1169
1170 if ( bnd_edges.insert(z).second ) /* true if bnd_edges didn't contain z yet */
1171 {
1172 /* Insert the line into the boundary data structure */
1173 CreateLine(min, max, Heap, thePatch, corners, lines, sides, &nlines, &err);
1174 }
1175 else
1176 bnd_edges.erase(z);
1177
1178 }
1179
1180 } else if (PATCH_TYPE (sides[s]) == PARAMETRIC_PATCH_TYPE) {
1181
1182 /* Handle edges of parametric patches */
1183
1184 /* Determine whether the boundary segment is a triangle or a quadrilateral.
1185 We assume that it is a triangle if the fourth vertex is invalid, and
1186 a quadrilateral otherwise. This is the best we can do.
1187 */
1188
1189 int cornersOfParametricPatch = (PARAM_PATCH_POINTS (sides[s], 3) < 0 || PARAM_PATCH_POINTS (sides[s], 3) > ncorners) ? 3 : 4;
1190
1191 for (nn = 0; nn < cornersOfParametricPatch; nn++)
1192 {
1193 i = PARAM_PATCH_POINTS (sides[s], nn);
1194 j = PARAM_PATCH_POINTS (sides[s], (nn+1)%cornersOfParametricPatch);
1195
1196 long max = std::max(i,j);
1197 long min = i + j - max;
1198 std::pair<long, long> z(min,max);
1199
1200 if ( bnd_edges.insert(z).second ) /* true if bnd_edges didn't contain z yet */
1201 {
1202 /* Insert the line into the boundary data structure */
1203 CreateLine(min, max, Heap, thePatch, corners, lines, sides, &nlines, &err);
1204 }
1205 else
1206 bnd_edges.erase(z);
1207
1208 }
1209
1210 } else
1211 UserWrite("Error: unknown PATCH_TYPE found for a boundary side!\n");
1212
1213 }
1214 ASSERT (err == 0);
1215 #endif
1216
1217 m = ncorners + nlines;
1218 theBVP->sideoffset = m;
1219 n = m + nsides;
1220 theBVP->patches = (PATCH **) GetFreelistMemory (Heap, n * sizeof (PATCH *));
1221 n = 0;
1222 for (i = 0; i < ncorners; i++)
1223 {
1224 thePatch = corners[i];
1225 for (j = 0; j < POINT_PATCH_N (thePatch); j++)
1226 POINT_PATCH_PID (thePatch, j) += m;
1227 theBVP->patches[n++] = thePatch;
1228 }
1229 #ifdef __THREEDIM__
1230 for (i = 0; i < nlines; i++)
1231 {
1232 thePatch = lines[i];
1233 PATCH_ID (thePatch) = n;
1234 for (j = 0; j < LINE_PATCH_N (thePatch); j++)
1235 LINE_PATCH_PID (thePatch, j) += m;
1236 theBVP->patches[n++] = thePatch;
1237 }
1238 #endif
1239 for (i = 0; i < nsides; i++)
1240 {
1241 thePatch = sides[i];
1242 PATCH_ID (thePatch) = n;
1243 theBVP->patches[n++] = thePatch;
1244 }
1245
1246 PRINTDEBUG (dom, 1, ("ncorners %d\n", theBVP->ncorners));
1247 for (i = 0; i < theBVP->ncorners; i++)
1248 PRINTDEBUG (dom, 1, (" id %d\n", PATCH_ID (theBVP->patches[i])));
1249
1250 if (Mesh != NULL)
1251 {
1252 Mesh->mesh_status = MESHSTAT_CNODES;
1253 Mesh->nBndP = theBVP->ncorners;
1254 Mesh->nInnP = 0;
1255 Mesh->nElements = NULL;
1256 Mesh->VertexLevel = NULL;
1257 Mesh->VertexPrio = NULL;
1258 Mesh->ElementLevel = NULL;
1259 Mesh->ElementPrio = NULL;
1260 Mesh->ElemSideOnBnd = NULL;
1261 Mesh->theBndPs =
1262 (BNDP **) GetTmpMem (Heap, n * sizeof (BNDP *), MarkKey);
1263 if (Mesh->theBndPs == NULL)
1264 return (NULL);
1265
1266 if (CreateCornerPoints (Heap, theBVP, Mesh->theBndPs))
1267 return (NULL);
1268
1269 PRINTDEBUG (dom, 1, ("mesh n %d\n", Mesh->nBndP));
1270 for (i = 0; i < theBVP->ncorners; i++)
1271 PRINTDEBUG (dom, 1, (" id %d\n",
1272 BND_PATCH_ID ((BND_PS *) (Mesh->theBndPs[i]))));
1273 }
1274
1275 /* allocate s2p table */
1276 STD_BVP_NDOMPART (theBVP) = DOMAIN_NPARTS (theDomain);
1277 STD_BVP_S2P_PTR (theBVP) =
1278 (INT *) GetFreelistMemory (Heap,
1279 (1 + STD_BVP_NSUBDOM (theBVP)) * sizeof (INT));
1280 if (STD_BVP_S2P_PTR (theBVP) == NULL)
1281 return (NULL);
1282
1283 /* fill number of parts */
1284 if (DOMAIN_NPARTS (theDomain) > 1)
1285 {
1286 const DOMAIN_PART_INFO *dpi;
1287
1288 if (DOMAIN_NPARTS (theDomain) > (1 << VPART_LEN))
1289 {
1290 UserWriteF("Too many parts for control entry in vector\n");
1291 UserWriteF("Domain requests %d parts, but only %d are possible!\n",
1292 DOMAIN_NPARTS (theDomain), (1 << VPART_LEN));
1293 ASSERT (false);
1294 return (NULL);
1295 }
1296
1297 /* transfer from part info (NB: STD_BVP_NSUBDOM only counts inner subdomains) */
1298 dpi = DOMAIN_PARTINFO (theDomain);
1299 for (i = 0; i <= STD_BVP_NSUBDOM (theBVP); i++)
1300 STD_BVP_S2P (theBVP, i) = DPI_SD2P (dpi, i);
1301 }
1302 else
1303 {
1304 /* 0 for each subdomnain by default */
1305 for (i = 0; i < STD_BVP_NSUBDOM (theBVP); i++)
1306 STD_BVP_S2P (theBVP, i) = 0;
1307 }
1308
1309 return ((BVP *) theBVP);
1310 }
1311
1312 /* domain interface function: for description see domain.h */
1313 INT NS_DIM_PREFIX
BVP_Dispose(BVP * theBVP)1314 BVP_Dispose (BVP * theBVP)
1315 {
1316 /* Deallocate the patch pointers
1317 The memory pointed to by theBVP->patches should also be freed when
1318 not using the system heap. However the UG heap data structure is not
1319 available here, and for this I don't know how to do the proper deallocation. */
1320 STD_BVP* stdBVP = (STD_BVP *) theBVP;
1321 /* npatches is the number of corners plus the number of lines plus the number of sides.
1322 * You apparently can't access nlines directly here, but sideoffset should be ncorners + nlines. */
1323 int npatches = stdBVP->sideoffset + stdBVP->nsides;
1324 for (int i=0; i<npatches; i++)
1325 free ( stdBVP->patches[i] );
1326
1327 free ( stdBVP->patches );
1328
1329 free ( STD_BVP_S2P_PTR (stdBVP) );
1330
1331 /* Unlock the item so it can be deleted from the environment tree */
1332 ((ENVITEM*)theBVP)->d.locked = 0;
1333
1334 if (ChangeEnvDir("/BVP")==NULL)
1335 return (1);
1336 if (RemoveEnvItem((ENVITEM *)theBVP))
1337 return (1);
1338
1339 return (0);
1340 }
1341
1342 void NS_DIM_PREFIX
Set_Current_BVP(BVP * theBVP)1343 Set_Current_BVP(BVP* theBVP)
1344 {
1345 currBVP = (STD_BVP*)theBVP;
1346 }
1347
1348 /* domain interface function: for description see domain.h */
1349 BVP *NS_DIM_PREFIX
BVP_GetByName(const char * name)1350 BVP_GetByName (const char *name)
1351 {
1352 return ((BVP *) SearchEnv (name, "/BVP", theBVPDirID, theBVPDirID));
1353 }
1354
1355 INT NS_DIM_PREFIX
BVP_SetBVPDesc(BVP * aBVP,BVP_DESC * theBVPDesc)1356 BVP_SetBVPDesc (BVP * aBVP, BVP_DESC * theBVPDesc)
1357 {
1358 STD_BVP *theBVP;
1359
1360 if (aBVP == NULL)
1361 return (1);
1362
1363 /* cast */
1364 theBVP = GetSTD_BVP (aBVP);
1365
1366 /* general part */
1367 strcpy (BVPD_NAME (theBVPDesc), ENVITEM_NAME (theBVP));
1368
1369 /* the domain part */
1370 BVPD_NSUBDOM (theBVPDesc) = theBVP->numOfSubdomains;
1371 BVPD_NPARTS (theBVPDesc) = theBVP->nDomainParts;
1372 BVPD_S2P_PTR (theBVPDesc) = STD_BVP_S2P_PTR (theBVP);
1373 BVPD_NCOEFFF (theBVPDesc) = theBVP->numOfCoeffFct;
1374 BVPD_NUSERF (theBVPDesc) = theBVP->numOfUserFct;
1375 BVPD_CONFIG (theBVPDesc) = theBVP->ConfigProc;
1376
1377 currBVP = theBVP;
1378
1379 return (0);
1380 }
1381
1382 /* domain interface function: for description see domain.h */
1383 INT NS_DIM_PREFIX
BVP_SetCoeffFct(BVP * aBVP,INT n,CoeffProcPtr * CoeffFct)1384 BVP_SetCoeffFct (BVP * aBVP, INT n, CoeffProcPtr * CoeffFct)
1385 {
1386 STD_BVP *theBVP;
1387 INT i;
1388
1389 theBVP = GetSTD_BVP (aBVP);
1390
1391 /* check */
1392 if (n < -1 || n >= theBVP->numOfCoeffFct)
1393 return (1);
1394
1395 if (n == -1)
1396 for (i = 0; i < theBVP->numOfCoeffFct; i++)
1397 CoeffFct[i] = (CoeffProcPtr) theBVP->CU_ProcPtr[i];
1398 else
1399 CoeffFct[0] = (CoeffProcPtr) theBVP->CU_ProcPtr[n];
1400
1401 return (0);
1402 }
1403
1404 /* domain interface function: for description see domain.h */
1405 INT NS_DIM_PREFIX
BVP_SetUserFct(BVP * aBVP,INT n,UserProcPtr * UserFct)1406 BVP_SetUserFct (BVP * aBVP, INT n, UserProcPtr * UserFct)
1407 {
1408 STD_BVP *theBVP;
1409 INT i;
1410
1411 theBVP = GetSTD_BVP (aBVP);
1412
1413 /* check */
1414 if (n < -1 || n >= theBVP->numOfUserFct)
1415 return (1);
1416
1417 if (n == -1)
1418 for (i = 0; i < theBVP->numOfUserFct; i++)
1419 UserFct[i] =
1420 (UserProcPtr) theBVP->CU_ProcPtr[i + theBVP->numOfCoeffFct];
1421 else
1422 UserFct[0] = (UserProcPtr) theBVP->CU_ProcPtr[n + theBVP->numOfCoeffFct];
1423
1424 return (0);
1425 }
1426
1427 /* domain interface function: for description see domain.h */
1428 INT NS_DIM_PREFIX
BVP_Check(BVP * aBVP)1429 BVP_Check (BVP * aBVP)
1430 {
1431 UserWrite ("BVP_Check: not implemented\n");
1432
1433 return (0);
1434 }
1435
1436 static INT
GetNumberOfCommonPatches(PATCH * p0,PATCH * p1,INT * Pid)1437 GetNumberOfCommonPatches (PATCH * p0, PATCH * p1, INT * Pid)
1438 {
1439 INT i, j, cnt, np0, np1, pid;
1440
1441 cnt = 0;
1442 np0 = GetNumberOfPatches (p0);
1443 np1 = GetNumberOfPatches (p1);
1444 for (i = 0; i < np0; i++)
1445 {
1446 pid = GetPatchId (p0, i);
1447 for (j = 0; j < np1; j++)
1448 if (pid == GetPatchId (p1, j))
1449 {
1450 if (!cnt)
1451 *Pid = pid;
1452 cnt++;
1453 }
1454 }
1455
1456 return (cnt);
1457 }
1458
1459 #ifdef __THREEDIM__
1460 static INT
GetCommonPatchId(PATCH * p0,PATCH * p1,INT k)1461 GetCommonPatchId (PATCH * p0, PATCH * p1, INT k)
1462 {
1463 INT i, j, cnt;
1464
1465 cnt = 0;
1466 for (i = 0; i < GetNumberOfPatches (p0); i++)
1467 for (j = 0; j < GetNumberOfPatches (p1); j++)
1468 if (GetPatchId (p0, i) == GetPatchId (p1, j))
1469 if (k == cnt++)
1470 return (GetPatchId (p1, j));
1471
1472 return (-1);
1473 }
1474
1475 static INT
GetCommonLinePatchId(PATCH * p0,PATCH * p1)1476 GetCommonLinePatchId (PATCH * p0, PATCH * p1)
1477 {
1478 INT i, k, l, cnt, cnt1;
1479 PATCH *p;
1480
1481 if (PATCH_TYPE (p0) == LINE_PATCH_TYPE)
1482 return (PATCH_ID (p0));
1483 else if (PATCH_TYPE (p1) == LINE_PATCH_TYPE)
1484 return (PATCH_ID (p1));
1485
1486 cnt = GetNumberOfCommonPatches (p0, p1, &i);
1487
1488 if (cnt < 1)
1489 return (-1);
1490
1491 for (k = currBVP->ncorners; k < currBVP->sideoffset; k++)
1492 {
1493 p = currBVP->patches[k];
1494 if (LINE_PATCH_N (p) != cnt)
1495 continue;
1496 cnt1 = 0;
1497 for (i = 0; i < cnt; i++)
1498 for (l = 0; l < LINE_PATCH_N (p); l++)
1499 if (GetCommonPatchId (p0, p1, i) == LINE_PATCH_PID (p, l))
1500 cnt1++;
1501 if (cnt == cnt1)
1502 return (k);
1503 }
1504
1505 return (-1);
1506 }
1507
1508 static BNDP *
CreateBndPOnLine(HEAP * Heap,PATCH * p0,PATCH * p1,DOUBLE lcoord)1509 CreateBndPOnLine (HEAP * Heap, PATCH * p0, PATCH * p1, DOUBLE lcoord)
1510 {
1511 BND_PS *bp;
1512 PATCH *p, *pp;
1513 DOUBLE local0[DIM_OF_BND];
1514 DOUBLE local1[DIM_OF_BND];
1515 INT k, l, cnt;
1516
1517 if (PATCH_TYPE (p0) != POINT_PATCH_TYPE)
1518 return (NULL);
1519 if (PATCH_TYPE (p1) != POINT_PATCH_TYPE)
1520 return (NULL);
1521
1522 PRINTDEBUG (dom, 1, (" p0 p1 %d %d\n", PATCH_ID (p0), PATCH_ID (p1)));
1523 for (l = 0; l < GetNumberOfPatches (p0); l++)
1524 PRINTDEBUG (dom, 1, (" bp pid %d\n", GetPatchId (p0, l)));
1525 for (l = 0; l < GetNumberOfPatches (p1); l++)
1526 PRINTDEBUG (dom, 1, (" bp pid %d\n", GetPatchId (p1, l)));
1527
1528 cnt = GetNumberOfCommonPatches (p0, p1, &k);
1529
1530 if (cnt < 2)
1531 return (NULL);
1532
1533 bp =
1534 (BND_PS *) GetFreelistMemory (Heap,
1535 (cnt - 1) * sizeof (COORD_BND_VECTOR) +
1536 sizeof (BND_PS));
1537 if (bp == NULL)
1538 return (NULL);
1539 bp->n = cnt;
1540
1541 k = GetCommonLinePatchId (p0, p1);
1542 if ((k < currBVP->ncorners) || (k >= currBVP->sideoffset))
1543 return (NULL);
1544 p = currBVP->patches[k];
1545 bp->patch_id = k;
1546
1547 PRINTDEBUG (dom, 1, (" Create BNDP line %d cnt %d\n", k, cnt));
1548 for (l = 0; l < GetNumberOfPatches (p0); l++)
1549 PRINTDEBUG (dom, 1, (" bp pid %d\n", GetPatchId (p0, l)));
1550 for (l = 0; l < GetNumberOfPatches (p1); l++)
1551 PRINTDEBUG (dom, 1, (" bp pid %d\n", GetPatchId (p1, l)));
1552
1553 for (l = 0; l < LINE_PATCH_N (p); l++)
1554 {
1555 pp = currBVP->patches[LINE_PATCH_PID (p, l)];
1556 switch (LINE_PATCH_CID0 (p, l))
1557 {
1558 case 0 :
1559 local0[0] = PARAM_PATCH_RANGE (pp)[0][0];
1560 local0[1] = PARAM_PATCH_RANGE (pp)[0][1];
1561 break;
1562 case 1 :
1563 local0[0] = PARAM_PATCH_RANGE (pp)[1][0];
1564 local0[1] = PARAM_PATCH_RANGE (pp)[0][1];
1565 break;
1566 case 2 :
1567 local0[0] = PARAM_PATCH_RANGE (pp)[1][0];
1568 local0[1] = PARAM_PATCH_RANGE (pp)[1][1];
1569 break;
1570 case 3 :
1571 local0[0] = PARAM_PATCH_RANGE (pp)[0][0];
1572 local0[1] = PARAM_PATCH_RANGE (pp)[1][1];
1573 break;
1574 }
1575 switch (LINE_PATCH_CID1 (p, l))
1576 {
1577 case 0 :
1578 local1[0] = PARAM_PATCH_RANGE (pp)[0][0];
1579 local1[1] = PARAM_PATCH_RANGE (pp)[0][1];
1580 break;
1581 case 1 :
1582 local1[0] = PARAM_PATCH_RANGE (pp)[1][0];
1583 local1[1] = PARAM_PATCH_RANGE (pp)[0][1];
1584 break;
1585 case 2 :
1586 local1[0] = PARAM_PATCH_RANGE (pp)[1][0];
1587 local1[1] = PARAM_PATCH_RANGE (pp)[1][1];
1588 break;
1589 case 3 :
1590 local1[0] = PARAM_PATCH_RANGE (pp)[0][0];
1591 local1[1] = PARAM_PATCH_RANGE (pp)[1][1];
1592 break;
1593 }
1594
1595 /* TODO: Why this? */
1596 if ((LINE_PATCH_CID0 (p, l) == 2) && (LINE_PATCH_CID1 (p, l) == 3))
1597 lcoord = 1.0 - lcoord;
1598 if ((LINE_PATCH_CID0 (p, l) == 1) && (LINE_PATCH_CID1 (p, l) == 0))
1599 lcoord = 1.0 - lcoord;
1600 if ((LINE_PATCH_CID0 (p, l) == 3) && (LINE_PATCH_CID1 (p, l) == 0))
1601 lcoord = 1.0 - lcoord;
1602 if ((LINE_PATCH_CID0 (p, l) == 2) && (LINE_PATCH_CID1 (p, l) == 1))
1603 lcoord = 1.0 - lcoord;
1604
1605 bp->local[l][0] = (1.0 - lcoord) * local0[0] + lcoord * local1[0];
1606 bp->local[l][1] = (1.0 - lcoord) * local0[1] + lcoord * local1[1];
1607
1608 if ((LINE_PATCH_CID0 (p, l) == 2) && (LINE_PATCH_CID1 (p, l) == 3))
1609 lcoord = 1.0 - lcoord;
1610 if ((LINE_PATCH_CID0 (p, l) == 1) && (LINE_PATCH_CID1 (p, l) == 0))
1611 lcoord = 1.0 - lcoord;
1612 if ((LINE_PATCH_CID0 (p, l) == 3) && (LINE_PATCH_CID1 (p, l) == 0))
1613 lcoord = 1.0 - lcoord;
1614 if ((LINE_PATCH_CID0 (p, l) == 2) && (LINE_PATCH_CID1 (p, l) == 1))
1615 lcoord = 1.0 - lcoord;
1616
1617
1618 PRINTDEBUG (dom, 1,
1619 (" Create bndp %d line %d C0 %d C1 %d l %d %f %f\n",
1620 bp->patch_id, LINE_PATCH_PID (p, l), LINE_PATCH_CID0 (p,
1621 l),
1622 LINE_PATCH_CID1 (p, l), l, bp->local[l][0],
1623 bp->local[l][1]));
1624 PRINTDEBUG (dom, 1, (" lcoord %f\n", lcoord));
1625 }
1626
1627 if (!PATCH_IS_FIXED (p))
1628 {
1629 /* store global coordinates */
1630 BND_DATA (bp) = GetFreelistMemory (Heap, DIM * sizeof (DOUBLE));
1631 if (BND_DATA (bp) == NULL)
1632 return (NULL);
1633
1634 if (BndPointGlobal ((BNDP *) bp, (DOUBLE *) BND_DATA (bp)))
1635 return (NULL);
1636 }
1637
1638 return ((BNDP *) bp);
1639 }
1640 #endif
1641
1642 #define BN_RES 100
1643
1644 static int
DropPerpendicularOnSegment(PATCH * patch,double range[][DIM_OF_BND],const double global[],double local[],double * mindist2)1645 DropPerpendicularOnSegment (PATCH * patch, double range[][DIM_OF_BND],
1646 const double global[], double local[],
1647 double *mindist2)
1648 {
1649 double sa = range[0][0];
1650 double sb = range[1][0];
1651 double ds = (sb - sa) / ((double) BN_RES);
1652 int k;
1653 #if (DIM==3)
1654 double ta = range[0][1];
1655 double tb = range[1][1];
1656 double dt = (tb - ta) / ((double) BN_RES);
1657 int l;
1658 #endif
1659
1660 for (k = 0; k <= BN_RES; k++)
1661 {
1662 DOUBLE lambda[DIM_OF_BND];
1663 DOUBLE_VECTOR point, diff;
1664 double dist2;
1665
1666 lambda[0] = (k == BN_RES) ? sb : sa + k * ds;
1667
1668 #if (DIM==3)
1669 for (l = 0; l <= BN_RES; l++)
1670 {
1671
1672 lambda[1] = (l == BN_RES) ? tb : ta + l * dt;
1673 #endif
1674
1675 if (PatchGlobal (patch, lambda, point))
1676 return 1;
1677 V_DIM_SUBTRACT (point, global, diff);
1678 dist2 = V_DIM_SCAL_PROD (diff, diff);
1679 if (*mindist2 > dist2)
1680 {
1681 *mindist2 = dist2;
1682 V_BDIM_COPY (lambda, local);
1683 }
1684 #if (DIM==3)
1685 }
1686 #endif
1687 }
1688 return 0;
1689 }
1690
1691 static int
ResolvePointOnSegment(PATCH * patch,int depth,double resolution2,double range[][DIM_OF_BND],const double global[],double local[])1692 ResolvePointOnSegment (PATCH * patch, int depth, double resolution2,
1693 double range[][DIM_OF_BND], const double global[],
1694 double local[])
1695 {
1696 double ds = (range[1][0] - range[0][0]) / ((double) BN_RES);
1697 #if (DIM==3)
1698 double dt = (range[1][1] - range[0][1]) / ((double) BN_RES);
1699 #endif
1700 double new_range[2][DIM_OF_BND];
1701 DOUBLE mindist2 = MAX_D;
1702
1703 new_range[0][0] = local[0] - ds;
1704 new_range[1][0] = local[0] + ds;
1705 #if (DIM==3)
1706 new_range[0][1] = local[1] - dt;
1707 new_range[1][1] = local[1] + dt;
1708 #endif
1709
1710 if (DropPerpendicularOnSegment (patch, new_range, global, local, &mindist2))
1711 return 1;
1712 if (mindist2 > resolution2)
1713 {
1714 if (depth > 0)
1715 {
1716 /* recursive call */
1717 if (ResolvePointOnSegment
1718 (patch, depth - 1, resolution2, new_range, global, local))
1719 return 1;
1720 }
1721 else
1722 return 2;
1723 }
1724 return 0;
1725 }
1726
1727 /* domain interface function: for description see domain.h */
1728 /* TODO: syntax for manpages??? */
1729 BNDP *NS_DIM_PREFIX
BVP_InsertBndP(HEAP * Heap,BVP * aBVP,INT argc,char ** argv)1730 BVP_InsertBndP (HEAP * Heap, BVP * aBVP, INT argc, char **argv)
1731 {
1732 using std::abs;
1733
1734 STD_BVP *theBVP;
1735 BND_PS *ps;
1736 PATCH *p;
1737 INT j, pid;
1738 int i;
1739 DOUBLE pos[2];
1740 # ifdef __THREEDIM__
1741 DOUBLE lc;
1742 # endif
1743
1744 theBVP = GetSTD_BVP (aBVP);
1745
1746 if (ReadArgvOption ("g", argc, argv))
1747 {
1748 DOUBLE resolution2;
1749 /* double global[DIM]; extension to constant 3 components for sscanf; Ch. Wrobel 981002 */
1750 double global[3];
1751
1752 /* insert bn from global coordinates */
1753 if (sscanf (argv[0], "bn %lf %lf %lf", global, global +1, global +2) !=
1754 DIM)
1755 {
1756 PrintErrorMessageF ('E', "BVP_InsertBndP",
1757 "g option specified but could not scan\n"
1758 "global coordinates from '%s'\n", argv[0]);
1759 return (NULL);
1760 }
1761
1762 if (ReadArgvDOUBLE ("r", &resolution2, argc, argv) != 0)
1763 resolution2 = 1e-2; /* default */
1764 resolution2 *= resolution2;
1765
1766 /* find segment id and local coordinates (on any segment) */
1767 {
1768 DOUBLE mindist2 = MAX_D;
1769 int seg;
1770 DOUBLE lambda[DIM_OF_BND];
1771
1772 for (seg = 0; seg < STD_BVP_NSIDES (theBVP); seg++)
1773 {
1774 PATCH *patch =
1775 STD_BVP_PATCH (theBVP, seg + STD_BVP_SIDEOFFSET (theBVP));
1776 DOUBLE seg_mindist2 = mindist2;
1777
1778 if (DropPerpendicularOnSegment
1779 (patch, PARAM_PATCH_RANGE (patch), global, lambda,
1780 &seg_mindist2))
1781 return NULL;
1782
1783 if (mindist2 > seg_mindist2)
1784 {
1785 mindist2 = seg_mindist2;
1786 i = seg;
1787 V_BDIM_COPY (lambda, pos);
1788 }
1789 if (mindist2 <= resolution2)
1790 break;
1791 }
1792 if (mindist2 > resolution2)
1793 {
1794 /* refine search */
1795 PATCH *patch =
1796 STD_BVP_PATCH (theBVP, i + STD_BVP_SIDEOFFSET (theBVP));
1797
1798 V_BDIM_COPY (pos, lambda);
1799 if (ResolvePointOnSegment
1800 (patch, 2, resolution2, PARAM_PATCH_RANGE (patch), global,
1801 lambda))
1802 return NULL;
1803 V_BDIM_COPY (lambda, pos);
1804 }
1805 }
1806
1807 }
1808 else
1809 {
1810 if (sscanf (argv[0], "bn %d %lf %lf", &i, pos, pos + 1) != DIM_OF_BND + 1)
1811 {
1812 PrintErrorMessageF ('E', "BVP_InsertBndP",
1813 "could not scan segment id and\n"
1814 "local coordinates on segment from '%s'\n",
1815 argv[0]);
1816 return (NULL);
1817 }
1818 }
1819 pid = i + theBVP->sideoffset;
1820 p = theBVP->patches[pid];
1821
1822 #ifdef __THREEDIM__
1823 /* check point on line or on point patch */
1824 if (abs(pos[0] - PARAM_PATCH_RANGE (p)[0][0]) < SMALL_DIFF)
1825 {
1826 lc = (pos[1] - PARAM_PATCH_RANGE (p)[0][1])
1827 / (PARAM_PATCH_RANGE (p)[1][1] - PARAM_PATCH_RANGE (p)[0][1]);
1828 if (abs(lc) < SMALL_DIFF)
1829 return (CreateBndPOnPoint
1830 (Heap, currBVP->patches[PARAM_PATCH_POINTS (p, 0)]));
1831 else if (abs(lc - 1.) < SMALL_DIFF)
1832 return (CreateBndPOnPoint
1833 (Heap, currBVP->patches[PARAM_PATCH_POINTS (p, 3)]));
1834 return (CreateBndPOnLine
1835 (Heap, currBVP->patches[PARAM_PATCH_POINTS (p, 0)],
1836 currBVP->patches[PARAM_PATCH_POINTS (p, 3)], lc));
1837 }
1838 else if (abs(pos[0] - PARAM_PATCH_RANGE (p)[1][0]) < SMALL_DIFF)
1839 {
1840 lc = (pos[1] - PARAM_PATCH_RANGE (p)[0][1])
1841 / (PARAM_PATCH_RANGE (p)[1][1] - PARAM_PATCH_RANGE (p)[0][1]);
1842 if (abs(lc) < SMALL_DIFF)
1843 return (CreateBndPOnPoint
1844 (Heap, currBVP->patches[PARAM_PATCH_POINTS (p, 1)]));
1845 else if (abs(lc - 1.) < SMALL_DIFF)
1846 return (CreateBndPOnPoint
1847 (Heap, currBVP->patches[PARAM_PATCH_POINTS (p, 2)]));
1848 return (CreateBndPOnLine
1849 (Heap, currBVP->patches[PARAM_PATCH_POINTS (p, 1)],
1850 currBVP->patches[PARAM_PATCH_POINTS (p, 2)], lc));
1851 }
1852 else if (abs(pos[1] - PARAM_PATCH_RANGE (p)[0][1]) < SMALL_DIFF)
1853 {
1854 lc = (pos[0] - PARAM_PATCH_RANGE (p)[0][0])
1855 / (PARAM_PATCH_RANGE (p)[1][0] - PARAM_PATCH_RANGE (p)[0][0]);
1856 if (abs(lc) < SMALL_DIFF)
1857 return (CreateBndPOnPoint
1858 (Heap, currBVP->patches[PARAM_PATCH_POINTS (p, 0)]));
1859 else if (abs(lc - 1.) < SMALL_DIFF)
1860 return (CreateBndPOnPoint
1861 (Heap, currBVP->patches[PARAM_PATCH_POINTS (p, 1)]));
1862 return (CreateBndPOnLine
1863 (Heap, currBVP->patches[PARAM_PATCH_POINTS (p, 0)],
1864 currBVP->patches[PARAM_PATCH_POINTS (p, 1)], lc));
1865 }
1866 else if (abs(pos[1] - PARAM_PATCH_RANGE (p)[1][1]) < SMALL_DIFF)
1867 {
1868 lc = (pos[0] - PARAM_PATCH_RANGE (p)[0][0])
1869 / (PARAM_PATCH_RANGE (p)[1][0] - PARAM_PATCH_RANGE (p)[0][0]);
1870 if (abs(lc) < SMALL_DIFF)
1871 return (CreateBndPOnPoint
1872 (Heap, currBVP->patches[PARAM_PATCH_POINTS (p, 3)]));
1873 else if (abs(lc - 1.) < SMALL_DIFF)
1874 return (CreateBndPOnPoint
1875 (Heap, currBVP->patches[PARAM_PATCH_POINTS (p, 2)]));
1876 return (CreateBndPOnLine
1877 (Heap, currBVP->patches[PARAM_PATCH_POINTS (p, 3)],
1878 currBVP->patches[PARAM_PATCH_POINTS (p, 2)], lc));
1879 }
1880 #endif
1881 #ifdef __TWODIM__
1882 /* check point on point patch */
1883 if (abs(pos[0] - PARAM_PATCH_RANGE (p)[0][0]) < SMALL_DIFF)
1884 return (CreateBndPOnPoint
1885 (Heap, currBVP->patches[PARAM_PATCH_POINTS (p, 0)]));
1886 else if (abs(pos[0] - PARAM_PATCH_RANGE (p)[0][1]) < SMALL_DIFF)
1887 return (CreateBndPOnPoint
1888 (Heap, currBVP->patches[PARAM_PATCH_POINTS (p, 1)]));
1889 #endif
1890
1891 if (PATCH_TYPE (p) == PARAMETRIC_PATCH_TYPE)
1892 {
1893 PRINTDEBUG (dom, 1, (" id %d i %d ns %d %s pos %f\n",
1894 pid, i, currBVP->sideoffset, argv[0], pos[0]));
1895 ps = (BND_PS *) GetFreelistMemory (Heap, sizeof (BND_PS));
1896 if (ps == NULL)
1897 return (NULL);
1898 ps->patch_id = pid;
1899 ps->n = 1;
1900 for (j = 0; j < DIM_OF_BND; j++)
1901 ps->local[0][j] = pos[j];
1902 }
1903 else
1904 return (NULL);
1905
1906 if (!PATCH_IS_FIXED (p))
1907 {
1908 /* store global coordinates */
1909 BND_DATA (ps) = GetFreelistMemory (Heap, DIM * sizeof (DOUBLE));
1910 if (BND_DATA (ps) == NULL)
1911 return (NULL);
1912
1913 if (BndPointGlobal ((BNDP *) ps, (DOUBLE *) BND_DATA (ps)))
1914 return (NULL);
1915 }
1916
1917 return ((BNDP *) ps);
1918 }
1919
1920 /****************************************************************************/
1921 /** \brief Decompose a strip into triangles
1922 *
1923 * \todo Which function does this doc block belong to?
1924 *
1925 * @param n,m - stripe with n+1 nodes on the bottom and m+1 nodes on the top
1926 * @param c0,c1,c2,c3 - corner node ids
1927 * @param s0,s1,s2,s3 - side node ids
1928 *
1929 * This function splits a stripe into triangles and calls AddBoundaryElements().
1930 *
1931 * @return <ul>
1932 * <li> 0 if ok </li>
1933 * <li> 1 if error occured </li>
1934 * </ul>
1935 */
1936 /****************************************************************************/
1937
1938 static INT
AddBoundaryElement(INT n,INT * nodelist,INT left,INT right,INT *** corners,INT * sides)1939 AddBoundaryElement (INT n, INT * nodelist,
1940 INT left, INT right, INT *** corners, INT * sides)
1941 {
1942 if (left > 0)
1943 {
1944 if (corners != NULL)
1945 {
1946 corners[left][sides[left]][0] = nodelist[0];
1947 corners[left][sides[left]][1] = nodelist[1];
1948 corners[left][sides[left]][2] = nodelist[2];
1949 }
1950 sides[left]++;
1951 }
1952 if (right > 0)
1953 {
1954 if (corners != NULL)
1955 {
1956 corners[right][sides[right]][0] = nodelist[0];
1957 corners[right][sides[right]][1] = nodelist[2];
1958 corners[right][sides[right]][2] = nodelist[1];
1959 }
1960 sides[right]++;
1961 }
1962
1963 return (0);
1964 }
1965
1966 static INT
AddBoundaryElements(INT n,INT m,INT c0,INT c1,INT c2,INT c3,INT s0,INT s1,INT s2,INT s3,INT left,INT right,INT *** corners,INT * sides)1967 AddBoundaryElements (INT n, INT m,
1968 INT c0, INT c1, INT c2, INT c3,
1969 INT s0, INT s1, INT s2, INT s3,
1970 INT left, INT right, INT *** corners, INT * sides)
1971 {
1972 INT nodelist[3];
1973
1974 PRINTDEBUG (dom, 1, (" add n %d m %d c0 %d c1 %d s0 %d s1 %d\n",
1975 n, m, c0, c1, s0, s1));
1976 if (m < n)
1977 {
1978 if (n == 1)
1979 {
1980 nodelist[0] = c0;
1981 nodelist[1] = c2;
1982 nodelist[2] = c1;
1983 AddBoundaryElement (3, nodelist, left, right, corners, sides);
1984 }
1985 else
1986 {
1987 nodelist[0] = c0;
1988 nodelist[1] = c2;
1989 nodelist[2] = s0;
1990 AddBoundaryElement (3, nodelist, left, right, corners, sides);
1991 c0 = s0;
1992 if (s0 < s1)
1993 s0++;
1994 else
1995 s0--;
1996 AddBoundaryElements (n - 1, m, c0, c1, c2, c3, s0, s1, s2, s3,
1997 left, right, corners, sides);
1998 }
1999 }
2000 else
2001 {
2002 if (m == 1)
2003 {
2004 nodelist[0] = c0;
2005 nodelist[1] = c2;
2006 nodelist[2] = c1;
2007 AddBoundaryElement (3, nodelist, left, right, corners, sides);
2008 nodelist[0] = c1;
2009 nodelist[1] = c2;
2010 nodelist[2] = c3;
2011 AddBoundaryElement (3, nodelist, left, right, corners, sides);
2012 }
2013 else
2014 {
2015 nodelist[0] = c2;
2016 nodelist[1] = s2;
2017 nodelist[2] = c0;
2018 AddBoundaryElement (3, nodelist, left, right, corners, sides);
2019 c2 = s2;
2020 if (s2 < s3)
2021 s2++;
2022 else
2023 s2--;
2024 AddBoundaryElements (n, m - 1, c0, c1, c2, c3, s0, s1, s2, s3,
2025 left, right, corners, sides);
2026 }
2027 }
2028
2029 return (0);
2030 }
2031
2032 static INT nodeid;
2033
2034 /****************************************************************************/
2035 /** \brief Decompose a patch into stripes
2036 *
2037 * @param h - maximal length of an edge of the boundary triangles
2038 * @param thePatch - poiter to a patch
2039 * @param npc - number of corner nodes of the patch
2040 * @param cornerid - ids of the corner nodes
2041 * @param local - local coordinates of the corners
2042 * @param sideid - ids of the nodes on the edges of the patch
2043 * @param siden - number of nodes on the edges of the patch
2044 *
2045 * This function splits the patch into stripes and calls
2046 * AddBoundaryElements() to decompose the stripes into triangles.
2047 *
2048 * @return <ul>
2049 * <li> 0 if ok </li>
2050 * <li> 1 if error occured </li>
2051 * </ul>
2052 */
2053 /****************************************************************************/
2054
2055
2056 static INT
TriangulatePatch(HEAP * Heap,PATCH * p,BNDP ** bndp,INT * sides,INT *** corners,DOUBLE h,INT npc,INT * cornerid,DOUBLE local[CORNERS_OF_BND_SEG][DIM-1],INT sideid[CORNERS_OF_BND_SEG][2],INT * siden)2057 TriangulatePatch (HEAP * Heap, PATCH * p, BNDP ** bndp,
2058 INT * sides, INT *** corners, DOUBLE h,
2059 INT npc, INT * cornerid,
2060 DOUBLE local[CORNERS_OF_BND_SEG][DIM - 1],
2061 INT sideid[CORNERS_OF_BND_SEG][2], INT * siden)
2062 {
2063 BND_PS *ps;
2064 INT i, k, left, right;
2065 DOUBLE gvect[DIM], gvect1[DIM];
2066 DOUBLE next_local[CORNERS_OF_BND_SEG][DIM - 1];
2067 DOUBLE lambda, dist, step;
2068 INT next_cornerid[CORNERS_OF_BND_SEG];
2069 INT next_sideid[CORNERS_OF_BND_SEG][2];
2070 INT next_siden[CORNERS_OF_BND_SEG];
2071
2072 left = PARAM_PATCH_LEFT (p);
2073 right = PARAM_PATCH_RIGHT (p);
2074
2075
2076 PRINTDEBUG (dom, 1, ("Triang nid %d sn %d %d %d %d\n", nodeid,
2077 siden[0], siden[1], siden[2], siden[3]));
2078
2079 if ((siden[npc - 1] > 1) && (siden[1] > 1))
2080 {
2081 for (k = 2; k < npc; k++)
2082 for (i = 0; i < DIM - 1; i++)
2083 next_local[k][i] = local[k][i];
2084 next_siden[0] = siden[0];
2085 next_siden[1] = siden[1] - 1;
2086 next_siden[2] = siden[2];
2087 next_siden[npc - 1] = siden[npc - 1] - 1;
2088 next_sideid[2][0] = sideid[2][0];
2089 next_sideid[2][1] = sideid[2][1];
2090 for (k = 2; k < npc; k++)
2091 next_cornerid[k] = cornerid[k];
2092 next_cornerid[0] = sideid[npc - 1][1];
2093 next_sideid[npc - 1][0] = sideid[npc - 1][0];
2094 if (sideid[npc - 1][0] < sideid[npc - 1][1])
2095 next_sideid[npc - 1][1] = sideid[npc - 1][1] - 1;
2096 else
2097 next_sideid[npc - 1][1] = sideid[npc - 1][1] + 1;
2098 lambda = (siden[npc - 1] - 1.0) / siden[npc - 1];
2099 V2_LINCOMB (lambda, local[0], (1.0 - lambda), local[npc - 1],
2100 next_local[0]);
2101 next_cornerid[1] = sideid[1][0];
2102 next_sideid[1][1] = sideid[1][1];
2103 if (sideid[1][0] < sideid[1][1])
2104 next_sideid[1][0] = sideid[1][0] + 1;
2105 else
2106 next_sideid[1][0] = sideid[1][0] - 1;
2107 lambda = (siden[1] - 1.0) / siden[1];
2108 V2_LINCOMB (lambda, local[1], (1.0 - lambda), local[2], next_local[1]);
2109 next_siden[0] = siden[0];
2110 if (h > 0)
2111 {
2112 if ((*PARAM_PATCH_BS (p))
2113 (PARAM_PATCH_BSD (p), next_local[0], gvect))
2114 return (-1);
2115 if ((*PARAM_PATCH_BS (p))
2116 (PARAM_PATCH_BSD (p), next_local[1], gvect1))
2117 return (-1);
2118 V_DIM_EUKLIDNORM_OF_DIFF (gvect, gvect1, dist);
2119 if (((INT) (dist / h)) < siden[0])
2120 next_siden[0] -= 1;
2121 else if (((INT) (dist / h)) > siden[0])
2122 next_siden[0] += 1;
2123 }
2124 next_sideid[0][0] = nodeid;
2125 next_sideid[0][1] = nodeid + next_siden[0] - 2;
2126 step = 1.0 / next_siden[0];
2127 if (bndp != NULL)
2128 for (i = 1; i < next_siden[0]; i++)
2129 {
2130 ps = (BND_PS *) GetFreelistMemory (Heap, sizeof (BND_PS));
2131 if (ps == NULL)
2132 return (0);
2133 ps->n = 1;
2134 ps->patch_id = PATCH_ID (p);
2135 lambda = i * step;
2136 V2_LINCOMB (lambda, next_local[1], (1.0 - lambda), next_local[0],
2137 ps->local[0]);
2138 if (!PATCH_IS_FIXED (p))
2139 {
2140 /* store global coordinates */
2141 BND_DATA (ps) =
2142 GetFreelistMemory (Heap, DIM * sizeof (DOUBLE));
2143 if (BND_DATA (ps) == NULL)
2144 return (1);
2145
2146 if (BndPointGlobal ((BNDP *) ps, (DOUBLE *) BND_DATA (ps)))
2147 return (1);
2148 }
2149 bndp[nodeid++] = (BNDP *) ps;
2150
2151
2152 PRINTDEBUG (dom, 1, (" lambda nid %d %f %f\n",
2153 nodeid - 1, ps->local[0][0],
2154 ps->local[0][1]));
2155
2156 }
2157 else
2158 nodeid += next_siden[0] - 1;
2159
2160 AddBoundaryElements (siden[0], next_siden[0],
2161 cornerid[0], cornerid[1],
2162 next_cornerid[0], next_cornerid[1],
2163 sideid[0][0], sideid[0][1],
2164 next_sideid[0][0], next_sideid[0][1],
2165 left, right, corners, sides);
2166
2167 return (TriangulatePatch
2168 (Heap, p, bndp, sides, corners, h, npc, next_cornerid,
2169 next_local, next_sideid, next_siden));
2170 }
2171 else if ((siden[npc - 1] == 1) && (siden[1] == 1))
2172 {
2173 if (npc == 3)
2174 return (AddBoundaryElements (siden[0], 0,
2175 cornerid[0], cornerid[1], cornerid[2], 0,
2176 sideid[0][0], sideid[0][1], 0, 0,
2177 left, right, corners, sides));
2178 else
2179 return (AddBoundaryElements (siden[0], siden[2],
2180 cornerid[0], cornerid[1],
2181 cornerid[3], cornerid[2],
2182 sideid[0][0], sideid[0][1],
2183 sideid[2][1], sideid[2][0],
2184 left, right, corners, sides));
2185 }
2186 else if (((siden[npc - 1] > 1) && (siden[1] == 1)) ||
2187 ((siden[npc - 1] == 1) && (siden[1] > 1)))
2188 {
2189 if ((siden[npc - 2] > 1) && (siden[1] == 1) && (siden[0] == 1))
2190 {
2191 UserWrite ("TriangulatePatch: this case is not implemented\n");
2192 return (1);
2193 }
2194
2195 /* swap sides */
2196 for (k = 0; k < npc; k++)
2197 {
2198 for (i = 0; i < DIM - 1; i++)
2199 next_local[k][i] = local[(k + 1) % npc][i];
2200 next_siden[k] = siden[(k + 1) % npc];
2201 next_sideid[k][0] = sideid[(k + 1) % npc][0];
2202 next_sideid[k][1] = sideid[(k + 1) % npc][1];
2203 next_cornerid[k] = cornerid[(k + 1) % npc];
2204 }
2205 return (TriangulatePatch
2206 (Heap, p, bndp, sides, corners, h, npc, next_cornerid,
2207 next_local, next_sideid, next_siden));
2208 }
2209
2210 return (0);
2211 }
2212
2213 #define INDEX_IN_LIST(from,to,nc) (from < to ? from*nc+to : to*nc+from)
2214
2215 /****************************************************************************/
2216 /****************************************************************************/
2217 /* */
2218 /* functions for BNDS */
2219 /* */
2220 /****************************************************************************/
2221 /****************************************************************************/
2222
2223 static INT
PatchGlobal(const PATCH * p,DOUBLE * lambda,DOUBLE * global)2224 PatchGlobal (const PATCH * p, DOUBLE * lambda, DOUBLE * global)
2225 {
2226 if (PATCH_TYPE (p) == PARAMETRIC_PATCH_TYPE)
2227 return ((*PARAM_PATCH_BS (p))(PARAM_PATCH_BSD (p), lambda, global));
2228 else if (PATCH_TYPE (p) == LINEAR_PATCH_TYPE)
2229 {
2230 #ifdef __TWODIM__
2231 global[0] = (1 - lambda[0]) * LINEAR_PATCH_POS (p, 0)[0]
2232 + lambda[0] * LINEAR_PATCH_POS (p, 1)[0];
2233 global[1] = (1 - lambda[0]) * LINEAR_PATCH_POS (p, 0)[1]
2234 + lambda[0] * LINEAR_PATCH_POS (p, 1)[1];
2235 #endif
2236 #ifdef __THREEDIM__
2237
2238 if (LINEAR_PATCH_N(p) ==3) {
2239 /* Linear interpolation for a triangle boundary segment */
2240 int i;
2241 for (i=0; i<3; i++)
2242 global[i] = (1 - lambda[0] - lambda[1]) * LINEAR_PATCH_POS (p, 0)[i]
2243 + lambda[0] * LINEAR_PATCH_POS (p, 1)[i]
2244 + lambda[1] * LINEAR_PATCH_POS (p, 2)[i];
2245
2246 } else {
2247 /* Bilinear interpolation for a quadrilateral boundary segment */
2248 int i;
2249 for (i=0; i<3; i++)
2250 global[i] = LINEAR_PATCH_POS (p, 0)[i]
2251 + lambda[0]*(LINEAR_PATCH_POS (p, 1)[i] - LINEAR_PATCH_POS (p, 0)[i])
2252 + lambda[1]*(LINEAR_PATCH_POS (p, 3)[i] - LINEAR_PATCH_POS (p, 0)[i])
2253 + lambda[0]*lambda[1]*(LINEAR_PATCH_POS (p, 0)[i]+LINEAR_PATCH_POS (p, 2)[i]-LINEAR_PATCH_POS (p, 1)[i]-LINEAR_PATCH_POS (p, 3)[i]);
2254 }
2255 #endif
2256 return (0);
2257 }
2258
2259 return (1);
2260 }
2261
2262 static INT
FreeBNDS_Global(BND_PS * ps,DOUBLE * local,DOUBLE * global)2263 FreeBNDS_Global (BND_PS * ps, DOUBLE * local, DOUBLE * global)
2264 {
2265 BND_PS **ppt;
2266 PATCH *p;
2267 DOUBLE *pos[4];
2268 INT i;
2269
2270 p = currBVP->patches[ps->patch_id];
2271 if (p == NULL)
2272 return (1);
2273
2274 /* get corner coordinates */
2275 ppt = (BND_PS **) BND_DATA (ps);
2276 ASSERT (BND_N (ps) <= 4);
2277 for (i = 0; i < BND_N (ps); i++)
2278 {
2279 ASSERT (ppt != NULL);
2280 pos[i] = (DOUBLE *) BND_DATA (*(ppt++));
2281 }
2282
2283 /* claculate global coordinates */
2284 #ifdef __TWODIM__
2285 for (i = 0; i < DIM; i++)
2286 global[i] = (1.0 - local[0]) * pos[0][i] + local[0] * pos[1][i];
2287 #endif
2288 #ifdef __THREEDIM__
2289 switch (ps->n)
2290 {
2291 case 3 :
2292 for (i = 0; i < DIM; i++)
2293 global[i] = (1.0 - local[0] - local[1]) * pos[0][i]
2294 + local[0] * pos[1][i] + local[1] * pos[2][i];
2295 break;
2296 case 4 :
2297 for (i = 0; i < DIM; i++)
2298 global[i] = (1.0 - local[0]) * (1.0 - local[1]) * pos[0][i]
2299 + local[0] * (1.0 - local[1]) * pos[1][i]
2300 + local[0] * local[1] * pos[2][i]
2301 + (1.0 - local[0]) * local[1] * pos[3][i];
2302 break;
2303 }
2304 #endif
2305
2306 return (0);
2307 }
2308
2309 static INT
local2lambda(BND_PS * ps,DOUBLE local[],DOUBLE lambda[])2310 local2lambda (BND_PS * ps, DOUBLE local[], DOUBLE lambda[])
2311 {
2312 PATCH *p;
2313
2314 p = currBVP->patches[ps->patch_id];
2315
2316 if ((PATCH_TYPE (p) == PARAMETRIC_PATCH_TYPE)
2317 || (PATCH_TYPE (p) == LINEAR_PATCH_TYPE))
2318 #ifdef __TWODIM__
2319 lambda[0] =
2320 (1.0 - local[0]) * ps->local[0][0] + local[0] * ps->local[1][0];
2321 #endif
2322 #ifdef __THREEDIM__
2323 switch (ps->n)
2324 {
2325 case 3 :
2326 lambda[0] = (1.0 - local[0] - local[1]) * ps->local[0][0]
2327 + local[0] * ps->local[1][0] + local[1] * ps->local[2][0];
2328 lambda[1] = (1.0 - local[0] - local[1]) * ps->local[0][1]
2329 + local[0] * ps->local[1][1] + local[1] * ps->local[2][1];
2330 break;
2331 case 4 :
2332 lambda[0] = (1.0 - local[0]) * (1.0 - local[1]) * ps->local[0][0]
2333 + local[0] * (1.0 - local[1]) * ps->local[1][0]
2334 + local[0] * local[1] * ps->local[2][0]
2335 + (1.0 - local[0]) * local[1] * ps->local[3][0];
2336 lambda[1] = (1.0 - local[0]) * (1.0 - local[1]) * ps->local[0][1]
2337 + local[0] * (1.0 - local[1]) * ps->local[1][1]
2338 + local[0] * local[1] * ps->local[2][1]
2339 + (1.0 - local[0]) * local[1] * ps->local[3][1];
2340 break;
2341 }
2342 #endif
2343 else
2344 return (1);
2345
2346 return (0);
2347 }
2348
2349 /* domain interface function: for description see domain.h */
2350 INT NS_DIM_PREFIX
BNDS_Global(BNDS * aBndS,DOUBLE * local,DOUBLE * global)2351 BNDS_Global (BNDS * aBndS, DOUBLE * local, DOUBLE * global)
2352 {
2353 BND_PS *ps;
2354 PATCH *p;
2355 DOUBLE lambda[DIM_OF_BND];
2356
2357 ps = (BND_PS *) aBndS;
2358 p = currBVP->patches[ps->patch_id];
2359 if (p == NULL)
2360 return (1);
2361
2362 PRINTDEBUG (dom, 1, (" Bnds global loc %f pid %d\n",
2363 local[0], PATCH_ID (p)));
2364
2365 if (PATCH_IS_FREE (p))
2366 return (FreeBNDS_Global (ps, local, global));
2367
2368 if (local2lambda (ps, local, lambda))
2369 return (1);
2370
2371 return (PatchGlobal (p, lambda, global));
2372 }
2373
2374 static INT
SideIsCooriented(BND_PS * ps)2375 SideIsCooriented (BND_PS * ps)
2376 {
2377 # ifdef __TWODIM__
2378 if (BND_LOCAL (ps, 1)[0] > BND_LOCAL (ps, 0)[0])
2379 return (YES);
2380 else
2381 return (NO);
2382 # endif
2383
2384 # ifdef __THREEDIM__
2385 DOUBLE vp, x0[2], x1[2];
2386
2387 ASSERT (BND_N (ps) >= 3);
2388
2389 /* check whether an (arbitrary) angle of the side is > 180 degree */
2390 V2_SUBTRACT (BND_LOCAL (ps, 1), BND_LOCAL (ps, 0), x0);
2391 V2_SUBTRACT (BND_LOCAL (ps, 2), BND_LOCAL (ps, 0), x1);
2392 V2_VECTOR_PRODUCT (x1, x0, vp);
2393
2394 ASSERT (fabs (vp) > SMALL_C);
2395
2396 if (vp > SMALL_C)
2397 return (YES);
2398 else
2399 return (NO);
2400 # endif
2401 }
2402
2403 /* domain interface function: for description see domain.h */
2404 INT NS_DIM_PREFIX
BNDS_BndSDesc(BNDS * theBndS,INT * id,INT * nbid,INT * part)2405 BNDS_BndSDesc (BNDS * theBndS, INT * id, INT * nbid, INT * part)
2406 {
2407 BND_PS *ps;
2408 PATCH *p;
2409 INT left, right;
2410
2411 ps = (BND_PS *) theBndS;
2412 p = currBVP->patches[ps->patch_id];
2413
2414 /* fill part from segment */
2415 if (STD_BVP_NDOMPART (currBVP) > 1)
2416 {
2417 *part = DPI_SG2P (DOMAIN_PARTINFO (STD_BVP_DOMAIN (currBVP)),
2418 PATCH_ID (p) - STD_BVP_SIDEOFFSET (currBVP));
2419 /* this expression yields the segment id */
2420 }
2421 else
2422 /* default is 0 */
2423 *part = 0;
2424
2425 if (PATCH_TYPE (p) == PARAMETRIC_PATCH_TYPE)
2426 {
2427 left = PARAM_PATCH_LEFT (p);
2428 right = PARAM_PATCH_RIGHT (p);
2429 }
2430 else if (PATCH_TYPE (p) == LINEAR_PATCH_TYPE)
2431 {
2432 left = LINEAR_PATCH_LEFT (p);
2433 right = LINEAR_PATCH_RIGHT (p);
2434 }
2435 else
2436 return (1);
2437
2438 /* check orientation */
2439 if (SideIsCooriented (ps))
2440 {
2441 /* patch and side are co-oriented */
2442 *id = left;
2443 *nbid = right;
2444 }
2445 else
2446 {
2447 /* patch and side are anti-oriented */
2448 *id = right;
2449 *nbid = left;
2450 }
2451
2452 return (0);
2453 }
2454
2455 /* domain interface function: for description see domain.h */
2456 BNDP *NS_DIM_PREFIX
BNDS_CreateBndP(HEAP * Heap,BNDS * aBndS,DOUBLE * local)2457 BNDS_CreateBndP (HEAP * Heap, BNDS * aBndS, DOUBLE * local)
2458 {
2459 BND_PS *ps, *pp;
2460 PATCH *p;
2461
2462 if (aBndS == NULL)
2463 return (NULL);
2464
2465 ps = (BND_PS *) aBndS;
2466 p = currBVP->patches[ps->patch_id];
2467
2468 pp = (BND_PS *) GetFreelistMemory (Heap, sizeof (BND_PS));
2469 if (pp == NULL)
2470 return (NULL);
2471
2472 pp->patch_id = ps->patch_id;
2473 pp->n = 1;
2474
2475 if (local2lambda (ps, local, pp->local[0]))
2476 return (NULL);
2477
2478 if (!PATCH_IS_FIXED (p))
2479 {
2480 /* store global coordinates */
2481 BND_DATA (pp) = GetFreelistMemory (Heap, DIM * sizeof (DOUBLE));
2482 if (BND_DATA (pp) == NULL)
2483 return (NULL);
2484
2485 if (FreeBNDS_Global (ps, pp->local[0], (DOUBLE *) BND_DATA (pp)))
2486 return (NULL);
2487 }
2488
2489 PRINTDEBUG (dom, 1, (" BNDP s %d\n", pp->patch_id));
2490
2491 return ((BNDP *) pp);
2492 }
2493
2494 /****************************************************************************/
2495 /****************************************************************************/
2496 /* */
2497 /* functions for BNDP */
2498 /* */
2499 /****************************************************************************/
2500 /****************************************************************************/
2501
2502 static INT
BndPointGlobal(const BNDP * aBndP,DOUBLE * global)2503 BndPointGlobal (const BNDP * aBndP, DOUBLE * global)
2504 {
2505 using std::abs;
2506 BND_PS *ps;
2507 PATCH *p, *s;
2508 INT j, k;
2509 DOUBLE pglobal[DIM];
2510
2511 ps = (BND_PS *) aBndP;
2512 p = currBVP->patches[ps->patch_id];
2513
2514 PRINTDEBUG (dom, 1, (" bndp pid %d %d %d\n", ps->patch_id,
2515 PATCH_ID (p), PATCH_TYPE (p)));
2516
2517 switch (PATCH_TYPE (p))
2518 {
2519 case PARAMETRIC_PATCH_TYPE :
2520 case LINEAR_PATCH_TYPE :
2521 return (PatchGlobal (p, ps->local[0], global));
2522 case POINT_PATCH_TYPE :
2523
2524 s = currBVP->patches[POINT_PATCH_PID (p, 0)];
2525
2526 PRINTDEBUG (dom, 1, (" bndp n %d %d loc %f %f gl \n",
2527 POINT_PATCH_N (p),
2528 POINT_PATCH_PID (p, 0),
2529 ps->local[0][0], ps->local[0][1]));
2530
2531 PatchGlobal(s, ps->local[0], global);
2532
2533 for (j = 1; j < POINT_PATCH_N (p); j++)
2534 {
2535 s = currBVP->patches[POINT_PATCH_PID (p, j)];
2536
2537 if (PatchGlobal(s, ps->local[j], pglobal))
2538 REP_ERR_RETURN (1);
2539
2540 PRINTDEBUG (dom,1, (" bndp j %d %d loc %f %f gl %f %f %f\n", j,
2541 POINT_PATCH_PID (p, j),
2542 ps->local[j][0],
2543 ps->local[j][1],
2544 pglobal[0], pglobal[1], pglobal[2]));
2545 for (k = 0; k < DIM; k++)
2546 if (abs(pglobal[k] - global[k]) > SMALL_DIFF)
2547 REP_ERR_RETURN (1);
2548 }
2549
2550 return (0);
2551 #ifdef __THREEDIM__
2552 case LINE_PATCH_TYPE :
2553 s = currBVP->patches[LINE_PATCH_PID (p, 0)];
2554
2555 if (PatchGlobal(s, ps->local[0], global))
2556 REP_ERR_RETURN (1);
2557
2558 PRINTDEBUG (dom, 1, (" bndp n %d %d loc %f %f gl %f %f %f\n",
2559 POINT_PATCH_N (p),
2560 LINE_PATCH_PID (p, 0),
2561 ps->local[0][0],
2562 ps->local[0][1], global[0], global[1], global[2]));
2563 for (j = 1; j < LINE_PATCH_N (p); j++)
2564 {
2565 s = currBVP->patches[LINE_PATCH_PID (p, j)];
2566
2567 if (PatchGlobal(s, ps->local[j], pglobal))
2568 REP_ERR_RETURN (1);
2569
2570 PRINTDEBUG (dom, 1, (" bndp j %d %d loc %f %f gl %f %f %f\n", j,
2571 LINE_PATCH_PID (p, j),
2572 ps->local[j][0],
2573 ps->local[j][1],
2574 pglobal[0], pglobal[1], pglobal[2]));
2575 for (k = 0; k < DIM; k++)
2576 if (abs(pglobal[k] - global[k]) > SMALL_DIFF)
2577 REP_ERR_RETURN (1);
2578 }
2579 return (0);
2580 #endif
2581 }
2582
2583 REP_ERR_RETURN (1);
2584 }
2585
2586 /* domain interface function: for description see domain.h */
2587 INT NS_DIM_PREFIX
BNDP_Global(BNDP * aBndP,DOUBLE * global)2588 BNDP_Global (BNDP * aBndP, DOUBLE * global)
2589 {
2590 BND_PS *ps = (BND_PS *) aBndP;
2591 DOUBLE *pos;
2592 INT i;
2593
2594 if (!PATCH_IS_FIXED (currBVP->patches[ps->patch_id]))
2595 {
2596 /* copy globals */
2597 pos = (DOUBLE *) BND_DATA (ps);
2598 for (i = 0; i < DIM; i++)
2599 global[i] = pos[i];
2600 return (0);
2601 }
2602
2603 return (BndPointGlobal (aBndP, global));
2604 }
2605
2606 /* domain interface function: for description see domain.h */
2607 INT NS_DIM_PREFIX
BNDP_BndPDesc(BNDP * theBndP,INT * move,INT * part)2608 BNDP_BndPDesc (BNDP * theBndP, INT * move, INT * part)
2609 {
2610 BND_PS *ps;
2611 PATCH *p;
2612
2613 ps = (BND_PS *) theBndP;
2614 p = STD_BVP_PATCH (currBVP, ps->patch_id);
2615
2616 /* default part is 0 */
2617 *part = 0;
2618
2619 switch (PATCH_TYPE (p))
2620 {
2621 case PARAMETRIC_PATCH_TYPE :
2622 case LINEAR_PATCH_TYPE :
2623 if (STD_BVP_NDOMPART (currBVP) > 1)
2624 *part = DPI_SG2P (DOMAIN_PARTINFO (STD_BVP_DOMAIN (currBVP)), PATCH_ID (p) - STD_BVP_SIDEOFFSET (currBVP)); /* <-- this expression yields the segment id */
2625 *move = PATCH_IS_FREE (p) ? DIM : DIM_OF_BND;
2626 return (0);
2627
2628 case POINT_PATCH_TYPE :
2629 if (STD_BVP_NDOMPART (currBVP) > 1)
2630 *part = DPI_PT2P (DOMAIN_PARTINFO (STD_BVP_DOMAIN (currBVP)), PATCH_ID (p)); /* <-- this expression yields the corner id */
2631 *move = PATCH_IS_FREE (p) ? DIM : 0;
2632 return (0);
2633
2634 # ifdef __THREEDIM__
2635 case LINE_PATCH_TYPE :
2636 if (STD_BVP_NDOMPART (currBVP) > 1)
2637 *part = DPI_LN2P (DOMAIN_PARTINFO (STD_BVP_DOMAIN (currBVP)),
2638 LINE_PATCH_C0 (p), LINE_PATCH_C1 (p));
2639 *move = PATCH_IS_FREE (p) ? DIM : 1;
2640 return (0);
2641 # endif
2642 }
2643
2644 return (1);
2645 }
2646
2647 /* domain interface function: for description see domain.h */
2648 INT NS_DIM_PREFIX
BNDP_BndEDesc(BNDP * aBndP0,BNDP * aBndP1,INT * part)2649 BNDP_BndEDesc (BNDP * aBndP0, BNDP * aBndP1, INT * part)
2650 {
2651 BND_PS *bp0, *bp1;
2652 PATCH *p, *p0, *p1;
2653 INT pid, cnt;
2654
2655 bp0 = (BND_PS *) aBndP0;
2656 p0 = currBVP->patches[bp0->patch_id];
2657 bp1 = (BND_PS *) aBndP1;
2658 p1 = currBVP->patches[bp1->patch_id];
2659
2660 if ((bp0 == NULL) || (bp1 == NULL))
2661 REP_ERR_RETURN (1);
2662
2663 /* default part is 0 */
2664 *part = 0;
2665
2666 if (STD_BVP_NDOMPART (currBVP) == 1)
2667 return (0);
2668
2669 /* find common patch(es) of boundary points */
2670 cnt = GetNumberOfCommonPatches (p0, p1, &pid);
2671 if (cnt == 0)
2672 return (1);
2673
2674 #ifdef __THREEDIM__
2675 if (cnt > 1)
2676 {
2677 pid = GetCommonLinePatchId (p0, p1);
2678 ASSERT ((pid < currBVP->ncorners) || (pid >= currBVP->sideoffset));
2679 p = STD_BVP_PATCH (currBVP, pid);
2680
2681 *part = DPI_LN2P (DOMAIN_PARTINFO (STD_BVP_DOMAIN (currBVP)),
2682 LINE_PATCH_C0 (p), LINE_PATCH_C1 (p));
2683 return (0);
2684 }
2685 #endif
2686
2687 p = STD_BVP_PATCH (currBVP, pid);
2688 switch (PATCH_TYPE (p))
2689 {
2690 case PARAMETRIC_PATCH_TYPE :
2691 case LINEAR_PATCH_TYPE :
2692 *part = DPI_SG2P (DOMAIN_PARTINFO (STD_BVP_DOMAIN (currBVP)), PATCH_ID (p) - STD_BVP_SIDEOFFSET (currBVP)); /* <-- this expression yields the segment id */
2693 break;
2694 default :
2695 REP_ERR_RETURN (1);
2696 }
2697
2698 return (0);
2699 }
2700
2701 /* domain interface function: for description see domain.h */
2702 BNDS *NS_DIM_PREFIX
BNDP_CreateBndS(HEAP * Heap,BNDP ** aBndP,INT n)2703 BNDP_CreateBndS (HEAP * Heap, BNDP ** aBndP, INT n)
2704 {
2705 BND_PS *bp[4], *bs, **pps;
2706 PATCH *p[4];
2707 DOUBLE *lambda[4];
2708 INT i, j, l, pid;
2709 # ifdef __THREEDIM__
2710 INT k;
2711 # endif
2712
2713 PRINTDEBUG (dom, 1, ("Create BNDS:\n"));
2714 for (i = 0; i < n; i++)
2715 {
2716 bp[i] = (BND_PS *) aBndP[i];
2717 p[i] = currBVP->patches[bp[i]->patch_id];
2718
2719 PRINTDEBUG (dom, 1, (" bp %d p %d n %d\n",
2720 bp[i]->patch_id, PATCH_ID (p[i]), n));
2721 for (l = 0; l < GetNumberOfPatches (p[i]); l++)
2722 PRINTDEBUG (dom, 1, (" bp pid %d\n", GetPatchId (p[i], l)));
2723 }
2724
2725 pid = -1;
2726 #ifdef __TWODIM__
2727 if (n != 2)
2728 return (NULL);
2729 /* find common patch of boundary points */
2730 for (i = 0; i < GetNumberOfPatches (p[0]); i++)
2731 for (j = 0; j < GetNumberOfPatches (p[1]); j++)
2732 {
2733 PRINTDEBUG (dom, 1, (" pid i, j %d %d %d %d\n", i, j,
2734 GetPatchId (p[0], i), GetPatchId (p[1], j)));
2735 if (GetPatchId (p[0], i) == GetPatchId (p[1], j))
2736 {
2737 pid = GetPatchId (p[0], i);
2738 lambda[0] = bp[0]->local[i];
2739 lambda[1] = bp[1]->local[j];
2740 break;
2741 PRINTDEBUG (dom, 1, (" pid %d \n", pid));
2742 }
2743 }
2744 #endif
2745 #ifdef __THREEDIM__
2746 switch (n)
2747 {
2748 case 3 :
2749 for (i = 0; i < GetNumberOfPatches (p[0]); i++)
2750 for (j = 0; j < GetNumberOfPatches (p[1]); j++)
2751 if (GetPatchId (p[0], i) == GetPatchId (p[1], j))
2752 for (k = 0; k < GetNumberOfPatches (p[2]); k++)
2753 if (GetPatchId (p[0], i) == GetPatchId (p[2], k))
2754 {
2755 pid = GetPatchId (p[0], i);
2756 lambda[0] = bp[0]->local[i];
2757 lambda[1] = bp[1]->local[j];
2758 lambda[2] = bp[2]->local[k];
2759 break;
2760 }
2761 break;
2762 case 4 :
2763 for (i = 0; i < GetNumberOfPatches (p[0]); i++)
2764 for (j = 0; j < GetNumberOfPatches (p[1]); j++)
2765 if (GetPatchId (p[0], i) == GetPatchId (p[1], j))
2766 for (k = 0; k < GetNumberOfPatches (p[2]); k++)
2767 if (GetPatchId (p[0], i) == GetPatchId (p[2], k))
2768 for (l = 0; l < GetNumberOfPatches (p[3]); l++)
2769 if (GetPatchId (p[0], i) == GetPatchId (p[3], l))
2770 {
2771 pid = GetPatchId (p[0], i);
2772 lambda[0] = bp[0]->local[i];
2773 lambda[1] = bp[1]->local[j];
2774 lambda[2] = bp[2]->local[k];
2775 lambda[3] = bp[3]->local[l];
2776 break;
2777 }
2778 break;
2779 }
2780 #endif
2781
2782 if (pid == -1)
2783 return (NULL);
2784
2785 bs = (BND_PS *) GetFreelistMemory (Heap, (n - 1) * sizeof (COORD_BND_VECTOR)
2786 + sizeof (BND_PS));
2787 if (bs == NULL)
2788 return (NULL);
2789 bs->n = n;
2790 bs->patch_id = pid;
2791 for (i = 0; i < n; i++)
2792 for (j = 0; j < DIM_OF_BND; j++)
2793 bs->local[i][j] = lambda[i][j];
2794
2795 for (i = 0; i < n; i++)
2796 for (j = 0; j < DIM_OF_BND; j++)
2797 PRINTDEBUG (dom, 1, (" bnds i, j %d %d %f\n", i, j, lambda[i][j]));
2798
2799 PRINTDEBUG (dom, 1, (" Create BNDS %x %d\n", bs, pid));
2800
2801 if (!PATCH_IS_FIXED (currBVP->patches[pid]))
2802 {
2803 /* store corner patch pointers */
2804 BND_DATA (bs) = GetFreelistMemory (Heap, n * sizeof (BND_PS *));
2805 if (BND_DATA (bs) == NULL)
2806 return (NULL);
2807
2808 pps = (BND_PS **) BND_DATA (bs);
2809 for (i = 0; i < n; i++)
2810 *(pps++) = bp[i];
2811 }
2812
2813 return ((BNDS *) bs);
2814 }
2815
2816 /* domain interface function: for description see domain.h */
2817 BNDP *NS_DIM_PREFIX
BNDP_CreateBndP(HEAP * Heap,BNDP * aBndP0,BNDP * aBndP1,DOUBLE lcoord)2818 BNDP_CreateBndP (HEAP * Heap, BNDP * aBndP0, BNDP * aBndP1, DOUBLE lcoord)
2819 {
2820 BND_PS *bp0, *bp1, *bp;
2821 PATCH *p0, *p1;
2822 INT i, j, k, l, cnt;
2823
2824 bp0 = (BND_PS *) aBndP0;
2825 bp1 = (BND_PS *) aBndP1;
2826
2827 if ((bp0 == NULL) || (bp1 == NULL))
2828 return (NULL);
2829
2830 p0 = currBVP->patches[bp0->patch_id];
2831 p1 = currBVP->patches[bp1->patch_id];
2832
2833 PRINTDEBUG (dom, 1, (" bp0 %d pid %d\n", bp0->patch_id, PATCH_ID (p0)));
2834 for (l = 0; l < GetNumberOfPatches (p0); l++)
2835 PRINTDEBUG (dom, 1, (" bp pid %d\n", GetPatchId (p0, l)));
2836 PRINTDEBUG (dom, 1, (" bp1 %d pid %d\n", bp1->patch_id, PATCH_ID (p1)));
2837 for (l = 0; l < GetNumberOfPatches (p1); l++)
2838 PRINTDEBUG (dom, 1, (" bp pid %d\n", GetPatchId (p1, l)));
2839
2840 cnt = GetNumberOfCommonPatches (p0, p1, &i);
2841 if (cnt == 0)
2842 return (NULL);
2843
2844 bp =
2845 (BND_PS *) GetFreelistMemory (Heap,
2846 (cnt - 1) * sizeof (COORD_BND_VECTOR) +
2847 sizeof (BND_PS));
2848 if (bp == NULL)
2849 return (NULL);
2850 bp->n = cnt;
2851
2852 #ifdef __THREEDIM__
2853 if (cnt > 1)
2854 {
2855 PATCH *p;
2856 k = GetCommonLinePatchId (p0, p1);
2857 if ((k < currBVP->ncorners) || (k >= currBVP->sideoffset))
2858 return (NULL);
2859 p = currBVP->patches[k];
2860 bp->patch_id = k;
2861
2862 PRINTDEBUG (dom, 1, (" Create BNDP line %d cnt %d\n", k, cnt));
2863 PRINTDEBUG (dom, 1, (" bp0 %d pid %d\n",
2864 bp0->patch_id, PATCH_ID (p0)));
2865 for (l = 0; l < GetNumberOfPatches (p0); l++)
2866 PRINTDEBUG (dom, 1, (" bp pid %d\n", GetPatchId (p0, l)));
2867 PRINTDEBUG (dom, 1, (" bp1 %d pid %d\n",
2868 bp1->patch_id, PATCH_ID (p1)));
2869 for (l = 0; l < GetNumberOfPatches (p1); l++)
2870 PRINTDEBUG (dom, 1, (" bp pid %d\n", GetPatchId (p1, l)));
2871
2872 for (l = 0; l < LINE_PATCH_N (p); l++)
2873 {
2874 for (i = 0; i < GetNumberOfPatches (p0); i++)
2875 if (GetPatchId (p0, i) == LINE_PATCH_PID (p, l))
2876 for (j = 0; j < GetNumberOfPatches (p1); j++)
2877 if (GetPatchId (p1, j) == LINE_PATCH_PID (p, l))
2878 for (k = 0; k < DIM_OF_BND; k++)
2879 bp->local[l][k] = (1.0 - lcoord) * bp0->local[i][k]
2880 + lcoord * bp1->local[j][k];
2881 }
2882
2883 for (l = 0; l < LINE_PATCH_N (p); l++)
2884 PRINTDEBUG (dom, 1, (" Create BNDP line %d l %d %f %f\n",
2885 LINE_PATCH_PID (p, l), l,
2886 bp->local[l][0], bp->local[l][1]));
2887
2888 if (!PATCH_IS_FIXED (p))
2889 {
2890 /* store global coordinates */
2891 BND_DATA (bp) = GetFreelistMemory (Heap, DIM * sizeof (DOUBLE));
2892 if (BND_DATA (bp) == NULL)
2893 return (NULL);
2894
2895 if (BNDP_Global ((BNDP *) bp, (DOUBLE *) BND_DATA (bp)))
2896 return (NULL);
2897 }
2898 return ((BNDP *) bp);
2899 }
2900 #endif
2901
2902 for (i = 0; i < GetNumberOfPatches (p0); i++)
2903 for (j = 0; j < GetNumberOfPatches (p1); j++)
2904 if (GetPatchId (p0, i) == GetPatchId (p1, j))
2905 {
2906 bp->patch_id = GetPatchId (p0, i);
2907 for (k = 0; k < DIM_OF_BND; k++)
2908 bp->local[0][k] = (1.0 - lcoord) * bp0->local[i][k]
2909 + lcoord * bp1->local[j][k];
2910 PRINTDEBUG (dom, 1,
2911 (" Create BNDP param %d \n", GetPatchId (p0, i)));
2912
2913 break;
2914 }
2915
2916 if (!PATCH_IS_FIXED (currBVP->patches[bp->patch_id]))
2917 {
2918 DOUBLE *x, *a, *b;
2919
2920 /* store global coordinates */
2921 BND_DATA (bp) = GetFreelistMemory (Heap, DIM * sizeof (DOUBLE));
2922 if (BND_DATA (bp) == NULL)
2923 return (NULL);
2924
2925 x = (DOUBLE *) BND_DATA (bp);
2926 a = (DOUBLE *) BND_DATA (bp0);
2927 b = (DOUBLE *) BND_DATA (bp1);
2928 ASSERT (a != NULL);
2929 ASSERT (b != NULL);
2930 for (i = 0; i < DIM; i++)
2931 x[i] = a[i] * (1.0 - lcoord) + b[i] * lcoord;
2932 }
2933
2934 return ((BNDP *) bp);
2935 }
2936
2937 /* domain interface function: for description see domain.h */
2938 INT NS_DIM_PREFIX
BNDP_Move(BNDP * aBndP,const DOUBLE global[])2939 BNDP_Move (BNDP * aBndP, const DOUBLE global[])
2940 {
2941 BND_PS *ps;
2942 DOUBLE *pos;
2943 INT i;
2944
2945 #ifdef ModelP
2946 /* TODO: parallel version */
2947 PrintErrorMessage ('E', "BNDP_Move", "parallel not implemented");
2948 ASSERT (false);
2949 #endif
2950
2951 ps = (BND_PS *) aBndP;
2952 if (PATCH_IS_FREE (currBVP->patches[ps->patch_id]))
2953 {
2954 /* copy globals */
2955 pos = (DOUBLE *) BND_DATA (ps);
2956 for (i = 0; i < DIM; i++)
2957 pos[i] = global[i];
2958 return (0);
2959 }
2960
2961 return (1);
2962 }
2963
2964 /* domain interface function: for description see domain.h */
2965 INT NS_DIM_PREFIX
BNDP_SaveInsertedBndP(BNDP * theBndP,char * data,INT max_data_size)2966 BNDP_SaveInsertedBndP (BNDP * theBndP, char *data, INT max_data_size)
2967 {
2968 BND_PS *bp;
2969 PATCH *p;
2970 INT pid;
2971
2972 bp = (BND_PS *) theBndP;
2973
2974 if (bp == NULL)
2975 return (1);
2976
2977 pid = bp->patch_id;
2978 p = currBVP->patches[pid];
2979
2980 switch (PATCH_TYPE (p))
2981 {
2982 case PARAMETRIC_PATCH_TYPE :
2983 case LINEAR_PATCH_TYPE :
2984 pid -= currBVP->sideoffset;
2985 break;
2986 case POINT_PATCH_TYPE :
2987 pid = POINT_PATCH_PID (p, 0) - currBVP->sideoffset;
2988 break;
2989 #ifdef __THREEDIM__
2990 case LINE_PATCH_TYPE :
2991 pid = LINE_PATCH_PID (p, 0) - currBVP->sideoffset;
2992 break;
2993 #endif
2994 }
2995
2996 PRINTDEBUG (dom, 1, (" Insert pid %d %d\n", bp->patch_id, pid));
2997
2998 #ifdef __TWODIM__
2999 if (sprintf (data, "bn %d %f", pid, (float) bp->local[0][0]) >
3000 max_data_size)
3001 return (1);
3002 #endif
3003
3004 #ifdef __THREEDIM__
3005 if (sprintf (data, "bn %d %f %f", (int) pid,
3006 (float) bp->local[0][0],
3007 (float) bp->local[0][1]) > max_data_size)
3008 return (1);
3009 #endif
3010
3011 return (0);
3012 }
3013
3014 /* domain interface function: for description see domain.h */
3015 INT NS_DIM_PREFIX
BNDP_SurfaceId(BNDP * aBndP,INT * n,INT i)3016 BNDP_SurfaceId (BNDP * aBndP, INT * n, INT i)
3017 {
3018 BND_PS *ps;
3019
3020 if (i < 0)
3021 return (1);
3022
3023 ps = (BND_PS *) aBndP;
3024 if (ps == NULL)
3025 return (-1);
3026
3027 return ps->patch_id;
3028 }
3029
3030 /* domain interface function: for description see domain.h */
3031 INT NS_DIM_PREFIX
BNDP_Dispose(HEAP * Heap,BNDP * theBndP)3032 BNDP_Dispose (HEAP * Heap, BNDP * theBndP)
3033 {
3034 BND_PS *ps;
3035
3036 if (theBndP == NULL)
3037 return (0);
3038
3039 ps = (BND_PS *) theBndP;
3040 if (!PATCH_IS_FIXED (currBVP->patches[ps->patch_id]))
3041 DisposeMem(Heap, BND_DATA (ps));
3042 DisposeMem(Heap, ps);
3043 return 0;
3044 }
3045
3046 /* domain interface function: for description see domain.h */
3047 INT NS_DIM_PREFIX
BNDS_Dispose(HEAP * Heap,BNDS * theBndS)3048 BNDS_Dispose (HEAP * Heap, BNDS * theBndS)
3049 {
3050 BND_PS *ps;
3051
3052 if (theBndS == NULL)
3053 return (0);
3054
3055 ps = (BND_PS *) theBndS;
3056 if (!PATCH_IS_FIXED (currBVP->patches[ps->patch_id]))
3057 DisposeMem(Heap, BND_DATA (ps));
3058 DisposeMem(Heap, ps);
3059 return 0;
3060 }
3061
3062 /* domain interface function: for description see domain.h */
3063 INT NS_DIM_PREFIX
BNDP_SaveBndP(BNDP * BndP)3064 BNDP_SaveBndP (BNDP * BndP)
3065 {
3066 BND_PS *bp;
3067 INT i, j;
3068 int iList[2];
3069 double dList[DIM];
3070
3071 iList[0] = BND_PATCH_ID (BndP);
3072 iList[1] = BND_N (BndP);
3073 if (Bio_Write_mint (2, iList))
3074 return (1);
3075
3076 bp = (BND_PS *) BndP;
3077 for (i = 0; i < BND_N (BndP); i++)
3078 {
3079 for (j = 0; j < DIM - 1; j++)
3080 dList[j] = bp->local[i][j];
3081 if (Bio_Write_mdouble (DIM - 1, dList))
3082 return (1);
3083 }
3084 if (!PATCH_IS_FIXED (currBVP->patches[BND_PATCH_ID (BndP)]))
3085 {
3086 DOUBLE *pos = (DOUBLE *) BND_DATA (bp);
3087
3088 /* save free boundary point coordinates */
3089 for (j = 0; j < DIM; j++)
3090 dList[j] = pos[j];
3091 if (Bio_Write_mdouble (DIM, dList))
3092 return (1);
3093 }
3094
3095 return (0);
3096 }
3097
3098
3099 INT NS_DIM_PREFIX
BNDP_SaveBndP_Ext(BNDP * BndP)3100 BNDP_SaveBndP_Ext (BNDP * BndP)
3101 {
3102 BND_PS *bp;
3103 INT i, j;
3104 int iList[2];
3105 double dList[DIM];
3106
3107 /* TODO: save free boundary points */
3108 iList[0] = BND_PATCH_ID (BndP);
3109 iList[1] = BND_N (BndP);
3110 if (Bio_Write_mint (2, iList))
3111 return (1);
3112
3113 bp = (BND_PS *) BndP;
3114 for (i = 0; i < BND_N (BndP); i++)
3115 {
3116 for (j = 0; j < DIM - 1; j++)
3117 dList[j] = bp->local[i][j];
3118 if (Bio_Write_mdouble (DIM - 1, dList))
3119 return (1);
3120 }
3121 if (!PATCH_IS_FIXED (currBVP->patches[BND_PATCH_ID (BndP)]))
3122 {
3123 DOUBLE *pos = (DOUBLE *) BND_DATA (bp);
3124
3125 /* save free boundary point coordinates */
3126 for (j = 0; j < DIM; j++)
3127 dList[j] = pos[j];
3128 if (Bio_Write_mdouble (DIM, dList))
3129 return (1);
3130 }
3131
3132 return (0);
3133 }
3134
3135 /* domain interface function: for description see domain.h */
3136 BNDP *NS_DIM_PREFIX
BNDP_LoadBndP(BVP * theBVP,HEAP * Heap)3137 BNDP_LoadBndP (BVP * theBVP, HEAP * Heap)
3138 {
3139 BND_PS *bp;
3140 int i, j, pid, n;
3141 int iList[2];
3142 double dList[DIM];
3143
3144 if (Bio_Read_mint (2, iList))
3145 return (NULL);
3146 pid = iList[0];
3147 n = iList[1];
3148 bp =
3149 (BND_PS *) GetFreelistMemory (Heap,
3150 (n - 1) * sizeof (COORD_BND_VECTOR) +
3151 sizeof (BND_PS));
3152 bp->n = n;
3153 bp->patch_id = pid;
3154 for (i = 0; i < n; i++)
3155 {
3156 if (Bio_Read_mdouble (DIM - 1, dList))
3157 return (NULL);
3158 for (j = 0; j < DIM - 1; j++)
3159 bp->local[i][j] = dList[j];
3160 }
3161 /* load free boundary points properly */
3162 if (!PATCH_IS_FIXED (currBVP->patches[pid]))
3163 {
3164 DOUBLE *pos;
3165
3166 /* read global coordinates */
3167 BND_DATA (bp) = GetFreelistMemory (Heap, DIM * sizeof (DOUBLE));
3168 if (BND_DATA (bp) == NULL)
3169 return (NULL);
3170
3171 if (Bio_Read_mdouble (DIM, dList))
3172 return (NULL);
3173 pos = (DOUBLE *) BND_DATA (bp);
3174 for (j = 0; j < DIM; j++)
3175 pos[j] = dList[j];
3176 }
3177
3178 return ((BNDP *) bp);
3179 }
3180
3181 BNDP *NS_DIM_PREFIX
BNDP_LoadBndP_Ext(void)3182 BNDP_LoadBndP_Ext (void)
3183 {
3184 BND_PS *bp;
3185 int i, j, pid, n;
3186 int iList[2];
3187 double dList[DIM - 1];
3188
3189 if (Bio_Read_mint (2, iList))
3190 return (NULL);
3191 pid = iList[0];
3192 n = iList[1];
3193 bp =
3194 (BND_PS *) malloc ((n - 1) * sizeof (COORD_BND_VECTOR) + sizeof (BND_PS));
3195 bp->n = n;
3196 bp->patch_id = pid;
3197 for (i = 0; i < n; i++)
3198 {
3199 if (Bio_Read_mdouble (DIM - 1, dList))
3200 return (NULL);
3201 for (j = 0; j < DIM - 1; j++)
3202 bp->local[i][j] = dList[j];
3203 }
3204
3205 return ((BNDP *) bp);
3206 }
3207
3208 /****************************************************************************/
3209 /** \brief Create and initialize the std_domain
3210 *
3211 * This function creates the environments domain, problem and BVP.
3212 *
3213 * @return <ul>
3214 * <li> 0 if ok
3215 * <li> 1 when error occured.
3216 * </ul>
3217 */
3218 /****************************************************************************/
3219
3220 INT NS_DIM_PREFIX
InitDom(void)3221 InitDom (void)
3222 {
3223
3224 /* change to root directory */
3225 if (ChangeEnvDir ("/") == NULL)
3226 {
3227 PrintErrorMessage ('F', "InitDom", "could not changedir to root");
3228 return (__LINE__);
3229 }
3230
3231 /* get env dir/var IDs for the problems */
3232 theProblemDirID = GetNewEnvDirID ();
3233 theBdryCondVarID = GetNewEnvVarID ();
3234
3235 /* install the /Domains directory */
3236 theDomainDirID = GetNewEnvDirID ();
3237 if (MakeEnvItem ("Domains", theProblemDirID, sizeof (ENVDIR)) == NULL)
3238 {
3239 PrintErrorMessage ('F', "InitDom", "could not install '/Domains' dir");
3240 return (__LINE__);
3241 }
3242 theBdrySegVarID = GetNewEnvVarID ();
3243 theLinSegVarID = GetNewEnvVarID ();
3244
3245 /* install the /BVP directory */
3246 theBVPDirID = GetNewEnvDirID ();
3247 if (MakeEnvItem ("BVP", theBVPDirID, sizeof (ENVDIR)) == NULL)
3248 {
3249 PrintErrorMessage ('F', "InitDom", "could not install '/BVP' dir");
3250 return (__LINE__);
3251 }
3252
3253 return (0);
3254 }
3255
3256 /** @} */
3257