1 /*
2 * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
3 * Applied Mathematics, Norway.
4 *
5 * Contact information: E-mail: tor.dokken@sintef.no
6 * SINTEF ICT, Department of Applied Mathematics,
7 * P.O. Box 124 Blindern,
8 * 0314 Oslo, Norway.
9 *
10 * This file is part of SISL.
11 *
12 * SISL is free software: you can redistribute it and/or modify
13 * it under the terms of the GNU Affero General Public License as
14 * published by the Free Software Foundation, either version 3 of the
15 * License, or (at your option) any later version.
16 *
17 * SISL is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU Affero General Public License for more details.
21 *
22 * You should have received a copy of the GNU Affero General Public
23 * License along with SISL. If not, see
24 * <http://www.gnu.org/licenses/>.
25 *
26 * In accordance with Section 7(b) of the GNU Affero General Public
27 * License, a covered work must retain the producer line in every data
28 * file that is created or manipulated using SISL.
29 *
30 * Other Usage
31 * You can be released from the requirements of the license by purchasing
32 * a commercial license. Buying such a license is mandatory as soon as you
33 * develop commercial activities involving the SISL library without
34 * disclosing the source code of your own applications.
35 *
36 * This file may be used in accordance with the terms contained in a
37 * written agreement between you and SINTEF ICT.
38 */
39
40 #include "sisl-copyright.h"
41
42 /*
43 *
44 * $Id: s1162.c,v 1.3 2001-03-19 15:58:41 afr Exp $
45 *
46 */
47
48
49 #define S1162
50
51 #include "sislP.h"
52
53 /*
54 * Forward declarations.
55 * ---------------------
56 */
57
58 #if defined(SISLNEEDPROTOTYPES)
59 static void s1162_s9mic(SISLObject *,SISLObject *,SISLIntdat **,
60 SISLEdge *[],int *);
61 static void s1162_s9num(SISLObject *,int *,int *);
62 static void s1162_s9edge(SISLObject *[],SISLObject *[],int,int,
63 SISLIntdat *,SISLEdge *[],int *);
64 static void s1162_s9con(SISLObject *,double *,double,SISLIntdat **,
65 SISLEdge *[],int *,int *,int *);
66 static void s1162_s9update(SISLObject *,double *,double,SISLIntdat **,
67 SISLEdge *[2],int *);
68 static void s1162_s9div(SISLObject *,double *,double,int,int,int,
69 SISLObject *[],SISLIntdat **,SISLEdge *[2],int,int *);
70 #else
71 static void s1162_s9mic();
72 static void s1162_s9num();
73 static void s1162_s9edge();
74 static void s1162_s9con();
75 static void s1162_s9update();
76 static void s1162_s9div();
77 #endif
78
79 #if defined(SISLNEEDPROTOTYPES)
80 void
s1162(SISLObject * po1,double * cmax,double aepsge,SISLIntdat ** pintdat,SISLEdge * vedge[2],int ilevel,int inum,int * jstat)81 s1162(SISLObject *po1,double *cmax,double aepsge,
82 SISLIntdat **pintdat,SISLEdge *vedge[2],
83 int ilevel,int inum,int *jstat)
84 #else
85 void s1162(po1,cmax,aepsge,pintdat,vedge,ilevel,inum,jstat)
86 SISLObject *po1;
87 double *cmax;
88 double aepsge;
89 SISLIntdat **pintdat;
90 SISLEdge *vedge[2];
91 int ilevel;
92 int inum;
93 int *jstat;
94 #endif
95 /*
96 *********************************************************************
97 *
98 *********************************************************************
99 *
100 * PURPOSE : SISLObject - object maximum. Treat the inner of the
101 * object.
102 *
103 *
104 *
105 * INPUT : po1 - Pointer to object
106 * aepsge - Geometry resolution.
107 * vedge - Pointers to structure of edge-maximums.
108 * vedge[1] must be SISL_NULL
109 * ilevel - Debt in recursion with inumb(>2) max. on the edges(if bezier case).
110 * inum - Number of max. on the edges.
111 *
112 * INPUT/OUTPUT : pintdat - Pointer to maximum data.
113 * cmax - Level value.
114 *
115 *
116 * OUTPUT : jstat - status messages
117 * = 2 : maximum found equal to level value
118 * = 1 : maximum found over to level value
119 * = 0 : no maximum
120 * < 0 : error
121 *
122 * METHOD :
123 *
124 *
125 * REFERENCES :
126 *
127 *-
128 * CALLS : s6err - Treating error situation.
129 * s1119 - Simple Case test for maximums.
130 * s1190 - Box test.
131 * s1791 - Test if possible to subdivide
132 * s1792 - Computing midpoint of parameter intervall.
133 * s1770 - Curve/curve iteration.
134 * s1770p - Point/curve iteration.
135 * s1771 - Curve/surface iteration.
136 * s1771p - Point/surface iteration.
137 * s1231 - Subdivide curve.
138 * s1711 - Subdivide surface.
139 * s1161 - Object/object maximum.
140 * s1435 - Pick an edge curve from a surface.
141 * s1438 - Pick an end point from a curve.
142 * s6idnpt -
143 * s6idkpt -
144 * s6idcpt -
145 * s6idcon -
146 * s6idput -
147 * s6idedg -
148 * s6idint -
149 * newPoint -
150 * newObject -
151 * newIntpt -
152 * newEdge -
153 * newIntdat -
154 * freeObject - Free space occupied by a given object
155 * freeCurve -
156 * freeEdge -
157 * freeIntdat -
158 *
159 * WRITTEN BY : Ulf J. Krystad, SI, 89-05.
160 *
161 *********************************************************************
162 */
163 {
164 int klevel; /* Local - Debt in recursion with. */
165 int knumedge; /* Local - Number of max. on the edges*/
166 int kpos = 0; /* Position of error. */
167 int kstat = 0; /* Local error status. */
168 int ksimple = 0; /* Local simple case status. */
169 int kdiv = 0; /* Parameter direction of subdivsion. */
170 int knum; /* Number of edges in subproblems. */
171 int ki; /* Counter. */
172 int kind1,kind2; /* Index two knots with multiplicity. */
173 SISLObject *uob1[4]; /* Pointers to subdivided object. */
174 SISLObject *qdum; /* Pointer to dummy object. */
175 SISLEdge **uedge=SISL_NULL; /* Pointer to array (to be allocated)
176 of edges to use in subproblems. */
177 SISLIntpt *up[2];
178 SISLPtedge *qpt0,*qpt1;
179
180 /*Init*/
181 knumedge = inum;
182 klevel = ilevel;
183
184 for (ki=0;ki<4;ki++) uob1[ki] = SISL_NULL;
185 if ((qdum = newObject(SISLPOINT)) == SISL_NULL) goto err101;
186
187 /* Initiate no maximum.*/
188 *jstat = 0;
189
190 /* Test if maximum is possible (perform box-test). */
191 s1190(po1,cmax,aepsge,&kstat);
192 if (kstat < 0) goto error;
193
194 /* We may have four different values on kstat.
195 kstat = 1 : The SISLbox is beyond level value.
196 kstat = 2 : The object is of constant value.
197 kstat = 3 : The object is beyond one of its corners.
198 kstat = 0 : No conclusion.*/
199
200 if (kstat == 1);
201
202 /* No max is possible */
203
204 else if (kstat == 2)
205 {
206 /* The geometry is of constant value. Since it is not taken by
207 the SISLbox test and the edges already are treated in s1161,
208 we just connect the point on the edges. */
209
210
211 if (vedge[0] != SISL_NULL && vedge[0]->iedge == 2)
212 {
213 /* Only curves has to do connect */
214
215 qpt0=vedge[0]->prpt[0];
216 qpt1=vedge[0]->prpt[1];
217 if (qpt0 != SISL_NULL && qpt1 != SISL_NULL)
218 {
219
220 up[0] = qpt0->ppt;
221 up[1] = qpt1->ppt;
222 s6idcon(pintdat,&up[0],&up[1],&kstat);
223 if (kstat<0) goto error;
224 }
225 }
226 }
227
228 else if (kstat == 3);
229
230
231 /* Maximum for the object is a corner value, it has been found
232 while treating the edges. */
233
234
235 else
236 {
237 /* Simple Case test (more than one maximum possible?) */
238 if(po1->iobj ==SISLCURVE)
239
240 s1119(po1->c1->ecoef,po1->c1->et,po1->c1->et,
241 po1->c1->ik,po1->c1->in,
242 1,1,&ksimple,&kind1,&kind2,&kstat);
243
244 else
245 s1119(po1->s1->ecoef,po1->s1->et1,po1->s1->et2,
246 po1->s1->ik1,po1->s1->in1,
247 po1->s1->ik2,po1->s1->in2,&ksimple,&kind1,&kind2,&kstat);
248 if (kstat < 0) goto error;
249
250 /* We may have three different values on ksimple.
251 ksimple = 0 : Not possible with interior max.
252 ksimple = 1 : Simpel case
253 ksimple = 2 : Not simpel case.*/
254
255 if (ksimple == 0)
256 *jstat = 0;
257
258 else if (ksimple == 1)
259 {
260 /* Simple Case, uppdate maximum list. */
261
262 s1162_s9update(po1,cmax,aepsge,pintdat,vedge,&kstat);
263 if (kstat < 0) goto error;
264 *jstat = kstat;
265
266 }
267 else
268 {
269 /* Check for interval maximum.*/
270
271 s1162_s9con(po1,cmax,aepsge,pintdat,vedge,&klevel,&knumedge,&kstat);
272 if (kstat < 0) goto error;
273
274 /* We may have two different values on kstat.
275 kstat = 0 : No intervall maximum.
276 kstat = 1 : More than 2 maximum found on the edges.
277 (bezier case only).
278 kstat = 2 : Intervall maximum found.
279 kstat = 3 : Simple case */
280
281 if (kstat == 3)
282 /* Simple Case, uppdate maximum list. */
283 {
284
285 s1162_s9update(po1,cmax,aepsge,pintdat,vedge,&kstat);
286 if (kstat < 0) goto error;
287 *jstat = kstat;
288 }
289
290 else if (kstat == 2)
291
292 *jstat = kstat; /*Uppdating maximum found. */
293
294 else
295 {
296 /* Find number of possible subdivision directions.
297 kdiv may have 4 difference values :
298 kdiv = 0 : Subdivision not possible.
299 kdiv = 1 : Subdivision in first parameter direction.
300 kdiv = 2 : Subdivision in second parameter direction.
301 kdiv = 3 : Subdivision in both parameter directions. */
302
303 s1162_s9num(po1,&kdiv,&kstat);
304 if (kstat < 0) goto error;
305
306
307 if(kdiv == 0)
308 {
309 /* Microcase in parameter plane.*/
310
311 s1162_s9mic(po1,qdum,pintdat,vedge,&kstat);
312 if (kstat < 0) goto error;
313 else *jstat = kstat;
314 }
315 else
316 {
317 /* We do not have simpel case and it is possible to
318 subdivide. We therfor subdivide and uppdate the
319 edge maximum and then do a recurcive call
320 to treat the sub problems. Curves are subdivided
321 into two, surfaces into four. We can therfor get
322 up to four recursive calls.*/
323
324 /* Computing total number of subobjects in sub problems. */
325 knum = (kdiv<3 ? 2:4);
326
327 /***** Treating objects on sub problems. *****/
328
329 if (kdiv > 0) /* New objects for subdivision of po1. */
330 {
331 for (ki=0;ki<knum;ki++)
332 {
333 if ((uob1[ki] = newObject(po1->iobj)) == SISL_NULL)
334 goto err101;
335
336 /*Initiate o1 pointer to point to top level object.*/
337
338 uob1[ki]->o1 = po1->o1;
339 }
340
341 /* Subdivide the po1 object. */
342
343 s1162_s9div(po1,cmax,aepsge,kdiv,kind1,kind2,
344 uob1,pintdat,vedge,klevel,&kstat);
345 if (kstat < 0) goto error;
346 *jstat = max(*jstat,kstat);
347
348 }
349
350 /***** Treating edges on sub problems. *****/
351
352
353 /* Making array of pointers to edge object
354 to the sub problems. */
355 if ((uedge = new0array(2*knum,SISLEdge *)) == SISL_NULL)
356 goto err101;
357
358 /* Making new edge object to sub problems. */
359 for (ki=0; ki<2*knum; ki+=2)
360 {
361
362 if ((uedge[ki] = newEdge(vedge[0]->iedge)) == SISL_NULL)
363 goto err101;
364 /* No edge for the dummy point: */
365 uedge[ki+1] = SISL_NULL;
366
367 }
368
369
370 /***** Recursion. *****/
371 for (ki=0;ki<knum;ki+=1)
372 {
373
374 /* Uppdate edge maximum on sub problems. */
375 s1162_s9edge(uob1+ki, &qdum, 1, 1, *pintdat,
376 uedge+2*ki, &kstat);
377 if (kstat < 0) goto error;
378
379 s1162(uob1[ki],cmax,aepsge,pintdat,
380 uedge+2*ki,klevel,knumedge,&kstat);
381 if (kstat < 0) goto error;
382 else *jstat = max(*jstat, kstat);
383 }
384 }
385 }
386 }
387 }
388
389 /* Intersections in the inner found. */
390
391 goto out;
392
393 /* Error in space allocation. */
394 err101: *jstat = -101;
395 s6err("s1162",*jstat,kpos);
396 goto out;
397
398 /* Error in lower level routine. */
399 error : *jstat = kstat;
400 s6err("s1162",*jstat,kpos);
401 goto out;
402
403 /* Free the space that is allocated. */
404
405 out:
406 if (qdum != SISL_NULL) freeObject(qdum);
407
408 for (ki=0;ki<4;ki++)
409 if (uob1[ki] != SISL_NULL) freeObject(uob1[ki]);
410
411 if (uedge != SISL_NULL)
412 {
413 /* 26.10.92 UJK/ BEOrd13969 */
414 /* for (ki=0;ki<knum;ki++) */
415 for (ki=0;ki<2*knum;ki++)
416 if (uedge[ki] != SISL_NULL) freeEdge(uedge[ki]);
417
418 freearray(uedge);
419 }
420 }
421
422 #if defined(SISLNEEDPROTOTYPES)
423 static void
s1162_s9mic(SISLObject * po1,SISLObject * po2,SISLIntdat ** rintdat,SISLEdge * vedge[],int * jstat)424 s1162_s9mic(SISLObject *po1,SISLObject *po2,SISLIntdat **rintdat,SISLEdge *vedge[],int *jstat)
425 #else
426 static void s1162_s9mic(po1,po2,rintdat,vedge,jstat)
427 SISLObject *po1;
428 SISLObject *po2;
429 SISLIntdat **rintdat;
430 SISLEdge *vedge[];
431 int *jstat;
432 #endif
433 /*
434 *********************************************************************
435 *
436 *********************************************************************
437 *
438 * PURPOSE : Treat intersection when it is not possible to
439 * subdivide any futher, and it is not simple case.
440 *
441 *
442 *
443 * INPUT : vedge[2] - SISLEdge intersection objects to the two
444 * objects in intersection problem.
445 *
446 *
447 *
448 * OUTPUT : rintdat - intersection data.
449 * jstat - status messages
450 * = 1 : Intersection found.
451 * = 0 : Intersection not found.
452 * < 0 : error.
453 *
454 *
455 * METHOD :
456 *
457 *
458 * REFERENCES :
459 *
460 *
461 * WRITTEN BY : Arne Laksaa, SI, 89-04.
462 *
463 *********************************************************************
464 */
465 {
466 int kpos = 0; /* Position of error. */
467 int kstat=0; /* Local error status. */
468 int kpoint; /* Number of intpt on edges. */
469 double *spar = SISL_NULL; /* Array to store parameter values. */
470 SISLIntpt **up = SISL_NULL; /* Array of poiners to intersection point. */
471
472
473 /* Initiate to now new intersection point. */
474
475
476 *jstat = 0;
477
478
479 /* Compute number of intersection points on edges. */
480
481 if (vedge[0] == SISL_NULL )
482 kpoint = 0;
483 else
484 kpoint = vedge[0]->ipoint;
485
486 if (vedge[1] != SISL_NULL )
487 kpoint += vedge[1]->ipoint;
488
489
490 if (kpoint == 0 )
491 {
492 int kpar = 0;
493 SISLIntpt *qt;
494
495
496 /* There is not any intersection points on the edges.
497 We therfor make one new intersection point with parameter
498 values in senter of each object. */
499
500
501 /* Number of parameter values of object 1. */
502
503 if (po1->iobj == SISLCURVE) kpar = 1;
504 else if (po1->iobj == SISLSURFACE) kpar = 2;
505
506
507 /* Number of parameter values of object 2. */
508
509 if (po2->iobj == SISLCURVE) kpar++;
510 else if (po2->iobj == SISLSURFACE) kpar += 2;
511
512
513 /* Allocate array to store midpoint parameter values. */
514
515 if ((spar = newarray(kpar,double)) == SISL_NULL)
516 goto err101;
517
518
519 /* Compute midpoint parameter values. */
520
521 if (po1->iobj == SISLCURVE)
522 {
523 spar[0] = (po1->c1->et[po1->c1->ik - 1] +
524 po1->c1->et[po1->c1->in])*(double)0.5;
525 kpar = 1;
526 }
527 else if (po1->iobj == SISLSURFACE)
528 {
529 spar[0] = (po1->s1->et1[po1->s1->ik1 - 1] +
530 po1->s1->et1[po1->s1->in1])*(double)0.5;
531 spar[1] = (po1->s1->et2[po1->s1->ik2 - 1] +
532 po1->s1->et2[po1->s1->in2])*(double)0.5;
533 kpar = 2;
534 }
535
536 if (po2->iobj == SISLCURVE)
537 {
538 spar[kpar] = (po2->c1->et[po2->c1->ik - 1] +
539 po2->c1->et[po2->c1->in])*(double)0.5;
540 kpar++;
541 }
542 else if (po2->iobj == SISLSURFACE)
543 {
544 spar[kpar] = (po2->s1->et1[po2->s1->ik1 - 1] +
545 po2->s1->et1[po2->s1->in1])*(double)0.5;
546 spar[kpar+1] = (po2->s1->et2[po2->s1->ik2 - 1] +
547 po2->s1->et2[po2->s1->in2])*(double)0.5;
548 kpar += 2;
549 }
550
551 *jstat = 1; /* Mark intersection found. */
552
553
554 /* Makeing intersection point. */
555
556 qt = newIntpt(kpar,spar,DZERO);
557 if (qt == SISL_NULL) goto err101;
558
559 /* Uppdating pintdat. */
560
561 s6idnpt(rintdat,&qt,1,&kstat);
562 if (kstat < 0) goto error;
563 }
564 else if (kpoint > 1)
565 {
566 int kn,kn1,ki,kj;
567 SISLPtedge *qpt;
568
569
570 /* We have more than one intersection point on the edges,
571 we therfor conect these points to each other. */
572
573 /* Allacate array of pointers to these points. */
574
575 if ((up = newarray(kpoint,SISLIntpt *)) == SISL_NULL) goto err101;
576
577
578 /* Uppdate the array. */
579
580 for (kn=0,kn1=0; kn<2; kn++)
581 if (vedge[kn] != SISL_NULL && vedge[kn]->ipoint > 0)
582 for(kj=0; kj<vedge[kn]->iedge; kj++)
583 for(qpt=vedge[kn]->prpt[kj]; qpt != SISL_NULL; qpt=qpt->pnext,kn1++)
584 up[kn1] = qpt->ppt;
585
586
587 /* Connect the points to each other. */
588
589 for (ki=1; ki<kpoint; ki++)
590 {
591 s6idcon(rintdat,&up[ki-1],&up[ki],&kstat);
592 if (kstat<0) goto error;
593 }
594 }
595
596 goto out;
597
598 /* Error in space allocation. */
599
600 err101: *jstat = -101;
601 s6err("s1162_s9mic",*jstat,kpos);
602 goto out;
603
604 /* Error in lower level routine. */
605
606 error : *jstat = kstat;
607 s6err("s1162_s9mic",*jstat,kpos);
608 goto out;
609
610 out: if (spar != SISL_NULL) freearray(spar);
611 if (up != SISL_NULL) freearray(up);
612 }
613
614 #if defined(SISLNEEDPROTOTYPES)
615 static void
s1162_s9num(SISLObject * po,int * jdiv,int * jstat)616 s1162_s9num(SISLObject *po,int *jdiv,int *jstat)
617 #else
618 static void s1162_s9num(po,jdiv,jstat)
619 SISLObject *po;
620 int *jdiv;
621 int *jstat;
622 #endif
623 /*
624 *********************************************************************
625 *
626 *********************************************************************
627 *
628 * PURPOSE : Find number of possible subdivisions of an object.
629 *
630 *
631 *
632 * INPUT : po - SISLObject to subdevide.
633 *
634 *
635 *
636 * OUTPUT : jdiv - Possible subdivisions of object.
637 * = 0 : No subdivision.
638 * = 1 : Subdivision in first parameter direction.
639 * = 2 : Subdivision in second parameter direction.
640 * = 3 : Subdivision in both parameter direction.
641 * jstat - status messages
642 * = 0 : no error.
643 * < 0 : error
644 *
645 *
646 * METHOD :
647 *
648 *
649 * REFERENCES :
650 *
651 *
652 * WRITTEN BY : Arne Laksaa, SI, 89-04.
653 * Revised BY : Christophe Rene Birkeland, SINTEF Oslo, May 93
654 * *jstat = 0 initialization
655 *
656 *********************************************************************
657 */
658 {
659 *jstat = 0;
660 if (po->iobj == SISLPOINT) *jdiv = 0;
661 else if (po->iobj == SISLCURVE)
662 {
663 if(s1791(po->c1->et,po->c1->ik,po->c1->in)) *jdiv = 1;
664 else *jdiv = 0;
665 }
666 else if (po->iobj == SISLSURFACE)
667 {
668 if(s1791(po->s1->et1,po->s1->ik1,po->s1->in1)) *jdiv = 1;
669 else *jdiv = 0;
670
671 if(s1791(po->s1->et2,po->s1->ik2,po->s1->in2)) *jdiv += 2;
672 }
673 else
674 {
675
676 /* Error. Kind of object does not exist. */
677
678 *jstat = -121;
679 s6err("s1162_s9num",*jstat,0);
680 }
681 }
682
683 #if defined(SISLNEEDPROTOTYPES)
684 static void
s1162_s9edge(SISLObject * vob1[],SISLObject * vob2[],int iobj1,int iobj2,SISLIntdat * pintdat,SISLEdge * wedge[],int * jstat)685 s1162_s9edge(SISLObject *vob1[],SISLObject *vob2[],
686 int iobj1,int iobj2,SISLIntdat *pintdat,SISLEdge *wedge[],int *jstat)
687 #else
688 static void s1162_s9edge(vob1,vob2,iobj1,iobj2,pintdat,wedge,jstat)
689 SISLObject *vob1[];
690 SISLObject *vob2[];
691 int iobj1;
692 int iobj2;
693 SISLIntdat *pintdat;
694 SISLEdge *wedge[];
695 int *jstat;
696 #endif
697 /*
698 *********************************************************************
699 *
700 *********************************************************************
701 *
702 * PURPOSE : Uppdate edge structure. It may be up to four object
703 * in vob1 and vob2, and wedge may therfor contain
704 * seexteen elements to be uppdated.
705 *
706 *
707 *
708 * INPUT : vob1[] - Array of pointers to first objects.
709 * vob2[] - Array of pointers to second objects.
710 * iobj1 - Number of elements in vob1.
711 * iobj2 - Number of elements in vob2.
712 * pintdat - intersection data.
713 *
714 *
715 *
716 * OUTPUT : wedge[] - SISLEdge structures to uppdate.
717 * jstat - status messages
718 * = 0 : OK!
719 * < 0 : error.
720 *
721 *
722 * METHOD :
723 *
724 *
725 * REFERENCES :
726 *
727 *
728 * WRITTEN BY : Arne Laksaa, SI, 89-05.
729 *
730 *********************************************************************
731 */
732 {
733 int kpos = 0; /* Position of error. */
734 int kstat=0; /* Local error status. */
735 int ki1,ki2,kj,kn; /* Counters. */
736 int kedg; /* Number of edges. */
737 int kpar; /* Parameter number. */
738 double tpar; /* Parameter value at edge. */
739
740
741 for (kn=0,ki1=0; ki1<iobj1; ki1++)
742 for (ki2=0; ki2<iobj2; ki2++,kn+=2)
743 {
744 kedg = (vob1[ki1]->iobj == SISLPOINT ?0:(vob1[ki1]->iobj == SISLCURVE ?2:4));
745
746 for (kj=0; kj<kedg; kj++)
747 {
748 if (vob1[ki1]->iobj == SISLCURVE)
749 {
750 tpar = (kj == 0 ? vob1[ki1]->c1->et[vob1[ki1]->c1->ik-1] :
751 vob1[ki1]->c1->et[vob1[ki1]->c1->in]);
752 kpar = 1;
753 }
754 else if (kj == 0)
755 {
756 tpar = vob1[ki1]->s1->et2[vob1[ki1]->s1->ik2-1];
757 kpar = 2;
758 }
759 else if (kj == 1)
760 {
761 tpar = vob1[ki1]->s1->et1[vob1[ki1]->s1->in1];
762 kpar = 1;
763 }
764 else if (kj == 2)
765 {
766 tpar = vob1[ki1]->s1->et2[vob1[ki1]->s1->in2];
767 kpar = 2;
768 }
769 else
770 {
771 tpar = vob1[ki1]->s1->et1[vob1[ki1]->s1->ik1-1];
772 kpar = 1;
773 }
774
775
776 s6idedg(vob1[ki1],vob2[ki2],1,kpar,tpar,pintdat,
777 &(wedge[kn]->prpt[kj]),&(wedge[kn]->ipoint),&kstat);
778 if (kstat < 0) goto error;
779 }
780
781 kedg = (vob2[ki2]->iobj == SISLPOINT ?0:(vob2[ki2]->iobj == SISLCURVE ?2:4));
782
783 for (kj=0; kj<kedg; kj++)
784 {
785 if (vob2[ki2]->iobj == SISLCURVE)
786 {
787 tpar = (kj == 0 ? vob2[ki2]->c1->et[vob2[ki2]->c1->ik-1] :
788 vob2[ki2]->c1->et[vob2[ki2]->c1->in]);
789 kpar = 1;
790 }
791 else if (kj == 0)
792 {
793 tpar = vob2[ki2]->s1->et2[vob2[ki2]->s1->ik2-1];
794 kpar = 2;
795 }
796 else if (kj == 1)
797 {
798 tpar = vob2[ki2]->s1->et1[vob2[ki2]->s1->in1];
799 kpar = 1;
800 }
801 else if (kj == 2)
802 {
803 tpar = vob2[ki2]->s1->et2[vob2[ki2]->s1->in2];
804 kpar = 2;
805 }
806 else
807 {
808 tpar = vob2[ki2]->s1->et1[vob2[ki2]->s1->ik1-1];
809 kpar = 1;
810 }
811
812
813 s6idedg(vob1[ki1],vob2[ki2],2,kpar,tpar,pintdat,
814 &(wedge[kn+1]->prpt[kj]),&(wedge[kn+1]->ipoint),&kstat);
815 if (kstat < 0) goto error;
816 }
817 }
818
819 *jstat = 0;
820
821 goto out;
822
823 /* Error in lower level routine. */
824
825 error : *jstat = kstat;
826 s6err("s1162_s9edge",*jstat,kpos);
827 goto out;
828
829 out: ;
830 }
831
832 #if defined(SISLNEEDPROTOTYPES)
833 static void
s1162_s9con(SISLObject * po1,double * cmax,double aepsge,SISLIntdat ** pintdat,SISLEdge * vedge[],int * jlevel,int * jnum,int * jstat)834 s1162_s9con(SISLObject *po1,double *cmax,double aepsge,SISLIntdat **pintdat,
835 SISLEdge *vedge[],int *jlevel,int *jnum,int *jstat)
836 #else
837 static void s1162_s9con(po1,cmax,aepsge,pintdat,vedge,jlevel,jnum,jstat)
838 SISLObject *po1;
839 double *cmax;
840 double aepsge;
841 SISLIntdat **pintdat;
842 SISLEdge *vedge[];
843 int *jlevel;
844 int *jnum;
845 int *jstat;
846 #endif
847 /*
848 *********************************************************************
849 *
850 *********************************************************************
851 *
852 * PURPOSE : Check if we are able to connect the maxpoints of the object
853 * found on the edges.
854 *
855 *
856 *
857 * INPUT : po1 - The first SISLObject to check.
858 * cmax - The max value found.
859 * aepsge - Geometrical resolution.
860 * vedge[] - SISLEdge intersection.
861 *
862 *
863 * INPUT/OUTPUT:pintdat - Intersection data.
864 * jlevel - Debt in recursion with inumb(>2) max. on the edges(if bezier case).
865 * jnum - Number of max. on the edges.
866 *
867 *
868 * OUTPUT :
869 * jstat - status messages
870 * = 3 : Simple Case.
871 * = 2 : Connection done.
872 * = 1 : Level value set to one or increased by one.
873 * = 0 : Level value set to zero.
874 * < 0 : error
875 *
876 *
877 * METHOD : When - object is a surface
878 * - the surf is a bezier patch
879 * - the patch has three (or more) max on the edges
880 * If jlevel = 4 and number of max on the edges are equal to jnum,
881 * we connect the points treating one of them as singulear.
882 * If jlevel in {1,2,3} and number of max
883 * on the edges are equal to jnum,
884 * we increase the jlevel by one.
885 * If jlevel is 0 we set jlevel to one and
886 * jnum to number of max on edge.
887 * If jnum not equal to number of max on edge,
888 * we set jlevel to one and
889 * jnum to number of max on edge..
890 *
891 * If the number of points on the edge is > 10, we do nothing.
892 * REFERENCES :
893 *
894 *
895 * WRITTEN BY : Ulf J Krystad, SI, 89-05.
896 *
897 *********************************************************************
898 */
899 {
900 SISLIntpt *qintpt,*up[10];
901 SISLPtedge *qpt;
902
903 int kstat; /* Local status. */
904 /*guen int kpos; */ /* Local status counter. */
905 /*guen changed into:*/
906 int kpos = 0; /* Local status counter. */
907 int kk1; /* Local SURFACE attribute. */
908 int kk2; /* Local SURFACE attribute. */
909 int kn1; /* Local SURFACE attribute. */
910 int kn2; /* Local SURFACE attribute. */
911 int ki,kj; /* Local counter. */
912 int kfound; /* Local flag in loop. */
913 int knum = 0; /* Local number of max on edge. */
914 int klevel = 0; /* Local level. */
915 int kleft1 = 0; /* Local input parameter s1421 */
916 int kleft2 = 0; /* Local input parameter s1421 */
917 int kder = 1; /* Local input parameter s1421 */
918
919 double spar[2]; /* Parameter value */
920 double spar1[2]; /* Parameter value */
921 double smidle[2];/* middle parameter value */
922 double *sval=SISL_NULL;/* Values from s1421. */
923 double *snorm=SISL_NULL;/* Values from s1421. */
924
925
926 kstat = 0;
927
928 if (po1->iobj == SISLSURFACE)
929 {
930 if ((po1->s1->in1 == po1->s1->ik1) && (po1->s1->in2 == po1->s1->ik2))
931 /* Bezier case for surface */
932 {
933
934 /*-------------------------------------------------------*/
935 /* Count number of max on the edges. */
936 kk1 = po1->s1->ik1;
937 kk2 = po1->s1->ik2;
938 kn1 = po1->s1->in1;
939 kn2 = po1->s1->in2;
940
941 for (kj=0,knum=0;kj<vedge[0]->iedge;kj++)
942 {
943 qpt = vedge[0]->prpt[kj];
944 while(qpt != SISL_NULL)
945 {
946 qintpt = qpt->ppt;
947 for (ki=0,kfound=0;ki<knum && kfound == 0;ki++)
948 if (qintpt == up[ki]) kfound = 1;
949
950 if (kfound == 0)
951 {
952 if (knum > 9) goto out;
953 up[knum]=qintpt;
954 knum++;
955 }
956
957 qpt = qpt->pnext;
958 }
959
960 }
961
962 /*---------------------------------------------------------*/
963
964 if (knum > 0 )
965 /* Number of max on the edges more than 1. */
966 {
967 klevel = *jlevel;
968
969 if (klevel == 0 || knum !=*jnum)
970 /* No continuation of suspected singulear point,
971 start a new one. */
972 {
973 kstat = 1;
974 klevel = 1;
975 }
976 else if (klevel < 2)
977 /* Continuation of suspected singulear point. */
978 {
979 kstat = 1;
980 klevel += 1;
981 }
982 else if (knum < 2 )
983 /* Simple Case */
984 {
985 kstat = 3;
986 klevel += 1;
987 }
988 else
989 {
990
991 /*--------------------------------------------------*/
992 /* Connection case. */
993
994 /* Allocate local used memory */
995
996 sval = newarray(4,double);
997 if (sval == SISL_NULL) goto err101;
998 snorm = sval + 3;
999
1000 for (kj=0;kj<knum-1;kj++)
1001 {
1002 spar[0] = up[kj]->epar[0];
1003 spar[1] = up[kj]->epar[1];
1004
1005 for (ki=kj+1;ki<knum;ki++)
1006 {
1007 /* First we linearize. */
1008 spar1[0] = up[ki]->epar[0];
1009 spar1[1] = up[ki]->epar[1];
1010 smidle[0] = (spar[0] + spar1[0])/(double)2.0;
1011 smidle[1] = (spar[1] + spar1[1])/(double)2.0;
1012
1013 /* Evaluate 0-1.st derivatives of surface */
1014
1015 s1421(po1->s1,kder,smidle,&kleft1,&kleft2,
1016 sval,snorm,&kstat);
1017 if (kstat < 0) goto error;
1018 if (fabs(sval[0]-*cmax) < aepsge)
1019 {
1020 /* Connect. */
1021 s6idcon(pintdat,&up[kj],&up[ki],&kstat);
1022 if (kstat<0) goto error;
1023 }
1024
1025 }
1026 }
1027
1028
1029 kstat = 2;
1030 /*------------------------------------------*/
1031
1032
1033 }
1034 }
1035 }
1036 }
1037
1038 goto out;
1039
1040 /* Error in allocation */
1041 err101: kstat = -101;
1042 s6err("s1162_s9con",kstat,kpos);
1043 goto out;
1044
1045 /* Error in lower level routine. */
1046 error : s6err("s1162_s9con",kstat,kpos);
1047 goto out;
1048
1049 out: if (sval != SISL_NULL) freearray(sval);
1050 *jlevel = klevel;
1051 *jnum = knum;
1052 *jstat = kstat;
1053
1054 }
1055
1056 #if defined(SISLNEEDPROTOTYPES)
1057 static void
s1162_s9update(SISLObject * po1,double * cmax,double aepsge,SISLIntdat ** pintdat,SISLEdge * vedge[2],int * jstat)1058 s1162_s9update(SISLObject *po1,double *cmax,double aepsge,
1059 SISLIntdat **pintdat,SISLEdge *vedge[2],int *jstat)
1060 #else
1061 static void s1162_s9update(po1,cmax,aepsge,pintdat,vedge,jstat)
1062 SISLObject *po1;
1063 double *cmax;
1064 double aepsge;
1065 SISLIntdat **pintdat;
1066 SISLEdge *vedge[2];
1067 int *jstat;
1068 #endif
1069 /*
1070 *********************************************************************
1071 *
1072 *********************************************************************
1073 *
1074 * PURPOSE : To find an inner maxima in an object when
1075 * it has no more then one.
1076 *
1077 * INPUT : po1 - The object.
1078 * aepsge - Geometry resolution.
1079 * vedge - Pointer to edge maximum.
1080 *
1081 *
1082 * INPUT/OUTPUT : cmax - The level value.
1083 *
1084 * OUTPUT : pintdat - Maximum data.
1085 * jstat - Status messages
1086 * = 2 : New maximum found inside edges.
1087 * = 1 : Maximum found inside object.
1088 * = 0 : no maximum found.
1089 * < 0 : error
1090 *
1091 *
1092 * METHOD :
1093 *
1094 *
1095 * REFERENCES :
1096 *
1097 *
1098 * WRITTEN BY : Ulf J. Krystad, SI, 89-05.
1099 *
1100 *********************************************************************
1101 */
1102 {
1103 int i, kj, ki; /* Counters. */
1104 int kpos = 0; /* Position of error. */
1105 int kstat= 0; /* Local status */
1106 int kk1, kk2, kn1, kn2; /* Local number of knots and vertices. */
1107 int kmax, kind1, kind2; /* Indexes of the maximum vertice. */
1108 int kleft = 0; /* Used in s1221 . */
1109 int kleft2 = 0; /* Used in s1424 . */
1110 int kconn = 0; /* Connection flag. */
1111 int knum = 0; /* Number of max on the edge. */
1112 int kfound = 0; /* Flag. */
1113
1114 double tstart, tend; /* Start, end values for curve parameter. */
1115 double sstart[2], send[2]; /* Start, end values for surface parameter. */
1116 double tpar; /* The parameter vallue for
1117 subdivision of a curve. */
1118 double spar[2]; /* The parameter vallue for subdivision
1119 of a surface. */
1120 double tmax, tmin; /* Local max and min value for the
1121 vertices of object. */
1122 double tval; /* The value of the geometry at the found point.*/
1123
1124
1125 SISLObject *qop=SISL_NULL,*qcuo = SISL_NULL;/* Help pointers */
1126 SISLIntdat *qintdat=SISL_NULL; /* Local max data. */
1127 SISLIntdat *qintdat1=SISL_NULL; /* Local for double upgrading. */
1128 SISLIntpt *qintpt,*up[3];
1129 SISLPtedge *qpt;
1130
1131 /* Init */
1132 *jstat = 0;
1133 if (po1 == SISL_NULL || po1->iobj == SISLPOINT) goto out;
1134 if ((qop = newObject(SISLPOINT)) == SISL_NULL) goto err101;
1135
1136 if (po1->iobj == SISLCURVE)
1137 {
1138 kk1 = po1->c1->ik;
1139 kn1 = po1->c1->in;
1140 kmax = po1->c1->pbox->imax;
1141 tmax = po1->c1->pbox->emax[0];
1142 tmin = po1->c1->pbox->emin[0];
1143
1144 tstart = po1->c1->et[kk1-1];
1145 tend = po1->c1->et[kn1];
1146
1147 /* Try to find an inner ekstremal point by iteration. */
1148
1149
1150 /* First get a good starting point for the iteration. */
1151 tpar = 0;
1152 for (i=kmax+1;i<kmax+kk1;i++)
1153 tpar += po1->c1->et[i];
1154
1155 tpar /= kk1 - 1;
1156
1157 s1252(po1->c1,aepsge,tpar,&tpar,&kstat);
1158 if (kstat < 0) goto error;
1159
1160 /* Test if the found point is at start or end. */
1161 if(DEQUAL(tpar,tstart) || DEQUAL(tpar,tend)) goto out;
1162
1163
1164 /* Evaluate curve at parameter value. */
1165 kleft = 0;
1166 s1221(po1->o1->c1,0,tpar,&kleft,&tval,&kstat);
1167 if (kstat < 0) goto error;
1168
1169 /* Here we are ready to examine if we really found a new max point.*/
1170 if ((qop->p1 = newPoint(&tval,1,1)) == SISL_NULL) goto err101;
1171
1172 s1161(qop,cmax,aepsge,&qintdat,&kstat);
1173 if (kstat < 0) goto error;
1174
1175 if (kstat == 2)
1176 /* New maximum found, delete old ones */
1177 if (*pintdat != SISL_NULL)
1178 {
1179 freeIntdat(*pintdat);
1180 *pintdat = SISL_NULL;
1181 }
1182
1183 if ( kstat )
1184 {
1185 /* Maximum found, add it to the list */
1186 *jstat = max(*jstat,kstat); /* Mark maximum found. */
1187
1188 /* Set parameter parameter value of curve. */
1189 s6idput(pintdat,qintdat,0,tpar,&kstat);
1190 if (kstat < 0) goto error;
1191 }
1192
1193 }
1194 else if (po1->iobj == SISLSURFACE)
1195 {
1196
1197
1198 kk1 = po1->s1->ik1;
1199 kn1 = po1->s1->in1;
1200 kk2 = po1->s1->ik2;
1201 kn2 = po1->s1->in2;
1202 kmax = po1->s1->pbox->imax;
1203 tmax = po1->s1->pbox->emax[0];
1204 tmin = po1->s1->pbox->emin[0];
1205
1206 sstart[0] = po1->s1->et1[kk1-1];
1207 sstart[1] = po1->s1->et2[kk2-1];
1208
1209 send[0] = po1->s1->et1[kn1];
1210 send[1] = po1->s1->et2[kn2];
1211
1212
1213 /* Get the two dimensional index of the greatest vertice. */
1214 kind2 = kmax/kk1;
1215 kind1 = kmax - kind2*kk1;
1216
1217
1218 /*-----------------------------------------------------------*/
1219 /* Count number of max on the edges. */
1220
1221 for (kj=0,knum=0;kj<vedge[0]->iedge&&knum<3;kj++)
1222 {
1223 qpt = vedge[0]->prpt[kj];
1224 while(qpt != SISL_NULL && knum<3)
1225 {
1226 qintpt = qpt->ppt;
1227 for (ki=0,kfound=0;ki<knum && kfound == 0;ki++)
1228 if (qintpt == up[ki]) kfound = 1;
1229
1230 if (kfound == 0)
1231 {
1232 up[knum]=qintpt;
1233 knum++;
1234 }
1235
1236 qpt = qpt->pnext;
1237 }
1238
1239 }
1240
1241 /* Try if connection is possible.*/
1242 if (knum == 2)
1243 {
1244 /* if on same edge, they are be connected before
1245 (when in simple case.)*/
1246 if ((DEQUAL(up[0]->epar[0],sstart[0]) &&
1247 DEQUAL(up[1]->epar[0],sstart[0]))||
1248 (DEQUAL(up[0]->epar[0],send[0]) &&
1249 DEQUAL(up[1]->epar[0],send[0])) ||
1250 (DEQUAL(up[0]->epar[1],sstart[1]) &&
1251 DEQUAL(up[1]->epar[1],sstart[1]))||
1252 (DEQUAL(up[0]->epar[1],send[1]) &&
1253 DEQUAL(up[1]->epar[1],send[1])))
1254 kconn = 0;
1255
1256 else
1257 {
1258 /* Pick out two curves between the parameter
1259 value on the edges. */
1260 kconn = 0;
1261 ki = 0;
1262 if (fabs(up[0]->epar[0]-up[1]->epar[0]) <
1263 fabs(up[0]->epar[1]-up[1]->epar[1]))
1264 ki =1;
1265
1266 tpar = (double)0.25*up[0]->epar[ki] +
1267 (double)0.75*up[1]->epar[ki];
1268 if ((qcuo = newObject(SISLCURVE)) == SISL_NULL) goto err101;
1269 if (ki==0)
1270 s1437(po1->s1,tpar,&(qcuo->c1),&kstat);
1271 else
1272 s1436(po1->s1,tpar,&(qcuo->c1),&kstat);
1273 if (kstat < 0) goto error;
1274
1275 s1161(qcuo,cmax,aepsge,&qintdat,&kstat);
1276 if (kstat < 0) goto error;
1277
1278 if (kstat == 1)
1279 {
1280 freeCurve(qcuo->c1);
1281 qcuo->c1 = SISL_NULL;
1282
1283 tpar = (double)0.75*up[0]->epar[ki] +
1284 (double)0.25*up[1]->epar[ki];
1285 if (ki==0)
1286 s1437(po1->s1,tpar,&(qcuo->c1),&kstat);
1287 else
1288 s1436(po1->s1,tpar,&(qcuo->c1),&kstat);
1289 if (kstat < 0) goto error;
1290
1291 s1161(qcuo,cmax,aepsge,&qintdat,&kstat);
1292 if (kstat < 0) goto error;
1293
1294 if (kstat == 1)
1295 {
1296 /* Connect. */
1297 kconn = 1;
1298 s6idcon(pintdat,&up[0],&up[1],&kstat);
1299 if (kstat<0) goto error;
1300 }
1301 }
1302 }
1303 }
1304
1305
1306 if (kconn == 0)
1307 {
1308 /* No connection is done. */
1309
1310 /* Try to find an inner ekstremal point by iteration. */
1311
1312 /* First get a good starting point for the iteration. */
1313 spar[0] = 0;
1314 for (i=kind1+1;i<kind1+kk1;i++)
1315 spar[0] += po1->s1->et1[i];
1316
1317 spar[0] /= kk1 - 1;
1318
1319 spar[1] = 0;
1320 for (i=kind2+1;i<kind2+kk2;i++)
1321 spar[1] += po1->s1->et2[i];
1322
1323 spar[1] /= kk2 - 1;
1324
1325
1326 /* Create a point greater than the surface */
1327 if ((qop->p1 = newPoint(&tmax,1,1)) == SISL_NULL) goto err101;
1328
1329 /* Iterate using aepsge=tmax-tmin to ensure covergence. */
1330 s1173(qop->p1,po1->o1->s1,aepsge,sstart,send,spar,spar,&kstat);
1331 if (kstat < 0) goto error;
1332
1333 /* Test if the found point is at start or end. */
1334 if(DEQUAL(spar[0],sstart[0]) ||
1335 DEQUAL(spar[0],send[0]) ||
1336 DEQUAL(spar[1],sstart[1]) ||
1337 DEQUAL(spar[1],send[1])) goto out;
1338
1339 /* Evaluate surface at parameter value. */
1340 kleft = 0;
1341 kleft2 = 0;
1342 s1424(po1->o1->s1,0,0,spar,&kleft,&kleft2,&tval,&kstat);
1343 if (kstat < 0) goto error;
1344
1345 /* Here we are ready to examine if we really found a max point.*/
1346 freePoint(qop->p1);
1347 qop->p1 = SISL_NULL;
1348 if ((qop->p1 = newPoint(&tval,1,1)) == SISL_NULL) goto err101;
1349
1350 s1161(qop,cmax,aepsge,&qintdat,&kstat);
1351 if (kstat < 0) goto error;
1352
1353 if (kstat == 2)
1354 /* New maximum found, delete old ones */
1355 if (*pintdat != SISL_NULL)
1356 {
1357 freeIntdat(*pintdat);
1358 *pintdat = SISL_NULL;
1359 }
1360
1361 if ( kstat )
1362 {
1363 /* Maximum found, add them to the list */
1364
1365 *jstat = max(*jstat,kstat); /* Mark maximum found. */
1366
1367 /* Special treatment for putting two
1368 new parameters into pintdat from qintdat. */
1369 s6idput(&qintdat1,qintdat,0,spar[0],&kstat);
1370 if (kstat < 0) goto error;
1371 s6idput(pintdat,qintdat1,1,spar[1],&kstat);
1372 if (kstat < 0) goto error;
1373
1374 }
1375 }
1376 }
1377
1378
1379 goto out;
1380
1381 /* -------------------ERROR SECTION----------------------------*/
1382
1383 /* Error in space allocation. */
1384 err101: *jstat = -101;
1385 s6err("s1162_s9update",*jstat,kpos);
1386 goto out;
1387
1388 /* Error in lower level routine. */
1389 error : *jstat = kstat;
1390
1391 out:
1392 if (qcuo != SISL_NULL)
1393 {
1394 if (qcuo->c1 != SISL_NULL)
1395 {
1396 freeCurve(qcuo->c1);
1397 qcuo->c1 = SISL_NULL;
1398 }
1399 freeObject(qcuo);
1400 qcuo = SISL_NULL;
1401 }
1402
1403 if (qop != SISL_NULL)
1404 {
1405 if (qop->p1 != SISL_NULL)
1406 {
1407 freePoint(qop->p1);
1408 qop->p1 = SISL_NULL;
1409 }
1410 freeObject(qop);
1411 qop = SISL_NULL;
1412 }
1413 if (qintdat != SISL_NULL)
1414 {
1415 freeIntdat(qintdat);
1416 qintdat = SISL_NULL;
1417 }
1418 if (qintdat1 != SISL_NULL)
1419 {
1420 freeIntdat(qintdat1);
1421 qintdat1 = SISL_NULL;
1422 }
1423 }
1424
1425 #if defined(SISLNEEDPROTOTYPES)
1426 static void
s1162_s9div(SISLObject * po1,double * cmax,double aepsge,int idiv,int iind1,int iind2,SISLObject * wob[],SISLIntdat ** pintdat,SISLEdge * vedge[2],int ilevel,int * jstat)1427 s1162_s9div(SISLObject *po1,double *cmax,double aepsge,int idiv,int iind1,
1428 int iind2,SISLObject *wob[],SISLIntdat **pintdat,SISLEdge *vedge[2],
1429 int ilevel,int *jstat)
1430 #else
1431 static void s1162_s9div(po1,cmax,aepsge,idiv,iind1,iind2,
1432 wob,pintdat,vedge,ilevel,jstat)
1433 SISLObject *po1;
1434 double *cmax;
1435 double aepsge;
1436 int idiv;
1437 int iind1;
1438 int iind2;
1439 SISLObject *wob[];
1440 SISLIntdat **pintdat;
1441 SISLEdge *vedge[2];
1442 int ilevel;
1443 int *jstat;
1444 #endif
1445 /*
1446 *********************************************************************
1447 *
1448 *********************************************************************
1449 *
1450 * PURPOSE : Subdivide an object possible.
1451 *
1452 *
1453 *
1454 * INPUT : po1 - The object to subdivide.
1455 * aepsge - Geometry resolution.
1456 * idiv - Subdivision direction.
1457 * = 0 : No subdivision.
1458 * = 1 : Subdivision in first parameter
1459 * direction.
1460 * = 2 : Subdivision in second parameter
1461 * direction.
1462 * = 3 : Subdivision in first and second
1463 * parameter direction.
1464 *
1465 * iind1 - Index to first interior knot multiplicity in
1466 * first parameter direction.
1467 * iind2 - Index to first interior knot multiplicity in
1468 * second parameter direction.
1469 * ilevel - If > 0 we subdivide in middlepoint.(SISLSurface only)
1470 * vedge - Pointer to edge maximum.
1471 *
1472 * INPUT/OUTPUT : cmax - The level value.
1473 *
1474 * OUTPUT : pintdat - Maximum data.
1475 * wob[] - Pointers to the subdivided objects.
1476 * jstat - Status messages
1477 * = 2 : New maximum found on new edges.
1478 * = 1 : Maximum found on new edges.
1479 * = 0 : no maximum found.
1480 * < 0 : error
1481 *
1482 *
1483 * METHOD :
1484 *
1485 *
1486 * REFERENCES :
1487 *
1488 *
1489 * WRITTEN BY : Ulf J. Krystad, SI, 89-05.
1490 * REVISED BY : Paal Fugelli, SINTEF, Oslo, 94-07. Fixed memory leaks.
1491 *
1492 *********************************************************************
1493 */
1494 {
1495
1496 SISLPtedge *qpt;
1497
1498 int ki, kj,i; /* Counters. */
1499 int kpos = 0; /* Position of error. */
1500 int kstat= 0; /* Local status */
1501 int kk1, kk2, kn1, kn2; /* Local number of knots and vertices. */
1502 int kmax, kind1, kind2; /* Indexes of the maximum vertice. */
1503 double tstart, tend; /* Start,end and middle values for curve parameter.*/
1504 double sstart[2], send[2],tmidle;/* Start, end values for surface parameter*/
1505 double tpar, tparold; /* The parameter vallue for subdivision curve.*/
1506 double spar[2],sparold[2]; /* The parameter vallue for subdivision
1507 of a surface. */
1508 double *tmax, *tmin;/* Local max and min value for the vertices of object.*/
1509 double sdiff[2]; /* The length of parameter intervall for surface. */
1510 double smin[2]; /* The lower allowed limit in the prameter intervall
1511 for subdividing a surface. */
1512 double smax[2]; /* The upper allowed limit in the prameter intervall
1513 for subdividing a surface. */
1514 SISLSurf *qs1=SISL_NULL; /* Help pointers while subdividing */
1515 SISLSurf *qs2=SISL_NULL; /* Help pointers while subdividing */
1516 SISLObject *qop = SISL_NULL;/* Help pointers while subdividing */
1517 SISLObject *qoc = SISL_NULL;/* Help pointers while subdividing */
1518 SISLIntdat *qintdat=SISL_NULL;/* Local max data for the new edges. */
1519
1520
1521
1522 /* Init */
1523 *jstat = 0;
1524 if ((qop = newObject(SISLPOINT)) == SISL_NULL) goto err101;
1525
1526 if (po1 == SISL_NULL || po1->iobj == SISLPOINT)
1527 /* Nothing to do. */
1528 ;
1529 else if (po1->iobj == SISLCURVE)
1530 {
1531 kk1 = po1->c1->ik;
1532 kn1 = po1->c1->in;
1533 kmax = po1->c1->pbox->imax;
1534 tmax = po1->c1->pbox->emax;
1535 tmin = po1->c1->pbox->emin;
1536
1537 tstart = po1->c1->et[kk1-1];
1538 tend = po1->c1->et[kn1];
1539
1540 /* If we got problems with subdiv in max points, remove as comment: */
1541 /* kmax = 0; */
1542
1543
1544 /* ------------------Determination of sudiv parameter value-----------*/
1545 if (iind1 != 0)
1546 /* We subdivide in an interior knot with multiplicity. */
1547 tpar = po1->c1->et[iind1];
1548
1549 else if (kmax == 0 || kmax == kn1-1)
1550 /*The greatest coeff is the first or last, divide in middlepoint. */
1551 tpar = s1792(po1->c1->et,kk1, kn1);
1552
1553
1554 else
1555 /* Try to find an inner subdivision (ekstremal) point by iteration. */
1556 {
1557
1558 /* First get a good starting point for the iteration. */
1559 tpar = 0;
1560 for (i=kmax+1;i<kmax+kk1;i++)
1561 tpar += po1->c1->et[i];
1562
1563 tpar /= kk1 - 1;
1564 tparold = tpar;
1565
1566 /*Iterate using Newton. */
1567 s1252(po1->c1,aepsge,tpar,&tpar,&kstat);
1568 if (kstat < 0) goto error;
1569
1570 /* Test if the found point is at start or end. */
1571 if(DEQUAL(tpar,tstart) || DEQUAL(tpar,tend))
1572 /*Try Schoenbergs approximation to max vertice. */
1573 {
1574 tpar = tparold;
1575
1576 if(DEQUAL(tpar,tstart) || DEQUAL(tpar,tend))
1577 /*Divide in middlepoint. */
1578 tpar = s1792(po1->c1->et,kk1,kn1);
1579 }
1580 }
1581 /* ------------------Subdivision -------------------------------- */
1582
1583 /* Subdivide the curve at the given parameter value. */
1584 s1231(po1->c1,tpar,&(wob[0]->c1),&(wob[1]->c1),&kstat);
1585 if (kstat < 0) goto error;
1586
1587
1588 /* Pick out end point from a curve. */
1589 s1438(wob[0]->c1,1,&(qop->p1),&tpar,&kstat);
1590 if (kstat < 0) goto error;
1591
1592
1593 /* Examin if the subdividing point is a max. */
1594 s1161(qop,cmax,aepsge,&qintdat,&kstat);
1595 if (kstat < 0) goto error;
1596
1597 if (kstat == 2)
1598 /* New maximum found, delete old ones */
1599 if (*pintdat != SISL_NULL)
1600 {
1601
1602 freeIntdat(*pintdat);
1603 *pintdat = SISL_NULL;
1604 }
1605
1606 if (kstat)
1607 {
1608 /* Maximum found, add them to the list */
1609
1610 *jstat = max(*jstat,kstat); /* Mark maximum found. */
1611
1612 /* Put maximum found on edges into pintdat. */
1613
1614 /* Set parameter border values of object. */
1615 s6idput(pintdat,qintdat,0,tpar,&kstat);
1616 if (kstat < 0) goto error;
1617
1618 if (qintdat != SISL_NULL)
1619 {
1620 freeIntdat(qintdat);
1621 qintdat = SISL_NULL;
1622 }
1623 }
1624 }
1625 else if (po1->iobj == SISLSURFACE)
1626 {
1627 kk1 = po1->s1->ik1;
1628 kn1 = po1->s1->in1;
1629 kk2 = po1->s1->ik2;
1630 kn2 = po1->s1->in2;
1631 kmax = po1->s1->pbox->imax;
1632 tmax = po1->s1->pbox->emax;
1633 tmin = po1->s1->pbox->emin;
1634
1635 sstart[0] = po1->s1->et1[kk1-1];
1636 sstart[1] = po1->s1->et2[kk2-1];
1637
1638 send[0] = po1->s1->et1[kn1];
1639 send[1] = po1->s1->et2[kn2];
1640
1641 sdiff[0] = send[0] - sstart[0];
1642 sdiff[1] = send[1] - sstart[1];
1643 smin[0] = sstart[0] + (double)0.01*sdiff[0];
1644 smin[1] = sstart[1] + (double)0.01*sdiff[1];
1645 smax[0] = send[0] - (double)0.01*sdiff[0];
1646 smax[1] = send[1] - (double)0.01*sdiff[0];
1647
1648 kind2 = kmax/kn1;
1649 kind1 = kmax - kind2*kn1;
1650
1651
1652 /* If we got problems with subdiv in max points, remove as comment: */
1653 /* kind1 = 0; */
1654
1655 /* ------------------Determination of sudiv parameter value-------*/
1656 if (iind1 != 0 || iind2 != 0 || ilevel > 0)
1657 {
1658 if (ilevel > 0)
1659 /* We are forced to subdivide in middlepoint. */
1660 {
1661 spar[0] = s1792(po1->s1->et1,kk1, kn1);
1662 spar[1] = s1792(po1->s1->et2,kk2, kn2);
1663 }
1664
1665 else
1666 /*We have knot multiplicity at least in one parameter direction.
1667 Subdivide in interior knot multiplicity. If the other parameter
1668 direction is without multiplicities, subdivide in middlepoint.*/
1669 {
1670 if (iind1 != 0)
1671 spar[0] = po1->s1->et1[iind1];
1672 else
1673 spar[0] = s1792(po1->s1->et1,kk1, kn1);
1674
1675 if (iind2 != 0 )
1676 spar[1] = po1->s1->et2[iind2];
1677 else
1678 spar[1] = s1792(po1->s1->et2,kk2, kn2);
1679 }
1680 }
1681
1682
1683 else if (kind1 == 0 || kind1 == kn1-1 || kind2 == 0 || kind2 == kn2-1)
1684 {
1685
1686 /*The greatest coeff is on the edge.
1687 Examin the edge for max and divide
1688 in these parameter values. If more than one max,
1689 use the one closest to the middlepoint*/
1690
1691 tmidle = s1792(po1->s1->et1,kk1, kn1);
1692 spar[0] = sstart[0];
1693
1694 for (kj=0;kj<3;kj+=2)
1695 {
1696 qpt = vedge[0]->prpt[kj];
1697 while (qpt != SISL_NULL)
1698 {
1699 if (fabs(qpt->ppt->epar[0] - tmidle) <
1700 fabs(spar[0] - tmidle))
1701 spar[0] = qpt->ppt->epar[0];
1702 qpt = qpt->pnext;
1703 }
1704 }
1705 if (DEQUAL(spar[0],sstart[0]) || DEQUAL(spar[0],send[0]))
1706 spar[0] = tmidle;
1707
1708 tmidle = s1792(po1->s1->et2,kk2, kn2);
1709 spar[1] = sstart[1];
1710
1711 for (kj=1;kj<4;kj+=2)
1712 {
1713 qpt = vedge[0]->prpt[kj];
1714 while (qpt != SISL_NULL)
1715 {
1716 if (fabs(qpt->ppt->epar[0] - tmidle) <
1717 fabs(spar[1] - tmidle))
1718 spar[1] = qpt->ppt->epar[1];
1719 qpt = qpt->pnext;
1720 }
1721 }
1722 if (DEQUAL(spar[1],sstart[1]) || DEQUAL(spar[1],send[1]))
1723 spar[1] = tmidle;
1724 }
1725
1726
1727 else
1728 /* Try to find an inner subdivision (ekstremal) point by iteration. */
1729 {
1730
1731 /* First get a good starting point for the iteration. */
1732 spar[0] = 0;
1733 for (i=kind1+1;i<kind1+kk1;i++)
1734 spar[0] += po1->s1->et1[i];
1735
1736 spar[0] /= kk1 - 1;
1737 sparold[0] = spar[0];
1738
1739 spar[1] = 0;
1740 for (i=kind2+1;i<kind2+kk2;i++)
1741 spar[1] += po1->s1->et2[i];
1742
1743 spar[1] /= kk2 - 1;
1744 sparold[1] = spar[1];
1745
1746
1747 /* Create a point greater than the surface */
1748 if ((qop->p1 = newPoint(tmax,1,1)) == SISL_NULL) goto err101;
1749
1750 /* Iterate using Newton. */
1751 s1173(qop->p1,po1->o1->s1,aepsge,sstart,send,spar,spar,&kstat);
1752 freePoint(qop->p1);
1753 qop->p1 = SISL_NULL;
1754 if (kstat < 0) goto error;
1755
1756 /* Test if the found point is near one edge. */
1757 if(spar[0] < smin[0] ||spar[0] > smax[0]
1758 || spar[1] < smin[1] ||spar[1] > smax[1])
1759 {
1760 /*Try Schoenberg. */
1761 spar[0] = sparold[0];
1762 spar[1] = sparold[1];
1763
1764 if(spar[0] < smin[0] ||spar[0] > smax[0]
1765 || spar[1] < smin[1] ||spar[1] > smax[1])
1766 {
1767 /*Divide in middlepoint. */
1768 spar[0] = s1792(po1->s1->et1,kk1,kn1);
1769 spar[1] = s1792(po1->s1->et2,kk2,kn2);
1770 }
1771 }
1772
1773
1774 }
1775
1776
1777 /* ------------------Subdivision ------------------------------*/
1778 /* Now we have found the parameters for subdivision, divide! */
1779
1780 if ((qoc = newObject(SISLCURVE)) == SISL_NULL)
1781 goto err101;
1782
1783 for (ki=0; ki<(idiv<3 ? 1:3); ki++)
1784 {
1785
1786 if (idiv == 1)
1787 {
1788 s1711(po1->s1,1,spar[0],&(wob[0]->s1),&(wob[1]->s1),&kstat);
1789 if (kstat < 0) goto error;
1790
1791 /* Pick out edge curve from a surface. */
1792
1793 s1435(wob[0]->s1,1,&(qoc->c1),spar,&kstat);
1794 if (kstat < 0) goto error;
1795 }
1796 else if (idiv == 2)
1797 {
1798 s1711(po1->s1,2,spar[1],&(wob[0]->s1),&(wob[1]->s1),&kstat);
1799 if (kstat < 0) goto error;
1800
1801 /* Pick out edge curve from a surface. */
1802
1803 s1435(wob[0]->s1,2,&(qoc->c1),spar+1,&kstat);
1804 if (kstat < 0) goto error;
1805 }
1806 else if (ki == 0)
1807 {
1808 s1711(po1->s1,1,spar[0],&qs1,&qs2,&kstat);
1809 if (kstat < 0) goto error;
1810
1811 /* Pick out edge curve from a surface. */
1812
1813 s1435(qs1,1,&(qoc->c1),spar,&kstat);
1814 if (kstat < 0) goto error;
1815 }
1816 else if (ki == 1)
1817 {
1818 s1711(qs1,2,spar[1],&(wob[0]->s1),&(wob[1]->s1),&kstat);
1819 if (kstat < 0) goto error;
1820
1821 /* Pick out edge curve from a surface. */
1822
1823 s1435(wob[0]->s1,2,&(qoc->c1),spar+1,&kstat);
1824 if (kstat < 0) goto error;
1825 }
1826 else /* if (ki == 2) */
1827 {
1828 s1711(qs2,2,spar[1],&(wob[2]->s1),&(wob[3]->s1),&kstat);
1829 if (kstat < 0) goto error;
1830
1831 /* Pick out edge curve from a surface. */
1832
1833 s1435(wob[2]->s1,2,&(qoc->c1),spar+1,&kstat);
1834 if (kstat < 0) goto error;
1835 }
1836
1837
1838 /* Examine the new edge for max. */
1839
1840 s1161(qoc, cmax, aepsge, &qintdat, &kstat);
1841 if (kstat < 0) goto error;
1842
1843 freeCurve(qoc->c1);
1844 qoc->c1 = SISL_NULL;
1845
1846
1847 if (kstat == 2)
1848 /* New maximum found, delete old ones */
1849 if (*pintdat != SISL_NULL)
1850 {
1851 freeIntdat(*pintdat);
1852 *pintdat = SISL_NULL;
1853 }
1854
1855 if (kstat)
1856 {
1857 /* Maximum found, add them to the list */
1858
1859 *jstat = max(kstat,*jstat); /* Mark maximum found. */
1860
1861 /* Put maximum found on edges into pintdat. */
1862
1863 /* Test if we can pick the second subdivision parameter
1864 from a max on subdiv curve.*/
1865 if(ki==0 && qintdat->vpoint[0]->epar[0] > smin[1]
1866 && qintdat->vpoint[0]->epar[0] < smax[1] )
1867 spar[1]=qintdat->vpoint[0]->epar[0];
1868
1869
1870 /* Set parameter border values of object. */
1871 s6idput(pintdat,qintdat,(ki==0 ? 0:1),spar[(ki==0 ? 0:1)],&kstat);
1872 if (kstat < 0) goto error;
1873
1874 if (qintdat != SISL_NULL)
1875 {
1876 freeIntdat(qintdat);
1877 qintdat = SISL_NULL;
1878 }
1879
1880 }
1881
1882 /* End of for (ki=/..............) */
1883 }
1884
1885 }
1886 goto out;
1887
1888 /* -------------------ERROR SECTION------------------------------------*/
1889
1890 /* Error in space allocation. */
1891 err101: *jstat = -101;
1892 s6err("s1162_s9div",*jstat,kpos);
1893 goto out;
1894
1895 /* Error in lower level routine. */
1896 error : *jstat = kstat;
1897 s6err("s1162_s9div",*jstat,kpos);
1898 goto out;
1899 /* -------------------END OF ERROR SECTION----------------------------*/
1900
1901 out:
1902 if (qop != SISL_NULL) freeObject(qop);
1903 if (qoc != SISL_NULL) freeObject(qoc);
1904 if (qs1 != SISL_NULL) freeSurf(qs1); /* PFU 15/07-94 */
1905 if (qs2 != SISL_NULL) freeSurf(qs2); /* PFU 15/07-94 */
1906 }
1907