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: sh1780.c,v 1.3 2005-02-28 09:04:50 afr Exp $
45 *
46 */
47
48
49 #define SH1780
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
sh1780(SISLObject * po1,SISLObject * po2,double aepsge,SISLIntdat ** rintdat,SISLIntpt * pintpt,int * jnewpt,int * jstat)55 sh1780 (SISLObject * po1, SISLObject * po2, double aepsge,
56 SISLIntdat ** rintdat, SISLIntpt * pintpt, int *jnewpt, int *jstat)
57 #else
58 void
59 sh1780 (po1, po2, aepsge, rintdat, pintpt, jnewpt, jstat)
60 SISLObject *po1;
61 SISLObject *po2;
62 double aepsge;
63 SISLIntdat **rintdat;
64 SISLIntpt *pintpt;
65 int *jnewpt;
66 int *jstat;
67 #endif
68
69 /* UPDATE: (ujk) : Must tune the test on when to use local
70 info (=subtask no ?),
71 As in sh1779 it is necessary to determine when a helppoint
72 is close enough to a mainpoint to be neglected. */
73 /*
74 *********************************************************************
75 *
76 *********************************************************************
77 *
78 * PURPOSE : Generate help points and set pre-topology data in
79 * curve-curve intersection.
80 *
81 *
82 * INPUT : po1 - Pointer to the first object in the intersection.
83 * po2 - Pointer to the second object in the intersection.
84 * aepsge - Geometry resolution.
85 * rintdat - Intersection data structure.
86 * pintpt - Current intersection point.
87 *
88 *
89 * OUTPUT : pintpt - Intersection point with updated pre-topology info.
90 * jnewpt - Number of new int.pt. created.
91 * jstat - status messages
92 * > 0 : Warning.
93 * = 0 : Ok.
94 * < 0 : Error.
95 *
96 *
97 * METHOD :
98 *
99 * CALLS : shevalc - ? (s1221 used )Evaluate curve using filtered coefficients.
100 * s1421 - Evaluate surface using filtered coefficients.
101 * sh1783 - March along two curves as long as they coincide.
102 * s6ang - Angle between two vectors.
103 * s6length - Length of vector.
104 * s6dist - Distance between two points.
105 * sh6idcon - Connect intersection points.
106 * sh6genhbrs - Get neighbours.
107 * sh6getgeo - Get geometric info from intersection point.
108 * sh6settop - Set topology of intersection point.
109 * hp_newIntpt - Create new intersection point.
110 *
111 * REFERENCES :
112 *
113 * WRITTEN BY : Vibeke Skytt, SI, Oslo, Norway. 04.91.
114 * Ulf J. Krystad, SI, 06.91
115 *********************************************************************
116 */
117 {
118 int kstat = 0; /* Status variable. */
119 int ki; /* Counters. */
120 int kleft1 = 0; /* Parameters to the evaluator. */
121 int kdim; /* Dimension of geometry space. */
122 int kpos = 0; /* Current position in output array. */
123 int kdir1, kdir2; /* Directions in which to march the curves.*/
124 int kk1, kk2; /* Orders of the two curves. */
125 int kn1, kn2; /* Number of vertices in the curves. */
126 int lleft[2]; /* Array storing pre-topology information. */
127 int lright[2]; /* Array storing pre-topology information. */
128 double tref; /* Reference value in equality test. */
129 double *st1, *st2; /* Pointers to knot vectors of curves. */
130 double sder[6]; /* Result of curve evaluation. */
131 double stang1[3]; /* Tangent vector of curve. */
132 double stang2[3]; /* Tangent vector of level value. */
133 double slast[3]; /* Last parameter value of coincidence. */
134 double snext[3]; /* First parameter value outside interval
135 of coincidence. */
136 double *ret_val; /* Pointer to geo data from sh6getgeom */
137 double *ret_norm; /* Pointer to geo data from sh6getgeom */
138 double *sptpar = pintpt->epar;/* Parameter array of int.pt. */
139 SISLIntpt *uintpt[2]; /* Pointer to new intersection points. */
140 double *nullp = SISL_NULL;
141
142 /* Don't make pretop for help points ! */
143 if (sh6ishelp (pintpt))
144 {
145 *jstat = 0;
146 goto out;
147 }
148
149 /* Test dimension of geometry space. */
150
151 kdim = po1->c1->idim;
152 if (kdim > 3)
153 goto err108;
154 if (kdim != po2->c1->idim)
155 goto err106;
156
157 /* Express the curve by local parameters. */
158
159 kn1 = po1->c1->in;
160 kk1 = po1->c1->ik;
161 st1 = po1->c1->et;
162 kn2 = po2->c1->in;
163 kk2 = po2->c1->ik;
164 st2 = po2->c1->et;
165 tref = MAX (st1[kn1] - st1[kk1 - 1], st2[kn2] - st2[kk2 - 1]);
166
167 /* Fetch already existing topology. */
168 sh6gettop (pintpt, -1, lleft, lright, lleft + 1, lright + 1, &kstat);
169
170 /* Fetch geometry information, first curve. */
171 sh6getgeom (po1, 1, pintpt, &ret_val, &ret_norm, aepsge, &kstat);
172 if (kstat < 0)
173 goto error;
174
175 /* Local copy of curve tangent */
176 memcopy (stang1, ret_val + kdim, kdim, DOUBLE);
177
178 /* Fetch geometry information,second curve. */
179 sh6getgeom (po2, 2, pintpt, &ret_val, &ret_norm, aepsge, &kstat);
180 if (kstat < 0)
181 goto error;
182
183 /* Local copy of curve tangent */
184 memcopy (stang2, ret_val + kdim, kdim, DOUBLE);
185
186 /* Compute the angle between the tangent vectors of the curves
187 in the current intersection point, and check if marching is
188 necessary to compute the pre-topology information. */
189
190 /* UPDATE (ujk) : tune */
191 if (s6ang (stang1, stang2, kdim) <= ANGULAR_TOLERANCE)
192 {
193 /* Perform marching in positive direction of the first curve. */
194
195 kdir1 = 1;
196 kdir2 = (s6scpr (stang1, stang2, kdim) >= DZERO) ? 1 : -1;
197
198 /* Check if the intersection point is situated at the endpoint
199 of a curve. */
200
201 if (DEQUAL (sptpar[0] + tref, st1[kn1] + tref) ||
202 (kdir2 == 1 && DEQUAL (sptpar[1] + tref, st2[kn2] + tref)) ||
203 (kdir2 == -1 && DEQUAL (sptpar[1] + tref, st2[kk2 - 1] + tref)))
204 {
205 }
206 else
207 {
208 /* Perform marching. */
209
210 sh1783 (po1->c1, po2->c1, aepsge, sptpar, kdir1, kdir2, slast,
211 snext, &kstat);
212 if (kstat < 0)
213 goto error;
214
215 if (kstat > 0)
216 {
217 /* An intersection interval is found. */
218 /* Set pre-topology */
219
220 lright[0] = SI_ON;
221 if (kdir2 == 1)
222 lright[1] = SI_ON;
223 else
224 lleft[1] = SI_ON;
225 }
226 else
227 {
228 /* Create help point. First fetch geometry information. */
229
230 s1221 (po1->c1, 0, slast[0], &kleft1, sder, &kstat);
231 if (kstat < 0)
232 goto error;
233
234 s1221 (po1->c1, 0, snext[0], &kleft1, sder + kdim, &kstat);
235 if (kstat < 0)
236 goto error;
237 s6diff (sder + kdim, sder, kdim, stang1);
238
239 s1221 (po2->c1, 0, slast[1], &kleft1, sder, &kstat);
240 if (kstat < 0)
241 goto error;
242
243 s1221 (po2->c1, 0, snext[1], &kleft1, sder + kdim, &kstat);
244 if (kstat < 0)
245 goto error;
246 s6diff (sder + kdim, sder, kdim, stang2);
247
248 /* Discuss directions of vectors and set up pre-topology
249 information in one direction of the curves. */
250
251 if ((stang1[0] * stang2[1] - stang1[1] * stang2[0]) * (double) kdir2
252 < DZERO)
253 lright[0] = SI_OUT;
254 else
255 lright[0] = SI_IN;
256
257 if (kdir2 == 1)
258 lright[1] = (lright[0] == SI_IN) ? SI_OUT : SI_IN;
259 else
260 lleft[1] = (lright[0] == SI_OUT) ? SI_OUT : SI_IN;
261
262 /* UPDATE (ujk) : tune */
263 if (s6dist (sptpar, slast, 2) > (double) 0.05 * tref)
264 {
265 /* Create help point. Set pre-topology data as SI_UNDEF. */
266
267 uintpt[kpos] = SISL_NULL;
268 if ((uintpt[kpos] = hp_newIntpt (2, slast, DZERO, -SI_ORD,
269 SI_UNDEF, SI_UNDEF, SI_UNDEF, SI_UNDEF,
270 0, 0, nullp, nullp)) == SISL_NULL)
271 goto err101;
272
273 kpos++;
274 }
275 }
276 }
277
278 /* Perform marching in negative direction of the first curve. */
279
280 kdir1 = -1;
281 kdir2 = -kdir2;
282
283 /* Check if the intersection point is situated at the endpoint
284 of a curve. */
285
286 if (DEQUAL (sptpar[0] + tref, st1[kk1 - 1] + tref) ||
287 (kdir2 == 1 && DEQUAL (sptpar[1] + tref, st2[kn2] + tref)) ||
288 (kdir2 == -1 && DEQUAL (sptpar[1] + tref, st2[kk2 - 1] + tref)))
289 {
290 }
291 else
292 {
293 /* Perform marching. */
294
295 sh1783 (po1->c1, po2->c1, aepsge, sptpar, kdir1, kdir2, slast,
296 snext, &kstat);
297 if (kstat < 0)
298 goto error;
299
300 if (kstat > 0)
301 {
302 /* An intersection interval is found. Set pre-topology. */
303
304 lleft[0] = SI_ON;
305 if (kdir2 == 1)
306 lright[1] = SI_ON;
307 else
308 lleft[1] = SI_ON;
309 }
310 else
311 {
312 /* Create help point. First fetch geometry information. */
313
314 s1221 (po1->c1, 0, slast[0], &kleft1, sder, &kstat);
315 if (kstat < 0)
316 goto error;
317
318 s1221 (po1->c1, 0, snext[0], &kleft1, sder + kdim, &kstat);
319 if (kstat < 0)
320 goto error;
321 s6diff (sder + kdim, sder, kdim, stang1);
322
323 s1221 (po2->c1, 0, slast[1], &kleft1, sder, &kstat);
324 if (kstat < 0)
325 goto error;
326
327 s1221 (po2->c1, 0, snext[1], &kleft1, sder + kdim, &kstat);
328 if (kstat < 0)
329 goto error;
330 s6diff (sder + kdim, sder, kdim, stang2);
331
332 /* Discuss directions of vectors and set up pre-topology
333 information in one direction of the curves. */
334
335 if ((stang1[0] * stang2[1] - stang1[1] * stang2[0]) * (double) kdir2
336 < DZERO)
337 lleft[0] = SI_OUT;
338 else
339 lleft[0] = SI_IN;
340
341 if (kdir2 == -1)
342 lleft[1] = (lleft[0] == SI_IN) ? SI_OUT : SI_IN;
343 else
344 lright[1] = (lleft[0] == SI_OUT) ? SI_OUT : SI_IN;
345
346 /* UPDATE (ujk) : tune */
347 if (s6dist (sptpar, slast, 2) > (double) 0.05 * tref)
348 {
349 /* Create help point. Set pre-topology data as SI_UNDEF. */
350
351 uintpt[kpos] = SISL_NULL;
352 if ((uintpt[kpos] = hp_newIntpt (2, slast, DZERO, -SI_ORD,
353 SI_UNDEF, SI_UNDEF, SI_UNDEF, SI_UNDEF,
354 0, 0, nullp, nullp)) == SISL_NULL)
355 goto err101;
356
357 kpos++;
358 }
359 }
360 }
361 }
362 else
363 {
364 /* The pretopology may be computed using local information. */
365
366 if (stang1[0] * stang2[1] - stang1[1] * stang2[0] < DZERO)
367 {
368 lleft[0] = SI_IN;
369 lright[0] = SI_OUT;
370 lleft[1] = SI_OUT;
371 lright[1] = SI_IN;
372 }
373 else
374 {
375 lleft[0] = SI_OUT;
376 lright[0] = SI_IN;
377 lleft[1] = SI_IN;
378 lright[1] = SI_OUT;
379 }
380
381
382 }
383
384 /* Update pre-topology of intersection point. */
385 /* UPDATE (ujk), index = -1 ?? */
386 sh6settop (pintpt, -1, lleft[0], lright[0], lleft[1], lright[1], &kstat);
387
388 /* Join intersection points, and set pretopology of help points. */
389
390 for (ki = 0; ki < kpos; ki++)
391 {
392 sh6idnpt (rintdat, &uintpt[ki], 1, &kstat);
393 if (kstat < 0)
394 goto error;
395
396 if (sh6ishelp (uintpt[ki]) && uintpt[ki]->no_of_curves == 0)
397 {
398 sh6settop (uintpt[ki], -1, *(pintpt->left_obj_1), *(pintpt->right_obj_1),
399 *(pintpt->left_obj_2), *(pintpt->right_obj_2), &kstat);
400
401 /* UPDATE (ujk) : Transfer pintpt to main point ?? */
402 /* Mark that an intersection interval is found. */
403 sh6idcon (rintdat, &uintpt[ki], &pintpt, &kstat);
404 if (kstat < 0)
405 goto error;
406 }
407 }
408
409 /* Pre-topology information computed. */
410
411 *jnewpt = kpos;
412 *jstat = 0;
413 goto out;
414
415 /* Error in scratch allocation. */
416
417 err101:*jstat = -101;
418 goto out;
419
420 /* Error in input. Conflicting dimensions. */
421
422 err106:*jstat = -106;
423 goto out;
424
425 /* Error in input. Dimension not equal to 2. */
426
427 err108:*jstat = -108;
428 goto out;
429
430 /* Error lower level routine. */
431
432 error:*jstat = kstat;
433 goto out;
434
435 out:
436 return;
437 }
438
439