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 #define SH1794
43
44 #include "sislP.h"
45
46
47 /*
48 * Forward declarations.
49 * ---------------------
50 */
51
52 #if defined(SISLNEEDPROTOTYPES)
53 static int sh1794_s9findedg(SISLSurf *ps1, SISLSurf *ps2, SISLIntpt **up,
54 int nmb_pt, SISLIntpt **pt1, SISLIntpt **pt2,
55 int *jstat);
56 #else
57 static int sh1794_s9findedg();
58 #endif
59
60 #if defined(SISLNEEDPROTOTYPES)
sh1794(SISLObject * po1,SISLObject * po2,SISLIntpt ** up,int nmb_pt,double aepsge,int * jstat)61 void sh1794(SISLObject *po1, SISLObject *po2, SISLIntpt **up, int nmb_pt,
62 double aepsge, int *jstat)
63 #else
64 void sh1794(po1, po2, up, nmb_pt, aepsge, jstat)
65 SISLObject *po1;
66 SISLObject *po2;
67 SISLIntpt **up;
68 int nmb_pt;
69 double aepsge;
70 int *jstat;
71 #endif
72 /*
73 *********************************************************************
74 *
75 *********************************************************************
76 *
77 * PURPOSE : Check if an intersection curve following a constant
78 * parameter direction between two specified intersection
79 * points represent a tangential intersection curve
80 * between two corners in the intersection problem
81 *
82 *
83 *
84 * INPUT : po1 - First surface in intersection.
85 * po2 - Second surface in intersection
86 * pt1 - Intersection point in the start of the curve
87 * pt2 - Intersection point in the end of the curve
88 * idir - Constant parameter direction (0 and 1 corresponds
89 * aepsge - Geometry resolution.
90 *
91 *
92 * OUTPUT : jstat - status messages
93 * = 2 : Entire surface withing tangential zone
94 * = 1 : Tangential curve found
95 * = 0 : The curve is not tangential or
96 * do not end in corners
97 * < 0 : error
98 *
99 *
100 * METHOD :
101 *
102 *
103 * REFERENCES :
104 *
105 *-
106 * CALLS :
107 *
108 * WRITTEN BY : Vibeke Skytt, SINTEF, 2018-02
109 *
110 *********************************************************************
111 */
112 {
113 int kstat = 0; /* Local status variable */
114 int ki, kj, kr; /* Counters */
115 int kmax = 7; /* Maximum number of iterations to find par. val.*/
116 int kdim; /* Dimension of geometry space */
117 SISLSurf *qs1; /* Surface with constant parameter intersection */
118 SISLSurf *qs2; /* The other surface */
119 int nsample; /* Number of sampling points */
120 int kleft1=0, kleft2=0, kleft3=0, kleft4=0; /* Indices in knot vector */
121 double spar1[2], spar2[2]; /* Parameter value of point on intersection
122 curve with respect to surface */
123 double tstart, tend; /* Endparameters of intersection curve in surface */
124 double tdel, tdel2; /* Step length */
125 int kdir; /* Parameter direction of constant parameter in surface */
126 int kd; /* Parameter direction */
127 double scratch[60]; /* Storage for results of surface evaluation */
128 double *sder1, *sder2, *snorm1, *snorm2; /* Pointers to position and
129 surface normals */
130 double *cvder1, *cvder2;
131 double tang; /* Angle between surface normals */
132 double sstart[2], send[2]; /* Parameter limits of surface */
133 int kk, kn; /* Order and number of coefficients of boundary
134 intersection curve */
135 double *st; /* Knot vector along boundary intersection curve */
136 double angtol = 2.0*ANGULAR_TOLERANCE; /* Needs tuning */
137 double sdist[2]; /* Specified distances between surfaces */
138 double dtol = 0.1*aepsge;/* Tolerance in search for zone par. val.*/
139 double *zonepar = NULL; /* Parameter values corresponding to distances */
140 double *samplepar = NULL; /* Parameter values at samples */
141 double *seg = NULL; /* Segmentation parameters along intersection curve */
142 int nseg = 0; /* Number of segmentation parameters */
143 double ptol = 100*REL_PAR_RES;
144 double tpar; /* Parameter value at the end of the zone */
145 double tprev; /* Previous estimate for parameter value */
146 double tdistprev; /* Distance at previous parameter value */
147 double initdist; /* Distance at initial parameter value */
148 double tdum; /* Estimate for distance at given parameter value */
149
150 double par1[2], par2[2], pos1[3], pos2[3];
151 double tstart2, tend2, dist2;
152 int sgn;
153 double tmin, tmax; /* Allowed endparameters of tangential zone */
154 double td1; /* Distance between intersection points */
155 int failure;
156 int idir;
157 SISLIntpt *pt1, *pt2;
158
159 sdist[0] = 1.2*aepsge;
160 sdist[1] = 2.0*aepsge;
161
162 /* Test input */
163 if (po1->iobj != SISLSURFACE || po2->iobj != SISLSURFACE)
164 goto warn0; /* Not a surface-surface intersection */
165
166 /* Identify edge intersection curve */
167 idir = sh1794_s9findedg(po1->s1, po2->s1, up, nmb_pt, &pt1, &pt2, &kstat);
168 if (kstat < 0)
169 goto error;
170 if (idir < 0 || idir > 3)
171 goto warn0; /* No constant parameter direction specified */
172
173 /* Initiate */
174 *jstat = 0;
175
176 /* Test for corner configuration */
177 sh6isinside(po1, po2, pt1, &kstat);
178 if (kstat < 0)
179 goto error;
180 if (kstat < 3)
181 goto warn0; /* Not a corner */
182
183 sh6isinside(po1, po2, pt2, &kstat);
184 if (kstat < 0)
185 goto error;
186 if (kstat < 3)
187 goto warn0; /* Not a corner */
188
189 /* Evaluate the intersection curve in a number of sampling points
190 to check whether it is tangential */
191 qs1 = (idir <= 1) ? po1->s1 : po2->s1;
192 qs2 = (idir <= 1) ? po2->s1 : po1->s1;
193 kdim = qs1->idim;
194 if (kdim != 3)
195 goto err104;
196
197 /* Prepare for checking */
198 kdir = (idir <= 1) ? idir : idir-2;
199 kd = (idir <= 1) ? 1-idir : 3-kdir;
200 spar1[kdir] = pt1->epar[idir];
201 tstart = pt1->epar[kd];
202 tend = pt2->epar[kd];
203 if (tstart > tend)
204 {
205 tdum = tstart;
206 tstart = tend;
207 tend = tdum;
208 }
209
210 sstart[0] = qs2->et1[qs2->ik1-1];
211 send[0] = qs2->et1[qs2->in1];
212 sstart[1] = qs2->et2[qs2->ik2-1];
213 send[1] = qs2->et2[qs2->in2];
214 spar2[0] = (idir <= 1) ? pt1->epar[2] : pt1->epar[0];
215 spar2[1] = (idir <= 1) ? pt1->epar[3] : pt1->epar[1];
216
217 if (kdir == 0)
218 {
219 tstart2 = qs1->et1[qs1->ik1-1];
220 tend2 = qs1->et1[qs1->in1];
221 }
222 else
223 {
224 tstart2 = qs1->et2[qs1->ik2-1];
225 tend2 = qs1->et1[qs2->in2];
226 }
227 sgn = (spar1[kdir]-tstart < tend-spar1[kdir]) ? 1 : -1;
228
229 /* Set pointers to evaluation results */
230 sder1 = scratch;
231 snorm1 = sder1 + 18;
232 sder2 = snorm1 + 3;
233 snorm2 = sder2 + 18;
234 cvder1 = snorm2 + 3;
235 cvder2 = cvder1 + 9;
236
237 /* Compute number of sampling points */
238 kn = (kdir == 0) ? qs1->in2 : qs1->in1;
239 kk = (kdir == 0) ? qs1->ik2 : qs1->ik1;
240 st = (kdir == 0) ? qs1->et2 : qs1->et1;
241 s1219(st, kk, kn, &kleft1, tstart, &kstat);
242 if (kstat < 0)
243 goto error;
244 s1219(st, kk, kn, &kleft2, tend, &kstat);
245 if (kstat < 0)
246 goto error;
247 nsample = 3*(kleft2 - kleft1)*kk;
248 nsample = max(3*kk, min(nsample, 20*kk));
249 tdel = (tend - tstart)/(double)(nsample-1);
250
251 /* Initial check to avoid overemphasising tiny curves */
252 s1421(qs1, 0, pt1->epar+2*(idir>0), &kleft1, &kleft2, sder1, snorm1, &kstat);
253 if (kstat < 0)
254 goto error;
255 s1421(qs1, 0, pt2->epar+2*(idir>0), &kleft1, &kleft2, sder2, snorm2, &kstat);
256 if (kstat < 0)
257 goto error;
258 td1 = s6dist(sder1, sder2, kdim);
259 if (td1 < aepsge && (DNEQUAL(tstart,st[kk-1]) || DNEQUAL(tend,st[kn])))
260 goto warn0;
261
262 nsample = max(2, min(nsample, (int)(td1/aepsge)));
263
264 zonepar = new0array(2*nsample, DOUBLE);
265 samplepar = newarray(nsample, DOUBLE);
266 seg = newarray(2+nsample, DOUBLE);
267 if (zonepar == NULL || samplepar == NULL || seg == NULL)
268 goto err101;
269
270 for (ki=0, spar1[1-kdir]=tstart; ki<nsample; ++ki, spar1[1-kdir]+=tdel)
271 {
272 samplepar[ki] = spar1[1-kdir];
273
274 /* Test if the two surfaces intersect tangentially in the
275 current sampling point.
276 First evaluate the surface in which there is an intersection
277 curve along the boundary */
278 s1421(qs1, 2, spar1, &kleft1, &kleft2, sder1, snorm1, &kstat);
279 if (kstat < 0)
280 goto error;
281
282 /* Find the closest point in the other surface. */
283 s1775(qs2, sder1, kdim, aepsge, sstart, send, spar2, spar2, &kstat);
284 if (kstat < 0)
285 goto error;
286
287 /* Evaluate the second surface. */
288 s1421(qs2, 2, spar2, &kleft3, &kleft4, sder2, snorm2, &kstat);
289 if (kstat < 0)
290 goto error;
291
292 // Test if the current point is a touch point
293 tang = s6ang(snorm1, snorm2, kdim);
294 if (tang > angtol)
295 break;
296
297 /* Compute 0-2. derivative of corresponding curves in surfaces */
298 /* Surface with boundary intersection */
299 memmove(cvder1, sder1, kdim*sizeof(double));
300 memmove(cvder1+kdim, sder1+(1+kdir)*kdim, kdim*sizeof(double));
301 memmove(cvder1+2*kdim, sder1+(3+2*kdir)*kdim, kdim*sizeof(double));
302
303 /* The other surface */
304 s1291(cvder1, sder2, kdim, cvder2, &kstat);
305 if (kstat < 0)
306 goto error;
307
308 /* Compute the parameter value where two curves have the
309 specified distances (approximation) */
310 tprev = 0.0;
311 initdist = tdistprev = s6dist(cvder1, cvder2, kdim);
312 tdum = s6dist(cvder1+2*kdim, cvder2+2*kdim, kdim);
313
314 for (kj=0; kj<2; ++kj)
315 {
316 if (tdum > REL_PAR_RES)
317 {
318 tpar = 2*sdist[kj]/tdum;
319 tpar = sqrt(tpar);
320 }
321 else
322 tpar = 0.0;
323
324 /* Iterate to a better position */
325 /* Could also take the parameter value of the previous
326 sample point as input */
327 failure = 0;
328 for (kr=0; kr<kmax; ++kr)
329 {
330 /* Check accuracy of current estimate */
331 par1[1-kdir] = spar1[1-kdir];
332 par1[kdir] = spar1[kdir] + sgn*tpar;
333 par2[0] = spar2[0];
334 par2[1] = spar2[1];
335 s1424(qs1, 0, 0, par1, &kleft1, &kleft2, pos1, &kstat);
336 if (kstat < 0)
337 goto error;
338 s1775(qs2, pos1, kdim, aepsge, sstart, send, par2, par2, &kstat);
339 if (kstat < 0)
340 goto error;
341 s1424(qs2, 0, 0, par2, &kleft3, &kleft4, pos2, &kstat);
342 if (kstat < 0)
343 goto error;
344 dist2 = s6dist(pos1, pos2, kdim);
345 /*printf("Dist: %7.13f, expected: %7.13f \n",dist2, sdist[kj]); */
346 if (fabs(dist2 - sdist[kj]) < dtol)
347 {
348 tprev = tpar;
349 break; /* The estimate for the tangential zone limit
350 is close enough */
351 }
352 if (DEQUAL(dist2, tdistprev) || dist2 < initdist)
353 {
354 failure = 1;
355 break;
356 }
357 if (fabs(dist2 - sdist[kj]) > fabs(tdistprev - sdist[kj]))
358 {
359 tpar = 0.5*(tprev+tpar); // We are loosing accuracy
360 continue;
361 }
362
363 /* Find new position */
364 tdel2 = (sdist[kj]-dist2)*(tpar-tprev)/(dist2-tdistprev);
365 tprev = tpar;
366 tdistprev = dist2;
367 tpar += tdel2;
368 }
369 if (failure)
370 break;
371 zonepar[2*ki+kj] = tprev;
372 }
373
374 if (failure)
375 break;
376
377 if (zonepar[2*ki] > zonepar[2*ki+1])
378 {
379 tdum = zonepar[2*ki];
380 zonepar[2*ki] = zonepar[2*ki+1];
381 zonepar[2*ki+1] = tdum;
382 }
383 }
384
385 if (ki < nsample)
386 {
387 /* Not a tangential intersection curve */
388 *jstat = 0;
389 goto out;
390 }
391
392 /* Identify non-overlapping pairs of tangential zone end parameters */
393 /* Check also whether the tangential intersection curve follows the
394 entire surface boundary */
395 if (tstart > st[kk-1] + ptol)
396 seg[nseg++] = tstart;
397
398 for (ki=0; ki<nsample; ki=kj)
399 {
400 for (kj=ki+1; kj<nsample; ++kj)
401 {
402 if (zonepar[2*ki] > zonepar[2*kj+1] ||
403 zonepar[2*ki+1] < zonepar[2*kj])
404 {
405 /* No overlap. Define segmentation parameter */
406 /* First look for a suitable knot value */
407 s1219(st, kk, kn, &kleft1, samplepar[ki], &kstat);
408 if (kstat < 0)
409 goto error;
410 s1219(st, kk, kn, &kleft2, samplepar[kj], &kstat);
411 if (kstat < 0)
412 goto error;
413 if (samplepar[kj] > st[kleft2]+ptol && kleft2 > kleft1)
414 {
415 seg[nseg++] = st[kleft2];
416 for (kr=kj-1; kr>ki && samplepar[kr]>st[kleft2]; --kr);
417 kj = max(kr, ki+1);
418 }
419 else if (kj > ki+1)
420 {
421 seg[nseg++] = samplepar[kj-1];
422 --kj;
423 }
424 else
425 seg[nseg++] = 0.5*(samplepar[ki]+samplepar[kj]);
426
427 break;
428 }
429 }
430 }
431
432 if (tend < st[kn] - ptol)
433 seg[nseg++] = tend;
434
435 if (nseg > 0)
436 {
437 /* Define segmentation parameters along the tangential
438 intersection curve */
439 sh6setseg(qs1, 1-kdir, seg, nseg, LIMITING_SEG, &kstat);
440 if (kstat < 0)
441 goto error;
442 }
443 else
444 {
445 /* Set width of tangential zone */
446 tmin = zonepar[0];
447 tmax = zonepar[1];
448 for (ki=1; ki<nsample; ++ki)
449 {
450 tmin = max(tmin, zonepar[2*ki]);
451 tmax = min(tmax, zonepar[2*ki+1]);
452 }
453 tpar = spar1[kdir] + 0.5*sgn*(tmin + tmax);
454
455 if (tpar < tstart2 || tpar > tend2)
456 {
457 /* The entire surface is within the tangential intersectin
458 zone. Set output status */
459 *jstat = 2;
460 goto out;
461 }
462
463 /* Define segmentation parameter perpendicular to the tangential
464 intersection curve */
465 seg[nseg++] = tpar;
466 sh6setseg(qs1, kdir, seg, nseg, (sgn == 1) ? TANGENTIAL_BELT_LEFT :
467 TANGENTIAL_BELT_RIGHT, &kstat);
468 if (kstat < 0)
469 goto error;
470 }
471
472 *jstat = 1;
473 goto out;
474
475 warn0:
476 /* Not a tangential curve */
477 *jstat = 0;
478 goto out;
479
480 error:
481 /* Error in lower order function */
482 *jstat = kstat;
483 goto out;
484
485 err101:
486 /* Error in scratch allocation */
487 *jstat = -101;
488 goto out;
489
490 err104:
491 /* Dimension of geometry space not equal to 3 */
492 *jstat = -104;
493 goto out;
494
495 out:
496 if (zonepar != NULL) freearray(zonepar);
497 if (samplepar != NULL) freearray(samplepar);
498 if (seg != NULL) freearray(seg);
499
500 return;
501 }
502
503
504 #if defined(SISLNEEDPROTOTYPES)
sh1794_s9findedg(SISLSurf * ps1,SISLSurf * ps2,SISLIntpt ** up,int nmb_pt,SISLIntpt ** pt1,SISLIntpt ** pt2,int * jstat)505 static int sh1794_s9findedg(SISLSurf *ps1, SISLSurf *ps2, SISLIntpt **up,
506 int nmb_pt, SISLIntpt **pt1, SISLIntpt **pt2,
507 int *jstat)
508 #else
509 static int sh1794_s9findedg(ps1, ps2, up, nmb_pt, pt1, pt2, jstat)
510 SISLSurf *ps1;
511 SISLSurf *ps2;
512 SISLIntpt **up;
513 int nmb_pt;
514 SISLIntpt **pt1;
515 SISLIntpt **pt2;
516 int *jstat;
517 #endif
518
519 /*
520 *********************************************************************
521 *
522 *********************************************************************
523 *
524 * PURPOSE : Check if a chain of linked intersection points
525 * contain an intersection curve along some surface edge
526 * and select the most signicant candidate
527 *
528 *
529 *
530 * INPUT : ps1 - First surface in intersection.
531 * ps2 - Second surface in intersection
532 * up - Chain of intersection points
533 * nmb_pt - Number of intersection points
534 *
535 *
536 * OUTPUT : return value - Parameter direction of detected edge (both sfs)
537 * pt1 - Intersection point in the start of the curve
538 * pt2 - Intersection point in the end of the curve
539 * jstat - status messages
540 * = 0 : OK
541 * < 0 : error
542 *
543 *
544 * METHOD :
545 *
546 *
547 * REFERENCES :
548 *
549 *-
550 * CALLS :
551 *
552 * WRITTEN BY : Vibeke Skytt, SINTEF, 2018-02
553 *
554 *********************************************************************
555 */
556 {
557 int kstat = 0; /* Local status variable */
558 int kdir; /* Parameter direction of possible tangential intersection curve */
559 int kdir2;
560 int pdir[4]; /* Maximum possible parameter directions */
561 int ix[8]; /* Maximum possible end parameters of constant chain */
562 int ki, kj; /* Counters */
563 int kidx = 0; /* Current parameter direction */
564 int klist1=0, klist2=0;
565 double parlen1, parlen2; /* Distance between intersection points in
566 parameter domain along constant direction */
567 int p1, p2; /* Parameter directions */
568 double del1, del2; /* Edge length of surface in parameter domain */
569
570 pdir[0] = pdir[1] = pdir[2] = pdir[3] = -1;
571 ix[2*kidx] = 0;
572 for (kj=0; kj<nmb_pt; kj=max(kj+1,ki-1))
573 {
574 kdir = -1;
575 for (ki=kj+1; ki<nmb_pt; ++ki)
576 {
577 sh6getlist(up[ki-1], up[ki], &klist1, &klist2, &kstat);
578 if (kstat < 0)
579 goto error;
580
581 for (kdir2=0; kdir2<up[ki-1]->ipar; ++kdir2)
582 {
583 if (up[ki-1]->curve_dir[klist1] & (1 << (kdir2+1)))
584 break;
585 }
586 ix[2*kidx] = kj;
587 if (kdir2 == up[ki-1]->ipar)
588 {
589 if (kdir >= 0)
590 {
591 pdir[kidx] = kdir;
592 ix[2*kidx+1] = ki-1;
593 kidx++;
594 kdir = 0;
595 }
596 kdir2 = -1;
597 break;
598 }
599 else if (kdir >= 0 && kdir2 != kdir)
600 {
601 pdir[kidx] = kdir;
602 ix[2*kidx+1] = ki;
603 kidx++;
604 kdir = kdir2;
605 break;
606 }
607 else
608 {
609 kdir = kdir2;
610 }
611 }
612 }
613 if (kdir2 >= 0 && (kidx == 0 || kdir2 != pdir[kidx-1]))
614 {
615 pdir[kidx] = kdir2;
616 ix[2*kidx+1] = ki;
617 kidx++;
618 }
619
620 kdir = -1;
621 if (kidx > 0)
622 {
623 kdir = pdir[0];
624 p1 = (kdir <= 1) ? 1 - kdir : 5 - kdir;
625 parlen1 = fabs(up[ix[1]-1]->epar[p1]-up[ix[0]]->epar[p1]);
626 if (kdir <= 0)
627 del1 = (p1 == 0) ? ps1->et1[ps1->in1] - ps1->et1[ps1->ik1-1] :
628 ps1->et2[ps1->in2] - ps1->et2[ps1->ik2-1];
629 else
630 del1 = (p1 == 0) ? ps2->et1[ps2->in1] - ps2->et1[ps2->ik1-1] :
631 ps2->et2[ps2->in2] - ps2->et2[ps2->ik2-1];
632
633 (*pt1) = up[ix[0]];
634 (*pt2) = up[ix[1]-1];
635 for (ki=1; ki<kidx; ++ki)
636 {
637 /* Check if other edge intersections are more significant */
638 kdir2 = pdir[ki];
639 p2 = (kdir2 <= 1) ? 1 - kdir2 : 5 - kdir2;
640 parlen2 = fabs(up[ix[2*ki+1]-1]->epar[p2]-up[ix[2*ki]]->epar[p2]);
641 if (kdir2 <= 0)
642 del2 = (p2 == 0) ? ps1->et1[ps1->in1] - ps1->et1[ps1->ik1-1] :
643 ps1->et2[ps1->in2] - ps1->et2[ps1->ik2-1];
644 else
645 del2 = (p2 == 0) ? ps2->et1[ps2->in1] - ps2->et1[ps2->ik1-1] :
646 ps2->et2[ps2->in2] - ps2->et2[ps2->ik2-1];
647 if (parlen2/del2 > parlen1/del1)
648 {
649 kdir = kdir2;
650 (*pt1) = up[ix[2*ki]];
651 (*pt2) = up[ix[2*ki+1]-1];
652 parlen1 = parlen2;
653 del1 = del2;
654 }
655 }
656 }
657
658 *jstat = 0;
659 goto out;
660
661 error:
662 *jstat = kstat;
663 goto out;
664
665 out:
666 return kdir;
667 }
668
669