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: s1935.c,v 1.4 2001-03-19 15:58:56 afr Exp $
45 *
46 */
47
48
49 #define S1935
50
51 #include "sislP.h"
52
53
54 #if defined(SISLNEEDPROTOTYPES)
s1935(double et1[],int in1,double et2[],int in2,double * knt[],int * in,int ik,int * jstat)55 void s1935 (double et1[], int in1, double et2[], int in2,
56 double *knt[], int *in, int ik, int *jstat)
57 #else
58 void
59 s1935 (et1, in1, et2, in2, knt, in, ik, jstat)
60 double *et1;
61 int in1;
62 double *et2;
63 int in2;
64 double *knt[];
65 int *in;
66 int ik;
67 int *jstat;
68
69 #endif
70 /*
71 *********************************************************************
72 *
73 *********************************************************************
74 *
75 * PURPOSE: To produce a knot vector that is the union of two knot
76 * vectors. The original knot vectors and the produced knot
77 * vector are of the same order.
78 *
79 *
80 * INPUT: et1 - The first knot vector.
81 * in1 - The number of degrees of freedom in the
82 * B-basis given by the first knot vector.
83 * et2 - The second knot vector.
84 * in2 - The number of degrees of freedom in the
85 * B-basis given by the second knot vector.
86 * ik - The order of the knot vector to be produced.
87 *
88 * OUTPUT: jstat - Output status:
89 * < 0: Error.
90 * = 0: Ok.
91 * > 0: Warning.
92 *
93 * METHOD:
94 *
95 * REFERENCES: Fortran version:
96 * T.Dokken, SI, 1981-11
97 *
98 * CALLS: s6err.
99 *
100 * WRITTEN BY: Christophe R. Birkeland, SI, 1991-07
101 * CORRECTED BY: Vibeke Skytt, SI, 92-10. Removed exact test on equality
102 * of knots.
103 * CORRECTED BY: Christophe R. Birkeland, SI, 93-05.
104 * Reinstalled exact test on equality of knots. Necessary
105 * for correct working of routine s1931-s1937.
106 * CORRECTED BY: Paal Fugelli, SINTEF, 1994-07.
107 * Changed the setting of 'curr', for the while loop, to avoid
108 * over running array bounds (address error) and added DEQUAL in
109 * tests for knot equality.
110 *
111 *********************************************************************
112 */
113 {
114 int kpos = 0; /* Error position indicator. */
115 int pek1, pek2; /* Index for array et1 and et2. */
116 int stop1; /* Length of et1, equals in1+ik. */
117 int stop2; /* Length of et2, equals in2+ik. */
118 double curr; /* Parameter used in calculation of
119 new knots. */
120
121 *jstat = 0;
122
123 /* Test if legal input */
124
125 if (ik < 1) goto err110;
126 if ((in1 < ik) || (in2 < ik)) goto err111;
127
128
129 /* Allocate array for new knot vector */
130
131 if((*knt = newarray (2 * ik + in1 + in2, DOUBLE))==SISL_NULL) goto err101;
132
133 /* Test if input knot vectors degenerate */
134
135 if (et1[ik - 1] >= et1[in1]) goto err112;
136
137 if (et2[ik - 1] >= et2[in2]) goto err112;
138
139 /* PRODUCTION OF KNOTS */
140
141 *in = 0;
142 /* PFU 05/07-94 curr = MIN (et1[0], et2[0]); */
143 pek1 = 0;
144 pek2 = 0;
145 stop1 = in1 + ik;
146 stop2 = in2 + ik;
147
148 while ((pek1 < stop1) && (pek2 < stop2))
149 {
150 /* Test if error in knot vector */
151
152 curr = MIN (et1[pek1], et2[pek2]);
153
154 if ((et1[pek1] < curr) || (et2[pek2] < curr)) goto err112;
155
156 if ( DEQUAL(et1[pek1],curr) ) pek1++;
157 if ( DEQUAL(et2[pek2],curr) ) pek2++;
158 (*knt)[*in] = curr;
159 (*in) ++;
160 }
161
162 /* Some knots may remain in one of the arrays */
163
164 if ((pek1 < stop1) || (pek2 < stop2))
165 {
166 if (pek1 >= stop1)
167 {
168 for (; pek2 < stop2; pek2++, (*in) ++)
169 (*knt)[*in] = et2[pek2];
170 }
171 else
172 for (; pek1 < stop1; pek1++, (*in) ++)
173 (*knt)[*in] = et1[pek1];
174 }
175
176 /* Knots produced */
177
178 *in -=ik;
179 *knt = increasearray (*knt, *in +ik, DOUBLE);
180 if (*knt == SISL_NULL) goto err101;
181 goto out;
182
183
184 /* Memory error */
185
186 err101:
187 *jstat = -101;
188 s6err ("s1935", *jstat, kpos);
189 goto out;
190
191 /* Error in description of B-spline curve: Order less than 1.*/
192
193 err110:
194 *jstat = -110;
195 s6err ("s1935", *jstat, kpos);
196 goto out;
197
198 /* No. of vertices less than order. */
199
200 err111:
201 *jstat = -111;
202 s6err ("s1935", *jstat, kpos);
203 goto out;
204
205 /* Error in knot vector */
206
207 err112:
208 *jstat = -112;
209 s6err ("s1935", *jstat, kpos);
210 goto out;
211
212 out:
213 return;
214 }
215