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: sh1992.c,v 1.4 2007-08-06 13:09:13 vsk Exp $
45  *
46  */
47 
48 
49 #define SH1992
50 
51 #include "sislP.h"
52 
53 /*
54 * Forward declarations.
55 * ---------------------
56 */
57 
58 #if defined(SISLNEEDPROTOTYPES)
59 static void sh1992_s9mbox3(double [],int,int,double,double,double [],double []);
60 static void sh1992_s9mbox2(double [],int,int,double,double,double [],double []);
61 static void sh1992_s9mbox(double [],int,int,int,double,double,double[],double [],int *);
62 #else
63 static void sh1992_s9mbox3();
64 static void sh1992_s9mbox2();
65 static void sh1992_s9mbox();
66 #endif
67 
68 #if defined(SISLNEEDPROTOTYPES)
sh1992(SISLObject * po,int itype,double aepsge,int * jstat)69 void sh1992(SISLObject *po,int itype,double aepsge,int *jstat)
70 #else
71 void sh1992(po,itype,aepsge,jstat)
72      SISLObject *po;
73      int    itype;
74      double aepsge;
75      int    *jstat;
76 #endif
77 /*
78 *********************************************************************
79 *
80 *********************************************************************
81 *
82 * PURPOSE    : Make a SISLbox on the control-polygons given by
83 *              ecoef to the object. If dimension is 2 then one
84 *	       SISLbox rotated 45 degree is also made. If dimension
85 *	       is 3 then one SISLbox roteted 45 degree around each
86 *	       main axes, in all 4 boxes is made.
87 *
88 *
89 *
90 * INPUT      : po     - Object to treat.
91 *              itype  - Kind of box to make.
92 *                       = 0 : Do not expand box.
93 *                       = 1 : Make a totally expanded box.
94 *                       = 2 : Make a box expanded in the inner of the
95 *                             object, and reduced along the edges/endpoints.
96 *                       If itype>=10, it is interpreted as itype-10, exept
97 *                       for that no rotation of the box is performed.
98 *              aepsge - Geometry resolution.
99 *
100 * OUTPUT     : jstat  - status messages
101 *                                         > 0      : warning
102 *                                         = 0      : ok
103 *                                         < 0      : error
104 *
105 *
106 *
107 * METHOD     :
108 *
109 *
110 * REFERENCES :
111 *
112 *-
113 * CALLS      :  s6existbox - Check if the wanted box exist already.
114 *               s6newbox   - Make a box of a given type.
115 *
116 * WRITTEN BY : Arne Laksaa, SI, 89-02.
117 * MODIFIED BY : Vibeke Skytt, SI, 91-01. Tolerance dependant boxes.
118 * MODIFIED BY : Vibeke Skytt, SI, 92-10. Points are connected
119 *                                        non-expanded boxes. Exept for
120 *                                        1D, boxes are expanded with
121 *                                        eps/2 and eps.
122 *
123 *********************************************************************
124 */
125 {
126    int kstat = 0;                       /* Status variable.        */
127    int kdim;                       /* Dimension of geometry space. */
128    int ktype = itype % 10;              /* Kind of box.            */
129    int knum;                            /* Number of sides of box. */
130    int k2;                              /* Other box type.         */
131    int kbez = 0;                        /* Indicates if Bezier case. */
132    double teps_inner;     /* Tolerance with which to expand in the inner. */
133    double teps_edge;      /* Tolerance with which to expand at the edge.  */
134 
135    /* Set correct tolerances.  */
136 
137    teps_inner = (ktype == 0) ? DZERO : (double)0.5*aepsge;
138    teps_edge = (ktype == 2) ? -teps_inner : teps_inner;
139 
140    if (po -> iobj == SISLPOINT)
141    {
142       if (po->p1->pbox == SISL_NULL)
143 	 if ((po->p1->pbox = newbox(po->p1->idim)) == SISL_NULL) goto err101;
144 
145       if (s6existbox(po->p1->pbox,ktype,aepsge) < 1)
146       {
147      	 kdim = po->p1->idim;
148 	 if (itype < 10 && kdim == 3) knum = 9;
149 	 else if (itype < 10 && kdim == 2) knum = 4;
150 	 else knum = kdim;
151 
152 	 /* The box do not exist already. For a point we always
153 	    use non-expanded boxes.  */
154 
155 	 /* Create the box.  */
156 
157 	 s6newbox(po->p1->pbox,knum,ktype,aepsge,&kstat);
158 	 if (kstat < 0) goto error;
159 
160 	 teps_inner = teps_edge = DZERO;
161 
162 	 k2 = (ktype == 0) ? 0 : ((ktype == 1) ? 2 : 1);
163 	 if (ktype > 0 && s6existbox(po->p1->pbox,k2,aepsge))
164 	    {
165 	       memcopy(po->p1->pbox->e2min[ktype],po->p1->pbox->e2min[k2],
166 		       (1+(kdim!=1))*knum,double);
167 	       memcopy(po->p1->pbox->e2max[ktype],po->p1->pbox->e2max[k2],
168 		       (1+(kdim!=1))*knum,double);
169 	    }
170 	    else
171 	    {
172 	       /* Make the requested box. */
173 
174 	       if (knum == 9)
175 		  sh1992_s9mbox3(po->p1->ecoef,1,1,teps_inner,teps_edge,
176 			  po->p1->pbox->e2max[ktype],po->p1->pbox->e2min[ktype]);
177 	       else if (knum == 4)
178 		  sh1992_s9mbox2(po->p1->ecoef,1,1,teps_inner,teps_edge,
179 			  po->p1->pbox->e2max[ktype],po->p1->pbox->e2min[ktype]);
180 	       else
181 	       {
182 		  sh1992_s9mbox(po->p1->ecoef,1,1,kdim,teps_inner,teps_edge,
183 			 po->p1->pbox->e2max[ktype],po->p1->pbox->e2min[ktype],
184 			 &kstat);
185 		  if (kstat < 0) goto error;
186 	       }
187 	    }
188       }
189    }
190    else if (po -> iobj == SISLCURVE)
191    {
192       if (po->c1->pbox == SISL_NULL)
193 	 if ((po->c1->pbox = newbox(po->c1->idim)) == SISL_NULL) goto err101;
194 
195       if (s6existbox(po->c1->pbox,ktype,aepsge) < 1)
196       {
197      	 kdim = po->c1->idim;
198 	 if (itype < 10 && kdim == 3) knum = 9;
199 	 else if (itype < 10 && kdim == 2) knum = 4;
200 	 else knum = kdim;
201 
202 	 /* The box do not exist already. In the Bezier case,
203 	    it is not necessary to expand in the inner of the curve.  */
204 
205 	 /* Create the box.  */
206 
207 	 s6newbox(po->c1->pbox,knum,ktype,aepsge,&kstat);
208 	 if (kstat < 0) goto error;
209 
210 	 if (kdim != 1 && po->c1->ik == po->c1->in)
211          {
212             teps_inner = DZERO;
213             kbez = 1;
214 	    }
215 
216 	 /* Make the requested box. First allocate scratch for
217 	    box arrays.  */
218 
219 	 if (knum == 9)
220 	    sh1992_s9mbox3(po->c1->ecoef,po->c1->in,1,teps_inner,teps_edge,
221 		    po->c1->pbox->e2max[ktype],po->c1->pbox->e2min[ktype]);
222 	 else if (knum == 4)
223 	    sh1992_s9mbox2(po->c1->ecoef,po->c1->in,1,teps_inner,teps_edge,
224 		    po->c1->pbox->e2max[ktype],po->c1->pbox->e2min[ktype]);
225 	 else
226 	 {
227 	    sh1992_s9mbox(po->c1->ecoef,po->c1->in,1,kdim,teps_inner,
228 		   teps_edge,po->c1->pbox->e2max[ktype],
229 		   po->c1->pbox->e2min[ktype],&kstat);
230 	    if (kstat < 0) goto error;
231          }
232       }
233    }
234    else if (po -> iobj == SISLSURFACE)
235    {
236       if (po->s1->pbox == SISL_NULL)
237 	 if ((po->s1->pbox = newbox(po->s1->idim)) == SISL_NULL) goto err101;
238 
239       if (s6existbox(po->s1->pbox,ktype,aepsge) < 1)
240       {
241      	 kdim = po->s1->idim;
242 	 if (itype < 10 && kdim == 3) knum = 9;
243 	 else if (itype < 10 && kdim == 2) knum = 4;
244 	 else knum = kdim;
245 
246 	 /* The box do not exist already. In the Bezier case, it
247 	    is not necessary to expand in the inner of the surface.  */
248 
249 	 /* Create the box.  */
250 
251 	 s6newbox(po->s1->pbox,knum,ktype,aepsge,&kstat);
252 	 if (kstat < 0) goto error;
253 
254 	 if (kdim != 1 &&
255 	     po->s1->ik1 == po->s1->in1 && po->s1->ik2 == po->s1->in2)
256          {
257 	    teps_inner = DZERO;
258             kbez = 1;
259 	 }
260 
261 	 /* Make the requested box. First allocate scratch for
262 	    box arrays.  */
263 
264 	 if (knum == 9)
265 	    sh1992_s9mbox3(po->s1->ecoef,po->s1->in1,po->s1->in2,teps_inner,
266 		    teps_edge,po->s1->pbox->e2max[ktype],
267 		    po->s1->pbox->e2min[ktype]);
268 	 else if (knum == 4)
269 	    sh1992_s9mbox2(po->s1->ecoef,po->s1->in1,po->s1->in2,teps_inner,
270 		    teps_edge,po->s1->pbox->e2max[ktype],
271 		    po->s1->pbox->e2min[ktype]);
272 	 else
273 	 {
274 	    sh1992_s9mbox(po->s1->ecoef,po->s1->in1,po->s1->in2,kdim,
275 		   teps_inner,teps_edge,po->s1->pbox->e2max[ktype],
276 		   po->s1->pbox->e2min[ktype],&kstat);
277 	    if (kstat < 0) goto error;
278 	 }
279       }
280    }
281 
282   *jstat = kbez;
283   goto out;
284 
285   /* Error in space allocation.  */
286 
287   err101 : *jstat = -101;
288   goto out;
289 
290   /* Error in lower level routine.  */
291 
292   error : *jstat = kstat;
293   goto out;
294 
295  out:
296     return;
297 }
298 
299 #if defined(SISLNEEDPROTOTYPES)
sh1992cu(SISLCurve * pc,int itype,double aepsge,int * jstat)300 void sh1992cu(SISLCurve *pc,int itype,double aepsge,int *jstat)
301 #else
302 void sh1992cu(pc,itype,aepsge,jstat)
303      SISLCurve *pc;
304      int   itype;
305      double aepsge;
306      int   *jstat;
307 #endif
308 /*
309 *********************************************************************
310 *
311 *********************************************************************
312 *
313 * PURPOSE    : Make a SISLbox on the control-polygons given by
314 *              ecoef to the object. If dimension is 2 then one
315 *	       SISLbox rotated 45 degree is also made. If dimension
316 *	       is 3 then one SISLbox roteted 45 degree around each
317 *	       main axes, in all 4 boxes is made.
318 *
319 *
320 *
321 * INPUT      : pc     - SISLObject to treat.
322 *              itype  - Kind of box to make.
323 *                       = 0 : Do not expand box.
324 *                       = 1 : Make a totally expanded box.
325 *                       = 2 : Make a box expanded in the inner of the
326 *                             object, and reduced along the edges/endpoints.
327 *                       If itype>=10, it is interpreted as itype-10, exept
328 *                       for that no rotation of the box is performed.
329 *              aepsge - Geometry resolution.
330 *
331 *
332 * OUTPUT     : jstat  - status messages
333 *                                         > 0      : warning
334 *                                         = 0      : ok
335 *                                         < 0      : error
336 *
337 *
338 *
339 * METHOD     :
340 *
341 *
342 * REFERENCES :
343 *
344 *-
345 * CALLS      :  s6existbox - Check if the wanted box exist already.
346 *               s6newbox   - Make a box of a given type.
347 *
348 *
349 * WRITTEN BY : Arne Laksaa, SI, 89-02.
350 * MODIFIED BY : Vibeke Skytt, SI, 91-01. Tolerance dependant boxes.
351 *
352 *********************************************************************
353 */
354 {
355    int kstat = 0;                       /* Status variable.        */
356    int kdim = pc->idim;                 /* Dimension of geometry space. */
357    int ktype = itype % 10;              /* Kind of box.            */
358    int knum;                            /* Number of sides of box. */
359    int kbez = 0;                        /* Indicates if Bezier case. */
360    double teps_inner;     /* Tolerance with which to expand in the inner. */
361    double teps_edge;      /* Tolerance with which to expand at the edge.  */
362 
363    /* Set number of box sides.  */
364 
365    if (itype < 10 && kdim == 3) knum = 9;
366    else if (itype < 10 && kdim == 2) knum = 4;
367    else knum = kdim;
368 
369    /* Set correct tolerances.  */
370 
371    teps_inner = (ktype == 0) ? DZERO : (double)0.5*aepsge;
372    teps_edge = (ktype == 2) ? -teps_inner : teps_inner;
373 
374    if (pc->pbox == SISL_NULL)
375       if ((pc->pbox = newbox(pc->idim)) == SISL_NULL) goto err101;
376 
377    if (s6existbox(pc->pbox,ktype,aepsge) < 1)
378    {
379       /* The box do not exist already. In the Bezier case,
380 	 it is not necessary to expand in the inner of the curve.  */
381 
382       /* Create the box.  */
383 
384       s6newbox(pc->pbox,knum,ktype,aepsge,&kstat);
385       if (kstat < 0) goto error;
386 
387       if (pc->ik == pc->in)
388       {
389           teps_inner = DZERO;
390           kbez = 1;
391       }
392 
393       /* Make the requested box. First allocate scratch for
394 	 box arrays.  */
395 
396       if (knum == 9)
397 	 sh1992_s9mbox3(pc->ecoef,pc->in,1,teps_inner,teps_edge,
398 		 pc->pbox->e2max[ktype],pc->pbox->e2min[ktype]);
399       else if (knum == 4)
400 	 sh1992_s9mbox2(pc->ecoef,pc->in,1,teps_inner,teps_edge,
401 		 pc->pbox->e2max[ktype],pc->pbox->e2min[ktype]);
402       else
403       {
404 	 sh1992_s9mbox(pc->ecoef,pc->in,1,kdim,teps_inner,teps_edge,
405 		pc->pbox->e2max[ktype],pc->pbox->e2min[ktype],&kstat);
406 	 if (kstat < 0) goto error;
407        }
408    }
409 
410   *jstat = kbez;
411   goto out;
412 
413   /* Error in space allocation.  */
414 
415   err101 : *jstat = -101;
416   goto out;
417 
418   /* Error in lower level routine.  */
419 
420   error : *jstat = kstat;
421   goto out;
422 
423  out:
424     return;
425 }
426 
427 #if defined(SISLNEEDPROTOTYPES)
sh1992su(SISLSurf * ps,int itype,double aepsge,int * jstat)428 void sh1992su(SISLSurf *ps,int itype,double aepsge,int *jstat)
429 #else
430 void sh1992su(ps,itype,aepsge,jstat)
431      SISLSurf *ps;
432      int  itype;
433      double aepsge;
434      int  *jstat;
435 #endif
436 /*
437 *********************************************************************
438 *
439 *********************************************************************
440 *
441 * PURPOSE    : Make a SISLbox on the control-polygons given by
442 *              ecoef to the object. If dimension is 2 then one
443 *	       SISLbox rotated 45 degree is also made. If dimension
444 *	       is 3 then one SISLbox roteted 45 degree around each
445 *	       main axes, in all 4 boxes is made.
446 *
447 *
448 *
449 * INPUT      : ps     - SISLObject to treat.
450 *              itype  - Kind of box to make.
451 *                       = 0 : Do not expand box.
452 *                       = 1 : Make a totally expanded box.
453 *                       = 2 : Make a box expanded in the inner of the
454 *                             object, and reduced along the edges/endpoints.
455 *                       If itype>=10, it is interpreted as itype-10, exept
456 *                       for that no rotation of the box is performed.
457 *              aepsge - Geometry resolution.
458 *
459 *
460 * OUTPUT     : jstat  - status messages
461 *                                         > 0      : warning
462 *                                         = 0      : ok
463 *                                         < 0      : error
464 *
465 *
466 *
467 * METHOD     :
468 *
469 *
470 * REFERENCES :
471 *
472 *-
473 * CALLS      :  s6existbox - Check if the wanted box exist already.
474 *               s6newbox   - Make a box of a given type.
475 *
476 *
477 * WRITTEN BY : Arne Laksaa, SI, 89-02.
478 * MODIFIED BY : Vibeke Skytt, SI, 91-01. Tolerance dependant boxes.
479 *
480 *********************************************************************
481 */
482 {
483    int kstat = 0;                       /* Status variable.        */
484    int kdim = ps->idim;                 /* Dimension of geometry space. */
485    int ktype = itype % 10;              /* Kind of box.            */
486    int knum;                            /* Number of sides of box. */
487    int kbez = 0;                        /* Indicates if Bezier.    */
488    double teps_inner;     /* Tolerance with which to expand in the inner. */
489    double teps_edge;      /* Tolerance with which to expand at the edge.  */
490 
491    /* Set correct tolerances.  */
492 
493    teps_inner = (ktype == 0) ? DZERO : (double)0.5*aepsge;
494    teps_edge = (ktype == 2) ? -teps_inner : teps_inner;
495 
496    /* Set number of box sides.  */
497 
498    if (itype < 10 && kdim == 3) knum = 9;
499    else if (itype < 10 && kdim == 2) knum = 4;
500    else knum = kdim;
501 
502    if (ps->pbox == SISL_NULL)
503       if ((ps->pbox = newbox(ps->idim)) == SISL_NULL) goto err101;
504 
505    if (s6existbox(ps->pbox,ktype,aepsge) < 1)
506    {
507       /* The box do not exist already. In the Bezier case, it
508 	 is not necessary to expand in the inner of the surface.  */
509 
510       /* Create the box.  */
511 
512       s6newbox(ps->pbox,knum,ktype,aepsge,&kstat);
513       if (kstat < 0) goto error;
514 
515       if (ps->ik1 == ps->in1 && ps->ik2 == ps->in2)
516       {
517 	 teps_inner = DZERO;
518          kbez = 1;
519       }
520 
521       /* Make the requested box. First allocate scratch for
522 	 box arrays.  */
523 
524       if (knum == 9)
525 	 sh1992_s9mbox3(ps->ecoef,ps->in1,ps->in2,teps_inner,teps_edge,
526 		 ps->pbox->e2max[ktype],ps->pbox->e2min[ktype]);
527       else if (knum == 4)
528 	 sh1992_s9mbox2(ps->ecoef,ps->in1,ps->in2,teps_inner,teps_edge,
529 		 ps->pbox->e2max[ktype],ps->pbox->e2min[ktype]);
530       else
531       {
532 	 sh1992_s9mbox(ps->ecoef,ps->in1,ps->in2,kdim,
533 		teps_inner,teps_edge,ps->pbox->e2max[ktype],
534 		ps->pbox->e2min[ktype],&kstat);
535 	 if (kstat < 0) goto error;
536       }
537    }
538 
539   *jstat = 0;
540   goto out;
541 
542   /* Error in space allocation.  */
543 
544   err101 : *jstat = -101;
545   goto out;
546 
547   /* Error in lower level routine.  */
548 
549   error : *jstat = kstat;
550   goto out;
551 
552  out:
553     return;
554 }
555 
556 #if defined(SISLNEEDPROTOTYPES)
557 static void
sh1992_s9mbox3(double ecoef[],int icoef1,int icoef2,double aeps1,double aeps2,double e2max[],double e2min[])558   sh1992_s9mbox3(double ecoef[],int icoef1,int icoef2,double aeps1,
559 		 double aeps2,double e2max[],double e2min[])
560 #else
561 static void sh1992_s9mbox3(ecoef,icoef1,icoef2,aeps1,aeps2,e2max,e2min)
562      double ecoef[];
563      int    icoef1;
564      int    icoef2;
565      double aeps1;
566      double aeps2;
567      double e2max[];
568      double e2min[];
569 #endif
570 /*
571 *********************************************************************
572 *
573 *********************************************************************
574 *
575 * PURPOSE    : Make 4 boxes on the control-polygons given by
576 *              ecoef to the object.
577 *
578 *
579 * INPUT      : ecoef   - Control-polygon.
580 *              icoef1  - Number of vertices in control-polygon in
581 *                        first parameter direction.
582 *              icoef2  - Number of vertices in control-polygon in
583 *                        second parameter direction.
584 *              aeps1   - Size of expansion to use in the inner of
585 *                        the object.
586 *              aeps2   - Size of expansion to use at the edge of the object.
587 *
588 *
589 *
590 * OUTPUT     : e2max   - Array to contain maximum values of the box.
591 *	       e2min   - Array to contain minimum values of the box.
592 *
593 *
594 * METHOD     : Make 4 boxes. One ordinary and tree boxes rotated
595 *	       45 degree around the main axes.
596 *
597 *
598 * REFERENCES :
599 *
600 *-
601 * CALLS      :
602 *
603 * WRITTEN BY : Arne Laksaa, SI, 89-02.
604 * MODIFIED BY : Vibeke Skytt, SI, 91-01.
605 *
606 *********************************************************************
607 */
608 {
609   int ki,kj;             /* Counters.                                 */
610   int kant = 9;          /* Number of box sides.                      */
611   int kinset = 0;        /* Indicates if inner box is set.            */
612   double teps1 = aeps1+aeps1; /* Double tolerance in the inner.       */
613   double teps2;               /* Tolerance at edge.                   */
614   double teps3;               /* Double tolerance at edge.            */
615   double t1,t2,t3,t4;    /* To store elements of the rotation matrix. */
616   double *tmin,*tmax;    /* Pointers used to traverse e2min and e2max.  */
617   double sminin[9],smaxin[9];   /* Box boundaries in the inner.       */
618   double sminedg[9],smaxedg[9]; /* Box boundaries at the edge.        */
619   double *sc1,*sc2,*sc3; /* Pointers to coefficients.                 */
620 
621   /* Set tolerances at edge. If the tolerance is positive or dimension
622      is 1D, the input tolerance is used, otherwise we must make sure
623      that the maximum distance from the total box at the edges to the
624      reduced box is aeps2.                                             */
625 
626   if (aeps2 >= DZERO)
627      teps2 = aeps2;
628   else
629      teps2 = (double)0.2767326953*aeps2;
630   teps3 = teps2 + teps2;
631 
632   /* Initiate box boundaries of inner box.  */
633 
634   for (ki=0; ki<kant; ki++)
635   {
636      sminin[ki] = HUGE;
637      smaxin[ki] = -HUGE;
638   }
639 
640   /* Fetch value of first vertex.  */
641 
642   sc1 = ecoef; sc2 = sc1+1;  sc3 = sc2 + 1;
643   t1= ROTM * sc1[0];
644   t2= ROTM * sc2[0];
645   t3= ROTM * sc3[0];
646 
647   tmin = sminedg;
648   tmax = smaxedg;
649   *tmin = *tmax = *sc1;
650   tmin++; tmax++;
651   *tmin = *tmax = *sc2;
652   tmin++; tmax++;
653   *tmin = *tmax = *sc3;
654   tmin++; tmax++;
655   *tmin = *tmax = t2-t3;
656   tmin++; tmax++;
657   *tmin = *tmax = t2+t3;
658   tmin++; tmax++;
659   *tmin = *tmax = t1-t3;
660   tmin++; tmax++;
661   *tmin = *tmax = t1+t3;
662   tmin++; tmax++;
663   *tmin = *tmax = t1-t2;
664   tmin++; tmax++;
665   *tmin = *tmax = t1+t2;
666 
667   /* For each vertice at the edge check and corrigate the box.  */
668 
669   for (ki=0,sc1+=3,sc2+=3,sc3+=3; ki<icoef2; ki++)
670      for (kj=(ki==0); kj<icoef1; kj++,sc1+=3,sc2+=3,sc3+=3)
671      {
672 
673 	/* Set correct pointers.  */
674 
675 	if (((ki==0 || ki==icoef2-1) && icoef2>1) ||
676 		  ((kj==0 || kj==icoef1-1) && icoef1>1))
677 	   tmin = sminedg, tmax = smaxedg;
678 	else
679 	   kinset = 1,  tmin = sminin,  tmax = smaxin;
680 
681 	t1= ROTM * sc1[0];
682 	t2= ROTM * sc2[0];
683 	t3= ROTM * sc3[0];
684 
685 	if(*sc1 < *tmin) *tmin = *sc1;
686 	if(*sc1 > *tmax) *tmax = *sc1;
687 	tmin++; tmax++;
688 	if(*sc2 < *tmin) *tmin = *sc2;
689 	if(*sc2 > *tmax) *tmax = *sc2;
690 	tmin++; tmax++;
691 	if(*sc3 < *tmin) *tmin = *sc3;
692 	if(*sc3 > *tmax) *tmax = *sc3;
693 	tmin++; tmax++;
694 	t4= t2 - t3;
695 	if(t4 < *tmin) *tmin = t4;
696 	if(t4 > *tmax) *tmax = t4;
697 	tmin++; tmax++;
698 	t4= t2 + t3;
699 	if(t4 < *tmin) *tmin = t4;
700 	if(t4 > *tmax) *tmax = t4;
701 	tmin++; tmax++;
702 	t4= t1 - t3;
703 	if(t4 < *tmin) *tmin = t4;
704 	if(t4 > *tmax) *tmax = t4;
705 	tmin++; tmax++;
706 	t4= t1 + t3;
707 	if(t4 < *tmin) *tmin = t4;
708 	if(t4 > *tmax) *tmax = t4;
709 	tmin++; tmax++;
710 	t4= t1 - t2;
711 	if(t4 < *tmin) *tmin = t4;
712 	if(t4 > *tmax) *tmax = t4;
713 	tmin++; tmax++;
714 	t4= t1 + t2;
715 	if(t4 < *tmin) *tmin = t4;
716 	if(t4 > *tmax) *tmax = t4;
717      }
718 
719   /* Merge the inner and the outer box, and adjust with the
720      tolerance.  */
721 
722   if (!kinset)
723   {
724      memcopy(sminin,sminedg,kant,DOUBLE);
725      memcopy(smaxin,smaxedg,kant,DOUBLE);
726   }
727   for (ki=0; ki<kant; ki++)
728   {
729      e2min[ki] = MIN(sminin[ki]-aeps1,sminedg[ki]-teps2);
730      e2max[ki] = MAX(smaxin[ki]+aeps1,smaxedg[ki]+teps2);
731      e2min[kant+ki] = MIN(sminin[ki]-teps1,sminedg[ki]-teps3);
732      e2max[kant+ki] = MAX(smaxin[ki]+teps1,smaxedg[ki]+teps3);
733   }
734 }
735 
736 #if defined(SISLNEEDPROTOTYPES)
737 static void
sh1992_s9mbox2(double ecoef[],int icoef1,int icoef2,double aeps1,double aeps2,double e2max[],double e2min[])738   sh1992_s9mbox2(double ecoef[],int icoef1,int icoef2,double aeps1,
739 		 double aeps2,double e2max[],double e2min[])
740 #else
741 static void sh1992_s9mbox2(ecoef,icoef1,icoef2,aeps1,aeps2,e2max,e2min)
742      double ecoef[];
743      int    icoef1;
744      int    icoef2;
745      double aeps1;
746      double aeps2;
747      double e2max[];
748      double e2min[];
749 #endif
750 /*
751 *********************************************************************
752 *
753 *********************************************************************
754 *
755 * PURPOSE    : Make 2 boxes on the control-polygons given by
756 *              ecoef to the object.
757 *
758 *
759 * INPUT      : ecoef  - Control-polygon.
760 *              icoef1  - Number of vertices in control-polygon in
761 *                        first parameter direction.
762 *              icoef2  - Number of vertices in control-polygon in
763 *                        second parameter direction.
764 *              aeps1   - Size of expansion to use in the inner of
765 *                        the object.
766 *              aeps2   - Size of expansion to use at the edge of the object.
767 *
768 *
769 *
770 * OUTPUT     : e2max   - Array to contain maximum values of the box.
771 *	       e2min   - Array to contain minimum values of the box.
772 *
773 *
774 * METHOD     : Make two boxes. One ordinary and one rotated 45 degree.
775 *
776 *
777 * REFERENCES :
778 *
779 *-
780 * CALLS      :
781 *
782 * WRITTEN BY : Arne Laksaa, SI, 89-02.
783 * MODIFIED BY : Vibeke Skytt, SI, 91-01.
784 *
785 *********************************************************************
786 */
787 {
788   int ki,kj;             /* Counters.                                 */
789   int kant = 4;          /* Number of box sides.                      */
790   int kinset = 0;        /* Indicates if an inner box is found.       */
791   double teps1 = aeps1+aeps1; /* Double tolerance in the inner.       */
792   double teps2;               /* Tolerance at edge.                   */
793   double teps3;               /* Double tolerance at edge.            */
794   double t1,t2,t3;       /* To store elements of the rotation matrix. */
795   double *tmin,*tmax;    /* Pointers used to traverse e2min and e2max.  */
796   double sminin[4],smaxin[4];   /* Box boundaries in the inner.       */
797   double sminedg[4],smaxedg[4]; /* Box boundaries at the edge.        */
798   double *sc1,*sc2;      /* Pointers into coefficient array.          */
799 
800   /* Set tolerances at edge. If the tolerance is positive or dimension
801      is 1D, the input tolerance is used, otherwise we must make sure
802      that the maximum distance from the total box at the edges to the
803      reduced box is aeps2.                                             */
804 
805   if (aeps2 >= DZERO)
806      teps2 = aeps2;
807   else
808      teps2 = (double)0.38268343*aeps2;   /* aeps2 * sin(PI/8).   */
809   teps3 = teps2 + teps2;
810 
811   /* Initiate box boundaries of inner box.  */
812 
813   for (ki=0; ki<kant; ki++)
814   {
815      sminin[ki] = HUGE;
816      smaxin[ki] = -HUGE;
817   }
818 
819   /* Fetch value of first vertex.  */
820 
821   sc1 = ecoef;  sc2 = sc1 + 1;
822   t1= ROTM * sc1[0];
823   t2= ROTM * sc2[0];
824 
825   tmin = sminedg;
826   tmax = smaxedg;
827   *tmin = *tmax = *sc1;
828   tmin++; tmax++;
829   *tmin = *tmax = *sc2;
830   tmin++; tmax++;
831   *tmin = *tmax = t1-t2;
832   tmin++; tmax++;
833   *tmin = *tmax = t1+t2;
834 
835   /* For each vertex check and corrigate the box.  */
836 
837   for (ki=0,sc1+=2,sc2+=2; ki<icoef2; ki++)
838      /* UJK, writing error */
839      /*for (kj=(ki==1); kj<icoef1; kj++,sc1+=2,sc2+=2) */
840 
841      for (kj=(ki==0); kj<icoef1; kj++,sc1+=2,sc2+=2)
842      {
843 	/* Set correct box boundaries.  */
844 
845 	if (((ki==0 || ki==icoef2-1) && icoef2>1) ||
846 		  ((kj==0 || kj==icoef1-1) && icoef1>1))
847 	   tmin = sminedg,  tmax = smaxedg;
848 	else
849 	   kinset = 1,  tmin = sminin,  tmax = smaxin;
850 
851 	t1= ROTM * sc1[0];
852 	t2= ROTM * sc2[0];
853 
854 	if(*sc1 < *tmin) *tmin = *sc1;
855 	if(*sc1 > *tmax) *tmax = *sc1;
856 	tmin++; tmax++;
857 	if(*sc2 < *tmin) *tmin = *sc2;
858 	if(*sc2 > *tmax) *tmax = *sc2;
859 	tmin++; tmax++;
860 	t3= t1 - t2;
861 	if(t3 < *tmin) *tmin = t3;
862 	if(t3 > *tmax) *tmax = t3;
863 	tmin++; tmax++;
864 	t3= t1 + t2;
865 	if(t3 < *tmin) *tmin = t3;
866 	if(t3 > *tmax) *tmax = t3;
867      }
868 
869   /* Merge the inner and the outer box, and adjust with the
870      tolerance.  */
871 
872   if (!kinset)
873   {
874      memcopy(sminin,sminedg,kant,DOUBLE);
875      memcopy(smaxin,smaxedg,kant,DOUBLE);
876   }
877   for (ki=0; ki<kant; ki++)
878   {
879      e2min[ki] = MIN(sminin[ki]-aeps1,sminedg[ki]-teps2);
880      e2max[ki] = MAX(smaxin[ki]+aeps1,smaxedg[ki]+teps2);
881      e2min[kant+ki] = MIN(sminin[ki]-teps1,sminedg[ki]-teps3);
882      e2max[kant+ki] = MAX(smaxin[ki]+teps1,smaxedg[ki]+teps3);
883   }
884 }
885 
886 #if defined(SISLNEEDPROTOTYPES)
887 static void
sh1992_s9mbox(double ecoef[],int icoef1,int icoef2,int idim,double aeps1,double aeps2,double e2max[],double e2min[],int * jstat)888   sh1992_s9mbox(double ecoef[],int icoef1,int icoef2,int idim,
889 		double aeps1,double aeps2,double e2max[],
890 		double e2min[],int *jstat)
891 #else
892 static void sh1992_s9mbox(ecoef,icoef1,icoef2,idim,aeps1,aeps2,
893 			  e2max,e2min,jstat)
894      double ecoef[];
895      int    icoef1;
896      int    icoef2;
897      int    idim;
898      double aeps1;
899      double aeps2;
900      double e2max[];
901      double e2min[];
902      int *jstat;
903 #endif
904 /*
905 *********************************************************************
906 *
907 *********************************************************************
908 *
909 * PURPOSE    : Make a SISLbox on the control-polygons given by
910 *              ecoef to the object.
911 *
912 *
913 * INPUT      : ecoef  - Control-polygon.
914 *              icoef1  - Number of vertices in control-polygon in
915 *                        first parameter direction.
916 *              icoef2  - Number of vertices in control-polygon in
917 *                       second parameter direction.
918 *              aeps1   - Size of expansion to use in the inner of
919 *                        the object.
920 *              aeps2   - Size of expansion to use at the edge of the object.
921 *
922 *
923 *
924 * OUTPUT     : e2max   - Array to contain maximum values of the box.
925 *	       e2min   - Array to contain minimum values of the box.
926 *              jstat  - status messages
927 *                                         > 0      : warning
928 *                                         = 0      : ok
929 *                                         < 0      : error
930 *
931 *
932 * METHOD     :
933 *
934 *
935 * REFERENCES :
936 *
937 *-
938 * CALLS      :
939 *
940 * WRITTEN BY : Arne Laksaa, SI, 89-02.
941 * MODIFIED BY : Vibeke Skytt, SI, 91-01.
942 *
943 *********************************************************************
944 */
945 {
946   int ki,ki1,kj;       /* Counters.  */
947   int kant = idim;     /* Number of box sides.                        */
948   int kinset = 0;      /* Indicates if the inner box is set.          */
949   double noice = (double)100.0*REL_COMP_RES;   /* Noice killer.       */
950   double teps1 = aeps1+aeps1; /* Double tolerance in the inner.       */
951   double teps2;               /* Tolerance at edge.                   */
952   double teps3;               /* Double tolerance at edge.            */
953   double *tmin,*tmax;  /* Pointers into box boundary arrays.          */
954   double *sc;          /* Pointer into coefficient array.             */
955   double *sminin=SISL_NULL,*smaxin=SISL_NULL;  /* Box boundaries of the inner.  */
956   double *sminedg=SISL_NULL,*smaxedg=SISL_NULL; /* Box boundaries of the edge.  */
957 
958   /* Set tolerances at edge. If the tolerance is positive or dimension
959      is 1D, the input tolerance is used, otherwise we must make sure
960      that the maximum distance from the total box at the edges to the
961      reduced box is aeps2.                                             */
962 
963   if (idim == 1 || aeps2 >= DZERO)
964      teps2 = aeps2;
965   else
966      teps2 = aeps2/sqrt((double)idim);
967   teps3 = teps2 + teps2;
968 
969   /* Allocate scratch for intermediate box arrays.  */
970 
971   if ((sminin = newarray(kant,double)) == SISL_NULL) goto err101;
972   if ((smaxin = newarray(kant,double)) == SISL_NULL) goto err101;
973   if ((sminedg = newarray(kant,double)) == SISL_NULL) goto err101;
974   if ((smaxedg = newarray(kant,double)) == SISL_NULL) goto err101;
975 
976   /* Initiate box boundaries of inner box.  */
977 
978   for (ki=0; ki<kant; ki++)
979   {
980      sminin[ki] = HUGE;
981      smaxin[ki] = -HUGE;
982   }
983 
984   /* Fetch value of first vertex.  */
985 
986   for (ki = 0; ki < idim; ki++)
987      sminedg[ki] = smaxedg[ki] = ecoef[ki];
988 
989   /* For each vertice check and corrigate the box.  */
990 
991   for (sc=ecoef+idim, ki=0; ki<icoef2; ki++)
992      for (kj=(ki==0); kj<icoef1; kj++)
993      {
994 	/* Set correct box.  */
995 
996 	if (((ki==0 || ki==icoef2-1) && icoef2>1) ||
997 		  ((kj==0 || kj==icoef1-1) && icoef1>1))
998 	   tmin = sminedg,  tmax = smaxedg;
999 	else
1000 	    kinset = 1,  tmin = sminin,  tmax = smaxin;
1001 
1002 	for (ki1=0; ki1<idim; ki1++,sc++,tmin++,tmax++)
1003 	{
1004 	   if(*sc < *tmin) *tmin = *sc;
1005 	   if(*sc > *tmax) *tmax = *sc;
1006 	}
1007      }
1008 
1009   /* Merge the inner and the outer box, and adjust with the
1010      tolerance.  */
1011 
1012   if (!kinset)
1013   {
1014      memcopy(sminin,sminedg,kant,DOUBLE);
1015      memcopy(smaxin,smaxedg,kant,DOUBLE);
1016   }
1017   for (ki=0; ki<kant; ki++)
1018   {
1019      e2min[ki] = MIN(sminin[ki]-aeps1,sminedg[ki]-teps2);
1020      e2max[ki] = MAX(smaxin[ki]+aeps1,smaxedg[ki]+teps2);
1021      if (idim > 1)
1022      {
1023 	e2min[kant+ki] = MIN(sminin[ki]-teps1,sminedg[ki]-teps3);
1024 	e2max[kant+ki] = MAX(smaxin[ki]+teps1,smaxedg[ki]+teps3);
1025      }
1026   }
1027 
1028   /* ALA and UJK 30.10.90, remove noice near by zero.  */
1029 
1030   if (idim == 1)
1031   {
1032      if (fabs(e2max[0]) < noice) e2max[0] = DZERO;
1033      if (fabs(e2min[0]) < noice) e2min[0] = DZERO;
1034   }
1035 
1036   *jstat = 0;
1037   goto out;
1038 
1039   /* Error in scratch allocation. */
1040 
1041   err101 : *jstat = -101;
1042   goto out;
1043 
1044   out :
1045   if (sminin != SISL_NULL) freearray(sminin);
1046   if (smaxin != SISL_NULL) freearray(smaxin);
1047   if (sminedg != SISL_NULL) freearray(sminedg);
1048   if (smaxedg != SISL_NULL) freearray(smaxedg);
1049 }
1050 
1051 
1052