21#include <polysolve/nonlinear/Solver.hpp>
23#include <spdlog/sinks/stdout_color_sinks.h>
24#include <spdlog/sinks/basic_file_sink.h>
25#include <spdlog/sinks/ostream_sink.h>
42 spdlog::level::level_enum,
43 {{spdlog::level::level_enum::trace,
"trace"},
44 {spdlog::level::level_enum::debug,
"debug"},
45 {spdlog::level::level_enum::info,
"info"},
46 {spdlog::level::level_enum::warn,
"warning"},
47 {spdlog::level::level_enum::err,
"error"},
48 {spdlog::level::level_enum::critical,
"critical"},
49 {spdlog::level::level_enum::off,
"off"},
50 {spdlog::level::level_enum::trace, 0},
51 {spdlog::level::level_enum::debug, 1},
52 {spdlog::level::level_enum::info, 2},
53 {spdlog::level::level_enum::warn, 3},
54 {spdlog::level::level_enum::err, 3},
55 {spdlog::level::level_enum::critical, 4},
56 {spdlog::level::level_enum::off, 5}})
63 bool is_failed_status(
const polysolve::nonlinear::Status status)
65 using Status = polysolve::nonlinear::Status;
66 return status == Status::NanEncountered
67 || status == Status::NotDescentDirection
68 || status == Status::LineSearchFailed
69 || status == Status::UpdateDirectionFailed
70 || status == Status::NotStarted
71 || status == Status::Continue;
74 std::string parse_remeshing_trigger(
const json &args)
76 const json &remeshing =
args[
"remeshing"];
77 if (remeshing.is_null())
80 const json &trigger = remeshing[
"trigger"];
81 if (trigger.is_null())
84 const bool periodic = trigger.contains(
"periodic") && trigger[
"periodic"].is_object();
85 const bool scaled_jacobian = trigger.contains(
"scaled_jacobian") && trigger[
"scaled_jacobian"].is_object();
86 if (periodic && scaled_jacobian)
91 return "scaled_jacobian";
95 std::vector<json> load_state_jsons(
96 const io::ResourceIO &resources,
98 std::vector<std::unique_ptr<const io::ResourceIO>> &state_resources)
100 std::vector<json> result;
101 state_resources.clear();
102 for (
int i = 0; i <
args.size(); ++i)
104 const std::string
path =
args[i][
"path"].get<std::string>();
105 if (!resources.exists(path))
107 result.push_back(json::parse(resources.read_string(path)));
108 const std::string parent = std::filesystem::path(path).parent_path().generic_string();
109 state_resources.push_back(resources.with_root(parent));
114 void validate_remeshing_json(
const json &args)
116 const std::string trigger = parse_remeshing_trigger(args);
117 if (trigger ==
"multiple")
118 log_and_throw_error(
"Optimization remeshing requires exactly one trigger: periodic or scaled_jacobian.");
119 const json &remeshing =
args[
"remeshing"];
120 const int state = remeshing[
"state"];
121 if (state >= args[
"states"].size())
124 for (
const auto &variable :
args[
"variable_to_simulation"])
126 if (variable[
"type"] !=
"shape")
129 const json &variable_state = variable[
"state"];
130 const bool targets_state = variable_state.is_array()
131 ? std::find(variable_state.begin(), variable_state.end(), state) != variable_state.end()
132 : variable_state.get<
int>() == state;
136 const json &selection = variable[
"active_geometry_nodes"];
137 const bool selects_all_nodes = selection.is_array() && selection.empty();
138 if (selects_all_nodes || selection.is_object())
142 "Optimization remeshing changes mesh vertex numbering. "
143 "Shape variable {} targeting state {} must use a geometry-based "
144 "active_geometry_nodes selector (interior, boundary, or "
145 "boundary_excluding_surface), or an empty selector for all nodes.",
146 variable.value(
"name",
"shape"), state);
150 void validate_remeshing_geometry(
const json &state_args,
const int geometry)
152 if (geometry < 0 || geometry >= state_args[
"geometry"].size())
153 log_and_throw_error(
"Optimization remeshing geometry index {} is out of range.", geometry);
154 if (state_args[
"geometry"][geometry].value(
"is_obstacle",
false))
159 double minimum_scaled_jacobian(
160 const Eigen::MatrixXd &vertices,
161 const Eigen::MatrixXi &elements)
163 Eigen::VectorXd quality;
165 return quality.minCoeff();
169 double minimum_scaled_jacobian(
const varform::DifferentiableVarForm &varform)
172 Eigen::MatrixXi elements;
173 varform.get_vertices(vertices);
174 varform.get_elements(elements);
175 return minimum_scaled_jacobian(vertices, elements);
184 bool remesh_and_write(
185 const OptState &opt_state,
188 const std::filesystem::path &path,
189 const std::optional<double> quality_threshold)
191 mesh::MmgOptions options;
192 options.optim =
true;
195 Eigen::MatrixXi elements;
196 opt_state.varforms[state]->get_vertices(vertices);
197 opt_state.varforms[state]->get_elements(elements);
199 Eigen::MatrixXd remeshed_vertices;
200 Eigen::MatrixXi remeshed_boundary;
201 Eigen::MatrixXi remeshed_elements;
202 bool has_shape_variable =
false;
203 for (
const auto &var2sim : opt_state.variable_to_simulations.data)
205 const auto shape_v2s = std::dynamic_pointer_cast<solver::ShapeVariableToSimulation>(var2sim);
206 if (!shape_v2s || !shape_v2s->affects_varform(*opt_state.varforms[state]))
208 has_shape_variable =
true;
210 for (
int other_state = 0; other_state < opt_state.varforms.size(); ++other_state)
212 if (other_state != state
213 && shape_v2s->affects_varform(*opt_state.varforms[other_state]))
216 "Optimization remeshing cannot remesh state {} while its shape variable also affects state {}. This is not supported currently.",
221 if (!has_shape_variable)
223 logger().error(
"Can not remesh state {} because it's not affected by shape variables. No reason to remesh if you are not doing shape optimization.", state);
227 const mesh::Mesh &mesh = opt_state.varforms[state]->get_mesh();
228 bool success =
false;
229 if (mesh.dimension() == 2)
231 success = mesh::remesh_2d(
232 vertices, elements, remeshed_vertices, remeshed_elements,
237 success = mesh::remesh_3d(
238 vertices, elements, remeshed_vertices, remeshed_boundary,
239 remeshed_elements, options);
246 const double input_quality = minimum_scaled_jacobian(vertices, elements);
247 const double output_quality = minimum_scaled_jacobian(remeshed_vertices, remeshed_elements);
249 "Optimization remeshing minimum scaled Jacobian: {} -> {}.",
250 input_quality, output_quality);
252 if (quality_threshold.has_value() && output_quality <= *quality_threshold)
255 "MMG remeshing failed to satisfy minimum scaled Jacobian threshold: "
256 "input={}, output={}, required > {}. This is because MMG use a different convergence criteria, please switch remeshing trigger to periodic.",
257 input_quality, output_quality, *quality_threshold);
261 path.string(), remeshed_vertices, remeshed_elements,
262 std::vector<int>(remeshed_elements.rows(), body_id),
263 opt_state.varforms[state]->get_mesh().is_volume(),
false);
279 std::filesystem::path root = std::filesystem::current_path();
281 root = input_args[
"root_path"].get<std::string>();
282 resources_ = std::make_unique<io::FileSystemIO>(root);
283 return run(std::move(input_args), *
resources_, strict_validation);
289 const bool strict_validation)
294 std::string mode = parse_remeshing_trigger(input_args);
295 bool remeshing_enabled = mode !=
"none";
296 if (remeshing_enabled)
297 validate_remeshing_json(input_args);
299 std::filesystem::path output_root = input_args[
"output"][
"directory"].get<std::string>();
300 if (!output_root.is_absolute())
301 output_root =
resources_->host_directory() / output_root;
302 int max_restarts = remeshing_enabled
303 ? input_args[
"remeshing"][
"max_restarts"].get<
int>()
310 if (remeshing_enabled)
312 state = input_args[
"remeshing"][
"state"].get<
int>();
313 geometry = input_args[
"remeshing"][
"geometry"].get<
int>();
314 validate_remeshing_geometry(
state_args[state], geometry);
315 body_id =
state_args[state][
"geometry"][geometry].value(
"volume_selection", 0);
318 int remesh_count = 0;
324 json round_args = input_args;
325 if (remeshing_enabled)
327 const std::filesystem::path round_dir = output_root / fmt::format(
"remesh_round_{:d}", remesh_count);
328 round_args[
"output"][
"directory"] = round_dir.string();
329 state_args[state][
"output"][
"directory"] = round_dir.string();
340 if (
args[
"compute_objective"].get<bool>())
346 const polysolve::nonlinear::Status status =
solve(
x);
347 if (is_failed_status(status))
349 logger().error(
"Optimization failed: {}.", polysolve::nonlinear::status_message(status));
356 if (remesh_count >= max_restarts)
358 logger().info(
"Reached the optimization remeshing limit of {} restart(s).", max_restarts);
362 const std::filesystem::path round_dir = output_root / fmt::format(
"remesh_round_{:d}", remesh_count);
363 const std::filesystem::path remeshed_path = std::filesystem::absolute(round_dir /
"remeshed.msh");
364 const std::optional<double> quality_threshold = mode ==
"scaled_jacobian"
365 ? std::optional<double>(
args[
"remeshing"][
"trigger"][
"scaled_jacobian"][
"quality_threshold"].get<double>())
367 if (!remesh_and_write(
368 *
this, state, body_id,
369 remeshed_path, quality_threshold))
371 logger().error(
"MMG failed to produce a valid optimization restart mesh.");
375 state_args[state][
"geometry"][geometry][
"mesh"] = remeshed_path.string();
376 logger().info(
"Restarting optimization from remeshed mesh {}.", remeshed_path.string());
382 const std::string &log_file,
383 const spdlog::level::level_enum log_level,
384 const spdlog::level::level_enum file_log_level,
387 std::vector<spdlog::sink_ptr> sinks;
391 console_sink_ = std::make_shared<spdlog::sinks::stdout_color_sink_mt>();
395 if (!log_file.empty())
397 file_sink_ = std::make_shared<spdlog::sinks::basic_file_sink_mt>(log_file,
true);
404 spdlog::flush_every(std::chrono::seconds(3));
409 std::vector<spdlog::sink_ptr> sinks;
410 sinks.emplace_back(std::make_shared<spdlog::sinks::ostream_sink_mt>(os,
false));
415 const std::vector<spdlog::sink_ptr> &sinks,
416 const spdlog::level::level_enum log_level)
418 set_adjoint_logger(std::make_shared<spdlog::logger>(
"adjoint-polyfem", sinks.begin(), sinks.end()));
434 std::filesystem::path root = std::filesystem::current_path();
436 root = p_args_in[
"root_path"].get<std::string>();
437 resources_ = std::make_unique<io::FileSystemIO>(root);
442 const json &p_args_in,
444 const bool strict_validation)
448 json args_in = p_args_in;
452 std::filesystem::path output_path =
args[
"output"][
"directory"].get<std::string>();
453 if (!output_path.is_absolute())
454 output_path =
resources_->host_directory() / output_path;
455 const std::string
output_dir = output_path.lexically_normal().string();
458 std::filesystem::create_directories(
output_dir);
462 std::string out_path_log =
args[
"output"][
"log"][
"path"];
463 if (!out_path_log.empty())
465 out_path_log = (std::filesystem::path(out_path_log).is_absolute()
466 ? std::filesystem::path(out_path_log)
467 : std::filesystem::path(
output_dir) / out_path_log)
474 args[
"output"][
"log"][
"level"],
475 args[
"output"][
"log"][
"file_level"],
476 args[
"output"][
"log"][
"quiet"]);
480 const int thread_in =
args[
"solver"][
"max_threads"];
489 size_t threads = max_threads <= 0
490 ? std::numeric_limits<unsigned int>::max()
496 if (!
args[
"output"][
"log"].empty())
497 cur_args[
"output"][
"log"].merge_patch(
args[
"output"][
"log"]);
504 diff_cache = std::make_shared<DiffCache>();
514 for (
int i = 0; i <
varforms.size(); ++i)
520 "varform::DifferentiableVarForm {} ({}) does not expose solve data required by optimization.",
528 "varform::DifferentiableVarForm {}: transient linear problem is not supported in optimization.", i);
534 if (!varform.
get_args()[
"contact"][
"use_gcp_formulation"].get<
bool>()
535 && !varform.
get_args()[
"contact"][
"use_convergent_formulation"].get<
bool>())
538 "varform::DifferentiableVarForm {}: non-convergent contact formulation is not supported in optimization.", i);
542 if (varform.
get_args()[
"/solver/contact/barrier_stiffness"_json_pointer].is_string())
545 "varform::DifferentiableVarForm {}: only constant barrier stiffness is supported in optimization.", i);
550 if (varform.
get_args().contains(
"boundary_conditions") && varform.
get_args()[
"boundary_conditions"].contains(
"rhs"))
552 const json &rhs = varform.
get_args()[
"boundary_conditions"][
"rhs"];
553 if (rhs.is_string() || (rhs.is_array() && rhs.size() > 0 && rhs[0].is_string()))
556 "varform::DifferentiableVarForm {}: only constant rhs over space is supported in optimization.", i);
563 for (
const auto &basis : element_bases.bases)
565 if (basis.order() > 1)
568 "varform::DifferentiableVarForm {}: high-order geometry basis is not supported in optimization.", i);
577 std::vector<std::reference_wrapper<const io::ResourceIO>> resource_refs;
580 resource_refs.emplace_back(*resources);
582 const json ¶meters =
args[
"parameters"];
583 bool is_auto = parameters.is_string() && parameters.get<std::string>() ==
"auto";
589 if (
args[
"variable_to_simulation"].size() != 1)
592 "Auto parameters are only supported with a single variable to simulation.");
597 if (composition[
"type"].get<std::string>() ==
"slice")
615 for (
const auto &arg :
args[
"parameters"])
630 int inv_dof = var2sim->inverse_dof();
634 "VariableToSimulation {} (type {}) expects {} DOF, but parameters define {} DOF.",
635 i, var2sim->name(), inv_dof,
ndof);
642 std::vector<std::reference_wrapper<const io::ResourceIO>> resource_refs;
645 resource_refs.emplace_back(*resources);
652 std::vector<std::shared_ptr<solver::AdjointForm>> stopping_conditions;
653 for (
const auto &arg :
args[
"stopping_conditions"])
654 stopping_conditions.push_back(
657 std::function<bool()> remeshing_trigger;
658 const std::string mode = parse_remeshing_trigger(
args);
659 if (mode ==
"periodic")
661 int period =
args[
"remeshing"][
"trigger"][
"periodic"][
"period"];
663 remeshing_trigger = [
this, period, iter = 0]()
mutable {
668 "Periodic optimization remeshing triggered after {} accepted iteration(s).", iter);
672 else if (mode ==
"scaled_jacobian")
674 int state =
args[
"remeshing"][
"state"];
675 double threshold =
args[
"remeshing"][
"trigger"][
"scaled_jacobian"][
"quality_threshold"];
678 remeshing_trigger = [
this, varform, threshold]() {
679 double quality = minimum_scaled_jacobian(*varform);
681 "Minimum scaled Jacobian: {} (remesh threshold: {}).",
688 nl_problem = std::make_unique<solver::AdjointNLProblem>(
690 args, std::move(remeshing_trigger));
709 args[
"solver"][
"nonlinear"],
710 args[
"solver"][
"linear"],
711 args[
"solver"][
"advanced"][
"characteristic_length"],
715 return nl_solver->status();
void initial_guess(Eigen::VectorXd &x)
void set_log_level(const spdlog::level::level_enum log_level)
change log level
json args
main input arguments containing all defaults
std::vector< std::shared_ptr< varform::DifferentiableVarForm > > varforms
Variational formulations used by the optimization.
void init(const json &args, const bool strict_validation)
initialize the polyfem solver with a json settings
int run(json args, const bool strict_validation)
Run optimization, including remeshing restarts when the selected trigger requests them.
std::vector< std::shared_ptr< DiffCache > > diff_caches
spdlog::sink_ptr file_sink_
void check_unsupported() const
Check and throw if any forward simulation varform::DifferentiableVarForm is not supported.
double eval(Eigen::VectorXd &x) const
std::vector< std::unique_ptr< const io::ResourceIO > > state_resources
std::vector< int > variable_sizes
variables
std::unique_ptr< solver::AdjointNLProblem > nl_problem
std::string output_dir
Directory for output files.
std::vector< json > state_args
spdlog::sink_ptr console_sink_
logger sink to stdout
void init_logger(const std::string &log_file, const spdlog::level::level_enum log_level, const spdlog::level::level_enum file_log_level, const bool is_quiet)
initializing the logger
std::unique_ptr< const io::ResourceIO > resources_
solver::VariableToSimulationGroup variable_to_simulations
void create_varforms(const int max_threads=-1)
Create the optimization variational formulations.
polysolve::nonlinear::Status solve(Eigen::VectorXd &x)
void init_variables()
init variables
bool remeshing_requested_
virtual bool is_time_dependent() const
static void write(const std::string &path, const mesh::Mesh &mesh, const bool binary)
saves the mesh
Read-only interface to a hierarchical collection of simulation resources.
virtual std::unique_ptr< const ResourceIO > with_root(const std::string &root) const =0
void update(const Eigen::VectorXd &x)
std::vector< std::shared_ptr< VariableToSimulation > > data
void set_logger(spdlog::logger &logger)
static GeogramUtils & instance()
void set_num_threads(const int max_threads)
std::shared_ptr< varform::DifferentiableVarForm > build_differentiable_varform(const json &args, const io::ResourceIO &resources, const size_t max_threads)
solver::VariableToSimulationGroup build_variable_to_simulation_group(const json &args, const std::vector< std::shared_ptr< varform::DifferentiableVarForm > > &varforms, const std::vector< std::shared_ptr< DiffCache > > &diff_caches, const std::vector< int > &variable_sizes, const std::vector< std::reference_wrapper< const io::ResourceIO > > &resources)
std::shared_ptr< solver::AdjointForm > build_form(const json &args, const solver::VariableToSimulationGroup &var2sim, const std::vector< std::shared_ptr< varform::DifferentiableVarForm > > &varforms, const std::vector< std::shared_ptr< DiffCache > > &diff_caches, const std::vector< std::reference_wrapper< const io::ResourceIO > > &resources)
bool scaled_jacobian(Mesh3DStorage &hmi, Mesh_Quality &mq)
NLOHMANN_JSON_SERIALIZE_ENUM(CollisionProxyTessellation, {{CollisionProxyTessellation::REGULAR, "regular"}, {CollisionProxyTessellation::IRREGULAR, "irregular"}})
std::vector< T > json_as_array(const json &j)
Return the value of a json object as an array.
bool is_param_valid(const json ¶ms, const std::string &key)
Determine if a key exists and is non-null in a json object.
spdlog::logger & logger()
Retrieves the current logger.
spdlog::logger & adjoint_logger()
Retrieves the current logger for adjoint.
void set_adjoint_logger(std::shared_ptr< spdlog::logger > p_logger)
Setup a logger object to be used by adjoint Polyfem.
void log_and_throw_adjoint_error(const std::string &msg)
void log_and_throw_error(const std::string &msg)
static int compute_variable_size(const json &args, const std::vector< std::shared_ptr< varform::DifferentiableVarForm > > &varforms)
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 std::shared_ptr< polysolve::nonlinear::Solver > make_nl_solver(const json &solver_params, const json &linear_solver_params, const double characteristic_length, const bool strict_validation=true)