1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <stdint.h>
9 #include "curvature.h"
10
11 #define tolerance 0.1e-20
12
Inv4x4ColumnMajor(double m[16],double inv[16],double * det)13 static bool Inv4x4ColumnMajor(double m[16], double inv[16], double *det)
14 {
15 int i;
16
17 inv[0] = m[5] * m[10] * m[15] - m[5] * m[11] * m[14] - m[9] * m[6] * m[15] +
18 m[9] * m[7] * m[14] + m[13] * m[6] * m[11] - m[13] * m[7] * m[10];
19 inv[4] = -m[4] * m[10] * m[15] + m[4] * m[11] * m[14] + m[8] * m[6] * m[15] -
20 m[8] * m[7] * m[14] - m[12] * m[6] * m[11] + m[12] * m[7] * m[10];
21 inv[8] = m[4] * m[9] * m[15] - m[4] * m[11] * m[13] - m[8] * m[5] * m[15] +
22 m[8] * m[7] * m[13] + m[12] * m[5] * m[11] - m[12] * m[7] * m[9];
23 inv[12] = -m[4] * m[9] * m[14] + m[4] * m[10] * m[13] + m[8] * m[5] * m[14] -
24 m[8] * m[6] * m[13] - m[12] * m[5] * m[10] + m[12] * m[6] * m[9];
25 inv[1] = -m[1] * m[10] * m[15] + m[1] * m[11] * m[14] + m[9] * m[2] * m[15] -
26 m[9] * m[3] * m[14] - m[13] * m[2] * m[11] + m[13] * m[3] * m[10];
27 inv[5] = m[0] * m[10] * m[15] - m[0] * m[11] * m[14] - m[8] * m[2] * m[15] +
28 m[8] * m[3] * m[14] + m[12] * m[2] * m[11] - m[12] * m[3] * m[10];
29 inv[9] = -m[0] * m[9] * m[15] + m[0] * m[11] * m[13] + m[8] * m[1] * m[15] -
30 m[8] * m[3] * m[13] - m[12] * m[1] * m[11] + m[12] * m[3] * m[9];
31 inv[13] = m[0] * m[9] * m[14] - m[0] * m[10] * m[13] - m[8] * m[1] * m[14] +
32 m[8] * m[2] * m[13] + m[12] * m[1] * m[10] - m[12] * m[2] * m[9];
33 inv[2] = m[1] * m[6] * m[15] - m[1] * m[7] * m[14] - m[5] * m[2] * m[15] +
34 m[5] * m[3] * m[14] + m[13] * m[2] * m[7] - m[13] * m[3] * m[6];
35 inv[6] = -m[0] * m[6] * m[15] + m[0] * m[7] * m[14] + m[4] * m[2] * m[15] -
36 m[4] * m[3] * m[14] - m[12] * m[2] * m[7] + m[12] * m[3] * m[6];
37 inv[10] = m[0] * m[5] * m[15] - m[0] * m[7] * m[13] - m[4] * m[1] * m[15] +
38 m[4] * m[3] * m[13] + m[12] * m[1] * m[7] - m[12] * m[3] * m[5];
39 inv[14] = -m[0] * m[5] * m[14] + m[0] * m[6] * m[13] + m[4] * m[1] * m[14] -
40 m[4] * m[2] * m[13] - m[12] * m[1] * m[6] + m[12] * m[2] * m[5];
41 inv[3] = -m[1] * m[6] * m[11] + m[1] * m[7] * m[10] + m[5] * m[2] * m[11] -
42 m[5] * m[3] * m[10] - m[9] * m[2] * m[7] + m[9] * m[3] * m[6];
43 inv[7] = m[0] * m[6] * m[11] - m[0] * m[7] * m[10] - m[4] * m[2] * m[11] +
44 m[4] * m[3] * m[10] + m[8] * m[2] * m[7] - m[8] * m[3] * m[6];
45 inv[11] = -m[0] * m[5] * m[11] + m[0] * m[7] * m[9] + m[4] * m[1] * m[11] -
46 m[4] * m[3] * m[9] - m[8] * m[1] * m[7] + m[8] * m[3] * m[5];
47 inv[15] = m[0] * m[5] * m[10] - m[0] * m[6] * m[9] - m[4] * m[1] * m[10] +
48 m[4] * m[2] * m[9] + m[8] * m[1] * m[6] - m[8] * m[2] * m[5];
49
50 *det = m[0] * inv[0] + m[1] * inv[4] + m[2] * inv[8] + m[3] * inv[12];
51
52 if(*det == 0) return false;
53
54 double invDet = 1.0 / *det;
55
56 for(i = 0; i < 16; i++) inv[i] *= invDet;
57
58 return true;
59 }
60
solveEig(double A,double B,double C,double D,double * lambda1,double * v1x,double * v1y,double * lambda2,double * v2x,double * v2y)61 static void solveEig(double A, double B, double C, double D, double *lambda1,
62 double *v1x, double *v1y, double *lambda2, double *v2x,
63 double *v2y)
64 {
65 if(B * C <= tolerance) {
66 *lambda1 = A;
67 *v1x = 1;
68 *v1y = 0;
69 *lambda2 = D;
70 *v2x = 0;
71 *v2y = 1;
72 return;
73 }
74
75 double tr = A + D;
76 double det = A * D - B * C;
77 double S = sqrt((tr * tr / 4) - det);
78 *lambda1 = tr / 2 + S;
79 *lambda2 = tr / 2 - S;
80
81 double temp = ((A - D) * (A - D) / 4) + B * C;
82
83 double SS = sqrt(temp > 0 ? temp : 0.0);
84 if(A - D < 0) {
85 *v1x = C;
86 *v1y = -(A - D) / 2 + SS;
87 *v2x = +(A - D) / 2 - SS;
88 *v2y = B;
89 }
90 else {
91 *v2x = C;
92 *v2y = -(A - D) / 2 - SS;
93 *v1x = +(A - D) / 2 + SS;
94 *v1y = B;
95 }
96
97 double n1 = sqrt((*v1x) * (*v1x) + (*v1y) * (*v1y));
98 *v1x /= n1;
99 *v1y /= n1;
100 double n2 = sqrt((*v2x) * (*v2x) + (*v2y) * (*v2y));
101 *v2x /= n2;
102 *v2y /= n2;
103 }
104
node2trianglescmp(const void * p0,const void * p1)105 static inline int node2trianglescmp(const void *p0, const void *p1)
106 {
107 const std::size_t *e0 = (const std::size_t *)p0;
108 const std::size_t *e1 = (const std::size_t *)p1;
109
110 if(e0[0] < e1[0]) return -1;
111 if(e0[0] > e1[0]) return 1;
112 return 0;
113 }
114
normalize(double * n)115 static inline double normalize(double *n)
116 {
117 double d = sqrt(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
118 if(d != 0.0) {
119 n[0] /= d;
120 n[1] /= d;
121 n[2] /= d;
122 }
123 return fabs(d);
124 }
125
makevector(const double * a,const double * b,double * ba)126 static inline void makevector(const double *a, const double *b, double *ba)
127 {
128 ba[0] = b[0] - a[0];
129 ba[1] = b[1] - a[1];
130 ba[2] = b[2] - a[2];
131 }
132
dotprod(double * a,double * b)133 static inline double dotprod(double *a, double *b)
134 {
135 return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
136 }
137
crossprod(double * a,double * b,double * n)138 static inline void crossprod(double *a, double *b, double *n)
139 {
140 n[0] = a[1] * b[2] - a[2] * b[1];
141 n[1] = -(a[0] * b[2] - a[2] * b[0]);
142 n[2] = a[0] * b[1] - a[1] * b[0];
143 }
144
unitNormal2Triangle(const double * x1,const double * x2,const double * x3,double * n,double * s)145 static inline void unitNormal2Triangle(const double *x1, const double *x2,
146 const double *x3, double *n, double *s)
147 {
148 double a[3] = {x2[0] - x1[0], x2[1] - x1[1], x2[2] - x1[2]};
149 double b[3] = {x3[0] - x1[0], x3[1] - x1[1], x3[2] - x1[2]};
150
151 crossprod(a, b, n);
152 *s = 2 * normalize(n);
153 }
154
computeLocalFrame(double * n,double * u,double * v)155 static void computeLocalFrame(double *n, double *u, double *v)
156 {
157 if(fabs(n[0]) > fabs(n[1]) && fabs(n[0]) > fabs(n[2])) {
158 u[0] = 0;
159 u[1] = 0;
160 u[2] = 1;
161 }
162 else {
163 u[0] = 1;
164 u[1] = 0;
165 u[2] = 0;
166 }
167 crossprod(n, u, v);
168 normalize(v);
169 crossprod(v, n, u);
170 }
171
CurvatureRusinkiewicz(const std::vector<int> & triangles,const std::vector<SPoint3> & nodes,std::vector<std::pair<SVector3,SVector3>> & nodalCurvatures)172 bool CurvatureRusinkiewicz(
173 const std::vector<int> &triangles, const std::vector<SPoint3> &nodes,
174 std::vector<std::pair<SVector3, SVector3> > &nodalCurvatures)
175 {
176 std::size_t nTriangles = triangles.size() / 3;
177 std::size_t nVertices = nodes.size();
178
179 nodalCurvatures.resize(nVertices);
180
181 // for(int i = 0; i < nodes.size(); ++i){
182 // nodalCurvatures[i] = std::make_pair(0.0,0.0);
183 // }
184
185 std::vector<std::size_t> node2tri(3 * 2 * nTriangles);
186
187 // first compute node normals and node-to-triangle connectivity
188
189 std::vector<double> nodeNormals(3 * nVertices, 0.);
190 std::size_t counter = 0;
191 double n[3], surf;
192
193 for(std::size_t i = 0; i < nTriangles; i++) {
194 node2tri[counter++] = triangles[3 * i + 0];
195 node2tri[counter++] = i;
196 node2tri[counter++] = triangles[3 * i + 1];
197 node2tri[counter++] = i;
198 node2tri[counter++] = triangles[3 * i + 2];
199 node2tri[counter++] = i;
200 unitNormal2Triangle(nodes[triangles[3 * i + 0]].data(),
201 nodes[triangles[3 * i + 1]].data(),
202 nodes[triangles[3 * i + 2]].data(), n, &surf);
203 double *n0 = &nodeNormals[3 * triangles[3 * i + 0]];
204 double *n1 = &nodeNormals[3 * triangles[3 * i + 1]];
205 double *n2 = &nodeNormals[3 * triangles[3 * i + 2]];
206 for(std::size_t i1 = 0; i1 < 3; i1++) {
207 n0[i1] += n[i1];
208 n1[i1] += n[i1];
209 n2[i1] += n[i1];
210 }
211 }
212 for(std::size_t i = 0; i < nVertices; i++) normalize(&nodeNormals[3 * i]);
213
214 qsort(&node2tri[0], 3 * nTriangles, 2 * sizeof(std::size_t),
215 node2trianglescmp);
216
217 // compute the second fundamental tensor on each triangle using least squares
218
219 std::vector<double> CURV(4 * nTriangles);
220
221 for(std::size_t i = 0; i < nTriangles; i++) {
222 double *n0 = &nodeNormals[3 * triangles[3 * i + 0]];
223 double *n1 = &nodeNormals[3 * triangles[3 * i + 1]];
224 double *n2 = &nodeNormals[3 * triangles[3 * i + 2]];
225
226 double e0[3], e1[3], e2[3];
227 makevector(nodes[triangles[3 * i + 2]].data(),
228 nodes[triangles[3 * i + 1]].data(), e0);
229 makevector(nodes[triangles[3 * i + 0]].data(),
230 nodes[triangles[3 * i + 2]].data(), e1);
231 makevector(nodes[triangles[3 * i + 1]].data(),
232 nodes[triangles[3 * i + 0]].data(), e2);
233 unitNormal2Triangle(nodes[triangles[3 * i + 0]].data(),
234 nodes[triangles[3 * i + 1]].data(),
235 nodes[triangles[3 * i + 2]].data(), n, &surf);
236 double u[3], v[3];
237 makevector(nodes[triangles[3 * i + 0]].data(),
238 nodes[triangles[3 * i + 1]].data(), u);
239 normalize(u);
240 crossprod(n, u, v);
241
242 double sys[6][4], rhs[6], temp[3], invA[16], A[16], B[4];
243 sys[0][0] = sys[1][2] = dotprod(e0, u);
244 sys[0][1] = sys[1][3] = dotprod(e0, v);
245 sys[0][2] = sys[0][3] = sys[1][0] = sys[1][1] = 0;
246 makevector(n2, n1, temp);
247 rhs[0] = dotprod(temp, u);
248 rhs[1] = dotprod(temp, v);
249
250 sys[2][0] = sys[3][2] = dotprod(e1, u);
251 sys[2][1] = sys[3][3] = dotprod(e1, v);
252 sys[2][2] = sys[2][3] = sys[3][0] = sys[3][1] = 0;
253 makevector(n0, n2, temp);
254 rhs[2] = dotprod(temp, u);
255 rhs[3] = dotprod(temp, v);
256
257 sys[4][0] = sys[5][2] = dotprod(e2, u);
258 sys[4][1] = sys[5][3] = dotprod(e2, v);
259 sys[4][2] = sys[4][3] = sys[5][0] = sys[5][1] = 0;
260 makevector(n1, n0, temp);
261 rhs[4] = dotprod(temp, u);
262 rhs[5] = dotprod(temp, v);
263
264 for(std::size_t i1 = 0; i1 < 4; i1++) {
265 B[i1] = 0.0;
266 for(std::size_t i3 = 0; i3 < 6; i3++) { B[i1] += sys[i3][i1] * rhs[i3]; }
267 for(std::size_t i2 = 0; i2 < 4; i2++) {
268 A[i1 + 4 * i2] = 0.0;
269 for(std::size_t i3 = 0; i3 < 6; i3++) {
270 A[i1 + 4 * i2] += sys[i3][i2] * sys[i3][i1];
271 }
272 }
273 }
274 double det;
275 Inv4x4ColumnMajor(A, invA, &det);
276 for(std::size_t i1 = 0; i1 < 4; i1++) {
277 CURV[4 * i + i1] = 0.0;
278 for(std::size_t i2 = 0; i2 < 4; i2++) {
279 CURV[4 * i + i1] += invA[i1 + 4 * i2] * B[i2];
280 }
281 }
282 CURV[4 * i + 1] = CURV[4 * i + 2] =
283 0.5 * (CURV[4 * i + 1] + CURV[4 * i + 2]);
284 }
285
286 uint64_t count = 0;
287 double uP[3] = {0., 0., 0.}, vP[3] = {0., 0., 0.};
288 double A = 0., B = 0., D = 0.;
289 uint64_t iTriangle;
290 uint64_t ind = 0;
291 for(uint64_t iVert = 0; iVert < nVertices; ++iVert) {
292 count = 0;
293 A = 0.0;
294 B = 0.0;
295 D = 0.0;
296 while(node2tri[2 * ind + 0] == iVert) {
297 iTriangle = node2tri[2 * ind + 1];
298 computeLocalFrame(&nodeNormals[3 * iVert], uP, vP);
299 // computing each curvature around a vertex
300 unitNormal2Triangle(nodes[triangles[3 * iTriangle + 0]].data(),
301 nodes[triangles[3 * iTriangle + 1]].data(),
302 nodes[triangles[3 * iTriangle + 2]].data(), n, &surf);
303 double uF[3], vF[3];
304 makevector(nodes[triangles[3 * iTriangle + 0]].data(),
305 nodes[triangles[3 * iTriangle + 1]].data(), uF);
306 normalize(uF);
307 crossprod(n, uF, vF);
308 double *c = &CURV[4 * iTriangle];
309 double UP[3] = {dotprod(uP, uF), dotprod(uP, vF), 0};
310 normalize(UP);
311 double VP[3] = {dotprod(vP, uF), dotprod(vP, vF), 0};
312 normalize(VP);
313 A += (UP[0] * UP[0] * c[0] + 2 * UP[0] * UP[1] * c[1] +
314 UP[1] * UP[1] * c[3]);
315 D += (VP[0] * VP[0] * c[0] + 2 * VP[0] * VP[1] * c[1] +
316 VP[1] * VP[1] * c[3]);
317 B += (VP[0] * UP[0] * c[0] + (VP[1] * UP[0] + VP[0] * UP[1]) * c[1] +
318 VP[1] * UP[1] * c[3]);
319 ++count;
320 ++ind;
321 if(ind >= 3 * nTriangles) break;
322 }
323
324 A /= (double)count;
325 B /= (double)count;
326 D /= (double)count;
327 double lambda1, lambda2, v1x, v1y, v2x, v2y;
328 solveEig(A, B, B, D, &lambda1, &v1x, &v1y, &lambda2, &v2x, &v2y);
329 SVector3 cMax(fabs(lambda1) * (v1x * uP[0] + v1y * vP[0]),
330 fabs(lambda1) * (v1x * uP[1] + v1y * vP[1]),
331 fabs(lambda1) * (v1x * uP[2] + v1y * vP[2]));
332 SVector3 cMin(fabs(lambda2) * (v2x * uP[0] + v2y * vP[0]),
333 fabs(lambda2) * (v2x * uP[1] + v2y * vP[1]),
334 fabs(lambda2) * (v2x * uP[2] + v2y * vP[2]));
335 nodalCurvatures[iVert] = std::make_pair(cMax, cMin);
336 } // for iVert
337
338 return true;
339 }
340
CurvatureRusinkiewicz(const std::vector<int> & triangles,const std::vector<SPoint3> & nodes,std::vector<std::pair<SVector3,SVector3>> & nodalCurvatures,std::vector<double> & nodeNormals)341 bool CurvatureRusinkiewicz(
342 const std::vector<int> &triangles, const std::vector<SPoint3> &nodes,
343 std::vector<std::pair<SVector3, SVector3> > &nodalCurvatures,
344 std::vector<double> &nodeNormals)
345 {
346 std::size_t nTriangles = triangles.size() / 3;
347 std::size_t nVertices = nodes.size();
348
349 nodalCurvatures.resize(nVertices);
350
351 // for(int i = 0; i < nodes.size(); ++i){
352 // nodalCurvatures[i] = std::make_pair(0.0,0.0);
353 // }
354
355 std::vector<std::size_t> node2tri(3 * 2 * nTriangles);
356
357 // first compute node normals and node-to-triangle connectivity
358
359 // std::vector<double> nodeNormals(3 * nVertices, 0.);
360 std::size_t counter = 0;
361 double n[3], surf;
362
363 for(std::size_t i = 0; i < nTriangles; i++) {
364 node2tri[counter++] = triangles[3 * i + 0];
365 node2tri[counter++] = i;
366 node2tri[counter++] = triangles[3 * i + 1];
367 node2tri[counter++] = i;
368 node2tri[counter++] = triangles[3 * i + 2];
369 node2tri[counter++] = i;
370 unitNormal2Triangle(nodes[triangles[3 * i + 0]].data(),
371 nodes[triangles[3 * i + 1]].data(),
372 nodes[triangles[3 * i + 2]].data(), n, &surf);
373 double *n0 = &nodeNormals[3 * triangles[3 * i + 0]];
374 double *n1 = &nodeNormals[3 * triangles[3 * i + 1]];
375 double *n2 = &nodeNormals[3 * triangles[3 * i + 2]];
376 for(std::size_t i1 = 0; i1 < 3; i1++) {
377 n0[i1] += n[i1];
378 n1[i1] += n[i1];
379 n2[i1] += n[i1];
380 }
381 }
382 for(std::size_t i = 0; i < nVertices; i++) normalize(&nodeNormals[3 * i]);
383
384 qsort(&node2tri[0], 3 * nTriangles, 2 * sizeof(std::size_t),
385 node2trianglescmp);
386
387 // compute the second fundamental tensor on each triangle using least squares
388
389 std::vector<double> CURV(4 * nTriangles);
390
391 for(std::size_t i = 0; i < nTriangles; i++) {
392 double *n0 = &nodeNormals[3 * triangles[3 * i + 0]];
393 double *n1 = &nodeNormals[3 * triangles[3 * i + 1]];
394 double *n2 = &nodeNormals[3 * triangles[3 * i + 2]];
395
396 double e0[3], e1[3], e2[3];
397 makevector(nodes[triangles[3 * i + 2]].data(),
398 nodes[triangles[3 * i + 1]].data(), e0);
399 makevector(nodes[triangles[3 * i + 0]].data(),
400 nodes[triangles[3 * i + 2]].data(), e1);
401 makevector(nodes[triangles[3 * i + 1]].data(),
402 nodes[triangles[3 * i + 0]].data(), e2);
403 unitNormal2Triangle(nodes[triangles[3 * i + 0]].data(),
404 nodes[triangles[3 * i + 1]].data(),
405 nodes[triangles[3 * i + 2]].data(), n, &surf);
406 double u[3], v[3];
407 makevector(nodes[triangles[3 * i + 0]].data(),
408 nodes[triangles[3 * i + 1]].data(), u);
409 normalize(u);
410 crossprod(n, u, v);
411
412 double sys[6][4], rhs[6], temp[3], invA[16], A[16], B[4];
413 sys[0][0] = sys[1][2] = dotprod(e0, u);
414 sys[0][1] = sys[1][3] = dotprod(e0, v);
415 sys[0][2] = sys[0][3] = sys[1][0] = sys[1][1] = 0;
416 makevector(n2, n1, temp);
417 rhs[0] = dotprod(temp, u);
418 rhs[1] = dotprod(temp, v);
419
420 sys[2][0] = sys[3][2] = dotprod(e1, u);
421 sys[2][1] = sys[3][3] = dotprod(e1, v);
422 sys[2][2] = sys[2][3] = sys[3][0] = sys[3][1] = 0;
423 makevector(n0, n2, temp);
424 rhs[2] = dotprod(temp, u);
425 rhs[3] = dotprod(temp, v);
426
427 sys[4][0] = sys[5][2] = dotprod(e2, u);
428 sys[4][1] = sys[5][3] = dotprod(e2, v);
429 sys[4][2] = sys[4][3] = sys[5][0] = sys[5][1] = 0;
430 makevector(n1, n0, temp);
431 rhs[4] = dotprod(temp, u);
432 rhs[5] = dotprod(temp, v);
433
434 for(std::size_t i1 = 0; i1 < 4; i1++) {
435 B[i1] = 0.0;
436 for(std::size_t i3 = 0; i3 < 6; i3++) { B[i1] += sys[i3][i1] * rhs[i3]; }
437 for(std::size_t i2 = 0; i2 < 4; i2++) {
438 A[i1 + 4 * i2] = 0.0;
439 for(std::size_t i3 = 0; i3 < 6; i3++) {
440 A[i1 + 4 * i2] += sys[i3][i2] * sys[i3][i1];
441 }
442 }
443 }
444 double det;
445 Inv4x4ColumnMajor(A, invA, &det);
446 for(std::size_t i1 = 0; i1 < 4; i1++) {
447 CURV[4 * i + i1] = 0.0;
448 for(std::size_t i2 = 0; i2 < 4; i2++) {
449 CURV[4 * i + i1] += invA[i1 + 4 * i2] * B[i2];
450 }
451 }
452 CURV[4 * i + 1] = CURV[4 * i + 2] =
453 0.5 * (CURV[4 * i + 1] + CURV[4 * i + 2]);
454 }
455
456 uint64_t count = 0;
457 double uP[3] = {0., 0., 0.}, vP[3] = {0., 0., 0.};
458 double A = 0., B = 0., D = 0.;
459 uint64_t iTriangle;
460 uint64_t ind = 0;
461 for(uint64_t iVert = 0; iVert < nVertices; ++iVert) {
462 count = 0;
463 A = 0.0;
464 B = 0.0;
465 D = 0.0;
466 while(node2tri[2 * ind + 0] == iVert) {
467 iTriangle = node2tri[2 * ind + 1];
468 computeLocalFrame(&nodeNormals[3 * iVert], uP, vP);
469 // computing each curvature around a vertex
470 unitNormal2Triangle(nodes[triangles[3 * iTriangle + 0]].data(),
471 nodes[triangles[3 * iTriangle + 1]].data(),
472 nodes[triangles[3 * iTriangle + 2]].data(), n, &surf);
473 double uF[3], vF[3];
474 makevector(nodes[triangles[3 * iTriangle + 0]].data(),
475 nodes[triangles[3 * iTriangle + 1]].data(), uF);
476 normalize(uF);
477 crossprod(n, uF, vF);
478 double *c = &CURV[4 * iTriangle];
479 double UP[3] = {dotprod(uP, uF), dotprod(uP, vF), 0};
480 normalize(UP);
481 double VP[3] = {dotprod(vP, uF), dotprod(vP, vF), 0};
482 normalize(VP);
483 A += (UP[0] * UP[0] * c[0] + 2 * UP[0] * UP[1] * c[1] +
484 UP[1] * UP[1] * c[3]);
485 D += (VP[0] * VP[0] * c[0] + 2 * VP[0] * VP[1] * c[1] +
486 VP[1] * VP[1] * c[3]);
487 B += (VP[0] * UP[0] * c[0] + (VP[1] * UP[0] + VP[0] * UP[1]) * c[1] +
488 VP[1] * UP[1] * c[3]);
489 ++count;
490 ++ind;
491 if(ind >= 3 * nTriangles) break;
492 }
493
494 A /= (double)count;
495 B /= (double)count;
496 D /= (double)count;
497 double lambda1, lambda2, v1x, v1y, v2x, v2y;
498 solveEig(A, B, B, D, &lambda1, &v1x, &v1y, &lambda2, &v2x, &v2y);
499 SVector3 cMax(fabs(lambda1) * (v1x * uP[0] + v1y * vP[0]),
500 fabs(lambda1) * (v1x * uP[1] + v1y * vP[1]),
501 fabs(lambda1) * (v1x * uP[2] + v1y * vP[2]));
502 SVector3 cMin(fabs(lambda2) * (v2x * uP[0] + v2y * vP[0]),
503 fabs(lambda2) * (v2x * uP[1] + v2y * vP[1]),
504 fabs(lambda2) * (v2x * uP[2] + v2y * vP[2]));
505 nodalCurvatures[iVert] = std::make_pair(cMax, cMin);
506 } // for iVert
507
508 return true;
509 }