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: s2543.c,v 1.7 1996-08-02 07:29:05 jka Exp $
45 *
46 */
47
48
49 #define S2543
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s2543(SISLSurf * surf,int ider,double derive[],double normal[],double * k1,double * k2,double d1[],double d2[],int * jstat)55 s2543(SISLSurf *surf, int ider, double derive[], double normal[], double *k1,
56 double *k2, double d1[], double d2[], int *jstat)
57 #else
58 void s2543(surf, ider, derive, normal, k1, k2, d1, d2, jstat)
59 SISLSurf *surf;
60 int ider;
61 double derive[];
62 double normal[];
63 double *k1;
64 double *k2;
65 double d1[];
66 double d2[];
67 int *jstat;
68 #endif
69 /*
70 ***************************************************************************
71 *
72 ***************************************************************************
73 * PURPOSE : To compute principal curvature (k1,k2) with corresponding
74 * principal directions (d1,d2) of a surface for
75 * given values (u,v). This is a lower level routine,
76 * used for evaluation of many T(u,v)'s.
77 *
78 * INPUT :
79 * surf - Pointer to the surface to evaluate.
80 * ider - Number of derivatives to calculate.
81 * Only implemented for ider=0.
82 * < 0 : No derivative calculated.
83 * = 0 : Principal curvature calculated.
84 * = 1 : Principal curvature and its first derivative
85 * calculated.
86 * derive - Array containing derivatives from routine s1421().
87 * Size = idim*6.
88 * normal - Array containing the normal from routine s1421().
89 * Size = 3.
90 *
91 * INPUT/OUTPUT :
92 *
93 * OUTPUT :
94 * k1 - Max. principal curvature.
95 * k2 - Min. principal curvature.
96 * d1 - Max. direction of the principal curvature k1, given
97 * in local coordiantes (with regard to Xu,Xv).
98 * Dimension = 2.
99 * d2 - Min. direction of the principal curvature k2, given
100 * in local coordiantes (with regard to Xu,Xv).
101 * Dimension = 2.
102 * jstat - Status messages
103 * = 0 : Ok.
104 * < 0 : error.
105 *
106 * METHOD : The princpal curvatures -k1 and -k2 are eigenvalues of
107 * dN, thus it turn out that we have to solve a
108 * eigenvalue/eigenvector problem, see references.
109 *
110 * REFERENCES : Differential Geometry of Curves and Surfaces,
111 * (Manfredo P. Do Carmo, Prentice Hall,
112 * ISBN: 0-13-212589-7).
113 *
114 * Elementary Linear Algebra 5e
115 * (Howard Anton, Wiley, ISBN:0-471-84819-0)
116 *-
117 * CALLS :
118 *
119 * LIMITATIONS :
120 * (i) If the surface is degenrated (not regular) at the point
121 * (u,v), it makes now sence to speak about curvature.
122 * (ii) If the surface is closed to degenrate, the resuts
123 * can be numerical unstable.
124 * (iv) The dimension of the space in which the surface lies must
125 * be 1,2 or 3. The routine return istat < 0.
126 *
127 *
128 * WRITTEN BY : Geir Westgaard, SINTEF, Oslo, Norway. Date: 1995-1
129 * REWRITTEN BY : Johannes Kaasa, SINTEF, Oslo, Norway. Date: 1995-9
130 *****************************************************************************
131 */
132 {
133 double denom; /* Denominator in a fraction. */
134 double a, b, c; /* Coefficients of second degree equation. */
135 double sqrt_arg; /* Square argument. */
136 double ratio; /* Ratio of principal direction. */
137 double length; /* Parameter length. */
138 double fundform[6]; /* The coefficients of the fundamental forms.
139 The sequence is: E, F, G, e, f, g. */
140 double transform[4]; /* Transformation matrix.
141 The sequence is a11, a12, a21, a22. */
142 double Su[3]; /* Tangent in first parameter direction. */
143 double Sv[3]; /* Tangent in second parameter direction. */
144
145 if (surf->idim == 1 || surf->idim == 3) /* 1D and 3D surface */
146 {
147
148 /* Set up the tangents. */
149
150 if (surf->idim == 1)
151 {
152 Su[0] = 1.;
153 Su[1] = 0.;
154 Su[2] = derive[1];
155 Sv[0] = 0.;
156 Sv[1] = 1.;
157 Sv[2] = derive[2];
158 }
159 else
160 {
161 Su[0] = derive[3];
162 Su[1] = derive[4];
163 Su[2] = derive[5];
164 Sv[0] = derive[6];
165 Sv[1] = derive[7];
166 Sv[2] = derive[8];
167 }
168
169 /* Calculate the fundamental forms. */
170
171 s2513(surf, ider, 2, 1, derive, normal, fundform, jstat);
172 if (*jstat < 0) goto error;
173
174 /* Calculate the transformation matrix. */
175
176 denom = fundform[0]*fundform[2] - fundform[1]*fundform[1];
177
178 transform[0] = (fundform[1]*fundform[4] - fundform[2]*fundform[3])/
179 denom;
180 transform[1] = (fundform[1]*fundform[5] - fundform[2]*fundform[4])/
181 denom;
182 transform[2] = (fundform[1]*fundform[3] - fundform[0]*fundform[4])/
183 denom;
184 transform[3] = (fundform[1]*fundform[4] - fundform[0]*fundform[5])/
185 denom;
186
187 /* Calculate the principal curvature. */
188
189 a = 1.;
190 b = transform[0] + transform[3];
191 c = transform[0]*transform[3] - transform[1]*transform[2];
192
193 sqrt_arg = b*b - 4.*a*c;
194 if (sqrt_arg < REL_PAR_RES)
195 goto war100;
196
197 *k1 = (- b + sqrt(sqrt_arg))/(2.*a);
198 *k2 = (- b - sqrt(sqrt_arg))/(2.*a);
199
200 /* Calculate the principal directions. */
201
202 /* Maximal curvature direction. */
203
204 if (fabs(transform[0] + *k1) < REL_PAR_RES &&
205 fabs(transform[1]) < REL_PAR_RES)
206 {
207
208 /* Parallel to the u direction. */
209
210 length = 1./sqrt(Su[0]*Su[0] + Su[1]*Su[1] + Su[2]*Su[2]);
211
212 d1[0] = length;
213 d1[1] = 0.;
214 }
215 else if (fabs(transform[3] + *k1) < REL_PAR_RES &&
216 fabs(transform[2]) < REL_PAR_RES)
217 {
218
219 /* Parallel to the v direction. */
220
221 length = 1./sqrt(Sv[0]*Sv[0] + Sv[1]*Sv[1] + Sv[2]*Sv[2]);
222
223 d1[0] = 0.;
224 d1[1] = length;
225 }
226 else if (fabs(transform[0] + *k1) < fabs(transform[1]))
227 {
228 ratio = (transform[0] + *k1)/transform[1];
229 length = 1./sqrt((Su[0] - ratio*Sv[0])*(Su[0] - ratio*Sv[0]) +
230 (Su[1] - ratio*Sv[1])*(Su[1] - ratio*Sv[1]) +
231 (Su[2] - ratio*Sv[2])*(Su[2] - ratio*Sv[2]));
232
233 d1[0] = length;
234 d1[1] = -ratio*length;
235 }
236 else
237 {
238 ratio = transform[1]/(transform[0] + *k1);
239 length = 1./sqrt((Sv[0] - ratio*Su[0])*(Sv[0] - ratio*Su[0]) +
240 (Sv[1] - ratio*Su[1])*(Sv[1] - ratio*Su[1]) +
241 (Sv[2] - ratio*Su[2])*(Sv[2] - ratio*Su[2]));
242
243 d1[0] = -ratio*length;
244 d1[1] = length;
245 }
246
247 /* Minimal curvature direction. */
248
249 if (fabs(transform[0] + *k2) < REL_PAR_RES &&
250 fabs(transform[1]) < REL_PAR_RES)
251 {
252
253 /* Parallel to the u direction. */
254
255 length = 1./sqrt(Su[0]*Su[0] + Su[1]*Su[1] + Su[2]*Su[2]);
256
257 d2[0] = length;
258 d2[1] = 0.;
259 }
260 else if (fabs(transform[3] + *k2) < REL_PAR_RES &&
261 fabs(transform[2]) < REL_PAR_RES)
262 {
263
264 /* Parallel to the v direction. */
265
266 length = 1./sqrt(Sv[0]*Sv[0] + Sv[1]*Sv[1] + Sv[2]*Sv[2]);
267
268 d2[0] = 0.;
269 d2[1] = length;
270 }
271 else if (fabs(transform[0] + *k2) < fabs(transform[1]))
272 {
273 ratio = (transform[0] + *k2)/transform[1];
274 length = 1./sqrt((Su[0] - ratio*Sv[0])*(Su[0] - ratio*Sv[0]) +
275 (Su[1] - ratio*Sv[1])*(Su[1] - ratio*Sv[1]) +
276 (Su[2] - ratio*Sv[2])*(Su[2] - ratio*Sv[2]));
277
278 d2[0] = length;
279 d2[1] = -ratio*length;
280 }
281 else
282 {
283 ratio = transform[1]/(transform[0] + *k2);
284 length = 1./sqrt((Sv[0] - ratio*Su[0])*(Sv[0] - ratio*Su[0]) +
285 (Sv[1] - ratio*Su[1])*(Sv[1] - ratio*Su[1]) +
286 (Sv[2] - ratio*Su[2])*(Sv[2] - ratio*Su[2]));
287
288 d2[0] = -ratio*length;
289 d2[1] = length;
290 }
291
292 }
293 else if (surf->idim == 2) /* 2D surface */
294 {
295 /* The surface lies in a plane => T(u,v) = 0 */
296
297 *k1 = 0.0;
298 *k2 = 0.0;
299 d1[0] = 1.0;
300 d1[1] = 0.0;
301 d2[0] = 0.0;
302 d2[1] = 1.0;
303 }
304 else /* When surf->idim != 1,2 or 3 */
305 {
306 goto err105;
307 }
308
309
310 /* Successful computations */
311
312 *jstat = 0;
313 goto out;
314
315
316 /* The surface does not have principal curvatures. */
317 war100:
318 if (fabs(sqrt_arg) < REL_PAR_RES)
319 {
320 *k1 = - b/(2.*a);
321 *k2 = *k1;
322 }
323 else
324 {
325 *k1 = 0.0;
326 *k2 = 0.0;
327 }
328 d1[0] = 1.0;
329 d1[1] = 0.0;
330 d2[0] = 0.0;
331 d2[1] = 1.0;
332 goto out;
333
334 /* Error in input, surf->idim != 1,2 or 3. */
335 err105:
336 *jstat = -105;
337 s6err("s2543",*jstat,0);
338 goto out;
339
340 error:
341 s6err("s2543",*jstat,0);
342 goto out;
343
344 out:
345
346 return;
347
348 }
349