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: sh1461.c,v 1.2 2001-03-19 15:59:03 afr Exp $ 45 * 46 */ 47 48 49 #define SH1461 50 51 #include "sislP.h" 52 53 #if defined(SISLNEEDPROTOTYPES) 54 typedef void (*fshapeProc)(double [],double [],int,int,int *); 55 typedef void (*initProc)(fshapeProc,SISLCurve *[],int,double [], 56 double [],double [],int *); 57 #else 58 typedef void (*fshapeProc)(); 59 typedef void (*initProc)(); 60 #endif 61 62 #if defined(SISLNEEDPROTOTYPES) 63 static void sh1461_s9coef(double [],double [],double [],double [],int, 64 double *,double *,double *,double *,int *); 65 static void sh1461_s9hermit(double [],int,int,int *); 66 static void sh1461_s9chcoor(double [],int,double [],int,double [], 67 int,double [],double [],double [],int, 68 double [],int *,int *); 69 static void sh1461_s9mult(double [],double [],int,int,double [],int *); 70 static void sh1461_s9comder(int,int,double [],int,double,double,double, 71 double,double [],int *); 72 static double sh1461_s9ang(double [],double [],int); 73 #else 74 static void sh1461_s9coef(); 75 static void sh1461_s9hermit(); 76 static void sh1461_s9chcoor(); 77 static void sh1461_s9mult(); 78 static void sh1461_s9comder(); 79 static double sh1461_s9ang(); 80 #endif 81 82 #if defined(SISLNEEDPROTOTYPES) 83 void sh1461(fshapeProc fshape,initProc f_initmid,SISLCurve * vboundc[],int icurv,SISLSurf * vsurf[],int * jstat)84 sh1461(fshapeProc fshape,initProc f_initmid, 85 SISLCurve *vboundc[],int icurv,SISLSurf *vsurf[],int *jstat) 86 #else 87 void sh1461(fshape,f_initmid,vboundc,icurv,vsurf,jstat) 88 fshapeProc fshape; 89 initProc f_initmid; 90 SISLCurve *vboundc[]; 91 SISLSurf *vsurf[]; 92 int icurv; 93 int *jstat; 94 #endif 95 /* 96 ********************************************************************* 97 * 98 * PURPOSE : Make a blend consisting of n tensor-product surfaces over 99 * a vertex region with n edges, n>=3. Along each edge position 100 * and cross-tangent conditions are given. The blend is supposed 101 * to be G1. 102 * 103 * 104 * 105 * INPUT : fshape - Application driven routine that gives the user an 106 * ability to change the middle point of the region 107 * (the vertex at which the blending surfaces meet), 108 * and the tangent vectors in the middle point along 109 * the curves which divedes the region. 110 * f_initmid - Function used to compute the position and 111 * derivatives of the first blending surface in 112 * the middle of the vertex region. This function 113 * is dependant on the number of edges of the region. 114 * vboundc - Pointers to curves describing edge-conditions. 115 * For each edge, two curves are given, one describing 116 * position and one cross-derivatives. The edges are 117 * sorted counter-clockwise around the region, and the 118 * curves are oriented counter-clockwise. The array has 119 * dimension 2*icurv, i.e. 6. 120 * icurv - Number of boundary curves. icurv >= 3. 121 * 122 * 123 * OUTPUT : vsurf - Pointers to blending surfaces. The array has dimension 124 * icurv. 125 * jstat - status messages 126 * > 0 : warning 127 * = 0 : ok 128 * < 0 : error 129 * 130 * 131 * METHOD : Hahn's method. Split the region in n 4-sided regions. Define 132 * position and cross tangent conditions along the inner edges 133 * such that G1-continuity is satisfied. Define the blending 134 * patches over the 4-sided regions as Coons patches. 135 * 136 * 137 * REFERENCES : Joerg Hahn : "Filling Polygonal Holes with Rectangular Patches" 138 * Theory and Practice of Geometric Modeling, Blackburn Oct 1988. 139 * 140 * 141 * USE : 3D geometry only. 142 * 143 *- 144 * CALLS : s1221 - Evaluate curve. 145 * s1706 - Turn orientation of curve. 146 * s1401 - Rectangular blending. 147 * newCurve - Create new curve. 148 * freeCurve - Free space occupied by a curve. 149 * sh1461_s9coef - Compute coefficients used for computing 150 * derivatives in the midpoint. 151 * sh1461_s9hermit - Compute vertices of Bezier curve. 152 * sh1461_s9chcoor - Perform coordinate change on derivative curve. 153 * sh1461_s9mult - Multiply two 4-order Bezier curves. 154 * sh1461_s9comder - Compute given 2. derivative of patch in midpoint 155 * of region. 156 * sh1461_s9ang - The angle between two vectors. 157 * 158 * 159 * WRITTEN BY : Vibeke Skytt, SI, 03.90. 160 * 161 ********************************************************************* 162 */ 163 { 164 int kstat = 0; /* Status variable. */ 165 int kdim = 3; /* Dimension of geometry space. */ 166 int ki,kj,kk,kh; /* Counters. */ 167 int kder = 2; /* Number of derivatives to 168 compute while evaluating curve. */ 169 int knmb = 0; /* Number of derivatives computed 170 in midpoint of patch. */ 171 int kleft = 0; /* Parameter of curve evaluation. */ 172 int kkcrt2; /* Order of cross tangent curve. */ 173 int lder[4]; 174 double tpar; /* Parameter value at which to 175 evaluate curve. */ 176 double tro00,tro01,tro10,tro11; /* Coefficients used to compute 177 derivatives at midpoint of region. */ 178 double *sder = SISL_NULL; /* Array of derivatives at midpoint */ 179 double *stwist = SISL_NULL; /* Twist vectors of corners of region. */ 180 double *stang = SISL_NULL; /* Tangent vectors at midpoint. */ 181 double spos[15]; /* Conditions and vertices of position 182 curve along inner edge of region. */ 183 double scrt1[21]; /* Conditions and vertices of cross tangent 184 curve along inner edge of region. */ 185 double scrt2[12]; /* Conditions and vertices of cross tangent 186 curve along inner edge of region. */ 187 double stpos[10]; /* Knot vector corresponding to spos. */ 188 double stcrt1[14]; /* Knot vector corresponding to scrt1. */ 189 double stcrt2[8]; /* Knot vector corresponding to scrt2. */ 190 double sblend[4]; /* Vertices of blending function. */ 191 double sdum[6]; /* Value and 1. derivative of position 192 boundary curve. */ 193 double sdum2[6]; /* Value and 1. derivative of cross 194 tangent boundary curve. */ 195 double *stwist2 = SISL_NULL; /* Twist vectors of 4-sided region. */ 196 double *st,*st2; /* Pointers into knot vectors. Used to 197 change parametrization of curve. */ 198 double tstart; /* Start parameter of knot vector. */ 199 SISLCurve *qc; /* Pointer to curve. */ 200 SISLCurve **qbound = SISL_NULL; /* Boundary conditions of 4-sided regions. */ 201 202 for (ki=0; ki<4; ki++) lder[ki] = 2; 203 204 for (ki=0; ki<=kder; ki++) knmb += ki + 1; 205 206 /* Initiate knot vectors of position curve and cross tangent 207 curves of inner edges. */ 208 209 for (ki=0; ki<5; ki++) 210 { 211 stpos[ki] = (double)0.0; 212 stpos[5+ki] = (double)1.0; 213 } 214 for (ki=0; ki<4; ki++) 215 { 216 stcrt2[ki] = (double)0.0; 217 stcrt2[4+ki] = (double)1.0; 218 } 219 for (ki=0; ki<7; ki++) 220 { 221 stcrt1[ki] = (double)0.0; 222 stcrt1[7+ki] = (double)1.0; 223 } 224 225 226 /* Allocate scratch for local arrays. */ 227 228 if ((sder = new0array(icurv*kdim*knmb,DOUBLE)) == SISL_NULL) goto err101; 229 if ((stwist = newarray(icurv*kdim,DOUBLE)) == SISL_NULL) goto err101; 230 if ((stwist2 = newarray(4*icurv*kdim,DOUBLE)) == SISL_NULL) goto err101; 231 if ((stang = newarray(icurv*kdim,DOUBLE)) == SISL_NULL) goto err101; 232 if ((qbound = newarray(8*icurv,SISLCurve*)) == SISL_NULL) goto err101; 233 234 for (ki=0; ki<icurv; ki++) 235 { 236 /* Test dimension of curve. */ 237 238 if (vboundc[ki]->idim != kdim) goto err102; 239 240 /* Initiate twist vector of n-sided region. */ 241 242 tpar = (double)0.0; 243 s1221(vboundc[2*ki+1],1,tpar,&kleft,sdum,&kstat); 244 if (kstat < 0) goto error; 245 246 memcopy(stwist+ki*kdim,sdum+kdim,kdim,double); 247 } 248 249 /* Initiate those twist vectors of the rectangular blending surfaces 250 that coincide with the twist vectors of the n-sided region. */ 251 252 for (ki=0; ki<kdim; ki++) 253 for (kj=0; kj<icurv; kj++) 254 { 255 kk = (kj + 1) % icurv; 256 stwist2[(kj*4+2)*kdim+ki] = stwist[kk*kdim+ki]/(double)4.0; 257 } 258 259 260 /* Set up blending function. */ 261 262 sblend[0] = (double)1.0; 263 sblend[1] = sblend[2] = sblend[3] = (double)0.0; 264 sh1461_s9hermit(sblend,4,1,&kstat); 265 266 /* Set up value and derivative of the 1. blending surface in the 267 midpoint of the vertex region. */ 268 269 f_initmid(fshape,vboundc,icurv,stwist,stang,sder,&kstat); 270 if (kstat < 0) goto error; 271 272 for (ki=1; ki<icurv; ki++) 273 { 274 /* Compute value and derivatives of the other blending surfaces 275 in the midpoint of the vertex region. */ 276 277 /* Copy the value in the midpoint. */ 278 279 memcopy(sder+ki*knmb*kdim,sder,kdim,DOUBLE); 280 281 /* Copy the first derivatives of the current patch. */ 282 283 memcopy(sder+(ki*knmb+1)*kdim,stang+(ki-1)*kdim,kdim,DOUBLE); 284 memcopy(sder+(ki*knmb+2)*kdim,stang+ki*kdim,kdim,DOUBLE); 285 286 /* The second derivative in the first parameter direction 287 is equal to the first derivative in the second parameter 288 direction of the previous surface. */ 289 290 memcopy(sder+(ki*knmb+3)*kdim,sder+((ki-1)*knmb+5)*kdim,kdim,DOUBLE); 291 292 /* Compute help variables used to define the rest of the derivatives. */ 293 294 sh1461_s9coef(stang+(icurv-1)*kdim,stang,stang+(ki-1)*kdim, 295 stang+ki*kdim,kdim,&tro00,&tro01,&tro10,&tro11,&kstat); 296 if (kstat < 0) goto error; 297 298 /* Compute mixed derivative of the patch. */ 299 300 sh1461_s9comder(1,1,sder+3*kdim,kdim,tro00,tro01,tro10,tro11, 301 sder+(ki*knmb+4)*kdim,&kstat); 302 303 /* Compute the second derivative in the second parameter direction. */ 304 305 sh1461_s9comder(0,2,sder+3*kdim,kdim,tro00,tro01,tro10,tro11, 306 sder+(ki*knmb+5)*kdim,&kstat); 307 } 308 309 /* Set up conditions for 4-edged blending. */ 310 311 for (ki=0; ki<icurv; ki++) 312 { 313 /* Copy value and derivatives of the current patch of the midpoint of 314 the vertex region into arrays containing interpolation conditions 315 and twist vector of current blending surface. */ 316 317 memcopy(spos,sder+ki*knmb*kdim,kdim,DOUBLE); 318 memcopy(spos+kdim,sder+(ki*knmb+2)*kdim,kdim,DOUBLE); 319 memcopy(spos+2*kdim,sder+(ki*knmb+5)*kdim,kdim,DOUBLE); 320 memcopy(scrt2,sder+(ki*knmb+1)*kdim,kdim,DOUBLE); 321 memcopy(scrt2+kdim,sder+(ki*knmb+4)*kdim,kdim,DOUBLE); 322 memcopy(stwist2+4*ki*kdim,sder+(ki*knmb+4)*kdim,kdim,DOUBLE); 323 324 /* Compute value and derivatives of the current inner edge curve at 325 the boundary of the vertex region. */ 326 327 kj = (ki < icurv-1) ? ki+1 : 0; 328 329 qc = vboundc[2*kj]; 330 tpar = (*(qc->et + qc->ik - 1) + *(qc->et + qc->in))/(double)2.0; 331 332 s1221(qc,1,tpar,&kleft,sdum,&kstat); 333 if (kstat < 0) goto error; 334 335 s1221(vboundc[2*kj+1],1,tpar,&kleft,sdum2,&kstat); 336 if (kstat < 0) goto error; 337 338 /* Copy value and derivatives at the current edge into arrays containing 339 interpolation conditions and twist vector of current blending surface. */ 340 341 memcopy(spos+4*kdim,sdum,kdim,DOUBLE); 342 for (kh=0; kh<kdim; kh++) 343 { 344 spos[3*kdim+kh] = sdum2[kh]/(double)2.0; 345 scrt2[3*kdim+kh] = -sdum[kdim+kh]/(double)2.0; 346 stwist2[(ki*4+3)*kdim+kh] = scrt2[2*kdim+kh] = -sdum2[kdim+kh]/(double)4.0; 347 stwist2[(kj*4+1)*kdim+kh] = -stwist2[(ki*4+3)*kdim+kh]; 348 } 349 350 /* Perform Hermit interpolation of position and first cross tangent curve. */ 351 352 sh1461_s9hermit(spos,5,kdim,&kstat); 353 if (kstat < 0) goto error; 354 355 sh1461_s9hermit(scrt2,4,kdim,&kstat); 356 if (kstat < 0) goto error; 357 358 /* Construct second cross tangent curve. */ 359 360 sh1461_s9chcoor(sblend,4,spos,5,scrt2,4,sder+(knmb*ki+1)*kdim,sder+(knmb*ki+2)*kdim, 361 sder+(kj*knmb+2)*kdim,kdim,scrt1,&kkcrt2,&kstat); 362 if (kstat < 0) goto error; 363 364 /* Represent inner boundary conditions as curves. */ 365 366 if ((qbound[8*kj] = newCurve(5,5,stpos,spos,1,kdim,1)) == SISL_NULL) 367 goto err101; 368 if ((qbound[8*kj+1] = newCurve(kkcrt2,kkcrt2,stcrt1,scrt1,1,kdim,1)) == SISL_NULL) 369 goto err101; 370 if ((qbound[8*ki+6] = newCurve(5,5,stpos,spos,1,kdim,1)) == SISL_NULL) 371 goto err101; 372 if ((qbound[8*ki+7] = newCurve(4,4,stcrt2,scrt2,1,kdim,1)) == SISL_NULL) 373 goto err101; 374 375 /* Split current boundary of the vertex region. */ 376 377 s1710(qc,tpar,qbound+8*ki+4,qbound+8*kj+2,&kstat); 378 if (kstat < 0) goto error; 379 380 s1710(vboundc[2*kj+1],tpar,qbound+8*ki+5,qbound+8*kj+3,&kstat); 381 if (kstat < 0) goto error; 382 383 /* Turn direction of curves corresponding to standard edge 3. */ 384 385 s1706(qbound[8*ki+4]); 386 s1706(qbound[8*ki+5]); 387 388 /* Adjust length of derivative curves. */ 389 390 for (kk=0; kk<kdim*(qbound[8*ki+5]->in); kk++) 391 *(qbound[8*ki+5]->ecoef+kk) *= (double)0.5; 392 393 for (kk=0; kk<kdim*(qbound[8*kj+3]->in); kk++) 394 *(qbound[8*kj+3]->ecoef+kk) *= (double)0.5; 395 396 /* Reparametrize boundary curves in order to get correct tangent length. */ 397 398 for (kh=2; kh<6; kh++) 399 { 400 kk = (kh > 3) ? 8*ki : 8*kj; 401 kk += kh; 402 403 for (st=qbound[kk]->et,tstart=st[0], 404 st2=st+qbound[kk]->in+qbound[kk]->ik; 405 st<st2; st++) 406 *st = (double)2.0*(*st - tstart); 407 } 408 409 } 410 411 for (ki=0; ki<icurv; ki++) 412 { 413 /* Perform rectangular blending. */ 414 415 s1401(qbound+8*ki,stwist2+4*ki*kdim,vsurf+ki,&kstat); 416 /* s1390(qbound+8*ki,vsurf+ki,lder,&kstat); */ 417 if (kstat < 0) goto error; 418 } 419 420 /* Blending performed. */ 421 422 *jstat = 0; 423 goto out; 424 425 /* Error in scratch allocation. */ 426 427 err101 : 428 *jstat = -101; 429 goto out; 430 431 /* Dimension not equal to 3. */ 432 433 err102 : 434 *jstat = -102; 435 goto out; 436 437 /* Error in lower level function. */ 438 439 error : 440 *jstat = kstat; 441 goto out; 442 443 out : 444 445 /* Free space occupied by local arrays. */ 446 447 if (sder) freearray(sder); 448 if (stwist) freearray(stwist); 449 if (stwist2) freearray(stwist2); 450 if (stang) freearray(stang); 451 if (qbound) 452 { 453 for (ki=0; ki<8*icurv; ki++) 454 if (qbound[ki]) freeCurve(qbound[ki]); 455 freearray(qbound); 456 } 457 458 return; 459 } 460 461 #if defined(SISLNEEDPROTOTYPES) 462 static void sh1461_s9hermit(double econd[],int icond,int idim,int * jstat)463 sh1461_s9hermit(double econd[],int icond,int idim,int *jstat) 464 #else 465 static void sh1461_s9hermit(econd,icond,idim,jstat) 466 int icond,idim,*jstat; 467 double econd[]; 468 #endif 469 /* 470 ********************************************************************* 471 * 472 * PURPOSE : Hermite interpolation of position and icond-3 derivatives 473 * in one endpoint and position and derivative in the other 474 * endpoint, represented as a Bezier curve on the interval [0,1]. 475 * 476 * 477 * 478 * INPUT : icond - Number of interpolation conditions. 479 * icond = 4 or icond = 5. 480 * idim - Dimension of geometry space. 481 * 482 * 483 * INPUT/OUTPUT : econd - Interpolation conditions as input, Bezier coefficients 484 * as output. The dimension is icond*idim. 485 * 486 * 487 * OUTPUT : jstat - status messages 488 * > 0 : warning 489 * = 0 : ok 490 * < 0 : error 491 * 492 * 493 ********************************************************************* 494 */ 495 { 496 int ki; /* Index. */ 497 498 /* Test input. The number of conditions has to be 4 or 5. */ 499 500 if (icond != 4 && icond != 5) goto err001; 501 502 if (icond == 4) 503 { 504 /* Hermit interpolation with Bezier curve of order 4. */ 505 506 for (ki=0; ki<idim; ki++) 507 { 508 econd[idim+ki] = ONE_THIRD*econd[idim+ki] + econd[ki]; 509 econd[2*idim+ki] = ONE_THIRD*econd[2*idim+ki] + econd[3*idim+ki]; 510 } 511 } 512 513 if (icond == 5) 514 { 515 /* Hermit interpolation with Bezier curve of order 5. */ 516 517 for (ki=0; ki<idim; ki++) 518 { 519 econd[idim+ki] = ONE_FOURTH*econd[idim+ki] + econd[ki]; 520 econd[2*idim+ki] = econd[2*idim+ki]/(double)12.0 521 + (double)2.0*econd[idim+ki] - econd[ki]; 522 econd[3*idim+ki] = ONE_FOURTH*econd[3*idim+ki] + econd[4*idim+ki]; 523 } 524 } 525 526 *jstat = 0; 527 goto out; 528 529 /* Error in input. Number of coefficients not 4 or 5. */ 530 531 err001 : 532 *jstat = -1; 533 goto out; 534 535 out : 536 return; 537 } 538 539 #if defined(SISLNEEDPROTOTYPES) 540 static void sh1461_s9coef(double evec1[],double evec2[],double evec3[],double evec4[],int idim,double * cro00,double * cro01,double * cro10,double * cro11,int * jstat)541 sh1461_s9coef(double evec1[],double evec2[],double evec3[], 542 double evec4[],int idim,double *cro00,double *cro01, 543 double *cro10,double *cro11,int *jstat) 544 #else 545 static void sh1461_s9coef(evec1,evec2,evec3,evec4,idim,cro00,cro01, 546 cro10,cro11,jstat) 547 int idim,*jstat; 548 double evec1[],evec2[],evec3[],evec4[]; 549 double *cro00,*cro01,*cro10,*cro11; 550 #endif 551 /* 552 ********************************************************************* 553 * 554 * PURPOSE : Compute factors of expression to find derivatives of 555 * patches in the midpoint of the vertex region. 556 * 557 * 558 * 559 * INPUT : evec1 - Derivative in first parameter direction of first patch. 560 * evec2 - Derivative in second parameter direction of first patch. 561 * evec3 - Derivative in first parameter direction of current patch. 562 * evec4 - Derivative in second parameter direction of current patch. 563 * NB! All these derivative vectors are expected to lie in 564 * the same plane. 565 * idim - Dimension of geometry space. 566 * 567 * 568 * OUTPUT : cro00 - First factor. 569 * cro01 - Second factor. 570 * cro10 - Third factor. 571 * cro11 - Fourth factor. 572 * jstat - status messages 573 * > 0 : warning 574 * = 0 : ok 575 * < 0 : error 576 * 577 * 578 * CALLS : s6length - Length of vector. 579 * s6scpr - Scalar product between two vectors. 580 * s6crss - Cross product between two vectors. 581 * 582 ********************************************************************* 583 */ 584 { 585 int kstat = 0; /* Status variable. */ 586 int ksin,ksin1,ksin2,ksin3,ksin4; /* Sign of sinus of angles. */ 587 double tang,tang1,tang2,tang3,tang4; /* Angles between input vectors. */ 588 double tl1,tl2,tl3,tl4; /* Lengths of input vectors. */ 589 double tsin,tsin1,tsin2,tsin3,tsin4; /* Sinus of angles. */ 590 double snorm[3]; /* Normal of plane spanned by 591 tangent vectors. */ 592 double svec[3]; /* Vector in tangent plane normal 593 to a specific tangent vector. */ 594 double tlvec; /* Length of the vector svec. */ 595 596 /* Compute the normal to the tangent plane spanned by the input vectors. */ 597 598 s6crss(evec1,evec2,snorm); 599 600 /* Compute the lengths of the vectors evec1 to evec4. */ 601 602 tl1 = s6length(evec1,idim,&kstat); 603 tl2 = s6length(evec2,idim,&kstat); 604 tl3 = s6length(evec3,idim,&kstat); 605 tl4 = s6length(evec4,idim,&kstat); 606 607 /* Compute the vector lying in the tangent plane normal to evec1. */ 608 609 s6crss(snorm,evec1,svec); 610 tlvec = s6length(svec,idim,&kstat); 611 612 /* Compute the sinus of angles between evec1 and the other input vectors. */ 613 614 tsin = s6scpr(svec,evec2,idim)/(tlvec*tl2); 615 tsin1 = s6scpr(svec,evec3,idim)/(tlvec*tl3); 616 tsin2 = s6scpr(svec,evec4,idim)/(tlvec*tl4); 617 618 /* Compute the vector lying in the tangent plane normal to evec2. */ 619 620 s6crss(snorm,evec2,svec); 621 tlvec = s6length(svec,idim,&kstat); 622 623 /* Compute the sinus of angles between evec2 and the later input vectors. */ 624 625 tsin3 = s6scpr(svec,evec3,idim)/(tlvec*tl3); 626 tsin4 = s6scpr(svec,evec4,idim)/(tlvec*tl4); 627 628 /* Fetch the sign of the sinuses of the angles. */ 629 630 ksin = (tsin < DZERO) ? -1 : 1; 631 ksin1 = (tsin1 < DZERO) ? -1 : 1; 632 ksin2 = (tsin2 < DZERO) ? -1 : 1; 633 ksin3 = (tsin3 < DZERO) ? -1 : 1; 634 ksin4 = (tsin4 < DZERO) ? -1 : 1; 635 636 /* Compute the angles in a more stable way. */ 637 638 tang = sh1461_s9ang(evec1,evec2,idim); 639 tang1 = sh1461_s9ang(evec1,evec3,idim); 640 tang2 = sh1461_s9ang(evec1,evec4,idim); 641 tang3 = sh1461_s9ang(evec2,evec3,idim); 642 tang4 = sh1461_s9ang(evec2,evec4,idim); 643 644 /* Compute the sinuses of the angles. */ 645 646 tsin = ksin*sin(tang); 647 tsin1 = ksin1*sin(tang1); 648 tsin2 = ksin2*sin(tang2); 649 tsin3 = ksin3*sin(tang3); 650 tsin4 = ksin4*sin(tang4); 651 652 /* Compute the output factors. */ 653 654 *cro00 = (tl3*tsin1)/(tl2*tsin); 655 *cro01 = (tl4*tsin2)/(tl2*tsin); 656 *cro10 = -(tl3*tsin3)/(tl1*tsin); 657 *cro11 = -(tl4*tsin4)/(tl1*tsin); 658 659 *jstat = 0; 660 661 return; 662 } 663 664 #if defined(SISLNEEDPROTOTYPES) 665 static void sh1461_s9chcoor(double eblend[],int iordblend,double epos[],int iordpos,double ecrt1[],int iordcrt1,double evec1[],double evec2[],double evec3[],int idim,double ecrt2[],int * jordcrt2,int * jstat)666 sh1461_s9chcoor(double eblend[],int iordblend,double epos[], 667 int iordpos,double ecrt1[],int iordcrt1, 668 double evec1[],double evec2[],double evec3[], 669 int idim,double ecrt2[],int *jordcrt2,int *jstat) 670 #else 671 static void sh1461_s9chcoor(eblend,iordblend,epos,iordpos,ecrt1,iordcrt1, 672 evec1,evec2,evec3,idim,ecrt2,jordcrt2,jstat) 673 int iordblend,iordpos,iordcrt1,idim,*jordcrt2,*jstat; 674 double evec1[],evec2[],evec3[]; 675 double eblend[],epos[],ecrt1[],ecrt2[]; 676 #endif 677 /* 678 ********************************************************************* 679 * 680 * PURPOSE : Compute that cross tangent curve belonging to one of the 681 * blending patches, that is a mapping of the derivative 682 * curves of the previous patch. 683 * 684 * 685 * 686 * INPUT : eblend - Vertices of blending function. The dimension 687 * of the array is iordblend. 688 * iordblend - Number of vertices of blending function. 689 * iordblend = 4. 690 * epos - Vertices of position curve. This curve is 691 * between the corners (0,0) and (0,1) of the 692 * previous patch, and the corners (0,0) and 693 * (1,0) of the current patch. 694 * iordpos - Number of vertices of the position curve. 695 * iordpos = 5. 696 * ecrt1 - Vertices of input cross tangent curve. This 697 * curve is the derivative curve in the 1. parameter 698 * direction along the edge from (0,0) to (0,1) 699 * of the previous patch. 700 * iordcrt1 - Number of vertices of input cross tangent curve. 701 * iordcrt1 = 4. 702 * evec1 - Tangent vector at the corner (0,0) in the 1. 703 * parameter direction of the previous patch. 704 * evec2 - Tangent vector at the corner (0,0) in the 2. 705 * parameter direction of the previous patch, and 706 * in the 1. parameter direction of the current patch. 707 * evec3 - Tangent vector at the corner (0,0) in the 2. 708 * parameter direction of the current patch. 709 * idim - Dimension of geometry space. 710 * 711 * 712 * OUTPUT : ecrt2 - Vertices of the produced cross tangent curve. 713 * This curve is the derivative curve in the 2. 714 * parameter direction along the edge from (0,0) to 715 * (1,0) of the current patch. 716 * jordcrt2 - Number of vertices of produced cross tangent curve. 717 * *jordcrt2 = 7. 718 * jstat - status messages 719 * > 0 : warning 720 * = 0 : ok 721 * < 0 : error 722 * 723 * 724 * METHOD : Let pc1 be the input cross tangent curve, pc2 be the 725 * derivative curve along the position curve and alpha the 726 * blending function. Then the output cross tangent curve is 727 * given by : 728 * pc3(s) = ((my+1)*alpha(s)-1)*pc1(s) + lambda*alpha(s)*pc2(s) 729 * my and lambda are constants depending on the lengths of 730 * input vectors, evec1, evec2 and evec3, and the angles 731 * between them. See the paper of Hahn. 732 * 733 * CALLS : s6length - Length of vector. 734 * 735 ********************************************************************* 736 */ 737 { 738 int kstat; /* Status variable. */ 739 int ki; /* Counters. */ 740 double tl1,tl2,tl3; /* Lengths of the input vectors. */ 741 double tang1,tang2; /* Angles between input vectors. */ 742 double tsin1,tsin2,tsin3; /* Sinus of angles between input 743 vectors. */ 744 double tmy,tmy1; /* Coefficients dependant on the 745 input vectors. */ 746 double tlambda; /* Coefficient dependant on the 747 input vectors. */ 748 double scoef[12]; /* Vertices of the derivative 749 curve along the position curve. */ 750 double sc1[21],sc2[21]; /* Coefficients of products of curves. */ 751 double sblend2[4]; /* Coefficients of curve equal to 752 a factor times the blending 753 curve minus 1. */ 754 755 if (iordblend != 4 || iordcrt1 != 4) goto err002; 756 757 *jordcrt2 = iordblend + iordcrt1 - 1; 758 759 /* Compute constants. */ 760 761 tl1 = s6length(evec1,idim,&kstat); 762 tl2 = s6length(evec2,idim,&kstat); 763 tl3 = s6length(evec3,idim,&kstat); 764 tang1 = sh1461_s9ang(evec1,evec2,idim); 765 tang2 = sh1461_s9ang(evec2,evec3,idim); 766 tsin1 = sin(tang1); 767 tsin2 = sin(tang2); 768 tsin3 = sin(tang1+tang2); 769 tmy = - (tl3*tsin2)/(tl1*tsin1); 770 tlambda = (tl3*tsin3)/(tl2*tsin1); 771 tmy1 = tmy + (double)1.0; 772 773 /* Compute coefficients of the curve (my+1)alpha(s)-1. */ 774 775 for (ki=0; ki<4; ki++) sblend2[ki] = tmy1*eblend[ki] - (double)1.0; 776 777 /* Compute coefficients of derivative curve along position curve. */ 778 779 for (ki=0; ki<4*idim; ki++) 780 scoef[ki] = (double)4.0*(epos[idim+ki] - epos[ki]); 781 782 /* Compute first part of the second cross tangent curve. */ 783 784 sh1461_s9mult(eblend,scoef,4,idim,sc1,&kstat); 785 786 /* Compute second part of the second cross tangent curve. */ 787 788 sh1461_s9mult(sblend2,ecrt1,4,idim,sc2,&kstat); 789 790 /* Compute cross tangent curve. */ 791 792 for (ki=0; ki<7*idim; ki++) 793 ecrt2[ki] = (tlambda*sc1[ki] + sc2[ki]); 794 795 *jstat = 0; 796 goto out; 797 798 /* Error in input. */ 799 800 err002 : 801 *jstat = -2; 802 goto out; 803 804 out : 805 return; 806 } 807 808 809 #if defined(SISLNEEDPROTOTYPES) 810 static void sh1461_s9mult(double eblend[],double ecoef[],int iord,int idim,double ecoefnew[],int * jstat)811 sh1461_s9mult(double eblend[],double ecoef[],int iord, 812 int idim,double ecoefnew[],int *jstat) 813 #else 814 static void sh1461_s9mult(eblend,ecoef,iord,idim,ecoefnew,jstat) 815 int iord,idim,*jstat; 816 double eblend[],ecoef[],ecoefnew[]; 817 #endif 818 /* 819 ********************************************************************* 820 * 821 * PURPOSE : Compute the product of two Bezier curves of order 4, 822 * one of which is a blending function of dimension 1. 823 * 824 * 825 * 826 * INPUT : eblend - Vertices of blending function. The dimension 827 * of the array is iord, i.e. 4. 828 * ecoef - Vertices of the other Bezier curve. The 829 * dimension of the array is iord*idim, i.e. 4*idim. 830 * iord - Order of the Bezier curves. iord = 4. 831 * idim - Dimension of geometry space. 832 * 833 * 834 * OUTPUT : ecoefnew - Vertices of the product curve. The array is 835 * allocated outside this routine, and has 836 * dimension (2*iord-1)*idim, i.e. 7*idim 837 * jstat - status messages 838 * > 0 : warning 839 * = 0 : ok 840 * < 0 : error 841 * 842 * USE : Used in Hahn's method when computing the cross tangent 843 * curve which is a mapping of the derivative curves of the 844 * previous patch. Called from s9chcoor. 845 * 846 ********************************************************************* 847 */ 848 { 849 int ki; /* Index. */ 850 851 /* Test if the order of the curves is equal to 4. */ 852 853 if (iord != 4) goto err001; 854 855 for (ki=0; ki<idim; ki++) 856 { 857 /* Compute the vertices of the product curve. */ 858 859 ecoefnew[ki] = eblend[0]*ecoef[ki]; 860 ecoefnew[idim+ki] = (eblend[0]*ecoef[idim+ki] 861 + eblend[1]*ecoef[ki])/(double)2.0; 862 ecoefnew[2*idim+ki] = (eblend[0]*ecoef[2*idim+ki] 863 + (double)3.0*eblend[1]*ecoef[idim+ki] 864 + eblend[2]*ecoef[ki])/(double)5.0; 865 ecoefnew[3*idim+ki] = (eblend[0]*ecoef[3*idim+ki] + eblend[3]*ecoef[ki] 866 + (double)9.0*(eblend[1]*ecoef[2*idim+ki] + 867 eblend[2]*ecoef[idim+ki]))/(double)20.0; 868 ecoefnew[4*idim+ki] = (eblend[1]*ecoef[3*idim+ki] 869 + (double)3.0*eblend[2]*ecoef[2*idim+ki] 870 + eblend[3]*ecoef[idim+ki])/(double)5.0; 871 ecoefnew[5*idim+ki] = (eblend[2]*ecoef[3*idim+ki] 872 + eblend[3]*ecoef[2*idim+ki])/(double)2.0; 873 ecoefnew[6*idim+ki] = eblend[3]*ecoef[3*idim+ki]; 874 } 875 876 *jstat = 0; 877 goto out; 878 879 /* Error in input. The order of the curves is not equal to 4. */ 880 881 err001 : 882 *jstat = -1; 883 goto out; 884 885 out : 886 return; 887 } 888 889 #if defined(SISLNEEDPROTOTYPES) 890 static double sh1461_s9ang(double evec1[],double evec2[],int idim)891 sh1461_s9ang(double evec1[],double evec2[],int idim) 892 #else 893 static double sh1461_s9ang(evec1,evec2,idim) 894 double evec1[]; 895 double evec2[]; 896 int idim; 897 #endif 898 /* 899 ********************************************************************* 900 * 901 ********************************************************************* 902 * 903 * PURPOSE : Compute the angle (in radians) between two vectors 904 * 905 * 906 * 907 * INPUT : evec1 - First vector 908 * evec2 - Second vector 909 * idim - Dimension of the space in which the vectors lie. 910 * 911 * 912 * 913 * OUTPUT : sh1461_s9ang - Angle in radians between vectors 914 * 915 * 916 * METHOD : Make cosine of the angle by computing the scalar product, 917 * then divide by the length of the two vectors. 918 * 919 * REFERENCES : 920 * 921 *- 922 * CALLS : s6scpr - Scalar product between two vectors. 923 * s6length - Length of vector. 924 * 925 * WRITTEN BY : Tor Dokken SI, 88-07. 926 * Arne Laksaa SI, 89-07. 927 * Vibeke Skytt SI, 90-04. 928 * 929 ********************************************************************* 930 */ 931 { 932 double tscpr,tang,tlength1,tlength2,tcos; 933 int kstat1,kstat2; 934 935 tscpr = s6scpr(evec1,evec2,idim); 936 937 tlength1 = s6length(evec1,idim,&kstat1); 938 tlength2 = s6length(evec2,idim,&kstat2); 939 940 if (!kstat1 || !kstat2) 941 tang = DZERO; 942 else 943 { 944 tcos = tscpr/(tlength1*tlength2); 945 tcos = MIN((double)1.0,tcos); 946 tcos = MAX(-(double)1.0,tcos); 947 tang = acos(tcos); 948 } 949 950 return(tang); 951 } 952 953 #if defined(SISLNEEDPROTOTYPES) 954 static void sh1461_s9comder(int ider1,int ider2,double ederprev[],int idim,double aro00,double aro01,double aro10,double aro11,double eder[],int * jstat)955 sh1461_s9comder(int ider1,int ider2,double ederprev[],int idim, 956 double aro00,double aro01,double aro10, 957 double aro11,double eder[],int *jstat) 958 #else 959 static void sh1461_s9comder(ider1,ider2,ederprev,idim,aro00,aro01,aro10, 960 aro11,eder,jstat) 961 int ider1,ider2,idim,*jstat; 962 double ederprev[],aro00,aro01,aro10,aro11,eder[]; 963 #endif 964 /* 965 ********************************************************************* 966 * 967 * PURPOSE : Compute given 2. derivative of a blending surface at the 968 * midpoint of the vertex region when the 2. derivatives of the 969 * first blending patch are known. 970 * 971 * 972 * 973 * INPUT : ider1 - Order of differentiation in first parameter direction. 974 * ider2 - Order of differentiation in second parameter direction. 975 * ederprev - 2. derivatives of 1. patch stored in the following 976 * order, (2,0)-, (1,1)- and (0,2)-derivative. 977 * idim - Dimension of geometry space. 978 * aro00 - First factor of coordinate transformation. 979 * aro01 - Second factor of coordinate transformation. 980 * aro10 - Third factor of coordinate transformation. 981 * aro11 - Fourth factor of coordinate transformation. 982 * 983 * 984 * OUTPUT : eder - Actual 2. derivative. 985 * jstat - status messages 986 * > 0 : warning 987 * = 0 : ok 988 * < 0 : error 989 * 990 * 991 ********************************************************************* 992 */ 993 { 994 int kj,kk,kh; /* Counters. */ 995 double t00,t01,t10,t11; /* Coefficients used to compute 996 derivatives at midpoint of region. */ 997 double tfac1,tfac2; /* Factor used to compute derivatives 998 at midpoint of region. */ 999 1000 /* Test input. */ 1001 1002 if (ider1 + ider2 != 2) goto err001; 1003 1004 /* Compute requested 2. derivative. */ 1005 1006 t00 = t10 = (double)1.0; 1007 for (kj=0; kj<ider1; kj++) t00 *= aro00; 1008 1009 for (kj=0; kj<=ider1; kj++) 1010 { 1011 tfac1 = (ider1 > 0) ? (double)((kj % ider1) + 1) : (double)1.0; 1012 1013 t01 = t11 = (double)1.0; 1014 for (kk=0; kk<ider2; kk++) t01*=aro01; 1015 1016 for (kk=0; kk<=ider2; kk++) 1017 { 1018 tfac2 = (ider2 > 0) ? (double)((kk % ider2) + 1) : (double)1.0; 1019 1020 for (kh=0; kh<idim; kh++) 1021 eder[kh] += tfac1*tfac2*t00*t01*t10*t11*ederprev[(2-kj-kk)*idim+kh]; 1022 1023 if (aro01 != DZERO) t01 /= aro01; 1024 else if (kk == ider2-1) t01 = (double)1.0; 1025 t11 *= aro11; 1026 } 1027 if (aro00 != DZERO) t00 /= aro00; 1028 else if (kj == ider1-1) t00 = (double)1.0; 1029 t10 *= aro10; 1030 } 1031 1032 *jstat = 0; 1033 goto out; 1034 1035 /* Error in input. Wrong order of differentiation. */ 1036 1037 err001 : 1038 *jstat = -1; 1039 goto out; 1040 1041 out : 1042 return; 1043 } 1044