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: s6strider.c,v 1.3 2001-03-19 15:59:02 afr Exp $
45 *
46 */
47
48
49 #define S6STRIDER
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s6strider(double eder[],int idim,int ider,double gder[],int * jstat)55 s6strider(double eder[],int idim,int ider,double gder[],int *jstat)
56 #else
57 void s6strider(eder,idim,ider,gder,jstat)
58 double eder[];
59 int idim;
60 int ider;
61 double gder[];
62 int *jstat;
63 #endif
64 /*
65 *********************************************************************
66 *
67 * PURPOSE : To calculate the value and ider*ider derivatives of
68 * a rational B-spline surface.
69 *
70 * INPUT : eder - Double array of dimenson [(ider+1)*(ider+2)*(idim+1)/2]
71 * containing the position and the derivative vectors
72 * of the homogeneous surface at the point with parameter value
73 * (epar[0],epar[1]).
74 * (idim+1 is the number of components of each B-spline
75 * coefficient, i.e. the dimension of the homogemous
76 * space in which the surface lies.)
77 * These vectors are stored in the following order:
78 * First the idim+1 components of the position vector,
79 * then the idim+1 components of the D(1,0) vector,
80 * then the idim+1 components of the D(0,1) vector,
81 * then the idim+1 components of the D(2,0) vector,
82 * followed by D(1,1), D(0,2)
83 * and so on up to the idim+1 components of the D(0,ider).
84 * idim - The dimension of the non homogenous space
85 * ider - The number of input derivatives with respect
86 * to both parameter directions.
87 *
88 *
89 * OUTPUT : jstat - Status message
90 * >0 : Warning
91 * =0 : ok
92 * <0 : Error
93 * gder - Double array of dimension [(ider+1)*(ider+2)*idim/2]
94 * containing the position and the derivative vectors
95 * of the surface at the point with parameter value
96 * (epar[0],epar[1]).
97 * (idim is the number of components of each B-spline
98 * coefficient, i.e. the dimension of the Euclidean
99 * space in which the surface lies.)
100 * These vectors are stored in the following order:
101 * First the idim components of the position vector,
102 * then the idim components of the D(1,0) vector,
103 * then the idim components of the D(0,1) vector,
104 * then the idim components of the D(2,0) vector,
105 * followed by D(1,1), D(0,2)
106 * and so on up to the idim components of the D(0,ider).
107 *
108 *
109 * METHOD : The surface P(u,v) can be written as the quotient
110 * P(u,v) = T(u,v) / w(u,v) where T and w are ordinary splines.
111 * The dimensions of T and w are idim and 1
112 * respectively. The array eder contains position
113 * and derivatives of the idim+1 dimensional surface
114 * (T(u,v),w(u,v)).
115 *
116 * Now, since wP = T, we find, by the Leibnitz formula,
117 *
118 * k l
119 * k! l! (k-i,l-j) (i,j) (k,l)
120 * sum sum -------- -------- w P = T .
121 * i!(k-i)! j!(l-j)!
122 * i=0 j=0
123 *
124 * Therefore
125 *
126 *
127 * -- k l --
128 * (k,l) | (k,l) k! l! (k-i,l-j) (i,j) |
129 * P = | T - sum sum -------- -------- w P | / w .
130 * | i!(k-i)! j!(l-j)! |
131 * -- i=0 j=0 --
132 * i+j<k+l
133 *
134 * This formula is applied recursively to evaluate P's derivatives.
135 *
136 * MF.
137 *
138 *
139 *
140 * CALLS :
141 *
142 * WRITTEN BY : Michael Floater, SI, 3.9.92.
143 * Essentially the same as s6sratder
144 * except that we work with triangular matrices
145 * ((0,0), (1,0), (0,1), (2,0), (1,1), ...)
146 * instead of rectangular ones.
147 * Revised by : Christophe Rene Birkeland, SINTEF Oslo, May 1993.
148 * Error message corrected
149 *
150 *********************************************************************
151 */
152 {
153 int kpos=0; /* Position of error. */
154 double w0; /* The denominator. */
155 int ki; /* Count through dimensions. */
156 int idu; /* Count through derivatives in u. */
157 int idv; /* Count through derivatives in v. */
158 int *binom=SISL_NULL; /* Array for binomial coefficients. */
159 int *binomu=SISL_NULL; /* Pointer to binomial coefficients in u. */
160 int *binomv=SISL_NULL; /* Pointer to binomial coefficients in v. */
161 double *sum1=SISL_NULL; /* Leibnitz expansion in u */
162 double *sum2=SISL_NULL; /* Leibnitz expansion in u and v. */
163 double sumdum1[4]; /* Fixed space for sum1. */
164 double sumdum2[4]; /* Fixed space for sum2. */
165 int idimp1; /* idim + 1. */
166 int iw; /* Pointer to a weight. */
167 int igder; /* Pointer to already calculated derivs. */
168 int i,iu,iv,j,k; /* Counters. */
169 int iderp1; /* ider + 1. */
170 int igrow; /* (ider+1) * idim. */
171 int iwrow; /* (ider+1) * idimp1. */
172 int iutemp,ivtemp; /* Used to find next weight in the sum. */
173 int tot,temp1; /* Temporary variables. */
174 int bidum[10]; /* Array for storing binomial coeffs. */
175 double temp; /* Temporary multiple. */
176
177 if (ider<0) goto err178;
178 if (idim<1) goto err102;
179
180 *jstat = 0;
181
182 /* Find denominator. */
183
184 w0 = eder[idim];
185 if (DEQUAL(w0,DZERO)) w0 = (double)1.0;
186
187 /* If we're only asked for position, we'll do it
188 now and exit for the sake of speed. */
189
190 if(ider == 0)
191 {
192 for(ki=0; ki<idim; ki++)
193 gder[ki] = eder[ki] / w0;
194
195 goto out;
196 }
197
198 /* Set up some constants. */
199
200 idimp1 = idim + 1;
201 iderp1 = ider + 1;
202 igrow = iderp1 * idim;
203 iwrow = igrow + iderp1; /* = iderp1 * idimp1 */
204
205 /* Set up binomial coefficients.
206 Use new array only when ider > 3. */
207
208 if (ider > 3)
209 {
210 binom = newarray((iderp1*(iderp1+1)) >> 1, INT);
211 if(binom == SISL_NULL) goto err179;
212 }
213 else
214 {
215 binom = bidum;
216 }
217
218 for(j=0,k=0; j<=ider; j++,k+=j)
219 {
220 /* Calculate the new row of binomial coefficients. */
221
222 binom[k] = 1;
223
224 for(i=k+1; i<k+j; i++)
225 {
226 binom[i] = binom[i-j-1] + binom[i-j];
227 }
228
229 binom[k+j] = 1;
230 }
231
232 /* Set up space for sum1 and sum2 if necessary.
233 Use new arrays only when idim > 4. */
234
235 if (idim > 4)
236 {
237 sum1 = newarray(idim, DOUBLE);
238 if(sum1 == SISL_NULL) goto err179;
239 sum2 = newarray(idim, DOUBLE);
240 if(sum2 == SISL_NULL) goto err179;
241 }
242 else
243 {
244 sum1=sumdum1;
245 sum2=sumdum2;
246 }
247
248 /* Loop through derivatives in u and v. */
249
250 for(idv=0,binomv=binom; idv<=ider; idv++,binomv+=idv)
251 {
252 for(idu=0,binomu=binom; idu<=ider-idv; idu++,binomu+=idu)
253 {
254 if(idu == 0 && idv == 0)
255 {
256 /* Position is a special case. */
257
258 for(ki=0; ki<idim; ki++)
259 gder[ki] = eder[ki] / w0;
260 }
261 else
262 {
263 /* Calculate indices in eder and gder. */
264
265 tot = idu + idv;
266 temp1 = ((tot * (tot+1)) >> 1) + idv;
267
268 j = temp1 * idim;
269 k = j + temp1;
270
271 /* Calculating each coefficient of the (idu,idv)'th
272 derivative of the rational surface (in gder).
273
274 This requires calculating the Liebnitz sum from
275 the subarray of gder (0,..,idu, 0,...,idv) and
276 the subarray of eder (0,..,idu, 0,...,idv). */
277
278 /* Calculate the Leibnitz sum. */
279
280 for(ki=0; ki<idim; ki++)
281 sum2[ki] = (double)0.0;
282
283 for(iv=0; iv<=idv; iv++)
284 {
285 for(ki=0; ki<idim; ki++)
286 sum1[ki] = (double)0.0;
287 ivtemp = idv-iv;
288
289 for(iu=0; iu<=idu; iu++)
290 {
291 tot = iu + iv;
292 temp1 = ((tot * (tot+1)) >> 1) + iv;
293
294 igder = temp1 * idim;
295 iutemp = idu-iu;
296
297 tot = iutemp + ivtemp;
298 temp1 = ((tot * (tot+1)) >> 1) + ivtemp;
299
300 iw = temp1 * idimp1 + idim;
301
302 /* Add the next Leibnitz term unless we
303 have reached the last one (the unknown). */
304
305 if(iu<idu || iv<idv)
306 {
307 /* If iu=0 or iu=idu, the u binomial
308 coefficient is 1 so don't multiply. */
309
310 if(iu>0 && iu<idu)
311 {
312 temp = (double)binomu[iu] * eder[iw];
313 for(ki=0; ki<idim; ki++,igder++)
314 sum1[ki] += temp * gder[igder];
315 }
316 else
317 for(ki=0; ki<idim; ki++,igder++)
318 sum1[ki] += eder[iw] * gder[igder];
319 }
320 }
321
322 /* If iv=0 or iv=idv, the v binomial
323 coefficient is 1 so don't multiply. */
324
325 if(iv>0 && iv<idv)
326 for(ki=0; ki<idim; ki++)
327 sum2[ki] += (double)binomv[iv] * sum1[ki];
328 else
329 for(ki=0; ki<idim; ki++)
330 sum2[ki] += sum1[ki];
331 }
332 for(ki=0; ki<idim; ki++,j++,k++)
333 gder[j] = (eder[k] - sum2[ki]) / w0;
334 }
335 }
336 }
337
338 /* Free arrays. */
339
340 if (ider > 3 && binom != SISL_NULL)
341 freearray(binom);
342
343 if (idim > 4)
344 {
345 if(sum1 != SISL_NULL) freearray(sum1);
346 if(sum2 != SISL_NULL) freearray(sum2);
347 }
348
349 /* Done. */
350
351 goto out;
352
353 /* idim less than 1. */
354
355 err102:
356 *jstat = -102;
357 s6err("s6strider",*jstat,kpos);
358 goto out;
359
360 /* Derivative negative */
361
362 err178:
363 *jstat = -178;
364 s6err("s6strider",*jstat,kpos);
365 goto out;
366
367 /* Not enough memory */
368
369 err179:
370 *jstat = -179;
371 s6err("s6strider",*jstat,kpos);
372 goto out;
373
374 out:
375 return;
376 }
377
378