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: sh1779.c,v 1.2 2001-03-19 15:59:05 afr Exp $
45 *
46 */
47
48
49 #define SH1779
50
51 #include "sislP.h"
52
53
54 #if defined(SISLNEEDPROTOTYPES)
55 void
sh1779(SISLObject * po1,SISLObject * po2,double aepsge,SISLIntdat ** rintdat,SISLIntpt * pintpt,int * jnewpt,int * jstat)56 sh1779 (SISLObject * po1, SISLObject * po2, double aepsge,
57 SISLIntdat ** rintdat, SISLIntpt * pintpt,
58 int *jnewpt, int *jstat)
59 #else
60 void
61 sh1779 (po1, po2, aepsge, rintdat, pintpt, jnewpt, jstat)
62 SISLObject *po1;
63 SISLObject *po2;
64 double aepsge;
65 SISLIntdat **rintdat;
66 SISLIntpt *pintpt;
67 int *jnewpt;
68 int *jstat;
69 #endif
70
71 /* UPDATE: (ujk) : Must tune the test on when to use local
72 info (=subtask no ?),
73 As in sh1780 it is necessary to determine when a helppoint
74 is close enough to a mainpoint to be neglected.
75 sh1784 does not tell if it march outside the surface,
76 this might be a problem. */
77
78 /*
79 *********************************************************************
80 *
81 *********************************************************************
82 *
83 * PURPOSE : Set pre-topology data in 3-dimensional curve-surface
84 * intersection. Also find help points if necessary.
85 *
86 *
87 * INPUT : po1 - Pointer to the first object in the intersection.
88 * po2 - Pointer to the second object in the intersection.
89 * aepsge - Geometry resolution.
90 * rintdat - Intersection data structure.
91 * pintpt - Current intersection point.
92 *
93 *
94 * OUTPUT : pintpt - Intersection point after updating pre-topology.
95 * jnewpt - Number of new int.pt. created.
96 * jstat - status messages
97 * > 0 : Warning.
98 * = 0 : Ok.
99 * < 0 : Error.
100 *
101 *
102 * METHOD :
103 *
104 * CALLS : shevalc - ? (s1221 used )Evaluate curve using filtered coefficients.
105 * s1421 - Evaluate surface using filtered coefficients.
106 * sh1784 - March along curve as long as it coincide
107 * with a surface.
108 * s6ang - Angle between two vectors.
109 * s6length - Length of vector.
110 * s6dist - Distance between two points.
111 * sh6idcon - Connect intersection points.
112 * sh6genhbrs - Get neighbours.
113 * sh6getgeo - Get geometric info from intersection point.
114 * sh6settop - Set topology of intersection point.
115 * hp_newIntpt - Create new intersection point.
116 *
117 * REFERENCES :
118 *
119 * WRITTEN BY : Vibeke Skytt, SI, Oslo, Norway. 04.91.
120 * Ulf J. Krystad, SI, 06.91.
121 *********************************************************************
122 */
123 {
124 int kstat = 0; /* Status variable. */
125 int ki; /* Counters. */
126 int kleft1 = 0, kleft2 = 0; /* Parameters to the evaluator. */
127 int kdim; /* Dimension of geometry space. */
128 int kpos = 0; /* Current position in int.pt. array. */
129 int kpar1, kpar2; /* Index of parameter value of object. */
130 int kn; /* Number of vertices of curve. */
131 int kk; /* Order of curve. */
132 int kmarch = 0; /* Indicates if marching is necessary. */
133 int lleft[2]; /* Array storing pre-topology information. */
134 int lright[2]; /* Array storing pre-topology information. */
135 int *ll1, *ll2, *lr1, *lr2; /* Pointers into pre-topology arrays. */
136 double tref; /* Referance value in equality test. */
137 double *st; /* Knot vector of curve. */
138 double sder[9]; /* Result of curve evaluation. */
139 double stang[3]; /* Tangent vector of curve. */
140 double snorm[3]; /* Normal vector of surface. */
141 double slast[3]; /* Last parameter value of coincidence. */
142 double snext[3]; /* First parameter value outside interval
143 of coincidence. */
144 double *ret_val; /* Pointer to geo data from sh6getgeom */
145 double *ret_norm; /* Pointer to geo data from sh6getgeom */
146 double *sptpar = pintpt->epar;/* Pointer to parameter values of int.pt. */
147 SISLCurve *qc; /* Pointer to the curve. */
148 SISLSurf *qs; /* Pointer to the surface. */
149 SISLIntpt *uintpt[2]; /* Array containing new intersection points. */
150 SISLIntpt *qpt1, *qpt2; /* Intersection points in list. */
151 double *nullp = SISL_NULL;
152 double sf_low_lim[2];
153 double sf_high_lim[2];
154
155 /* Don't make pretop for help points ! */
156 if (sh6ishelp (pintpt))
157 {
158 *jstat = 0;
159 goto out;
160 }
161
162 /* Set pointers into the arrays storing pre-topology information. */
163
164 if (po1->iobj == SISLCURVE)
165 {
166 qc = po1->c1;
167 qs = po2->s1;
168
169 kpar1 = 0;
170 kpar2 = 1;
171 ll1 = lleft;
172 lr1 = lright;
173 ll2 = lleft + 1;
174 lr2 = lright + 1;
175 }
176 else
177 {
178 qc = po2->c1;
179 qs = po1->s1;
180
181 kpar1 = 2;
182 kpar2 = 0;
183 ll1 = lleft + 1;
184 lr1 = lright + 1;
185 ll2 = lleft;
186 lr2 = lright;
187 }
188
189 /* Get pre-topology of intersection point. */
190 sh6gettop (pintpt, -1, lleft, lright, lleft + 1, lright + 1, &kstat);
191
192 /* Describe curve partly by local parameters. */
193
194 kdim = qc->idim;
195 kk = qc->ik;
196 kn = qc->in;
197 st = qc->et;
198 tref = st[kn] - st[kk - 1];
199
200 sf_low_lim[0] = qs->et1[qs->ik1 - 1] + REL_COMP_RES;
201 sf_low_lim[1] = qs->et2[qs->ik2 - 1] + REL_COMP_RES;
202 sf_high_lim[0] = qs->et1[qs->in1] - REL_COMP_RES;
203 sf_high_lim[1] = qs->et2[qs->in2] - REL_COMP_RES;
204
205 /* Fetch geometry information, curve. */
206 sh6getgeom ((po1->iobj == SISLCURVE) ? po1 : po2,
207 (po1->iobj == SISLCURVE) ? 1 : 2,
208 pintpt, &ret_val, &ret_norm, aepsge, &kstat);
209 if (kstat < 0)
210 goto error;
211
212 /* Local copy of curve tangent */
213 memcopy (stang, ret_val + kdim, kdim, DOUBLE);
214
215 /* Fetch geometry information, surface. */
216 sh6getgeom ((po1->iobj == SISLSURFACE) ? po1 : po2,
217 (po1->iobj == SISLSURFACE) ? 1 : 2,
218 pintpt, &ret_val, &ret_norm, aepsge, &kstat);
219 if (kstat < 0)
220 goto error;
221
222 /* Local copy of surface normal */
223 memcopy (snorm, ret_norm, kdim, DOUBLE);
224
225
226 /* (ALA) Test if local information may be used to compute pre-topology. */
227 s6length (snorm, kdim, &kstat);
228 s6length (snorm, kdim, &ki);
229
230 if (!kstat || !ki || fabs (PIHALF - s6ang (snorm, stang, kdim)) < 0.05)
231 {
232 /* Check if the intersection point lies at the start point of
233 the curve. */
234
235 if (DEQUAL (sptpar[kpar1] + tref, st[kn] + tref))
236 ;
237 else
238 {
239 /* Check if the intersection point is member of a list
240 in this parameter direction of the curve. */
241
242 qpt1 = qpt2 = SISL_NULL;
243 kmarch = 1;
244
245 /* UPDATE (ujk) : only one list ? */
246 sh6getnhbrs (pintpt, &qpt1, &qpt2, &kstat);
247 if (kstat < 0)
248 goto error;
249
250 kmarch = 0;
251 if (qpt1 != SISL_NULL && qpt1->epar[kpar1] > sptpar[kpar1])
252 *lr1 = SI_ON;
253 else if (qpt2 != SISL_NULL && qpt2->epar[kpar1] > sptpar[kpar1])
254 *lr1 = SI_ON;
255 else
256 kmarch = 1;
257 }
258
259 if (kmarch)
260 {
261 /* Perform marching to compute pre-topology. March first in the
262 positive direction of the curve. */
263
264 sh1784 (qc, qs, aepsge, sptpar, (kpar1 == 0), 1, slast, snext, &kstat);
265 if (kstat < 0)
266 goto error;
267
268 if (kstat == 1)
269 {
270 /* The endpoint of the curve is reached. */
271 ;
272 }
273 else if (kstat == 2)
274 ;
275 else
276 {
277
278 if (slast[kpar2] > sf_high_lim[0] ||
279 slast[kpar2 + 1] > sf_high_lim[1] ||
280 slast[kpar2] < sf_low_lim[0] ||
281 slast[kpar2 + 1] < sf_low_lim[1])
282 ;
283 else
284 {
285
286
287 /* Create help point. First fetch geometry information. */
288
289 s1221 (qc, 0, slast[kpar1], &kleft1, sder, &kstat);
290 if (kstat < 0)
291 goto error;
292
293 s1221 (qc, 0, snext[kpar1], &kleft1, sder + kdim, &kstat);
294 if (kstat < 0)
295 goto error;
296 s6diff (sder + kdim, sder, kdim, stang);
297
298 s1421 (qs, 1, slast + kpar2, &kleft1, &kleft2, sder, snorm, &kstat);
299 if (kstat < 0)
300 goto error;
301
302 /* Discuss tangent- and normal vector, and set up pre-topology
303 in one direction of the curve. */
304
305 if (s6scpr (snorm, stang, kdim) > DZERO)
306 *lr1 = SI_OUT;
307 else
308 *lr1 = SI_IN;
309
310 /* UPDATE (ujk) : Tuning on distance */
311 if (s6dist (sptpar, slast, 3) > (double) 0.05 * tref)
312 {
313 /* Create help point. Set pre-topology data as undefined. */
314 /* UPDATE (ujk) : If calculated values is stored, kder must
315 be 1 for curve and 2 for surface (sh6getgeom). Should
316 shevalc be used in stead of s1221 ? */
317
318 uintpt[kpos] = SISL_NULL;
319 if ((uintpt[kpos] = hp_newIntpt (3, slast, DZERO, -SI_ORD,
320 lleft[0], lright[0], lleft[1],
321 lright[1], 0, 0, nullp, nullp)) == SISL_NULL)
322 goto err101;
323
324 kpos++;
325 }
326 }
327 }
328 }
329
330 /* Check if the intersection point lies at the end point of
331 the curve. */
332
333 kmarch = 0;
334 if (DEQUAL (sptpar[kpar1] + tref, st[kk - 1] + tref))
335 ;
336 else
337 {
338 /* Check if the intersection point is member of a list
339 in this parameter direction of the curve. */
340
341 qpt1 = qpt2 = SISL_NULL;
342 kmarch = 1;
343
344 /* UPDATE (ujk) : only one list ? */
345 /* UPDATE (ujk) : only one list ? */
346 sh6getnhbrs (pintpt, &qpt1, &qpt2, &kstat);
347 if (kstat < 0)
348 goto error;
349
350 kmarch = 0;
351 if (qpt1 != SISL_NULL && qpt1->epar[kpar1] < sptpar[kpar1])
352 *ll1 = SI_ON;
353 else if (qpt2 != SISL_NULL && qpt2->epar[kpar1] < sptpar[kpar1])
354 *ll1 = SI_ON;
355 else
356 kmarch = 1;
357
358 }
359
360 if (kmarch)
361 {
362 /* March in the negative direction of the curve. */
363
364 sh1784 (qc, qs, aepsge, sptpar, (kpar1 == 0), -1, slast, snext, &kstat);
365 if (kstat < 0)
366 goto error;
367
368 if (kstat == 1)
369 {
370 /* The endpoint of the curve is reached. */
371 ;
372 }
373 else if (kstat == 2)
374 ;
375 else
376 {
377 if (slast[kpar2] > sf_high_lim[0] ||
378 slast[kpar2 + 1] > sf_high_lim[1] ||
379 slast[kpar2] < sf_low_lim[0] ||
380 slast[kpar2 + 1] < sf_low_lim[1])
381 ;
382 else
383 {
384
385 /* Create help point. First fetch geometry information. */
386
387 s1221 (qc, 0, slast[kpar1], &kleft1, sder, &kstat);
388 if (kstat < 0)
389 goto error;
390
391 s1221 (qc, 0, snext[kpar1], &kleft1, sder + kdim, &kstat);
392 if (kstat < 0)
393 goto error;
394 s6diff (sder + kdim, sder, kdim, stang);
395
396 s1421 (qs, 1, slast + kpar2, &kleft1, &kleft2, sder, snorm, &kstat);
397 if (kstat < 0)
398 goto error;
399
400 /* Discuss tangent- and normal vector, and set up pre-topology
401 in one direction of the curve. */
402
403 if (s6scpr (snorm, stang, kdim) > DZERO)
404 *ll1 = SI_OUT;
405 else
406 *ll1 = SI_IN;
407
408 /* UPDATE (ujk) : Tuning on distance */
409 if (s6dist (sptpar, slast, 3) > (double) 0.05 * tref)
410 {
411 /* Create help point. Set pre-topology data as undefined. */
412 /* UPDATE (ujk) : If calculated values is stored, kder must
413 be 1 for curve and 2 for surface (sh6getgeom). Should
414 shevalc be used in stead of s1221 ? */
415
416 uintpt[kpos] = SISL_NULL;
417 if ((uintpt[kpos] = hp_newIntpt (3, slast, DZERO, -SI_ORD,
418 lleft[0], lright[0], lleft[1],
419 lright[1], 0, 0, nullp, nullp)) == SISL_NULL)
420 goto err101;
421
422 kpos++;
423 }
424 }
425 }
426 }
427 }
428 else
429 {
430 /* Pre-topology data of the curve may be computed from
431 local information. */
432
433 if (s6scpr (snorm, stang, kdim) > DZERO)
434 {
435 *ll1 = SI_IN;
436 *lr1 = SI_OUT;
437 }
438 else
439 {
440 *ll1 = SI_OUT;
441 *lr1 = SI_IN;
442 }
443
444 }
445
446 /* Update pre-topology of intersection point. */
447 /* UPDATE (ujk), index = -1 ?? */
448 sh6settop (pintpt, -1, lleft[0], lright[0], lleft[1], lright[1], &kstat);
449
450 /* Join intersection points, and set pretopology of help points. */
451
452 for (ki = 0; ki < kpos; ki++)
453 {
454 /* Help point ? */
455 if (sh6ishelp (uintpt[ki]))
456 sh6settop (uintpt[ki], -1, *(pintpt->left_obj_1), *(pintpt->right_obj_1),
457 *(pintpt->left_obj_2), *(pintpt->right_obj_2), &kstat);
458
459 sh6idcon (rintdat, &uintpt[ki], &pintpt, &kstat);
460 if (kstat < 0)
461 goto error;
462 }
463
464 /* Pre-topology information computed. */
465
466 *jnewpt = kpos;
467 *jstat = 0;
468 goto out;
469
470 /* Error in scratch allocation. */
471
472 err101:*jstat = -101;
473 goto out;
474
475
476 /* Error lower level routine. */
477
478 error:*jstat = kstat;
479 goto out;
480
481 out:
482 return;
483 }
484