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: s1343.c,v 1.2 2001-03-19 15:58:46 afr Exp $
45 *
46 */
47
48 #define S1343
49
50 #include "sislP.h"
51
52 #if defined(SISLNEEDPROTOTYPES)
s1343(SISLCurve * pc,double eeps[],int ilend,int irend,double aepsco,int itmax,SISLCurve ** rc,int * jstat)53 void s1343(SISLCurve *pc,double eeps[],int ilend,int irend,double aepsco,
54 int itmax, SISLCurve **rc,int *jstat)
55 #else
56 void s1343(pc,eeps,ilend,irend,aepsco,itmax,rc,jstat)
57 SISLCurve *pc;
58 double eeps[];
59 int ilend;
60 int irend;
61 double aepsco;
62 int itmax;
63 SISLCurve **rc;
64 int *jstat;
65 #endif
66 /*
67 ***********************************************************
68 *
69 * Purpose: To approximate the input spline-curve by a cubic spline-
70 * curve with error less than eeps in each of the kdim components.
71 *
72 * Input :
73 * Pc - Pointer to curve.
74 * Eeps - Array (length kdim) giving the desired accuracy of
75 * the spline-approximation in each component.
76 * Ilend - The no. of derivatives that are not allowed to change
77 * at the left end of the curve.
78 * The (0 - (ilend-1)) derivatives will be kept fixed.
79 * If ilend <0, this routine will set it to 0.
80 * If ilend<kord, this routine will set it to kord.
81 * Irend - The no. of derivatives that are not allowed to change
82 * at the right end of the curve.
83 * The (0 - (irend-1)) derivatives will be kept fixed.
84 * If irend <0, this routine will set it to 0.
85 * If irend<kord, this routine will set it to kord.
86 * Aepsco - Two numbers differing by a relative amount <aepsco,
87 * will in some cases be considered equal. A suitable
88 * value is just above the unit roundoff in the machine.
89 * On a VAX 1.0E-6 is a reasonable choice.
90 * The computations are not guaranteed to have relative
91 * accuracy less than aepsco.
92 * Itmax - Max. no. of iterations.
93 *
94 * Output:
95 * Jstat - Output status:
96 * < 0 : Error.
97 * = 0 : Ok.
98 * > 0 : Warning:
99 * Rc - Pointer to curve.
100 *-
101 * Method: First sampling points on the input-curve are determined such
102 * that the cubic hermite spline interpolant to these points will
103 * deviate with less than half the tolerance from the input-spline,
104 * and this hermite spline interpolant is run through data
105 * reduction to reduce the no. of knots.
106 *
107 * Calls: s1379, s1340, s1355, s6err
108 *
109 * The main routine, s1340, is written by Knut M|rken, Si.
110 * Written by: C.R. Birkeland, Si, April 1993.
111 **********************************************************
112 */
113 {
114 int i,j; /* Loop control variables */
115 int index = 0;
116 int idim = pc->idim; /* Space dimension */
117 int stat = 0; /* Error control parameters */
118 int kpos = 0;
119 int km;
120 int leftknot = 0;
121 double *error1 = SISL_NULL;
122 double *error2 = SISL_NULL;
123 double *epar = SISL_NULL;
124 double *derive = SISL_NULL;
125 double *kp = SISL_NULL;
126 double *kder = SISL_NULL;
127 SISLCurve *ocurve = SISL_NULL;
128 SISLCurve *qc_kreg = SISL_NULL; /* Non-periodic version of the input curve. */
129
130
131 /* Check input-curve. */
132
133 if (!pc) goto err150;
134 s1707(pc, &stat);
135 if (stat<0) goto error;
136
137 /* Make sure that the input curve is non-periodic. */
138
139 if (pc->cuopen == SISL_CRV_PERIODIC)
140 {
141 make_cv_kreg(pc, &qc_kreg, &stat);
142 if (stat < 0) goto error;
143
144 /* The input curve is closed and periodic. Make sure that the
145 endpoints of the curve is still matching by fixing the
146 position and the derivative in the endpoints of the curve.
147 The change made to startfix and endfix is only locallly. */
148
149 ilend = MAX(ilend, 2);
150 irend = MAX(irend, 2);
151 }
152 else
153 qc_kreg = pc;
154
155 /* Check input */
156
157 if (qc_kreg->in < qc_kreg->ik || qc_kreg->ik < 1 || idim < 1) goto err103;
158 if (qc_kreg->ik < 5)
159 {
160 /* Increase order of curve to 4 (cubic) */
161
162 s1750(qc_kreg, 4, &ocurve, &stat);
163 if (stat < 0) goto error;
164
165 /* Curve is now a cubic spline
166 * Call reduction routine */
167
168 if( (error2 = newarray( idim, DOUBLE )) == SISL_NULL) goto err101;
169 s1340( ocurve, eeps, ilend, irend, aepsco, itmax, rc,
170 error2, &stat);
171 if (stat < 0) goto error;
172
173 /* Success ! Go out ! */
174
175 goto out;
176 }
177
178 /* Set local tolerance */
179
180 if( (error1 = newarray(idim, DOUBLE)) == SISL_NULL) goto err101;
181
182 for (i=0; i<idim; i++)
183 error1[i] = 0.5*eeps[i];
184
185 /* Determine the sample points */
186
187 s1355( qc_kreg, error1, &epar, &km, &stat );
188 if (stat < 0) goto error;
189
190 /* Compute positions and derivatives of input spline
191 * at the given sampling points */
192
193 derive = newarray( idim * 2, DOUBLE );
194 kp = newarray( idim * km, DOUBLE );
195 kder = newarray( idim * km, DOUBLE );
196 if (derive == SISL_NULL || kp == SISL_NULL || kder == SISL_NULL) goto err101;
197
198 for(i=0; i<km; i++)
199 {
200 if(epar[i] == epar[i+1])
201 s1227( qc_kreg, 1, epar[i], &leftknot, derive, &stat );
202 else
203 s1221( qc_kreg, 1, epar[i], &leftknot, derive, &stat );
204
205 for(j=0; j<idim; j++, index++)
206 {
207 kp[index] = derive[j];
208 kder[index] = derive[j+idim];
209 }
210 }
211 if (stat < 0) goto error;
212
213 /* Compute Hermite interpolant */
214
215 s1379( kp, kder, epar, km, idim, &ocurve, &stat );
216 if (stat < 0) goto error;
217
218 /* Compute datareduction on the cubic hermite interpolant */
219
220 if( (error2 = newarray( idim, DOUBLE )) == SISL_NULL) goto err101;
221 s1340( ocurve, error1, ilend, irend, aepsco, itmax, rc,
222 error2, &stat);
223 if (stat < 0) goto error;
224
225 /* Set periodicity flag. */
226
227 if (pc->cuopen == SISL_CRV_CLOSED ||
228 pc->cuopen == SISL_CRV_PERIODIC)
229 (*rc)->cuopen = SISL_CRV_CLOSED;
230
231 /* Success ! */
232
233 *jstat = 0;
234 goto out;
235
236 /* Allocation error. */
237
238 err101:
239 *jstat = -101;
240 s6err("s1343",*jstat,kpos);
241 goto out;
242
243 /* Error in input */
244
245 err103:
246 *jstat = -103;
247 s6err("s1343",*jstat,kpos);
248 goto out;
249
250 /* Empty curve. */
251
252 err150:
253 *jstat = -150;
254 s6err("s1343",*jstat,kpos);
255 goto out;
256
257 /* Error in lower level routine. */
258
259 error:
260 *jstat = stat;
261 s6err("s1343",*jstat,kpos);
262 goto out;
263
264 /* Exit */
265
266 out:
267 /* Free allocated arrays */
268
269 if( error1 != SISL_NULL) freearray(error1);
270 if( error2 != SISL_NULL) freearray(error2);
271 if( epar != SISL_NULL) freearray(epar);
272 if( derive != SISL_NULL) freearray(derive);
273 if( kp != SISL_NULL) freearray(kp);
274 if( kder != SISL_NULL) freearray(kder);
275
276 /* Free local SISL-curves */
277
278 if( ocurve != SISL_NULL) freeCurve(ocurve);
279 if (qc_kreg != SISL_NULL && qc_kreg != pc) freeCurve(qc_kreg);
280
281 return;
282 }
283