17 return std::make_unique<ShapeVariableToSimulation>(states, parametrization);
18 else if (type ==
"elastic")
19 return std::make_unique<ElasticVariableToSimulation>(states, parametrization);
20 else if (type ==
"friction")
21 return std::make_unique<FrictionCoeffientVariableToSimulation>(states, parametrization);
22 else if (type ==
"damping")
23 return std::make_unique<DampingCoeffientVariableToSimulation>(states, parametrization);
24 else if (type ==
"initial")
25 return std::make_unique<InitialConditionVariableToSimulation>(states, parametrization);
26 else if (type ==
"dirichlet")
27 return std::make_unique<DirichletVariableToSimulation>(states, parametrization);
28 else if (type ==
"dirichlet-nodes")
29 return std::make_unique<DirichletNodesVariableToSimulation>(states, parametrization);
30 else if (type ==
"pressure")
31 return std::make_unique<PressureVariableToSimulation>(states, parametrization);
32 else if (type ==
"periodic-shape")
33 return std::make_unique<PeriodicShapeVariableToSimulation>(states, parametrization);
36 return std::unique_ptr<VariableToSimulation>();
41 const std::string composite_map_type = args[
"composite_map_type"];
43 if (composite_map_type ==
"none")
47 else if (composite_map_type ==
"indices")
49 if (args[
"composite_map_indices"].is_string())
55 else if (args[
"composite_map_indices"].is_array())
72 ind.setLinSpaced(out_size, 0, out_size - 1);
77 return Eigen::VectorXi();
88 return Eigen::VectorXd();
98 std::vector<ValueType>().swap(
L);
99 for (
const auto &arg : args)
105 Eigen::VectorXd adjoint_term = Eigen::VectorXd::Zero(
x.size());
106 for (
const auto &v2s :
L)
107 adjoint_term += v2s->compute_adjoint_term(
x);
113 for (
const auto &v2s :
L)
115 if (v2s->get_parameter_type() != type)
118 const Eigen::VectorXd var = v2s->get_parametrization().eval(
x);
119 for (
const auto &state : v2s->get_states())
121 if (state.get() != state_ptr)
124 state_variable(v2s->get_output_indexing(
x)) = var;
131 Eigen::VectorXd gradv = Eigen::VectorXd::Zero(
x.size());
132 for (
const auto &v2s :
L)
134 if (v2s->get_parameter_type() != type)
137 for (
const auto &state : v2s->get_states())
139 if (state.get() != state_ptr)
142 gradv += v2s->apply_parametrization_jacobian(grad(),
x);
152 const int dim = state->mesh->dimension();
155 for (
int i = 0; i < indices.size(); i += dim)
156 for (
int j = 0; j < dim; j++)
157 assert(indices(i + j) == indices(i) + j);
159 for (
int i = 0; i < indices.size(); i += dim)
160 state->set_mesh_vertex(indices(i) / dim, state_variable(Eigen::seqN(i, dim)));
165 Eigen::VectorXd term, cur_term;
168 if (state->problem->is_time_dependent())
172 if (!state->is_homogenization())
178 if (term.size() != cur_term.size())
187 const int dim =
states_[0]->mesh->dimension();
188 const int npts =
states_[0]->mesh->n_vertices();
193 if (indices.size() == 0)
194 indices.setLinSpaced(npts * dim, 0, npts * dim - 1);
198 if (indices.maxCoeff() >=
V.size())
206 const std::string composite_map_type = args[
"composite_map_type"];
208 if (composite_map_type ==
"interior")
213 else if (composite_map_type ==
"boundary")
218 else if (composite_map_type ==
"boundary_excluding_surface")
220 const std::vector<int> excluded_surfaces = args[
"surface_selection"];
232 const int n_elem = state->bases.size();
233 assert(n_elem * 2 == state_variable.size());
234 state->assembler->update_lame_params(state_variable.segment(0, n_elem), state_variable.segment(n_elem, n_elem));
239 Eigen::VectorXd term, cur_term;
242 if (state->problem->is_time_dependent())
247 if (term.size() != cur_term.size())
257 auto params_map = state.assembler->parameters();
259 auto search_lambda = params_map.find(
"lambda");
260 auto search_mu = params_map.find(
"mu");
261 if (search_lambda == params_map.end() || search_mu == params_map.end())
264 return Eigen::VectorXd();
267 Eigen::VectorXd lambdas(state.mesh->n_elements());
268 Eigen::VectorXd mus(state.mesh->n_elements());
269 for (
int e = 0; e < state.mesh->n_elements(); e++)
272 if (!state.mesh->is_volume())
274 const auto &mesh2d = *
dynamic_cast<mesh::Mesh2D *
>(state.mesh.get());
279 const auto &mesh3d = *
dynamic_cast<mesh::Mesh3D *
>(state.mesh.get());
282 lambdas(e) = search_lambda->second(RowVectorNd::Zero(state.mesh->dimension()), barycenter, 0., e);
283 mus(e) = search_mu->second(RowVectorNd::Zero(state.mesh->dimension()), barycenter, 0., e);
285 state.assembler->update_lame_params(lambdas, mus);
287 Eigen::VectorXd params(lambdas.size() + mus.size());
288 params << lambdas, mus;
295 assert(state_variable.size() == 1);
296 assert(state_variable(0) >= 0);
298 state->args[
"contact"][
"friction_coefficient"] = state_variable(0);
302 Eigen::VectorXd term, cur_term;
305 if (state->problem->is_time_dependent())
310 if (term.size() != cur_term.size())
320 return Eigen::VectorXd();
325 assert(state_variable.size() == 2);
326 json damping_param = {
327 {
"psi", state_variable(0)},
328 {
"phi", state_variable(1)},
332 if (!state->args[
"materials"].is_array())
334 state->args[
"materials"][
"psi"] = damping_param[
"psi"];
335 state->args[
"materials"][
"phi"] = damping_param[
"phi"];
339 for (
auto &arg : state->args[
"materials"])
341 arg[
"psi"] = damping_param[
"psi"];
342 arg[
"phi"] = damping_param[
"phi"];
346 if (state->damping_assembler)
347 state->damping_assembler->add_multimaterial(0, damping_param, state->units);
349 logger().info(
"[{}] Current params: {}, {}",
name(), state_variable(0), state_variable(1));
353 Eigen::VectorXd term, cur_term;
356 if (state->problem->is_time_dependent())
361 if (term.size() != cur_term.size())
371 return Eigen::VectorXd();
378 if (state_variable.size() != state->ndof() * 2)
380 log_and_throw_adjoint_error(
"[{}] Inconsistent number of parameters {} and number of dofs in forward {}!",
name(), state_variable.size(), state->ndof() * 2);
382 state->initial_sol_update = state_variable.head(state->ndof());
383 state->initial_vel_update = state_variable.tail(state->ndof());
388 Eigen::VectorXd term, cur_term;
391 if (state->problem->is_time_dependent())
396 if (term.size() != cur_term.size())
406 Eigen::MatrixXd sol, vel;
407 state.initial_solution(sol);
408 state.initial_velocity(vel);
410 Eigen::VectorXd
x(sol.size() + vel.size());
417 auto tensor_problem = std::dynamic_pointer_cast<polyfem::assembler::GenericTensorProblem>(
states_[0]->problem);
419 int dim =
states_[0]->mesh->dimension();
420 int num_steps = indices.size() / dim;
421 for (
int i = 0; i < num_steps; ++i)
423 tensor_problem->update_dirichlet_boundary(b, indices(i * dim) + 1, state_variable.segment(i * dim, dim));
430 Eigen::VectorXd term, cur_term;
433 if (state->problem->is_time_dependent())
438 if (term.size() != cur_term.size())
454 int dim =
states_[0]->mesh->dimension();
456 for (
const auto &b :
states_[0]->args[
"boundary_conditions"][
"dirichlet_boundary"])
459 auto value = b[
"value"];
460 if (value[0].is_array())
462 if (!
states_[0]->problem->is_time_dependent())
464 Eigen::MatrixXd dirichlet = value;
465 x.setZero(dirichlet.rows() * (dirichlet.cols() - 1));
466 for (
int j = 1; j < dirichlet.cols(); ++j)
467 x.segment((j - 1) * dim, dim) = dirichlet.col(j);
469 else if (value[0].is_number())
471 if (
states_[0]->problem->is_time_dependent())
476 else if (value.is_string())
485 const std::string composite_map_type = args[
"composite_map_type"];
487 if (composite_map_type ==
"time_step_indexing")
489 const int time_steps = state.
args[
"time"][
"time_steps"];
490 const int dim = state.
mesh->dimension();
493 for (
int i = 0; i < time_steps; ++i)
494 for (
int k = 0; k < dim; ++k)
507 assert(state_variable.size() == (state->mesh->dimension() *
dirichlet_nodes_.size()));
508 auto tensor_problem = std::dynamic_pointer_cast<polyfem::assembler::GenericTensorProblem>(state->problem);
509 assert(!state->problem->is_time_dependent());
512 logger().info(
"Updated dirichlet nodes");
518 Eigen::VectorXd term, cur_term;
521 if (state->problem->is_time_dependent())
526 if (term.size() != cur_term.size())
544 const std::string composite_map_type = args_[
"composite_map_type"];
545 if (composite_map_type !=
"indices")
553 std::vector<int> composite_map_indices = {};
555 for (
int k = 0; k < dim; ++k)
557 args_[
"composite_map_indices"] = composite_map_indices;
564 auto tensor_problem = std::dynamic_pointer_cast<polyfem::assembler::GenericTensorProblem>(
states_[0]->problem);
566 for (
int i = 0; i < indices.size(); ++i)
568 tensor_problem->update_pressure_boundary(b, indices(i) + 1, state_variable(i));
575 Eigen::VectorXd term, cur_term;
578 if (state->problem->is_time_dependent())
580 Eigen::MatrixXd adjoint_nu, adjoint_p;
581 adjoint_nu = state->get_adjoint_mat(1);
582 adjoint_p = state->get_adjoint_mat(0);
589 if (term.size() != cur_term.size())
608 for (
const auto &b :
states_[0]->args[
"boundary_conditions"][
"pressure_boundary"])
611 auto value = b[
"value"];
612 if (value.is_array())
614 if (!
states_[0]->problem->is_time_dependent())
616 Eigen::VectorXd pressures = value;
617 x = pressures.segment(1, pressures.size() - 1);
619 else if (value.is_number())
621 if (
states_[0]->problem->is_time_dependent())
626 else if (value.is_string())
636 const std::string composite_map_type = args[
"composite_map_type"];
638 if (composite_map_type ==
"time_step_indexing")
640 const int time_steps = state.
args[
"time"][
"time_steps"];
642 for (
int i = 0; i < time_steps; ++i)
653 Eigen::VectorXd term, cur_term;
656 if (state->problem->is_time_dependent())
664 if (term.size() != cur_term.size())
673 const int dim =
states_[0]->mesh->dimension();
679 const int n_verts = state->mesh->n_vertices();
681 for (
int i = 0; i < n_verts; i++)
682 state->set_mesh_vertex(i,
V.row(i));
687 const auto &state = *(
states_[0]);
690 state.get_vertices(
V);
692 if (!state.periodic_bc->all_direction_periodic())
main class that contains the polyfem solver and all its state
std::unique_ptr< mesh::Mesh > mesh
current mesh, it can be a Mesh2D or Mesh3D
json args
main input arguments containing all defaults
std::string resolve_input_path(const std::string &path, const bool only_if_exists=false) const
Resolve input path relative to root_path() if the path is not absolute.
RowVectorNd face_barycenter(const int index) const override
face barycenter
virtual RowVectorNd cell_barycenter(const int c) const =0
cell barycenter
int size(const int x_size) const override
Eigen::VectorXd apply_jacobian(const Eigen::VectorXd &grad_full, const Eigen::VectorXd &x) const override
Eigen::VectorXd inverse_eval(const Eigen::VectorXd &y) override
Eigen::VectorXd eval(const Eigen::VectorXd &x) const override
Eigen::VectorXd compute_adjoint_term(const Eigen::VectorXd &x) const override
Eigen::VectorXd inverse_eval() override
void update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices) override
std::string name() const override
Eigen::VectorXd inverse_eval() override
void update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices) override
std::string name() const override
Eigen::VectorXi dirichlet_nodes_
void set_dirichlet_nodes(const Eigen::VectorXi &dirichlet_nodes)
std::string variable_to_string(const Eigen::VectorXd &variable)
Eigen::VectorXd compute_adjoint_term(const Eigen::VectorXd &x) const override
void set_output_indexing(const json &args) override
Eigen::VectorXd compute_adjoint_term(const Eigen::VectorXd &x) const override
std::string variable_to_string(const Eigen::VectorXd &variable)
std::vector< int > dirichlet_boundaries_
std::string name() const override
void update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices) override
Eigen::VectorXd inverse_eval() override
void set_dirichlet_boundaries(const std::vector< int > &dirichlet_boundaries)
void set_output_indexing(const json &args) override
std::string name() const override
Eigen::VectorXd compute_adjoint_term(const Eigen::VectorXd &x) const override
void update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices) override
Eigen::VectorXd inverse_eval() override
Eigen::VectorXd compute_adjoint_term(const Eigen::VectorXd &x) const override
std::string name() const override
Eigen::VectorXd inverse_eval() override
void update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices) override
Eigen::VectorXd inverse_eval() override
Eigen::VectorXd compute_adjoint_term(const Eigen::VectorXd &x) const override
std::string name() const override
void update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices) override
void update(const Eigen::VectorXd &x) override
Eigen::VectorXd compute_adjoint_term(const Eigen::VectorXd &x) const override
Eigen::VectorXd apply_parametrization_jacobian(const Eigen::VectorXd &term, const Eigen::VectorXd &x) const override
Eigen::VectorXd periodic_mesh_representation
Eigen::VectorXd inverse_eval() override
std::unique_ptr< PeriodicMeshToMesh > periodic_mesh_map
std::vector< int > pressure_boundaries_
void set_pressure_boundaries(const std::vector< int > &pressure_boundaries)
Eigen::VectorXd inverse_eval() override
std::string variable_to_string(const Eigen::VectorXd &variable)
void update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices) override
void set_output_indexing(const json &args) override
Eigen::VectorXd compute_adjoint_term(const Eigen::VectorXd &x) const override
virtual void update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices) override
void set_output_indexing(const json &args) override
Eigen::VectorXd compute_adjoint_term(const Eigen::VectorXd &x) const override
Eigen::VectorXd inverse_eval() override
const Eigen::VectorXi & get_output_indexing() const
void compute_state_variable(const ParameterType type, const State *state_ptr, const Eigen::VectorXd &x, Eigen::VectorXd &state_variable) const
Evaluate the variable to simulations and overwrite the state_variable based on x.
Eigen::VectorXd compute_adjoint_term(const Eigen::VectorXd &x) const
Computes the sum of adjoint terms for all VariableToSimulation.
virtual 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 init(const json &args, const std::vector< std::shared_ptr< State > > &states, const std::vector< int > &variable_sizes)
std::vector< ValueType > L
virtual Eigen::VectorXd apply_parametrization_jacobian(const Eigen::VectorXd &term, const Eigen::VectorXd &x) const
CompositeParametrization parametrization_
virtual Eigen::VectorXd inverse_eval()
const std::vector< std::shared_ptr< State > > states_
virtual void update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices)
virtual void set_output_indexing(const json &args)
virtual std::string name() const =0
Eigen::VectorXi get_output_indexing(const Eigen::VectorXd &x) const
static std::unique_ptr< VariableToSimulation > create(const std::string &type, const std::vector< std::shared_ptr< State > > &states, CompositeParametrization &¶metrization)
Eigen::VectorXi output_indexing_
bool read_matrix(const std::string &path, Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &mat)
Reads a matrix from a file. Determines the file format based on the path's extension.
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.
spdlog::logger & logger()
Retrieves the current logger.
void log_and_throw_adjoint_error(const std::string &msg)
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > RowVectorNd
void log_and_throw_error(const std::string &msg)
static std::unique_ptr< VariableToSimulation > create_variable_to_simulation(const json &args, const std::vector< std::shared_ptr< State > > &states, const std::vector< int > &variable_sizes)