Loading [MathJax]/extensions/tex2jax.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
VariableToSimulation.cpp
Go to the documentation of this file.
2#include <polyfem/State.hpp>
6
8
11
12namespace polyfem::solver
13{
14 std::unique_ptr<VariableToSimulation> VariableToSimulation::create(const std::string &type, const std::vector<std::shared_ptr<State>> &states, CompositeParametrization &&parametrization)
15 {
16 if (type == "shape")
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);
34
35 log_and_throw_adjoint_error("Invalid type of VariableToSimulation!");
36 return std::unique_ptr<VariableToSimulation>();
37 }
38
40 {
41 const std::string composite_map_type = args["composite_map_type"];
42 const State &state = *(states_[0]);
43 if (composite_map_type == "none")
44 {
45 output_indexing_.resize(0);
46 }
47 else if (composite_map_type == "indices")
48 {
49 if (args["composite_map_indices"].is_string())
50 {
51 Eigen::MatrixXi tmp_mat;
52 polyfem::io::read_matrix(state.resolve_input_path(args["composite_map_indices"].get<std::string>()), tmp_mat);
54 }
55 else if (args["composite_map_indices"].is_array())
56 output_indexing_ = args["composite_map_indices"];
57 else
58 log_and_throw_adjoint_error("Invalid composite map indices type!");
59 }
60 else
61 log_and_throw_adjoint_error("Unknown composite_map_type!");
62 }
63
64 Eigen::VectorXi VariableToSimulation::get_output_indexing(const Eigen::VectorXd &x) const
65 {
66 const int out_size = parametrization_.size(x.size());
67 if (output_indexing_.size() == out_size || out_size == 0)
68 return output_indexing_;
69 else if (output_indexing_.size() == 0)
70 {
71 Eigen::VectorXi ind;
72 ind.setLinSpaced(out_size, 0, out_size - 1);
73 return ind;
74 }
75 else
76 log_and_throw_adjoint_error(fmt::format("[{}] Indexing size and output size of the Parametrization do not match! {} vs {}", name(), output_indexing_.size(), out_size));
77 return Eigen::VectorXi();
78 }
79
80 Eigen::VectorXd VariableToSimulation::apply_parametrization_jacobian(const Eigen::VectorXd &term, const Eigen::VectorXd &x) const
81 {
83 }
84
86 {
87 log_and_throw_adjoint_error("[{}] inverse_eval not implemented!", name());
88 return Eigen::VectorXd();
89 }
90
91 void VariableToSimulation::update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices)
92 {
93 log_and_throw_adjoint_error("[{}] update_state not implemented!", name());
94 }
95
96 void VariableToSimulationGroup::init(const json &args, const std::vector<std::shared_ptr<State>> &states, const std::vector<int> &variable_sizes)
97 {
98 std::vector<ValueType>().swap(L);
99 for (const auto &arg : args)
100 L.push_back(AdjointOptUtils::create_variable_to_simulation(arg, states, variable_sizes));
101 }
102
103 Eigen::VectorXd VariableToSimulationGroup::compute_adjoint_term(const Eigen::VectorXd &x) const
104 {
105 Eigen::VectorXd adjoint_term = Eigen::VectorXd::Zero(x.size());
106 for (const auto &v2s : L)
107 adjoint_term += v2s->compute_adjoint_term(x);
108 return adjoint_term;
109 }
110
111 void VariableToSimulationGroup::compute_state_variable(const ParameterType type, const State *state_ptr, const Eigen::VectorXd &x, Eigen::VectorXd &state_variable) const
112 {
113 for (const auto &v2s : L)
114 {
115 if (v2s->get_parameter_type() != type)
116 continue;
117
118 const Eigen::VectorXd var = v2s->get_parametrization().eval(x);
119 for (const auto &state : v2s->get_states())
120 {
121 if (state.get() != state_ptr)
122 continue;
123
124 state_variable(v2s->get_output_indexing(x)) = var;
125 }
126 }
127 }
128
129 Eigen::VectorXd VariableToSimulationGroup::apply_parametrization_jacobian(const ParameterType type, const State *state_ptr, const Eigen::VectorXd &x, const std::function<Eigen::VectorXd()> &grad) const
130 {
131 Eigen::VectorXd gradv = Eigen::VectorXd::Zero(x.size());
132 for (const auto &v2s : L)
133 {
134 if (v2s->get_parameter_type() != type)
135 continue;
136
137 for (const auto &state : v2s->get_states())
138 {
139 if (state.get() != state_ptr)
140 continue;
141
142 gradv += v2s->apply_parametrization_jacobian(grad(), x);
143 }
144 }
145 return gradv;
146 }
147
148 void ShapeVariableToSimulation::update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices)
149 {
150 for (auto state : states_)
151 {
152 const int dim = state->mesh->dimension();
153
154 // If indices include one vertex entry, we assume it include all entries of this vertex.
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);
158
159 for (int i = 0; i < indices.size(); i += dim)
160 state->set_mesh_vertex(indices(i) / dim, state_variable(Eigen::seqN(i, dim)));
161 }
162 }
163 Eigen::VectorXd ShapeVariableToSimulation::compute_adjoint_term(const Eigen::VectorXd &x) const
164 {
165 Eigen::VectorXd term, cur_term;
166 for (auto state : states_)
167 {
168 if (state->problem->is_time_dependent())
169 AdjointTools::dJ_shape_transient_adjoint_term(*state, state->get_adjoint_mat(1), state->get_adjoint_mat(0), cur_term);
170 else
171 {
172 if (!state->is_homogenization())
173 AdjointTools::dJ_shape_static_adjoint_term(*state, state->diff_cached.u(0), state->get_adjoint_mat(0), cur_term);
174 else
175 AdjointTools::dJ_shape_homogenization_adjoint_term(*state, state->diff_cached.u(0), state->get_adjoint_mat(0), cur_term);
176 }
177
178 if (term.size() != cur_term.size())
179 term = cur_term;
180 else
181 term += cur_term;
182 }
183 return apply_parametrization_jacobian(term, x);
184 }
186 {
187 const int dim = states_[0]->mesh->dimension();
188 const int npts = states_[0]->mesh->n_vertices();
189
190 Eigen::VectorXd x;
191 Eigen::VectorXi indices = get_output_indexing(x);
192
193 if (indices.size() == 0)
194 indices.setLinSpaced(npts * dim, 0, npts * dim - 1);
195
196 Eigen::MatrixXd V;
197 states_[0]->get_vertices(V);
198 if (indices.maxCoeff() >= V.size())
199 log_and_throw_adjoint_error("Output indices larger than DoFs of vertices!");
200 x = utils::flatten(V)(indices);
201
203 }
205 {
206 const std::string composite_map_type = args["composite_map_type"];
207 const State &state = *(states_[0]);
208 if (composite_map_type == "interior")
209 {
210 VariableToInteriorNodes variable_to_node(state, args["volume_selection"]);
211 output_indexing_ = variable_to_node.get_output_indexing();
212 }
213 else if (composite_map_type == "boundary")
214 {
215 VariableToBoundaryNodes variable_to_node(state, args["surface_selection"]);
216 output_indexing_ = variable_to_node.get_output_indexing();
217 }
218 else if (composite_map_type == "boundary_excluding_surface")
219 {
220 const std::vector<int> excluded_surfaces = args["surface_selection"];
221 VariableToBoundaryNodesExclusive variable_to_node(state, excluded_surfaces);
222 output_indexing_ = variable_to_node.get_output_indexing();
223 }
224 else
226 }
227
228 void ElasticVariableToSimulation::update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices)
229 {
230 for (auto state : states_)
231 {
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));
235 }
236 }
237 Eigen::VectorXd ElasticVariableToSimulation::compute_adjoint_term(const Eigen::VectorXd &x) const
238 {
239 Eigen::VectorXd term, cur_term;
240 for (auto state : states_)
241 {
242 if (state->problem->is_time_dependent())
243 AdjointTools::dJ_material_transient_adjoint_term(*state, state->get_adjoint_mat(1), state->get_adjoint_mat(0), cur_term);
244 else
245 AdjointTools::dJ_material_static_adjoint_term(*state, state->diff_cached.u(0), state->get_adjoint_mat(0), cur_term);
246
247 if (term.size() != cur_term.size())
248 term = cur_term;
249 else
250 term += cur_term;
251 }
252 return apply_parametrization_jacobian(term, x);
253 }
255 {
256 auto &state = *(states_[0]);
257 auto params_map = state.assembler->parameters();
258
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())
262 {
263 log_and_throw_adjoint_error("[{}] Failed to find Lame parameters!", name());
264 return Eigen::VectorXd();
265 }
266
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++)
270 {
271 RowVectorNd barycenter;
272 if (!state.mesh->is_volume())
273 {
274 const auto &mesh2d = *dynamic_cast<mesh::Mesh2D *>(state.mesh.get());
275 barycenter = mesh2d.face_barycenter(e);
276 }
277 else
278 {
279 const auto &mesh3d = *dynamic_cast<mesh::Mesh3D *>(state.mesh.get());
280 barycenter = mesh3d.cell_barycenter(e);
281 }
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);
284 }
285 state.assembler->update_lame_params(lambdas, mus);
286
287 Eigen::VectorXd params(lambdas.size() + mus.size());
288 params << lambdas, mus;
289
290 return parametrization_.inverse_eval(params);
291 }
292
293 void FrictionCoeffientVariableToSimulation::update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices)
294 {
295 assert(state_variable.size() == 1);
296 assert(state_variable(0) >= 0);
297 for (auto state : states_)
298 state->args["contact"]["friction_coefficient"] = state_variable(0);
299 }
300 Eigen::VectorXd FrictionCoeffientVariableToSimulation::compute_adjoint_term(const Eigen::VectorXd &x) const
301 {
302 Eigen::VectorXd term, cur_term;
303 for (auto state : states_)
304 {
305 if (state->problem->is_time_dependent())
306 AdjointTools::dJ_friction_transient_adjoint_term(*state, state->get_adjoint_mat(1), state->get_adjoint_mat(0), cur_term);
307 else
308 log_and_throw_adjoint_error("[{}] Gradient in static simulations not implemented!", name());
309
310 if (term.size() != cur_term.size())
311 term = cur_term;
312 else
313 term += cur_term;
314 }
315 return apply_parametrization_jacobian(term, x);
316 }
318 {
319 log_and_throw_adjoint_error("[{}] inverse_eval not implemented!", name());
320 return Eigen::VectorXd();
321 }
322
323 void DampingCoeffientVariableToSimulation::update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices)
324 {
325 assert(state_variable.size() == 2);
326 json damping_param = {
327 {"psi", state_variable(0)},
328 {"phi", state_variable(1)},
329 };
330 for (auto state : states_)
331 {
332 if (!state->args["materials"].is_array())
333 {
334 state->args["materials"]["psi"] = damping_param["psi"];
335 state->args["materials"]["phi"] = damping_param["phi"];
336 }
337 else
338 {
339 for (auto &arg : state->args["materials"])
340 {
341 arg["psi"] = damping_param["psi"];
342 arg["phi"] = damping_param["phi"];
343 }
344 }
345
346 if (state->damping_assembler)
347 state->damping_assembler->add_multimaterial(0, damping_param, state->units);
348 }
349 logger().info("[{}] Current params: {}, {}", name(), state_variable(0), state_variable(1));
350 }
351 Eigen::VectorXd DampingCoeffientVariableToSimulation::compute_adjoint_term(const Eigen::VectorXd &x) const
352 {
353 Eigen::VectorXd term, cur_term;
354 for (auto state : states_)
355 {
356 if (state->problem->is_time_dependent())
357 AdjointTools::dJ_damping_transient_adjoint_term(*state, state->get_adjoint_mat(1), state->get_adjoint_mat(0), cur_term);
358 else
359 log_and_throw_adjoint_error("[{}] Static simulation not supported!", name());
360
361 if (term.size() != cur_term.size())
362 term = cur_term;
363 else
364 term += cur_term;
365 }
366 return apply_parametrization_jacobian(term, x);
367 }
369 {
370 log_and_throw_adjoint_error("[{}] inverse_eval not implemented!", name());
371 return Eigen::VectorXd();
372 }
373
374 void InitialConditionVariableToSimulation::update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices)
375 {
376 for (auto state : states_)
377 {
378 if (state_variable.size() != state->ndof() * 2)
379 {
380 log_and_throw_adjoint_error("[{}] Inconsistent number of parameters {} and number of dofs in forward {}!", name(), state_variable.size(), state->ndof() * 2);
381 }
382 state->initial_sol_update = state_variable.head(state->ndof());
383 state->initial_vel_update = state_variable.tail(state->ndof());
384 }
385 }
386 Eigen::VectorXd InitialConditionVariableToSimulation::compute_adjoint_term(const Eigen::VectorXd &x) const
387 {
388 Eigen::VectorXd term, cur_term;
389 for (auto state : states_)
390 {
391 if (state->problem->is_time_dependent())
392 AdjointTools::dJ_initial_condition_adjoint_term(*state, state->get_adjoint_mat(1), state->get_adjoint_mat(0), cur_term);
393 else
394 log_and_throw_adjoint_error("[{}] Static simulation not supported!", name());
395
396 if (term.size() != cur_term.size())
397 term = cur_term;
398 else
399 term += cur_term;
400 }
401 return apply_parametrization_jacobian(term, x);
402 }
404 {
405 auto &state = *states_[0];
406 Eigen::MatrixXd sol, vel;
407 state.initial_solution(sol);
408 state.initial_velocity(vel);
409
410 Eigen::VectorXd x(sol.size() + vel.size());
411 x << sol, vel;
412 return x;
413 }
414
415 void DirichletVariableToSimulation::update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices)
416 {
417 auto tensor_problem = std::dynamic_pointer_cast<polyfem::assembler::GenericTensorProblem>(states_[0]->problem);
418 assert(dirichlet_boundaries_.size() > 0);
419 int dim = states_[0]->mesh->dimension();
420 int num_steps = indices.size() / dim;
421 for (int i = 0; i < num_steps; ++i)
422 for (const int &b : dirichlet_boundaries_)
423 tensor_problem->update_dirichlet_boundary(b, indices(i * dim) + 1, state_variable.segment(i * dim, dim));
424
425 logger().info("Current dirichlet boundary {} is {}.", dirichlet_boundaries_[0], state_variable.transpose());
426 }
427
428 Eigen::VectorXd DirichletVariableToSimulation::compute_adjoint_term(const Eigen::VectorXd &x) const
429 {
430 Eigen::VectorXd term, cur_term;
431 for (auto state : states_)
432 {
433 if (state->problem->is_time_dependent())
434 AdjointTools::dJ_dirichlet_transient_adjoint_term(*state, state->get_adjoint_mat(1), state->get_adjoint_mat(0), cur_term);
435 else
436 log_and_throw_adjoint_error("[{}] Static dirichlet boundary optimization not supported!", name());
437
438 if (term.size() != cur_term.size())
439 term = cur_term;
440 else
441 term += cur_term;
442 }
443 return apply_parametrization_jacobian(term, x);
444 }
445 std::string DirichletVariableToSimulation::variable_to_string(const Eigen::VectorXd &variable)
446 {
447 return "";
448 }
450 {
451 assert(dirichlet_boundaries_.size() > 0);
452 assert(states_.size() > 0);
453
454 int dim = states_[0]->mesh->dimension();
455 Eigen::VectorXd x;
456 for (const auto &b : states_[0]->args["boundary_conditions"]["dirichlet_boundary"])
457 if (b["id"].get<int>() == dirichlet_boundaries_[0])
458 {
459 auto value = b["value"];
460 if (value[0].is_array())
461 {
462 if (!states_[0]->problem->is_time_dependent())
463 log_and_throw_adjoint_error("Simulation must be time dependent for timestep wise dirichlet.");
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);
468 }
469 else if (value[0].is_number())
470 {
471 if (states_[0]->problem->is_time_dependent())
472 log_and_throw_adjoint_error("Simulation must be quasistatic for single value dirichlet.");
473 x.resize(dim);
474 x = value;
475 }
476 else if (value.is_string())
477 assert(false);
478 break;
479 }
480
482 }
484 {
485 const std::string composite_map_type = args["composite_map_type"];
486 const State &state = *(states_[0]);
487 if (composite_map_type == "time_step_indexing")
488 {
489 const int time_steps = state.args["time"]["time_steps"];
490 const int dim = state.mesh->dimension();
491
492 output_indexing_.setZero(time_steps * dim);
493 for (int i = 0; i < time_steps; ++i)
494 for (int k = 0; k < dim; ++k)
495 output_indexing_(i * dim + k) = i;
496 }
497 else
499
500 set_dirichlet_boundaries(args["surface_selection"]);
501 }
502
503 void DirichletNodesVariableToSimulation::update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices)
504 {
505 for (auto state : states_)
506 {
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());
510 tensor_problem->update_dirichlet_nodes(state->in_node_to_node, dirichlet_nodes_, utils::unflatten(state_variable, state->mesh->dimension()));
511
512 logger().info("Updated dirichlet nodes");
513 }
514 }
515
516 Eigen::VectorXd DirichletNodesVariableToSimulation::compute_adjoint_term(const Eigen::VectorXd &x) const
517 {
518 Eigen::VectorXd term, cur_term;
519 for (auto state : states_)
520 {
521 if (state->problem->is_time_dependent())
522 log_and_throw_adjoint_error("[{}] Transient dirichlet node optimization not supported!", name());
523 else
524 AdjointTools::dJ_dirichlet_static_adjoint_term(*state, state->get_adjoint_mat(0), cur_term);
525
526 if (term.size() != cur_term.size())
527 term = cur_term;
528 else
529 term += cur_term;
530 }
531 return apply_parametrization_jacobian(term, x);
532 }
533 std::string DirichletNodesVariableToSimulation::variable_to_string(const Eigen::VectorXd &variable)
534 {
535 return "";
536 }
538 {
539 log_and_throw_adjoint_error("Inverse eval not implemented!");
540 }
542 {
543 json args_ = args;
544 const std::string composite_map_type = args_["composite_map_type"];
545 if (composite_map_type != "indices")
546 {
547 log_and_throw_adjoint_error("Must specify indices on which the nodal dirichlet is applied!");
548 }
549 else
550 {
551 set_dirichlet_nodes(args_["composite_map_indices"]);
552 int dim = 3;
553 std::vector<int> composite_map_indices = {};
554 for (int i = 0; i < dirichlet_nodes_.size(); ++i)
555 for (int k = 0; k < dim; ++k)
556 composite_map_indices.push_back(dirichlet_nodes_[i] * dim + k);
557 args_["composite_map_indices"] = composite_map_indices;
558 }
560 }
561
562 void PressureVariableToSimulation::update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices)
563 {
564 auto tensor_problem = std::dynamic_pointer_cast<polyfem::assembler::GenericTensorProblem>(states_[0]->problem);
565 assert(pressure_boundaries_.size() > 0);
566 for (int i = 0; i < indices.size(); ++i)
567 for (const int &b : pressure_boundaries_)
568 tensor_problem->update_pressure_boundary(b, indices(i) + 1, state_variable(i));
569
570 logger().info("Current pressure boundary {} is {}.", pressure_boundaries_[0], state_variable.transpose());
571 }
572
573 Eigen::VectorXd PressureVariableToSimulation::compute_adjoint_term(const Eigen::VectorXd &x) const
574 {
575 Eigen::VectorXd term, cur_term;
576 for (auto state : states_)
577 {
578 if (state->problem->is_time_dependent())
579 {
580 Eigen::MatrixXd adjoint_nu, adjoint_p;
581 adjoint_nu = state->get_adjoint_mat(1);
582 adjoint_p = state->get_adjoint_mat(0);
583 AdjointTools::dJ_pressure_transient_adjoint_term(*state, pressure_boundaries_, adjoint_nu, adjoint_p, cur_term);
584 }
585 else
586 {
587 AdjointTools::dJ_pressure_static_adjoint_term(*state, pressure_boundaries_, state->diff_cached.u(0), state->get_adjoint_mat(0), cur_term);
588 }
589 if (term.size() != cur_term.size())
590 term = cur_term;
591 else
592 term += cur_term;
593 }
594 return apply_parametrization_jacobian(term, x);
595 }
596
597 std::string PressureVariableToSimulation::variable_to_string(const Eigen::VectorXd &variable)
598 {
599 return "";
600 }
601
603 {
604 assert(pressure_boundaries_.size() > 0);
605 assert(states_.size() > 0);
606
607 Eigen::VectorXd x;
608 for (const auto &b : states_[0]->args["boundary_conditions"]["pressure_boundary"])
609 if (b["id"].get<int>() == pressure_boundaries_[0])
610 {
611 auto value = b["value"];
612 if (value.is_array())
613 {
614 if (!states_[0]->problem->is_time_dependent())
615 log_and_throw_adjoint_error("Simulation must be time dependent for timestep wise pressure.");
616 Eigen::VectorXd pressures = value;
617 x = pressures.segment(1, pressures.size() - 1);
618 }
619 else if (value.is_number())
620 {
621 if (states_[0]->problem->is_time_dependent())
622 log_and_throw_adjoint_error("Simulation must be quasistatic for single value pressure.");
623 x.resize(1);
624 x(0) = value;
625 }
626 else if (value.is_string())
627 assert(false);
628 break;
629 }
630
632 }
633
635 {
636 const std::string composite_map_type = args["composite_map_type"];
637 const State &state = *(states_[0]);
638 if (composite_map_type == "time_step_indexing")
639 {
640 const int time_steps = state.args["time"]["time_steps"];
641 output_indexing_.setZero(time_steps);
642 for (int i = 0; i < time_steps; ++i)
643 output_indexing_(i) = i;
644 }
645 else
647
648 set_pressure_boundaries(args["surface_selection"]);
649 }
650
651 Eigen::VectorXd PeriodicShapeVariableToSimulation::compute_adjoint_term(const Eigen::VectorXd &x) const
652 {
653 Eigen::VectorXd term, cur_term;
654 for (auto state : states_)
655 {
656 if (state->problem->is_time_dependent())
657 {
658 log_and_throw_error("Not implemented!");
659 }
660 else
661 {
662 AdjointTools::dJ_periodic_shape_adjoint_term(*state, *periodic_mesh_map, periodic_mesh_representation, state->diff_cached.u(0), state->get_adjoint_mat(0), cur_term);
663 }
664 if (term.size() != cur_term.size())
665 term = cur_term;
666 else
667 term += cur_term;
668 }
670 }
671 void PeriodicShapeVariableToSimulation::update(const Eigen::VectorXd &x)
672 {
673 const int dim = states_[0]->mesh->dimension();
675 const Eigen::MatrixXd V = utils::unflatten(periodic_mesh_map->eval(periodic_mesh_representation), dim);
676
677 for (auto state : states_)
678 {
679 const int n_verts = state->mesh->n_vertices();
680
681 for (int i = 0; i < n_verts; i++)
682 state->set_mesh_vertex(i, V.row(i));
683 }
684 }
686 {
687 const auto &state = *(states_[0]);
688
689 Eigen::MatrixXd V;
690 state.get_vertices(V);
691
692 if (!state.periodic_bc->all_direction_periodic())
693 log_and_throw_error("Cannot inverse evaluate periodic shape!");
694
695 periodic_mesh_map = std::make_unique<PeriodicMeshToMesh>(V);
697
699 }
700 Eigen::VectorXd PeriodicShapeVariableToSimulation::apply_parametrization_jacobian(const Eigen::VectorXd &term, const Eigen::VectorXd &x) const
701 {
702 const Eigen::VectorXd mid = periodic_mesh_map->apply_jacobian(term, periodic_mesh_representation);
703 return parametrization_.apply_jacobian(mid, x);
704 }
705} // namespace polyfem::solver
int V
StiffnessMatrix tmp_mat
int x
main class that contains the polyfem solver and all its state
Definition State.hpp:79
std::unique_ptr< mesh::Mesh > mesh
current mesh, it can be a Mesh2D or Mesh3D
Definition State.hpp:471
json args
main input arguments containing all defaults
Definition State.hpp:101
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
Definition Mesh2D.cpp:69
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
void update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices) override
void update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices) override
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
Eigen::VectorXd compute_adjoint_term(const Eigen::VectorXd &x) const 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_dirichlet_boundaries(const std::vector< int > &dirichlet_boundaries)
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 compute_adjoint_term(const Eigen::VectorXd &x) const override
void update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices) 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 compute_adjoint_term(const Eigen::VectorXd &x) const override
Eigen::VectorXd apply_parametrization_jacobian(const Eigen::VectorXd &term, const Eigen::VectorXd &x) const override
std::unique_ptr< PeriodicMeshToMesh > periodic_mesh_map
void set_pressure_boundaries(const std::vector< int > &pressure_boundaries)
std::string variable_to_string(const Eigen::VectorXd &variable)
void update_state(const Eigen::VectorXd &state_variable, const Eigen::VectorXi &indices) 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
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)
virtual Eigen::VectorXd apply_parametrization_jacobian(const Eigen::VectorXd &term, const Eigen::VectorXd &x) const
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 &&parametrization)
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.
Definition MatrixIO.cpp:18
void dJ_pressure_static_adjoint_term(const State &state, const std::vector< int > &boundary_ids, const Eigen::MatrixXd &sol, const Eigen::MatrixXd &adjoint, Eigen::VectorXd &one_form)
void dJ_periodic_shape_adjoint_term(const State &state, const PeriodicMeshToMesh &periodic_mesh_map, const Eigen::VectorXd &periodic_mesh_representation, const Eigen::MatrixXd &sol, const Eigen::MatrixXd &adjoint, Eigen::VectorXd &one_form)
void dJ_dirichlet_static_adjoint_term(const State &state, const Eigen::MatrixXd &adjoint, Eigen::VectorXd &one_form)
void dJ_friction_transient_adjoint_term(const State &state, const Eigen::MatrixXd &adjoint_nu, const Eigen::MatrixXd &adjoint_p, Eigen::VectorXd &one_form)
void dJ_pressure_transient_adjoint_term(const State &state, const std::vector< int > &boundary_ids, const Eigen::MatrixXd &adjoint_nu, const Eigen::MatrixXd &adjoint_p, Eigen::VectorXd &one_form)
void dJ_initial_condition_adjoint_term(const State &state, const Eigen::MatrixXd &adjoint_nu, const Eigen::MatrixXd &adjoint_p, Eigen::VectorXd &one_form)
void dJ_material_transient_adjoint_term(const State &state, const Eigen::MatrixXd &adjoint_nu, const Eigen::MatrixXd &adjoint_p, Eigen::VectorXd &one_form)
void dJ_shape_homogenization_adjoint_term(const State &state, const Eigen::MatrixXd &sol, const Eigen::MatrixXd &adjoint, Eigen::VectorXd &one_form)
void dJ_shape_transient_adjoint_term(const State &state, const Eigen::MatrixXd &adjoint_nu, const Eigen::MatrixXd &adjoint_p, Eigen::VectorXd &one_form)
void dJ_material_static_adjoint_term(const State &state, const Eigen::MatrixXd &sol, const Eigen::MatrixXd &adjoint, Eigen::VectorXd &one_form)
void dJ_dirichlet_transient_adjoint_term(const State &state, const Eigen::MatrixXd &adjoint_nu, const Eigen::MatrixXd &adjoint_p, Eigen::VectorXd &one_form)
void dJ_shape_static_adjoint_term(const State &state, const Eigen::MatrixXd &sol, const Eigen::MatrixXd &adjoint, Eigen::VectorXd &one_form)
void dJ_damping_transient_adjoint_term(const State &state, const Eigen::MatrixXd &adjoint_nu, const Eigen::MatrixXd &adjoint_p, Eigen::VectorXd &one_form)
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.
Definition Logger.cpp:44
nlohmann::json json
Definition Common.hpp:9
void log_and_throw_adjoint_error(const std::string &msg)
Definition Logger.cpp:79
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > RowVectorNd
Definition Types.hpp:13
void log_and_throw_error(const std::string &msg)
Definition Logger.cpp:73
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)