1 #ifndef KATOPTRON_STRESSCOMPUTATION_H
2 #define KATOPTRON_STRESSCOMPUTATION_H
6 #include "wGaussHex8.h"
7 #include "wCacheHex8.h"
10 #include "wSfTetra4.h"
11 #include "wGaussTetra4.h"
12 #include "wCacheTetra4.h"
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)
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)
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);
31 Strain += (J * du_d_Phi).transpose() * (J * du_d_Phi) * 0.5;
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)
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);
44 for (
size_t i = 0; i < 3; ++i)
45 for (
size_t j = 0; j < 3; ++j)
48 for (
size_t k = 0; k < n_n; ++k)
49 tmp += D(i, j) * T[k];
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)
57 tbox::Cache &cache = elem->getVCache();
59 size_t n_gp = cache.getVGauss().getN();
60 size_t n_n = elem->nodes.size();
62 Eigen::MatrixXd u_e(3, n_n);
64 for (
size_t i = 0; i < n_n; ++i)
71 Strain = Eigen::Matrix3d::Zero();
72 for (
size_t a = 0; a < n_gp; ++a)
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)
83 tbox::Cache &cache = elem->getVCache();
85 size_t n_gp = cache.getVGauss().getN();
86 size_t n_n = elem->nodes.size();
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);
93 for (
auto &k : Strain_at_gp)
95 for (
auto &k : Stress_at_gp)
98 for (
size_t i = 0; i < n_n; ++i)
105 for (
size_t a = 0; a < n_gp; ++a)
107 Strain = Eigen::Matrix3d::Zero();
109 Strain_at_gp[a] = Strain;
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)
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);
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)
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);
133 #endif //KATOPTRON_STRESSCOMPUTATION_H