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: s1753.c,v 1.2 2005-02-28 09:04:48 afr Exp $
45 *
46 */
47
48
49 #define S1753
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1753(double et[],double ecf[],int in,int ik,int idim,double etr[],double ecfr[],int inr,double ecc[],double ecw[],int * jstat)55 s1753 (double et[], double ecf[], int in, int ik, int idim, double etr[],
56 double ecfr[], int inr, double ecc[], double ecw[], int *jstat)
57
58 #else
59 void
60 s1753 (et, ecf, in, ik, idim, etr, ecfr, inr, ecc, ecw, jstat)
61 double et[];
62 double ecf[];
63 int in;
64 int ik;
65 int idim;
66 double etr[];
67 double ecfr[];
68 int inr;
69 double ecc[];
70 double ecw[];
71 int *jstat;
72
73 #endif
74 /*
75 *********************************************************************
76 *
77 *********************************************************************
78 *
79 * PURPOSE : To raise the description of a B-spline curve one order.
80 *
81 *
82 * INPUT : et - Description of knot vector of original description
83 * ecf - Coefficients of original description
84 * in - Number of vertices of original description
85 * ik - Order of original description
86 * idim - The dimension of the space in which the curve lies
87 * etr - Knot vector of the raised basis
88 * inr - Number of vertices in the raised curve
89 * ecc - Array for internal use only
90 * ecw - ---- " ----
91 *
92 * OUTPUT : ecfr - Knots of the raised curve
93 * jstat - Status variable:
94 * < 0 : error
95 * = 0 : OK.
96 *
97 * METHOD : The order raising algorithm of Cohen, Lyche and Schumaker
98 * is used.
99 *
100 *
101 * REFERENCES : Fortran version:
102 * T.Dokken, SI, 1984-06
103 *
104 *
105 * CALLS : s6err.
106 *
107 *
108 * WRITTEN BY : Christophe R. Birkeland, SI, 1991-07
109 * REWRITTEN BY :
110 * REVISED BY :
111 *
112 *********************************************************************
113 */
114 {
115 int ki, kj, kk, kl, kr, kstop;/* Loop control variables */
116 int kjmid, ikmid; /* kjmid=(kj-1)*idim ikmid=(ik-1)*idim */
117 int kpos = 0; /* Error position indicator */
118 double ty1, ty2, tyi, tyik; /* Parameters used in Main Loop */
119 double dummy;
120 double tden;
121
122 *jstat = 0;
123
124
125 /* Check input values. */
126
127 if ((ik < 1) || (in <ik) ||(inr < (ik + 1)))
128 goto err112;
129
130
131 /* Initiate local variables. */
132
133 kr = 1;
134 for (kj = 1; kj <= inr; kj++)
135 {
136
137 /* Find kr, such that et[kr-1]<=etr[kj-1]<et[kr] */
138
139 for (kr--; et[kr] <= etr[kj - 1]; kr++) ;
140
141
142 /* Set ecc and ecw to zero. */
143
144 for (ki = 0; ki < ik * idim; ki++)
145 {
146 ecc[ki] = (double) 0.0;
147 ecw[ki] = (double) 0.0;
148 }
149
150 /* Initialize the remaining ecc and ecw entries. */
151
152 kstop = MIN (ik, in +ik - kr);
153 for (ki = MAX (0, ik - kr); ki < kstop; ki++)
154 for (kl = 0; kl < idim; kl++)
155 {
156 dummy = ecf[(ki + kr - ik) * idim + kl];
157 ecc[ki * idim + kl] = dummy;
158 ecw[ki * idim + kl] = dummy;
159 }
160
161 /* MAIN LOOP. */
162
163 for (kk = ik - 1; kk > 0; kk--)
164 {
165 ty1 = etr[kj + kk - 1];
166 ty2 = etr[kj + kk];
167 kstop = MAX (ik - kk, ik - kr);
168
169 for (ki = MIN (ik - 1, in +2 * ik - kk - kr - 1); ki >= kstop; ki--)
170 {
171 tyi = et[kr + ki - ik];
172 tyik = et[kr + ki + kk - ik];
173 tden = tyik - tyi;
174
175 for (kl = 0; kl < idim; kl++)
176 {
177 ecc[ki * idim + kl] = ((ty2 - tyi) * ecc[ki * idim + kl] +
178 (tyik - ty2) * ecc[(ki - 1) * idim + kl]) / tden;
179 ecw[ki * idim + kl] = ((ty1 - tyi) * ecw[ki * idim + kl] +
180 (tyik - ty1) * ecw[(ki - 1) * idim + kl]) / tden +
181 ecc[ki * idim + kl];
182 }
183 }
184 }
185 kjmid = (kj - 1) * idim;
186 ikmid = (ik - 1) * idim;
187
188 for (kl = 0; kl < idim; kl++)
189 ecfr[kjmid + kl] = ecw[ikmid + kl] / ik;
190 }
191
192 goto out;
193
194
195 /* Error in description of bases */
196
197 err112:
198 *jstat = -112;
199 s6err ("s1753", *jstat, kpos);
200 goto out;
201
202 out:
203 return;
204 }
205