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: s1252.c,v 1.3 2001-03-19 15:58:43 afr Exp $ 45 * 46 */ 47 48 49 #define S1252 50 51 #include "sislP.h" 52 53 /* 54 * Forward declarations. 55 * --------------------- 56 */ 57 #if defined (SISLNEEDPROTOTYPES) 58 static void s1252_s6corr(double *,double,double [],int,int,int *,int *); 59 static void s1252_s6dir(double *,double,double [],double,double); 60 #else 61 static void s1252_s6corr(); 62 static void s1252_s6dir(); 63 #endif 64 65 #if defined (SISLNEEDPROTOTYPES) 66 void 67 s1252(SISLCurve * pcurve,double aepsge,double astart,double * cpos,int * jstat)68 s1252(SISLCurve *pcurve,double aepsge,double astart,double *cpos,int *jstat) 69 #else 70 void s1252(pcurve,aepsge,astart,cpos,jstat) 71 SISLCurve *pcurve; 72 double aepsge; 73 double astart; 74 double *cpos; 75 int *jstat; 76 #endif 77 /* 78 ********************************************************************* 79 * 80 ********************************************************************* 81 * 82 * PURPOSE : Newton iteration to a local maximum on a function 83 * in one variable. 84 * 85 * INPUT : pcurve - Pointer to the first curve in the intersection. 86 * aepsge - Geometry resolution. 87 * astart - Start value of the first curve to the iteration. 88 * 89 * 90 * 91 * OUTPUT : cpos - Parameter value of of first curve in intersection 92 * point. 93 * jstat - status messages 94 * = 2 : Divergence or approximative 95 * intersection found. 96 * = 1 : Intersection found. 97 * < 0 : error. 98 * 99 * 100 * METHOD : Newton iteration in one parameter direction. 101 * 102 * 103 * REFERENCES : 104 * 105 *- 106 * CALLS : s1221 - Evaluate expression of curve in given 107 * parameter values. 108 * 109 * WRITTEN BY : Tor Dokken, SI, Mars 1989 110 * 111 ********************************************************************* 112 */ 113 { 114 int kstat = 0; /* Local status variable. */ 115 int kpos = 0; /* Position of error. */ 116 int kleft=0; /* Variables used in the evaluator. */ 117 int kder=3; /* Order of derivatives to be calulated */ 118 int kdim; /* Dimension of space the curves lie in */ 119 int knbit; /* Number of iterations */ 120 int kn,kk; /* Number of vertices and order */ 121 int kdir=1; /* Direction of derivative to be calculated */ 122 double tstart,tend; /* Ends of parameter interval of first curve. */ 123 double tdelta; /* Parameter interval of the curves. */ 124 double tdist=DZERO; /* Distance between position and origo. */ 125 double td; /* Distances between old and new parameter value */ 126 double tnext; /* Parameter-value of expression in first curve. */ 127 double tprev; /* Previous difference between the curves. */ 128 double sval[4]; /* Value ,first and second derivative on function */ 129 double *st; /* Knot vector */ 130 double ref; /* Refferance value for equality test. */ 131 132 /* Test input. */ 133 134 if (pcurve->idim != 1) goto err106; 135 136 kdim = pcurve -> idim; 137 138 /* Fetch endpoints and the intervals of parameter interval of curves. */ 139 140 st = pcurve->et; 141 kn = pcurve->in; 142 kk = pcurve->ik; 143 144 tstart = *(pcurve->et + pcurve->ik - 1); 145 tend = *(pcurve->et + pcurve->in); 146 tdelta = tend - tstart; 147 if (tdelta == DZERO) tdelta = fabs(tend); 148 if (tdelta == DZERO) tdelta = (double)1.0; 149 150 /* Initiate variables. */ 151 152 tnext = astart; 153 154 /* Evaluate 0-1.st derivatives of function */ 155 156 s1221(pcurve,kder,tnext,&kleft,sval,&kstat); 157 if (kstat < 0) goto error; 158 159 tprev = sval[0]; 160 161 /* Evaluate step */ 162 163 s1252_s6dir(&td,tnext,sval,tstart,tend); 164 165 /* Correct if we not are inside the parameter intervall. */ 166 167 s1252_s6corr(&td,tnext,st,kn,kk,&kleft,&kdir); 168 169 /* Iterate to find the intersection point. */ 170 171 for (knbit = 0; knbit < 20; knbit++) 172 { 173 174 /* If the tnext is a break point test if it is a local maximum */ 175 176 if (kdir == -2 || kdir == 2) 177 { 178 double tder1,tder2; 179 /* Break point, test if local maximum */ 180 181 s1221(pcurve,kder,tnext,&kleft,sval,&kstat); 182 if (kstat < 0) goto error; 183 tder2 = sval[1]; 184 185 s1227(pcurve,kder,tnext,&kleft,sval,&kstat); 186 if (kstat < 0) goto error; 187 tder1 = sval[1]; 188 189 /* Test if top point */ 190 191 if (tder1>=DZERO && tder2<=DZERO) break; 192 193 /* Not a top point */ 194 } 195 196 197 /* Evaluate 0-1.st derivatives of both curves, dependent of the 198 sign of td we calculate derivatives from the right or the left */ 199 200 if (kdir>=1) 201 { 202 s1221(pcurve,kder,tnext+td,&kleft,sval,&kstat); 203 if (kstat < 0) goto error; 204 } 205 else 206 { 207 s1227(pcurve,kder,tnext+td,&kleft,sval,&kstat); 208 if (kstat < 0) goto error; 209 } 210 211 tdist = sval[0]; 212 if (fabs(tdist) < (double)1.0) ref = (double)2.0; 213 else ref = DZERO; 214 215 if (tdist >= tprev || DEQUAL(ref+tdist,ref+tprev)) 216 { 217 tnext += td; 218 219 /* Evaluate step */ 220 s1252_s6dir(&td,tnext,sval,tstart,tend); 221 s1252_s6corr(&td,tnext,st,kn,kk,&kleft,&kdir); 222 223 if (fabs(td/tdelta) <= REL_COMP_RES) break; 224 225 tprev = tdist; 226 227 } 228 229 /* Not converging, correct and try again. */ 230 231 else 232 { 233 234 td /= (double)2; 235 if (fabs(td/tdelta) <= REL_COMP_RES) break; 236 } 237 238 239 } 240 241 242 /* Iteration stopped, test if point founds found is within resolution */ 243 244 if (tdist <= aepsge) 245 *jstat = 1; 246 else 247 *jstat = 2; 248 249 /*ujk,july 89:*/ 250 /* Test if the iteration is close to a knot */ 251 if (DEQUAL(tnext,pcurve->et[kleft])) 252 *cpos = pcurve->et[kleft]; 253 else if (DEQUAL(tnext,pcurve->et[kleft+1])) 254 *cpos = pcurve->et[kleft+1]; 255 else 256 *cpos = tnext; 257 258 /* Iteration completed. */ 259 260 goto out; 261 262 263 /* Error in input. Conflicting dimensions. */ 264 265 err106: *jstat = -106; 266 s6err("S1252",*jstat,kpos); 267 goto out; 268 269 /* Error in lower level routine. */ 270 271 error : *jstat = kstat; 272 s6err("S1252",*jstat,kpos); 273 goto out; 274 275 out:; 276 } 277 278 #if defined (SISLNEEDPROTOTYPES) 279 static void s1252_s6corr(double * gdn,double acoef,double et[],int in,int ik,int * ileft,int * jdir)280 s1252_s6corr(double *gdn,double acoef,double et[], 281 int in,int ik,int *ileft,int *jdir) 282 #else 283 static void s1252_s6corr(gdn,acoef,et,in,ik,ileft,jdir) 284 double *gdn; 285 double acoef; 286 double et[]; 287 int in; 288 int ik; 289 int *ileft; 290 int *jdir; 291 #endif 292 /* 293 ********************************************************************* 294 * 295 ********************************************************************* 296 * 297 * PURPOSE : Adjust the step to not cross knot values or out 298 * of the curve 299 * 300 * 301 * INPUT : acoef - Current parameter value 302 * st - knots 303 * in - number of vertices 304 * ik - polynomial order 305 * 306 * INPUT/OUTPUT : 307 * ileft - Pointer to the interval in the knot vector 308 * where ax is located, check the relations above. 309 * 310 * 311 * OUTPUT : gdn - Old and new step value. 312 * jdir - Direction of derivative to be calculated 313 * -2 Negative direction acoef at break point 314 * -1 Negative direction 315 * +1 Positive direction 316 * +2 Positive direction acoef at break point 317 * 318 * METHOD : We are making the step keep inside the parameter interval. 319 * 320 * REFERENCES : 321 * 322 *- 323 * CALLS : 324 * 325 * 326 * WRITTEN BY : Tor Dokken, SI, Mars 1989 327 * 328 ********************************************************************* 329 */ 330 { 331 int kmult,kstat; 332 333 /* Make sure the point is inside the interval */ 334 335 *gdn = MAX(et[ik-1]-acoef,*gdn); 336 *gdn = MIN(et[in] -acoef,*gdn); 337 338 if (acoef+*gdn<et[*ileft] && acoef>et[*ileft]) 339 { 340 *gdn = MAX(et[*ileft]-acoef,*gdn); 341 } 342 343 else if(acoef<et[*ileft+1] && acoef+*gdn>et[*ileft+1]) 344 { 345 /* We cross a knot value */ 346 347 *gdn = MIN(et[*ileft+1]-acoef,*gdn); 348 } 349 350 /* Make sure that we calculate the left or right handed derivatives */ 351 352 if (*gdn>=0) 353 { 354 *jdir = 1; 355 } 356 else 357 { 358 *jdir = -1; 359 } 360 361 kmult = s6knotmult(et,ik,in,ileft,acoef,&kstat); 362 363 if (acoef==et[*ileft]) 364 { 365 366 if(kmult>ik-2) 367 { 368 if (*jdir == -1) 369 { 370 *jdir = -2; 371 } 372 else 373 { 374 *jdir = 2; 375 } 376 } 377 } 378 } 379 380 #if defined (SISLNEEDPROTOTYPES) 381 static void s1252_s6dir(double * cdiff,double acoef,double eval[],double astart,double aend)382 s1252_s6dir(double *cdiff,double acoef,double eval[],double astart, 383 double aend) 384 #else 385 static void s1252_s6dir(cdiff,acoef,eval,astart,aend) 386 double *cdiff; 387 double acoef; 388 double eval[]; 389 double astart; 390 double aend; 391 #endif 392 /* 393 ********************************************************************* 394 * 395 ********************************************************************* 396 * 397 * PURPOSE : To compute the next iteration step 398 * 399 * INPUT : eval - Value and derivative 400 * astart - The lower interval end 401 * aend - The higher interval end 402 * 403 * 404 * OUTPUT : cdiff - Iteration step 405 *- 406 * 407 * WRITTEN BY : Tor Dokken, SI, Mars 1989 408 * 409 ********************************************************************* 410 */ 411 { 412 double t1,t2,t3,t4,t5,t6; /* Constants in equation. */ 413 double tmax; /* Max values in equation. */ 414 double ttol=(double)1e-10; /* Relative tolerance in equation. */ 415 416 /* Dummy statements to avoid warning. */ 417 t1=acoef; 418 t2=astart; 419 t3=aend; 420 421 422 t1 = eval[1]; 423 t2 = eval[2]; 424 t3 = eval[3]/(double)2.0; 425 426 tmax = max(fabs(t1),fabs(t2)); 427 tmax = max(fabs(t3),tmax); 428 429 if (DEQUAL(tmax,DZERO)) *cdiff = DZERO; 430 else if (fabs(t3)/tmax < ttol) /* The second degree part is degenerated. */ 431 { 432 if (fabs(t2) == DZERO ) *cdiff = DZERO; 433 else *cdiff = (-t1/t2); 434 } 435 else 436 { 437 /* An ordinary second degree equation. */ 438 t4 = t2*t2 - (double)4*t3*t1; 439 if (t4 < DZERO) 440 { 441 /* Use linear equation. */ 442 if (fabs(t2) == DZERO ) *cdiff = DZERO; 443 else *cdiff = (-t1/t2); 444 } 445 446 else 447 { 448 t6 = sqrt(t4); 449 t5 = (-t2 + t6)/((double)2*t3); 450 t6 = (-t2 - t6)/((double)2*t3); 451 t4 = min(fabs(t5),fabs(t6)); 452 453 /* We have two solutions and we want to use the one 454 with the one with smallest value. */ 455 456 if (t4 == DZERO) 457 { 458 /* Use linear equation. */ 459 if (fabs(t2) == DZERO ) *cdiff = DZERO; 460 else *cdiff = (-t1/t2); 461 } 462 else if (fabs(t5) <= fabs(t6)) *cdiff = t5; 463 else *cdiff = t6; 464 } 465 } 466 } 467