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 S1535
43
44 #include "sislP.h"
45
46 #if defined(SISLNEEDPROTOTYPES)
s1535(double points[],double der10[],double der01[],double der11[],int im1,int im2,int idim,double par1[],double par2[],int con1,int con2,int con3,int con4,int order1,int order2,int iopen1,int iopen2,SISLSurf ** rsurf,int * jstat)47 void s1535(double points[],double der10[],double der01[],double der11[],
48 int im1,int im2,int idim,double par1[],
49 double par2[],int con1,int con2,int con3,int con4,
50 int order1, int order2, int iopen1, int iopen2,
51 SISLSurf **rsurf,int *jstat)
52 #else
53 void s1535(points,der10, der01,der11,im1,im2,idim,par1, par2,
54 con1,con2,con3,con4,order1, order2,iopen1,iopen2,rsurf,jstat)
55 double points[];
56 double der10[];
57 double der01[];
58 double der11[];
59 int im1;
60 int im2;
61 int idim;
62 double par1[];
63 double par2[];
64 int con1;
65 int con2;
66 int con3;
67 int con4;
68 int order1;
69 int order2;
70 int iopen1;
71 int iopen2;
72 SISLSurf **rsurf;
73 int *jstat;
74 #endif
75 /*
76 ************************************************************************
77 *
78 * PURPOSE: To compute a B-spline tensor surface interpolating a set
79 * of points.
80 *
81 * INPUT:
82 * points - Array of dimension idim*im1*im2 containing
83 * the positions of the nodes (using the same ordering
84 * as ecoef in the SISLSurf structure).
85 *
86 * der10 - Array of dimension idim*im1*im2 containing the first
87 * derivatives in the first parameter direction.
88 *
89 * der01 - Array of dimension idim*im1*im2 containing the first
90 * derivatives in the second parameter direction.
91 *
92 * der11 - Array of dimension idim*im1*im2 containing the cross
93 * derivatives (the twists).
94 *
95 * im1 - The number of interpolation points in the
96 * first parameter direction.
97 *
98 * im2 - The number of interpolation points in the
99 * second parameter direction.
100 *
101 * idim - Dimension of the space we are working in.
102 *
103 * par1 - Parametrization in first parameter direction.
104 * For closed curves, one additional parameter value
105 * must be spesified. The last entry contains
106 * the parametrization of the repeted start point.
107 * (if the endpoint is equal to the startpoint of
108 * the interpolation the lenght of the array could
109 * be equal to im1 also in the closed case).
110 *
111 * par2 - Parametrization in second parameter direction.
112 * For closed curves, one additional parameter value
113 * must be spesified. The last entry contains
114 * the parametrization of the repeted start point.
115 * (if the endpoint is equal to the startpoint of
116 * the interpolation the lenght of the array could
117 * be equal to im2 also in the closed case).
118 *
119 *
120 * ^ Second par. direction
121 * |
122 * | (2.)
123 * |-----------|
124 * | |
125 * (3.) | | (4.)
126 * | |
127 * | |
128 * |-----------|-> First par. direction
129 * (1.)
130 *
131 * con1 - Additional condition along edge 1:
132 * = 0: No additional condition.
133 * = 1: Zero curvature.
134 *
135 * con2 - Additional condition along edge 2:
136 * = 0: No additional condition.
137 * = 1: Zero curvature.
138 *
139 * con3 - Additional condition along edge 3:
140 * = 0: No additional condition.
141 * = 1: Zero curvature.
142 *
143 * con4 - Additional condition along edge 4:
144 * = 0: No additional condition.
145 * = 1: Zero curvature.
146 *
147 * order1 - Order of surface in first parameter direction.
148 *
149 * order2 - Order of surface in second parameter direction.
150 *
151 * iopen1 - Open/close parameter in first parameter direction.
152 * = 1 : open surface.
153 * = 0 : closed, non-periodic surface.
154 * = -1 : periodic surface
155 *
156 * iopen2 - Open/close parameter in second parameter direction.
157 * = 1 : open surface.
158 * = 0 : closed, non-periodic surface.
159 * = -1 : periodic surface
160 *
161 *
162 * Output:
163 * rsurf - Pointer to the surf produced
164 * jstat - Status variable
165 * < 0 - Error.
166 *
167 * Method:
168 * The interpolation is accomplished by using a one dimensional
169 * routine for spline interpolation called several times.
170 * First, the datapoints
171 * are considered to be idim*im1 dimentional and so on...
172 *
173 *
174 * REFERENCES :
175 *
176 * CALLS : s1357, s6chpar.
177 *
178 * WRITTEN BY : Christophe Rene Birkeland, SINTEF, June 1993.
179 * CHANGED BY : Vibeke Skytt, SINTEF, 0394. Introduced iopen1 and iopen2.
180 *
181 *********************************************************************
182 */
183 {
184 int i, j, k, len; /* Loop control parameter */
185 int kpek1, kpek2, kpek3;
186 int idimm1, newindim;
187 int maxim; /* Max (im1, im2) */
188 int kstat=0; /* Status variable */
189 int kpos=0; /* Position of error */
190 int newin1, newin2; /* Number of vertices along par. dir. 1 & 2 */
191 int numpt; /* Needed in call to s1357 */
192 double start=0.; /* Needed in call to s1357 */
193 double end;
194 int *typept=SISL_NULL; /* Array needed for call to s1357 */
195 double *pointpar=SISL_NULL; /* Array needed for call to s1357 */
196 double *coeffpos=SISL_NULL; /* Array needed for call to s1357 */
197 double *coeffder=SISL_NULL; /* Array needed for call to s1357 */
198 double *coeffposder=SISL_NULL; /* Array needed for call to s1357 */
199 double *newpoint=SISL_NULL;
200 double *newder=SISL_NULL;
201 SISLCurve *curve1a=SISL_NULL, *curve1b=SISL_NULL;
202 SISLCurve *curve2=SISL_NULL;
203
204 /* Allocate and initialize necessary arrays for call to s1357 */
205
206 maxim = 2 * MAX(im1, im2);
207 if((typept = newarray(maxim, INT))==SISL_NULL) goto err101;
208 for(i=0; i<maxim; i+=2)
209 {
210 typept[i] = 1;
211 typept[i+1] = 4;
212 }
213
214 idimm1 = idim*im1;
215 len = 2*im2*idimm1;
216 if((newpoint = newarray(len, DOUBLE))==SISL_NULL) goto err101;
217 if((newder = newarray(len, DOUBLE))==SISL_NULL) goto err101;
218 for(i=0, kpek1=0, kpek2=0, kpek3 = idimm1;
219 i<im2;
220 i++, kpek1+=(2*idimm1), kpek3+=(2*idimm1), kpek2+=idimm1)
221 {
222 for(k=0; k<idimm1; k++)
223 {
224 newpoint[kpek1+k] = points[kpek2+k];
225 newpoint[kpek3+k] = der01[kpek2+k];
226 newder[kpek1+k] = der10[kpek2+k];
227 newder[kpek3+k] = der11[kpek2+k];
228 }
229 }
230
231 /* INTERPOLATION in SECOND direction : position. */
232
233 s1357(newpoint, im2*2, idimm1, typept, par2, con1, con2, iopen2, order2,
234 start, &end, &curve1a, &pointpar, &numpt, &kstat);
235 if(kstat < 0) goto error;
236 if(pointpar != SISL_NULL)
237 {
238 freearray(pointpar);
239 pointpar = SISL_NULL;
240 }
241 newin2 = curve1a->in;
242
243 /* INTERPOLATION in SECOND direction : derivative. */
244
245 s1357(newder, im2*2, idimm1, typept, par2, con1, con2, iopen2, order2,
246 start, &end, &curve1b, &pointpar, &numpt, &kstat);
247 if(kstat < 0) goto error;
248 if(pointpar != SISL_NULL)
249 {
250 freearray(pointpar);
251 pointpar = SISL_NULL;
252 }
253 if(curve1b->in != newin2) goto err116;
254
255 /* Transpose results, store new coefficients in
256 * arrays coeffpos and coeffder */
257
258 newindim = newin2 * idim;
259 if((coeffpos = newarray(im1*newindim, DOUBLE)) == SISL_NULL)
260 goto err101;
261 if((coeffder = newarray(im1*newindim, DOUBLE)) == SISL_NULL)
262 goto err101;
263 s6chpar(curve1a->ecoef, im1, newin2, idim, coeffpos);
264 s6chpar(curve1b->ecoef, im1, newin2, idim, coeffder);
265
266 if((coeffposder = newarray(2*im1*newindim, DOUBLE)) == SISL_NULL)
267 goto err101;
268 for(j=0, kpek1=0, kpek2=newindim, kpek3=0;
269 j<im1;
270 j++, kpek1+=(2*newindim), kpek2+=(2*newindim), kpek3+=newindim)
271 {
272 for(k=0; k<newindim; k++)
273 {
274 coeffposder[kpek1+k] = coeffpos[kpek3+k];
275 coeffposder[kpek2+k] = coeffder[kpek3+k];
276 }
277 }
278
279 /* Interpolation in FIRST parameter direction */
280
281 s1357(coeffposder, 2*im1, idim*newin2, typept, par1, con3, con4, iopen1,
282 order1, start, &end, &curve2, &pointpar, &numpt, &kstat);
283 if(kstat < 0) goto error;
284 if(pointpar != SISL_NULL)
285 {
286 freearray(pointpar);
287 pointpar = SISL_NULL;
288 }
289 newin1 = curve2->in;
290
291 /* Transpose back coefficients */
292
293 if((coeffposder=increasearray(coeffposder, idim*newin1*newin2, DOUBLE))
294 == SISL_NULL) goto err101;
295 s6chpar(curve2->ecoef, newin2, newin1, idim, coeffposder);
296
297 /* Create instance of surface */
298
299 if (((*rsurf) = newSurf(newin1, newin2, order1, order2, curve2->et,
300 curve1a->et, coeffposder, 1, idim, 1)) == SISL_NULL)
301 goto err101;
302
303 /* Set periodicity flag. */
304
305 (*rsurf)->cuopen_1 = curve2->cuopen;
306 (*rsurf)->cuopen_2 = curve1a->cuopen;
307
308 /* Success */
309
310 *jstat = 0;
311 goto out;
312
313 /* Allocation error. */
314
315 err101:
316 *jstat = -101;
317 s6err("s1535",*jstat,kpos);
318 goto out;
319
320 /* Error. */
321
322 err116:
323 *jstat = -116;
324 s6err("s1535",*jstat,kpos);
325 goto out;
326
327 /* Error in lower level routine. */
328
329 error: *jstat =kstat;
330 s6err("s1535",*jstat,kpos);
331 goto out;
332
333 out:
334 /* Free arrays */
335
336 if (typept != SISL_NULL) freearray(typept);
337 if(newpoint != SISL_NULL) freearray(newpoint);
338 if(newder != SISL_NULL) freearray(newder);
339 if(coeffpos != SISL_NULL) freearray(coeffpos);
340 if(coeffder != SISL_NULL) freearray(coeffder);
341 if(coeffposder != SISL_NULL) freearray(coeffposder);
342
343 /* Free local SISL-curve objects */
344
345 if(curve1a != SISL_NULL) freeCurve(curve1a);
346 if(curve1b != SISL_NULL) freeCurve(curve1b);
347 if(curve2 != SISL_NULL) freeCurve(curve2);
348
349 return;
350 }
351