GCC Code Coverage Report

Directory: ./
File: Geometry/curvatures.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 167 202 82.7%
Branches: 198 426 46.5%
Line Branch Exec Source
1 // ______ ______ _ _ _____ ______
2 // | ____| ____| | (_)/ ____| | ____|
3 // | |__ | |__ | | _| (___ ___| |__
4 // | __| | __| | | | |\___ \ / __| __|
5 // | | | |____| |____| |____) | (__| |____
6 // |_| |______|______|_|_____/ \___|______|
7 // Finite Elements for Life Sciences and Engineering
8 //
9 // License: LGL2.1 License
10 // FELiScE default license: LICENSE in root folder
11 //
12 // Main authors: Matteo
13 //
14
15 // System includes
16 #include <vector>
17
18 // External includes
19
20 // Project includes
21 #include "Geometry/curvatures.hpp"
22
23 namespace felisce
24 {
25
9/18
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 16 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 16 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 16 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 16 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 16 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 16 times.
✗ Branch 21 not taken.
✓ Branch 23 taken 16 times.
✗ Branch 24 not taken.
✓ Branch 26 taken 16 times.
✗ Branch 27 not taken.
 16 Curvatures::Curvatures() {
26
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 m_minProportion=FelisceParam::instance().fiberDispersionParameter;
27 16 }
28
29 215856 void Curvatures::update(const ElementField& normal, const UBlasMatrix& firstFundForm, const UBlasMatrix& covBasis) {
30
31 215856 m_firstFundForm = firstFundForm;
32
33 215856 computeSecondFundamentalForm(covBasis,normal);
34 215856 computeShapeOperator();
35 215856 computeCurvatures();
36 215856 computePrincipalDirections();
37 215856 computeArcLengths();
38 215856 computeProportions();
39 215856 computeProjectors();
40 215856 computeEigInCartCoordinate(covBasis,normal);
41
42 215856 }
43
44 215856 void Curvatures::computeEigInCartCoordinate(const UBlasMatrix& covBasis, const ElementField& normal ) {
45
1/2
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
 215856 m_eigMax.resize(3);
46
1/2
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
 215856 m_eigMin.resize(3);
47
2/2
✓ Branch 0 taken 647568 times.
✓ Branch 1 taken 215856 times.
 863424 for (int c(0); c<3; c++) {
48
5/10
✓ Branch 1 taken 647568 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 647568 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 647568 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 647568 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 647568 times.
✗ Branch 14 not taken.
 647568 m_eigMax(c) = m_maxEigVec(0)*covBasis(0,c) + m_maxEigVec(1)*covBasis(1,c);
49
5/10
✓ Branch 1 taken 647568 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 647568 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 647568 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 647568 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 647568 times.
✗ Branch 14 not taken.
 647568 m_eigMin(c) = m_minEigVec(0)*covBasis(0,c) + m_minEigVec(1)*covBasis(1,c);
50 }
51
3/6
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 215856 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 215856 times.
✗ Branch 8 not taken.
 215856 m_eigMax = m_eigMax/norm_2(m_eigMax);
52
3/6
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 215856 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 215856 times.
✗ Branch 8 not taken.
 215856 m_eigMin = m_eigMin/norm_2(m_eigMin);
53
54
1/2
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
 215856 UBlasVector ones(3);
55
56
1/2
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
 215856 ones(0) = 1;
57
1/2
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
 215856 ones(1) = 1;
58
1/2
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
 215856 ones(2) = 0;
59
3/4
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 108420 times.
✓ Branch 4 taken 107436 times.
 215856 if ( inner_prod ( m_eigMin, ones ) < 0 )
60
2/4
✓ Branch 1 taken 108420 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 108420 times.
✗ Branch 5 not taken.
 108420 m_eigMin = -1*m_eigMin;
61
62
1/2
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
 215856 UBlasVector n(3);
63
2/4
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 215856 times.
✗ Branch 5 not taken.
 215856 n(0)=normal.val(0,0);
64
2/4
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 215856 times.
✗ Branch 5 not taken.
 215856 n(1)=normal.val(1,0);
65
2/4
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 215856 times.
✗ Branch 5 not taken.
 215856 n(2)=normal.val(2,0);
66
67
1/2
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
 215856 UBlasVector crossProd(3);
68
5/10
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 215856 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 215856 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 215856 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 215856 times.
✗ Branch 14 not taken.
 215856 crossProd(0) = m_eigMax(1)*m_eigMin(2) - m_eigMax(2)*m_eigMin(1);
69
5/10
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 215856 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 215856 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 215856 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 215856 times.
✗ Branch 14 not taken.
 215856 crossProd(1) = m_eigMax(2)*m_eigMin(0) - m_eigMax(0)*m_eigMin(2);
70
5/10
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 215856 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 215856 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 215856 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 215856 times.
✗ Branch 14 not taken.
 215856 crossProd(2) = m_eigMax(0)*m_eigMin(1) - m_eigMax(1)*m_eigMin(0);
71
72
3/4
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 107628 times.
✓ Branch 4 taken 108228 times.
 215856 if ( inner_prod(crossProd, n) < 0 )
73
2/4
✓ Branch 1 taken 107628 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 107628 times.
✗ Branch 5 not taken.
 107628 m_eigMax = -1*m_eigMax;
74 215856 }
75
76 215856 void Curvatures::computeSecondFundamentalForm( const UBlasMatrix& covBasis, const ElementField& normal ) {
77
78 215856 felInt dimension = 3;
79 215856 std::vector<double> d10,d20;
80
1/2
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
 215856 d10.resize(dimension,0.0);
81
1/2
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
 215856 d20.resize(dimension,0.0);
82
83 // gradient of the normal
84
2/2
✓ Branch 0 taken 647568 times.
✓ Branch 1 taken 215856 times.
 863424 for ( int c(0.0); c<dimension; ++c) {
85
2/4
✓ Branch 1 taken 647568 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 647568 times.
✗ Branch 5 not taken.
 647568 d10[c] = normal.val(c,1) - normal.val(c,0);
86
2/4
✓ Branch 1 taken 647568 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 647568 times.
✗ Branch 5 not taken.
 647568 d20[c] = normal.val(c,2) - normal.val(c,0);
87 }
88
89 215856 double e(0),f1(0),f2(0),g(0);
90 // second fundamental form b_{\alpha,\beta} = - \partial_{\alpha}a_3 \dot a_{\beta}
91
2/2
✓ Branch 0 taken 647568 times.
✓ Branch 1 taken 215856 times.
 863424 for ( int c(0); c < dimension; ++c) {
92
1/2
✓ Branch 2 taken 647568 times.
✗ Branch 3 not taken.
 647568 e += -d10[c]*covBasis(0,c);
93
1/2
✓ Branch 2 taken 647568 times.
✗ Branch 3 not taken.
 647568 f1 += -d10[c]*covBasis(1,c);
94
1/2
✓ Branch 2 taken 647568 times.
✗ Branch 3 not taken.
 647568 f2 += -d20[c]*covBasis(0,c);
95
1/2
✓ Branch 2 taken 647568 times.
✗ Branch 3 not taken.
 647568 g += -d20[c]*covBasis(1,c);
96 }
97
98
1/2
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
 215856 m_secondFundForm.resize(2,2);
99 // numerically it is not symmetric -> we take its symmetric part
100
1/2
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
 215856 m_secondFundForm(0,0) = e;
101
1/2
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
 215856 m_secondFundForm(1,1) = g;
102
1/2
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
 215856 m_secondFundForm(0,1) = (f1 + f2)*0.5;
103
1/2
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
 215856 m_secondFundForm(1,0) = (f1 + f2)*0.5;
104 215856 }
105
106 215856 void Curvatures::computeShapeOperator() {
107
108 215856 m_invFFT.resize(2,2);
109
110 215856 m_invFFT(0,0) = m_firstFundForm(1,1);
111 215856 m_invFFT(1,1) = m_firstFundForm(0,0);
112
113 215856 m_invFFT(0,1) = - m_firstFundForm(0,1);
114 215856 m_invFFT(1,0) = - m_firstFundForm(1,0);
115
116 431712 m_invFFT = m_invFFT/( m_firstFundForm(0,0)*m_firstFundForm(1,1)
117
2/4
✓ Branch 3 taken 215856 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 215856 times.
✗ Branch 7 not taken.
 215856 -m_firstFundForm(0,1)*m_firstFundForm(1,0) );
118
119
1/2
✓ Branch 2 taken 215856 times.
✗ Branch 3 not taken.
 215856 m_shapeOperator = prod( m_invFFT, m_secondFundForm);
120
121 215856 }
122
123 215856 void Curvatures::computeCurvatures() {
124
125 215856 double tr(0.0),det(1.0);
126 215856 tr = ( m_shapeOperator(0,0) + m_shapeOperator(1,1) );
127 215856 det = m_shapeOperator(0,0)*m_shapeOperator(1,1) - m_shapeOperator(0,1)*m_shapeOperator(1,0);
128
129 // max and min here refer to the case of normal going inside.
130 // since the normal is going outside the domain they are swapped
131 215856 m_maxCurv = (tr - std::sqrt( tr*tr - 4*det ) )/2;
132 215856 m_minCurv = (tr + std::sqrt( tr*tr - 4*det ) )/2;
133
134 215856 m_meanCurv = tr/2;
135 215856 m_gaussCurv = det;
136
137 215856 }
138
139 ✗ void Curvatures::print() const {
140 ✗ std::cout<<"================CURVATURES================="<<std::endl;
141 ✗ std::cout<<" Kmax: "<<m_maxCurv<<std::endl;
142 ✗ std::cout<<" Kmin: "<<m_minCurv<<std::endl;
143 ✗ std::cout<<" GaussK: "<<m_gaussCurv<<std::endl;
144 ✗ std::cout<<" MeanK: "<<m_meanCurv<<std::endl;
145 ✗ std::cout<<" Principal direction ( max "<<m_maxEigVec<<std::endl;
146 ✗ std::cout<<" Principal direction ( min "<<m_minEigVec<<std::endl;
147 ✗ std::cout<<" Length Arc (lambda max) "<<m_rayleighMaxEigVec<<std::endl;
148 ✗ std::cout<<" Length Arc (lambda min) "<<m_rayleighMinEigVec<<std::endl;
149 ✗ std::cout<<" Proportions: max: "<<m_maxProp<<std::endl;
150 ✗ std::cout<<" Proportions: min: "<<m_minProp<<std::endl;
151 ✗ std::cout<<" Projectors: max: "<<m_projOnMaxEigVec<<std::endl;
152 ✗ std::cout<<" Projectors: min: "<<m_projOnMinEigVec<<std::endl;
153 ✗ std::cout<<"==========================================="<<std::endl;
154 }
155
156 215856 void Curvatures::computePrincipalDirections() {
157
158
2/4
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 215856 times.
✗ Branch 5 not taken.
 215856 UBlasMatrix Smin(2,2), Smax(2,2);
159
2/4
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 215856 times.
✗ Branch 5 not taken.
 215856 double toll( std::max( 1e-4 * norm_frobenius(m_shapeOperator), 1e-10 )); //it used to be std::set to 0.01
160
161
1/2
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
 215856 m_maxEigVec.resize(2);
162
1/2
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
 215856 m_minEigVec.resize(2);
163
164
4/8
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 215856 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 215856 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 215856 times.
✗ Branch 11 not taken.
 215856 Smin = m_shapeOperator - m_minCurv*( UBlasIdentityMatrix ( 2 ) );
165
4/8
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 215856 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 215856 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 215856 times.
✗ Branch 11 not taken.
 215856 Smax = m_shapeOperator - m_maxCurv*( UBlasIdentityMatrix ( 2 ) );
166
167
3/4
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 212427 times.
✓ Branch 4 taken 3429 times.
 215856 if ( std::fabs( Smin(0,0) ) > toll ) {
168
3/6
✓ Branch 1 taken 212427 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 212427 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 212427 times.
✗ Branch 8 not taken.
 212427 m_minEigVec(0) = -Smin(0,1)/Smin(0,0);
169
1/2
✓ Branch 1 taken 212427 times.
✗ Branch 2 not taken.
 212427 m_minEigVec(1) = 1;
170
3/4
✓ Branch 1 taken 3429 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 201 times.
✓ Branch 4 taken 3228 times.
 3429 } else if ( std::fabs( Smin(0,1) ) > toll ) {
171
1/2
✓ Branch 1 taken 201 times.
✗ Branch 2 not taken.
 201 m_minEigVec(0) = 1;
172
3/6
✓ Branch 1 taken 201 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 201 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 201 times.
✗ Branch 8 not taken.
 201 m_minEigVec(1) = -Smin(0,0)/Smin(0,1);
173
3/4
✓ Branch 1 taken 3228 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 27 times.
✓ Branch 4 taken 3201 times.
 3228 } else if ( std::fabs( Smin(1,1) ) > toll ) {
174
1/2
✓ Branch 1 taken 27 times.
✗ Branch 2 not taken.
 27 m_minEigVec(0) = 1;
175
3/6
✓ Branch 1 taken 27 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 27 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 27 times.
✗ Branch 8 not taken.
 27 m_minEigVec(1) = -Smin(1,0)/Smin(1,1);
176
2/4
✓ Branch 1 taken 3201 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 3201 times.
 3201 } else if ( std::fabs( Smin(1,0) ) > toll ) {
177 ✗ m_minEigVec(0) = -Smin(1,1)/Smin(1,0);
178 ✗ m_minEigVec(1) = 1;
179
1/2
✓ Branch 1 taken 3201 times.
✗ Branch 2 not taken.
 3201 } else if ( std::fabs( Smin(0,0) )
180
1/2
✓ Branch 1 taken 3201 times.
✗ Branch 2 not taken.
 3201 + std::fabs( Smin(1,0) )
181
1/2
✓ Branch 1 taken 3201 times.
✗ Branch 2 not taken.
 3201 + std::fabs( Smin(0,1) )
182
2/4
✓ Branch 1 taken 3201 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3201 times.
✗ Branch 4 not taken.
 3201 + std::fabs( Smin(1,1) ) < toll ) {
183
1/2
✓ Branch 1 taken 3201 times.
✗ Branch 2 not taken.
 3201 m_minEigVec(0) = 1;
184
1/2
✓ Branch 1 taken 3201 times.
✗ Branch 2 not taken.
 3201 m_minEigVec(1) = 0;
185 } else {
186 ✗ std::cout<<"toll : "<<toll<<std::endl;
187 ✗ std::cout<<"Shape operator: ";
188 ✗ std::cout<<m_shapeOperator<<std::endl;
189 ✗ std::cout<<"Shape operator - lambda I:";
190 ✗ std::cout<<Smin<<std::endl;
191 ✗ FEL_ERROR(" TODO in curvatures.cpp ");
192 }
193
194
3/6
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 215856 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 215856 times.
✗ Branch 8 not taken.
 215856 m_minEigVec = m_minEigVec/norm_2(m_minEigVec);
195
196
197
3/4
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 212469 times.
✓ Branch 4 taken 3387 times.
 215856 if ( std::fabs( Smax(0,0) ) > toll ) {
198
3/6
✓ Branch 1 taken 212469 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 212469 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 212469 times.
✗ Branch 8 not taken.
 212469 m_maxEigVec(0) = -Smax(0,1)/Smax(0,0);
199
1/2
✓ Branch 1 taken 212469 times.
✗ Branch 2 not taken.
 212469 m_maxEigVec(1) = 1;
200
3/4
✓ Branch 1 taken 3387 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 150 times.
✓ Branch 4 taken 3237 times.
 3387 } else if ( std::fabs( Smax(0,1) ) > toll ) {
201
1/2
✓ Branch 1 taken 150 times.
✗ Branch 2 not taken.
 150 m_maxEigVec(0) = 1;
202
3/6
✓ Branch 1 taken 150 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 150 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 150 times.
✗ Branch 8 not taken.
 150 m_maxEigVec(1) = -Smax(0,0)/Smax(0,1);
203
3/4
✓ Branch 1 taken 3237 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 36 times.
✓ Branch 4 taken 3201 times.
 3237 } else if ( std::fabs( Smax(1,1) ) > toll ) {
204
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 36 m_maxEigVec(0) = 1;
205
3/6
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 36 times.
✗ Branch 8 not taken.
 36 m_maxEigVec(1) = -Smax(1,0)/Smax(1,1);
206
2/4
✓ Branch 1 taken 3201 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 3201 times.
 3201 } else if ( std::fabs( Smax(1,0) ) > toll ) {
207 ✗ m_maxEigVec(0) = -Smax(1,1)/Smax(1,0);
208 ✗ m_maxEigVec(1) = 1;
209
1/2
✓ Branch 1 taken 3201 times.
✗ Branch 2 not taken.
 3201 } else if ( std::fabs( Smax(0,0) )
210
1/2
✓ Branch 1 taken 3201 times.
✗ Branch 2 not taken.
 3201 + std::fabs( Smax(1,0) )
211
1/2
✓ Branch 1 taken 3201 times.
✗ Branch 2 not taken.
 3201 + std::fabs( Smax(0,1) )
212
2/4
✓ Branch 1 taken 3201 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3201 times.
✗ Branch 4 not taken.
 3201 + std::fabs( Smax(1,1) ) < toll ) {
213
1/2
✓ Branch 1 taken 3201 times.
✗ Branch 2 not taken.
 3201 m_maxEigVec(0) = 0;
214
1/2
✓ Branch 1 taken 3201 times.
✗ Branch 2 not taken.
 3201 m_maxEigVec(1) = 1;
215 } else {
216 ✗ std::cout<<"toll : "<<toll<<std::endl;
217 ✗ std::cout<<"Shape operator: ";
218 ✗ std::cout<<m_shapeOperator<<std::endl;
219 ✗ std::cout<<"Shape operator - lambda I:";
220 ✗ std::cout<<Smin<<std::endl;
221 ✗ FEL_ERROR(" TODO in curvatures.cpp");
222 }
223
224
3/6
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 215856 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 215856 times.
✗ Branch 8 not taken.
 215856 m_maxEigVec = m_maxEigVec/norm_2(m_maxEigVec);
225
226 215856 }
227
228 215856 void Curvatures::computeArcLengths() {
229
230
1/2
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
 215856 m_rayleighMaxEigVec = std::sqrt( inner_prod( m_maxEigVec,
231 431712 prod( m_firstFundForm, m_maxEigVec )
232 ));
233
234
1/2
✓ Branch 1 taken 215856 times.
✗ Branch 2 not taken.
 215856 m_rayleighMinEigVec = std::sqrt( inner_prod( m_minEigVec,
235 431712 prod( m_firstFundForm, m_minEigVec )
236 ));
237 215856 }
238
239 215856 void Curvatures::computeProportions() {
240
241 215856 double c(0.5);
242
243 215856 double k1 = m_minCurv;
244 215856 double k2 = m_maxCurv;
245
246 215856 double toll = 1;
247
2/2
✓ Branch 0 taken 101595 times.
✓ Branch 1 taken 114261 times.
 215856 if( k1 * k2 < 0 ) { //k2 positive and k1 negative ( saddle point )
248 101595 c = 0;
249
2/2
✓ Branch 0 taken 3201 times.
✓ Branch 1 taken 111060 times.
 114261 } else if( k1 + k2 > - toll) {//both negative or both equal to zero
250 3201 c = 0.5;
251
1/2
✓ Branch 0 taken 111060 times.
✗ Branch 1 not taken.
 111060 } else if( k1 <= 0 ) {//both positive
252 111060 c = 0.5*k1/k2;
253 } else {
254 ✗ FEL_ERROR(" TODO in curvatures.cpp ");
255 }
256
257 215856 m_maxProp = (1-c)*(1-m_minProportion) + m_minProportion/2;
258 215856 m_minProp = (c)*(1-m_minProportion) + m_minProportion/2;
259 215856 }
260
261 215856 void Curvatures::computeProjectors() {
262
263
2/4
✓ Branch 2 taken 215856 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 215856 times.
✗ Branch 6 not taken.
 215856 m_projOnMaxEigVec = m_maxProp/(m_rayleighMaxEigVec*m_rayleighMaxEigVec)*outer_prod( m_maxEigVec, m_maxEigVec);
264
2/4
✓ Branch 2 taken 215856 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 215856 times.
✗ Branch 6 not taken.
 215856 m_projOnMinEigVec = m_minProp/(m_rayleighMinEigVec*m_rayleighMinEigVec)*outer_prod( m_minEigVec, m_minEigVec);
265
266 215856 }
267
268 20939787 UBlasVector Curvatures::testFuncGradient( int idLocDof ) {
269 20939787 UBlasVector testFuncGrad(2);
270
3/4
✓ Branch 0 taken 6979929 times.
✓ Branch 1 taken 6979929 times.
✓ Branch 2 taken 6979929 times.
✗ Branch 3 not taken.
 20939787 switch(idLocDof) {
271 6979929 case 0:
272
1/2
✓ Branch 1 taken 6979929 times.
✗ Branch 2 not taken.
 6979929 testFuncGrad(0) = -1;
273
1/2
✓ Branch 1 taken 6979929 times.
✗ Branch 2 not taken.
 6979929 testFuncGrad(1) = -1;
274 6979929 break;
275 6979929 case 1:
276
1/2
✓ Branch 1 taken 6979929 times.
✗ Branch 2 not taken.
 6979929 testFuncGrad(0) = 1;
277
1/2
✓ Branch 1 taken 6979929 times.
✗ Branch 2 not taken.
 6979929 testFuncGrad(1) = 0;
278 6979929 break;
279 6979929 case 2:
280
1/2
✓ Branch 1 taken 6979929 times.
✗ Branch 2 not taken.
 6979929 testFuncGrad(0) = 0;
281
1/2
✓ Branch 1 taken 6979929 times.
✗ Branch 2 not taken.
 6979929 testFuncGrad(1) = 1;
282 6979929 break;
283 ✗ default:
284 ✗ FEL_ERROR("you can not use testFuncGradient on elements with more than 3 vertices");
285 ✗ break;
286 }
287 20939787 return testFuncGrad;
288 }
289
290 2866107 UBlasVector Curvatures::elemFieldGradient( const ElementField& f ) {
291 2866107 UBlasVector fGrad(2);
292
3/6
✓ Branch 1 taken 2866107 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2866107 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2866107 times.
✗ Branch 8 not taken.
 2866107 fGrad(0) = f.val(0,1) - f.val(0,0);
293
3/6
✓ Branch 1 taken 2866107 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2866107 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2866107 times.
✗ Branch 8 not taken.
 2866107 fGrad(1) = f.val(0,2) - f.val(0,0);
294 2866107 return fGrad;
295 }
296
297 967185 double Curvatures::mNormOfGradient( const ElementField& f , const UBlasMatrix& M) {
298
3/6
✓ Branch 2 taken 967185 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 967185 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 967185 times.
✗ Branch 9 not taken.
 967185 return inner_prod( elemFieldGradient(f), prod(M,elemFieldGradient(f) ) );
299 }
300 }
301