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: s2542.c,v 1.3 2001-03-19 15:59:00 afr Exp $ 45 * 46 */ 47 48 49 #define S2542 50 51 #include "sislP.h" 52 53 #if defined(SISLNEEDPROTOTYPES) 54 void s2542(SISLSurf * surf,int ider,int iside1,int iside2,double parvalue[],int * leftknot1,int * leftknot2,double * k1,double * k2,double d1[],double d2[],int * jstat)55 s2542(SISLSurf *surf, int ider, int iside1, int iside2, double parvalue[], 56 int *leftknot1,int *leftknot2, double *k1, double *k2, 57 double d1[], double d2[],int *jstat) 58 #else 59 void s2542(surf, ider, iside1, iside2, parvalue, leftknot1, leftknot2, 60 k1, k2, d1, d2, jstat) 61 SISLSurf *surf; 62 int ider; 63 int iside1; 64 int iside2; 65 double parvalue[]; 66 int *leftknot1; 67 int *leftknot2; 68 double *k1; 69 double *k2; 70 double d1[]; 71 double d2[]; 72 int *jstat; 73 #endif 74 /* 75 *************************************************************************** 76 * 77 *************************************************************************** 78 * PURPOSE : To compute principal curvature (k1,k2) with corresponding 79 * principal directions (d1,d2) of a surface for 80 * given values (u,v) = (parvalue[0],parvalue[1]), where: 81 * 82 * etl[leftknot1] <= parvalue[0] < etl[leftknot1+1], 83 * et2[leftknot2] <= parvalue[1] < et2[leftknot2+1]. 84 * 85 * INPUT : 86 * surf - Pointer to the surface to evaluate. 87 * ider - Number of derivatives to calculate. 88 * Only implemented for ider=0. 89 * < 0 : No derivative calculated. 90 * = 0 : Principal curvature calculated. 91 * = 1 : Principal curvature and its first derivative 92 * calculated. 93 * etc. 94 * iside1 - Indicator telling if the principal curvature in the first 95 * parameter direction is to be calculated from the 96 * left or from the right: 97 * < 0 calculate principal curvature from the left hand 98 * side 99 * >= 0 calculate principal curvature from the right hand 100 * side. 101 * iside2 - Indicator telling if the principal curvature in the second 102 * parameter direction is to be calculated from the 103 * left or from the right: 104 * < 0 calculate principal curvature from the left hand 105 * side 106 * >= 0 calculate principal curvature from the right hand 107 * side. 108 * parvalue - Parameter-value at which to evaluate. Dimension of 109 * parvalue is 2. 110 * 111 * INPUT/OUTPUT : 112 * leftknot1 - Pointer to the interval in the knot vektor in the 113 * first parameter direction where parvalue[0] is found, 114 * that is: 115 * etl[leftknot1] <= parvalue[0] < etl[leftknot1+1]. 116 * leftknot1 should be set equal to zero at the first call 117 * to the rutine. 118 * 119 * leftknot1 - Pointer to the interval in the knot vektor in the 120 * second parameter direction where parvalue[1] is found, 121 * that is: 122 * et2[leftknot2] <= parvalue[1] < et2[leftknot2+1]. 123 * leftknot2 should be set equal to zero at the first call 124 * to the rutine. 125 * 126 * OUTPUT : 127 * k1 - Max. principal curvature. 128 * k2 - Min. principal curvature. 129 * d1 - Max. direction of the principal curvature k1, given 130 * in local coordiantes (with regard to Xu,Xv). 131 * Dimension = 2. 132 * d2 - Min. direction of the principal curvature k2, given 133 * in local coordiantes (with regard to Xu,Xv). 134 * Dimension = 2. 135 * jstat - Status messages 136 * = 2 : Surface is degenerate at the point, that is, 137 * the surface is not regular at this point. 138 * = 1 : Surface is close to degenerate at the point. 139 * Angle between tangents is less than the angular 140 * tolerance. 141 * = 0 : Ok. 142 * < 0 : Error. 143 * 144 * 145 * METHOD : The princpal curvatures -k1 and -k2 are eigenvalues of 146 * dN, thus it turn out that we have to solve a 147 * eigenvalue/eigenvector problem, see references. 148 * 149 * REFERENCES : Differential Geometry of Curves and Surfaces, 150 * (Manfredo P. Do Carmo, Prentice Hall, 151 * ISBN: 0-13-212589-7). 152 * 153 * Elementary Linear Algebra 5e 154 * (Howard Anton, Wiley, ISBN:0-471-84819-0) 155 *- 156 * CALLS : s1422() and s2543(). 157 * 158 * LIMITATIONS : 159 * (i) If the surface is degenrated (not regular) at the point 160 * (u,v), it makes now sence to speak about curvature. 161 * The routine return istat == 2. 162 * (ii) If the surface is close to degenrate, the resuts 163 * can be numerical unstable. The routine return 164 * istat == 1. 165 * (iii) The dimension of the space in which the surface lies must 166 * be 1,2 or 3. 167 * 168 * 169 * WRITTEN BY : Geir Westgaard, SINTEF, Oslo, Norway. Date: 1995-1 170 ****************************************************************************** 171 */ 172 { 173 int kwarn = 0; /* Local staus variable(warning). */ 174 int der = 0; /* (dummy) */ 175 double derive[18]; /* Array containing the computed derivatives. */ 176 double normal[3]; /* Array containing the computed normalvektor. */ 177 178 179 if (ider != 0) goto err178; 180 181 182 if (surf == SISL_NULL) goto err150; 183 else 184 { 185 /* Compute derivates and normal. */ 186 187 s1422(surf,2,iside1,iside2,parvalue,leftknot1,leftknot2,derive,normal, 188 jstat); 189 if (*jstat > 0) kwarn = *jstat; 190 191 if (*jstat < 0) /* Error in lower level routine. */ 192 { 193 goto error; 194 } 195 else if (*jstat != 2) /* The surface is not degenerate */ 196 { 197 s2543(surf, der, derive, normal, k1, k2, d1, d2, jstat); 198 199 if (*jstat < 0) 200 goto error; 201 } 202 else if (*jstat == 2) /* The surface is degenerated. */ 203 { 204 *k1 = 0.0; 205 *k2 = 0.0; 206 d1[0] = 1.0; 207 d1[1] = 0.0; 208 d2[0] = 0.0; 209 d2[1] = 1.0; 210 211 goto war002; 212 } 213 214 } 215 216 /* Successful computations */ 217 218 *jstat = kwarn; 219 goto out; 220 221 /* The surface is degenerated at (u,v) */ 222 war002: 223 *jstat = 2; 224 goto out; 225 226 /* Error. Input (surface) pointer is SISL_NULL. */ 227 err150: 228 *jstat = -150; 229 s6err("s2542", *jstat, 0); 230 goto out; 231 232 /* Illegal derivative requested. */ 233 err178: 234 *jstat = -178; 235 s6err("s2542",*jstat,0); 236 goto out; 237 /* Error in lower level routine. */ 238 239 error: 240 s6err("s2542",*jstat,0); 241 goto out; 242 243 244 out: 245 246 return; 247 248 } 249 250 251