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: s1896.c,v 1.3 2001-03-19 15:58:55 afr Exp $
45 *
46 */
47
48
49 #define S1896
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1896(SISLSurf * osurf,double earray[],int dimp1,int narr,int ders1[],int dert1[],int ders2[],int dert2[],SISLSurf ** nsurf,int * jstat)55 s1896 (SISLSurf * osurf, double earray[], int dimp1, int narr, int ders1[],
56 int dert1[], int ders2[], int dert2[], SISLSurf ** nsurf, int *jstat)
57
58 #else
59 void
60 s1896 (osurf, earray, dimp1, narr, ders1, dert1, ders2, dert2, nsurf, jstat)
61 SISLSurf *osurf;
62 double earray[];
63 int dimp1;
64 int narr;
65 int ders1[];
66 int dert1[];
67 int ders2[];
68 int dert2[];
69 SISLSurf **nsurf;
70 int *jstat;
71
72 #endif
73 /*
74 *********************************************************************
75 *
76 *********************************************************************
77 *
78 * PURPOSE : For each of the planes in the matrix we make the function
79 *
80 * iders1 idert1 iders2 idert2
81 * /D \ /D \ /D \ /D \
82 * F(s,t) = -- * -- P(s,t) earray -- * -- P(s,t)
83 * \DS/ \DT/ \DS/ \DT/
84 *
85 * Here:
86 * iders1=ders1[i], iders2=ders2[i]
87 * idert1=dert1[i], idert2=dert2[i]
88 *
89 * If earray is the identity matrix then we make the dor product
90 * of (iders1,idert1)-derivateive and (iders2,idert2)-derivative
91 * of the B-spline "surface" P(s,t) described by osurf.
92 *
93 * INPUT : osurf - The original B-pline surface.
94 * earray - The array to be used.
95 * dimp1 - The dimension idim+1. This is due to that the matrix
96 * earray is used for calculation with homogenous
97 * coordinates.
98 * narr - Number of parallel arrays.
99 * ders1 - The derivatives to be applied in the first parameter
100 * direction in the first occurence of the surface.
101 * dert1 - The derivatives to be applied in the second parameter
102 * direction in the first occurence of the surface.
103 * ders2 - The derivatives to be applied in the first parameter
104 * direction in the second occurence of the surface.
105 * dert2 - The derivatives to be applied in the second parameter
106 * direction in the second occurence of the surface.
107 *
108 * OUTPUT : nsurf - The new B-spline surface.
109 * jstat - Status variable
110 * > 0 : warning
111 * = 0 : ok
112 * < 0 : error
113 *
114 * METHOD : The product of the (iders1,idert1)-derivative of a B-spline
115 * curve, a matrix and the (iders2,idert2)-derivative of the same
116 * B-spline curve is a function of order 2*ik1-ders1-ders2-1 in
117 * the first parameter direction and 2*ik2-dert1-dert2-1 in the
118 * second parameter direction, if the order of the B-spline
119 * surface is ik1,ik2.
120 * ds1 = min( ders1[i] ), i=1,...,narr
121 * dt1 = min( dert1[i] ), i=1,...,narr
122 * ds2 = min( ders2[i] ), i=1,...,narr
123 * dt2 = min( dert2[i] ), i=1,...,narr
124 *
125 * The subroutine have the following steps.
126 * 1. Calculate knot vectors of B-spline functions.
127 * 2. Calculate parameter value pairs to be used for calculation
128 * of points on the B-spline function.
129 * 3. Interpolate the calculated values on the B-spline function.
130 *
131 * Note: The parameter values of the B-spline function values must
132 * be calculated to ensure reproduction of the B-spline function.
133 *
134 * REFERENCES : Fortran version:
135 * Tor Dokken, SI, 1983-05
136 *
137 * CALLS : s1894,s1890,s1421,s1891,s6err.
138 *
139 * WRITTEN BY : Trond Vidar Stensby, SI, 1991-06
140 * REVISED BY : Johannes Kaasa, SI, May 92 (Changed the order of the indexes
141 * in the multiplication of earray, val1 and val2. In addition I
142 * changed the sequence in the second assignment of ds2 and dt2.
143 * I also had to put on paranthesises in the expression for pos1
144 * and pos2).
145 * REVISED BY : Michael Floater, SI, June 92. The rational stuff
146 * was completely messed up. But it works now.
147 * Dimension of output curve was wrong -- now narr.
148 * Revised by : Paal Fugelli, SINTEF, Oslo, Norway, August 94. Fixed over-running
149 * of array indices (found with Purify).
150 *
151 *********************************************************************
152 */
153 {
154 int nik1; /* Order of new surface in
155 first parameter direction. */
156 int nin1; /* Order of new surface in
157 second parameter direction. */
158 int nik2; /* Number of vertices in first
159 parameter direction. */
160 int nin2; /* Number of vertices in second parameter direction. */
161 int lfs; /* Interval indicator. (left side) */
162 int lft; /* Interval indicator. (left side) */
163 int tpos; /* Used to index array tau. */
164 int epos; /* Used to index earray. */
165 int pos1; /* Position of values of first derivatives. */
166 int pos2; /* Position of values of second derivatives. */
167 int ds1; /* Order of derivatives. */
168 int dt1;
169 int ds2;
170 int dt2;
171 int mds1; /* Maximum order of derivatives. */
172 int mdt1;
173 int mds2;
174 int mdt2;
175 int nder1; /* Total order of derivatives.
176 (Both directions) */
177 int nder2;
178 int dim; /* Dimension of tau. */
179 int maxder; /* Largest total order of derivatives.
180 (Both functions.) */
181 int count1; /* Loop control variables. */
182 int kj, ki;
183 int kl, kr, kp;
184 double parval[2];
185 double sum; /* Used for calculation of P(s,t). */
186 double *nknots1 = SISL_NULL; /* New knots in first parameter direction. */
187 double *nknots2 = SISL_NULL; /* New knots in second parameter direction. */
188 double *coef1 = SISL_NULL; /* New coeficients */
189 double *coef2 = SISL_NULL; /* New coeficients */
190 double *par1 = SISL_NULL; /* Parameter values in first direction. */
191 double *par2 = SISL_NULL; /* Parameter values in second direction. */
192 int *der1 = SISL_NULL; /* Derivative indicators in first direction. */
193 int *der2 = SISL_NULL; /* Derivative indicators in second direction.*/
194 double *deriv = SISL_NULL; /* Derivatives returned by s1421. */
195 double *normal = SISL_NULL; /* Normal returned by s1421. (not used) */
196 double *val1 = SISL_NULL; /* Values extracted from deriv. */
197 double *val2 = SISL_NULL; /* Values extracted from deriv. */
198 double *tau = SISL_NULL; /* Interpolation points. */
199 int kstat = 0;
200 int kpos = 0;
201
202 *jstat = 0;
203
204 /* Test if legal input. */
205
206 if (osurf->ik1 <= 1 || osurf->in1 < osurf->ik1)
207 goto err112;
208 if (osurf->ik2 <= 1 || osurf->in2 < osurf->ik2)
209 goto err112;
210
211 /* Find minimal and maximal order of derivatives */
212
213 ds1 = mds1 = ders1[0];
214 dt1 = mdt1 = dert1[0];
215 ds2 = mds2 = ders2[0];
216 dt2 = mdt2 = dert2[0];
217
218 for (count1 = 1; count1 < narr; count1++)
219 {
220 if (ds1 > ders1[count1]) ds1 = ders1[count1];
221 if (dt1 > dert1[count1]) dt1 = dert1[count1];
222 if (ds2 > ders2[count1]) ds2 = ders2[count1];
223 if (dt2 > dert2[count1]) dt2 = dert2[count1];
224
225 if (mds1 < ders1[count1]) mds1 = ders1[count1];
226 if (mdt1 < dert1[count1]) mdt1 = dert1[count1];
227 if (mds2 < ders2[count1]) mds2 = ders2[count1];
228 if (mdt2 < dert2[count1]) mdt2 = dert2[count1];
229 }
230
231 /* Produce a knot vector in the first parameter direction. */
232
233 s1894 (osurf->et1, osurf->ik1, osurf->in1, ds1, ds2, earray, dimp1, narr,
234 &nknots1, &nik1, &nin1, &kstat);
235 if (kstat < 0) goto error;
236
237 /* Produce a knot vector in second parameter direction. */
238
239 s1894 (osurf->et2, osurf->ik2, osurf->in2, dt1, dt2, earray, dimp1, narr,
240 &nknots2, &nik2, &nin2, &kstat);
241 if (kstat < 0) goto error;
242
243 /* Produce parameter values and derivative indicators in first
244 * parameter direction. */
245
246 s1890 (nknots1, nik1, nin1, &par1, &der1, &kstat);
247 if (kstat < 0) goto error;
248
249 /* Produce parameter values and derivative indicators in second
250 * parameter direction. */
251
252 s1890 (nknots2, nik2, nin2, &par2, &der2, &kstat);
253 if (kstat < 0) goto error;
254
255 /* Allocate memory for point calculation. */
256
257 val1 = newarray (dimp1, DOUBLE);
258 if (val1 == SISL_NULL) goto err101;
259 val2 = newarray (dimp1, DOUBLE);
260 if (val2 == SISL_NULL) goto err101;
261 tau = newarray (narr * nin1 * nin2, DOUBLE);
262 if (tau == SISL_NULL) goto err101;
263 maxder = max (max (mds1, mds2), max (mdt1, mdt2));
264 deriv = newarray (osurf->idim * (maxder + 1) * (maxder + 2) / 2, DOUBLE);
265 if (deriv == SISL_NULL) goto err101;
266 normal = newarray (osurf->idim * (maxder + 1) * (maxder + 2) / 2, DOUBLE);
267 if (normal == SISL_NULL) goto err101;
268
269 /* Calculate interpolation points. */
270
271 lfs = 0;
272 lft = 0;
273 tpos = 0;
274 for (kj = 0; kj < nin2; kj++)
275 {
276 parval[1] = par2[kj];
277 for (ki = 0; ki < nin1; ki++)
278 {
279 parval[0] = par1[ki];
280 epos = 0;
281 for (kl = 0; kl < narr; kl++)
282 {
283 ds1 = ders1[kl];
284 dt1 = dert1[kl];
285 ds2 = ders2[kl];
286 dt2 = dert2[kl];
287
288 /* ds2 = dert2[kl];
289 dt2 = ders2[kl]; */
290
291 maxder = max (max (ds1, ds2), max (dt1, dt2));
292
293 s1421 (osurf, maxder, parval, &lfs, &lft, deriv, normal, &kstat);
294 if (kstat < 0) goto error;
295
296 nder1 = ds1 + dt1;
297 nder2 = ds2 + dt2;
298 pos1 = osurf->idim * (nder1 * (nder1 + 1) / 2 + dt1);
299 pos2 = osurf->idim * (nder2 * (nder2 + 1) / 2 + dt2);
300
301 for (count1 = 0; count1 < osurf->idim; count1++)
302 {
303 val1[count1] = deriv[pos1++];
304 val2[count1] = deriv[pos2++];
305 }
306 if (osurf->idim < dimp1)
307 {
308 val1[osurf->idim] = (double) 1.0;
309 val2[osurf->idim] = (double) 1.0;
310 if (ds1 > 0 || dt1 > 0)
311 val1[osurf->idim] = (double) 0.0;
312 if (ds2 > 0 || dt2 > 0)
313 val2[osurf->idim] = (double) 0.0;
314 }
315
316 /* Can now calculate a interpolation point. */
317
318 sum = (double) 0.0;
319 for (kr = 0; kr < dimp1; kr++, epos += dimp1)
320 {
321 for (kp = 0; kp < dimp1; kp++)
322 sum += earray[epos + kp] * val1[kr] * val2[kp];
323 /* sum += earray[epos + kp] * val1[kp] * val2[kr]; */
324 }
325 tau[tpos++] = sum;
326 }
327 }
328 }
329
330 /* Calculate new surface description. */
331
332 /* Interpolate in second parameter direction. */
333
334 dim = narr * nin1;
335
336 s1891 (par2, tau, dim, nin2, 1, der2, TRUE, nknots2, &coef1, &nin2,
337 nik2, 0, 0, &kstat);
338 if (kstat < 0) goto error;
339
340 /* Interpolate in first parameter direction. */
341
342 s1891 (par1, coef1, narr, nin1, nin2, der1, TRUE, nknots1, &coef2,
343 &nin1, nik1, 0, 0, &kstat);
344 if (kstat < 0) goto error;
345
346 /* OK */
347
348 *nsurf = newSurf (nin1, nin2, nik1, nik2, nknots1, nknots2,
349 coef2, osurf->ikind, narr, 2);
350 if (*nsurf == SISL_NULL) goto err171;
351
352 goto out;
353
354 /* Not enough memory. */
355
356 err101:
357 *jstat = -101;
358 s6err ("s1896", *jstat, kpos);
359 goto out;
360
361 /* Could not create surface, */
362
363 err171:
364 *jstat = -171;
365 s6err ("s1896", *jstat, kpos);
366 goto out;
367
368 /* Error in description of B-spline. */
369
370 err112:
371 *jstat = -112;
372 s6err ("s1896", *jstat, kpos);
373 goto out;
374
375 /* Error in lower level routine. */
376
377 error:
378 *jstat = kstat;
379 s6err ("s1896", *jstat, kpos);
380 goto out;
381
382 /* Free pointers. */
383
384 out:
385 if (coef1 != SISL_NULL) freearray (coef1);
386 if (val1 != SISL_NULL) freearray (val1);
387 if (val2 != SISL_NULL) freearray (val2);
388 if (par1 != SISL_NULL) freearray (par1);
389 if (par2 != SISL_NULL) freearray (par2);
390 if (der1 != SISL_NULL) freearray (der1);
391 if (der2 != SISL_NULL) freearray (der2);
392 if (normal != SISL_NULL) freearray (normal);
393 if (deriv != SISL_NULL) freearray (deriv);
394 if (tau != SISL_NULL) freearray (tau);
395
396 return;
397 }
398