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 #define S1966
44
45 #include "sislP.h"
46
47 #if defined(SISLNEEDPROTOTYPES)
s1966(double ep[],int im1,int im2,int idim,int ipar,double epar1[],double epar2[],double eeps[],int nend[],int iopen1,int iopen2,double edgeps[],double afctol,int iopt,int itmax,int ik1,int ik2,SISLSurf ** rs,double emxerr[],int * jstat)48 void s1966(double ep[],int im1,int im2,int idim,int ipar,double epar1[],
49 double epar2[],double eeps[],int nend[],int iopen1,int iopen2,
50 double edgeps[],double afctol,int iopt,int itmax,
51 int ik1,int ik2,SISLSurf **rs,double emxerr[],int *jstat)
52 #else
53 void s1966(ep,im1,im2,idim,ipar,epar1,epar2,eeps,nend,iopen1,iopen2,
54 edgeps,afctol,iopt,itmax,ik1,ik2,
55 rs,emxerr,jstat)
56 double ep[];
57 int im1;
58 int im2;
59 int idim;
60 int ipar;
61 double epar1[];
62 double epar2[];
63 double eeps[];
64 int nend[];
65 int iopen1;
66 int iopen2;
67 double edgeps[];
68 double afctol;
69 int iopt;
70 int itmax;
71 int ik1;
72 int ik2;
73 SISLSurf **rs;
74 double emxerr[];
75 int *jstat;
76 #endif
77 /*
78 ********************************************************************
79 *
80 * Purpose: To compute a tensor-product spline-approximation of order
81 * (ik1,ik2) to the rectangular array of idim-dimensional
82 * points given by ep.
83 *
84 * Input : Ep - Array (length idim*im1*im2) containing the points
85 * to be approximated.
86 * Im1 - The no. of points in first parameter-direction.
87 * Im2 - The no. of points in second parameter-direction.
88 * Idim - The no. of components of each input-point.
89 * The approximation will be a parametric surface
90 * situated in the idim-dimensional euclidean space
91 * (usually 3).
92 * Ipar - Flag determining the parametrization of the data points:
93 * = 1: Mean accumulated cord-length parameterization.
94 * = 2: Uniform parametrization.
95 * = 3: Parametrization given by epar1 and epar2.
96 * Epar1 - Array (length im1) containing a parametrization
97 * in the first parameter-direction. (Will only
98 * be used if ipar=3).
99 * Epar2 - Array (length im2) containing a parametrization
100 * in the second parameter-direction. (Will only
101 * be used if ipar=3).
102 * Eeps - Array (length idim) containing the max. permissible
103 * deviation of the approximation from the given data
104 * points, in each of the components. More specifically,
105 * the approximation will not deviate more than eeps(kdim)
106 * in component no. kdim, from the bilinear approximation
107 * to the data.
108 * Nend - Array (length 4) giving the no. of derivatives to be
109 * kept fixed along each edge of the bilinear interpolant.
110 * The numbering of the edges is the same as for edgeps below.
111 * All the derivatives of order < (nend(i)-1) will be kept
112 * fixed along the edge i. Hence nend(i)=0 indicates that
113 * nothing is to be kpet fixed along edge i. (Used by the
114 * data reduction routine.)
115 * To be kept fixed here means to have error less than edgeps.
116 * In general, it is impossible to remove any knots and keep
117 * an edge completely fixed.
118 * iopen1 - Open/closed parameter in first parameter direction.
119 * = 1 : Produce open surface.
120 * = 0 : Produce closed, non-periodic surface if possible.
121 * = -1 : Produce closed, periodic surface if possible.
122 * NB! The surface will be closed/periodic only if the first
123 * and last column of data points are (approximately) equal.
124 * iopen2 - Open/closed parameter in second parameter direction.
125 * = 1 : Produce open surface.
126 * = 0 : Produce closed, non-periodic surface if possible.
127 * = -1 : Produce closed, periodic surface if possible.
128 * NB! The surface will be closed/periodic only if the first
129 * and last row of data points are (approximately) equal.
130 * Edgeps - Array (length idim*4) containing the max. deviation from
131 * the bilinear interpolant which is acceptable along the
132 * edges of the surface.
133 * Edgeps(1,i):edgeps(idim,i) gives the tolerance along
134 * the edge corresponding to the i-th parameter having
135 * one of it`s extremal-values.
136 * i=1: min value of first parameter.
137 * i=2: max value of first parameter.
138 * i=3: min value of second parameter.
139 * i=4: max value of second parameter.
140 * (Used by the data-reduction routine.)
141 * Edgeps(kp,i) will only have significance if nend(i)>0.
142 * Afctol - 0.0 >= afctol <= 1.0.
143 * Afctol indicates how the tolerance is to be shared
144 * between the two data-reduction stages. For the linear
145 * reduction, a tolerance of afctol*eeps will be used,
146 * while a tolerance of (1.0-afctol)*eeps will be used
147 * during the final data reduction (similarly for edgeps.)
148 * Default is 0.
149 * Iopt - Flag indicating the order in which the data-reduction
150 * is to be performed:
151 * = 1: Remove knots inparameter-direction 1 only.
152 * = 2: Remove knots inparameter-direction 2 only.
153 * = 3: Remove knots first in parameter-direction 1 and
154 * then in parameter-direction 2.
155 * = 4: Remove knots first in parameter-direction 2 and
156 * then in parameter-direction 1.
157 * Itmax - Max. no. of iterations in the data-reduction.
158 * Ik1 - The order of the approximation in first
159 * parameter-directon.
160 * Ik2 - The order of the approximation in second
161 * parameter-directon.
162 *
163 * Output:
164 * Jstat - Output status:
165 * < 0 : Error.
166 * = 0 : Ok.
167 * > 0 : Warning:
168 * Rs - Pointer to surface.
169 * Emxerr - Array (length idim) (allocated outside this routine.)
170 * containing the error in the approximation to the data.
171 * This is guaranteed upper bound on the max. deviation
172 * in each component, between the final approximation
173 * and the bilinear spline-approximation to the original data.
174 *
175 * Method:
176 * First the bilinear interpolant to the data is computed, using the
177 * parameterization given by ipar, and knots are removed from this
178 * initial approximation by a call to the data-reduction routine for
179 * surfaces. Then the order is raised to (ik1,ik2) and the final data
180 * reduction is performed.
181 *-
182 * Calls: s1965, s1350, s6chpar, s6err.
183 *
184 * Written by: C.R.Birkeland, Si, April 1993.
185 * Changed by: Per OEyvind, SINTEF, 1994-11.
186 * Removed following memory leaks:
187 * 1) Improper use of copy flag to newSurf()
188 * 2) Forgetting to free temp array after using icopy == 1
189 * Changed and renamed by : Vibeke Skytt, SINTEF Oslo, 02.95. Introduced
190 * periodicity.
191 **********************************************************************
192 */
193 {
194 int in1,in2; /* Number of vertices */
195 int newin1, newin2;
196 int fouridim=4*idim;
197 int i; /* Loop control parameters */
198 int stat=0, kpos=0; /* Error message parameters */
199 double *par1 = SISL_NULL;
200 double *par2 = SISL_NULL;
201 double *knot1 = SISL_NULL; /* Knot vectors in 1 and 2. par.dir. */
202 double *knot2 = SISL_NULL;
203 double *error1 = SISL_NULL; /* Arrays for error storage */
204 double *error2 = SISL_NULL;
205 double *maxerr = SISL_NULL;
206 double *newcoeff = SISL_NULL; /* Coefficients array */
207 SISLCurve *ocurve1 = SISL_NULL; /* Used to store local curves */
208 SISLCurve *ocurve2 = SISL_NULL;
209 SISLSurf *osurf1 = SISL_NULL; /* Used to store local surfaces */
210 SISLSurf *osurf2 = SISL_NULL;
211
212 /* Check Input */
213
214 if (im1 < 2 || im2 < 2 || ik1 < 1 || ik2 < 1 || idim < 1)
215 goto err103;
216 if (ipar < 1 || ipar > 3) ipar = 1;
217
218 if (ipar != 3)
219 {
220 /* Generate parametrization */
221
222 s1528(idim, im1, im2, ep, ipar, SISL_CRV_OPEN, SISL_CRV_OPEN,
223 &par1, &par2, &stat);
224 if (stat<0) goto error;
225 }
226 else
227 {
228 /* Parametrization is passed as parameter */
229
230 par1 = epar1;
231 par2 = epar2;
232 }
233
234 /* Represent input (points) as a surface of
235 * order 2 (linear) in both directions.
236 * First, generate knot vectors */
237
238 knot1 = newarray(im1+2, DOUBLE);
239 knot2 = newarray(im2+2, DOUBLE);
240 if(knot1 == SISL_NULL || knot2 == SISL_NULL) goto err101;
241 memcopy(&knot1[1],par1,im1,DOUBLE);
242 memcopy(&knot2[1],par2,im2,DOUBLE);
243 knot1[0] = knot1[1];
244 knot2[0] = knot2[1];
245 knot1[im1+1] = knot1[im1];
246 knot2[im2+1] = knot2[im2];
247 osurf1 = newSurf(im1, im2, 2, 2, knot1, knot2, ep,
248 1,idim, 1);
249 if (osurf1 == SISL_NULL) goto err101;
250 if (knot1 != SISL_NULL) freearray(knot1); knot1 = SISL_NULL;
251 if (knot2 != SISL_NULL) freearray(knot2); knot2 = SISL_NULL;
252
253 /* Compute tolerance vectors for linear reduction
254 * Both max deviation of surface and max dev. of edges */
255
256 maxerr = newarray(idim, DOUBLE);
257 error1 = newarray(idim, DOUBLE);
258 error2 = newarray(fouridim, DOUBLE);
259 if (error1 == SISL_NULL || error2 == SISL_NULL || maxerr == SISL_NULL)
260 goto err101;
261 for (i=0; i<fouridim; i++)
262 {
263 edgeps[i] = MIN(edgeps[i], eeps[(i+idim)%idim]);
264 error2[i] = afctol * edgeps[i];
265 }
266 for (i=0; i<idim; i++)
267 error1[i] = afctol*eeps[i];
268
269 /* Perform datareduction on the bilinear interpolant */
270
271 s1965(osurf1, error1, nend, SISL_CRV_OPEN, SISL_CRV_OPEN, error2,
272 iopt, itmax, &osurf2, maxerr, &stat);
273 if (stat<0) goto error;
274
275 in1 = osurf2->in1;
276 in2 = osurf2->in2;
277
278 /* Free surface osurf1 */
279
280 if(osurf1 != SISL_NULL)
281 {
282 freeSurf(osurf1);
283 osurf1 = SISL_NULL;
284 }
285
286 /* Piecewise linear interpolant to the reduced
287 * bilinear interpolant expressed as a surface
288 * of orders ik1 and ik2 */
289
290 /* Second parameter direction */
291
292 s1350(osurf2->ecoef,&(osurf2->et2)[1], in2,
293 in1 * idim, ik2, &ocurve1, &stat);
294 if (stat<0) goto error;
295
296 newin2 = ocurve1->in;
297
298 /* Transpose result, store new coefficients in
299 * array newcoeff */
300
301 if( (newcoeff = newarray(idim * in1 * newin2, DOUBLE)) == SISL_NULL )
302 goto err101;
303 s6chpar(ocurve1->ecoef, in1, newin2, idim, newcoeff);
304
305 /* First parameter direction */
306
307 s1350(newcoeff, &(osurf2->et1)[1], in1,
308 idim*newin2, ik1, &ocurve2, &stat);
309 if (stat<0) goto error;
310 newin1 = ocurve2->in;
311
312 /* Free surface osurf2 */
313
314 if(osurf2 != SISL_NULL)
315 {
316 freeSurf(osurf2);
317 osurf2 = SISL_NULL;
318 }
319
320 /* Transpose back and get coefficients of bilinear
321 * approximatoin surface of orders ik1 and ik2 */
322
323 newcoeff = increasearray(newcoeff,
324 idim * newin1 * newin2, DOUBLE);
325 if (newcoeff == SISL_NULL) goto err101;
326 s6chpar(ocurve2->ecoef, newin2, newin1, idim, newcoeff);
327
328 /* Store results as a surface */
329
330 osurf1 = newSurf(newin1, newin2, ik1, ik2, ocurve2->et,
331 ocurve1->et, newcoeff, 1, idim, 1);
332
333 if (newcoeff != SISL_NULL) freearray(newcoeff); newcoeff = SISL_NULL;
334
335 if (osurf1 == SISL_NULL) goto err101;
336
337 /* Set periodicity flag. */
338
339 osurf1->cuopen_1 = ocurve2->cuopen;
340 osurf1->cuopen_2 = ocurve1->cuopen;
341
342 /* Compute tolerance for final datareduction */
343
344 for (i=0; i<fouridim; i++)
345 error2[i] = edgeps[i]-error2[i];
346 for (i=0; i<idim; i++)
347 error1[i] = eeps[i]-maxerr[i];
348
349 /* Perform final datareduction step */
350
351 s1965(osurf1, error1, nend, iopen1, iopen2, error2, iopt, itmax,
352 rs, emxerr, &stat);
353 if (stat<0) goto error;
354
355 /* Compute total (and final) error */
356
357 for(i=0; i<idim; i++)
358 emxerr[i] += maxerr[i];
359
360 /* Success */
361
362 *jstat = 0;
363 goto out;
364
365 /* Empty array. */
366
367 err101:
368 *jstat = -101;
369 s6err("s1966",*jstat,kpos);
370 goto out;
371
372 /* Error in input */
373
374 err103:
375 *jstat = -103;
376 s6err("s1966",*jstat,kpos);
377 goto out;
378
379 /* Error in lower level routine. */
380
381 error:
382 *jstat = stat;
383 s6err("s1966",*jstat,kpos);
384 goto out;
385
386 /* Exit */
387
388 out:
389 /* Free SISL-curves allocated in this routine */
390
391 if(ocurve1 != SISL_NULL) freeCurve(ocurve1);
392 if(ocurve2 != SISL_NULL) freeCurve(ocurve2);
393
394 /* Free SISL-surfaces allocated in this routine */
395
396 if(osurf1 != SISL_NULL) freeSurf(osurf1);
397
398 /* Free arrays */
399
400 if(error1 != SISL_NULL) freearray(error1);
401 if(error2 != SISL_NULL) freearray(error2);
402 if(maxerr != SISL_NULL) freearray(maxerr);
403
404 if (ipar != 3)
405 {
406 freearray(par1);
407 freearray(par2);
408 }
409
410 return;
411 }
412