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: s1333cycli.c,v 1.4 2005-02-28 09:04:48 afr Exp $
45 *
46 */
47
48
49 #define S1333_CYCLIC
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1333_cyclic(SISLSurf * vsurf,int icont,int * jstat)55 s1333_cyclic(SISLSurf *vsurf,int icont,int *jstat)
56 #else
57 void s1333_cyclic(vsurf,icont,jstat)
58 SISLSurf *vsurf;
59 int icont;
60 int *jstat;
61 #endif
62 /*
63 *********************************************************************
64 *
65 * PURPOSE : To describe the B-spline surface (i.e. not NURBS) vsurf with a
66 * cyclic basis of continuity icont in the first parameter direction.
67 *
68 * INPUT : vsurf - Pointer to the surface
69 * icont - The required continuity
70 *
71 * OUTPUT : jstat - status messages
72 * > 0 : warning
73 * = 0 : ok
74 * < 0 : error
75 * vsurf - The surface with modified description
76 *
77 * METHOD : 1. The cyclic knot vector with the right continuity is made
78 * 2. The transformation matrix for the ik2 first vertices between
79 * the cyclic and the old knot vector in second paramter direction
80 * is made.
81 * 3. The transformation matrix is inverted and used to update the
82 * ik2 first vertices.
83 * 4. The transformation matrix for the ik2 last vertices between
84 * the cyclic and the old knot vector in second paramter direction
85 * is made.
86 * 5. The transformation matrix is inverted and used to update the
87 * ik2 last vertices.
88 * 6. The knot vector is updated.
89 *-
90 * CALLS : s1701,s6err.
91 *
92 * Written by : Tor Dokken, SI, Oslo, Norway,feb-1992
93 *
94 *********************************************************************
95 */
96 {
97 double *scycl=SISL_NULL; /* Cyclic version of knot vector */
98 double *smatrix=SISL_NULL; /* Matrix converting between baes */
99 double *salloc=SISL_NULL; /* Matrix for memory allocation */
100 double *salfa=SISL_NULL; /* The values of a discrete B-spline
101 calculation */
102 double *spek=SISL_NULL; /* Pointer used in traversing arrays */
103 double *scoef=SISL_NULL; /* Copy of the vertices of the surface */
104 double *sb=SISL_NULL; /* Right hand side of equation */
105 double *sfrom,*sto;
106 double *sp; /* Hlep array for s1701 */
107 double *st1=SISL_NULL; /* Internal version of et1 */
108 double *stx=SISL_NULL; /* Knot vector after insertion of knots
109 at start */
110
111 int kdim = vsurf->idim;
112 int rat = (vsurf->ikind == 1 || vsurf->ikind == 3) ? 0 : 1;
113 double *sourcecoef = (rat) ? vsurf->rcoef : vsurf->ecoef;
114 int kdim2 = kdim + rat;
115 int kk1 = vsurf->ik1;
116 int kn1 = vsurf->in1;
117 int kn2 = vsurf->in2;
118 int knumb = kdim2*kn1;
119
120 int kcont; /* Checked continuity */
121 int kmult; /* Multiplicity of knot kk2-1 and kn2*/
122 int ki,kj,kl;
123 int kperlength;
124 int kant;
125 int kleft=0; /* Pointer into knot vector */
126 int kpl,kfi,kla; /* Pointers into conversion matrix */
127 int kstat;
128 int *mpiv=SISL_NULL; /* Pointer to pivotation array */
129 int kpos = 0;
130 int knst1; /* NUmber of basis functions in st1 */
131 int knstx; /* Number of basis functions in stx */
132
133
134
135
136 /* Test continuity */
137
138 if (icont < 0) goto finished;
139 kcont = icont;
140 if (icont >= kk1-2) icont = kk1-2;
141
142 /* Make multiplicty to be used at value et1[ik1-1] and et1[in1] */
143 kmult = kk1 - kcont - 1;
144
145 /* Make the number of knots to be changed at the start and the end, this
146 is also equal to extra knot to be inserted in internal version of array et */
147
148 kant = kk1-kmult;
149
150
151 /* Alloocate array for pivotation vector */
152
153 mpiv = new0array(2*kk1,INT);
154 if (mpiv == SISL_NULL) goto err101;
155
156 salloc = new0array(3*kn1+9*kk1+4*kk1*kk1+kdim2*kn1*kn2,DOUBLE);
157 if (salloc == SISL_NULL) goto err101;
158 scycl = salloc; /* Size kn1+kk1 */
159 smatrix = scycl + kn1 + kk1; /* Max size 4*kk1*kk1 */
160 salfa = smatrix + 4*kk1*kk1; /* Size kk1 */
161 scoef = salfa + kk1; /* Size kdim2*kn1*kn2 */
162 sb = scoef + kdim2*kn1*kn2; /* Size 2*kk1 */
163 sp = sb + 2*kk1; /* Size kk1 */
164 st1 = sp + kk1; /* Size kn1 + 2*kk1 */
165 stx = st1 + kn1 + 2*kk1; /* Size kn1 + 2*kk1 */
166
167
168
169 /* Copy vertices, to avoid destruction of surface */
170
171 memcopy(scoef,sourcecoef,kdim2*kn1*kn2,DOUBLE);
172
173
174
175 /* Make cyclic knot vector */
176
177
178 /* First copy all knots */
179
180 memcopy(scycl,vsurf->et1,kn1+kk1,DOUBLE);
181
182 /* The change the ik1 first and the ik1 last knots to make a cyclic basis */
183
184 kperlength = kn1 - kk1 + kmult;
185
186 /* Make knots 0 to ik - kmult - 1 */
187
188 for (ki=kk1-kmult-1 ; 0<=ki ; ki--)
189 {
190 scycl[ki] = scycl[kk1-1] - (scycl[kn1] - scycl[ki+kperlength]);
191 }
192
193
194 /* Make knots kn1 + kmult to kn1 + kk1 -1 */
195
196 for (ki=kmult ; ki < kk1 ; ki++)
197 {
198 scycl[kn1+ki] = scycl[kn1] + (scycl[kk1+ki-kmult] - scycl[kk1-1]);
199
200 }
201 /* s1701 expects et1 to be a refinement of scyclic, thus we have to make a new
202 version of et1 with the extra kk1-kmult new knots before the start and
203 after the end and one intermediate version with only kk1-kmult at the start */
204
205 memcopy(st1,scycl,kant,DOUBLE);
206 memcopy(st1+kant,vsurf->et1,kn1+kk1,DOUBLE);
207 memcopy(st1+kant+kk1+kn1,scycl+kn1+kk1-kant,kant,DOUBLE);
208 knst1 = kn1 + 2*kant;
209
210 memcopy(stx,scycl,kn1,DOUBLE);
211 memcopy(stx+kn1,st1+kant+kn1,kk1+kant,DOUBLE);
212 knstx = kn1 + kant;
213
214 /* STEP 2 Make matrix going between bases, only the kk2-kmult first and last knots
215 are to be changed. */
216
217
218 /* Now we have two cases. We know that only the kk1-kmult first and kk1-kmult
219 last vertices are to be changed. However 2*(kk1-kmult) might be a bigger
220 number than kn1. Thus we have to change all vertices if kn1<=2(kk1-kmult) */
221
222
223 /* Make two steps one for the start and one for the end of the surface */
224
225
226 /* Make matrix for the kk1 first vertices */
227
228 for (ki=kant,spek=smatrix ; ki <kk1+kant ; ki++, spek+=kk1)
229 {
230 /* we use kn1 instead of knstx since s1219 expects et[in-1] != et[in], we only
231 address vertices at the start so this does not matter*/
232
233 s1219(stx,kk1,kn1,&kleft,st1[ki],&kstat);
234 if (kstat<0) goto error;
235
236 s1701(ki,kleft,kk1,knstx,&kpl,&kfi,&kla,st1,stx,sp,salfa,&kstat);
237 if(kstat<0) goto error;
238
239 /* Copy the discrete B-splines into the right position */
240
241 memcopy(spek+kfi,salfa+kpl+kfi,kla-kfi+1,DOUBLE);
242 }
243
244
245
246 /* Do the factorisation of the matrix */
247
248 s6lufacp(smatrix,mpiv,kk1,&kstat);
249 if (kstat<0) goto error;
250
251 /* The vertices in the surface are ordered in the sequence
252 (x11,y11,z11),..,(xij,yij,zij), i=1,..,in1 and j=1,..,in2.
253 i is running faster than j. The only vertices
254 affected by this backsubstitution is the kant first rows.
255 We want to treat the back substitution
256 as idim(=3)*in2 backsubstitutions. Thus we have to copy the proper
257 parts of the vertices into a temporary array. Do backsubstitution and
258 copy back into the surface object */
259
260
261 for (ki=0 ; ki<kn2 ; ki++)
262 for (kl=0 ; kl<kdim2 ; kl++)
263 {
264 for (kj=0, sfrom=sourcecoef+ki*knumb+kl,sto=sb ;
265 kj<kk1 ; kj++,sfrom+=kdim2,sto++)
266 *sto = *sfrom;
267
268 /* sb now contains the vertices to be backsubsituted */
269
270 s6lusolp(smatrix,sb,mpiv,kk1,&kstat);
271 if (kstat<0) goto error;
272
273 /* Copy the backsubsituted vertices back into scoef */
274
275 for (kj=0, sto=scoef+ki*knumb+kl,sfrom=sb ;
276 kj<kk1 ; kj++,sfrom++,sto+=kdim2)
277 *sto = *sfrom;
278 }
279
280
281 /* Make matrix for the kk1 last vertices */
282
283
284 for (ki=0,spek=smatrix ; ki<kk1*kk1 ; ki++,spek++) *spek = DZERO;
285
286
287 for (ki=kn1-kk1 ,spek=smatrix ; ki <kn1 ; ki++, spek+=kk1)
288 {
289 s1219(scycl,kk1,kn1,&kleft,stx[ki],&kstat);
290 if (kstat<0) goto error;
291
292 s1701(ki,kleft,kk1,kn1,&kpl,&kfi,&kla,stx,scycl,sp,salfa,&kstat);
293 if(kstat<0) goto error;
294
295 /* Copy the discrete B-splines into the right position */
296
297 memcopy(spek+kfi-(kn1-kk1),salfa+kpl+kfi,kla-kfi+1,DOUBLE);
298 }
299
300
301
302 /* Do the factorisation of the matrix */
303
304 s6lufacp(smatrix,mpiv,kk1,&kstat);
305 if (kstat<0) goto error;
306
307 /* The vertices in the surface are ordered in the sequence
308 (x11,y11,z11),..,(xij,yij,zij), i=1,..,in1 and j=1,..,in2.
309 i is running faster than j The only vertices
310 affected by this backsubstitution is the kant last rows.
311 We want to treat the back substitution
312 as idim(=3)*in2 backsubstitutions. Thus we have to copy the proper
313 parts of the vertices into a temporary array. Do backsubstitution and
314 copy back into the surface object */
315
316 for (ki=0 ; ki<kn2 ; ki++)
317 for (kl=0 ; kl<kdim2 ; kl++)
318 {
319 for (kj=0, sfrom=scoef+ki*knumb+kdim2*(kn1-kk1)+kl,sto=sb ;
320 kj<kk1 ; kj++,sfrom+=kdim2,sto++)
321 *sto = *sfrom;
322
323 /* sb now contains the vertices to be backsubsituted */
324
325 s6lusolp(smatrix,sb,mpiv,kk1,&kstat);
326 if (kstat<0) goto error;
327
328 /* Copy the backsubsituted vertices back into scoef */
329
330 for (kj=0, sto=scoef+ki*knumb+kdim2*(kn1-kk1)+kl,sfrom=sb ;
331 kj<kk1 ; kj++,sto+=kdim2,sfrom++)
332 *sto = *sfrom;
333 }
334
335
336 /* Copy knots and vertices into the surface object */
337
338 memcopy(sourcecoef,scoef,kdim2*kn1*kn2,DOUBLE);
339 memcopy(vsurf->et1,scycl,kn1+kk1,DOUBLE);
340
341 /* Set periodicity flag */
342 vsurf->cuopen_1 = SISL_SURF_PERIODIC;
343
344 /* Update divided coefficients */
345 if (rat)
346 {
347 for (ki=0; ki<kn1*kn2; ++ki)
348 {
349 for (kj=0; kj<kdim; ++kj)
350 vsurf->ecoef[ki*kdim+kj] =
351 vsurf->rcoef[ki*kdim2+kj]/vsurf->rcoef[ki*kdim2+kdim];
352 }
353 }
354
355 /* Task done */
356
357 finished:
358
359 *jstat = 0;
360 goto out;
361
362 /* Error in allocation. */
363
364 err101:
365 *jstat = -101;
366 s6err("s1333_cyclic",*jstat,kpos);
367 goto out;
368
369
370
371 /* Error in lower level routine. */
372
373 error :
374 *jstat = kstat;
375 s6err("s1333_cyclic",*jstat,kpos);
376 goto out;
377 out:
378
379 /* Free allocated scratch */
380 if (salloc != SISL_NULL) freearray(salloc);
381 if (mpiv != SISL_NULL) freearray(mpiv);
382
383 return;
384
385 }
386