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: s1741.c,v 1.3 2001-03-19 15:58:53 afr Exp $ 45 * 46 */ 47 48 49 #define S1741 50 51 #include "sislP.h" 52 53 #if defined(SISLNEEDPROTOTYPES) 54 void s1741(SISLObject * po1,SISLObject * po2,double aepsge,int * jstat)55 s1741(SISLObject *po1,SISLObject *po2,double aepsge,int *jstat) 56 #else 57 void s1741(po1,po2,aepsge,jstat) 58 SISLObject *po1; 59 SISLObject *po2; 60 double aepsge; 61 int *jstat; 62 #endif 63 /* 64 ********************************************************************* 65 * 66 ********************************************************************* 67 * 68 * PURPOSE : Check if a object-object intersection is a simple case, 69 * i.e. the intersection will result in one single point. 70 * 71 * 72 * 73 * INPUT : po1 - First curve in the intersection problem. 74 * po2 - Second curve in the intersection problem. 75 * aepsge - Geometry resolution. 76 * 77 * 78 * OUTPUT : jstat - status messages 79 * = 1 : simpel case. 80 * = 0 : not simpel case. 81 * < 0 : error. 82 * 83 * 84 * METHOD : 85 * 86 * 87 * REFERENCES : 88 * 89 * CALLS : s1990 - Making the direction cone of surface. 90 * s1991 - Making the direction cone of curve. 91 * sh1993 - Simple case of curve in one dimention. 92 * sh1994 - Simple case of surface in one dimention. 93 * s6err - Gives error message. 94 * 95 * WRITTEN BY : Arne Laksaa, SI, 89-05. 96 * REVISED BY : Michael Floater, SI, May 1997 to agree with the HP version. 97 * 98 ********************************************************************* 99 */ 100 { 101 int kstat = 0; /* Local status variable. */ 102 int kpos = 0; /* Position of the error. */ 103 int k1; /* Control variable in loop. */ 104 double tang; /* Angel between two vectors. */ 105 double small_tang;/* Smallest angle between two vectors. */ 106 107 if (po1->iobj == SISLPOINT || po2->iobj == SISLPOINT) 108 { 109 SISLObject *qo1,*qo2; 110 111 if(po1->iobj == SISLPOINT) 112 { 113 qo1 = po1; 114 qo2 = po2; 115 } 116 else 117 { 118 qo1 = po2; 119 qo2 = po1; 120 } 121 122 if (qo2->iobj == SISLCURVE) 123 { 124 /* Test if the curve lies in the same space as the point. */ 125 126 if (qo1->p1->idim != qo2->c1->idim) goto err106; 127 128 if (qo2->c1->idim == 1) 129 { 130 sh1993(qo2->c1,aepsge,&kstat); 131 132 *jstat = kstat; 133 goto out; 134 } 135 136 /* Computing the direction cone of the curve. If the curve 137 have cones greater then pi we just return not a simple case. */ 138 139 s1991(qo2->c1,aepsge,&kstat); 140 if (kstat < 0) goto error; 141 else if (qo2->c1->pdir->igtpi != 0) goto out2;/* Not a simple case.*/ 142 143 144 /* Performing a simple case check. */ 145 146 if (qo2->c1->pdir->aang<PIHALF) 147 { 148 /* A simpel case. The iteration is able to 149 find intersection.*/ 150 151 *jstat = 1; 152 goto out; 153 } 154 } 155 else if (qo2->iobj == SISLSURFACE) 156 { 157 /* Test if the surface lies in the same space as the point. */ 158 159 if (qo1->p1->idim != qo2->s1->idim) goto err106; 160 161 162 if (qo2->s1->idim == 1) 163 { 164 sh1994(qo2->s1,aepsge,&kstat); 165 166 *jstat = kstat; 167 goto out; 168 } 169 else 170 { 171 /* Computing the direction cone of the surface. If the surface 172 have cones greater then pi we just return not a simple case.*/ 173 174 s1990(qo2->s1,aepsge,&kstat); 175 if (kstat < 0) goto error; 176 else if (qo2->s1->pdir->igtpi != 0) goto out2; /*No simple case*/ 177 178 179 /* Performing a simple case check. */ 180 181 if (qo2->s1->pdir->aang<PIHALF) 182 { 183 /* A simpel case. The iteration is able to 184 find intersection.*/ 185 186 187 *jstat = 1; 188 goto out; 189 } 190 } 191 } 192 } 193 else if (po1->iobj == SISLCURVE && po2->iobj == SISLCURVE) 194 { 195 /* Test if the curves lies in the same space. */ 196 197 if (po2->c1->idim != po1->c1->idim) goto err106; 198 199 200 201 /* Computing the direction cone of the two curves. If one of them 202 have cones greater then pi we just return not a simple case. */ 203 204 s1991(po1->c1,aepsge,&kstat); 205 if (kstat < 0) goto error; 206 207 s1991(po2->c1,aepsge,&kstat); 208 if (kstat < 0) goto error; 209 210 if (po1->c1->pdir->igtpi != 0) goto out2; /* Not a simple case.*/ 211 if (po2->c1->pdir->igtpi != 0) goto out2; /* Not a simple case.*/ 212 213 214 /* Computing the angle beetween the senters of the two cones. */ 215 216 for (tang=DZERO,k1=0;k1<po1->c1->idim;k1++) 217 tang += po1->c1->pdir->ecoef[k1]*po2->c1->pdir->ecoef[k1]; 218 219 if (tang >= DZERO) tang = min((double)1.0,tang); 220 else tang = max((double)-1.0,tang); 221 222 tang = acos(tang); 223 224 if (tang > PIHALF) 225 small_tang = PI - tang; 226 else 227 small_tang = tang; 228 229 /* Performing a simple case check. */ 230 231 if ((tang+po1->c1->pdir->aang+po2->c1->pdir->aang)<PI && 232 (po1->c1->pdir->aang+po2->c1->pdir->aang)<tang) 233 { 234 /* A simpel case. The two cones and their mirrors 235 are not intersecting.*/ 236 237 *jstat = 1; 238 goto out; 239 } 240 else if (po1->c1->idim == 2) 241 { 242 *jstat = 0; 243 goto out; 244 } 245 else if (tang < PI - ANGULAR_TOLERANCE && 246 tang > ANGULAR_TOLERANCE && 247 po1->c1->pdir->aang <= (double)1.3*small_tang && 248 po2->c1->pdir->aang <= (double)1.3*small_tang) 249 /*po1->c1->pdir->aang <= (double)1.3*tang && 250 po2->c1->pdir->aang <= (double)1.3*tang)*/ 251 { 252 s1796(po1->c1,po2->c1,aepsge,tang,&kstat); 253 if (kstat<0) goto error; 254 else *jstat = kstat; 255 goto out; 256 } 257 } 258 else if (po1->iobj == SISLSURFACE && po2->iobj == SISLSURFACE) 259 { 260 261 /* Test if the surfaces lies in the same space. */ 262 263 if (po2->s1->idim != po1->s1->idim) goto err106; 264 265 266 267 /* Computing the direction cone of the two surfaces. If one of them 268 have cones greater then pi we just return not a simple case. */ 269 270 s1990(po1->s1,aepsge,&kstat); 271 if (kstat < 0) goto error; 272 273 s1990(po2->s1,aepsge,&kstat); 274 if (kstat < 0) goto error; 275 276 if (po1->s1->pdir->igtpi != 0) goto out2; /* Not a simple case. */ 277 278 if (po2->s1->pdir->igtpi != 0) goto out2; /* Not a simple case. */ 279 280 /* Computing the angle beetween the senters of the two cones. */ 281 282 for (tang=DZERO,k1=0;k1<po1->s1->idim;k1++) 283 tang += po1->s1->pdir->ecoef[k1]*po2->s1->pdir->ecoef[k1]; 284 285 if (tang >= DZERO) tang = min((double)1.0,tang); 286 else tang = max((double)-1.0,tang); 287 288 tang = acos(tang); 289 290 291 /* Performing a simple case check. */ 292 293 if ((tang+po1->s1->pdir->aang+po2->s1->pdir->aang)<PI && 294 (po1->s1->pdir->aang+po2->s1->pdir->aang)<tang) 295 { 296 /* A simpel case. The two cones and their mirrors 297 are not intersecting.*/ 298 299 po1->psimple = po2; 300 *jstat = 1; 301 goto out; 302 } 303 else if (tang < PI - ANGULAR_TOLERANCE && 304 tang > ANGULAR_TOLERANCE && 305 po1->s1->pdir->aang <= (double)1.3*tang && 306 po2->s1->pdir->aang <= (double)1.3*tang) 307 { 308 s1795(po1->s1,po2->s1,aepsge,tang,&kstat); 309 if (kstat < 0) goto error; 310 if (kstat == 1) po1->psimple = po2; 311 *jstat = kstat; 312 goto out; 313 } 314 } 315 else if (po1->iobj == SISLCURVE || po2->iobj == SISLCURVE) 316 { 317 SISLObject *qo1,*qo2; 318 319 if(po1->iobj == SISLCURVE) 320 { 321 qo1 = po1; 322 qo2 = po2; 323 } 324 else 325 { 326 qo1 = po2; 327 qo2 = po1; 328 } 329 330 331 /* Test if the surface and curve lies in the same space. */ 332 333 if (qo2->s1->idim != qo1->c1->idim) goto err106; 334 335 336 337 /* Computing the direction cone of the curve and the surface. If one of 338 them have cones greater then pi we just return not a simple case. */ 339 340 341 s1990(qo2->s1,aepsge,&kstat); 342 if (kstat < 0) goto error; 343 344 s1991(qo1->c1,aepsge,&kstat); 345 if (kstat < 0) goto error; 346 347 if (qo1->c1->pdir->igtpi != 0) goto out2; /* Not a simple case. */ 348 if (qo2->s1->pdir->igtpi != 0) goto out2; /* Not a simple case. */ 349 350 351 352 /* Computing the angle beetween the senters of the two cones. */ 353 354 for (tang=DZERO,k1=0;k1<qo2->s1->idim;k1++) 355 tang += qo2->s1->pdir->ecoef[k1]*qo1->c1->pdir->ecoef[k1]; 356 357 if (tang >= DZERO) tang = min((double)1.0,tang); 358 else tang = max((double)-1.0,tang); 359 360 tang = acos(tang); 361 362 363 /* Performing a simple case check. */ 364 365 if (((tang + qo1->c1->pdir->aang) < (PIHALF - qo2->s1->pdir->aang)) || 366 ((tang - PIHALF - qo1->c1->pdir->aang) > qo2->s1->pdir->aang)) 367 { 368 /* A simpel case. The curve cone or the mirror cone 369 are tottally inside the inverted surface cone. */ 370 371 *jstat = 1; 372 goto out; 373 } 374 else if (tang < PI - ANGULAR_TOLERANCE && 375 tang > ANGULAR_TOLERANCE && 376 min(tang,fabs(PI-tang)) < 377 (double)0.8*(PIHALF - qo2->s1->pdir->aang) && 378 qo1->c1->pdir->aang < (double)0.8*(PIHALF-qo2->s1->pdir->aang)) 379 { 380 s1797(qo2->s1,qo1->c1,aepsge,tang,&kstat); 381 if (kstat<0) goto error; 382 else *jstat = kstat; 383 goto out; 384 } 385 } 386 387 388 /* Not a simple case. */ 389 390 out2: *jstat = 0; 391 goto out; 392 393 /* Error. Dimensions conflicting. */ 394 395 err106: *jstat = -106; 396 s6err("s1741",*jstat,kpos); 397 goto out; 398 399 /* Error in lower level routine. */ 400 401 error : *jstat = kstat; 402 s6err("s1741",*jstat,kpos); 403 goto out; 404 405 out: ; 406 } 407