Loading [MathJax]/extensions/tex2jax.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
IncompressibleLinElast.cpp
Go to the documentation of this file.
2
3namespace polyfem::assembler
4{
5 using namespace basis;
6
7 void IncompressibleLinearElasticityDispacement::add_multimaterial(const int index, const json &params, const Units &units)
8 {
9 assert(size() == 2 || size() == 3);
10
11 params_.add_multimaterial(index, params, size() == 3, units.stress());
12 }
13
14 Eigen::Matrix<double, Eigen::Dynamic, 1, 0, 9, 1>
16 {
17 // 2mu (epsi : epsj)
18 assert(size() == 2 || size() == 3);
19
20 Eigen::Matrix<double, Eigen::Dynamic, 1, 0, 9, 1> res(size() * size());
21 res.setZero();
22
23 const Eigen::MatrixXd &gradi = data.vals.basis_values[data.i].grad_t_m;
24 const Eigen::MatrixXd &gradj = data.vals.basis_values[data.j].grad_t_m;
25
26 Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> epsi(size(), size());
27 Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> epsj(size(), size());
28
29 for (long p = 0; p < gradi.rows(); ++p)
30 {
31
32 double lambda, mu;
33 params_.lambda_mu(data.vals.quadrature.points.row(p), data.vals.val.row(p), data.t, data.vals.element_id, lambda, mu);
34
35 for (int di = 0; di < size(); ++di)
36 {
37 epsi.setZero();
38 epsi.row(di) = gradi.row(p);
39 epsi = ((epsi + epsi.transpose()) / 2).eval();
40 for (int dj = 0; dj < size(); ++dj)
41 {
42 epsj.setZero();
43 epsj.row(dj) = gradj.row(p);
44 epsj = ((epsj + epsj.transpose()) / 2.0).eval();
45
46 res(dj * size() + di) += 2 * mu * (epsi.array() * epsj.array()).sum() * data.da(p);
47 }
48 }
49 }
50
51 return res;
52 }
53
55 const int all_size,
56 const ElasticityTensorType &type,
57 Eigen::MatrixXd &all,
58 const std::function<Eigen::MatrixXd(const Eigen::MatrixXd &)> &fun) const
59 {
60 const auto &displacement = data.fun;
61 const auto &local_pts = data.local_pts;
62 const auto &bs = data.bs;
63 const auto &gbs = data.gbs;
64 const auto el_id = data.el_id;
65
66 assert(size() == 2 || size() == 3);
67 all.resize(local_pts.rows(), all_size);
68 assert(displacement.cols() == 1);
69
70 Eigen::MatrixXd displacement_grad(size(), size());
71
73 vals.compute(el_id, size() == 3, local_pts, bs, gbs);
74
75 for (long p = 0; p < local_pts.rows(); ++p)
76 {
77 compute_diplacement_grad(size(), bs, vals, local_pts, p, displacement, displacement_grad);
78
79 if (type == ElasticityTensorType::F)
80 {
81 all.row(p) = fun(displacement_grad + Eigen::MatrixXd::Identity(size(), size()));
82 continue;
83 }
84
85 double lambda, mu;
86 params_.lambda_mu(local_pts.row(p), vals.val.row(p), data.t, vals.element_id, lambda, mu);
87
88 const Eigen::MatrixXd strain = (displacement_grad + displacement_grad.transpose()) / 2;
89 Eigen::MatrixXd stress = 2 * mu * strain + lambda * strain.trace() * Eigen::MatrixXd::Identity(size(), size());
90
91 if (type == ElasticityTensorType::PK1)
92 stress = pk1_from_cauchy(stress, displacement_grad + Eigen::MatrixXd::Identity(size(), size()));
93 else if (type == ElasticityTensorType::PK2)
94 stress = pk2_from_cauchy(stress, displacement_grad + Eigen::MatrixXd::Identity(size(), size()));
95
96 all.row(p) = fun(stress);
97 }
98 }
99
100 std::map<std::string, Assembler::ParamFunc> IncompressibleLinearElasticityDispacement::parameters() const
101 {
102 std::map<std::string, ParamFunc> res;
103 const auto &params = params_;
104 const int size = this->size();
105
106 res["lambda"] = [&params](const RowVectorNd &uv, const RowVectorNd &p, double t, int e) {
107 double lambda, mu;
108
109 params.lambda_mu(uv, p, t, e, lambda, mu);
110 return lambda;
111 };
112
113 res["mu"] = [&params](const RowVectorNd &uv, const RowVectorNd &p, double t, int e) {
114 double lambda, mu;
115
116 params.lambda_mu(uv, p, t, e, lambda, mu);
117 return mu;
118 };
119
120 res["E"] = [&params, size](const RowVectorNd &uv, const RowVectorNd &p, double t, int e) {
121 double lambda, mu;
122 params.lambda_mu(uv, p, t, e, lambda, mu);
123
124 if (size == 3)
125 return mu * (3.0 * lambda + 2.0 * mu) / (lambda + mu);
126 else
127 return 2 * mu * (2.0 * lambda + 2.0 * mu) / (lambda + 2.0 * mu);
128 };
129
130 res["nu"] = [&params, size](const RowVectorNd &uv, const RowVectorNd &p, double t, int e) {
131 double lambda, mu;
132
133 params.lambda_mu(uv, p, t, e, lambda, mu);
134
135 if (size == 3)
136 return lambda / (2.0 * (lambda + mu));
137 else
138 return lambda / (lambda + 2.0 * mu);
139 };
140
141 return res;
142 }
143
144 Eigen::Matrix<double, Eigen::Dynamic, 1, 0, 3, 1>
146 {
147 // (psii : div phij) = -psii * gradphij
148 assert(size() == 2 || size() == 3);
149 Eigen::Matrix<double, Eigen::Dynamic, 1, 0, 3, 1> res(rows() * cols());
150 res.setZero();
151
152 const Eigen::MatrixXd &psii = data.psi_vals.basis_values[data.i].val;
153 const Eigen::MatrixXd &gradphij = data.phi_vals.basis_values[data.j].grad_t_m;
154 assert(psii.size() == gradphij.rows());
155 assert(gradphij.cols() == rows());
156
157 for (int k = 0; k < gradphij.rows(); ++k)
158 {
159 res -= psii(k) * gradphij.row(k) * data.da(k);
160 }
161
162 return res;
163 }
164
165 void IncompressibleLinearElasticityPressure::add_multimaterial(const int index, const json &params, const Units &units)
166 {
167 assert(disp_size_ == 2 || disp_size_ == 3);
168
169 params_.add_multimaterial(index, params, disp_size_ == 3, units.stress());
170 }
171
172 Eigen::Matrix<double, Eigen::Dynamic, 1, 0, 9, 1>
174 {
175 // -1/lambda phi_ * phi_j
176
177 const Eigen::MatrixXd &phii = data.vals.basis_values[data.i].val;
178 const Eigen::MatrixXd &phij = data.vals.basis_values[data.j].val;
179
180 double res = 0;
181
182 for (long p = 0; p < data.da.size(); ++p)
183 {
184 double lambda, mu;
185 params_.lambda_mu(data.vals.quadrature.points.row(p), data.vals.val.row(p), data.t, data.vals.element_id, lambda, mu);
186
187 res += -phii(p) * phij(p) * data.da(p) / lambda;
188 }
189
190 // double res = (phii.array() * phij.array() * da.array()).sum();
191 // res *= -1. / lambda;
192
193 return Eigen::Matrix<double, 1, 1>::Constant(res);
194 }
195} // namespace polyfem::assembler
ElementAssemblyValues vals
Definition Assembler.cpp:22
std::string stress() const
Definition Units.hpp:27
stores per element basis values at given quadrature points and geometric mapping
void compute(const int el_index, const bool is_volume, const Eigen::MatrixXd &pts, const basis::ElementBases &basis, const basis::ElementBases &gbasis)
computes the per element values at the local (ref el) points (pts) sets basis_values,...
void add_multimaterial(const int index, const json &params, const Units &units) override
Eigen::Matrix< double, Eigen::Dynamic, 1, 0, 9, 1 > assemble(const LinearAssemblerData &data) const override
local assembly function that defines the bilinear form (LHS) computes and returns a single local stif...
std::map< std::string, ParamFunc > parameters() const override
void assign_stress_tensor(const OutputData &data, const int all_size, const ElasticityTensorType &type, Eigen::MatrixXd &all, const std::function< Eigen::MatrixXd(const Eigen::MatrixXd &)> &fun) const override
Eigen::Matrix< double, Eigen::Dynamic, 1, 0, 3, 1 > assemble(const MixedAssemblerData &data) const override
void add_multimaterial(const int index, const json &params, const Units &units) override
Eigen::Matrix< double, Eigen::Dynamic, 1, 0, 9, 1 > assemble(const LinearAssemblerData &data) const override
local assembly function that defines the bilinear form (LHS) computes and returns a single local stif...
void lambda_mu(double px, double py, double pz, double x, double y, double z, double t, int el_id, double &lambda, double &mu) const
void add_multimaterial(const int index, const json &params, const bool is_volume, const std::string &stress_unit)
const ElementAssemblyValues & vals
stores the evaluation for that element
const QuadratureVector & da
contains both the quadrature weight and the change of metric in the integral
const QuadratureVector & da
contains both the quadrature weight and the change of metric in the integral
const ElementAssemblyValues & phi_vals
stores the evaluation for that element
const ElementAssemblyValues & psi_vals
stores the evaluation for that element
const basis::ElementBases & bs
const Eigen::MatrixXd & fun
const Eigen::MatrixXd & local_pts
const basis::ElementBases & gbs
Used for test only.
Eigen::MatrixXd pk2_from_cauchy(const Eigen::MatrixXd &stress, const Eigen::MatrixXd &F)
nlohmann::json json
Definition Common.hpp:9
void compute_diplacement_grad(const int size, const ElementAssemblyValues &vals, const Eigen::MatrixXd &local_pts, const int p, const Eigen::MatrixXd &displacement, Eigen::MatrixXd &displacement_grad)
Eigen::MatrixXd pk1_from_cauchy(const Eigen::MatrixXd &stress, const Eigen::MatrixXd &F)
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > RowVectorNd
Definition Types.hpp:13