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: s1530.c,v 1.2 2001-03-19 15:58:50 afr Exp $
45 *
46 */
47
48
49 #define S1530
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
s1530(double ep[],double eder10[],double eder01[],double eder11[],double epar1[],double epar2[],int im1,int im2,int idim,SISLSurf ** rsurf,int * jstat)54 void s1530(double ep[],double eder10[],double eder01[],double eder11[],
55 double epar1[],double epar2[],
56 int im1,int im2,int idim,
57 SISLSurf **rsurf,int *jstat)
58 #else
59 void s1530(ep,eder10,eder01,eder11,epar1,epar2,
60 im1,im2,idim,rsurf,jstat)
61 double ep[];
62 double eder10[];
63 double eder01[];
64 double eder11[];
65 double epar1[];
66 double epar2[];
67 int im1;
68 int im2;
69 int idim;
70 SISLSurf **rsurf;
71 int *jstat;
72 #endif
73 /*
74 ************************************************************************
75 *
76 * PURPOSE: To compute the cubic Hermite interpolant to the data given.
77 * More specifically, given positions, 10, 01, and 11
78 * derivatives at points of a rectangular grid, the routine
79 * computes a cubic tensor-product B-spline interpolant to
80 * the given data with double knots at each data (the first
81 * knot vector will have double knots at all interior points
82 * in epar1, quadruple knots at the first and last points,
83 * and similarly for the second knot vector).
84 *
85 * INPUT:
86 * ep - Array of dimension idim*im1*im2 containing
87 * the positions of the nodes (using the same ordering
88 * as ecoef in the SISLSurf structure).
89 *
90 * eder10 - Array of dimension idim*im1*im2 containing the
91 * first derivative in the first parameter direction.
92 *
93 * eder01 - Array of dimension idim*im1*im2 containing the
94 * first derivative in the second parameter direction.
95 *
96 * eder11 - Array of dimension idim*im1*im2 containing
97 * the cross derivative (twist vector).
98 *
99 * epar1 - Array of dimension im1 containing the
100 * parametrization in the first direction.
101 *
102 * epar2 - Array of dimension im2 containing the
103 * parametrization in the first direction.
104 *
105 * im1 - The number of interpolation points in the
106 * first parameter direction.
107 *
108 * im2 - The number of interpolation points in the
109 * second parameter direction.
110 *
111 * idim - Dimension of the space we are working in.
112 *
113 * Output:
114 * rsurf - Pointer to the surf produced
115 * jstat - Status variable
116 * < 0 - Error.
117 *
118 * Method:
119 * The interpolation is accomplished by using a one dimensional
120 * routine for cubic Hermite spline interpolation. First, the data
121 * is considered to be im2 positional and derivative vectors on
122 * two curves in idim * im1 dimensional space sampled at the
123 * points of epar2.
124 * The first of these has position vectors given by ep and
125 * derivative vectors given by eder01, the second position vectors
126 * given by eder10 and derivative vectors given by eder11.
127 * Running these curves through the one dimensional cubic Hermite
128 * spline interpolation routine then produces two cubic splines rpos
129 * and rder with coefficients of dimension (idim * im1) * (2 * im2)
130 * on the knot vector et2 which is just the points of epar2 with
131 * multiplicity 2 for the interior points and 4 for the
132 * end points.
133 * These coefficients are then considered to be im1 position vectors
134 * and derivative vectors on a curve in 2*idim*im2 dimensional
135 * space (after an appropriate tranposition) sampled at the
136 * points of epar1. Running this data through the one dimensional
137 * cubic Hermite spline routine results in a cubic spline
138 * with coefficients of dimension (2 * idim * im2) * (2 * im1)
139 * with knot vector et1 similar to et2.
140 * A transposition of these coefficients yields the B-spline
141 * coefficients of the bicubic Hermite tensor-product spline
142 * interpolant.
143
144 * REFERENCES :
145 *
146 * CALLS :
147 *
148 * WRITTEN BY : Michael Floater, SI, June 1992.
149 *
150 *********************************************************************
151 */
152 {
153 SISLCurve *rpos; /* Curve of positions of dimension idim*im1 */
154 SISLCurve *rder; /* Curve of derivatives of dimension idim*im1 */
155 SISLCurve *rcurve; /* Curve of dimension idim*im2 */
156 int kstat=0; /* Status variable */
157 int kpos=0; /* Position of error */
158 double *ph=SISL_NULL; /* Transposed positions (in rpos) */
159 double *dh=SISL_NULL; /* Transposed derivatives (in rder) */
160 double *scoef=SISL_NULL; /* Transposed positions in rcurve */
161
162
163
164 /* Check input */
165
166 if (im1 < 2 || im2 < 2 || idim < 1) goto err102;
167
168
169 /* Interpolate position and derivative
170 in the second parameter direction. */
171
172 s1379(ep,eder01,epar2,im2,idim * im1,&rpos,&kstat);
173 if(kstat < 0) goto error;
174
175 s1379(eder10,eder11,epar2,im2,idim * im1,&rder,&kstat);
176 if(kstat < 0) goto error;
177
178
179
180 /* Transpose the matrices (ecoef arrays) of the curves. */
181
182 s1531(rpos->ecoef, idim, im1, rpos->in, &ph, &kstat);
183 if(kstat < 0) goto error;
184
185 s1531(rder->ecoef, idim, im1, rder->in, &dh, &kstat);
186 if(kstat < 0) goto error;
187
188
189
190 /* Interpolate in the first parameter direction. */
191
192 s1379(ph,dh,epar1,im1,idim * rpos->in,&rcurve,&kstat);
193 if(kstat < 0) goto error;
194
195 /* Transpose coefficient matrix of rcurve. */
196
197 s1531(rcurve->ecoef, idim, rpos->in, rcurve->in, &scoef, &kstat);
198 if(kstat < 0) goto error;
199
200 /* Create instance of surface. */
201
202 (*rsurf) = newSurf(rcurve->in,rpos->in,rcurve->ik,rpos->ik,
203 rcurve->et,rpos->et,scoef,1,idim,1);
204 if((*rsurf) == SISL_NULL) goto err101;
205
206 /* Set periodicity flag. */
207
208 (*rsurf)->cuopen_1 = rcurve->cuopen;
209 (*rsurf)->cuopen_2 = rpos->cuopen;
210
211 /* Bicubic surface made. */
212
213 *jstat = 0;
214 goto out;
215
216 /* Error in space allocation. */
217
218 err101: *jstat = -101;
219 s6err("s1530",*jstat,kpos);
220 goto out;
221
222
223 /* Error in input data. */
224
225 err102: *jstat = -102;
226 s6err("s1530",*jstat,kpos);
227 goto out;
228
229
230 /* Error in lower level routine. */
231
232 error: *jstat =kstat;
233 s6err("s1530",*jstat,kpos);
234 goto out;
235
236 out:
237 if (rpos != SISL_NULL) freeCurve(rpos);
238 if (rder != SISL_NULL) freeCurve(rder);
239 if (rcurve != SISL_NULL) freeCurve(rcurve);
240 if (scoef != SISL_NULL) freearray(scoef);
241 if (ph != SISL_NULL) freearray(ph);
242 if (dh != SISL_NULL) freearray(dh);
243
244 return;
245 }
246