Loading [MathJax]/extensions/tex2jax.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Optimizations.cpp
Go to the documentation of this file.
1#include "Optimizations.hpp"
2
4#include <jse/jse.h>
5
7
8#include <polyfem/State.hpp>
18
21
23
27
28#include <polysolve/nonlinear/BoxConstraintSolver.hpp>
29
30namespace spdlog::level
31{
33 spdlog::level::level_enum,
34 {{spdlog::level::level_enum::trace, "trace"},
35 {spdlog::level::level_enum::debug, "debug"},
36 {spdlog::level::level_enum::info, "info"},
37 {spdlog::level::level_enum::warn, "warning"},
38 {spdlog::level::level_enum::err, "error"},
39 {spdlog::level::level_enum::critical, "critical"},
40 {spdlog::level::level_enum::off, "off"},
41 {spdlog::level::level_enum::trace, 0},
42 {spdlog::level::level_enum::debug, 1},
43 {spdlog::level::level_enum::info, 2},
44 {spdlog::level::level_enum::warn, 3},
45 {spdlog::level::level_enum::err, 3},
46 {spdlog::level::level_enum::critical, 4},
47 {spdlog::level::level_enum::off, 5}})
48}
49
50namespace polyfem::solver
51{
52 namespace
53 {
54 bool load_json(const std::string &json_file, json &out)
55 {
56 std::ifstream file(json_file);
57
58 if (!file.is_open())
59 return false;
60
61 file >> out;
62
63 out["root_path"] = json_file;
64
65 return true;
66 }
67 } // namespace
68
69 std::shared_ptr<polysolve::nonlinear::Solver> AdjointOptUtils::make_nl_solver(const json &solver_params, const json &linear_solver_params, const double characteristic_length)
70 {
71 auto names = polysolve::nonlinear::Solver::available_solvers();
72 if (std::find(names.begin(), names.end(), solver_params["solver"]) != names.end())
73 return polysolve::nonlinear::Solver::create(solver_params, linear_solver_params, characteristic_length, adjoint_logger());
74
75 names = polysolve::nonlinear::BoxConstraintSolver::available_solvers();
76 if (std::find(names.begin(), names.end(), solver_params["solver"]) != names.end())
77 return polysolve::nonlinear::BoxConstraintSolver::create(solver_params, linear_solver_params, characteristic_length, adjoint_logger());
78
79 log_and_throw_adjoint_error("Invalid nonlinear solver name!");
80 }
81
82 std::shared_ptr<AdjointForm> AdjointOptUtils::create_simple_form(const std::string &obj_type, const std::string &param_type, const std::shared_ptr<State> &state, const json &args)
83 {
84 std::shared_ptr<AdjointForm> obj;
87 if (obj_type == "compliance")
88 obj = std::make_shared<ComplianceForm>(var2sim, *state, args);
89 else if (obj_type == "acceleration")
90 obj = std::make_shared<AccelerationForm>(var2sim, *state, args);
91 else if (obj_type == "kinetic")
92 obj = std::make_shared<AccelerationForm>(var2sim, *state, args);
93 else if (obj_type == "position")
94 obj = std::make_shared<PositionForm>(var2sim, *state, args);
95 else if (obj_type == "stress")
96 obj = std::make_shared<StressForm>(var2sim, *state, args);
97 else if (obj_type == "stress_norm")
98 obj = std::make_shared<StressNormForm>(var2sim, *state, args);
99 else if (obj_type == "dirichlet_energy")
100 obj = std::make_shared<DirichletEnergyForm>(var2sim, *state, args);
101 else if (obj_type == "elastic_energy")
102 obj = std::make_shared<ElasticEnergyForm>(var2sim, *state, args);
103 else if (obj_type == "max_stress")
104 obj = std::make_shared<MaxStressForm>(var2sim, *state, args);
105 else if (obj_type == "volume")
106 obj = std::make_shared<VolumeForm>(var2sim, *state, args);
107 else if (obj_type == "min_jacobian")
108 obj = std::make_shared<MinJacobianForm>(var2sim, *state);
109 else if (obj_type == "AMIPS")
110 obj = std::make_shared<AMIPSForm>(var2sim, *state);
111 else if (obj_type == "layer_thickness")
112 obj = std::make_shared<LayerThicknessForm>(var2sim, *state, args["boundary_ids"].get<std::vector<int>>(), args["dhat"]);
113 else if (obj_type == "layer_thickness_log")
114 obj = std::make_shared<LayerThicknessForm>(var2sim, *state, args["boundary_ids"].get<std::vector<int>>(), args["dhat"], true, args["dmin"]);
115 else if (obj_type == "function-target")
116 {
117 std::shared_ptr<TargetForm> tmp = std::make_shared<TargetForm>(var2sim, *state, args);
118 tmp->set_reference(args["target_function"], args["target_function_gradient"]);
119 obj = tmp;
120 }
121 else if (obj_type == "boundary_smoothing")
122 {
123 if (args["surface_selection"].is_array())
124 obj = std::make_shared<BoundarySmoothingForm>(var2sim, *state, args["scale_invariant"], args["power"], args["surface_selection"].get<std::vector<int>>(), args["dimensions"].get<std::vector<int>>());
125 else
126 obj = std::make_shared<BoundarySmoothingForm>(var2sim, *state, args["scale_invariant"], args["power"], std::vector<int>{args["surface_selection"].get<int>()}, args["dimensions"].get<std::vector<int>>());
127 }
128 else
129 log_and_throw_adjoint_error("Invalid simple form type {}!", obj_type);
130
131 return obj;
132 }
133
134 std::shared_ptr<AdjointForm> AdjointOptUtils::create_form(const json &args, const VariableToSimulationGroup &var2sim, const std::vector<std::shared_ptr<State>> &states)
135 {
136 std::shared_ptr<AdjointForm> obj;
137 if (args.is_array())
138 {
139 std::vector<std::shared_ptr<AdjointForm>> forms;
140 for (const auto &arg : args)
141 forms.push_back(create_form(arg, var2sim, states));
142
143 obj = std::make_shared<SumCompositeForm>(var2sim, forms);
144 }
145 else
146 {
147 const std::string type = args["type"];
148 if (type == "transient_integral")
149 {
150 std::shared_ptr<StaticForm> static_obj = std::dynamic_pointer_cast<StaticForm>(create_form(args["static_objective"], var2sim, states));
151 if (!static_obj)
152 log_and_throw_adjoint_error("Transient integral objective must have a static objective!");
153 const auto &state = states[args["state"]];
154 obj = std::make_shared<TransientForm>(var2sim, state->args["time"]["time_steps"], state->args["time"]["dt"], args["integral_type"], args["steps"].get<std::vector<int>>(), static_obj);
155 }
156 else if (type == "proxy_transient_integral")
157 {
158 std::shared_ptr<StaticForm> static_obj = std::dynamic_pointer_cast<StaticForm>(create_form(args["static_objective"], var2sim, states));
159 if (!static_obj)
160 log_and_throw_adjoint_error("Transient integral objective must have a static objective!");
161 if (args["steps"].size() == 0)
162 log_and_throw_adjoint_error("ProxyTransientForm requires non-empty \"steps\"!");
163 const auto &state = states[args["state"]];
164 obj = std::make_shared<ProxyTransientForm>(var2sim, state->args["time"]["time_steps"], state->args["time"]["dt"], args["integral_type"], args["steps"].get<std::vector<int>>(), static_obj);
165 }
166 else if (type == "power")
167 {
168 std::shared_ptr<AdjointForm> obj_aux = create_form(args["objective"], var2sim, states);
169 obj = std::make_shared<PowerForm>(obj_aux, args["power"]);
170 }
171 else if (type == "divide")
172 {
173 std::shared_ptr<AdjointForm> obj1 = create_form(args["objective"][0], var2sim, states);
174 std::shared_ptr<AdjointForm> obj2 = create_form(args["objective"][1], var2sim, states);
175 std::vector<std::shared_ptr<AdjointForm>> objs({obj1, obj2});
176 obj = std::make_shared<DivideForm>(objs);
177 }
178 else if (type == "plus-const")
179 {
180 obj = std::make_shared<PlusConstCompositeForm>(create_form(args["objective"], var2sim, states), args["value"]);
181 }
182 else if (type == "log")
183 {
184 obj = std::make_shared<LogCompositeForm>(create_form(args["objective"], var2sim, states));
185 }
186 else if (type == "compliance")
187 {
188 obj = std::make_shared<ComplianceForm>(var2sim, *(states[args["state"]]), args);
189 }
190 else if (type == "acceleration")
191 {
192 obj = std::make_shared<AccelerationForm>(var2sim, *(states[args["state"]]), args);
193 }
194 else if (type == "kinetic")
195 {
196 obj = std::make_shared<AccelerationForm>(var2sim, *(states[args["state"]]), args);
197 }
198 else if (type == "target")
199 {
200 std::shared_ptr<TargetForm> tmp = std::make_shared<TargetForm>(var2sim, *(states[args["state"]]), args);
201 auto reference_cached = args["reference_cached_body_ids"].get<std::vector<int>>();
202 tmp->set_reference(states[args["target_state"]], std::set(reference_cached.begin(), reference_cached.end()));
203 obj = tmp;
204 }
205 else if (type == "displacement-target")
206 {
207 std::shared_ptr<TargetForm> tmp = std::make_shared<TargetForm>(var2sim, *(states[args["state"]]), args);
208 Eigen::VectorXd target_displacement;
209 target_displacement.setZero(states[args["state"]]->mesh->dimension());
210 if (target_displacement.size() != args["target_displacement"].size())
211 log_and_throw_error("Target displacement shape must match the dimension of the simulation");
212 for (int i = 0; i < target_displacement.size(); ++i)
213 target_displacement(i) = args["target_displacement"][i].get<double>();
214 if (args["active_dimension"].size() > 0)
215 {
216 if (target_displacement.size() != args["active_dimension"].size())
217 log_and_throw_error("Active dimension shape must match the dimension of the simulation");
218 std::vector<bool> active_dimension_mask(args["active_dimension"].size());
219 for (int i = 0; i < args["active_dimension"].size(); ++i)
220 active_dimension_mask[i] = args["active_dimension"][i].get<bool>();
221 tmp->set_active_dimension(active_dimension_mask);
222 }
223 tmp->set_reference(target_displacement);
224 obj = tmp;
225 }
226 else if (type == "center-target")
227 {
228 obj = std::make_shared<BarycenterTargetForm>(var2sim, args, states[args["state"]], states[args["target_state"]]);
229 }
230 else if (type == "sdf-target")
231 {
232 std::shared_ptr<SDFTargetForm> tmp = std::make_shared<SDFTargetForm>(var2sim, *(states[args["state"]]), args);
233 double delta = args["delta"].get<double>();
234 if (!states[args["state"]]->mesh->is_volume())
235 {
236 int dim = 2;
237 Eigen::MatrixXd control_points(args["control_points"].size(), dim);
238 for (int i = 0; i < control_points.rows(); ++i)
239 for (int j = 0; j < control_points.cols(); ++j)
240 control_points(i, j) = args["control_points"][i][j].get<double>();
241 Eigen::VectorXd knots(args["knots"].size());
242 for (int i = 0; i < knots.size(); ++i)
243 knots(i) = args["knots"][i].get<double>();
244 tmp->set_bspline_target(control_points, knots, delta);
245 }
246 else
247 {
248 int dim = 3;
249 Eigen::MatrixXd control_points_grid(args["control_points_grid"].size(), dim);
250 for (int i = 0; i < control_points_grid.rows(); ++i)
251 for (int j = 0; j < control_points_grid.cols(); ++j)
252 control_points_grid(i, j) = args["control_points_grid"][i][j].get<double>();
253 Eigen::VectorXd knots_u(args["knots_u"].size());
254 for (int i = 0; i < knots_u.size(); ++i)
255 knots_u(i) = args["knots_u"][i].get<double>();
256 Eigen::VectorXd knots_v(args["knots_v"].size());
257 for (int i = 0; i < knots_v.size(); ++i)
258 knots_v(i) = args["knots_v"][i].get<double>();
259 tmp->set_bspline_target(control_points_grid, knots_u, knots_v, delta);
260 }
261
262 obj = tmp;
263 }
264 else if (type == "mesh-target")
265 {
266 std::shared_ptr<MeshTargetForm> tmp = std::make_shared<MeshTargetForm>(var2sim, *(states[args["state"]]), args);
267 double delta = args["delta"].get<double>();
268 Eigen::MatrixXd V;
269 Eigen::MatrixXi E, F;
270 bool read = polyfem::io::OBJReader::read(args["mesh_path"], V, E, F);
271 if (!read)
272 log_and_throw_error(fmt::format("Could not read mesh! {}", args["mesh"]));
273 tmp->set_surface_mesh_target(V, F, delta);
274 obj = tmp;
275 }
276 else if (type == "function-target")
277 {
278 std::shared_ptr<TargetForm> tmp = std::make_shared<TargetForm>(var2sim, *(states[args["state"]]), args);
279 tmp->set_reference(args["target_function"], args["target_function_gradient"]);
280 obj = tmp;
281 }
282 else if (type == "node-target")
283 {
284 obj = std::make_shared<NodeTargetForm>(*(states[args["state"]]), var2sim, args);
285 }
286 else if (type == "min-dist-target")
287 {
288 obj = std::make_shared<MinTargetDistForm>(var2sim, args["steps"], args["target"], args, states[args["state"]]);
289 }
290 else if (type == "position")
291 {
292 obj = std::make_shared<PositionForm>(var2sim, *(states[args["state"]]), args);
293 }
294 else if (type == "stress")
295 {
296 obj = std::make_shared<StressForm>(var2sim, *(states[args["state"]]), args);
297 }
298 else if (type == "stress_norm")
299 {
300 obj = std::make_shared<StressNormForm>(var2sim, *(states[args["state"]]), args);
301 }
302 else if (type == "dirichlet_energy")
303 {
304 obj = std::make_shared<DirichletEnergyForm>(var2sim, *(states[args["state"]]), args);
305 }
306 else if (type == "elastic_energy")
307 {
308 obj = std::make_shared<ElasticEnergyForm>(var2sim, *(states[args["state"]]), args);
309 }
310 else if (type == "quadratic_contact_force_norm")
311 {
312 obj = std::make_shared<ProxyContactForceForm>(var2sim, *(states[args["state"]]), args["dhat"], true, args);
313 }
314 else if (type == "log_contact_force_norm")
315 {
316 obj = std::make_shared<ProxyContactForceForm>(var2sim, *(states[args["state"]]), args["dhat"], false, args);
317 }
318 else if (type == "max_stress")
319 {
320 obj = std::make_shared<MaxStressForm>(var2sim, *(states[args["state"]]), args);
321 }
322 else if (type == "smooth_contact_force_norm")
323 {
324 // assert(states[args["state"]]->args["contact"]["use_gcp_formulation"]);
325 obj = std::make_shared<SmoothContactForceForm>(var2sim, *(states[args["state"]]), args);
326 }
327 else if (type == "volume")
328 {
329 obj = std::make_shared<VolumeForm>(var2sim, *(states[args["state"]]), args);
330 }
331 else if (type == "soft_constraint")
332 {
333 std::vector<std::shared_ptr<AdjointForm>> forms({create_form(args["objective"], var2sim, states)});
334 Eigen::VectorXd bounds = args["soft_bound"];
335 obj = std::make_shared<InequalityConstraintForm>(forms, bounds, args["power"]);
336 }
337 else if (type == "min_jacobian")
338 {
339 obj = std::make_shared<MinJacobianForm>(var2sim, *(states[args["state"]]));
340 }
341 else if (type == "AMIPS")
342 {
343 obj = std::make_shared<AMIPSForm>(var2sim, *(states[args["state"]]));
344 }
345 else if (type == "boundary_smoothing")
346 {
347 if (args["surface_selection"].is_array())
348 obj = std::make_shared<BoundarySmoothingForm>(var2sim, *(states[args["state"]]), args["scale_invariant"], args["power"], args["surface_selection"].get<std::vector<int>>(), args["dimensions"].get<std::vector<int>>());
349 else
350 obj = std::make_shared<BoundarySmoothingForm>(var2sim, *(states[args["state"]]), args["scale_invariant"], args["power"], std::vector<int>{args["surface_selection"].get<int>()}, args["dimensions"].get<std::vector<int>>());
351 }
352 else if (type == "collision_barrier")
353 {
354 obj = std::make_shared<CollisionBarrierForm>(var2sim, *(states[args["state"]]), args["dhat"]);
355 }
356 else if (type == "layer_thickness")
357 {
358 obj = std::make_shared<LayerThicknessForm>(var2sim, *(states[args["state"]]), args["boundary_ids"].get<std::vector<int>>(), args["dhat"]);
359 }
360 else if (type == "layer_thickness_log")
361 {
362 obj = std::make_shared<LayerThicknessForm>(var2sim, *(states[args["state"]]), args["boundary_ids"].get<std::vector<int>>(), args["dhat"], true, args["dmin"]);
363 }
364 else if (type == "deformed_collision_barrier")
365 {
366 obj = std::make_shared<DeformedCollisionBarrierForm>(var2sim, *(states[args["state"]]), args["dhat"]);
367 }
368 else if (type == "parametrized_product")
369 {
370 std::vector<std::shared_ptr<Parametrization>> map_list;
371 for (const auto &arg : args["parametrization"])
372 map_list.push_back(create_parametrization(arg, states, {}));
373 obj = std::make_shared<ParametrizedProductForm>(CompositeParametrization(std::move(map_list)));
374 }
375 else
376 log_and_throw_adjoint_error("Objective not implemented!");
377
378 obj->set_weight(args["weight"]);
379 if (args["print_energy"].get<std::string>() != "")
380 obj->enable_energy_print(args["print_energy"]);
381 }
382
383 return obj;
384 }
385
386 std::shared_ptr<Parametrization> AdjointOptUtils::create_parametrization(const json &args, const std::vector<std::shared_ptr<State>> &states, const std::vector<int> &variable_sizes)
387 {
388 std::shared_ptr<Parametrization> map;
389 const std::string type = args["type"];
390 if (type == "per-body-to-per-elem")
391 {
392 map = std::make_shared<PerBody2PerElem>(*(states[args["state"]]->mesh));
393 }
394 else if (type == "per-body-to-per-node")
395 {
396 map = std::make_shared<PerBody2PerNode>(*(states[args["state"]]->mesh), states[args["state"]]->bases, states[args["state"]]->n_bases);
397 }
398 else if (type == "E-nu-to-lambda-mu")
399 {
400 map = std::make_shared<ENu2LambdaMu>(args["is_volume"]);
401 }
402 else if (type == "slice")
403 {
404 if (args["from"] != -1 || args["to"] != -1)
405 map = std::make_shared<SliceMap>(args["from"], args["to"], args["last"]);
406 else if (args["parameter_index"] != -1)
407 {
408 int idx = args["parameter_index"].get<int>();
409 int from, to, last;
410 int cumulative = 0;
411 for (int i = 0; i < variable_sizes.size(); ++i)
412 {
413 if (i == idx)
414 {
415 from = cumulative;
416 to = from + variable_sizes[i];
417 }
418 cumulative += variable_sizes[i];
419 }
420 last = cumulative;
421 map = std::make_shared<SliceMap>(from, to, last);
422 }
423 else
424 log_and_throw_adjoint_error("Incorrect spec for SliceMap!");
425 }
426 else if (type == "exp")
427 {
428 map = std::make_shared<ExponentialMap>(args["from"], args["to"]);
429 }
430 else if (type == "scale")
431 {
432 map = std::make_shared<Scaling>(args["value"]);
433 }
434 else if (type == "power")
435 {
436 map = std::make_shared<PowerMap>(args["power"]);
437 }
438 else if (type == "append-values")
439 {
440 Eigen::VectorXd vals = args["values"];
441 map = std::make_shared<InsertConstantMap>(vals, args["start"]);
442 }
443 else if (type == "append-const")
444 {
445 map = std::make_shared<InsertConstantMap>(args["size"], args["value"], args["start"]);
446 }
447 else if (type == "linear-filter")
448 {
449 map = std::make_shared<LinearFilter>(*(states[args["state"]]->mesh), args["radius"]);
450 }
451 else if (type == "bounded-biharmonic-weights")
452 {
453 map = std::make_shared<BoundedBiharmonicWeights2Dto3D>(args["num_control_vertices"], args["num_vertices"], *states[args["state"]], args["allow_rotations"]);
454 }
455 else if (type == "scalar-velocity-parametrization")
456 {
457 map = std::make_shared<ScalarVelocityParametrization>(args["start_val"], args["dt"]);
458 }
459 else
460 log_and_throw_adjoint_error("Unkown parametrization!");
461
462 return map;
463 }
464
465 std::unique_ptr<VariableToSimulation> AdjointOptUtils::create_variable_to_simulation(const json &args, const std::vector<std::shared_ptr<State>> &states, const std::vector<int> &variable_sizes)
466 {
467 const std::string type = args["type"];
468
469 std::vector<std::shared_ptr<State>> cur_states;
470 if (args["state"].is_array())
471 for (int i : args["state"])
472 cur_states.push_back(states[i]);
473 else
474 cur_states.push_back(states[args["state"]]);
475
476 std::vector<std::shared_ptr<Parametrization>> map_list;
477 for (const auto &arg : args["composition"])
478 map_list.push_back(create_parametrization(arg, states, variable_sizes));
479
480 std::unique_ptr<VariableToSimulation> var2sim = VariableToSimulation::create(type, cur_states, CompositeParametrization(std::move(map_list)));
481 var2sim->set_output_indexing(args);
482
483 return var2sim;
484 }
485
486 Eigen::VectorXd AdjointOptUtils::inverse_evaluation(const json &args, const int ndof, const std::vector<int> &variable_sizes, VariableToSimulationGroup &var2sim)
487 {
488 Eigen::VectorXd x;
489 x.setZero(ndof);
490 int accumulative = 0;
491 int var = 0;
492 for (const auto &arg : args)
493 {
494 const auto &arg_initial = arg["initial"];
495 Eigen::VectorXd tmp(variable_sizes[var]);
496 if (arg_initial.is_array() && arg_initial.size() > 0)
497 {
498 tmp = arg_initial;
499 x.segment(accumulative, tmp.size()) = tmp;
500 }
501 else if (arg_initial.is_number())
502 {
503 tmp.setConstant(arg_initial.get<double>());
504 x.segment(accumulative, tmp.size()) = tmp;
505 }
506 else // arg["initial"] is empty array
507 x += var2sim[var]->inverse_eval();
508
509 accumulative += tmp.size();
510 var++;
511 }
512
513 return x;
514 }
515
516 std::shared_ptr<State> AdjointOptUtils::create_state(const json &args, CacheLevel level, const size_t max_threads)
517 {
518 std::shared_ptr<State> state = std::make_shared<State>();
519 state->set_max_threads(max_threads);
520
521 json in_args = args;
522 in_args["solver"]["max_threads"] = max_threads;
523 if (!args.contains("output") || !args["output"].contains("log") || !args["output"]["log"].contains("level"))
524 {
525 const json tmp = R"({
526 "output": {
527 "log": {
528 "level": "error"
529 }
530 }
531 })"_json;
532
533 in_args.merge_patch(tmp);
534 }
535
536 state->optimization_enabled = level;
537 state->init(in_args, true);
538 state->load_mesh();
539 Eigen::MatrixXd sol, pressure;
540 state->build_basis();
541 state->assemble_rhs();
542 state->assemble_mass_mat();
543
544 return state;
545 }
546
547 std::vector<std::shared_ptr<State>> AdjointOptUtils::create_states(const json &state_args, const CacheLevel &level, const size_t max_threads)
548 {
549 std::vector<std::shared_ptr<State>> states(state_args.size());
550 int i = 0;
551 for (const json &args : state_args)
552 {
553 json cur_args;
554 if (!load_json(args["path"], cur_args))
555 log_and_throw_adjoint_error("Can't find json for State {}", i);
556
557 states[i++] = AdjointOptUtils::create_state(cur_args, level, max_threads);
558 }
559 return states;
560 }
561
563 {
564 state.assemble_rhs();
565 state.assemble_mass_mat();
566 Eigen::MatrixXd sol, pressure;
567 state.solve_problem(sol, pressure);
568 }
569
570 void apply_objective_json_spec(json &args, const json &rules)
571 {
572 if (args.is_array())
573 {
574 for (auto &arg : args)
575 apply_objective_json_spec(arg, rules);
576 }
577 else if (args.is_object())
578 {
579 jse::JSE jse;
580 const bool valid_input = jse.verify_json(args, rules);
581
582 if (!valid_input)
583 {
584 logger().error("invalid objective json:\n{}", jse.log2str());
585 throw std::runtime_error("Invald objective json file");
586 }
587
588 args = jse.inject_defaults(args, rules);
589
590 for (auto &it : args.items())
591 {
592 if (it.key().find("objective") != std::string::npos)
593 apply_objective_json_spec(it.value(), rules);
594 }
595 }
596 }
597
598 json AdjointOptUtils::apply_opt_json_spec(const json &input_args, bool strict_validation)
599 {
600 json args_in = input_args;
601
602 // CHECK validity json
603 json rules;
604 jse::JSE jse;
605 {
606 jse.strict = strict_validation;
607 std::ifstream file(POLYFEM_OPT_INPUT_SPEC);
608
609 if (file.is_open())
610 file >> rules;
611 else
612 {
613 logger().error("unable to open {} rules", POLYFEM_OPT_INPUT_SPEC);
614 throw std::runtime_error("Invald spec file");
615 }
616
617 jse.include_directories.push_back(POLYFEM_JSON_SPEC_DIR);
618 jse.include_directories.push_back(POLYSOLVE_JSON_SPEC_DIR);
619 rules = jse.inject_include(rules);
620
621 // polysolve::linear::Solver::apply_default_solver(rules, "/solver/linear");
622 }
623
624 // polysolve::linear::Solver::select_valid_solver(args_in["solver"]["linear"], logger());
625
626 const bool valid_input = jse.verify_json(args_in, rules);
627
628 if (!valid_input)
629 {
630 logger().error("invalid input json:\n{}", jse.log2str());
631 throw std::runtime_error("Invald input json file");
632 }
633
634 json args = jse.inject_defaults(args_in, rules);
635
636 json obj_rules;
637 {
638 const std::string polyfem_objective_spec = POLYFEM_OBJECTIVE_INPUT_SPEC;
639 std::ifstream file(polyfem_objective_spec);
640
641 if (file.is_open())
642 file >> obj_rules;
643 else
644 {
645 logger().error("unable to open {} rules", polyfem_objective_spec);
646 throw std::runtime_error("Invald spec file");
647 }
648 }
649 apply_objective_json_spec(args["functionals"], obj_rules);
650
651 if (args.contains("stopping_conditions"))
652 apply_objective_json_spec(args["stopping_conditions"], obj_rules);
653
654 return args;
655 }
656
657 int AdjointOptUtils::compute_variable_size(const json &args, const std::vector<std::shared_ptr<State>> &states)
658 {
659 if (args["number"].is_number())
660 {
661 return args["number"].get<int>();
662 }
663 else if (args["number"].is_null() && args["initial"].size() > 0)
664 {
665 return args["initial"].size();
666 }
667 else if (args["number"].is_object())
668 {
669 auto selection = args["number"];
670 if (selection.contains("surface_selection"))
671 {
672 auto surface_selection = selection["surface_selection"].get<std::vector<int>>();
673 auto state_id = selection["state"];
674 std::set<int> node_ids = {};
675 for (const auto &surface : surface_selection)
676 {
677 std::vector<int> ids;
678 states[state_id]->compute_surface_node_ids(surface, ids);
679 for (const auto &i : ids)
680 node_ids.insert(i);
681 }
682 return node_ids.size() * states[state_id]->mesh->dimension();
683 }
684 else if (selection.contains("volume_selection"))
685 {
686 auto volume_selection = selection["volume_selection"].get<std::vector<int>>();
687 auto state_id = selection["state"];
688 std::set<int> node_ids = {};
689 for (const auto &volume : volume_selection)
690 {
691 std::vector<int> ids;
692 states[state_id]->compute_volume_node_ids(volume, ids);
693 for (const auto &i : ids)
694 node_ids.insert(i);
695 }
696
697 if (selection["exclude_boundary_nodes"])
698 {
699 std::vector<int> ids;
700 states[state_id]->compute_total_surface_node_ids(ids);
701 for (const auto &i : ids)
702 node_ids.erase(i);
703 }
704
705 return node_ids.size() * states[state_id]->mesh->dimension();
706 }
707 }
708
709 log_and_throw_adjoint_error("Incorrect specification for parameters.");
710 return -1;
711 }
712} // namespace polyfem::solver
int V
ElementAssemblyValues vals
Definition Assembler.cpp:22
int x
main class that contains the polyfem solver and all its state
Definition State.hpp:79
void assemble_mass_mat()
assemble mass, step 4 of solve build mass matrix based on defined basis modifies mass (and maybe more...
Definition State.cpp:1496
void solve_problem(Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
solves the problems
Definition State.cpp:1692
void assemble_rhs()
compute rhs, step 3 of solve build rhs vector based on defined basis and given rhs of the problem mod...
Definition State.cpp:1620
static bool read(const std::string obj_file_name, std::vector< std::vector< double > > &V, std::vector< std::vector< double > > &TC, std::vector< std::vector< double > > &N, std::vector< std::vector< int > > &F, std::vector< std::vector< int > > &FTC, std::vector< std::vector< int > > &FN, std::vector< std::vector< int > > &L)
Read a mesh from an ascii obj file.
Definition OBJReader.cpp:32
A collection of VariableToSimulation.
static std::unique_ptr< VariableToSimulation > create(const std::string &type, const std::vector< std::shared_ptr< State > > &states, CompositeParametrization &&parametrization)
bool load_json(const std::string &json_file, json &out)
Definition main.cpp:26
NLOHMANN_JSON_SERIALIZE_ENUM(CollisionProxyTessellation, {{CollisionProxyTessellation::REGULAR, "regular"}, {CollisionProxyTessellation::IRREGULAR, "irregular"}})
void apply_objective_json_spec(json &args, const json &rules)
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:44
spdlog::logger & adjoint_logger()
Retrieves the current logger for adjoint.
Definition Logger.cpp:30
nlohmann::json json
Definition Common.hpp:9
void log_and_throw_adjoint_error(const std::string &msg)
Definition Logger.cpp:79
void log_and_throw_error(const std::string &msg)
Definition Logger.cpp:73
static std::shared_ptr< AdjointForm > create_form(const json &args, const VariableToSimulationGroup &var2sim, const std::vector< std::shared_ptr< State > > &states)
static std::shared_ptr< AdjointForm > create_simple_form(const std::string &obj_type, const std::string &param_type, const std::shared_ptr< State > &state, const json &args)
static std::shared_ptr< polysolve::nonlinear::Solver > make_nl_solver(const json &solver_params, const json &linear_solver_params, const double characteristic_length)
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)
static json apply_opt_json_spec(const json &input_args, bool strict_validation)
static Eigen::VectorXd inverse_evaluation(const json &args, const int ndof, const std::vector< int > &variable_sizes, VariableToSimulationGroup &var2sim)
static void solve_pde(State &state)
static std::shared_ptr< Parametrization > create_parametrization(const json &args, const std::vector< std::shared_ptr< State > > &states, const std::vector< int > &variable_sizes)
static int compute_variable_size(const json &args, const std::vector< std::shared_ptr< State > > &states)
static std::vector< std::shared_ptr< State > > create_states(const json &state_args, const CacheLevel &level, const size_t max_threads)
static std::shared_ptr< State > create_state(const json &args, CacheLevel level, const size_t max_threads)