19 bool delta(
int i,
int j)
21 return (i == j) ? true :
false;
24 double dot(
const Eigen::MatrixXd &A,
const Eigen::MatrixXd &B) {
return (A.array() * B.array()).sum(); }
26 class LocalThreadScalarStorage
33 LocalThreadScalarStorage()
39 class LocalThreadVecStorage
46 LocalThreadVecStorage(
const int size)
73 return Eigen::VectorXd::Zero(state.
ndof());
89 j.
set_j([formulation,
this](
const Eigen::MatrixXd &local_pts,
const Eigen::MatrixXd &pts,
const Eigen::MatrixXd &u,
const Eigen::MatrixXd &grad_u,
const Eigen::VectorXd &lambda,
const Eigen::VectorXd &mu,
const Eigen::MatrixXd &reference_normals,
const assembler::ElementAssemblyValues &
vals,
const IntegrableFunctional::ParameterType ¶ms, Eigen::MatrixXd &
val) {
90 val.setZero(grad_u.rows(), 1);
91 const int dim = u.cols();
92 Eigen::MatrixXd grad_u_q;
93 for (
int q = 0; q < grad_u.rows(); q++)
96 if (formulation ==
"LinearElasticity")
98 grad_u_q = (grad_u_q + grad_u_q.transpose()).eval() / 2.;
99 val(q) = mu(q) * (grad_u_q.transpose() * grad_u_q).trace() + lambda(q) / 2. * grad_u_q.trace() * grad_u_q.trace();
101 else if (formulation ==
"NeoHookean")
103 Eigen::MatrixXd def_grad = grad_u_q + Eigen::MatrixXd::Identity(dim, dim);
104 double log_det_j = log(def_grad.determinant());
105 val(q) = mu(q) / 2 * ((def_grad.transpose() * def_grad).trace() - dim - 2 * log_det_j) + lambda(q) / 2 * log_det_j * log_det_j;
112 j.
set_dj_dgradu([formulation,
this](
const Eigen::MatrixXd &local_pts,
const Eigen::MatrixXd &pts,
const Eigen::MatrixXd &u,
const Eigen::MatrixXd &grad_u,
const Eigen::VectorXd &lambda,
const Eigen::VectorXd &mu,
const Eigen::MatrixXd &reference_normals,
const assembler::ElementAssemblyValues &
vals,
const IntegrableFunctional::ParameterType ¶ms, Eigen::MatrixXd &
val) {
113 val.setZero(grad_u.rows(), grad_u.cols());
114 Eigen::MatrixXd grad_u_q, def_grad, FmT, stress;
115 for (
int q = 0; q < grad_u.rows(); q++)
118 if (formulation ==
"LinearElasticity")
120 stress = mu(q) * (grad_u_q + grad_u_q.transpose()) + lambda(q) * grad_u_q.trace() * Eigen::MatrixXd::Identity(grad_u_q.rows(), grad_u_q.cols());
122 else if (formulation ==
"NeoHookean")
124 def_grad = Eigen::MatrixXd::Identity(grad_u_q.rows(), grad_u_q.cols()) + grad_u_q;
125 FmT = def_grad.inverse().transpose();
126 stress = mu(q) * (def_grad - FmT) + lambda(q) * std::log(def_grad.determinant()) * FmT;
142 return Eigen::VectorXd::Zero(0).eval();
153 j.
set_j([formulation, power, &state = std::as_const(
state_)](
const Eigen::MatrixXd &local_pts,
const Eigen::MatrixXd &pts,
const Eigen::MatrixXd &u,
const Eigen::MatrixXd &grad_u,
const Eigen::VectorXd &lambda,
const Eigen::VectorXd &mu,
const Eigen::MatrixXd &reference_normals,
const assembler::ElementAssemblyValues &
vals,
const IntegrableFunctional::ParameterType ¶ms, Eigen::MatrixXd &
val) {
154 val.setZero(grad_u.rows(), 1);
155 const double dt = state.problem->is_time_dependent() ? state.args[
"time"][
"dt"].get<
double>() : 0;
157 Eigen::MatrixXd grad_u_q, stress, grad_unused;
158 for (
int q = 0; q < grad_u.rows(); q++)
160 if (formulation ==
"Laplacian")
161 stress = grad_u.row(q);
165 state.assembler->compute_stress_grad_multiply_mat(
OptAssemblerData(params.
t, dt, params.
elem, local_pts.row(q), pts.row(q), grad_u_q), Eigen::MatrixXd::Zero(grad_u_q.rows(), grad_u_q.cols()), stress, grad_unused);
167 val(q) = pow(stress.squaredNorm(), power / 2.);
171 j.
set_dj_dgradu([formulation, power, &state = std::as_const(
state_)](
const Eigen::MatrixXd &local_pts,
const Eigen::MatrixXd &pts,
const Eigen::MatrixXd &u,
const Eigen::MatrixXd &grad_u,
const Eigen::VectorXd &lambda,
const Eigen::VectorXd &mu,
const Eigen::MatrixXd &reference_normals,
const assembler::ElementAssemblyValues &
vals,
const IntegrableFunctional::ParameterType ¶ms, Eigen::MatrixXd &
val) {
172 val.setZero(grad_u.rows(), grad_u.cols());
173 const double dt = state.problem->is_time_dependent() ? state.args[
"time"][
"dt"].get<
double>() : 0;
174 const int dim = state.mesh->dimension();
176 if (formulation ==
"Laplacian")
178 for (
int q = 0; q < grad_u.rows(); q++)
180 const double coef = power * pow(grad_u.row(q).squaredNorm(), power / 2. - 1.);
181 val.row(q) = coef * grad_u.row(q);
186 Eigen::MatrixXd grad_u_q, stress, stress_dstress;
187 for (
int q = 0; q < grad_u.rows(); q++)
190 state.assembler->compute_stress_grad_multiply_stress(
OptAssemblerData(params.
t, dt, params.
elem, local_pts.row(q), pts.row(q), grad_u_q), stress, stress_dstress);
192 const double coef = power * pow(stress.squaredNorm(), power / 2. - 1.);
206 return Eigen::VectorXd::Zero(0).eval();
218 val.setZero(grad_u.rows(), 1);
220 Eigen::MatrixXd grad_u_q, stress;
221 for (
int q = 0; q < grad_u.rows(); q++)
224 stress = mu(q) * (grad_u_q + grad_u_q.transpose()) + lambda(q) * grad_u_q.trace() * Eigen::MatrixXd::Identity(grad_u_q.rows(), grad_u_q.cols());
225 val(q) = (stress.array() * grad_u_q.array()).sum();
230 val.setZero(grad_u.rows(), grad_u.cols());
231 const int dim = u.cols();
233 Eigen::MatrixXd grad_u_q, stress;
234 for (
int q = 0; q < grad_u.rows(); q++)
237 stress = mu(q) * (grad_u_q + grad_u_q.transpose()) + lambda(q) * grad_u_q.trace() * Eigen::MatrixXd::Identity(grad_u_q.rows(), grad_u_q.cols());
248 const double t =
state_.
problem->is_time_dependent() ? dt * time_step +
state_.
args[
"time"][
"t0"].get<
double>() : 0;
253 Eigen::VectorXd term = Eigen::VectorXd::Zero(bases.size() * 2);
256 for (
int e = 0; e < bases.size(); e++)
264 Eigen::MatrixXd u, grad_u;
267 Eigen::MatrixXd grad_u_q;
276 Eigen::MatrixXd f_prime_dmu, f_prime_dlambda;
279 term(e + bases.size()) += dot(f_prime_dmu, grad_u_q) *
da(q);
280 term(e) += dot(f_prime_dlambda, grad_u_q) *
da(q);
290 const int dim =
dim_;
293 val = u.col(dim) + pts.col(dim);
297 val.setZero(u.rows(), u.cols());
298 val.col(dim).setOnes();
302 val.setZero(pts.rows(), pts.cols());
303 val.col(dim).setOnes();
312 const int dim = this->
dim_;
314 j.
set_j([dim, &state = std::as_const(this->
state_)](
const Eigen::MatrixXd &local_pts,
const Eigen::MatrixXd &pts,
const Eigen::MatrixXd &u,
const Eigen::MatrixXd &grad_u,
const Eigen::VectorXd &lambda,
const Eigen::VectorXd &mu,
const Eigen::MatrixXd &reference_normals,
const assembler::ElementAssemblyValues &
vals,
const IntegrableFunctional::ParameterType ¶ms, Eigen::MatrixXd &
val) {
315 Eigen::MatrixXd acc, grad_acc;
321 j.
set_dj_du([dim,
this](
const Eigen::MatrixXd &local_pts,
const Eigen::MatrixXd &pts,
const Eigen::MatrixXd &u,
const Eigen::MatrixXd &grad_u,
const Eigen::VectorXd &lambda,
const Eigen::VectorXd &mu,
const Eigen::MatrixXd &reference_normals,
const assembler::ElementAssemblyValues &
vals,
const IntegrableFunctional::ParameterType ¶ms, Eigen::MatrixXd &
val) {
322 val.setZero(u.rows(), u.cols());
334 Eigen::MatrixXd v, grad_v;
337 val.setZero(u.rows(), 1);
338 for (
int q = 0; q < v.rows(); q++)
341 val(q) = 0.5 * rho * v.row(q).squaredNorm();
346 val.setZero(u.rows(), u.cols());
359 return Eigen::VectorXd::Zero(0).eval();
370 j.
set_j([formulation, dimensions,
this](
const Eigen::MatrixXd &local_pts,
const Eigen::MatrixXd &pts,
const Eigen::MatrixXd &u,
const Eigen::MatrixXd &grad_u,
const Eigen::VectorXd &lambda,
const Eigen::VectorXd &mu,
const Eigen::MatrixXd &reference_normals,
const assembler::ElementAssemblyValues &
vals,
const IntegrableFunctional::ParameterType ¶ms, Eigen::MatrixXd &
val) {
371 val.setZero(grad_u.rows(), 1);
373 Eigen::MatrixXd grad_u_q, stress;
374 for (
int q = 0; q < grad_u.rows(); q++)
377 if (formulation ==
"LinearElasticity")
379 stress = mu(q) * (grad_u_q + grad_u_q.transpose()) + lambda(q) * grad_u_q.trace() * Eigen::MatrixXd::Identity(grad_u_q.rows(), grad_u_q.cols());
381 else if (formulation ==
"NeoHookean")
383 Eigen::MatrixXd def_grad = Eigen::MatrixXd::Identity(grad_u_q.rows(), grad_u_q.cols()) + grad_u_q;
384 Eigen::MatrixXd FmT = def_grad.inverse().transpose();
385 stress = mu(q) * (def_grad - FmT) + lambda(q) * std::log(def_grad.determinant()) * FmT;
389 val(q) = stress(dimensions[0], dimensions[1]);
393 j.
set_dj_dgradu([formulation, dimensions,
this](
const Eigen::MatrixXd &local_pts,
const Eigen::MatrixXd &pts,
const Eigen::MatrixXd &u,
const Eigen::MatrixXd &grad_u,
const Eigen::VectorXd &lambda,
const Eigen::VectorXd &mu,
const Eigen::MatrixXd &reference_normals,
const assembler::ElementAssemblyValues &
vals,
const IntegrableFunctional::ParameterType ¶ms, Eigen::MatrixXd &
val) {
394 val.setZero(grad_u.rows(), grad_u.cols());
397 Eigen::MatrixXd grad_u_q, stiffness, stress;
398 for (
int q = 0; q < grad_u.rows(); q++)
400 stiffness.setZero(1, dim * dim * dim * dim);
403 if (formulation ==
"LinearElasticity")
405 stress = mu(q) * (grad_u_q + grad_u_q.transpose()) + lambda(q) * grad_u_q.trace() * Eigen::MatrixXd::Identity(grad_u_q.rows(), grad_u_q.cols());
406 for (
int i = 0, idx = 0; i < dim; i++)
407 for (
int j = 0; j < dim; j++)
408 for (
int k = 0; k < dim; k++)
409 for (
int l = 0; l < dim; l++)
411 stiffness(idx++) = mu(q) * delta(i, k) * delta(j, l) + mu(q) * delta(i, l) * delta(j, k) + lambda(q) * delta(i, j) * delta(k, l);
414 else if (formulation ==
"NeoHookean")
416 Eigen::MatrixXd def_grad = Eigen::MatrixXd::Identity(grad_u_q.rows(), grad_u_q.cols()) + grad_u_q;
417 Eigen::MatrixXd FmT = def_grad.inverse().transpose();
418 stress = mu(q) * (def_grad - FmT) + lambda(q) * std::log(def_grad.determinant()) * FmT;
420 double J = def_grad.determinant();
421 double tmp1 = mu(q) - lambda(q) * std::log(J);
422 for (
int i = 0, idx = 0; i < dim; i++)
423 for (
int j = 0; j < dim; j++)
424 for (
int k = 0; k < dim; k++)
425 for (
int l = 0; l < dim; l++)
427 stiffness(idx++) = mu(q) * delta(i, k) * delta(j, l) + tmp1 * FmT(i, l) * FmT(k, j);
429 stiffness += lambda(q) *
utils::flatten(FmT_vec * FmT_vec.transpose()).transpose();
434 val.row(q) = stiffness.block(0, (dimensions[0] * dim + dimensions[1]) * dim * dim, 1, dim * dim);
446 val.setOnes(grad_u.rows(), 1);
ElementAssemblyValues vals
void set_dj_du(const functionalType &dj_du)
void set_dj_dgradu(const functionalType &dj_dgradu)
void set_j(const functionalType &j)
void set_dj_dx(const functionalType &dj_dx)
main class that contains the polyfem solver and all its state
std::string formulation() const
return the formulation (checks if the problem is scalar or not and deals with multiphysics)
assembler::AssemblyValsCache ass_vals_cache
used to store assembly values for small problems
const std::vector< basis::ElementBases > & geom_bases() const
Get a constant reference to the geometry mapping bases.
std::shared_ptr< assembler::Assembler > assembler
assemblers
std::shared_ptr< assembler::Mass > mass_matrix_assembler
std::unique_ptr< mesh::Mesh > mesh
current mesh, it can be a Mesh2D or Mesh3D
std::shared_ptr< assembler::Problem > problem
current problem, it contains rhs and bc
json args
main input arguments containing all defaults
solver::DiffCache diff_cached
std::vector< basis::ElementBases > bases
FE bases, the size is #elements.
void compute(const int el_index, const bool is_volume, const basis::ElementBases &basis, const basis::ElementBases &gbasis, ElementAssemblyValues &vals) const
retrieves cached basis evaluation and geometric for the given element if it doesn't exist,...
stores per element basis values at given quadrature points and geometric mapping
quadrature::Quadrature quadrature
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)
Eigen::VectorXd v(int step) const
Eigen::VectorXd u(int step) const
Eigen::VectorXd apply_parametrization_jacobian(const ParameterType type, const State *state_ptr, const Eigen::VectorXd &x, const std::function< Eigen::VectorXd()> &grad) const
Maps the partial gradient wrt.
void vector2matrix(const Eigen::VectorXd &vec, Eigen::MatrixXd &mat)
Eigen::MatrixXd unflatten(const Eigen::VectorXd &x, int dim)
Unflatten rowwises, so every dim elements in x become a row.
Eigen::VectorXd flatten(const Eigen::MatrixXd &X)
Flatten rowwises.
Eigen::Matrix< double, Eigen::Dynamic, 1, 0, MAX_QUAD_POINTS, 1 > QuadratureVector
void log_and_throw_adjoint_error(const std::string &msg)
Parameters for the functional evaluation.