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: s1891.c,v 1.2 2001-03-19 15:58:55 afr Exp $
45 *
46 */
47
48
49 #define S1891
50
51 #include "sislP.h"
52 #define MAX_SIZE 50
53
54
55 #if defined(SISLNEEDPROTOTYPES)
56 void
s1891(double etau[],double epoint[],int idim,int inbpnt,int iright,int eder[],int iopen,double et[],double * ebcoef[],int * in,int ik,int inlr,int inrc,int * jstat)57 s1891 (double etau[], double epoint[], int idim, int inbpnt, int iright,
58 int eder[], int iopen, double et[], double *ebcoef[], int *in,
59 int ik, int inlr, int inrc, int *jstat)
60
61 #else
62 void
63 s1891 (etau, epoint, idim, inbpnt, iright, eder, iopen, et, ebcoef,
64 in, ik, inlr, inrc, jstat)
65 double etau[];
66 double epoint[];
67 int idim;
68 int inbpnt;
69 int iright;
70 int eder[];
71 int iopen;
72 double et[];
73 double *ebcoef[];
74 int *in;
75 int ik;
76 int inlr;
77 int inrc;
78 int *jstat;
79 #endif
80 /*
81 *********************************************************************
82 *
83 *********************************************************************
84 *
85 * PURPOSE : To compute the B-spline coefficients required for
86 * interpolation of epoint with a B-spline.
87 *
88 * INPUT : etau - Parameter values (i.e. parametrization of
89 * points in epoint) calculated in s1890.
90 * epoint - Point array. Contains points/derivatives.
91 * idim - Dimension of the space in which the points
92 * lie.
93 * inbpnt - Number of points/derivatives
94 * iright - Number of right hand sides of dimension
95 * (idim,inbpnt)
96 * eder - Parametrization of derivatives in epoint
97 * (calculated in s1890)
98 * iopen - Open or closed curve.
99 * et - Knot vector of length (in+ik)
100 * Open curve: in = inbpnt
101 * Closed curve: in = inbpnt+ik-1
102 * ik - Order of B-spline basis to be used
103 * inlr - Indicates shape of matrix required for B-spline
104 * interpolation (see subroutines s1925 & s1926).
105 * inrc - Indicates shape of matrix required for B-spline
106 * interpolation (see subroutines s1925 & s1926).
107 *
108 * OUTPUT : ebcoef - Guiding points for B-spline curve.
109 * Dimension (1:idim*in*iright)
110 * in - Number of vertices in the B-spline curve
111 * produced. Open curve: in=inbpnt
112 * Closed curve: in=inbpnt+ik-1.
113 * jstat - Status variable:
114 * < 0 : error
115 * > 0 : warning
116 * = 0 : OK.
117 *
118 * METHOD : See: Carl De Boor: "A practical guide to splines"
119 *
120 * REFERENCES : Fortran version:
121 * E. Aarn[s, CP, 1979-01
122 *
123 * CALLS : s1925,s6err.
124 *
125 * WRITTEN BY : Christophe R. Birkeland 1991-06
126 *
127 *********************************************************************
128 */
129 {
130 int kstat = 0;
131 int kpos = 0; /* Position of error */
132 int ii; /* Loop control parameter */
133 int limit1, limit2; /* Loop parameters */
134 int kj, kl;
135 int kdum, stop;
136 int nur; /* Number of upper rows in W */
137 int inlx; /* Equal to inlr if inlr>0, else=1 */
138 int inrx; /* Equal to inrc if inrc>0, else=1 */
139 int edarray[MAX_SIZE]; /* Array for ed below */
140 int alloc_needed=FALSE;
141 int *ed = SISL_NULL; /* Arrays defining elements of W */
142 double *ewarray=SISL_NULL; /* Array for ew1, ew2 and ew3 */
143 double *ew1 = SISL_NULL; /* See subroutine s1926 */
144 double *ew2 = SISL_NULL;
145 double *ew3 = SISL_NULL;
146
147 *jstat = 0;
148
149
150 /* Test if legal input. */
151
152 if (ik < 1 || idim < 1) goto err112;
153
154 /* Indicate dimension of B-spline. */
155
156 *in = inbpnt;
157 if (iopen != SISL_CRV_OPEN) *in +=ik - 1;
158
159 *ebcoef = new0array (*in *idim * iright, DOUBLE);
160 if (*ebcoef == SISL_NULL) goto err101;
161
162 if ((nur = inbpnt - inlr) > MAX_SIZE)
163 alloc_needed = TRUE;
164
165 /* Allocate arrays ew1, ew2, ew3, ed. */
166
167 inlx = MAX (1, inlr);
168 inrx = MAX (1, inrc);
169 limit1 = (ik * nur) + (inrx * nur) + (inlx * inbpnt);
170
171 if ((ewarray = new0array(limit1 + 1,DOUBLE)) == SISL_NULL) goto err101;
172
173 ew1 = ewarray;
174 ew2 = ew1 + (ik * nur);
175 ew3 = ew2 + (inrx * nur);
176
177 if (alloc_needed)
178 {
179 if ((ed = new0array(nur,INT)) == SISL_NULL)
180 goto err101;
181 }
182 else
183 ed = edarray;
184
185 s1925 (etau, epoint, inbpnt, eder, et, *ebcoef,*in, ik, iright,
186 idim, ew1, nur, ed, ew2, inrc, ew3, inlr, &kstat);
187 if (kstat < 0) goto error;
188
189 /* For closed B-spline curves we have:
190 * ebcoef(i) = ebcoef(i+inbpnt) ; i=1,...,ik-1. */
191
192 if (iopen != SISL_CRV_OPEN)
193 {
194 stop = ik - 1;
195 for (kl = 0; kl < iright; kl++)
196 {
197 kdum = *in *kl;
198 for (kj = 0; kj < stop; kj++)
199 {
200 limit2 = (kj + kdum) * idim;
201 limit1 = inbpnt * idim + limit2;
202 for (ii = 0; ii < idim; ii++)
203 (*ebcoef)[limit1 + ii] = (*ebcoef)[limit2 + ii];
204 }
205 }
206 }
207
208 goto out;
209
210 /* Error in lower level routine */
211
212 error:
213 *jstat = kstat;
214 s6err ("s1891", *jstat, kpos);
215 goto out;
216
217 /* Error in array allocations */
218
219 err101:
220 *jstat = -101;
221 s6err ("s1891", *jstat, kpos);
222 goto out;
223
224 /* Error in description of B-spline */
225
226 err112:
227 *jstat = -112;
228 s6err ("s1891", *jstat, kpos);
229 goto out;
230
231 out:
232 if (alloc_needed) freearray (ed);
233 if (ewarray) freearray (ewarray);
234 return;
235 }
236 #undef MAX_SIZE
237