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: s9adsimp.c,v 1.2 2001-03-19 15:59:02 afr Exp $
45 *
46 */
47
48
49 #define S9ADSIMP
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 double
s9adsimp(double epnt1[],double epar1[],double eimpli[],int ideg,double egd1[],double epgd1[],double etang[],double eptan[],double astep,int * jstat)55 s9adsimp(double epnt1[],double epar1[],double eimpli[],int ideg,double egd1[],
56 double epgd1[],double etang[],double eptan[],double astep,int *jstat)
57 #else
58 double s9adsimp(epnt1,epar1,eimpli,ideg,egd1,epgd1,etang,eptan,astep,jstat)
59 double epnt1[];
60 double epar1[];
61 double eimpli[];
62 int ideg;
63 double egd1[];
64 double epgd1[];
65 double etang[];
66 double eptan[];
67 double astep;
68 int *jstat;
69 #endif
70 /*
71 *********************************************************************
72 *
73 * PURPOSE : To decide if we should step directly on to a guide point
74 *
75 * INPUT : epnt1 - start point, in first surface of iteration step
76 * epar1 - Parameter value of epnt1
77 * eimpli - Description of the implicit surface
78 * ideg - The degree of the implicit surface
79 * ideg=1: Plane
80 * ideg=2; Quadric surface
81 * ideg=1001: Torus surface
82 * ideg=1003: Silhouette line parallel projection
83 * ideg=1004: Silhouette line perspective projection
84 * ideg=1005: Silhouette line circular projection
85 * egd1 - Guide point in first surface
86 * For ideg=1,2 and 1001 the sequence is position,
87 * first derivative in first parameter direction,
88 * first derivative in second parameter direction,
89 * (2,0) derivative, (1,1) derivative, (0,2) derivative
90 * and normal. (21 numbers)
91 * For ideg=1003,1004,1005 the second derivatives are followed
92 * by the third derivatives and the normal (33 numbers)
93 * Compatible with output of s1421
94 * epgd1 - Parameter value of egd1
95 * etang - Tangent in step direction at epoint
96 * eptan - Tangent in parameter plane at current point
97 * idim - Dimension of space the vectors lie in
98 * astep - Current step length
99 *
100 *
101 * OUTPUT : s9adsimp - Distance to guide point
102 * jstat - Status variable
103 * = 0 : Don't step onto the guide point
104 * = 1 : Step through the guide point
105 *
106 * METHOD : The step is changed if
107 * - The guide point lies in the tangent direction
108 * - The distance between the guide point and start point
109 * of iteration step is less than or equal to astep
110 *
111 * USE : This routine is only working in 3-D
112 *
113 *
114 * REFERENCES :
115 *
116 *-
117 * CALLS : s6diff,s6scpr,s6length,s6crss
118 *
119 * WRITTEN BY : Tor Dokken, SI, 1988-june-11
120 * Revised by : Tor Dokken, SI, 1989-03
121 * Corrected stepping through singular points.
122 * Revised by : Mike Floater, SI, 1991-01
123 * Add perspective and circular silhouettes (ideg=1004,ideg=1005)
124 *
125 *********************************************************************
126 */
127 {
128 int kpos=1; /* Position indicator for errors */
129 int kstat; /* Dummy status variable */
130 int kdim=3; /* This routine is only working in 3-D */
131 int k2dim=2; /* Dimension of parameter plane */
132 int ksize; /* Number of doubles for storage of derivateves
133 and normal vector */
134 int ksizem3; /* ksize - 3 */
135 double tdum; /* Variable for storage of reals */
136 double tdist=(double)0.0;/* Distance between guide point and point */
137 double sdiff[3]; /* Vector for difference between two points*/
138 double scr1[3],scr2[3]; /* Normal vectors */
139 double snorm[3]; /* Normal vector */
140
141
142 /* If ideg=1,2 or 1001 then only derivatives up to second order
143 are calculated, then 18 doubles for derivatives and 3 for the
144 normal vector are to be used for calculation of points in the
145 spline surface. For ideg=1003,1004,1005 we have a silhouette curve and
146 derivatives up to the third are to be calculated,
147 thus 30 +3 a total of 33 doubles are to be calculated */
148
149 if (ideg==1003 || ideg==1004 || ideg==1005)
150 {
151 ksize = 33;
152 }
153 else
154 {
155 ksize = 21;
156 }
157 ksizem3 = ksize -3;
158
159 /* First see that we are not turning direction in the parameter plane */
160
161 s6diff(epgd1,epar1,k2dim,sdiff);
162 if (s6scpr(sdiff,eptan,k2dim) < DZERO) goto dontstepthrough;
163
164 s6diff(egd1,epnt1,kdim,sdiff);
165 tdum = s6scpr(sdiff,etang,kdim);
166 tdist = s6length(sdiff,kdim,&kstat);
167
168 /* Step onto point if it is within 2.0*aepsge */
169
170 if (tdum > DZERO)
171 {
172 if (DZERO < tdist && tdist <= (double)2.0*astep)
173 {
174 /* Guide point lies within step length and in step direction, test
175 if cross products of normal vectors at current point and guide point
176 point in the same direction, this is not possible to test on for
177 silhouette curves */
178
179 if (ideg < 1003)
180 {
181
182 /* Make normal to implicit surface at epnt1 */
183 s1308(epnt1,3,eimpli,ideg,snorm,&kstat);
184 if (kstat<0) goto error;
185
186 /* Make cross product of normals in start point */
187 s6crss(epnt1+ksizem3,snorm,scr1);
188
189
190 /* Make normal to implicit surface at epnt1 */
191 s1308(epnt1,3,eimpli,ideg,snorm,&kstat);
192 if (kstat<0) goto error;
193
194 /* Make cross product of normals in guide point */
195 s6crss(egd1+ksizem3,snorm,scr2);
196
197
198 /* Make scalar product of these two vectors */
199 tdum = s6scpr(scr1,scr2,kdim);
200 if (kstat<0) goto error;
201
202 /* If positive scalar product the curve at the two points point in
203 the same direction, step through point */
204 if (tdum > DZERO) goto stepthrough;
205 else if (tdum == DZERO)
206 {
207
208 double tl1,tl2;
209
210 /* Vectors orthogonal or at least one has length zero */
211
212 tl1 = s6length(scr1,kdim,&kstat);
213 tl2 = s6length(scr2,kdim,&kstat);
214
215 if (tl1 != DZERO && tl2 != DZERO)
216 goto dontstepthrough;
217 else if (tl1 == DZERO && tl2 == DZERO)
218 goto stepthrough;
219 else if (tl2 == DZERO)
220 goto stepthrough;
221 else if (tl2 != DZERO)
222 {
223 /* Test if scr2 points in the direction from start */
224
225 tl1 = s6scpr(sdiff,scr2,kdim);
226
227 if (tl1 < DZERO)
228 goto dontstepthrough;
229 else
230 goto stepthrough;
231 }
232 }
233 }
234 else
235 goto stepthrough;
236 }
237 }
238
239 dontstepthrough:
240
241 *jstat = 0;
242 goto out;
243
244 stepthrough:
245 *jstat = 1;
246 goto out;
247
248 /* Error in lower leve function */
249 error:
250 *jstat = kstat;
251 s6err("s9adsimp",*jstat,kpos);
252 goto out;
253
254 out:
255 return(tdist);
256 }
257