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: s1233.c,v 1.4 2005-02-28 09:04:48 afr Exp $
45 *
46 */
47
48
49 #define S1233
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1233(SISLCurve * pc,double afak1,double afak2,SISLCurve ** rc,int * jstat)55 s1233(SISLCurve *pc,double afak1,double afak2,SISLCurve **rc,int *jstat)
56 #else
57 void s1233(pc,afak1,afak2,rc,jstat)
58 SISLCurve *pc;
59 double afak1;
60 double afak2;
61 SISLCurve **rc;
62 int *jstat;
63 #endif
64 /*
65 *********************************************************************
66 *
67 * PURPOSE : To extend a B-spline curve (i.e. NOT rationals) at the start
68 * and/or the end of the curve by continuing the polynomial
69 * behaviour of the curve.
70 *
71 * INPUT : pc - Pointer to the curve that will be extended.
72 * afak1 - How much the curve is to be stretched at the
73 * start of the curve. The length of the stretched
74 * curve will be equal to (1+afak1) times the
75 * input curve. afak1 >= 0 and will be set to 0 if
76 * negative.
77 * afak2 - How much the curve is to be stretched at the
78 * end of the curve. The length of the stretched
79 * curve will be equal to (1+afak2) times the
80 * input curve. afak2 >= 0 and will be set to 0 if
81 * negative.
82 *
83 * OUTPUT : rc - Pointer to the extended curve.
84 * jstat - status messages
85 * > 0 : warning
86 * = 1 : Stretching factors less than 0 - readjusted
87 * factor(s) have been used.
88 * = 0 : ok
89 * < 0 : error
90 *
91 * METHOD : 1. The new knot vector is made
92 * 2. The transformation matrix for the ik first vertices between
93 * the new and the old knot vector is made.
94 * 3. The transformation matrix is inverted and used to update the
95 * ik first vertices.
96 * 4. The transformation matrix for the ik last vertices between
97 * the new and the old knot vector is made.
98 * 5. The transformation matrix is inverted and used to update the
99 * ik last vertices.
100 * 6. The knot vector is updated.
101 *-
102 * CALLS : make_cv_kreg,s1219,s1701,s6lufacp,s6lusolp,s6err.
103 *
104 * Written by : Paal Fugelli, SINTEF, Oslo, Norway, Sept-1992. Adapted from
105 * a FORTRAN based s1233() written by Vibeke Skytt and s1333_cyclic()
106 * written by Tor Dokken.
107 *
108 *********************************************************************
109 */
110 {
111 double *ext = SISL_NULL; /* Extended version of knot vector */
112 double *smatrix = SISL_NULL; /* Matrix converting between basises */
113 double *salloc = SISL_NULL; /* Matrix for memory allocation */
114 double *salfa = SISL_NULL; /* The values of a discrete B-spline
115 calculation */
116 double *spek = SISL_NULL; /* Pointer used in traversing arrays */
117 double *sb = SISL_NULL; /* Right hand side of equation */
118 double *sfrom, *sto;
119 double *sp; /* Help array for s1701 */
120 double *stx = SISL_NULL; /* Knot vector after insertion of knots
121 at start */
122 int *mpiv=SISL_NULL; /* Pointer to pivotation array */
123 SISLCurve *kreg; /* k-regular curve */
124
125 double *st = SISL_NULL; /* Internal version of et */
126 double *scoef = SISL_NULL; /* Copy of the vertices of the surface */
127 int kdim = pc->idim;
128 int kk = pc->ik;
129 int kn = pc->in;
130
131 double tlen; /* Length of curve parameterization */
132 double tstart; /* New start knot value */
133 double tend; /* New end knot value */
134 int ki, kj;
135 int knst;
136 int knstx;
137 int kleft = 0;
138 int kpl, kfi, kla;
139
140 int kstat; /* Error status from lower level */
141 int kpos = 0;
142
143
144
145 /* Ensure reasonable return value */
146 *rc = SISL_NULL;
147
148
149 /* Test input */
150
151 if ( kk < 1 ) goto err110;
152
153 if ( kn < kk ) goto err111;
154
155 if ( afak1 < DZERO || afak2 < DZERO )
156 {
157 /* Warning - so correct the factor(s) */
158 *jstat = 1;
159
160 if ( afak1 < DZERO ) afak1 = DZERO;
161 if ( afak2 < DZERO ) afak2 = DZERO;
162 }
163
164
165 /* Ensure k-regular curve */
166
167 make_cv_kreg(pc, &kreg, &kstat);
168 if ( kstat < 0 ) goto error;
169
170
171 /* Alloocate array for pivotation vector */
172
173 mpiv = new0array(2*kk, INT);
174 if ( mpiv == SISL_NULL ) goto err101;
175
176 /* Allocate space (en bloc) for the local vectors */
177
178 salloc = new0array(3*kn + 9*kk + 4*kk*kk + kdim*kn, DOUBLE);
179 if ( salloc == SISL_NULL ) goto err101;
180
181 ext = salloc; /* Size kn+kk */
182 smatrix = ext + kn + kk; /* Max size 4*kk*kk */
183 salfa = smatrix + 4*kk*kk; /* Size kk */
184 scoef = salfa + kk; /* Size kdim*kn */
185 sb = scoef + kdim*kn; /* Size 2*kk */
186 sp = sb + 2*kk; /* Size kk */
187 st = sp + kk; /* Size kn + 2*kk */
188 stx = st + kn + 2*kk; /* Size kn + kk */
189
190
191
192 /* Copy knots and vertices, to avoid destruction of curve */
193
194 memcopy(ext, kreg->et, kn + kk, DOUBLE);
195 memcopy(scoef, kreg->ecoef, kdim*kn, DOUBLE);
196
197
198 /* Make extended knot vector */
199
200 tlen = ext[kn] - ext[kk-1];
201
202 if ( afak1 > DZERO )
203 {
204 /* Extend the basis at the start of the curve */
205
206 tstart = ext[kk-1] - tlen*afak1;
207 for ( ki = 0; ki < kk; ki++ ) ext[ki] = tstart;
208 }
209
210 if ( afak2 > DZERO )
211 {
212 /* Extend the basis at the end of the curve */
213
214 tend = ext[kn] + tlen*afak2;
215 for ( ki = kn; ki < kn+kk; ki++ ) ext[ki] = tend;
216 }
217
218
219 /* s1701 expects et to be a refinement of ext, thus we have to make a new
220 version of e1 with the extra kk new knots before the start and after the
221 end and one intermediate version with only kk at the start */
222
223 memcopy(st, ext, kk - 1, DOUBLE);
224 memcopy(st + kk - 1, kreg->et, kn + kk, DOUBLE);
225 memcopy(st + 2*kk - 1 + kn, ext + kn + 1, kk - 1, DOUBLE);
226 knst = kn + 2*(kk - 1);
227
228 memcopy(stx, ext ,kn, DOUBLE);
229 memcopy(stx + kn, st + kk - 1 + kn, 2*kk - 1, DOUBLE);
230 knstx = kn + kk - 1;
231
232 /* STEP 2 Make matrix going between bases, only the kk first and last knots
233 are to be changed. */
234
235
236 /* Now we have two cases. We know that only the kk+1 first and kk+1
237 last vertices are to be changed. However 2*(kk+1) might be equal to kn.
238 Thus we have to change all vertices if kn<=2*(kk+1) */
239
240
241 /* Make two steps one for the start and one for the end of the curve */
242
243
244 /* Make matrix for the kk first vertices */
245
246 for ( ki=kk-1, spek=smatrix; ki < 2*kk-1; ki++, spek+=kk )
247 {
248 /* we use kn instead of knstx since s1219 expects et[in-1] != et[in], we only
249 address vertices at the start so this does not matter */
250
251 s1219(stx, kk, kn, &kleft, st[ki], &kstat);
252 if ( kstat < 0 ) goto error;
253
254 s1701(ki, kleft, kk, knstx, &kpl, &kfi, &kla, st, stx, sp, salfa, &kstat);
255 if( kstat < 0 ) goto error;
256
257 /* Copy the discrete B-splines into the right position */
258
259 memcopy(spek + kfi, salfa + kpl + kfi, kla - kfi + 1, DOUBLE);
260 }
261
262
263
264 /* Do the factorisation of the matrix */
265
266 s6lufacp(smatrix, mpiv, kk, &kstat);
267 if ( kstat < 0 ) goto error;
268
269 /* The vertices in the curve are ordered in the sequence
270 (x1,y1,z1),..,(xi,yi,zi), i=1,..,in.
271 The only vertices affected by this backsubstitution is the kk first.
272 We want to treat the back substitution as idim(=3) backsubstitutions.
273 Thus we have to copy the proper parts of the vertices into a temporary
274 array. Do backsubstitution and copy back into the curve object */
275
276
277 for ( ki=0; ki < kdim; ki++ )
278 {
279 for ( kj=0,sfrom=(kreg->ecoef)+ki,sto=sb; kj < kk; kj++,sfrom+=kdim,sto++ )
280 *sto = *sfrom;
281
282 /* sb now contains the parts of vertices to be backsubsituted */
283
284 s6lusolp(smatrix, sb, mpiv, kk, &kstat);
285 if ( kstat < 0 ) goto error;
286
287 /* Copy the backsubsituted vertices back into scoef */
288
289 for ( kj=0,sto=scoef+ki,sfrom=sb; kj < kk; kj++,sfrom++,sto+=kdim )
290 *sto = *sfrom;
291 }
292
293
294 /* Make matrix for the kk last vertices */
295
296 for ( ki=0,spek=smatrix; ki < kk*kk; ki++,spek++ ) *spek = DZERO;
297
298
299 for ( ki=kn-kk,spek=smatrix; ki < kn; ki++,spek+=kk )
300 {
301 s1219(ext, kk, kn, &kleft, stx[ki], &kstat);
302 if ( kstat < 0 ) goto error;
303
304 s1701(ki, kleft, kk, kn, &kpl, &kfi, &kla, stx, ext, sp, salfa, &kstat);
305 if ( kstat < 0 ) goto error;
306
307 /* Copy the discrete B-splines into the right position */
308
309 memcopy(spek+kfi-(kn-kk), salfa+kpl+kfi, kla-kfi+1, DOUBLE);
310 }
311
312
313 /* Do the factorisation of the matrix */
314
315 s6lufacp(smatrix, mpiv, kk, &kstat);
316 if ( kstat < 0 ) goto error;
317
318 /* The vertices in the curve are ordered in the sequence
319 (x1,y1,z1),..,(xi,yi,zi), i=1,..,in1. The only vertices
320 affected by this backsubstitution is the kk-1 last rows.
321 We want to treat the back substitution as idim(=3) backsubstitutions.
322 Thus we have to copy the proper parts of the vertices into a temporary
323 array. Do backsubstitution and copy back into the surface object */
324
325 for ( ki=0; ki < kdim; ki++ )
326 {
327 for ( kj=0,sfrom=scoef+kdim*(kn-kk)+ki,sto=sb; kj < kk; kj++,sfrom+=kdim,sto++ )
328 *sto = *sfrom;
329
330 /* sb now contains the vertices to be backsubsituted */
331
332 s6lusolp(smatrix, sb, mpiv, kk, &kstat);
333 if ( kstat < 0 ) goto error;
334
335 /* Copy the backsubsituted vertices back into scoef */
336
337 for ( kj=0,sto=scoef+kdim*(kn-kk)+ki,sfrom=sb; kj < kk; kj++,sto+=kdim,sfrom++ )
338 *sto = *sfrom;
339 }
340
341
342 /* Copy knots and vertices into the surface object */
343
344 memcopy(kreg->ecoef, scoef, kdim*kn, DOUBLE);
345 memcopy(kreg->et, ext, kn+kk, DOUBLE);
346
347 /* Set periodicity flag */
348 kreg->cuopen = SISL_CRV_OPEN;
349
350
351 /* Task done */
352
353 *rc = kreg;
354
355 *jstat = 0;
356 goto out;
357
358
359 /* Error in allocation. */
360
361 err101:
362 *jstat = -101;
363 s6err("s1233",*jstat,kpos);
364 goto out;
365
366
367 /* Error in curve description - order less than 1 */
368
369 err110:
370 *jstat = -110;
371 s6err("s1233",*jstat,kpos);
372 goto out;
373
374
375 /* Error in curve desctiption - number of vertices is less than the order */
376
377 err111:
378 *jstat = -111;
379 s6err("s1233",*jstat,kpos);
380 goto out;
381
382
383 /* Error in lower level routine. */
384
385 error:
386 *jstat = kstat;
387 s6err("s1233",*jstat,kpos);
388 goto out;
389
390
391 out:
392
393 /* Free allocated scratch */
394
395 if (salloc != SISL_NULL) freearray(salloc);
396 if (mpiv != SISL_NULL) freearray(mpiv);
397
398 return;
399
400 }
401