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: sh1781.c,v 1.2 2001-03-19 15:59:05 afr Exp $
45 *
46 */
47
48
49 #define SH1781
50
51 #include "sislP.h"
52
53
54 #if defined(SISLNEEDPROTOTYPES)
55 void
sh1781(SISLObject * po1,SISLObject * po2,double aepsge,SISLIntdat ** rintdat,SISLIntpt * pintpt,int * jnewpt,int * jstat)56 sh1781 (SISLObject * po1, SISLObject * po2, double aepsge,
57 SISLIntdat ** rintdat, SISLIntpt * pintpt, int *jnewpt,
58 int *jstat)
59 #else
60 void
61 sh1781 (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 /* UPDATE: (ujk) : Must tune the test on when to use local
71 info (=subtask no ?) */
72 /*
73 ******************************************************************
74 *
75 *********************************************************************
76 *
77 * PURPOSE : Set pre-topology data in 1-dimensional curve-point
78 * intersection. Also find help points if necessary.
79 *
80 *
81 * INPUT : po1 - Pointer to the first object in the intersection.
82 * po2 - Pointer to the second object in the intersection.
83 * aepsge - Geometry resolution.
84 * rintdat - Intersection data structure.
85 * pintpt - Current intersection point.
86 *
87 *
88 * OUTPUT : pintpt - Intersection point after updating pre-topology.
89 * jnewpt - Number of new int.pt. created.
90 * jstat - status messages
91 * > 0 : Warning.
92 * = 0 : Ok.
93 * < 0 : Error.
94 *
95 *
96 * METHOD :
97 *
98 * CALLS : shevalc - Evaluate curve using filtered coefficients.
99 * s6ang - Angle between two vectors.
100 * s6idcon - Connect two intersection points.
101 * hp_newIntpt - Create new intersection point.
102 *
103 * REFERENCES :
104 *
105 * WRITTEN BY : Vibeke Skytt, SI, Oslo, Norway. 04.91.
106 * Ulf J. Krystad, SI, 06.91.
107 * REWISED BY : VSK, 11-92. Make sure that a help point always is an
108 * intersection point.
109 *********************************************************************
110 */
111 {
112 int kstat = 0; /* Status variable. */
113 int ki, kj; /* Counters. */
114 int kleft = 0; /* Parameter to evaluator. */
115 int korgleft = 0; /* Knot index. */
116 int kdim; /* Dimension of geometry space. */
117 int kn; /* Number of vertices of curve. */
118 int kk; /* Order of curve. */
119 int kpos = 0; /* Current position in int.pt. array. */
120 int lleft[2]; /* Array storing pre-topology information. */
121 int lright[2]; /* Array storing pre-topology information. */
122 int *ll1, *ll2, *lr1, *lr2; /* Pointers into pre-topology arrays. */
123 double tpoint; /* Level value. */
124 double tpar0,tpar; /* Parameter value of point on curve. */
125 double spar[1]; /* Parameter value of endpoint of curve. */
126 double sder[2]; /* Result of curve evaluation. */
127 double stang1[2]; /* Tangent vector of curve. */
128 double stang2[2]; /* Tangent vector of level value. */
129 double *st; /* Pointer to knot vector of curve. */
130 double *sptpar = pintpt->epar;/* Pointer to parameter array of int.pt. */
131 double tref; /* Referance value in equality test. */
132 SISLCurve *qc; /* Pointer to current curve. */
133 SISLIntpt *uintpt[2]; /* Array storing new intersection points. */
134 double *ret_val; /* Pointer to geo data from sh6getgeom */
135 double *ret_norm; /* Pointer to geo data from sh6getgeom */
136 double *nullp = SISL_NULL;
137 int make_hp; /* Flag, make/not make help pt. */
138
139 /* Don't make pretop for help points ! */
140 if (sh6ishelp (pintpt))
141 {
142 *jstat = 0;
143 goto out;
144 }
145
146 /* Set pointers into the arrays storing pre-topology information. */
147
148 if (po1->iobj == SISLCURVE)
149 {
150 ll1 = lleft;
151 lr1 = lright;
152 ll2 = lleft + 1;
153 lr2 = lright + 1;
154 }
155 else
156 {
157 ll1 = lleft + 1;
158 lr1 = lright + 1;
159 ll2 = lleft;
160 lr2 = lright;
161 }
162
163 /* Get pre-topology information. */
164 sh6gettop (pintpt, -1, lleft, lright, lleft + 1, lright + 1, &kstat);
165 if (kstat < 0)
166 goto error;
167
168 /* Test dimension of geometry space. */
169 if (po1->iobj == SISLCURVE)
170 {
171 qc = po1->c1;
172 }
173 else
174 {
175 qc = po2->c1;
176 }
177
178 kdim = qc->idim;
179 if (kdim != 1)
180 goto err106;
181
182 /* Store curve information in local parameters. */
183
184 kn = qc->in;
185 kk = qc->ik;
186 st = qc->et;
187 tref = st[kn] - st[kk - 1];
188
189 /* Fetch geometry information, point. */
190 sh6getgeom ((po1->iobj == SISLPOINT) ? po1 : po2,
191 (po1->iobj == SISLPOINT) ? 1 : 2,
192 pintpt, &ret_val, &ret_norm, aepsge, &kstat);
193 if (kstat < 0)
194 goto error;
195
196 tpoint = ret_val[0];
197
198 /* Fetch geometry information, curve. */
199 sh6getgeom ((po1->iobj == SISLCURVE) ? po1 : po2,
200 (po1->iobj == SISLCURVE) ? 1 : 2,
201 pintpt, &ret_val, &ret_norm, aepsge, &kstat);
202 if (kstat < 0)
203 goto error;
204
205 s1219(st,kk,kn,&korgleft,sptpar[0],&kstat);
206 if (kstat < 0) goto error;
207
208 sder[0] = ret_val[0];
209 sder[1] = ret_val[1];
210
211 /* Set tangent vectors. */
212
213 stang1[0] = (double) 1.0;
214 stang1[1] = ret_val[1];
215 stang2[0] = (double) 1.0;
216 stang2[1] = DZERO;
217
218 /* UPDATE (ujk) : tune */
219 if (s6ang (stang1, stang2, 2) > 0.001*ANGULAR_TOLERANCE)
220 {
221 /* Compute pre-topology using local information. */
222
223 if (sder[1] > 0)
224 {
225 *ll1 = SI_IN;
226 *lr1 = SI_OUT;
227 *ll2 = SI_OUT;
228 *lr2 = SI_IN;
229 }
230 else
231 {
232 *ll1 = SI_OUT;
233 *lr1 = SI_IN;
234 *ll2 = SI_IN;
235 *lr2 = SI_OUT;
236 }
237
238 }
239 else
240 {
241 /* Test if the intersection point lies at the endpoint of
242 the curve. */
243
244 if (DEQUAL (sptpar[0] + tref, st[kn] + tref))
245 {
246
247 }
248 else
249 {
250 /* Find endpoint of coincidence interval in the positive
251 direction of the curve. */
252
253 ki = 0;
254 tpar = sptpar[0] + (double) 2.0 *sqrt (aepsge);
255 tpar = min (tpar, st[kn]);
256 tpar0 = tpar = min (tpar, st[korgleft+1]);
257 shevalc (qc, 0, tpar, aepsge, &kleft, sder, &kstat);
258 if (fabs (sder[0] - tpoint) <= aepsge)
259 {
260 make_hp = TRUE;
261 for (ki = kleft - kk + 1; ki < kn; ki++)
262 {
263 for (tpar = DZERO, kj = ki + 1; kj < ki + kk; kj++)
264 tpar += st[kj];
265 tpar /= (double) (kk - 1);
266
267 if (tpar > sptpar[0] && DNEQUAL(tpar,sptpar[0]))
268 {
269 shevalc (qc, 0, tpar, aepsge, &kleft, sder, &kstat);
270 if (fabs (sder[0] - tpoint) >= aepsge)
271 break;
272
273 tpar0 = tpar; /* Remember parameter value. */
274 }
275 }
276 }
277 /*UJK, sept 92, don't make help pt close to main */
278 else make_hp = FALSE;
279
280 /* Test if there is coincidence along the entire curve part. */
281
282 if (ki == kn)
283 {
284 /* Set right values of original point. */
285 *lr1 = *lr2 = SI_ON;
286 }
287 else
288 {
289 /* Compute right values of intersection point. */
290 *lr1 = (sder[0] > tpoint) ? SI_OUT : SI_IN;
291 *lr2 = (*lr1 == SI_IN) ? SI_OUT : SI_IN;
292
293 /*UJK, sept 92, don't make help pt close to main */
294 if (make_hp)
295 {
296 /* Create help point. */
297 if (sptpar[0] < st[kleft])
298 spar[0] = MIN(tpar0,st[kleft]);
299 else
300 spar[0] = tpar0;
301
302 uintpt[kpos] = SISL_NULL;
303 if ((uintpt[kpos] = hp_newIntpt (1, spar, DZERO, -SI_ORD,
304 SI_ON, lright[0], SI_ON,
305 lright[1], 0, 0, nullp, nullp)) == SISL_NULL)
306 goto err101;
307
308 /* Insert the point into the data structure. */
309
310 sh6idnpt (rintdat, &uintpt[kpos], 1, &kstat);
311 if (kstat < 0)
312 goto error;
313
314 kpos++;
315 }
316 }
317 }
318
319 /* Test if the intersection point lies at the startpoint
320 of the curve. */
321
322 if (DEQUAL (sptpar[0] + tref, st[kk - 1] + tref))
323 {
324 }
325 else
326 {
327 /* Find endpoint of coincidence interval in the negative
328 direction of the curve. */
329
330 ki = kn;
331 while (sptpar[0] == st[korgleft]) korgleft--;
332 tpar = sptpar[0] - (double) 2.0 *sqrt (aepsge);
333 tpar = max (tpar, st[kk - 1]);
334 tpar0 = tpar = max (tpar, st[korgleft]);
335 shevalc (qc, 0, tpar, aepsge, &kleft, sder, &kstat);
336 if (fabs (sder[0] - tpoint) <= aepsge)
337 {
338 make_hp = TRUE;
339 for (ki = kleft; ki >= 0; ki--)
340 {
341 for (tpar = DZERO, kj = ki + 1; kj < ki + kk; kj++)
342 tpar += st[kj];
343 tpar /= (double) (kk - 1);
344
345 if (tpar < sptpar[0] && DNEQUAL(tpar,sptpar[0]))
346 {
347 shevalc (qc, 0, tpar, aepsge, &kleft, sder, &kstat);
348 if (fabs (sder[0] - tpoint) >= aepsge)
349 break;
350
351 tpar0 = tpar;
352 }
353 }
354 }
355 /*UJK, sept 92, don't make help pt close to main */
356 else make_hp = FALSE;
357
358 /* Test if there is coincidence along the entire curve part. */
359 if (ki < 0)
360 {
361 /* Set left values of original point. */
362 *ll1 = *ll2 = SI_ON;
363 }
364 else
365 {
366 /* Compute left values of intersection point. */
367
368 *ll1 = (sder[0] > tpoint) ? SI_OUT : SI_IN;
369 *ll2 = (*ll1 == SI_IN) ? SI_OUT : SI_IN;
370
371 /*UJK, sept 92, don't make help pt close to main */
372 if (make_hp)
373 {
374 /* Create intersection point. */
375 if (sptpar[0] > st[kleft+1])
376 spar[0] = MAX(tpar0,st[kleft+1]);
377 else
378 spar[0] = tpar0;
379
380 uintpt[kpos] = SISL_NULL;
381 if ((uintpt[kpos] = hp_newIntpt (1, spar, DZERO, -SI_ORD,
382 lleft[0], SI_ON, lleft[1],
383 SI_ON, 0, 0, nullp, nullp)) == SISL_NULL)
384 goto err101;
385
386 /* Insert the point into the data structure. */
387
388 sh6idnpt (rintdat, &uintpt[kpos], 1, &kstat);
389 if (kstat < 0)
390 goto error;
391
392
393 kpos++;
394 }
395
396 }
397
398 }
399 }
400
401 /* Update pretopology of intersection point. */
402
403 sh6settop (pintpt, -1, lleft[0], lright[0], lleft[1], lright[1], &kstat);
404 if (kstat < 0)
405 goto error;
406 /* Change, if necessary, pintpt to mainpoint */
407 sh6tomain (pintpt, &kstat);
408
409 /* Join intersection points. (kpos=0,1,2)*/
410 for (ki = 0; ki < kpos; ki++)
411 {
412 sh6idnpt (rintdat, &uintpt[ki], 1, &kstat);
413 if (kstat < 0)
414 goto error;
415 /* Mark that an intersection interval is found. */
416 if (sh6ishelp (uintpt[ki]) && uintpt[ki]->no_of_curves == 0)
417 {
418 sh6idcon (rintdat, &uintpt[ki], &pintpt, &kstat);
419 if (kstat < 0)
420 goto error;
421 }
422 }
423
424 /* Pre-topology information computed. */
425
426 *jnewpt = kpos;
427 *jstat = 0;
428 goto out;
429
430 /* Error in scratch allocation. */
431
432 err101:*jstat = -101;
433 goto out;
434
435 /* Error in input. Incorrect dimension. */
436
437 err106:*jstat = -106;
438 goto out;
439
440 /* Error lower level routine. */
441
442 error:*jstat = kstat;
443 goto out;
444
445 out:
446 return;
447 }
448