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: sh1466.c,v 1.1 1994-04-21 12:10:42 boh Exp $ 45 * 46 */ 47 48 49 #define SH1466 50 51 #include "sislP.h" 52 53 #if defined(SISLNEEDPROTOTYPES) 54 void sh1466(SISLCurve * ecurve[],double etwist[],int ider,double ebar[],double eval[],int * jstat)55 sh1466(SISLCurve *ecurve[],double etwist[],int ider,double ebar[], 56 double eval[],int *jstat) 57 #else 58 void sh1466(ecurve,etwist,ider,ebar,eval,jstat) 59 double etwist[],ebar[],eval[]; 60 SISLCurve *ecurve[]; 61 int ider,*jstat; 62 #endif 63 /* 64 ********************************************************************* 65 * 66 * PURPOSE : Given the barycentric coordinates of a point in a 3-sided 67 * vertex region, evaluate the value of the ideal blending 68 * surface of the vertex region in this point. 69 * 70 * 71 * 72 * INPUT : ecurve - Position and cross-tangent curves around the vertex 73 * region. For each edge of the region position and cross- 74 * tangent curves are given. The curves follow each other 75 * around the region and are oriented counter-clock-wise. 76 * The dimension of the array is 6. 77 * etwist - Twist-vectors of the corners of the vertex region. The 78 * first element of the array is the twist in the corner 79 * before the first edge, etc. The dimension of the array 80 * is 3*kdim. 81 * ider - Number of derivatives to compute. Directions of 82 * differentiation is that of the two first barycentric 83 * coordinates. 0 <= ider <= 2. 84 * ebar - Barycentric coordinates of the point to be evaluated. 85 * The dimension of the array is 3. 86 * 87 * 88 * OUTPUT : eval - Value and derivatives of ideal blending surface in the 89 * given point. Dimension of the array is 3*(1+..+(ider+1)). 90 * jstat - status messages 91 * > 0 : warning 92 * = 0 : ok 93 * < 0 : error 94 * 95 * 96 * METHOD : A functional description of the ideal surface is given as 97 * a blend between three surfaces, each of which fulfill the 98 * continuity requirements over two edges. 99 * 100 * REFERENCES : Gregory and Charrot : A C1 Triangular Interpolation Patch for 101 * Computer Aided Geometric Design 102 * 103 * USE : 3D geometry only. 104 * 105 *- 106 * CALLS : s1221 - Evaluate curve at a given parameter value. 107 * 108 * WRITTEN BY : Vibeke Skytt, SI, Dec 89. 109 * 110 ********************************************************************* 111 */ 112 { 113 int kstat=0; /* Status variable. */ 114 int ki,kj,kk,kh; /* Counters. */ 115 int kder; /* Number of derivatives to evaluate. */ 116 int kleft = 0; /* Local parameter used in s1221. */ 117 int kdim = 3; /* Dimension of geometry. */ 118 int kwarn = 0; /* Indicates if a warning is to be sendt. */ 119 int knmb; /* Number of doubles pr derivative. */ 120 int kl = 0; /* Number of derivatives in the output array. */ 121 double tpar1; /* Parameter value of edge between actual corner 122 and next corner. */ 123 double tpar2; /* Parameter value of edge between actual corner 124 and previous corner. */ 125 double tlambi; /* First barycentric coordinate. */ 126 double tlambj; /* Second barycentric coordinate. */ 127 double tlambk; /* Third barycentric coordinate. */ 128 double salpha[3]; /* Weight of actual blending surface. */ 129 double sp[9]; /* Value of function interpolating two sides. */ 130 double sstart[3]; /* Start-parameters of edge-curves. */ 131 double send[3]; /* End-parameters of edge-curves. */ 132 double sint[3]; /* Parameter intervals of edge-curves. */ 133 double spos1[27]; /* Position of edge-curve in tpar1. */ 134 double sder1[27]; /* Cross tangent in tpar1. */ 135 double spos2[27]; /* Position of edge-curve in tpar2. */ 136 double sder2[27]; /* Cross tangent in tpar2. */ 137 double scorn[9]; /* Position of edge-curves in actual corner of 138 vertex region. */ 139 double scornder1[9]; /* Tangent in corner along next edge. */ 140 double scornder2[9]; /* Tangent in corner along previous edge. */ 141 142 /* Test input. */ 143 144 if (ider > 2) kwarn = 1; 145 146 /* Initialise. */ 147 148 kder = ider; 149 knmb = kdim*(ider+1); 150 for (ki=0; ki<ider; ki++) kl += ider + 1; 151 152 /* Initiate output array to zero. */ 153 154 for (kh=0; kh<kl*kdim; kh++) eval[kh] = (double)0.0; 155 156 /* Get endpoints of parameter intervals of edge curves. */ 157 158 for (ki=0; ki<3; ki++) 159 { 160 sstart[ki] = *(ecurve[2*ki] -> et + ecurve[2*ki] -> ik - 1); 161 send[ki] = *(ecurve[2*ki] -> et + ecurve[2*ki] -> in); 162 sint[ki] = send[ki] - sstart[ki]; 163 } 164 165 /* Evaluate position and cross-tangent curves at points on 166 the edges needed when evaluating surface. */ 167 168 for (ki=0; ki<3; ki++) 169 { 170 kj = (ki+1) % 3; 171 kk = (ki+2) % 3; 172 173 /* Copy barycentric coordinates to local variables. */ 174 175 tlambi = ebar[ki]; 176 tlambj = ebar[kj]; 177 tlambk = ebar[kk]; 178 179 /* Find parameter values of points on edges used to evaluate 180 actual blending surface. */ 181 182 tpar1 = ((double)1.0 - tlambj)*sstart[ki] + tlambj*send[ki]; 183 tpar2 = tlambk*sstart[kk] + ((double)1.0 - tlambk)*send[kk]; 184 185 /* Evaluate position and cross-tangent curves at first 186 found parameter value. */ 187 188 s1221(ecurve[2*ki],kder,tpar1,&kleft,spos1+knmb*ki,&kstat); 189 if (kstat < 0) goto error; 190 191 s1221(ecurve[2*ki+1],kder,tpar1,&kleft,sder1+knmb*ki,&kstat); 192 if (kstat < 0) goto error; 193 194 /* Evaluate position and cross-tangent curves at second 195 found parameter value. */ 196 197 s1221(ecurve[2*kk],kder,tpar2,&kleft,spos2+knmb*ki,&kstat); 198 if (kstat < 0) goto error; 199 200 s1221(ecurve[2*kk+1],kder,tpar2,&kleft,sder2+knmb*ki,&kstat); 201 if (kstat < 0) goto error; 202 203 /* Evaluate position and both cross-tangents at the corner nr ki. */ 204 205 s1221(ecurve[2*ki],0,sstart[ki],&kleft,scorn+ki*kdim,&kstat); 206 if (kstat < 0) goto error; 207 208 s1221(ecurve[2*kk+1],0,send[kk],&kleft,scornder1+ki*kdim,&kstat); 209 if (kstat < 0) goto error; 210 211 s1221(ecurve[2*ki+1],0,sstart[ki],&kleft,scornder2+ki*kdim,&kstat); 212 if (kstat < 0) goto error; 213 214 /* Compute the weigth of the actual blending surface. */ 215 216 salpha[ki] = tlambi*tlambi*((double)3.0 - (double)2.0*tlambi + 217 (double)6.0*tlambj*tlambk); 218 219 /* Add the contribution of the value of this blending surface to the 220 value of the ideal surface. */ 221 222 for (kh=0; kh<kdim; kh++) 223 { 224 sp[ki*kdim+kh] = spos1[knmb*ki+kh] + tlambk*sder1[knmb*ki+kh] + 225 spos2[knmb*ki+kh] + tlambj*sder2[knmb*ki+kh] - 226 scorn[ki*kdim+kh] - tlambj*scornder1[ki*kdim+kh] 227 - tlambk*scornder2[ki*kdim+kh] - 228 tlambj*tlambk*etwist[ki*kdim+kh]; 229 230 eval[kh] += salpha[ki]*sp[ki*kdim+kh]; 231 } 232 } 233 234 if (ider >= 1) 235 { 236 /* Compute first derivatives of the Gregory Charrot function. */ 237 238 double tl1,tl2,tl3; /* Barycentric coordinates. */ 239 double sd1alpha[3],sd2alpha[3]; /* 1. derivative of weight functions. */ 240 double sd1p[9],sd2p[9]; /* 1. derivative of blending functions. */ 241 242 /* Copy barycentric coordinates to local variables. */ 243 244 tl1 = ebar[0]; 245 tl2 = ebar[1]; 246 tl3 = ebar[2]; 247 248 /* Compute the 1. derivatives of the weight functions. */ 249 250 sd1alpha[0] = (double)6.0*tl1*((double)1.0 - tl1 - tl1*tl2 + 251 (double)2.0*tl2*tl3); 252 sd2alpha[0] = (double)6.0*tl1*tl1*(tl3 - tl2); 253 254 sd1alpha[1] = (double)6.0*tl2*tl2*(tl3 - tl1); 255 sd2alpha[1] = (double)6.0*tl2*((double)1.0 - tl2 - tl1*tl2 + 256 (double)2.0*tl1*tl3); 257 258 sd1alpha[2] = (double)6.0*tl3*(-(double)1.0 + tl3 + tl2*tl3 - 259 (double)2.0*tl1*tl2); 260 sd2alpha[2] = (double)6.0*tl3*(-(double)1.0 + tl3 + tl1*tl3 - 261 (double)2.0*tl1*tl2); 262 263 /* Compute 1. derivatives of the functions which blends two sides 264 of the region. */ 265 266 for (kh=0; kh<kdim; kh++) 267 { 268 sd1p[kh] = (spos2[kdim+kh] + tl2*sder2[kdim+kh])*sint[2] 269 - sder1[kh] + scornder2[kh] + tl2*etwist[kh]; 270 sd2p[kh] = (spos2[kdim+kh] + tl2*sder2[kdim+kh])*sint[2] 271 + (spos1[kdim+kh] + tl3*sder1[kdim+kh])*sint[0] 272 + sder2[kh] - sder1[kh] - scornder1[kh] 273 + scornder2[kh] + (tl2 - tl3)*etwist[kh]; 274 275 sd1p[kdim+kh] = -sder2[knmb+kh] + sder1[knmb+kh] 276 - (spos2[knmb+kdim+kh] + tl3*sder2[knmb+kdim+kh])*sint[0] 277 - (spos1[knmb+kdim+kh] + tl1*sder1[knmb+kdim+kh])*sint[1] 278 + scornder1[kdim+kh] - scornder2[kdim+kh] 279 + (tl1 - tl3)*etwist[kdim+kh]; 280 sd2p[kdim+kh] = - sder2[knmb+kh] 281 - (spos1[knmb+kdim+kh] + tl1*sder1[knmb+kdim+kh])*sint[1] 282 + scornder1[kdim+kh] + tl1*etwist[kdim+kh]; 283 284 sd1p[2*kdim+kh] = sder2[2*knmb+kh] 285 + (spos1[2*knmb+kdim+kh] + tl2*sder1[2*knmb+kdim+kh])*sint[2] 286 - scornder1[2*kdim+kh] - tl2*etwist[2*kdim+kh]; 287 sd2p[2*kdim+kh] = sder1[2*knmb+kh] 288 - (spos2[2*knmb+kdim+kh] + tl1*sder2[2*knmb+kdim+kh])*sint[1] 289 - scornder2[2*kdim+kh] - tl1*etwist[2*kdim+kh]; 290 291 /* Compute the first derivative of the Gregory Charrot function. */ 292 293 for (ki=0; ki<3; ki++) 294 { 295 eval[kdim+kh] += sd1alpha[ki]*sp[ki*kdim+kh] 296 + salpha[ki]*sd1p[ki*kdim+kh]; 297 298 eval[2*kdim+kh] += sd2alpha[ki]*sp[ki*kdim+kh] 299 + salpha[ki]*sd2p[ki*kdim+kh]; 300 } 301 } 302 303 if (ider >= 2) 304 { 305 double sd11alpha[3],sd12alpha[3],sd22alpha[3]; /* 2. derivatives of 306 blending patches. */ 307 double sd11p[9],sd12p[9],sd22p[9]; /* 2. derivatives of weight function. */ 308 309 /* Compute second derivatives of the Gregory Charrot function. */ 310 311 /* Compute the 2. derivatives of the weight functions. */ 312 313 sd11alpha[0] = (double)6.0 - (double)12.0*tl1 314 - (double)24.0*tl1*tl2 + (double)12.0*tl2*tl3; 315 sd12alpha[0] = tl1*((double)12.0 - (double)18.0*tl1 316 - (double)24.0*tl2); 317 sd22alpha[0] = -(double)12.0*tl1*tl1; 318 319 sd11alpha[1] = -(double)12.0*tl2*tl2; 320 sd12alpha[1] = tl2*((double)12.0 - (double)18.0*tl2 321 - (double)24.0*tl1); 322 sd22alpha[1] = (double)6.0 - (double)12.0*tl2 323 - (double)24.0*tl1*tl2 + (double)12.0*tl1*tl3; 324 325 sd11alpha[2] = (double)6.0 - (double)12.0*tl3 326 - (double)24.0*tl2*tl3 + (double)12.0*tl1*tl2; 327 sd12alpha[2] = (double)6.0 + (double)6.0*tl3*tl3 328 + (double)12.0*(-tl3 + tl1*tl2 - tl1*tl3 - tl2*tl3); 329 sd22alpha[2] = (double)6.0 - (double)12.0*tl3 330 - (double)24.0*tl1*tl3 + (double)12.0*tl1*tl2; 331 332 /* Compute 2. derivatives of the functions which blends two sides 333 of the region. */ 334 335 for (kh=0; kh<kdim; kh++) 336 { 337 sd11p[kh] = (spos2[2*kdim+kh]+ tl2*sder2[2*kdim+kh])*sint[2]*sint[2]; 338 sd12p[kh] = ((spos2[2*kdim+kh] + tl2*sder2[2*kdim+kh])*sint[2] 339 + sder2[kdim+kh]*sint[2] - sder1[kdim+kh])*sint[0] 340 + etwist[kh]; 341 sd22p[kh] = ((spos2[2*kdim+kh] + tl2*sder2[2*kdim+kh])*sint[2] 342 + (double)2.0*sder2[kdim+kh])*sint[2] 343 + ((spos1[2*kdim+kh] + tl3*sder1[2*kdim+kh])*sint[0] 344 - (double)2.0*sder1[kdim+kh])*sint[0] 345 + (double)2.0*etwist[kh]; 346 347 sd11p[kdim+kh] = ((spos2[knmb+2*kdim+kh] 348 + tl3*sder2[knmb+2*kdim+kh])*sint[0] 349 + (double)2.0*sder2[knmb+kdim+kh])*sint[0] 350 + ((spos1[knmb+2*kdim+kh] 351 + tl1*sder1[knmb+2*kdim+kh])*sint[1] 352 - (double)2.0*sder1[knmb+kdim+kh])*sint[1] 353 + (double)2.0*etwist[kdim+kh]; 354 sd12p[kdim+kh] = ((spos1[knmb+2*kdim+kh] 355 + tl1*sder1[knmb+2*kdim+kh])*sint[1] 356 - sder1[knmb+kdim+kh])*sint[1] 357 + sder2[knmb+kdim+kh]*sint[0] 358 + etwist[kdim+kh]; 359 sd22p[kdim+kh] = (spos1[knmb+2*kdim+kh] 360 + tl1*sder1[knmb+2*kdim+kh])*sint[1]*sint[1]; 361 362 sd11p[2*kdim+kh] = (spos1[2*knmb+2*kdim+kh] 363 + tl2*sder1[2*knmb+2*kdim+kh])*sint[2]*sint[2]; 364 sd12p[2*kdim+kh] = -sder2[2*knmb+kdim+kh]*sint[1] 365 + sder1[2*knmb+kdim+kh]*sint[2] - etwist[2*kdim+kh]; 366 sd22p[2*kdim+kh] = (spos2[2*knmb+2*kdim+kh] 367 + tl1*sder2[2*knmb+2*kdim+kh])*sint[1]*sint[1]; 368 369 /* Compute the 2. derivative of the Gregory Charrot function. */ 370 371 for (ki=0; ki<3; ki++) 372 { 373 eval[3*kdim+kh] += sd11alpha[ki]*sp[ki*kdim+kh] 374 + (double)2.0*sd1alpha[ki]*sd1p[ki*kdim+kh] 375 + salpha[ki]*sd11p[ki*kdim+kh]; 376 377 eval[4*kdim+kh] += sd12alpha[ki]*sp[ki*kdim+kh] 378 + sd1alpha[ki]*sd2p[ki*kdim+kh] 379 + sd2alpha[ki]*sd1p[ki*kdim+kh] 380 + salpha[ki]*sd12p[ki*kdim+kh]; 381 382 eval[5*kdim+kh] += sd22alpha[ki]*sp[ki*kdim+kh] 383 + (double)2.0*sd2alpha[ki]*sd2p[ki*kdim+kh] 384 + salpha[ki]*sd22p[ki*kdim+kh]; 385 } 386 } 387 } 388 } 389 390 /* Ideal surface evaluated. */ 391 392 *jstat = kwarn; 393 goto out; 394 395 396 /* Error in lower level function. */ 397 398 error : 399 *jstat = kstat; 400 goto out; 401 402 out : 403 return; 404 } 405 406 407 408 409 410 411