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: sh6evalint.c,v 1.3 2001-03-19 15:59:07 afr Exp $
45 *
46 */
47
48
49 #define SH6EVALINT
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
sh6evalint(SISLObject * ob1,SISLObject * ob2,double eimpli[],int ideg,SISLIntpt * pt,double aepsge,double * curve_3d[],double * curve_2d_1[],double * curve_2d_2[],int * jstat)55 sh6evalint (SISLObject * ob1, SISLObject * ob2, double eimpli[], int ideg,
56 SISLIntpt * pt, double aepsge,
57 double *curve_3d[], double *curve_2d_1[], double *curve_2d_2[],
58 int *jstat)
59 #else
60 void
61 sh6evalint (ob1, ob2, eimpli, ideg,
62 pt, aepsge,
63 curve_3d, curve_2d_1, curve_2d_2,
64 jstat)
65 SISLObject *ob1;
66 SISLObject *ob2;
67 double eimpli[];
68 int ideg;
69 SISLIntpt *pt;
70 double aepsge;
71 double *curve_3d[];
72 double *curve_2d_1[];
73 double *curve_2d_2[];
74 int *jstat;
75 #endif
76 /*
77 *********************************************************************
78 *
79 *********************************************************************
80 *
81 * PURPOSE : Given an Intpt, get the tangent and curvature of intersection
82 *
83 *
84 * INPUT : ob1 - Pointer to geometric object(first surf).
85 * ob2 - Pointer to geometric object(second surf).
86 * eimpli - Array containing descr. of implicit surf
87 * ideg - Type of impl surf.
88 * pt - Pointer to the Intpt.
89 * aepsge - Absolute tolerance
90 *
91 *
92 * OUTPUT : curve_3d - Interscetion data, object space (as in s1304, s1306)
93 * curve_2d_1 - Interscetion data, param space (as in s1304, s1306)
94 * curve_2d_2 - Interscetion data, param space (as in s1304)
95 *
96 *
97 * METHOD :
98 *
99 *
100 * REFERENCES :
101 *
102 * WRITTEN BY : Ulf J. Krystad, SI, Oslo, Norway. July 91.
103 * REVICED BY : UJK, Changed to deal with SISL_Curves also.
104 * Revised by : Paal Fugelli, SINTEF, Oslo, Norway, Sept. 1994. Avoid array bound
105 * over-run in memcopy.
106 *********************************************************************
107 */
108 {
109 int dim; /* Geometric dimension. */
110 int kstat;
111 int kpos = 1; /* Position indicator ofr errors */
112 int left1 = 0, left2 = 0; /* Knot navigators in s1421 */
113 int kder = 2; /* Numb of derivatives */
114 int silhouett; /* Flag silhouett case */
115 int ki; /* Variable used in loop */
116 int ksize; /* Size of output from s1421 or getgeom */
117 double *geom1 = SISL_NULL; /* Output values from s1421 or getgeom */
118 double con_tang[3]; /* Constant tangent. */
119 double *norm1 = SISL_NULL; /* Output values from s1421 or getgeom */
120 double *geom2 = SISL_NULL; /* Output values from s1421 or getgeom */
121 double *norm2 = SISL_NULL; /* Output values from s1421 or getgeom */
122 double normimpl[3]; /* Normal of impl surf */
123 double right_dir[3]; /* Right direction of 3D intersect. curve */
124 double dot; /* Scalar product */
125 double dummy[6];
126 double ang;
127 double min_hp_ang = 0.00000000001;
128 *jstat = 0;
129 con_tang[0] = (double) 1.0;
130 con_tang[1] = DZERO;
131 con_tang[2] = DZERO;
132
133
134 if (ob1->iobj != SISLSURFACE && ob1->iobj != SISLCURVE)
135 goto errinp;
136 if (!pt)
137 goto errinp;
138 if (ideg < 0)
139 goto errinp;
140
141 if (ob1->iobj == SISLSURFACE )
142 dim = ob1->s1->idim;
143 else
144 dim = ob1->c1->idim;
145
146 if (dim > 3 || dim < 1)
147 goto errinp;
148
149 *curve_3d = pt->geo_track_3d;
150 *curve_2d_1 = pt->geo_track_2d_1;
151 *curve_2d_2 = pt->geo_track_2d_2;
152
153 if (pt->evaluated)
154 goto out;
155
156 if (ideg == 0)
157 {
158 /* No implicit geometry involved */
159 kpos = 1;
160 if (ob2->iobj != SISLSURFACE && ob2->iobj != SISLCURVE)
161 goto errinp;
162
163 if (ob2->iobj == SISLCURVE)
164 {
165 /* At least the second object is a spline curve,
166 use this one */
167 kpos = 2;
168 if (ob2->c1->idim > 3) goto errinp;
169
170 /* Get geometry of first surface */
171 sh6getgeom (ob1, 1, pt, &geom1, &norm1, aepsge, &kstat);
172 if (kstat < 0)
173 goto error;
174
175 /* Get geometry of objects */
176 sh6getgeom (ob2, 2, pt, &geom2, &norm2, aepsge, &kstat);
177 if (kstat < 0)
178 goto error;
179
180 /* The number of elements to copy is given by pt->size_<obnr>
181 and we have obnr=2 (PFU 05/09-94) */
182 memcopy(*curve_3d,geom2,pt->size_2,double);
183
184 }
185 else
186 {
187
188
189 /* Two 3d surfaces */
190 kpos = 3;
191 if (ob2->iobj != SISLSURFACE)
192 goto errinp;
193 if ((dim = ob2->s1->idim) != 3)
194 goto errinp;
195
196 /* Get geometry of first surface */
197 sh6getgeom (ob1, 1, pt, &geom1, &norm1, aepsge, &kstat);
198 if (kstat < 0)
199 goto error;
200
201 /* Get geometry of second surface */
202 sh6getgeom (ob2, 2, pt, &geom2, &norm2, aepsge, &kstat);
203 if (kstat < 0)
204 goto error;
205
206 /* Get normal direction */
207 s6crss (norm1, norm2, right_dir);
208 if (kstat < 0)
209 goto error;
210
211 /* Compute angle. */
212 ang = s6ang(norm1, norm2,3);
213 if (ang < min_hp_ang)
214 {
215 /* The point is a singular meeting point.*/
216 if (pt->iinter == SI_ORD) pt->iinter = SI_SING;
217 }
218
219 /* Get tangent and curvature */
220 s1304 (geom1, geom2, pt->epar, pt->epar + 2,
221 *curve_3d, *curve_2d_1, *curve_2d_2, &kstat);
222 if (kstat < 0)
223 goto error;
224
225 if ((dot = s6scpr (right_dir, *curve_3d + 3, 3)) < DZERO)
226 {
227 /* Change direction for tangent */
228 for (ki = 0; ki < 3; ki++)
229 (*curve_3d)[ki + 3] *= -(double) 1;
230 for (ki = 0; ki < 2; ki++)
231 {
232 (*curve_2d_1)[ki + 2] *= -(double) 1;
233 (*curve_2d_2)[ki + 2] *= -(double) 1;
234 }
235 }
236 }
237
238 pt->evaluated = TRUE;
239 }
240 else
241 {
242 /* Implicit cases */
243 if (ideg == 2000)
244 {
245 /* Here we treat the cases
246 spline surf vs implicit analytic curve
247 spline curve vs implicit analytic curve
248 spline curve vs implicit analytic surf
249 in all these cases only 3D posisition is necessary */
250
251 /* Clean up from 1D or 2D result */
252 if (pt->geo_data_1)
253 freearray (pt->geo_data_1);
254 if (pt->geo_data_2)
255 freearray (pt->geo_data_2);
256 pt->geo_data_1 = SISL_NULL;
257 pt->size_1 = 0;
258 pt->geo_data_2 = SISL_NULL;
259 pt->size_2 = 0;
260
261 /* Get the right values are computed */
262 sh6getgeom (ob1, 1, pt, &geom1, &norm1, aepsge, &kstat);
263 if (kstat < 0)
264 goto error;
265
266 memcopy(*curve_3d,geom1,dim,double);
267 memcopy((*curve_3d)+dim,con_tang,dim,double);
268
269 }
270 else
271 {
272 if (ideg == 1003 || ideg == 1004 || ideg == 1005)
273 {
274 /* Silhouette cases, B-spline surface */
275 kpos = 3;
276 ksize = 33;
277 silhouett = TRUE;
278 kder = 3;
279
280 }
281 else
282 {
283 /* Analytic surf vs B-spline surface */
284
285 kpos = 4;
286 ksize = 21;
287 silhouett = FALSE;
288 kder = 2;
289 }
290
291 if (pt->size_1 != ksize)
292 {
293 /* Clean up from 1D result */
294 if (pt->geo_data_1)
295 freearray (pt->geo_data_1);
296 if (pt->geo_data_2)
297 freearray (pt->geo_data_2);
298 pt->geo_data_1 = SISL_NULL;
299 pt->size_1 = 0;
300 pt->geo_data_2 = SISL_NULL;
301 pt->size_2 = 0;
302
303
304 if ((pt->geo_data_1 = newarray (ksize, DOUBLE))
305 == SISL_NULL)
306 goto err101;
307 pt->size_1 = ksize;
308 geom1 = pt->geo_data_1;
309 norm1 = pt->geo_data_1 + ksize - 3;
310
311 s1422 (ob1->s1, kder, pt->iside_1, pt->iside_2,
312 pt->epar, &left1, &left2, geom1,
313 norm1, &kstat);
314 if (kstat < 0)
315 goto error;
316 }
317 else
318 {
319 /* The right values are computed */
320 sh6getgeom (ob1, 1, pt, &geom1, &norm1, aepsge, &kstat);
321 if (kstat < 0)
322 goto error;
323
324 }
325
326
327 /* Get normal of implicit surface */
328 s1331 (geom1, eimpli, ideg, kder = -1, dummy, normimpl, &kstat);
329 if (kstat < 0)
330 goto error;
331
332 /* Get the right direction of the intersection curve */
333 if (silhouett)
334 {
335 ang = 1.5; /* Not used */
336 memcopy (right_dir, normimpl, 3, DOUBLE);
337 for (ki=0;ki<3;ki++) right_dir[ki] *= -(double)1.0;
338 }
339 else
340 {
341 /* Compute angle. */
342 ang = s6ang(norm1, normimpl,3);
343 s6crss (norm1, normimpl, right_dir);
344 }
345
346 /* Get tangent and curvature to the real intersection. */
347 s1306 (geom1, pt->epar,
348 eimpli, ideg, *curve_3d, *curve_2d_1, &kstat);
349 if (kstat < 0)
350 goto error;
351 if (kstat == 2)
352 {
353 /* The point is a singular meeting point.*/
354 if (pt->iinter == SI_ORD) pt->iinter = SI_SING;
355 }
356 else if (kstat == 10)
357 {
358 /* The point is a singular non-meeting point.
359 Tangent found, but sign might be wrong. */
360 if (pt->iinter == SI_ORD || pt->iinter == SI_SING )
361 pt->iinter = SI_TOUCH;
362 }
363 else if (ang < min_hp_ang)
364 {
365 /* The point is a singular meeting point.*/
366 if (pt->iinter == SI_ORD) pt->iinter = SI_SING;
367 }
368 else
369 if ((dot = s6scpr (right_dir, *curve_3d + 3, 3)) < DZERO)
370 {
371 /* Change direction for tangent */
372 for (ki = 0; ki < 3; ki++)
373 (*curve_3d)[ki + 3] *= -(double) 1;
374 for (ki = 0; ki < 2; ki++)
375 (*curve_2d_1)[ki + 2] *= -(double) 1;
376
377 }
378 }
379
380
381 pt->evaluated = TRUE;
382
383 }
384
385 *jstat = 0;
386 goto out;
387
388 /* ---------- ERROR EXITS --------------------------- */
389 /* Error in alloc */
390 err101:
391 *jstat = -101;
392 s6err ("shevalint", *jstat, kpos);
393 goto out;
394
395 /* Error in lower level */
396 error:
397 *jstat = kstat;
398 s6err ("shevalint", *jstat, kpos);
399 goto out;
400
401 /* Error in input */
402 errinp:
403 *jstat = -200;
404 s6err ("shevalint", *jstat, kpos);
405 goto out;
406
407
408 out:;
409 }
410