MAST3
Multidisciplinary-design Adaptation and Sensitivity Toolkit (MAST)
von_mises_stress.hpp
Go to the documentation of this file.
1 /*
2 * MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
3 * Copyright (C) 2013-2020 Manav Bhatia and MAST authors
4 *
5 * This library is free software; you can redistribute it and/or
6 * modify it under the terms of the GNU Lesser General Public
7 * License as published by the Free Software Foundation; either
8 * version 2.1 of the License, or (at your option) any later version.
9 *
10 * This library is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * Lesser General Public License for more details.
14 *
15 * You should have received a copy of the GNU Lesser General Public
16 * License along with this library; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18 */
19 
20 #ifndef __mast_linear_continuum_von_mises_stress_h__
21 #define __mast_linear_continuum_von_mises_stress_h__
22 
23 // MAST includes
24 #include <mast/physics/elasticity/linear_elastic_strain_operator.hpp>
25 
26 namespace MAST {
27 namespace Physics {
28 namespace Elasticity {
29 namespace LinearContinuum {
30 
31 
32 template <typename ScalarType, uint_t Dim>
33 using stress_vec_t = Eigen::Matrix<ScalarType,
34  MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
35  1>;
36 
37 template <typename ScalarType, uint_t Dim>
38 using stress_adjoint_mat_t = Eigen::Matrix<ScalarType,
39  MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
40  Eigen::Dynamic>;
41 
42 
43 
44 template <typename ScalarType, uint_t Dim, typename VecType>
45 inline
46 typename std::enable_if<Dim==2, ScalarType>::type
47 vonMises_stress(const VecType& stress) {
48 
49  Assert1(stress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
50  stress.size(),
51  "Incorrect stress vector dimension");
52 
53  return
54  pow(0.5 * (pow(stress(0)-stress(1),2) + //(((sigma_xx - sigma_yy)^2 +
55  pow(stress(1),2) + // (sigma_yy)^2 +
56  pow(stress(0),2)) + // (sigma_xx)^2)/2 +
57  3.0 * (pow(stress(2), 2)), 0.5); // 3.0 * tau_xy^2)^.5
58 }
59 
60 
61 
62 template <typename ScalarType, uint_t Dim, typename VecType>
63 inline
64 typename std::enable_if<Dim==3, ScalarType>::type
65 vonMises_stress(const VecType& stress) {
66 
67  Assert1(stress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
68  stress.size(),
69  "Incorrect stress vector dimension");
70 
71  return
72  pow(0.5 * (pow(stress(0)-stress(1),2) + //(((sigma_xx - sigma_yy)^2 +
73  pow(stress(1)-stress(2),2) + // (sigma_yy - sigma_zz)^2 +
74  pow(stress(2)-stress(0),2)) + // (sigma_zz - sigma_xx)^2)/2 +
75  3.0 * (pow(stress(3), 2) + // 3* (tau_xy^2 +
76  pow(stress(4), 2) + // tau_yz^2 +
77  pow(stress(5), 2)), 0.5); // tau_zx^2))^.5
78 }
79 
80 
81 
82 template <typename ScalarType, uint_t Dim, typename Vec1Type, typename Vec2Type>
83 inline
84 typename std::enable_if<Dim==2, ScalarType>::type
85 vonMises_stress_derivative(const Vec1Type& stress,
86  const Vec2Type& dstress_dp) {
87 
88  Assert1(stress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
89  stress.size(),
90  "Incorrect stress vector dimension");
91  Assert1(dstress_dp.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
92  dstress_dp.size(),
93  "Incorrect stress vector dimension");
94 
95  ScalarType
96  p =
97  0.5 * (pow(stress(0)-stress(1),2) + //((sigma_xx - sigma_yy)^2 +
98  pow(stress(1),2) + // (sigma_yy)^2 +
99  pow(stress(0),2)) + // (sigma_xx)^2)/2 +
100  3.0 * (pow(stress(2), 2)), // 3.0 * tau_xy^2)
101  dp = 0.;
102 
103  // if p == 0, then the sensitivity returns nan
104  // Hence, we are avoiding this by setting it to zero whenever p = 0.
105  if (fabs(p) > 0.)
106  dp =
107  (((dstress_dp(0) - dstress_dp(1)) * (stress(0) - stress(1)) +
108  (dstress_dp(1)) * (stress(1)) +
109  (dstress_dp(0)) * (stress(0))) +
110  6.0 * (dstress_dp(2) * stress(2))) * 0.5 * pow(p, -0.5);
111 
112  return dp;
113 }
114 
115 
116 
117 
118 template <typename ScalarType, uint_t Dim, typename Vec1Type, typename Vec2Type>
119 inline
120 typename std::enable_if<Dim==3, ScalarType>::type
121 vonMises_stress_derivative(const Vec1Type& stress,
122  const Vec2Type& dstress_dp) {
123 
124  Assert1(stress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
125  stress.size(),
126  "Incorrect stress vector dimension");
127  Assert1(dstress_dp.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
128  dstress_dp.size(),
129  "Incorrect stress vector dimension");
130 
131  ScalarType
132  p =
133  0.5 * (pow(stress(0)-stress(1),2) + //((sigma_xx - sigma_yy)^2 +
134  pow(stress(1)-stress(2),2) + // (sigma_yy - sigma_zz)^2 +
135  pow(stress(2)-stress(0),2)) + // (sigma_zz - sigma_xx)^2)/2 +
136  3.0 * (pow(stress(3), 2) + // 3* (tau_xy^2 +
137  pow(stress(4), 2) + // tau_yz^2 +
138  pow(stress(5), 2)), // tau_zx^2)
139  dp = 0.;
140 
141  // if p == 0, then the sensitivity returns nan
142  // Hence, we are avoiding this by setting it to zero whenever p = 0.
143  if (fabs(p) > 0.)
144  dp =
145  (((dstress_dp(0) - dstress_dp(1)) * (stress(0) - stress(1)) +
146  (dstress_dp(1) - dstress_dp(2)) * (stress(1) - stress(2)) +
147  (dstress_dp(2) - dstress_dp(0)) * (stress(2) - stress(0))) +
148  6.0 * (dstress_dp(3) * stress(3)+
149  dstress_dp(4) * stress(4)+
150  dstress_dp(5) * stress(5))) * 0.5 * pow(p, -0.5);
151 
152  return dp;
153 }
154 
155 
156 
157 template <typename ScalarType,
158  uint_t Dim,
159  typename VecType,
160  typename AdjMatType>
161 inline
162 typename std::enable_if<Dim==2, void>::type
163 vonMises_stress_dX(const VecType &stress,
164  const AdjMatType &dstress_dX,
165  Eigen::Matrix<ScalarType, Eigen::Dynamic, 1> &vm_adjoint) {
166 
167  Assert1(stress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
168  stress.size(),
169  "Incorrect stress vector dimension");
170  Assert1(dstress_dX.rows() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
171  dstress_dX.rows(),
172  "Incorrect stress adjoint matrix rows");
173  Assert2(vm_adjoint.size() == dstress_dX.cols(),
174  vm_adjoint.size(), dstress_dX.cols(),
175  "Incorrect von Mises adjoint vector dimension");
176 
177  ScalarType
178  p =
179  0.5 * (pow(stress(0)-stress(1),2) + //((sigma_xx - sigma_yy)^2 +
180  pow(stress(1),2) + // (sigma_yy)^2 +
181  pow(stress(0),2)) + // (sigma_xx)^2)/2 +
182  3.0 * (pow(stress(2), 2)); // 3* (tau_xy^2)
183 
184 
185  // if p == 0, then the sensitivity returns nan
186  // Hence, we are avoiding this by setting it to zero whenever p = 0.
187  if (fabs(p) > 0.)
188  vm_adjoint =
189  (((dstress_dX.row(0) - dstress_dX.row(1)) * (stress(0) - stress(1)) +
190  dstress_dX.row(1) * stress(1) +
191  dstress_dX.row(0) * stress(0)) +
192  6.0 * dstress_dX.row(2) * stress(2)) * 0.5 * pow(p, -0.5);
193 }
194 
195 
196 
197 template <typename ScalarType,
198  uint_t Dim,
199  typename VecType,
200  typename AdjMatType>
201 inline
202 typename std::enable_if<Dim==3, void>::type
203 vonMises_stress_dX(const VecType &stress,
204  const AdjMatType &dstress_dX,
205  Eigen::Matrix<ScalarType, Eigen::Dynamic, 1> &vm_adjoint) {
206 
207  Assert1(stress.size() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
208  stress.size(),
209  "Incorrect stress vector dimension");
210  Assert1(dstress_dX.rows() == MAST::Physics::Elasticity::LinearContinuum::NStrainComponents<Dim>::value,
211  dstress_dX.rows(),
212  "Incorrect stress adjoint matrix rows");
213  Assert2(vm_adjoint.size() == dstress_dX.cols(),
214  vm_adjoint.size(), dstress_dX.cols(),
215  "Incorrect von Mises adjoint vector dimension");
216 
217  ScalarType
218  p =
219  0.5 * (pow(stress(0)-stress(1),2) + //((sigma_xx - sigma_yy)^2 +
220  pow(stress(1)-stress(2),2) + // (sigma_yy - sigma_zz)^2 +
221  pow(stress(2)-stress(0),2)) + // (sigma_zz - sigma_xx)^2)/2 +
222  3.0 * (pow(stress(3), 2) + // 3* (tau_xx^2 +
223  pow(stress(4), 2) + // tau_yy^2 +
224  pow(stress(5), 2)); // tau_zz^2)
225 
226 
227  // if p == 0, then the sensitivity returns nan
228  // Hence, we are avoiding this by setting it to zero whenever p = 0.
229  if (fabs(p) > 0.)
230  vm_adjoint =
231  (((dstress_dX.row(0) - dstress_dX.row(1)) * (stress(0) - stress(1)) +
232  (dstress_dX.row(1) - dstress_dX.row(2)) * (stress(1) - stress(2)) +
233  (dstress_dX.row(2) - dstress_dX.row(0)) * (stress(2) - stress(0))) +
234  6.0 * (dstress_dX.row(3) * stress(3)+
235  dstress_dX.row(4) * stress(4)+
236  dstress_dX.row(5) * stress(5))) * 0.5 * pow(p, -0.5);
237 }
238 
239 
240 } // namespace LinearContinuum
241 } // namespace Elasticity
242 } // namespace Physics
243 } // namespace MAST
244 
245 #endif // __mast_linear_continuum_von_mises_stress_h__
MAST::Physics::Elasticity::LinearContinuum::stress_adjoint_mat_t
Eigen::Matrix< ScalarType, MAST::Physics::Elasticity::LinearContinuum::NStrainComponents< Dim >::value, Eigen::Dynamic > stress_adjoint_mat_t
Definition: von_mises_stress.hpp:40
Assert1
#define Assert1(cond, v1, msg)
Definition: exceptions.hpp:143
MAST::Physics::Elasticity::LinearContinuum::vonMises_stress
std::enable_if< Dim==2, ScalarType >::type vonMises_stress(const VecType &stress)
Definition: von_mises_stress.hpp:47
MAST::Physics::Elasticity::LinearContinuum::vonMises_stress_derivative
std::enable_if< Dim==2, ScalarType >::type vonMises_stress_derivative(const Vec1Type &stress, const Vec2Type &dstress_dp)
Definition: von_mises_stress.hpp:85
uint_t
unsigned int uint_t
Definition: mast_data_types.h:43
Assert2
#define Assert2(cond, v1, v2, msg)
Definition: exceptions.hpp:152
MAST::Physics::Elasticity::LinearContinuum::stress_vec_t
Eigen::Matrix< ScalarType, MAST::Physics::Elasticity::LinearContinuum::NStrainComponents< Dim >::value, 1 > stress_vec_t
Definition: von_mises_stress.hpp:35
linear_elastic_strain_operator.hpp
MAST::Physics::Elasticity::LinearContinuum::vonMises_stress_dX
std::enable_if< Dim==2, void >::type vonMises_stress_dX(const VecType &stress, const AdjMatType &dstress_dX, Eigen::Matrix< ScalarType, Eigen::Dynamic, 1 > &vm_adjoint)
Definition: von_mises_stress.hpp:163