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 #define S1967
43
44 #include "sislP.h"
45
46 #if defined(SISLNEEDPROTOTYPES)
s1967(double ep[],double etang1[],double etang2[],double eder11[],int im1,int im2,int idim,int ipar,double epar1[],double epar2[],double eeps[],int nend[],int iopen1,int iopen2,double edgeps[],int iopt,int itmax,SISLSurf ** rs,double emxerr[],int * jstat)47 void s1967(double ep[],double etang1[],double etang2[],double eder11[],
48 int im1,int im2,int idim,int ipar,double epar1[],double epar2[],
49 double eeps[],int nend[],int iopen1,int iopen2,double edgeps[],
50 int iopt,int itmax,SISLSurf **rs,double emxerr[],
51 int *jstat)
52 #else
53 void s1967(ep,etang1,etang2,eder11,im1,im2,idim,ipar,epar1,epar2,
54 eeps,nend,iopen1,iopen2,edgeps,iopt,itmax,rs,emxerr,jstat)
55 double ep[];
56 double etang1[];
57 double etang2[];
58 double eder11[];
59 int im1;
60 int im2;
61 int idim;
62 int ipar;
63 double epar1[];
64 double epar2[];
65 double eeps[];
66 int nend[];
67 int iopen1;
68 int iopen2;
69 double edgeps[];
70 int iopt;
71 int itmax;
72 SISLSurf **rs;
73 double emxerr[];
74 int *jstat;
75 #endif
76 /*
77 ********************************************************************
78 *
79 * Purpose: To compute a bicubic hermite spline-approximation to the
80 * position and derivative data given by ep,etang1,etang2
81 * and eder11.
82 *
83 * Input : Ep - Array (length idim*im1*im2) containing the points
84 * to be approximated.
85 * Etang1 - Array (length idim*im1*im2) containing the derivatives
86 * (tangents) in the first parameter-direction at the
87 * data-points.
88 * Etang2 - Array (length idim*im1*im2) containing the derivatives
89 * (tangents) in the second parameter-direction at the
90 * data-points.
91 * Eder11 - Array (length idim*im1*im2) containing the cross (twist)
92 * derivatives at the data-points.
93 * Im1 - The no. of points in the first parameter-direction.
94 * Im2 - The no. of points in the second parameter-direction.
95 * Ipar - Flag determining the parametrization of the data-points.
96 * = 1: Mean accumulated cord length parametrization.
97 * = 2: Uniform parametrization.
98 * = 3: Parametrization given by epar1 and epar2.
99 * Epar1 - Array (length im1) containing a parametrization in the
100 * first parameter-direction. (Will only be used if ipar=3.)
101 * Epar2 - Array (length im2) containing a parametrization in the
102 * surface lies.)
103 * Eeps - Array (length idim) containing the maximum deviation
104 * which is acceptable in each of the idim components of
105 * the surface (except possibly along the edges).
106 * Nend - Array (length 4) containing the no. of derivatives
107 * to be kept fixed along each edge of the surface.
108 * The numbering of the edges is the same as for edeps below.
109 * All the derivatives of order < nend(i)-1 will be kept fixed
110 * along edge no. i. Hence nend(i)=0 indicates that nothing
111 * is to be kpet fixed along edge no. i.
112 * To be kept fixed here means to have error less than edgeps.
113 * In general, it is impossible to remove any knots and
114 * keep an edge completely fixed.
115 * iopen1 - Open/closed parameter in first parameter direction.
116 * = 1 : Produce open surface.
117 * = 0 : Produce closed, non-periodic surface if possible.
118 * = -1 : Produce closed, periodic surface if possible.
119 * NB! The surface will be closed/periodic only if the first
120 * and last column of data are (approximately) equal.
121 * iopen2 - Open/closed parameter in second parameter direction.
122 * = 1 : Produce open surface.
123 * = 0 : Produce closed, non-periodic surface if possible.
124 * = -1 : Produce closed, periodic surface if possible.
125 * NB! The surface will be closed/periodic only if the first
126 * and last row of data are (approximately) equal.
127 * Edgeps - Array (length idim*4) containing the max. deviation
128 * which is acceptable along the edges of the surfaces.
129 * Edgeps(1,i):edgeps(idim,i) gives the tolerance along
130 * the edge corresponding to the i-th parameter having
131 * its min. or max value.
132 * i=1 : min value of first parameter.
133 * i=2 : max value of first parameter.
134 * i=3 : min value of second parameter.
135 * i=4 : max value of second parameter.
136 * Edgeps(kp,i) will only have any significance if nend(i)>0.
137 * Iopt - Flag indicationg the order in which tha data-reduction
138 * is to be performed.
139 * = 1 : Remove knots in parameter-direction 1 only.
140 * = 2 : Remove knots in parameter-direction 2 only.
141 * = 3 : Remove knots in parameter-direction 1 and
142 * and then in parameter-direction 2.
143 * = 4 : Remove knots in parameter-direction 2 and
144 * and then in parameter-direction 1.
145 * Itmax - Max. no. of iteration.
146 *
147 * Ouput: Jstat - Output status:
148 * < 0 : Error.
149 * = 0 : Ok.
150 * > 0 : Warning.
151 * Rs - Pointer to surface.
152 * Emxerr - Array (length idim) (allocated outside this routine.)
153 * containing an upper bound for the error comitted in
154 * each component during the data reduction.
155 *
156 * Method: First the bicubic hermite spline-interpolant is computed
157 * using the appropriate parametrization, and then knots
158 * are removed from this approximation by a call to the
159 * data-reduction routine for surfaces.
160 *
161 * Calls: s1530, s1965, s6err.
162 *
163 * Written by: C.R.Birkeland, Si, Oslo, Norway, May 1993.
164 * Changed and renamed by : Vibeke Skytt, SINTEF Oslo, 02.95. Introduced
165 * periodicity.
166 **********************************************************************
167 */
168 {
169 int stat=0, kpos=0; /* Error message parameters */
170 double *par1 = SISL_NULL; /* Used to store parametrizations */
171 double *par2 = SISL_NULL;
172 SISLSurf *osurf = SISL_NULL; /* Hermite interp. surface */
173
174 /* Check Input */
175
176 if (im1 < 2 || im2 < 2 || idim < 1)
177 goto err103;
178 if (ipar < 1 || ipar > 3) ipar = 1;
179
180 /* Generate parametrization */
181
182 if (ipar != 3)
183 {
184 /* Generate parametrization */
185
186 s1528(idim, im1, im2, ep, ipar, SISL_CRV_OPEN, SISL_CRV_OPEN,
187 &par1, &par2, &stat);
188 if (stat<0) goto error;
189 }
190 else
191 {
192 /* Parametrization is passed as parameter */
193
194 par1 = epar1;
195 par2 = epar2;
196 }
197
198 /* Perform bicubic hermite spline interpolation */
199
200 s1530(ep, etang1, etang2, eder11, par1, par2,
201 im1, im2, idim, &osurf, &stat);
202 if (stat<0) goto error;
203
204 /* Perform final datareduction step */
205
206 s1965(osurf, eeps, nend, iopen1, iopen2, edgeps, iopt, itmax,
207 rs, emxerr, &stat);
208 if (stat<0) goto error;
209
210 /* Success */
211
212 *jstat = 0;
213 goto out;
214
215
216 /* Error in input */
217
218 err103:
219 *jstat = -103;
220 s6err("s1967",*jstat,kpos);
221 goto out;
222
223 /* Error in lower level routine. */
224
225 error:
226 *jstat = stat;
227 s6err("s1967",*jstat,kpos);
228 goto out;
229
230 /* Exit. */
231
232 out:
233 /* Free SISL-surface allocated in this routine */
234
235 if(osurf != SISL_NULL) freeSurf(osurf);
236
237 if (ipar != 3)
238 {
239 if(par1 != SISL_NULL) freearray(par1);
240 if(par2 != SISL_NULL) freearray(par2);
241 }
242 return;
243 }
244