waves
Basic FE playground
StressComputation.h
Go to the documentation of this file.
1 #ifndef KATOPTRON_STRESSCOMPUTATION_H
2 #define KATOPTRON_STRESSCOMPUTATION_H
3 
4 #include "wHex8.h"
5 #include "wSfHex8.h"
6 #include "wGaussHex8.h"
7 #include "wCacheHex8.h"
8 
9 #include "wTetra4.h"
10 #include "wSfTetra4.h"
11 #include "wGaussTetra4.h"
12 #include "wCacheTetra4.h"
13 
20 void strain_at_gp(size_t n_n, tbox::Element &el, tbox::Cache &cache, size_t a, const Eigen::MatrixXd &u_e, Eigen::MatrixXd &Strain)
21 {
22  Eigen::MatrixXd const &J = el.getJinv(a);
23  Eigen::MatrixXd const &dffi = cache.getDsf(a);
24  Eigen::Matrix3d du_d_Phi = Eigen::Matrix3d::Zero();
25  for (size_t i = 0; i < n_n; ++i)
26  {
27  for (auto k = 0; k < 3; ++k)
28  for (auto j = 0; j < 3; ++j)
29  du_d_Phi(j, k) += dffi(j, i) * u_e(k, i);
30  }
31  Strain += (J * du_d_Phi).transpose() * (J * du_d_Phi) * 0.5;
32 }
33 
34 void compute_stress(size_t n_n, const Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress, const Eigen::MatrixXd &H, const Eigen::MatrixXd &D, const std::vector<double> &T, bool thermal)
35 {
36  Stress = Eigen::Matrix3d::Zero();
37  for (size_t i = 0; i < 3; ++i)
38  for (size_t j = 0; j < 3; ++j)
39  for (size_t k = 0; k < 3; ++k)
40  for (size_t l = 0; l < 3; ++l)
41  Stress(i, j) += H(i * 3 + j, k * 3 + l) * Strain(k, l);
42 
43  if (thermal)
44  for (size_t i = 0; i < 3; ++i)
45  for (size_t j = 0; j < 3; ++j)
46  {
47  double tmp = 0.;
48  for (size_t k = 0; k < n_n; ++k)
49  tmp += D(i, j) * T[k];
50  tmp /= n_n;
51  Stress(i, j) += tmp;
52  }
53 }
54 
55 void strain_stress_x_y_z(tbox::Element *elem, Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress, const Eigen::MatrixXd &H, const Eigen::MatrixXd &D, const std::vector<double> &x, const std::vector<double> &y, const std::vector<double> &z, const std::vector<double> &T, bool thermal = false)
56 {
57  tbox::Cache &cache = elem->getVCache();
58 
59  size_t n_gp = cache.getVGauss().getN();
60  size_t n_n = elem->nodes.size();
61 
62  Eigen::MatrixXd u_e(3, n_n);
63 
64  for (size_t i = 0; i < n_n; ++i)
65  {
66  u_e(0, i) = x[i];
67  u_e(1, i) = y[i];
68  u_e(2, i) = z[i];
69  }
70 
71  Strain = Eigen::Matrix3d::Zero();
72  for (size_t a = 0; a < n_gp; ++a)
73  {
74  strain_at_gp(n_n, *elem, cache, a, u_e, Strain);
75  }
76  Strain /= n_n;
77 
78  compute_stress(n_n, Strain, Stress, H, D, T, thermal);
79 }
80 
81 void strain_stress_x_y_z_at_node(tbox::Element *elem, std::vector<Eigen::MatrixXd> &Strain_at_node, std::vector<Eigen::MatrixXd> &Stress_at_node, const Eigen::MatrixXd &H, const Eigen::MatrixXd &D, const std::vector<double> &x, const std::vector<double> &y, const std::vector<double> &z, const std::vector<double> &T, bool thermal = false)
82 {
83  tbox::Cache &cache = elem->getVCache();
84 
85  size_t n_gp = cache.getVGauss().getN();
86  size_t n_n = elem->nodes.size();
87 
88  std::vector<Eigen::MatrixXd> Strain_at_gp(n_gp);
89  std::vector<Eigen::MatrixXd> Stress_at_gp(n_gp);
90  Eigen::MatrixXd u_e(3, n_n);
91  Eigen::MatrixXd Strain(3, 3);
92 
93  for (auto &k : Strain_at_gp)
94  k.resize(3, 3);
95  for (auto &k : Stress_at_gp)
96  k.resize(3, 3);
97 
98  for (size_t i = 0; i < n_n; ++i)
99  {
100  u_e(0, i) = x[i];
101  u_e(1, i) = y[i];
102  u_e(2, i) = z[i];
103  }
104 
105  for (size_t a = 0; a < n_gp; ++a)
106  {
107  Strain = Eigen::Matrix3d::Zero();
108  strain_at_gp(n_n, *elem, cache, a, u_e, Strain);
109  Strain_at_gp[a] = Strain;
110  compute_stress(n_n, Strain, Stress_at_gp[a], H, D, T, thermal);
111  }
112 
113  // Extrapolation at node points
114  Eigen::MatrixXd values_at_gp(n_gp, 18);
115  for (size_t i = 0; i < 3; ++i)
116  for (size_t j = 0; j < 3; ++j)
117  for (size_t a = 0; a < n_gp; ++a)
118  {
119  values_at_gp(a, i * 3 + j) = Strain_at_gp[a](i, j);
120  values_at_gp(a, 9 + i * 3 + j) = Stress_at_gp[a](i, j);
121  }
122 
123  Eigen::MatrixXd values_at_node = elem->gp2node(values_at_gp);
124  for (size_t k = 0; k < n_n; ++k)
125  for (size_t i = 0; i < 3; ++i)
126  for (size_t j = 0; j < 3; ++j)
127  {
128  Strain_at_node[k](i, j) = values_at_node(k, i * 3 + j);
129  Stress_at_node[k](i, j) = values_at_node(k, 9 + i * 3 + j);
130  }
131 }
132 
133 #endif //KATOPTRON_STRESSCOMPUTATION_H
strain_at_gp
void strain_at_gp(size_t n_n, tbox::Element &el, tbox::Cache &cache, size_t a, const Eigen::MatrixXd &u_e, Eigen::MatrixXd &Strain)
Those functions compute the stress and strain at nodes and in average per element....
Definition: StressComputation.h:20
compute_stress
void compute_stress(size_t n_n, const Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress, const Eigen::MatrixXd &H, const Eigen::MatrixXd &D, const std::vector< double > &T, bool thermal)
Definition: StressComputation.h:34
strain_stress_x_y_z_at_node
void strain_stress_x_y_z_at_node(tbox::Element *elem, std::vector< Eigen::MatrixXd > &Strain_at_node, std::vector< Eigen::MatrixXd > &Stress_at_node, const Eigen::MatrixXd &H, const Eigen::MatrixXd &D, const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &z, const std::vector< double > &T, bool thermal=false)
Definition: StressComputation.h:81
strain_stress_x_y_z
void strain_stress_x_y_z(tbox::Element *elem, Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress, const Eigen::MatrixXd &H, const Eigen::MatrixXd &D, const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &z, const std::vector< double > &T, bool thermal=false)
Definition: StressComputation.h:55