Loading [MathJax]/extensions/tex2jax.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
LocalRelaxation.cpp
Go to the documentation of this file.
6
7#include <polysolve/nonlinear/Solver.hpp>
8
9namespace polyfem::mesh
10{
12 decltype(Remesher::timings) &timings,
13 const polyfem::json &solver_info)
14 {
15 // Copy over timing data
16 const int solver_iters = solver_info["iterations"];
17 for (const auto &[key, value] : solver_info.items())
18 {
19 if (utils::StringUtils::startswith(key, "time_")
20 && key != "total_time"
21 && key != "time_line_search")
22 {
23 // The solver reports the time per iteration, so we need to
24 // multiply by the number of iterations.
25 const std::string new_key = "NonlinearSolver::" + key.substr(5, key.size() - 5);
26 timings[new_key] += solver_iters * value.get<double>();
27 }
28 }
29 }
30
31 template <class WMTKMesh>
33 const std::vector<Tuple> &local_mesh_tuples,
34 const double acceptance_tolerance)
35 {
36 // --------------------------------------------------------------------
37 // 1. Get the n-ring of elements around the vertex.
38
39 const bool include_global_boundary =
40 state.is_contact_enabled() && std::any_of(local_mesh_tuples.begin(), local_mesh_tuples.end(), [&](const Tuple &t) {
41 const size_t tid = this->element_id(t);
42 for (int i = 0; i < Super::FACETS_PER_ELEMENT; ++i)
43 if (this->is_boundary_facet(this->tuple_from_facet(tid, i)))
44 return true;
45 return false;
46 });
47
48 LocalMesh<Super> local_mesh(*this, local_mesh_tuples, include_global_boundary);
49
50 // --------------------------------------------------------------------
51 // 2. Perform "relaxation" by minimizing the elastic energy of the
52 // n-ring with the internal boundary edges fixed.
53
54 LocalRelaxationData data(this->state, local_mesh, this->current_time, include_global_boundary);
55 solver::SolveData &solve_data = data.solve_data;
56
57 const int n_free_dof = data.n_free_dof();
58 if (n_free_dof <= 0)
59 {
60 assert(n_free_dof == 0);
61 return false;
62 }
63
64 this->total_ndofs += n_free_dof;
65 this->num_solves++;
66
67 // Nonlinear solver
68 auto nl_solver = state.make_nl_solver(/*for_al=*/false); // TODO: Use Eigen::LLT
69 nl_solver->stop_criteria().iterations = args["local_relaxation"]["max_nl_iterations"];
70 if (this->is_boundary_op())
71 nl_solver->stop_criteria().iterations = std::max(nl_solver->stop_criteria().iterations, size_t(5));
72 nl_solver->allow_out_of_iterations = true;
73
74 Eigen::VectorXd reduced_sol = solve_data.nl_problem->full_to_reduced(data.sol());
75
76 const auto level_before = logger().level();
77 logger().set_level(spdlog::level::warn);
78 try
79 {
80 POLYFEM_REMESHER_SCOPED_TIMER("Local relaxation solve");
81 solve_data.nl_problem->normalize_forms();
82 nl_solver->minimize(*(solve_data.nl_problem), reduced_sol);
83 }
84 catch (const std::runtime_error &e)
85 {
86 logger().set_level(level_before);
87 assert(false);
88 return false;
89 }
90 logger().set_level(level_before);
91
92 // Copy over timing data
93 add_solver_timings(this->timings, nl_solver->info());
94
95 Eigen::VectorXd sol = solve_data.nl_problem->reduced_to_full(reduced_sol);
96
97 // --------------------------------------------------------------------
98 // 3. Determine if we should accept the operation based on a decrease in
99 // energy.
100
101 const double local_energy_after = solve_data.nl_problem->value(sol);
102 assert(std::isfinite(local_energy_before()));
103 assert(std::isfinite(local_energy_after));
104 const double abs_diff = local_energy_before() - local_energy_after; // > 0 if energy decreased
105 // TODO: compute global_energy_before
106 // Right now using: starting_energy = state.solve_data.nl_problem->value(sol)
107 // const double global_energy_before = abs(starting_energy);
108 // const double rel_diff = abs_diff / global_energy_before;
109
110 // NOTE: only use abs_diff
111 // accept = rel_diff >= energy_relative_tolerance && abs_diff >= energy_absolute_tolerance;
112 // NOTE: account for Δt² in energy by multiplying acceptance tol by Δt²
113 const double dt_sqr = solve_data.time_integrator ? solve_data.time_integrator->acceleration_scaling() : 1.0;
114 const bool accept = abs_diff >= dt_sqr * acceptance_tolerance;
115
116 // Update positions only on acceptance
117 if (accept)
118 {
119 static int save_i = 0;
120 // local_mesh.write_mesh(state.resolve_output_path(fmt::format("local_mesh_{:04d}.vtu", save_i)), target_x);
121 // write_mesh(state.resolve_output_path(fmt::format("relaxation_{:04d}.vtu", save_i++)));
122
123 // Re-solve with more iterations
124 if (!is_converged_status(nl_solver->status()))
125 {
126 nl_solver->stop_criteria().iterations = 100;
127
128 const auto level_before = logger().level();
129 logger().set_level(spdlog::level::warn);
130 try
131 {
132 POLYFEM_REMESHER_SCOPED_TIMER("Local relaxation resolve");
133 solve_data.nl_problem->normalize_forms();
134 nl_solver->minimize(*(solve_data.nl_problem), reduced_sol);
135 }
136 catch (const std::runtime_error &e)
137 {
138 logger().set_level(level_before);
139 assert(false);
140 return false;
141 }
142 logger().set_level(level_before);
143
144 // Copy over timing data
145 add_solver_timings(this->timings, nl_solver->info());
146
147 sol = solve_data.nl_problem->reduced_to_full(reduced_sol);
148 }
149
150 for (const auto &[glob_vi, loc_vi] : local_mesh.global_to_local())
151 {
152 const auto u = sol.middleRows(this->dim() * loc_vi, this->dim());
153 const auto u_old = vertex_attrs[glob_vi].displacement();
154 vertex_attrs[glob_vi].position = vertex_attrs[glob_vi].rest_position + u;
155 }
156
157 // local_mesh.write_mesh(state.resolve_output_path(fmt::format("local_mesh_{:04d}.vtu", save_i)), sol);
158 // write_mesh(state.resolve_output_path(fmt::format("relaxation_{:04d}.vtu", save_i++)));
159 }
160
161 static const std::string accept_str =
162 fmt::format(fmt::fg(fmt::terminal_color::green), "accept");
163 static const std::string reject_str =
164 fmt::format(fmt::fg(fmt::terminal_color::yellow), "reject");
165 logger().debug(
166 "[{:s}] E0={:<10g} E1={:<10g} (E0-E1)={:<11g} tol={:g} local_ndof={:d} n_iters={:d}",
167 accept ? accept_str : reject_str, local_energy_before(),
168 local_energy_after, abs_diff, dt_sqr * acceptance_tolerance,
169 n_free_dof, nl_solver->current_criteria().iterations);
170
171 return accept;
172 }
173
174 // ----------------------------------------------------------------------------------------------
175 // Template specializations
176 template class PhysicsRemesher<wmtk::TriMesh>;
177 template class PhysicsRemesher<wmtk::TetMesh>;
178} // namespace polyfem::mesh
#define POLYFEM_REMESHER_SCOPED_TIMER(name)
Definition Remesher.hpp:15
bool local_relaxation(const Tuple &t, const double acceptance_tolerance)
Relax a local n-ring around a vertex.
static std::unordered_map< std::string, utils::Timing > timings
Timings for the remeshing operations.
Definition Remesher.hpp:229
class to store time stepping data
Definition SolveData.hpp:59
std::shared_ptr< solver::NLProblem > nl_problem
std::shared_ptr< time_integrator::ImplicitTimeIntegrator > time_integrator
void add_solver_timings(decltype(Remesher::timings) &timings, const polyfem::json &solver_info)
bool startswith(const std::string &str, const std::string &prefix)
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:42
nlohmann::json json
Definition Common.hpp:9