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 S1537
43
44 #include "sislP.h"
45
46 #if defined(SISLNEEDPROTOTYPES)
s1537(double points[],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 s1537(double points[],int im1,int im2,int idim,double par1[],
48 double par2[],int con1,int con2,int con3,int con4,
49 int order1, int order2,int iopen1,int iopen2,
50 SISLSurf **rsurf,int *jstat)
51 #else
52 void s1537(points,im1,im2,idim,par1, par2,con1,con2,con3,
53 con4,order1, order2,iopen1,iopen2,rsurf,jstat)
54 double points[];
55 int im1;
56 int im2;
57 int idim;
58 double par1[];
59 double par2[];
60 int con1;
61 int con2;
62 int con3;
63 int con4;
64 int order1;
65 int order2;
66 int iopen1;
67 int iopen2;
68 SISLSurf **rsurf;
69 int *jstat;
70 #endif
71 /*
72 ************************************************************************
73 *
74 * PURPOSE: To compute a B-spline tensor surface interpolating a set
75 * of points.
76 *
77 * INPUT:
78 * points - Array of dimension idim*im1*im2 containing
79 * the positions of the nodes (using the same ordering
80 * as ecoef in the SISLSurf structure).
81 *
82 * im1 - The number of interpolation points in the
83 * first parameter direction.
84 *
85 * im2 - The number of interpolation points in the
86 * second parameter direction.
87 *
88 * idim - Dimension of the space we are working in.
89 *
90 * par1 - Parametrization in first parameter direction.
91 * For closed curves, one additional parameter value
92 * must be spesified. The last entry contains
93 * the parametrization of the repeted start point.
94 * (if the endpoint is equal to the startpoint of
95 * the interpolation the lenght of the array could
96 * be equal to im1 also in the closed case).
97 *
98 * par2 - Parametrization in second parameter direction.
99 * For closed curves, one additional parameter value
100 * must be spesified. The last entry contains
101 * the parametrization of the repeted start point.
102 * (if the endpoint is equal to the startpoint of
103 * the interpolation the lenght of the array could
104 * be equal to im2 also in the closed case).
105 *
106 *
107 * ^ Second par. direction
108 * |
109 * | (2.)
110 * |-----------|
111 * | |
112 * (3.) | | (4.)
113 * | |
114 * | |
115 * |-----------|-> First par. direction
116 * (1.)
117 *
118 * con1 - Additional condition along edge 1:
119 * = 0: No additional condition.
120 * = 1: Zero curvature.
121 *
122 * con2 - Additional condition along edge 2:
123 * = 0: No additional condition.
124 * = 1: Zero curvature.
125 *
126 * con3 - Additional condition along edge 3:
127 * = 0: No additional condition.
128 * = 1: Zero curvature.
129 *
130 * con4 - Additional condition along edge 4:
131 * = 0: No additional condition.
132 * = 1: Zero curvature.
133 *
134 * order1 - Order of surface in first parameter direction.
135 *
136 * order2 - Order of surface in second parameter direction.
137 *
138 * iopen1 - Open/close parameter in first parameter direction.
139 * = 1 : open surface.
140 * = 0 : closed, non-periodic surface.
141 * = -1 : periodic surface.
142 *
143 * iopen2 - Open/close parameter in second parameter direction.
144 * = 1 : open surface.
145 * = 0 : closed, non-periodic surface.
146 * = -1 : periodic surface.
147 *
148 *
149 * Output:
150 * rsurf - Pointer to the surf produced
151 * jstat - Status variable
152 * < 0 - Error.
153 *
154 * Method:
155 * The interpolation is accomplished by using a one dimensional
156 * routine for spline interpolation called several times.
157 * First, the datapoints
158 * are considered to be idim*im1 dimentional and so on...
159 *
160 *
161 * REFERENCES :
162 *
163 * CALLS : s1357
164 *
165 * WRITTEN BY : Christophe Rene Birkeland, SINTEF, May 1993.
166 * REWISED BY : Vibeke Skytt, SINTEF, 0394. Introduced iopen1, iopen2.
167 *
168 *********************************************************************
169 */
170 {
171 int i; /* Loop control parameter */
172 int maxim; /* Max (im1, im2) */
173 int kstat=0; /* Status variable */
174 int kpos=0; /* Position of error */
175 int newin1, newin2; /* Number of vertices along par. dir. 1 & 2 */
176 int numpt; /* Needed in call to s1357 */
177 double start=0; /* Needed in call to s1357 */
178 double end;
179 int *typept=SISL_NULL; /* Array needed for call to s1357 */
180 double *pointpar=SISL_NULL; /* Array needed for call to s1357 */
181 double *newcoeff=SISL_NULL; /* Array needed for call to s1357 */
182 SISLCurve *curve1=SISL_NULL, *curve2=SISL_NULL;
183
184 /* Allocate necessary array for call to s1357 */
185
186 maxim = MAX( im1, im2 );
187 if((typept = newarray(maxim, INT))==SISL_NULL) goto err101;
188 for(i=0; i<maxim; i++)
189 typept[i] = 1;
190
191 /* Interpolation in second direction */
192
193 s1357(points, im2, idim*im1, typept, par2, con1, con2, iopen2, order2,
194 start, &end, &curve1, &pointpar, &numpt, &kstat);
195 if(kstat < 0) goto error;
196 if(pointpar != SISL_NULL)
197 {
198 freearray(pointpar);
199 pointpar = SISL_NULL;
200 }
201
202 newin2 = curve1->in;
203
204 /* Transpose result, store new coefficients in
205 * array newcoeff */
206
207 if( (newcoeff = newarray(idim * im1 * newin2, DOUBLE)) == SISL_NULL )
208 goto err101;
209 s6chpar(curve1->ecoef, im1, newin2, idim, newcoeff);
210
211 /* Interpolation in first parameter direction */
212
213 s1357(newcoeff, im1, idim*newin2, typept, par1, con3, con4, iopen1, order1,
214 start, &end, &curve2, &pointpar, &numpt, &kstat);
215 if(kstat < 0) goto error;
216 if(pointpar != SISL_NULL)
217 {
218 freearray(pointpar);
219 pointpar = SISL_NULL;
220 }
221
222 newin1 = curve2->in;
223
224 /* Transpose back coefficients */
225
226 if( (newcoeff=increasearray(newcoeff, idim*newin1*newin2, DOUBLE))
227 == SISL_NULL ) goto err101;
228 s6chpar(curve2->ecoef, newin2, newin1, idim, newcoeff);
229
230 /* Create instance of surface */
231
232 if (((*rsurf) = newSurf(newin1, newin2, order1, order2, curve2->et,
233 curve1->et, newcoeff, 1, idim, 1)) == SISL_NULL)
234 goto err101;
235
236 /* Set periodicity flag. */
237
238 (*rsurf)->cuopen_1 = curve2->cuopen;
239 (*rsurf)->cuopen_2 = curve1->cuopen;
240
241 /* Success */
242
243 *jstat = 0;
244 goto out;
245
246 /* Allocation error. */
247
248 err101:
249 *jstat = -101;
250 s6err("s1537",*jstat,kpos);
251 goto out;
252
253 /* Error in lower level routine. */
254
255 error: *jstat =kstat;
256 s6err("s1537",*jstat,kpos);
257 goto out;
258
259 out:
260 /* Free arrays */
261
262 if (newcoeff != SISL_NULL) freearray(newcoeff);
263 if (typept != SISL_NULL) freearray(typept);
264
265 /* Free local SISL-curve objects */
266
267 if (curve1 != SISL_NULL) freeCurve(curve1);
268 if (curve2 != SISL_NULL) freeCurve(curve2);
269
270 return;
271 }
272