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: s1018.c,v 1.3 2005-02-28 09:04:48 afr Exp $
45 *
46 */
47
48
49 #define S1018
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
s1018(SISLCurve * pc,double epar[],int inpar,SISLCurve ** rcnew,int * jstat)54 void s1018 (SISLCurve *pc, double epar[], int inpar,
55 SISLCurve **rcnew, int *jstat)
56 #else
57 void
58 s1018 (pc, epar, inpar, rcnew, jstat)
59 SISLCurve *pc;
60 double epar[];
61 int inpar;
62 SISLCurve **rcnew;
63 int *jstat;
64 #endif
65 /*
66 *********************************************************************
67 *
68 *********************************************************************
69 *
70 * PURPOSE : Insert a given set of knots into the description of
71 * a B-spline curve.
72 * NOTE : When the curve is periodic, the input parameter values
73 * must lie in the HALFOPEN [et[kk-1], et[kn), the function
74 * will automatically update the extra knots and
75 * coeffisients.
76 * rcnew->in is still eq to pc->in + inpar !
77 *
78 * INPUT : pc - SISLCurve to be refined.
79 * epar - Parameter values of knots to be inserted.
80 * The values are stored in increasing order.
81 * inpar - Number of parameter values in epar.
82 *
83 *
84 *
85 * OUTPUT : rcnew - The new, refined curve .
86 * jstat - status messages
87 * > 0 : warning
88 * = 0 : ok
89 * < 0 : error
90 *
91 *
92 * METHOD :
93 *
94 *
95 * REFERENCES :
96 *
97 *-
98 * CALLS : newCurve - Allocate space for a new curve-object.
99 * freeCurve - Free space occupied by given curve-object.
100 * S1701.C - Making the knot-inserten-transformation matrix.
101 *
102 * WRITTEN BY : Arne Laksaa, SI, 88-11.
103 * CHANGED BY : Ulf J. Krystad, SI, 92-01
104 * Treatment of periodic curves.
105 **********************************************************************/
106 {
107 register int ki, ki2; /* Control variable in loop. */
108 register int kj, kj1, kj2; /* Control variable in loop. */
109 register double *s1; /* Pointers used in loop. */
110 int kstat; /* Local status variable. */
111 int k1, k2, k3; /* Index to array of parameter values. */
112 int kmy; /* An index to the knot-vector. */
113 int kpl, kfi, kla; /* To posisjon elements in trans.-matrix. */
114 int kk = pc->ik; /* Order of the input curve. */
115 int kn = pc->in; /* Number of the vertices in input curves. */
116 int kdim = pc->idim; /* Dimensjon of the space in whice curve lies.*/
117 int kn1; /* Number of vertices in the new curve . */
118 double tstart, tend; /* Endparameters of curve. */
119 double tpar; /* Parameter value of knot to insert. */
120 double *st = SISL_NULL; /* The new knot-vector. */
121 double *sp = SISL_NULL; /* To use in s1701.c */
122 double *salfa = SISL_NULL; /* A line of the trans.-matrix. */
123 double *scoef = SISL_NULL; /* The new vertices. */
124 double *coef = SISL_NULL; /* The old vertices. */
125 SISLCurve *q1 = SISL_NULL; /* Pointer to new curve-object. */
126 /* Periodicity treatment -------------------------- */
127 double *st2 = SISL_NULL; /* The new knot-vector. */
128 double *scoef2 = SISL_NULL; /* The new vertices. */
129 double t1; /* Start of knot vector */
130 double t2; /* Start of full basis part of knot vector */
131 double t3; /* End of full basis part of knot vector */
132 double t4; /* End of knot vector */
133 double tmod; /* t3 - t2, period size */
134 double mod_neg; /* parametervalue shifted tmod left */
135 double mod_pos; /* parametervalue shifted tmod right */
136 int no_neg, no_pos; /* Nmb of Xtra parval. to insert (lft,rght)*/
137 double *negpar = SISL_NULL; /* Array of parameter values shift lft */
138 double *pospar = SISL_NULL; /* Array of parameter values shift rght */
139 double *pararr = SISL_NULL; /* Pointer to parameter array treated */
140 double *periodarr = SISL_NULL; /* Total parameter array, periodic case*/
141 /* --------------------------------------------------------- */
142
143 /* Check that we have a curve to treat. */
144
145 if (!pc)
146 goto err150;
147 if (pc->ikind == 2 || pc->ikind == 4)
148 {
149 kdim++;
150 coef = pc->rcoef;
151 }
152 else
153 {
154 coef = pc->ecoef;
155 }
156
157
158 /* Periodicity treatment -------------------------- */
159 if (pc->cuopen == SISL_CRV_PERIODIC)
160 {
161 /* Get values in knotvector for start, start of full basis,
162 end of full basis, end. */
163 t1 = *(pc->et);
164 t2 = *(pc->et + kk - 1);
165 t3 = *(pc->et + kn);
166 t4 = *(pc->et + kn + kk - 1);
167 tmod = t3 - t2;
168
169 /* Check that the input parameter values lie in the legal parameter
170 interval of the curve, ie in the HALFOPEN [et[kk-1], et[kn) */
171 for (k1 = 0; k1 < inpar; k1++)
172 if (epar[k1] < t2 ||epar[k1] >= t3)
173 goto err158;
174
175 no_neg = 0;
176 no_pos = 0;
177 if ((negpar = newarray (inpar, double)) == SISL_NULL)
178 goto err101;
179 if ((pospar = newarray (inpar, double)) == SISL_NULL)
180 goto err101;
181
182
183 for (k1 = 0; k1 < inpar; k1++)
184 {
185 mod_neg = epar[k1] - tmod;
186 mod_pos = epar[k1] + tmod;
187 if (mod_neg <= t2 &&
188 mod_neg > t1)
189 {
190 negpar[no_neg] = mod_neg;
191 no_neg++;
192 }
193 if (mod_pos >= t3 &&
194 mod_pos <t4)
195 {
196 pospar[no_pos] = mod_pos;
197 no_pos++;
198 }
199 }
200
201 if ((periodarr = newarray (inpar + no_neg + no_pos, double)) == SISL_NULL)
202 goto err101;
203 memcopy (periodarr, negpar, no_neg, double);
204 memcopy (periodarr + no_neg, epar, inpar, double);
205 memcopy (periodarr + no_neg + inpar, pospar, no_pos, double);
206 inpar += no_neg + no_pos;
207
208 pararr = periodarr;
209 }
210 else
211 pararr = epar;
212
213
214
215 /* Check that the input parameter values lie in the parameter
216 interval of the curve. */
217
218 tstart = *(pc->et);
219 tend = *(pc->et + kn + kk - 1);
220 for (k1 = 0; k1 < inpar; k1++)
221 if (pararr[k1] < tstart || pararr[k1] > tend)
222 goto err158;
223
224 /* Allocate space for the kk elements which may not be zero in eache
225 line of the basic transformation matrix, and space for new knots
226 to use in s1701.c */
227
228 if ((salfa = newarray (kk, double)) == SISL_NULL)
229 goto err101;
230 if ((sp = newarray (kk, double)) == SISL_NULL)
231 goto err101;
232
233 /* Find the number of vertices in the new curve. */
234
235 kn1 = kn + inpar;
236
237 /* Allocating the new arrays to the new curve. */
238
239 if ((st = newarray (kn1 + kk, double)) == SISL_NULL)
240 goto err101;
241 if ((scoef = new0array (kn1 * kdim, double)) == SISL_NULL)
242 goto err101;
243
244 /* Making the new knotvectors. */
245
246 for (k2 = 0, k3 = 0, k1 = 0; k1 < inpar; k1++)
247 {
248 tpar = pararr[k1];
249 for (; k2 < kn + kk; k2++)
250 {
251 if (pc->et[k2] <= tpar)
252 st[k3++] = pc->et[k2];
253 else
254 break;
255 }
256 st[k3++] = tpar;
257 }
258 for (; k2 < kn + kk; k2++)
259 st[k3++] = pc->et[k2];
260
261 /* Updating the coefisientvector to the new curve.*/
262
263 for (s1 = scoef, ki = 0, kmy = 0; ki < kn1; ki++)
264 {
265 /* Here we compute a new line with line number ki of
266 the knot inserten matrix. */
267
268 while (pc->et[kmy + 1] <= st[ki])
269 kmy++;
270 s1701 (ki, kmy, kk, kn, &kpl, &kfi, &kla, st, pc->et, sp, salfa, &kstat);
271 if (kstat)
272 goto err153;
273
274 /* Compute the kdim vertices with the same "index". */
275
276 for (kj = 0; kj < kdim; kj++, s1++)
277 for (*s1 = 0, kj1 = kfi, kj2 = kfi + kpl; kj1 <= kla; kj1++, kj2++)
278 {
279 ki2 = kj1 * kdim + kj;
280 *s1 += salfa[kj2] * coef[ki2];
281 }
282 }
283
284 /* Periodicity treatment -------------------------- */
285
286 if (pc->cuopen == SISL_CRV_PERIODIC)
287 {
288 kn1 -= no_pos + no_neg;
289 /* Allocating the new arrays to the new curve. */
290
291 if ((st2 = newarray (kn1 + kk, double)) == SISL_NULL)
292 goto err101;
293 if ((scoef2 = new0array (kn1 * kdim, double)) == SISL_NULL)
294 goto err101;
295 memcopy (st2, st + no_neg, kn1 + kk, double);
296 memcopy (scoef2, scoef + no_neg * kdim, kn1 * kdim, double);
297
298 if (st)
299 freearray (st);
300 if (scoef)
301 freearray (scoef);
302
303 /* Allocating new curve-objects.*/
304
305 if ((q1 = newCurve(kn1,kk,st2,scoef2,pc->ikind,pc->idim,2)) == SISL_NULL)
306 goto err101;
307 }
308 else
309 {
310 /* Allocating new curve-objects.*/
311 if ((q1 = newCurve (kn1, kk, st, scoef, pc->ikind, pc->idim, 2)) == SISL_NULL)
312 goto err101;
313 }
314
315 q1->cuopen = pc->cuopen;
316
317 /* Updating output. */
318
319 *rcnew = q1;
320 *jstat = 0;
321 goto out;
322
323
324 /* Error. Subrutine error. */
325
326 err153:*jstat = kstat;
327 goto outfree;
328
329
330 /* Error. No curve to treat. */
331
332 err150:*jstat = -150;
333 goto out;
334
335 /* Error. Parameter value to insert is outside the curve. */
336
337 err158:*jstat = -158;
338 goto out;
339
340 /* Error. Allocation error, not enough memory. */
341
342 err101:*jstat = -101;
343 goto outfree;
344
345
346 outfree:
347 if (q1)
348 freeCurve (q1);
349 else
350 {
351 if (st)
352 freearray (st);
353 if (scoef)
354 freearray (scoef);
355 if (st2)
356 freearray (st2);
357 if (scoef2)
358 freearray (scoef2);
359 }
360
361
362 /* Free local used memory. */
363
364 out:if (salfa)
365 freearray (salfa);
366 if (sp)
367 freearray (sp);
368 if (periodarr)
369 freearray (periodarr);
370 if (pospar)
371 freearray (pospar);
372 if (negpar)
373 freearray (negpar);
374
375 }
376