9#include <igl/boundary_facets.h>
10#include <igl/writeOBJ.h>
30 vector<list<int>>
adj;
33 void topologicalSortUtil(
int v, vector<bool> &visited, stack<int> &Stack);
39 void addEdge(
int v,
int w);
42 vector<int> topologicalSort();
51 void Graph::addEdge(
int v,
int w)
57 void Graph::topologicalSortUtil(
int v, vector<bool> &visited,
64 list<int>::iterator i;
65 for (i =
adj[v].begin(); i !=
adj[v].end(); ++i)
67 topologicalSortUtil(*i, visited, Stack);
75 vector<int> Graph::topologicalSort()
80 vector<bool> visited(
V,
false);
84 for (
int i = 0; i <
V; i++)
85 if (visited[i] ==
false)
86 topologicalSortUtil(i, visited, Stack);
90 while (Stack.empty() ==
false)
92 sorted.push_back(Stack.top());
103 variables_to_simulation_(variables_to_simulation),
104 all_states_(all_states),
105 save_freq(args[
"output"][
"save_frequency"]),
106 solve_in_parallel(args[
"solver"][
"advanced"][
"solve_in_parallel"])
110 if (args[
"output"][
"solution"] !=
"")
112 solution_ostream.open(args[
"output"][
"solution"].get<std::string>(), std::ofstream::out);
113 if (!solution_ostream.is_open())
117 solve_in_order.clear();
119 Graph G(all_states.size());
120 for (
int k = 0; k < all_states.size(); k++)
122 auto &arg = args[
"states"][k];
123 if (arg[
"initial_guess"].get<int>() >= 0)
124 G.addEdge(arg[
"initial_guess"].get<int>(), k);
127 solve_in_order = G.topologicalSort();
130 active_state_mask.assign(all_states_.size(),
false);
131 for (
int i = 0; i < all_states_.size(); i++)
133 for (
const auto &v2sim : variables_to_simulation_)
135 for (
const auto &state : v2sim->get_states())
137 if (all_states_[i].get() == state.get())
139 active_state_mask[i] =
true;
168 gradv.setZero(
x.size());
178 form_->first_derivative(
x, gradv);
187 return form_->is_step_valid(x0, x1);
192 return form_->is_step_collision_free(x0, x1);
197 return form_->max_step_size(x0, x1);
202 form_->line_search_begin(x0, x1);
207 form_->line_search_end();
214 form_->post_step(data);
223 adjoint_logger().debug(
"Save solution at iteration {} to file...", iter_num);
224 solution_ostream << iter_num <<
": " << std::setprecision(16) << x0.transpose() << std::endl;
233 bool save_vtu =
true;
234 bool save_rest_mesh =
true;
236 std::string vis_mesh_path = state->resolve_output_path(fmt::format(
"opt_state_{:d}_iter_{:d}.vtu",
id, iter_num));
237 std::string rest_mesh_path = state->resolve_output_path(fmt::format(
"opt_state_{:d}_iter_{:d}.obj",
id, iter_num));
242 adjoint_logger().debug(
"Save final vtu to file {} ...", vis_mesh_path);
244 double tend = state->args.value(
"tend", 1.0);
246 if (!state->args[
"time"].is_null())
247 dt = state->args[
"time"][
"dt"];
249 Eigen::MatrixXd sol = state->diff_cached.u(-1);
251 state->out_geom.save_vtu(
255 Eigen::MatrixXd::Zero(state->n_pressure_bases, 1),
258 state->is_contact_enabled(),
259 state->solution_frames);
263 adjoint_logger().debug(
"Save rest mesh to file {} ...", rest_mesh_path);
268 state->get_vertices(
V);
269 state->get_elements(
F);
270 if (state->mesh->dimension() == 3)
271 F = igl::boundary_facets<Eigen::MatrixXi, Eigen::MatrixXi>(
F);
279 bool need_rebuild_basis =
false;
286 need_rebuild_basis =
true;
289 if (need_rebuild_basis)
292 state->build_basis();
298 form_->solution_changed(newX);
308 for (int i = start; i < end; i++)
310 auto state = all_states_[i];
311 if (active_state_mask[i] || state->diff_cached.size() == 0)
313 state->assemble_rhs();
314 state->assemble_mass_mat();
315 Eigen::MatrixXd sol, pressure;
316 state->solve_problem(sol, pressure);
325 Eigen::MatrixXd sol, pressure;
326 for (
int i : solve_in_order)
328 auto state = all_states_[i];
329 if (active_state_mask[i] || state->diff_cached.size() == 0)
331 state->assemble_rhs();
332 state->assemble_mass_mat();
334 state->solve_problem(sol, pressure);
349 obj->solution_changed(
x);
350 if (obj->value(
x) > 0)
vector< list< int > > adj
#define POLYFEM_SCOPED_TIMER(...)
static bool write(const std::string &path, const Eigen::MatrixXd &v, const Eigen::MatrixXi &e, const Eigen::MatrixXi &f)
double max_step_size(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) override
void gradient(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) override
std::ofstream solution_ostream
void save_to_file(const int iter_num, const Eigen::VectorXd &x0)
bool stop(const TVector &x) override
bool is_step_valid(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) override
std::vector< std::shared_ptr< AdjointForm > > stopping_conditions_
void hessian(const Eigen::VectorXd &x, StiffnessMatrix &hessian) override
void post_step(const polysolve::nonlinear::PostStepData &data) override
AdjointNLProblem(std::shared_ptr< AdjointForm > form, const VariableToSimulationGroup &variables_to_simulation, const std::vector< std::shared_ptr< State > > &all_states, const json &args)
void line_search_end() override
void solution_changed(const Eigen::VectorXd &new_x) override
bool is_step_collision_free(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) override
double value(const Eigen::VectorXd &x) override
std::vector< std::shared_ptr< State > > all_states_
std::shared_ptr< AdjointForm > form_
const VariableToSimulationGroup variables_to_simulation_
void line_search_begin(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) override
const bool solve_in_parallel
A collection of VariableToSimulation.
void maybe_parallel_for(int size, const std::function< void(int, int, int)> &partial_for)
spdlog::logger & adjoint_logger()
Retrieves the current logger for adjoint.
void log_and_throw_adjoint_error(const std::string &msg)
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix