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: s1119.c,v 1.2 2001-03-19 15:58:41 afr Exp $
45 *
46 */
47
48
49 #define S1119
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1119(double * ecoef,double * et1,double * et2,int ik1,int in1,int ik2,int in2,int * jsimple,int * jind1,int * jind2,int * jstat)55 s1119(double *ecoef,double *et1,double *et2,int ik1,int in1,int ik2,
56 int in2,int *jsimple,int *jind1,int *jind2,int *jstat)
57 #else
58 void s1119(ecoef,et1,et2,ik1,in1,ik2,in2,jsimple,jind1,jind2,jstat)
59 double *ecoef;
60 double *et1;
61 double *et2;
62 int ik1;
63 int in1;
64 int ik2;
65 int in2;
66 int *jsimple;
67 int *jind1;
68 int *jind2;
69 int *jstat;
70 #endif
71 /*
72 *********************************************************************
73 *
74 *********************************************************************
75 *
76 * PURPOSE : Check if a one-dimensional surface can have only one
77 * single maximal-point.
78 *
79 *
80 *
81 * INPUT : ecoef - Vertices of surface.
82 * et1 - Knots, parameterdirection one.
83 * et2 - Knots, parameterdirection two.
84 * ik1 - Order of surface in first parameter direction.
85 * in1 - Number of vertices in first parameter direction.
86 * ik2 - Order of surface in second parameter direction.
87 * in2 - Number of vertices in second parameter direction.
88 *
89 *
90 *
91 * OUTPUT : jsimple - Indicates if single case
92 * = 0 : No interior max possible.
93 * = 1 : Only one maximal point or curve possible.
94 * = 2 : More than one maximal point possible.
95 * jind1 - Index to an interior knot with multiplicity ik1 in et1.
96 * = 0 : No knot with interior multiplicity.
97 * jind1 - Index to an interior knot with multiplicity ik2 in et2.
98 * = 0 : No knot with interior multiplicity.
99 * jstat - status messages
100 * > 0 : warning
101 * = 0 : ok
102 * < 0 : error
103 *
104 *
105 * METHOD : Count number of times the control-polygon changes
106 * direction in each parameter direction. If maximum
107 * times of changing is greater than one, this is not
108 * a simple case (possibility of several maxima).
109 *
110 *
111 * REFERENCES :
112 *
113 *-
114 * CALLS :
115 *
116 * WRITTEN BY : UJK, SI, 89-06.
117 *
118 *********************************************************************
119 */
120 {
121 int ki,kj; /* Counters. */
122 int ksimple; /* Indicates if simple case. */
123 int ksimple1; /* Indicates if simple case. */
124 int ksimple2; /* Indicates if simple case. */
125 int ksign; /* Number of direction changes in line/column. */
126 int kconvex1; /* Flag, if true, we have no interior min
127 in first direction.*/
128 int kconcav1; /* Flag, if true, we have no interior max
129 in first direction.*/
130 int kconvex2; /* Flag, if true, we have no interior min
131 in second direction.*/
132 int kconcav2; /* Flag, if true, we have no interior max
133 in second direction.*/
134 int kbez; /* Flagging for bezier case. */
135 double tfirst; /* First non-zero difference between two adjacent vertices. */
136 double tprev; /* Previous difference between two adjacent vertices. */
137 double tdiff; /* Current difference between two adjacent vertices. */
138 double *s1; /* Pointer used to traverse array of vertices. */
139
140 /* First we search for interior knotmultiplicity in
141 both parameter directions*/
142 *jind1 = 0;
143 ksimple1 = 1;
144 if (in1 > 1)
145 for (ki=ik1+1; ki<in1 && ksimple1; ki++)
146 {
147 if (et1[ki] == et1[ki+ik1-1])
148 {
149 *jind1 = ki;
150 ksimple1 = 0;
151 }
152 }
153
154 *jind2 = 0;
155 ksimple2 = 1;
156 if (in2 > 1)
157 for (ki=ik2+1; ki<in2 && ksimple2; ki++)
158 {
159 if (et2[ki] == et2[ki+ik2-1])
160 {
161 *jind2 = ki;
162 ksimple2 = 0;
163 }
164 }
165
166
167 ksimple = ksimple1 && ksimple2;
168 kbez = (((ik1 == in1) && (ik2 == in2)) ? 1 : 0);
169
170 /* Count number of direction changes in first parameter direction. */
171 /* Notify that we cannot accept equal coeffisient neighbours when
172 we are in a none-bezier case. */
173
174 kconcav1 = kconvex1 = 1;
175
176 if (in1 > 1)
177 for (s1=ecoef,kj=0; kj<in2 && ksimple; kj++)
178 {
179 ksign = 0;
180 tfirst = DZERO;
181
182 for (ki=0; ki<in1-1 && ksimple; ki++,s1++)
183 {
184 tdiff = *(s1+1) - *s1;
185 if (DEQUAL(tdiff,DZERO) )
186 {
187 if (kbez == 0) ksimple = 0;
188 }
189 else if (DEQUAL(tfirst,DZERO) )
190 {
191 /* First none-zero vector, save it. */
192 tfirst = tdiff;
193 tprev = tdiff;
194 }
195
196 else if (tprev*tdiff < DZERO)
197 {
198 tprev = tdiff;
199 ksign++;
200 if (ksign > 1) ksimple = 0;
201 }
202 }
203
204
205 if (kbez == 0)
206 {
207 /* We permit status simple case only in bezier case.
208 However, if the surface is strictly concav in one
209 parameter direction, we have found
210 the max on the edges. */
211 kconvex1 = 0;
212 kconcav1 = ((tfirst < DZERO) && kconcav1);
213 }
214 else
215 {
216
217 kconvex1 = (((ksign == 0) ||
218 (ksign == 1 && tfirst >= DZERO)) && kconvex1);
219 kconcav1 = (((ksign == 0) ||
220 (ksign == 1 && tfirst <= DZERO)) && kconcav1);
221 }
222
223 ksimple = ((kconvex1 || kconcav1) && ksimple);
224 s1++;
225
226 }
227
228 /* Count number of direction changes in second parameter direction. */
229 kconcav2 = kconvex2 = 1;
230 if (in2 > 1)
231 for (kj=0; kj<in1 && ksimple; kj++)
232 {
233 ksign = 0;
234 tfirst = DZERO;
235 s1 = ecoef + kj;
236
237 for (ki=0; ki<in2-1 && ksimple; ki++,s1+=in1)
238 {
239 tdiff = *(s1+in1) - *s1;
240 if (DEQUAL(tdiff,DZERO) )
241 {
242 if (kbez == 0) ksimple = 0;
243 }
244 else if (DEQUAL(tfirst,DZERO) )
245 {
246 /* First none-zero vector, save it. */
247 tfirst = tdiff;
248 tprev = tdiff;
249 }
250
251 else if (tprev*tdiff < DZERO)
252 {
253 tprev = tdiff;
254 ksign++;
255 if (ksign > 1) ksimple = 0;
256 }
257 }
258
259 if (kbez == 0)
260 {
261 /* We permit status simple case only in bezier case.
262 However, if the surface is strictly concav in one
263 parameter direction, we have found
264 the max on the edges. */
265 kconvex2 = 0;
266 kconcav2 = ((tfirst < DZERO) && kconcav2);
267 }
268 else
269 {
270
271 kconvex2 = (((ksign == 0) ||
272 (ksign == 1 && tfirst >= DZERO)) && kconvex2);
273 kconcav2 = (((ksign == 0) ||
274 (ksign == 1 && tfirst <= DZERO)) && kconcav2);
275 }
276 ksimple = ((kconvex2 || kconcav2) && ksimple);
277 }
278
279 /* Simple case test performed. */
280
281 if (ksimple)
282 {
283 if (kconvex1 && kconvex2)
284 *jsimple = 1;
285 else
286 *jsimple = 0;
287 }
288 else
289 *jsimple = 2;
290 *jstat = 0;
291
292 return;
293 }
294
295