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: s1894.c,v 1.3 2005-02-28 09:04:49 afr Exp $
45 *
46 */
47
48
49 #define S1894
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1894(double oknots[],int oik,int oin,int der1,int der2,double earray[],int dimp1,int narr,double * nknots[],int * nik,int * nin,int * jstat)55 s1894 (double oknots[], int oik, int oin, int der1, int der2, double earray[],
56 int dimp1, int narr, double *nknots[], int *nik, int *nin, int *jstat)
57 #else
58 void
59 s1894 (oknots, oik, oin, der1, der2, earray, dimp1, narr, nknots,
60 nik, nin, jstat)
61 double oknots[];
62 int oik;
63 int oin;
64 int der1;
65 int der2;
66 double earray[];
67 int dimp1;
68 int narr;
69 double *nknots[];
70 int *nik;
71 int *nin;
72 int *jstat;
73 #endif
74 /*
75 *********************************************************************
76 *
77 *********************************************************************
78 *
79 * PURPOSE : To produce a knot vector for the dot product of two derivates
80 * of a B-spline curve having the input knot vector.
81 *
82 * INPUT : oknots - The original knot vector.
83 * oik - The order of the original knot vector.
84 * oin - The number of degrees of freedom in the original
85 * knot vector.
86 * der1 - The first of the derivatives.
87 * der2 - The secondof the derivatives.
88 * earray - Description of the array to be used.
89 * dimp1 - The dimension of the first matrix plane.
90 * narr - Number of parallel matrix planes.
91 *
92 * OUTPUT : nknots - The new knot vector.
93 * nik - The order of the new knot vector.
94 * nin - The number of degrees of freedom in the new
95 * knot vector.
96 * jstat - Status variable:
97 * > 0 : warning
98 * = 0 : ok
99 * < 0 : error
100 *
101 * METHOD : The order nk of the new new basis is determined.
102 * The multiplicity of the knots are counted and the right number
103 * of knots inserted according to the order of the original basis
104 * and the derivates involved. At et[ik] and et[in+1] nk knots
105 * are inserted.
106 *
107 * REFERENCES : Fortran version:
108 * Tor Dokken, SI, 1982-10
109 *
110 * CALLS : s6err.
111 *
112 * WRITTEN BY : Trond Vidar Stensby, SI, 1991-06
113 * REVISED BY : Johannes Kaasa, SI, May 92. (Taking the number of derivatives
114 * into account when calculating the multiplicity of the interior
115 * knots, to reduce the continuity in accordance with the
116 * derivatives. I also changed minimum allowed order from 1 to 2,
117 * to avoid errors in s1890).
118 *
119 **************************************************************** */
120 {
121 int size; /* The total size of earray. */
122 int mult; /* Multiplicity of knots */
123 int numb; /* Number of new knots. */
124 int kdim; /* dimp1 -1 (sub-matrix dimension) */
125 int empty; /* Used to check if sub-matrix of earray
126 is zero. */
127 int kl; /* Loop control varibles. */
128 int count1;
129 int count2;
130 int count3;
131 int start;
132 int stop;
133
134 double eps; /* Resolution. */
135 double maximum; /* The maximum value in et. */
136 double prev; /* Knot value. (extracted from orig) */
137 double curr; /* Knot value. (extracted from orig) */
138 int kpos = 0;
139 int der = max(der1, der2);
140
141 *jstat = 0;
142
143
144 /* Test if legal input. */
145
146 if (oik <= 1 || oin < oik)
147 goto err112;
148
149
150 /* Test if knot vector degenerate. */
151
152 if (oknots[oik - 1] >= oknots[oin])
153 goto err112;
154
155
156 /* The maximal number of knots to be produced at a specified knot value
157 * is the order of the B-spline basis produced. */
158
159 /* Allocate space for new knot vector */
160
161 (*nknots) = newarray ((oin + oik) * oik, DOUBLE);
162 if (*nknots == SISL_NULL)
163 goto err101;
164
165
166 /* Check if sub-matrix is zero. */
167
168 kdim = dimp1 - 1;
169 size = dimp1 * dimp1;
170 empty = TRUE;
171
172 for (count1 = 0; count1 < narr && empty; count1++)
173 for (count2 = 0; count2 < kdim && empty; count2++)
174 for (count3 = 0; count3 < kdim && empty; count3++)
175 if (earray[count1 * size + count2 * dimp1 + count3] != (double)0.)
176 empty = FALSE;
177
178
179 /* Assign value to nk. */
180
181 if (empty)
182 (*nik) = oik - min (der1, der2);
183 else
184 (*nik) = 2 * oik - der1 - der2 - 1;
185 if ((*nik) < 2)
186 (*nik) = 2;
187 *nin = 0;
188
189
190 /* Make resolution to be used for testing of knot value equalness. */
191
192 eps = fabs (oknots[oin] - oknots[oik - 1]) * 1.0e-11;
193
194
195 /* Production of knots. Initiate for calculation of knots.
196 Find first knot not equal to start of curve. */
197
198 maximum = oknots[oin];
199 prev = oknots[oik - 1];
200 for (kl = oik; prev >= oknots[kl]; kl++) ;
201
202 curr = oknots[kl];
203 for (mult = oik; curr < maximum; mult++)
204 {
205 if (curr < prev)
206 goto err112;
207
208 if (prev > curr || curr > prev + eps)
209 {
210
211 /* New knot value found. Fill in old value. */
212
213 /* numb = (*nik) - oik + mult; */
214 numb = (*nik) - oik + mult + der;
215 if (numb > (*nik))
216 numb = (*nik);
217
218
219 /* If numb >= nik, test if all the numb knots are equal
220 or if they only are equal within the resolution eps.
221 If not totally equal knumb=nik-1. */
222
223 if (numb == (*nik))
224 {
225 /* start = max (kl - oik, 1);
226 stop = kl - 2;
227 for (count1 = start; count1 <= stop; count1++)
228 if (oknots[count1 - 1] != oknots[count1])
229 numb = (*nik) - 1; */
230
231 start = kl - oik + der;
232 stop = kl - 2;
233 for (count1 = start; count1 <= stop; count1++)
234 if (oknots[count1] != oknots[count1 + 1])
235 numb = (*nik) - 1;
236 }
237
238 if (prev == oknots[oik - 1])
239 numb = (*nik);
240 for (count1 = 1; count1 <= numb; count1++)
241 (*nknots)[(*nin)++] = prev;
242
243
244 /* Initialize multiplicity. */
245
246 mult = 0;
247 prev = curr;
248 }
249 kl++;
250 curr = oknots[kl];
251 }
252
253 /* Knot for the next last knot value not produced. */
254
255 /* numb = min ((*nik) - oik + mult, (*nik)); */
256 numb = min ((*nik) - oik + mult + der, (*nik));
257
258
259 /* If numb >= nik, test if all the numb knots are equal or if they
260 * only are equal within the resolution eps. */
261
262 /* I not totally equal numb=nik-1. */
263
264 if (numb >= (*nik))
265 {
266 /* start = max (kl - oik, 1);
267 stop = kl - 2;
268 for (count1 = start; count1 <= stop; count1++)
269 if (oknots[count1 - 1] != oknots[count1])
270 numb = (*nik) - 1; */
271
272 start = kl - oik + der;
273 stop = kl - 2;
274 for (count1 = start; count1 <= stop; count1++)
275 if (oknots[count1] != oknots[count1 + 1])
276 numb = (*nik) - 1;
277 }
278
279 for (count1 = 1; count1 <= numb; count1++)
280 (*nknots)[(*nin)++] = prev;
281
282
283 /* Knot at et[oin+1] not produced. */
284
285 for (count1 = 1; count1 <= (*nik); count1++)
286 (*nknots)[(*nin)++] = maximum;
287
288
289 /* Knots produced. Correct nin and length of nknots. */
290
291 (*nin) -= (*nik);
292 *nknots = increasearray (*nknots, (*nik) + (*nin), DOUBLE);
293 if (*nknots == SISL_NULL)
294 goto err101;
295
296 goto out;
297
298 /* Not enough memory. */
299
300 err101:
301 *jstat = -101;
302 s6err ("s1894", *jstat, kpos);
303 goto out;
304
305 /* Error in description of B-spline. */
306
307 err112:
308 *jstat = -112;
309 s6err ("s1894", *jstat, kpos);
310 goto out;
311
312 out:
313 return;
314 }
315