GCC Code Coverage Report

Directory: ./
File: HyperelasticityLaws/ciarletGeymonat.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 79 83 95.2%
Branches: 9 26 34.6%
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: S. Gilles
13 //
14
15 // System includes
16 #include <string>
17 #include <limits>
18
19 // External includes
20
21 // Project includes
22 #include "Core/felisceTools.hpp"
23 #include "Core/felisceParam.hpp"
24 #include "HyperelasticityLaws/ciarletGeymonat.hpp"
25
26 namespace felisce
27 {
28 577 CiarletGeymonat::CiarletGeymonat()
29 : BaseClassType(),
30 577 m_kappa1(1.e20), // silly value on purpose
31 577 m_kappa2(1.e20)
32 {
33 577 }
34
35 /***********************************************************************************/
36 /***********************************************************************************/
37
38 25 void CiarletGeymonat::InitParameters()
39 {
40 // Call base class
41 25 BaseClassType::InitParameters();
42
43 25 const auto& r_kappas = FelisceParam::instance().CiarletGeymonat;
44 25 m_kappa1 = r_kappas.kappa1;
45 25 m_kappa2 = r_kappas.kappa2;
46
47 // Check the values were correctly initialized in the input file
48 {
49 25 std::bitset<3> are_parameters_unused;
50
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 25 times.
 25 if (m_kappa1 < 0.)
51 ✗ are_parameters_unused.set(0);
52
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 25 times.
 25 if (m_kappa2 < 0.)
53 ✗ are_parameters_unused.set(1);
54
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 25 times.
 25 if (m_bulk < 0.)
55 ✗ are_parameters_unused.set(2);
56
57
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 25 times.
 25 if (are_parameters_unused.any())
58 ✗ throw ErrorInvalidParameter(are_parameters_unused);
59 }
60 25 }
61
62 /***********************************************************************************/
63 /***********************************************************************************/
64
65 3 double CiarletGeymonat::W() const
66 {
67 3 BaseClassType::W();
68
69 3 const double I1 = HyperElasticLaw::m_invariants[0];
70 3 const double I2 = HyperElasticLaw::m_invariants[1];
71 3 const double I3 = HyperElasticLaw::m_invariants[2];
72 3 const double I3_pow_minus_one_third = std::pow(I3, -1.0/3.0);
73
74 3 return m_kappa1 * (I1 * I3_pow_minus_one_third - 3.0)
75 3 + m_kappa2 * (I2 * std::pow(I3_pow_minus_one_third, 2) - 3.0);
76 }
77
78 /***********************************************************************************/
79 /***********************************************************************************/
80
81 798000 double CiarletGeymonat::FirstDerivativeWFirstInvariant() const
82 {
83 798000 BaseClassType::FirstDerivativeWFirstInvariant();
84 798000 return m_kappa1 * std::pow(HyperElasticLaw::m_invariants[2], -1.0/3.0);
85 }
86
87 /***********************************************************************************/
88 /***********************************************************************************/
89
90 798003 double CiarletGeymonat::FirstDerivativeWSecondInvariant() const
91 {
92 798003 BaseClassType::FirstDerivativeWSecondInvariant();
93 798003 return m_kappa2 * std::pow(HyperElasticLaw::m_invariants[2], -2.0/3.0);
94 }
95
96 /***********************************************************************************/
97 /***********************************************************************************/
98
99 798002 double CiarletGeymonat::SecondDerivativeWFirstInvariant() const
100 {
101 798002 return 0.0;
102 }
103
104 /***********************************************************************************/
105 /***********************************************************************************/
106
107 798002 double CiarletGeymonat::SecondDerivativeWSecondInvariant() const
108 {
109 798002 return 0.0;
110 }
111
112 /***********************************************************************************/
113 /***********************************************************************************/
114
115 798002 double CiarletGeymonat::SecondDerivativeWFirstAndSecondInvariant() const
116 {
117 798002 return 0.0;
118 }
119
120 /***********************************************************************************/
121 /***********************************************************************************/
122
123 798003 double CiarletGeymonat::FirstDerivativeWThirdInvariant() const
124 {
125 798003 BaseClassType::FirstDerivativeWThirdInvariant();
126
127 798003 const double I1 = HyperElasticLaw::m_invariants[0];
128 798003 const double I2 = HyperElasticLaw::m_invariants[1];
129 798003 const double I3 = HyperElasticLaw::m_invariants[2];
130
131
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 798003 times.
 798003 FEL_ASSERT("I3 should be strictly positive! (and more precisely somewhere in the vicinity of 1.)" && std::abs(I3) > std::numeric_limits<double>::epsilon());
132
133 798003 const double minus_one_third = -1.0/3.0;
134
135 798003 return m_kappa1 * I1 * minus_one_third * std::pow(I3, 4. * minus_one_third)
136 798003 + m_kappa2 * I2 * 2. * minus_one_third * std::pow(I3, 5. * minus_one_third);
137 }
138
139 /***********************************************************************************/
140 /***********************************************************************************/
141
142 798003 double CiarletGeymonat::SecondDerivativeWThirdInvariant() const
143 {
144 798003 BaseClassType::SecondDerivativeWThirdInvariant();
145
146 798003 const double I1 = HyperElasticLaw::m_invariants[0];
147 798003 const double I2 = HyperElasticLaw::m_invariants[1];
148 798003 const double I3 = HyperElasticLaw::m_invariants[2];
149 798003 const double minus_one_third = -1.0/3.0;
150 798003 const double one_ninth = 1.0/9.0;
151
152
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 798003 times.
 798003 FEL_ASSERT("I3 should be strictly positive!" && std::abs(I3) > std::numeric_limits<double>::epsilon());
153
154 798003 return m_kappa1 * I1 * 4.0 * one_ninth * std::pow(I3, 7.0 * minus_one_third)
155 798003 + m_kappa2 * I2 * 10.0 * one_ninth * std::pow(I3, 8.0 * minus_one_third);
156 }
157
158 /***********************************************************************************/
159 /***********************************************************************************/
160
161 798003 double CiarletGeymonat::SecondDerivativeWFirstAndThirdInvariant() const
162 {
163 798003 BaseClassType::SecondDerivativeWFirstAndThirdInvariant();
164
165 798003 const double I3 = HyperElasticLaw::m_invariants[2];
166 798003 const double minus_one_third = -1.0/3.0;
167
168 798003 return m_kappa1 * minus_one_third * std::pow(I3, 4.0 * minus_one_third);
169 }
170
171 /***********************************************************************************/
172 /***********************************************************************************/
173
174 798003 double CiarletGeymonat::SecondDerivativeWSecondAndThirdInvariant() const
175 {
176 798003 BaseClassType::SecondDerivativeWSecondAndThirdInvariant();
177 798003 const double I3 = HyperElasticLaw::m_invariants[2];
178 798003 const double minus_one_third = -1.0/3.0;
179
180 798003 return m_kappa2 * 2. * minus_one_third * std::pow(I3, 5.0 * minus_one_third);
181 }
182
183 /***********************************************************************************/
184 /***********************************************************************************/
185
186 3 double CiarletGeymonat::Wc(const double pressure) const
187 {
188 3 BaseClassType::Wc(pressure);
189 3 const double I3 = HyperElasticLaw::m_invariants[2];
190
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
 3 FEL_ASSERT("I3 should be strictly positive! (and more precisely somewhere in the vicinity of 1.)" && std::abs(I3) > std::numeric_limits<double>::epsilon());
191
192 3 return pressure * (std::sqrt(I3) - 1.0 - 0.5 * std::log(I3));
193 }
194
195 /***********************************************************************************/
196 /***********************************************************************************/
197
198 3 double CiarletGeymonat::FirstDerivativeWcThirdInvariant(const double pressure) const
199 {
200 3 BaseClassType::FirstDerivativeWcThirdInvariant(pressure);
201 3 const double I3 = HyperElasticLaw::m_invariants[2];
202
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
 3 FEL_ASSERT("I3 should be strictly positive! (and more precisely somewhere in the vicinity of 1.)" && std::abs(I3) > std::numeric_limits<double>::epsilon());
203
204 3 return pressure * (0.5/std::sqrt(I3) - 0.5/I3);
205 }
206
207 /***********************************************************************************/
208 /***********************************************************************************/
209
210 3 double CiarletGeymonat::SecondDerivativeWcThirdInvariant(const double pressure) const
211 {
212 3 BaseClassType::SecondDerivativeWcThirdInvariant(pressure);
213 3 const double I3 = HyperElasticLaw::m_invariants[2];
214
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
 3 FEL_ASSERT("I3 should be strictly positive! (and more precisely somewhere in the vicinity of 1.)" && std::abs(I3) > std::numeric_limits<double>::epsilon());
215
216 3 return pressure * (-0.25 * std::pow(I3, -1.5) + 0.5/std::pow(I3, 2));
217 }
218
219 } // namespace felisce
220
221