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: s1787.c,v 1.3 2005-02-28 09:04:49 afr Exp $ 45 * 46 */ 47 48 49 #define S1787 50 51 #include "sislP.h" 52 53 #if defined(SISLNEEDPROTOTYPES) 54 void s1787(SISLSurf * ps,double alevel,double aepsge,double epar[],double gpar1[],double gpar2[],int * jstat)55 s1787(SISLSurf *ps,double alevel,double aepsge,double epar[], 56 double gpar1[],double gpar2[],int *jstat) 57 #else 58 void s1787(ps,alevel,aepsge,epar,gpar1,gpar2,jstat) 59 SISLSurf *ps; 60 double alevel; 61 double aepsge; 62 double epar[]; 63 double gpar1[]; 64 double gpar2[]; 65 int *jstat; 66 #endif 67 /* 68 ********************************************************************* 69 * 70 ********************************************************************* 71 * 72 * PURPOSE : In surface-point intersection in dimention one. 73 * To search for endpoints on an intersection curve. 74 * The intersection curve must go through the parameter 75 * pair epar. If the function find a closed curve 76 * it will return a status message 2. 77 * Else if epar is an edge point the function 78 * will return parameter values for the other end point 79 * on the intersection curve in gpar1 and a status message 80 * 1. If epar is an internal point the function will 81 * return the parameter values for the two endpoints on 82 * the intersection curve in gpar1 and gpar2 and a status 83 * message 0. 84 * 85 * 86 * INPUT : ps - The surface in intersection. 87 * alevel - The constant value the surface is intersected with. 88 * aepsge - Geometry resolution. 89 * epar[2] - Parameter values for the intersection point. 90 * 91 * 92 * 93 * OUTPUT : gpar1[2] - Parameter values for the first intersection point. 94 * One of the ends are within computer tolerance from 95 * epar, then this point is put in this variable and 96 * the status 1 is returned. 97 * gpar2[2] - Parameter values for the second intersection point 98 * on the edges. If closed curve found with no singular 99 * point this point is equal to the first point. 100 * jstat - status messages 101 * < 0 : Error. 102 * = 0 : The marching did not succeed. 103 * = 11 : epar == gpar1. Open curve. 104 * = 12 : epar == gpar1. Open curve. gpar2 inside. 105 * = 13 : epar == gpar1. Open curve. gpar1 inside. 106 * = 14 : epar == gpar1. Open curve. Both ends inside. 107 * = 16 : epar == gpar1. Closed curve. gpar2 == gpar1 108 * = 17 : epar == gpar1. Closed curve. gpar2 singular. 109 * = 21 : epar != output points. Open curve. 110 * = 22 : epar != output points. Open curve. gpar2 inside. 111 * = 24 : epar != output points. Open curve. 112 * Output points inside. 113 * = 26 : epar != output points. Closed curve. gpar2 == gpar1 114 * = 27 : epar != output points. Closed curve. 115 * gpar2 singular. 116 * 117 * 118 * METHOD : 119 * 120 * 121 * REFERENCES : 122 * 123 * WRITTEN BY : Tor Dokken, SI, Oslo, Norway. July 1989 124 * 125 ********************************************************************* 126 */ 127 { 128 int kdeg=1; /* Indicate that a plane is used */ 129 int kk1,kk2,kn1,kn2; /* Orders and numbers of vertices */ 130 int kstat; /* Status variable */ 131 int kkm1,kkm2; /* Orders minus 1 */ 132 int kincre; /* Number of doubles in first vertex direction */ 133 int kpos=0; /* Position of error */ 134 int ki,kj,kl,kstop; 135 int kcur,kgraph; /* Indicators telling to control type of output 136 from marching */ 137 int kmark1,kmark2,kclose,kmatch1,kmatch2; /* Flags */ 138 double tmax; /* Box size */ 139 double tstart,tlength;/* Variables used in Marsdens identity */ 140 double tfak; 141 double tdum1; /* Max knot value used in DEQUAL comparing. */ 142 double tdum2; /* Max knot value used in DEQUAL comparing. */ 143 double tsum,*sp,*sq; 144 double simpli[4]; /* Description of plane */ 145 double *st1,*st2,*scoef; /* Knots and vertices of input surface */ 146 double *s3coef=SISL_NULL; /* 3-D coeff */ 147 148 double tepsco = REL_COMP_RES; 149 double tepsge; 150 double sval1[2]; /* Limits of parameter plane in first SISLdir */ 151 double sval2[2]; /* Limits of parameter plane in second SISLdir */ 152 double *spar1,*spar2; /* Pointers to arrays */ 153 double *spar=SISL_NULL; /* Pointer to allocated values for parameter values*/ 154 SISLSurf *qs=SISL_NULL; /* 3-D version of surface */ 155 SISLCurve *qcrv; /* Curve in parameter plane */ 156 SISLIntcurve *qintcr=SISL_NULL;/* Intersection curve object */ 157 kk1 = ps -> ik1; 158 kk2 = ps -> ik2; 159 kn1 = ps -> in1; 160 kn2 = ps -> in2; 161 st1 = ps -> et1; 162 st2 = ps -> et2; 163 scoef = ps -> ecoef; 164 sval1[0] = st1[kk1-1]; 165 sval1[1] = st1[kn1]; 166 sval2[0] = st2[kk2-1]; 167 sval2[1] = st2[kn2]; 168 169 /* Allocate array for 3-D representation of surface */ 170 171 if((s3coef = newarray(kn1*kn2*3,DOUBLE)) == SISL_NULL) goto err101; 172 173 sh1992su(ps,0,aepsge,&kstat); 174 if (kstat < 0) goto error; 175 176 tmax = ps->pbox->e2max[0][0] - ps->pbox->e2min[0][0]; 177 178 /* Make description of plane */ 179 180 simpli[0] = (double)0.0; 181 simpli[1] = (double)0.0; 182 simpli[2] = (double)1.0; 183 simpli[3] = -alevel; 184 185 /* Make 3-D description of the surface */ 186 187 188 /* Make representation of coefficients from Marsdens identity for the 189 * function f(t) = t, with the knot vector in first parameter direction 190 * scaled to [0,tmax]. This will be used as the x-coordinate in the 3-D 191 * representation */ 192 193 tstart = st1[kk1-1]; 194 tlength = st1[kn1] - tstart; 195 tfak = tmax/tlength; 196 kkm1 = kk1 - 1; 197 kincre = 3*kn1; 198 199 for (ki=0,kl=0,sp=s3coef ; ki<kn1 ; ki++,kl+=3,sp+=3) 200 { 201 tsum = (double)0.0; 202 kstop = ki+kk1; 203 for (kj=ki+1;kj<kstop;kj++) 204 tsum +=st1[kj]; 205 206 tsum = (tsum/kkm1-tstart)*tfak; 207 208 209 /* Copy x-coordinate to the other vertex rows */ 210 /*UJK,changed from kj<kn to kj<kn2.*/ 211 for (kj=0,sq=sp ; kj<kn2 ; kj++,sq+=kincre) *sq = tsum; 212 213 } 214 215 /* Make representation of coefficients from Marsdens identity for the 216 * function f(t) = t, with the knot vector in second parameter direction 217 * scaled to [0,tfak]. This will be used as the x-coordinate in the 3-D 218 * representation */ 219 220 kkm2 = kk2 - 1; 221 tstart = st2[kk2-1]; 222 tlength = st2[kn2] - tstart; 223 tfak = tmax/tlength; 224 for (ki=0,sp=s3coef+1 ; ki< kn2 ; ki++) 225 { 226 tsum = (double)0.0; 227 kstop = ki+kk2; 228 for (kj=ki+1;kj<kstop;kj++) 229 tsum +=st2[kj]; 230 231 tsum = (tsum/kkm2-tstart)*tfak; 232 233 /* Copy to remaining y-coordinates in first vertex row */ 234 235 for (kj=0 ; kj<kn1 ; kj++,sp+=3) *sp = tsum; 236 237 } 238 239 /* Copy z-coordinates */ 240 241 for (kj=0,sp=s3coef+2,sq=scoef ; kj < kn2 ; kj++) 242 for (ki=0 ; ki<kn1 ; ki++,sp+=3,sq++) 243 *sp = *sq; 244 245 /* Make 3-D surface */ 246 247 if((qs = newSurf(kn1,kn2,kk1,kk2,st1,st2,s3coef,1,3,1)) == SISL_NULL) goto err101; 248 249 kgraph = 0; 250 kcur = 3; 251 252 /* Make an intersection curve object with the parameter value */ 253 254 if ((spar=newarray(2,DOUBLE))==SISL_NULL) goto err101; 255 memcopy(spar,epar,2,DOUBLE); 256 257 if((qintcr = newIntcurve(1,2,0,spar,SISL_NULL,0)) == SISL_NULL) goto err101; 258 259 kcur = 2; 260 kgraph = 0; 261 tepsge = tmax*(double)0.01; 262 s1313(qs,simpli,kdeg,tepsco,tepsge,tmax,qintcr,kcur,kgraph,&kstat); 263 if (kstat==-185) goto war00; 264 if (kstat<0) goto error; 265 266 /* Identify first and last parameter pair in the intersection curve */ 267 268 qcrv = qintcr -> ppar1; 269 if (qcrv == SISL_NULL) goto war00; 270 271 spar1 = qcrv -> ecoef; 272 spar2 = spar1 + 2*(qcrv->in)-2; 273 /* Check if any of the points lie on the boundary */ 274 275 kmark1 = 0; 276 tdum1 = (double)2.0*max(fabs(sval1[0]),fabs(sval1[1])); 277 tdum2 = (double)2.0*max(fabs(sval2[0]),fabs(sval2[1])); 278 279 if (DEQUAL(spar1[0]+tdum1,sval1[0]+tdum1) || DEQUAL(spar1[0]+tdum1,sval1[1]+tdum1) || 280 DEQUAL(spar1[1]+tdum2,sval2[0]+tdum2) || DEQUAL(spar1[1]+tdum2,sval2[1]+tdum2) ) 281 kmark1 = 1; 282 283 kmark2 = 0; 284 285 if (DEQUAL(spar2[0]+tdum1,sval1[0]+tdum1) || DEQUAL(spar2[0]+tdum1,sval1[1]+tdum1) || 286 DEQUAL(spar2[1]+tdum2,sval2[0]+tdum2) || DEQUAL(spar2[1]+tdum2,sval2[1]+tdum2) ) 287 kmark2 = 1; 288 289 /* Check if closed */ 290 291 kclose = 0; 292 if (spar1[0] == spar2[0] && spar1[1] == spar2[1]) kclose = 1; 293 294 /* Check if first points matches start point */ 295 296 kmatch1 = 0; 297 if (DEQUAL(epar[0]+tdum1,spar1[0]+tdum1) && DEQUAL(epar[1]+tdum2,spar1[1]+tdum2) ) 298 kmatch1 = 1; 299 300 /* Check if second points matches start point */ 301 302 kmatch2 = 0; 303 if (DEQUAL(epar[0]+tdum1,spar2[0]+tdum1) && DEQUAL(epar[1]+tdum2,spar2[1]+tdum2) ) 304 kmatch2 = 1; 305 306 /* Check if any point matches start point */ 307 308 if (kmatch1 == 1 || kmatch2 == 1) 309 { 310 /* Start point matches one of the end points, status values in 311 the range 11-19*/ 312 313 if (kmark1 == 1 && kmark2 == 1 && kclose == 0) 314 { 315 /* Open curve, status 11 */ 316 *jstat = 11; 317 if(kmatch1==1) 318 goto copy; 319 else 320 goto invcopy; 321 } 322 else if (kmark1 ==1 || (kmark2 == 1 && kclose == 0)) 323 { 324 /* Open curve one point inside status 12 or 13 */ 325 326 if (kmark1 == 1 && kmatch1 == 1) 327 { 328 *jstat = 12; 329 goto copy; 330 } 331 else if (kmark2 == 1 && kmatch2 == 1) 332 { 333 *jstat = 12; 334 goto invcopy; 335 } 336 if (kmark1 == 1 && kmatch2 == 1) 337 { 338 *jstat = 13; 339 goto invcopy; 340 } 341 if (kmark2 == 1 && kmatch1 == 1) 342 { 343 *jstat = 13; 344 goto copy; 345 } 346 } 347 else if (kclose == 0) 348 { 349 /* Both ends inside */ 350 *jstat = 14; 351 if(kmatch1==1) 352 goto copy; 353 else 354 goto invcopy; 355 } 356 else if(kmatch1 == 1) 357 { 358 /* Closed curve, no singularity */ 359 *jstat = 16; 360 memcopy(gpar1,spar1,2,DOUBLE); 361 memcopy(gpar2,spar1,2,DOUBLE); 362 goto out; 363 } 364 else 365 { 366 /* Closed curve, with singularity */ 367 *jstat=17; 368 memcopy(gpar1,epar ,2,DOUBLE); 369 memcopy(gpar2,spar1,2,DOUBLE); 370 goto out; 371 } 372 } 373 else 374 { 375 /* epar does not match produced end points, status messages in 376 21-29 the range */ 377 378 if (kmark1 ==1 && kmark2 ==1 && kclose == 0) 379 { 380 /* Open curve, status 11 */ 381 *jstat = 21; 382 memcopy(gpar1,spar1,2,DOUBLE); 383 memcopy(gpar2,spar2,2,DOUBLE); 384 goto out; 385 } 386 else if (kmark1 ==1 && kclose == 0) 387 { 388 /* Open curve one point inside status 12 */ 389 *jstat=22; 390 goto copy; 391 } 392 else if (kmark2 ==1 && kclose == 0) 393 { 394 /* Open curve one point inside status 12 */ 395 *jstat=22; 396 goto invcopy; 397 } 398 else if (kclose == 0) 399 { 400 /* Both ends inside */ 401 *jstat=24; 402 goto copy; 403 } 404 else if(kmatch1 == 1) 405 { 406 /* Closed curve, no singularity */ 407 *jstat=26; 408 memcopy(gpar1,spar1,2,DOUBLE); 409 memcopy(gpar2,spar1,2,DOUBLE); 410 } 411 else 412 { 413 /* Closed curve, with singularity */ 414 *jstat = 27; 415 memcopy(gpar1,epar ,2,DOUBLE); 416 memcopy(gpar2,spar1,2,DOUBLE); 417 goto out; 418 } 419 } 420 /* Marching produced no curve */ 421 422 war00: *jstat = 0; 423 memcopy(gpar1,epar,2,DOUBLE); 424 memcopy(gpar2,epar,2,DOUBLE); 425 goto out; 426 427 copy: 428 memcopy(gpar1,spar1,2,DOUBLE); 429 memcopy(gpar2,spar2,2,DOUBLE); 430 goto out; 431 432 invcopy: 433 memcopy(gpar1,spar2,2,DOUBLE); 434 memcopy(gpar2,spar1,2,DOUBLE); 435 goto out; 436 437 /* Error in space allocation */ 438 err101: 439 *jstat = -101; 440 s6err("s1787",*jstat,kpos); 441 goto out; 442 443 /* Error in lower level function */ 444 error: 445 *jstat = kstat; 446 s6err("s1787",*jstat,kpos); 447 goto out; 448 449 out: 450 if (s3coef != SISL_NULL) freearray(s3coef); 451 if (qs != SISL_NULL) freeSurf (qs); 452 if (qintcr != SISL_NULL) freeIntcurve(qintcr); 453 } 454