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: s2510.c,v 1.4 2001-03-19 15:58:59 afr Exp $ 45 * 46 */ 47 48 49 #define S2510 50 51 #include "sislP.h" 52 53 #if defined(SISLNEEDPROTOTYPES) 54 void s2510(SISLSurf * surf,int ider,int iside1,int iside2,double parvalue[],int * leftknot1,int * leftknot2,double * mehlum,int * jstat)55 s2510(SISLSurf *surf, int ider, int iside1, int iside2, double parvalue[], 56 int *leftknot1,int *leftknot2, double *mehlum, int *jstat) 57 #else 58 void s2510(surf, ider, iside1, iside2, parvalue, leftknot1, leftknot2, mehlum, jstat) 59 SISLSurf *surf; 60 int ider; 61 int iside1; 62 int iside2; 63 double parvalue[]; 64 int *leftknot1; 65 int *leftknot2; 66 double *mehlum; 67 int *jstat; 68 #endif 69 /* 70 *************************************************************************** 71 * 72 *************************************************************************** 73 * PURPOSE : To compute the third order Mehlum curvature M(u,v) of a 74 * surface for given values (u,v) = (parvalue[0],parvalue[1]), 75 * where: 76 * 77 * et1[leftknot1] <= parvalue[0] < et1[leftknot1+1], 78 * et2[leftknot2] <= parvalue[1] < et2[leftknot2+1]. 79 * 80 * INPUT : 81 * surf - Pointer to the surface to evaluate. 82 * ider - Number of derivatives to calculate. 83 * Only implemented for ider=0. 84 * < 0 : No derivative calculated. 85 * = 0 : Position calculated. 86 * = 1 : Position and first derivative calculated. 87 * etc. 88 * iside1 - Indicator telling if the derivatives in the first 89 * parameter direction is to be calculated from the 90 * left or from the right: 91 * < 0 calculate derivative from the left hand side 92 * >= 0 calculate derivative from the right hand side. 93 * iside2 - Indicator telling if the derivatives in the second 94 * parameter direction is to be calculated from the 95 * left or from the right: 96 * < 0 calculate derivative from the left hand side 97 * >= 0 calculate derivative from the right hand side. 98 * parvalue - Parameter-value at which to evaluate. Dimension of 99 * parvalue is 2. 100 * 101 * INPUT/OUTPUT : 102 * leftknot1 - Pointer to the interval in the knot vector in the 103 * first parameter direction where parvalue[0] is found, 104 * that is: 105 * et1[leftknot1] <= parvalue[0] < et1[leftknot1+1]. 106 * leftknot1 should be set equal to zero at the first call 107 * to the routine. 108 * 109 * leftknot2 - Pointer to the interval in the knot vector in the 110 * second parameter direction where parvalue[1] is found, 111 * that is: 112 * et2[leftknot2] <= parvalue[1] < et2[leftknot2+1]. 113 * leftknot2 should be set equal to zero at the first call 114 * to the routine. 115 * 116 * OUTPUT : 117 * mehlum - Third order Mehlum curvature of the surface in (u,v) = 118 * (parvalue[0],parvalue[1]). 119 * jstat - Status messages 120 * = 2 : Surface is degenerate at the point, that is, 121 * the surface is not regular at this point. 122 * = 1 : Surface is close to degenerate at the point. 123 * Angle between tangents is less than the angular 124 * tolerance. 125 * = 0 : Ok. 126 * < 0 : Error. 127 * 128 * METHOD : The third order Mehlum curvature is given by 129 * 130 * M(u,v) = (5G^3P^2 131 * + (EG + 4F^2) 132 * *(9GQ^2 + 9ES^2 + 6GPS + 6EQT) 133 * + 5E^3T^2 134 * - 2F(3EG + 2F^2)(PT + 9QS) 135 * - 30F(G^2PQ + E^2ST)) 136 * /(16(EG - F^2)^3). 137 * 138 * The variables E,F,G,P,Q,S and T 139 * are the coefficients of the first and third fundamental form. 140 * They are given by: 141 * E = <Xu,Xu>, F = <Xu,Xv> and G = <Xv,Xv>. 142 * P = <N,Xuuu> + 3(a*alpha + b*beta), 143 * Q = <N,Xuuv> + c*alpha + d*beta + 2a*gamma + 2b*delta, 144 * S = <N,Xuvv> + 2c*gamma + 2d*delta + a*epsilon + b*mu, 145 * T = <N,Xvvv> + 3(c*epsilon + d*mu), 146 * where N is normalized, and 147 * a = Ff - Ge, 148 * b = Fe - Ef, 149 * c = Fg - Gf, 150 * d = Ff - Eg, 151 * e, f and g being the second fundamental form coefficients 152 * (e = <N,Xuu>, f = <N,Xuv> and g = <N,Xvv>), and 153 * alpha = <Xuu,Xu>/||N||^2, 154 * beta = <Xuu,Xv>/||N||^2, 155 * gamma = <Xuv,Xu>/||N||^2, 156 * delta = <Xuv,Xv>/||N||^2, 157 * epsilon = <Xvv,Xu>/||N||^2, 158 * mu = <Xvv,Xv>/||N||^2. 159 * 160 * REFERENCES : Differential Geometry of Curves and Surfaces, 161 * (Manfredo P. Do Carmo, Prentice Hall, 162 * ISBN: 0-13-212589-7). 163 *- 164 * CALLS : s1422() and s2511(). 165 * 166 * LIMITATIONS : 167 * (i) If the surface is degenerated (not regular) at the point 168 * (u,v), it makes no sense to speak about the Mehlum 169 * M(u,v). The routine returns jstat = 2. 170 * (ii) If the surface is closed to degenerate, the Mehlum 171 * M(u,v) can be numerical instable. The routine returns 172 * jstat = 1. 173 * (iii) The dimension of the space in which the surface lies must 174 * be 1,2 or 3. 175 * 176 * 177 * WRITTEN BY : Johannes Kaasa, SINTEF, Oslo, Norway. Date: 1995-9 178 ****************************************************************************** 179 */ 180 { 181 int kwarn = 0; /* Local staus variable(warning). */ 182 int kistat = 0; /* Local staus variable. */ 183 double derive[30]; /* Array containing the computed derivatives. */ 184 double normal[3]; /* Array containing the computed normalvektor. */ 185 186 187 if (ider != 0) goto err178; 188 189 190 if (surf == SISL_NULL) goto err150; 191 else 192 { 193 /* Compute derivates and normal. */ 194 195 s1422(surf,3,iside1,iside2,parvalue,leftknot1,leftknot2,derive,normal, 196 &kistat); 197 if (kistat > 0) kwarn=kistat; 198 199 if (kistat < 0) /* Error in lower level routine. */ 200 { 201 goto error; 202 } 203 else if (kistat != 2) /* The surface is not degenerate */ 204 { 205 s2511(surf, ider, derive, normal, mehlum, &kistat); 206 207 if (kistat < 0) 208 goto error; 209 } 210 else if (kistat == 2) /* The surface is degenerated. */ 211 { 212 *mehlum = 0.0; 213 goto war002; 214 } 215 216 } 217 218 219 /* Successful computations */ 220 221 *jstat = kwarn; 222 goto out; 223 224 225 /* The surface is degenerated at (u,v) */ 226 war002: 227 *jstat = 2; 228 goto out; 229 230 /* Error. Input (surface) pointer is SISL_NULL. */ 231 err150: 232 *jstat = -150; 233 s6err("s2510", *jstat, 0); 234 goto out; 235 236 /* Illegal derivative requested. */ 237 err178: 238 *jstat = -178; 239 s6err("s2510",*jstat,0); 240 goto out; 241 /* Error in lower level routine. */ 242 243 error: 244 *jstat = kistat; 245 s6err("s2510",*jstat,0); 246 goto out; 247 248 249 out: 250 251 return; 252 253 } 254