PolyFEM
Loading...
Searching...
No Matches
ElasticForceDerivative.cpp
Go to the documentation of this file.
1#include "ElasticForceDerivative.hpp"
2
3#include <Eigen/Core>
4#include <polyfem/solver/forms/ElasticForm.hpp>
5#include <polyfem/assembler/ElementAssemblyValues.hpp>
6#include <polyfem/assembler/MatParams.hpp>
7#include <polyfem/assembler/ViscousDamping.hpp>
8#include <polyfem/io/Evaluator.hpp>
9#include <polyfem/quadrature/Quadrature.hpp>
10#include <polyfem/utils/MatrixUtils.hpp>
11#include <polyfem/utils/MaybeParallelFor.hpp>
12#include <polyfem/utils/Types.hpp>
13
14using namespace polyfem::assembler;
15
16namespace polyfem::solver
17{
18 namespace
19 {
20 class LocalThreadVecStorage
21 {
22 public:
23 Eigen::MatrixXd vec;
24 assembler::ElementAssemblyValues vals;
25 QuadratureVector da;
26
27 LocalThreadVecStorage(const int size)
28 {
29 vec.resize(size, 1);
30 vec.setZero();
31 }
32 };
33 } // namespace
34
35 void ElasticForceDerivative::force_material_derivative(
36 ElasticForm &form,
37 const double t,
38 const Eigen::MatrixXd &x,
39 const Eigen::MatrixXd &x_prev,
40 const Eigen::MatrixXd &adjoint,
41 Eigen::VectorXd &term)
42 {
43 const int dim = form.is_volume_ ? 3 : 2;
44
45 const int n_elements = int(form.bases_.size());
46
47 if (form.assembler_.name() == "ViscousDamping")
48 {
49 term.setZero(2);
50
51 auto storage = utils::create_thread_storage(LocalThreadVecStorage(term.size()));
52
53 utils::maybe_parallel_for(n_elements, [&](int start, int end, int thread_id) {
54 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
55
56 for (int e = start; e < end; ++e)
57 {
58 assembler::ElementAssemblyValues &vals = local_storage.vals;
59 form.ass_vals_cache_.compute(e, form.is_volume_, form.bases_[e], form.geom_bases_[e], vals);
60
61 const quadrature::Quadrature &quadrature = vals.quadrature;
62 local_storage.da = vals.det.array() * quadrature.weights.array();
63
64 Eigen::MatrixXd u, grad_u, prev_u, prev_grad_u, p, grad_p;
65 io::Evaluator::interpolate_at_local_vals(e, dim, dim, vals, x, u, grad_u);
66 io::Evaluator::interpolate_at_local_vals(e, dim, dim, vals, x_prev, prev_u, prev_grad_u);
67 io::Evaluator::interpolate_at_local_vals(e, dim, dim, vals, adjoint, p, grad_p);
68
69 for (int q = 0; q < local_storage.da.size(); ++q)
70 {
71 Eigen::MatrixXd grad_p_i, grad_u_i, prev_grad_u_i;
72 utils::vector2matrix(grad_p.row(q), grad_p_i);
73 utils::vector2matrix(grad_u.row(q), grad_u_i);
74 utils::vector2matrix(prev_grad_u.row(q), prev_grad_u_i);
75
76 Eigen::MatrixXd f_prime_dpsi, f_prime_dphi;
77 assembler::ViscousDamping::compute_dstress_dpsi_dphi(
78 OptAssemblerData(t, form.dt_, e, quadrature.points.row(q), vals.val.row(q), grad_u_i),
79 prev_grad_u_i, f_prime_dpsi, f_prime_dphi);
80
81 // This needs to be a sum over material parameter basis.
82 local_storage.vec(0) += -utils::matrix_inner_product<double>(f_prime_dpsi, grad_p_i) * local_storage.da(q);
83 local_storage.vec(1) += -utils::matrix_inner_product<double>(f_prime_dphi, grad_p_i) * local_storage.da(q);
84 }
85 }
86 });
87
88 for (const LocalThreadVecStorage &local_storage : storage)
89 term += local_storage.vec;
90 }
91 else
92 {
93 term.setZero(n_elements * 2, 1);
94
95 auto storage = utils::create_thread_storage(LocalThreadVecStorage(term.size()));
96
97 utils::maybe_parallel_for(n_elements, [&](int start, int end, int thread_id) {
98 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
99
100 for (int e = start; e < end; ++e)
101 {
102 assembler::ElementAssemblyValues &vals = local_storage.vals;
103 form.ass_vals_cache_.compute(e, form.is_volume_, form.bases_[e], form.geom_bases_[e], vals);
104
105 const quadrature::Quadrature &quadrature = vals.quadrature;
106 local_storage.da = vals.det.array() * quadrature.weights.array();
107
108 Eigen::MatrixXd u, grad_u, p, grad_p;
109 io::Evaluator::interpolate_at_local_vals(e, dim, dim, vals, x, u, grad_u);
110 io::Evaluator::interpolate_at_local_vals(e, dim, dim, vals, adjoint, p, grad_p);
111
112 for (int q = 0; q < local_storage.da.size(); ++q)
113 {
114 Eigen::MatrixXd grad_p_i, grad_u_i;
115 utils::vector2matrix(grad_p.row(q), grad_p_i);
116 utils::vector2matrix(grad_u.row(q), grad_u_i);
117
118 Eigen::MatrixXd f_prime_dmu, f_prime_dlambda;
119 form.assembler_.compute_dstress_dmu_dlambda(
120 OptAssemblerData(t, form.dt_, e, quadrature.points.row(q), vals.val.row(q), grad_u_i),
121 f_prime_dmu, f_prime_dlambda);
122
123 // This needs to be a sum over material parameter basis.
124 local_storage.vec(e + n_elements) += -utils::matrix_inner_product<double>(f_prime_dmu, grad_p_i) * local_storage.da(q);
125 local_storage.vec(e) += -utils::matrix_inner_product<double>(f_prime_dlambda, grad_p_i) * local_storage.da(q);
126 }
127 }
128 });
129
130 for (const LocalThreadVecStorage &local_storage : storage)
131 term += local_storage.vec;
132 }
133 }
134
135 void ElasticForceDerivative::force_shape_derivative(
136 ElasticForm &form,
137 const double t,
138 const int n_verts,
139 const Eigen::MatrixXd &x,
140 const Eigen::MatrixXd &x_prev,
141 const Eigen::MatrixXd &adjoint,
142 Eigen::VectorXd &term)
143 {
144 const int dim = form.is_volume_ ? 3 : 2;
145 const int actual_dim = ((form.assembler_.name() == "Laplacian") || (form.assembler_.name() == "Electrostatics")) ? 1 : dim;
146
147 const int n_elements = int(form.bases_.size());
148 term.setZero(n_verts * dim, 1);
149
150 auto storage = utils::create_thread_storage(LocalThreadVecStorage(term.size()));
151
152 if (form.assembler_.name() == "ViscousDamping")
153 {
154 utils::maybe_parallel_for(n_elements, [&](int start, int end, int thread_id) {
155 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
156
157 for (int e = start; e < end; ++e)
158 {
159 assembler::ElementAssemblyValues &vals = local_storage.vals;
160 form.ass_vals_cache_.compute(e, form.is_volume_, form.bases_[e], form.geom_bases_[e], vals);
161 assembler::ElementAssemblyValues gvals;
162 gvals.compute(e, form.is_volume_, vals.quadrature.points, form.geom_bases_[e], form.geom_bases_[e]);
163
164 const quadrature::Quadrature &quadrature = vals.quadrature;
165 local_storage.da = vals.det.array() * quadrature.weights.array();
166
167 Eigen::MatrixXd u, grad_u, prev_u, prev_grad_u, p, grad_p;
168 io::Evaluator::interpolate_at_local_vals(e, dim, dim, vals, x, u, grad_u);
169 io::Evaluator::interpolate_at_local_vals(e, dim, dim, vals, x_prev, prev_u, prev_grad_u);
170 io::Evaluator::interpolate_at_local_vals(e, dim, dim, vals, adjoint, p, grad_p);
171
172 Eigen::MatrixXd grad_u_i, prev_grad_u_i;
173 Eigen::MatrixXd grad_v_i;
174 Eigen::MatrixXd stress_tensor;
175 Eigen::VectorXd f_prime_gradu_gradv, f_prev_prime_prev_gradu_gradv;
176
177 for (int q = 0; q < local_storage.da.size(); ++q)
178 {
179 utils::vector2matrix(grad_u.row(q), grad_u_i);
180 utils::vector2matrix(prev_grad_u.row(q), prev_grad_u_i);
181
182 for (auto &v : gvals.basis_values)
183 {
184 Eigen::MatrixXd stress_grad, stress_prev_grad;
185 form.assembler_.compute_stress_grad(
186 OptAssemblerData(t, form.dt_, e, quadrature.points.row(q), vals.val.row(q), grad_u_i),
187 prev_grad_u_i, stress_tensor, stress_grad);
188 form.assembler_.compute_stress_prev_grad(
189 OptAssemblerData(t, form.dt_, e, quadrature.points.row(q), vals.val.row(q), grad_u_i),
190 prev_grad_u_i, stress_prev_grad);
191 for (int d = 0; d < dim; d++)
192 {
193 grad_v_i.setZero(dim, dim);
194 grad_v_i.row(d) = v.grad_t_m.row(q);
195
196 f_prime_gradu_gradv = stress_grad * utils::flatten(grad_u_i * grad_v_i);
197 f_prev_prime_prev_gradu_gradv = stress_prev_grad * utils::flatten(prev_grad_u_i * grad_v_i);
198
199 Eigen::MatrixXd tmp = grad_v_i - grad_v_i.trace() * Eigen::MatrixXd::Identity(dim, dim);
200 local_storage.vec(v.global[0].index * dim + d) -= grad_p.row(q).dot(f_prime_gradu_gradv + f_prev_prime_prev_gradu_gradv + utils::flatten(stress_tensor * tmp.transpose())) * local_storage.da(q);
201 }
202 }
203 }
204 }
205 });
206 }
207 else
208 {
209 utils::maybe_parallel_for(n_elements, [&](int start, int end, int thread_id) {
210 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
211
212 for (int e = start; e < end; ++e)
213 {
214 assembler::ElementAssemblyValues &vals = local_storage.vals;
215 form.ass_vals_cache_.compute(e, form.is_volume_, form.bases_[e], form.geom_bases_[e], vals);
216 assembler::ElementAssemblyValues gvals;
217 gvals.compute(e, form.is_volume_, vals.quadrature.points, form.geom_bases_[e], form.geom_bases_[e]);
218
219 const quadrature::Quadrature &quadrature = vals.quadrature;
220 local_storage.da = vals.det.array() * quadrature.weights.array();
221
222 Eigen::MatrixXd u, grad_u, p, grad_p;
223 io::Evaluator::interpolate_at_local_vals(e, dim, actual_dim, vals, x, u, grad_u);
224 io::Evaluator::interpolate_at_local_vals(e, dim, actual_dim, vals, adjoint, p, grad_p);
225
226 for (int q = 0; q < local_storage.da.size(); ++q)
227 {
228 Eigen::MatrixXd grad_u_i, grad_p_i, stiffness_i;
229 if (actual_dim == 1)
230 {
231 grad_u_i = grad_u.row(q);
232 grad_p_i = grad_p.row(q);
233 }
234 else
235 {
236 utils::vector2matrix(grad_u.row(q), grad_u_i);
237 utils::vector2matrix(grad_p.row(q), grad_p_i);
238 }
239
240 for (auto &v : gvals.basis_values)
241 {
242 for (int d = 0; d < dim; d++)
243 {
244 Eigen::MatrixXd grad_v_i;
245 grad_v_i.setZero(dim, dim);
246 grad_v_i.row(d) = v.grad_t_m.row(q);
247
248 Eigen::MatrixXd stress_tensor, f_prime_gradu_gradv;
249 form.assembler_.compute_stress_grad_multiply_mat(
250 OptAssemblerData(t, form.dt_, e, quadrature.points.row(q), vals.val.row(q), grad_u_i),
251 grad_u_i * grad_v_i, stress_tensor, f_prime_gradu_gradv);
252
253 const Eigen::MatrixXd tmp = stress_tensor * grad_v_i.transpose() - grad_v_i.trace() * stress_tensor;
254 local_storage.vec(v.global[0].index * dim + d) -= utils::matrix_inner_product<double>(f_prime_gradu_gradv + tmp, grad_p_i) * local_storage.da(q);
255 }
256 }
257 }
258 }
259 });
260 }
261
262 for (const LocalThreadVecStorage &local_storage : storage)
263 term += local_storage.vec;
264 }
265} // namespace polyfem::solver
vec
Eigen::MatrixXd vec
Definition Assembler.cpp:72
da
QuadratureVector da
Definition Assembler.cpp:23
vals
ElementAssemblyValues vals
Definition Assembler.cpp:22
gvals
assembler::ElementAssemblyValues gvals
Definition BodyForceDerivative.cpp:26
quadrature
Quadrature quadrature
Definition MassMatrixAssembler.cpp:30
x
int x
Definition SplineBasis3d.cpp:55
polyfem::assembler::ElementAssemblyValues
stores per element basis values at given quadrature points and geometric mapping
Definition ElementAssemblyValues.hpp:14
polyfem::assembler::ElementAssemblyValues::compute
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,...
Definition ElementAssemblyValues.cpp:164
polyfem::assembler::ViscousDamping::compute_dstress_dpsi_dphi
static void compute_dstress_dpsi_dphi(const OptAssemblerData &data, const Eigen::MatrixXd &prev_grad_u_i, Eigen::MatrixXd &dstress_dpsi, Eigen::MatrixXd &dstress_dphi)
Definition ViscousDamping.cpp:380
polyfem::io::Evaluator::interpolate_at_local_vals
static void interpolate_at_local_vals(const mesh::Mesh &mesh, const bool is_problem_scalar, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const int el_index, const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &fun, Eigen::MatrixXd &result, Eigen::MatrixXd &result_grad)
interpolate solution and gradient at element (calls interpolate_at_local_vals with sol)
Definition Evaluator.cpp:634
polyfem::solver::ElasticForceDerivative::force_material_derivative
static void force_material_derivative(ElasticForm &form, const double t, const Eigen::MatrixXd &x, const Eigen::MatrixXd &x_prev, const Eigen::MatrixXd &adjoint, Eigen::VectorXd &term)
Definition ElasticForceDerivative.cpp:35
polyfem::solver::ElasticForceDerivative::force_shape_derivative
static void force_shape_derivative(ElasticForm &form, const double t, const int n_verts, const Eigen::MatrixXd &x, const Eigen::MatrixXd &x_prev, const Eigen::MatrixXd &adjoint, Eigen::VectorXd &term)
Definition ElasticForceDerivative.cpp:135
polyfem::assembler
Used for test only.
Definition AMIPSEnergy.cpp:6
polyfem::utils::vector2matrix
void vector2matrix(const Eigen::VectorXd &vec, Eigen::MatrixXd &mat)
Definition MatrixUtils.cpp:69
polyfem::utils::get_local_thread_storage
auto & get_local_thread_storage(Storages &storage, int thread_id)
polyfem::utils::create_thread_storage
auto create_thread_storage(const LocalStorage &initial_local_storage)
polyfem::utils::flatten
Eigen::VectorXd flatten(const Eigen::MatrixXd &X)
Flatten rowwises.
Definition MatrixUtils.cpp:34
polyfem::utils::maybe_parallel_for
void maybe_parallel_for(int size, const std::function< void(int, int, int)> &partial_for)
polyfem::QuadratureVector
Eigen::Matrix< double, Eigen::Dynamic, 1, 0, MAX_QUAD_POINTS, 1 > QuadratureVector
Definition Types.hpp:17