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: s1017.c,v 1.4 2005-02-28 09:04:48 afr Exp $
45 *
46 */
47
48
49 #define S1017
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1017(SISLCurve * pc,SISLCurve ** rc,double apar,int * jstat)55 s1017 (SISLCurve * pc, SISLCurve ** rc, double apar, int *jstat)
56 #else
57 void
58 s1017 (pc, rc, apar, jstat)
59 int *jstat;
60 double apar;
61 SISLCurve *pc, **rc;
62 #endif
63 /*
64 *********************************************************************
65 *
66 *********************************************************************
67 *
68 * PURPOSE : Insert a given knot into the description of
69 * a B-spline curve.
70 * NOTE : When the curve is periodic, the input parameter value
71 * must lie in the HALFOPEN [et[kk-1], et[kn), the function
72 * will automatically update the extra knots and
73 * coeffisients.
74 * rcnew->in is still eq to pc->in + 1!
75 *
76 *
77 * INPUT : pc - SISLCurve to be refined.
78 * apar - Parameter values of knot to be s1017ed.
79 *
80 *
81 *
82 * OUTPUT : rc - The new, refined curve.
83 * jstat - status messages
84 * > 0 : warning
85 * = 0 : ok
86 * < 0 : error
87 *
88 *
89 * METHOD :
90 *
91 *
92 * REFERENCES :
93 *
94 *-
95 * CALLS : newCurve - Allocate space for a new curve-object.
96 * freeCurve - Free space occupied by given curve-object.
97 * S1701.C - Making the knot-s1017en-transformation matrix.
98 * s1017knots in periodic case.
99 * WRITTEN BY : Arne Laksaa, SI, 88-11.
100 * CHANGED BY : Ulf J. Krystad, SI, 92-01
101 * Treatment of periodic curves.
102 * Revised by : Paal Fugelli, SINTEF, Oslo, Norway, Nov. 1994. Undefined
103 * symbol 's1017knots()' changed to 's1018()'. Closed ('cuopen'
104 * flag=0) and 'ikind' flag are now handled correctly. Fixed memory
105 * problem for rationals. Cleaned up a bit.
106 *
107 **********************************************************************/
108 {
109
110 int kstat; /* Local status variable. */
111 int kmy; /* An index to the knot-vector. */
112 int kpl, kfi, kla; /* To posisjon elements in trans.-matrix. */
113 int kk = pc->ik; /* Order of the input curve. */
114 int kn = pc->in; /* Number of the vertices in input curve. */
115 int kdim = pc->idim; /* Dimension of the space in whice curve lies.*/
116 int kn1; /* Number of vertices in the new curve. */
117 int kch; /* First vertice to be changes. */
118 int knum; /* Number of knots less and equal than
119 the intersection point. */
120 int ki; /* Control variable in loop. */
121 int kj, kj1, kj2; /* Control variable in loop. */
122 double *s1; /* Pointers used in loop. */
123 double *st = SISL_NULL; /* The first new knot-vector. */
124 double *salfa = SISL_NULL; /* A line of the trans.-matrix. */
125 double *scoef = SISL_NULL; /* The first new vertice. */
126 double *sp; /* Help array for use in s1701. */
127 double *ecoef; /* Pointer to input curve ecoef/rcoef vector. */
128 SISLCurve *qc = SISL_NULL; /* Pointer to new curve-object. */
129
130
131
132 /* Make sure the returned pointer is valid. */
133
134 *rc = SISL_NULL;
135
136
137 /* Check that we have a curve. */
138
139 if (!pc)
140 goto err150;
141
142
143 /* Find the number of vertices in the new curve. */
144
145 kn1 = kn + 1;
146
147 if ( kn1 <= 0 )
148 goto err150;
149
150
151
152 /* Periodicity treatment -------------------------- */
153 if (pc->cuopen == SISL_CRV_PERIODIC)
154 {
155 s1018(pc, &apar, 1, rc, &kstat);
156 if (kstat < 0)
157 goto err153;
158 goto out;
159 }
160
161
162
163 /* Check that the intersection point is an interior point. */
164
165 if (apar < *(pc->et) || apar > *(pc->et + kn + kk - 1))
166 goto err158;
167
168
169 /* Check if the curve is rational. */
170
171 if (pc->ikind == 2 || pc->ikind == 4)
172 {
173 ecoef = pc->rcoef;
174 kdim++;
175 }
176 else
177 ecoef = pc->ecoef;
178
179
180 /* Allocate space for the kk elements which may not be zero in eache
181 line of the basic transformation matrix and for a help array of
182 kk elements.*/
183
184 if ((salfa = newarray (2 * kk, double)) == SISL_NULL)
185 goto err101;
186
187 sp = salfa + kk;
188
189 /* Find the number of the knots which is smaller or like
190 the intersection point.*/
191
192 s1 = pc->et;
193
194 if (apar > pc->et[0] && apar < pc->et[kn + kk - 1])
195 {
196 /* Using binary search*/
197 kj1 = 0;
198 kj2 = kk + kn - 1;
199 knum = (kj1 + kj2) / 2;
200 while (knum != kj1)
201 {
202 if (s1[knum] < apar)
203 kj1 = knum;
204 else
205 kj2 = knum;
206 knum = (kj1 + kj2) / 2;
207 }
208 knum++; /* The smaller knots. */
209
210 while (s1[knum] == apar)
211 /* The knots that are equal to the intersection point. */
212 knum++;
213 }
214 else if (apar == pc->et[0])
215 {
216 knum = 0;
217 while (s1[knum] == apar)
218 /* The knots that are equal to the intersection point. */
219 knum++;
220 }
221 else if (apar == pc->et[kn + kk - 1])
222 {
223 knum = 0;
224 while (s1[knum] < apar)
225 /* The knots that are less than or equal to the intersection point. */
226 knum++;
227 }
228
229
230
231 /* Allocating the new arrays to the new curve. */
232
233 if ((scoef = newarray (kn1 * kdim, double)) == SISL_NULL)
234 goto err101;
235 if ((st = newarray (kn1 + kk, double)) == SISL_NULL)
236 goto err101;
237
238
239 /* Copying the knotvectors, all but the intersection point from
240 the old curve to the new curves */
241
242 memcopy (st, pc->et, knum, double);
243 st[knum] = apar;
244 if (knum < kn + kk)
245 memcopy (st + knum + 1, pc->et + knum, kn + kk - knum, double);
246
247
248 /* Copying the coefisientvector to the new curve. (Here 'ecoef' points
249 to 'pc->rcoef' or 'pc->ecoef' depening on if 'pc' is rational or not). */
250
251 kch = knum - kk + 1;
252
253 if (kch > 0)
254 memcopy (scoef, ecoef, kdim * kch, double);
255 if (knum < kn1)
256 memcopy (scoef + kdim * knum, ecoef + kdim * (knum - 1),
257 kdim * (kn1 - knum), double);
258
259
260 /* Updating the coefisientvectors the new curve.*/
261
262 /* Updating the first curve. */
263
264 for (ki = max (0, kch), kmy = 0, s1 = scoef + ki * kdim; ki < min (knum + 1, kn1); ki++)
265 {
266 /* Initialising:
267 ki = kch, Index of the vertices we are going to
268 change. Starting with kch, but if
269 kch is negativ we start at zero.
270 s1=scoef1+ki*kdim,Pointer at the first vertice to
271 change. */
272
273
274 /* Using the Oslo-algorithm to make a transformation-vector
275 from the old vertices to one new vertice. */
276
277 while (kmy < kn + kk && pc->et[kmy] <= st[ki])
278 kmy++;
279
280 s1701 (ki, kmy - 1, kk, kn, &kpl, &kfi, &kla, st, pc->et, sp, salfa, &kstat);
281 if (kstat)
282 goto err153;
283
284
285 /* Compute the kdim vertices with the same "index". */
286
287 for (kj = 0; kj < kdim; kj++, s1++)
288 for (*s1 = 0, kj1 = kfi, kj2 = kfi + kpl; kj1 <= kla; kj1++, kj2++)
289 *s1 += salfa[kj2] * ecoef[kj1 * kdim + kj];
290 }
291
292
293
294 /* Allocating new curve-objects.*/
295
296
297 if ((qc = newCurve (kn1, kk, st, scoef, pc->ikind, pc->idim, 2)) == SISL_NULL)
298 goto err101;
299
300
301 /* Updating output. */
302
303 qc->cuopen = pc->cuopen; /* Ok since only internal knots can be inserted. */
304
305 *rc = qc;
306
307 *jstat = 0;
308 goto out;
309
310
311 /* Error. Error in low level routine. */
312
313 err153:
314 *jstat = kstat;
315 goto outfree;
316
317
318 /* Error. No curve to s1017 a new knot. */
319
320 err150:
321 *jstat = -150;
322 goto out;
323
324
325 /* Error. The intersection-point is outside the curve. */
326
327 err158:
328 *jstat = -158;
329 goto out;
330
331
332 /* Error. Allocation error, not enough memory. */
333
334 err101:
335 *jstat = -101;
336 goto outfree;
337
338
339 outfree:
340 if (qc)
341 freeCurve (qc);
342 else
343 {
344 if (st)
345 freearray (st);
346 if (scoef)
347 freearray (scoef);
348 }
349
350 /* Free local used memory. */
351
352 out:if (salfa)
353 freearray (salfa);
354 }
355