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: sh1994.c,v 1.2 2001-03-19 15:59:07 afr Exp $
45 *
46 */
47
48
49 #define SH1994
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
sh1994(SISLSurf * s1,double aepsge,int * jstat)54 void sh1994(SISLSurf *s1,double aepsge,int *jstat)
55 #else
56 void sh1994(s1,aepsge,jstat)
57 SISLSurf *s1;
58 double aepsge;
59 int *jstat;
60 #endif
61 /*
62 *********************************************************************
63 *
64 *********************************************************************
65 *
66 * PURPOSE : Check if a point-surface intersection in one dimention
67 * is a simple case,
68 * i.e. the intersection will result in one single point.
69 *
70 *
71 *
72 * INPUT : s1 - Surface in the intersection problem.
73 * aepsge - Geometry resolution.
74 *
75 *
76 *
77 * OUTPUT : jstat - status messages
78 * = 1 : simpel case.
79 * = 0 : not simpel case.
80 * < 0 : error.
81 *
82 *
83 * METHOD :
84 *
85 *
86 * REFERENCES :
87 *
88 * CALLS :
89 *
90 * WRITTEN BY : TDO SI, 89-08.
91 * REWISED BY : Vibeke Skytt, SI, 91-02.
92 * UJK, SI,91-10 Bug in first direction when a row
93 * ends with 2 or more equal numbers.
94 * Also added a test to include
95 * tmax < tmin (both HUGE)
96 *********************************************************************
97 */
98 {
99 register int ki,kj,kh;
100 int kk1, kk2, kn1, kn2;
101 int kbez;
102
103 double tmaxt, tmaxs;
104 double tmint, tmins;
105 double tdiff;
106 double *scoef=SISL_NULL;
107
108 /* Init to simple case. */
109 *jstat = 1;
110
111 tmaxt = tmaxs = - HUGE;
112 tmint = tmins = HUGE;
113
114 /* Get surface attributes. */
115 kk1 = s1->ik1;
116 kk2 = s1->ik2;
117 kn1 = s1->in1;
118 kn2 = s1->in2;
119 kbez = (kk1 == kn1) && (kk2 == kn2);
120
121
122 /* If the surface is linear in some direction it is simpel case. */
123 if ((kk1 == 2 && kn1 == 2) || (kk2 == 2 && kn2 == 2)) goto out;
124
125
126 /* Run through vertices in first parameter direction to find
127 intervall of first derivative. */
128
129 /* UJK, 91-10 */
130 /* for (kj=0, scoef=s1->ecoef; kj<kn2; kj++,scoef++) */
131 for (kj=0, scoef=s1->ecoef; kj<kn2; kj++,scoef=s1->ecoef+kn1*kj)
132 for (tdiff=DZERO, ki=1; ki<kn1; ki+=kh, scoef+=kh)
133 {
134 for (kh=1; ki+kh<=kn1; kh++)
135 {
136 if (tdiff*(*(scoef+kh) - *(scoef+kh-1)) < DZERO)
137 {
138 scoef += (kh-1);
139 ki += (kh-1);
140 kh = 1;
141 }
142 tdiff = *(scoef + kh) - *scoef;
143 if (fabs(tdiff) >= aepsge) break;
144 }
145 if (ki+kh > kn1) break;
146
147 tmint = min(tmint,tdiff);
148 tmaxt = max(tmaxt,tdiff);
149 }
150
151 /* Run through vertices in second parameter direction to find
152 intervall of first derivative. */
153
154 for (ki=0; ki<kn1; ki++)
155 for (tdiff=DZERO, kj=1, scoef=s1->ecoef+ki; kj<kn2; kj+=kh, scoef+=kh*kn1)
156 {
157 for (kh=1; kj+kh<=kn2; kh++)
158 {
159 if (tdiff*(*(scoef+kh*kn1) - *(scoef+(kh-1)*kn1)) < DZERO)
160 {
161 scoef += (kh-1)*kn1;
162 kj += (kh-1);
163 kh = 1;
164 }
165 tdiff = *(scoef + kh*kn1) - *scoef;
166 if (fabs(tdiff) >= aepsge) break;
167 }
168 if (kj+kh > kn2) break;
169
170 tmins = min(tmins,tdiff);
171 tmaxs = max(tmaxs,tdiff);
172 }
173
174 /* UJK, 91-10, maybe parameters not set */
175 if (tmint > tmaxt || tmins > tmaxs)
176 {
177 *jstat = 1;
178 goto out;
179 }
180
181 /* The first derivatives decide directions of possible intersection curves. */
182 if (kbez && (tmint*tmaxt >=DZERO || tmins*tmaxs >=DZERO))
183 *jstat = 1;
184 else if (tmint*tmaxt > DZERO || tmins*tmaxs > DZERO)
185 *jstat = 1;
186 else if (tmint == tmaxt || tmins == tmaxs)
187 *jstat = 1;
188 else
189 /* Not a simple case. */
190 *jstat = 0;
191
192 goto out;
193 out: ;
194 }
195
196