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: s1920.c,v 1.3 2001-03-19 15:58:56 afr Exp $
45  *
46  */
47 
48 
49 #define S1920
50 
51 #include "sislP.h"
52 
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1920(SISLCurve * pc1,double edir[],int idim,double aepsco,double aepsge,int * jpt,double ** gpar,int * jcrv,SISLIntcurve *** wcurve,int * jstat)55 s1920(SISLCurve *pc1,double edir[],int idim,double aepsco,double aepsge,
56 	   int *jpt,double **gpar,int *jcrv,SISLIntcurve ***wcurve,int *jstat)
57 #else
58 void s1920(pc1,edir,idim,aepsco,aepsge,jpt,gpar,jcrv,wcurve,jstat)
59      SISLCurve    *pc1;
60      double   edir[];
61      int      idim;
62      double   aepsco;
63      double   aepsge;
64      int      *jpt;
65      double   **gpar;
66      int      *jcrv;
67      SISLIntcurve ***wcurve;
68      int      *jstat;
69 #endif
70 /*
71 *********************************************************************
72 *
73 *********************************************************************
74 *
75 * PURPOSE    : Find the extremal points/intervals of the curve pc1 in
76 *              the direction edir.
77 *
78 *
79 *
80 * INPUT      : pc1    - Pointer to the curve.
81 *              edir   - The direction in which the extremal point(s)
82 *                       and/or interval(s) are to be calculated. If
83 *                       idim=1 a positive value indicates the maximum
84 *                       of the B-spline function and a negative value
85 *                       the minimum. If the dimension is greater that
86 *                       1 the array contains the coordinates of the
87 *                       direction vector.
88 *              idim   - Dimension of the space in which the vector edir
89 *                       lies.
90 *              aepsco - Computational resolution.
91 *              aepsge - Geometry resolution.
92 *
93 *
94 *
95 * OUTPUT     : jpt    - Number of single extremal points.
96 *              gpar   - Array containing the parameter values of the
97 *                       single extremal points in the parameter
98 *                       interval of the curve. The points lie continuous.
99 *                       Extremal curves are stored in wcurve.
100 *              jcrv   - Number of extremal curves.
101 *              wcurve - Array containing descriptions of the extremal
102 *                       curves. The curves are only described by points
103 *                       in the parameter interval. The curve-pointers points
104 *                       to nothing. (See descrjption of Intcurve
105 *                       in intcurve.dcl).
106 *              jstat  - status messages
107 *                                         > 0      : warning
108 *                                         = 0      : ok
109 *                                         < 0      : error
110 *
111 *
112 * METHOD     : The scalar-product of the coefficients of the curve with
113 *              the vector edir is calculated giving a curve with dimension
114 *              1. Then the maxima of this curve is calculated.
115 *
116 *
117 *
118 * REFERENCES :
119 *
120 *-
121 * CALLS      : s1880 - Put extremal points/intervals into output format.
122 *              s1161 - Find maxima of one-dimensional curve.
123 *              s6scpr - Scalar-product between two vectors.
124 *              newCurve   - Create new curve.
125 *              newObject  - Create new object.
126 *              newIntpt   - Create new extremal point.
127 *              freeObject - Free the space occupied by a given object.
128 *              freeIntdat - Free space occupied by an extremal list.
129 *
130 * WRITTEN BY : Vibeke Skytt, SI, 88-10.
131 * REVISED BY : Mike Floater, SI, 1991-04 for a rational curve.
132 * Revised by : Paal Fugelli, SINTEF, Oslo, Norway, 1994-08.  Updated to handle
133 *              input curve with periodic basis correctly.
134 *********************************************************************
135 */
136 {
137   int ikind;               /* Type of curve pc1 is                     */
138   int kstat = 0;           /* Local status variable.                     */
139   int kpos = 0;            /* Position of error.                         */
140   int ki;                  /* Counter.                                   */
141   int kn;                  /* Number of vertices of curve.               */
142   int kk;                  /* Order of curve.                            */
143   double tmax;             /* Estimate of maximal value of 1-dim. curve. */
144   double *st;              /* Pointer to knotvector of curve.            */
145   double *scoef;           /* Pointer to vertices of curve.              */
146   double *sc = SISL_NULL;       /* Pointer to vertices of curve in maxima
147 			      calculation.                               */
148   double *spar = SISL_NULL;     /* Values of extremal points in the parameter
149 			      area of the second object. Empty in this case.*/
150   double *s1,*s2,*sstop;   /* Pointers used to traverse double-arrays.   */
151   SISLIntdat *qintdat = SISL_NULL;  /* Maximum results.                     */
152   SISLCurve *qc = SISL_NULL;        /* Pointer to curve in maxima calculation.    */
153   SISLObject *qo1 = SISL_NULL;      /* Pointer to object in maxima calculation. */
154   SISLCurve *qkreg = SISL_NULL;     /* Input curve with k-regularity ensured. */
155 
156 
157   *jpt = 0;
158   *jcrv = 0;
159 
160 
161   /* Ensure k-regular basis. */
162 
163   if ( pc1 -> cuopen == SISL_CRV_PERIODIC )
164   {
165     /* Cyclic (periodic) basis. */
166 
167     make_cv_kreg(pc1, &qkreg, &kstat);
168     if ( kstat < 0 )  goto error;
169   }
170   else
171     qkreg = pc1;
172 
173 
174   /* Check dimension.  */
175 
176   if ( qkreg -> idim != idim )  goto err106;
177 
178   /* Describe curve with local variables.  */
179 
180   kn = qkreg -> in;
181   kk = qkreg -> ik;
182   st = qkreg -> et;
183   ikind = qkreg -> ikind;
184 
185   if ( ikind == 2 || ikind == 4 )
186   {
187       scoef = qkreg -> rcoef;
188       /* Allocate space for coeffecients of new curve.  */
189 
190       if ( (sc = newarray(2*kn, DOUBLE)) == SISL_NULL )  goto err101;
191 
192       /* Compute scalar-product of curve-vertices and direction vector. */
193       /* Copy over weights. */
194 
195       for ( s1=scoef, s2=sc, sstop=s2+2*kn;  s2<sstop;  s1+=idim+1, s2+=2 )
196       {
197           *s2 = s6scpr(s1, edir, idim);
198           *(s2+1) = *(s1+idim);
199       }
200   }
201   else
202   {
203       scoef = qkreg -> ecoef;
204       /* Allocate space for coeffecients of new curve.  */
205 
206       if ( (sc = newarray(kn, DOUBLE)) == SISL_NULL )  goto err101;
207 
208       /* Compute scalar-product of curve-vertices and direction vector. */
209 
210       for ( s1=scoef, s2=sc, sstop=s2+kn;  s2<sstop;  s1+=idim, s2++ )
211       {
212           *s2 = s6scpr(s1, edir, idim);
213       }
214   }
215 
216 
217   /* Create new curve.  */
218 
219   qc = newCurve(kn, kk, st, sc, qkreg->ikind, 1, 1);
220   if ( qc == SISL_NULL )  goto err101;
221 
222   /* Create new object and connect curve to object.  */
223 
224   qo1 = newObject(SISLCURVE);
225   if ( qo1 == SISL_NULL )  goto err101;
226   qo1 -> c1 = qc;
227 
228   tmax = -(double)HUGE;
229 
230   /* Find maxima. */
231   s1161(qo1, &tmax, aepsge, &qintdat, &kstat);
232   if ( kstat < 0 )  goto error;
233 
234   if (qintdat)
235   {
236 
237     /* Express maximal points/intervals on output format.  */
238     s1880(1, 0, &qintdat->ipoint, qintdat->vpoint,
239 	  &qintdat->ilist, qintdat->vlist,
240 	  jpt, gpar, &spar, jcrv, wcurve, &kstat);
241     if ( kstat < 0 )  goto error;
242 
243     /* Handle periodicity (remove points at periodic curve ends). */
244     if ( *jpt > 1 && idim > 1 && pc1->cuopen == SISL_CRV_PERIODIC )
245     {
246       for ( ki=0; ki < (*jpt); ki++ )
247       {
248 	if ( (*gpar)[ki] == pc1->et[pc1->in] )
249 	{
250 	  (*jpt)--;
251 	  (*gpar)[ki] = (*gpar)[*jpt];
252 	  ki--;
253 	}
254       }
255     }
256   }
257 
258   /* Extremal points/intervals found.  */
259 
260   *jstat = 0;
261   goto out;
262 
263   /* Error in space allocation.  */
264 
265  err101:
266   *jstat = -101;
267   s6err("s1920", *jstat, kpos);
268   goto out;
269 
270   /* Dimensions conflicting.  */
271 
272  err106:
273   *jstat = -106;
274   s6err("s1920", *jstat, kpos);
275   goto out;
276 
277   /* Error in lower level routine.  */
278 
279  error:
280   *jstat = kstat;
281   s6err("s1920", *jstat, kpos);
282   goto out;
283 
284  out:
285 
286   /* Free allocated space.  */
287 
288   if ( qkreg && qkreg != pc1 ) freeCurve(qkreg);
289   if (sc) freearray(sc);
290   if (spar) freearray(spar);
291   if (qintdat) freeIntdat(qintdat);
292   if (qo1) freeObject(qo1);
293 
294 
295   return;
296 }
297