Loading [MathJax]/extensions/tex2jax.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
BarrierForms.cpp
Go to the documentation of this file.
1#include "BarrierForms.hpp"
2#include <polyfem/State.hpp>
3
4namespace polyfem::solver
5{
6 namespace
7 {
8 class QuadraticBarrier : public ipc::Barrier
9 {
10 public:
11 QuadraticBarrier(const double weight = 1) : weight_(weight) {}
12
13 double operator()(const double d, const double dhat) const override
14 {
15 if (d > dhat)
16 return 0;
17 else
18 return weight_ * (d - dhat) * (d - dhat);
19 }
20 double first_derivative(const double d, const double dhat) const override
21 {
22 if (d > dhat)
23 return 0;
24 else
25 return 2 * weight_ * (d - dhat);
26 }
27 double second_derivative(const double d, const double dhat) const override
28 {
29 if (d > dhat)
30 return 0;
31 else
32 return 2 * weight_;
33 }
34 double units(const double dhat) const override
35 {
36 return dhat * dhat;
37 }
38
39 private:
40 const double weight_;
41 };
42
43 } // namespace
44
45 CollisionBarrierForm::CollisionBarrierForm(const VariableToSimulationGroup &variable_to_simulation, const State &state, const double dhat, const double dmin)
46 : AdjointForm(variable_to_simulation), state_(state), dhat_(dhat), dmin_(dmin), barrier_potential_(dhat)
47 {
51 [this](const std::string &p) { return this->state_.resolve_input_path(p); },
53
54 Eigen::MatrixXd V;
57
58 broad_phase_method_ = ipc::BroadPhaseMethod::HASH_GRID;
59 }
60
61 double CollisionBarrierForm::value_unweighted(const Eigen::VectorXd &x) const
62 {
63 const Eigen::MatrixXd displaced_surface = collision_mesh_.vertices(utils::unflatten(get_updated_mesh_nodes(x), state_.mesh->dimension()));
64 return barrier_potential_(collision_set, collision_mesh_, displaced_surface);
65 }
66
67 void CollisionBarrierForm::compute_partial_gradient(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const
68 {
70 const Eigen::MatrixXd displaced_surface = collision_mesh_.vertices(utils::unflatten(get_updated_mesh_nodes(x), state_.mesh->dimension()));
71 const Eigen::VectorXd grad = collision_mesh_.to_full_dof(barrier_potential_.gradient(collision_set, collision_mesh_, displaced_surface));
73 });
74 }
75
76 void CollisionBarrierForm::solution_changed(const Eigen::VectorXd &x)
77 {
79
80 const Eigen::MatrixXd displaced_surface = collision_mesh_.vertices(utils::unflatten(get_updated_mesh_nodes(x), state_.mesh->dimension()));
81 build_collision_set(displaced_surface);
82 }
83
84 bool CollisionBarrierForm::is_step_collision_free(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const
85 {
86 const Eigen::MatrixXd V0 = utils::unflatten(get_updated_mesh_nodes(x0), state_.mesh->dimension());
87 const Eigen::MatrixXd V1 = utils::unflatten(get_updated_mesh_nodes(x1), state_.mesh->dimension());
88
89 // Skip CCD if the displacement is zero.
90 if ((V1 - V0).lpNorm<Eigen::Infinity>() == 0.0)
91 return true;
92
93 const ipc::TightInclusionCCD tight_inclusion_ccd(1e-6, 1e6);
94 bool is_valid = ipc::is_step_collision_free(
96 collision_mesh_.vertices(V0),
97 collision_mesh_.vertices(V1),
98 dmin_,
100 tight_inclusion_ccd);
101
102 return is_valid;
103 }
104
105 double CollisionBarrierForm::max_step_size(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const
106 {
107 const Eigen::MatrixXd V0 = utils::unflatten(get_updated_mesh_nodes(x0), state_.mesh->dimension());
108 const Eigen::MatrixXd V1 = utils::unflatten(get_updated_mesh_nodes(x1), state_.mesh->dimension());
109
110 const ipc::TightInclusionCCD tight_inclusion_ccd(1e-6, 1e6);
111 double max_step = ipc::compute_collision_free_stepsize(
113 collision_mesh_.vertices(V0),
114 collision_mesh_.vertices(V1),
115 dmin_,
117 tight_inclusion_ccd);
118
119 adjoint_logger().info("Objective {}: max step size is {}.", name(), max_step);
120
121 return max_step;
122 }
123
124 void CollisionBarrierForm::build_collision_set(const Eigen::MatrixXd &displaced_surface)
125 {
126 static Eigen::MatrixXd cached_displaced_surface;
127 if (cached_displaced_surface.size() == displaced_surface.size() && cached_displaced_surface == displaced_surface)
128 return;
129
130 collision_set.build(collision_mesh_, displaced_surface, dhat_, dmin_, broad_phase_method_);
131
132 cached_displaced_surface = displaced_surface;
133 }
134
135 Eigen::VectorXd CollisionBarrierForm::get_updated_mesh_nodes(const Eigen::VectorXd &x) const
136 {
137 Eigen::VectorXd X = X_init;
140 }
141
142 DeformedCollisionBarrierForm::DeformedCollisionBarrierForm(const VariableToSimulationGroup &variable_to_simulation, const State &state, const double dhat)
143 : AdjointForm(variable_to_simulation), state_(state), dhat_(dhat), barrier_potential_(dhat)
144 {
146 log_and_throw_adjoint_error("[{}] Should use linear FE basis!", name());
147
151 [this](const std::string &p) { return this->state_.resolve_input_path(p); },
153
154 Eigen::MatrixXd V;
157
158 broad_phase_method_ = ipc::BroadPhaseMethod::HASH_GRID;
159 }
160
161 double DeformedCollisionBarrierForm::value_unweighted(const Eigen::VectorXd &x) const
162 {
163 const Eigen::MatrixXd displaced_surface = collision_mesh_.vertices(utils::unflatten(get_updated_mesh_nodes(x), state_.mesh->dimension()));
164
165 return barrier_potential_(collision_set, collision_mesh_, displaced_surface);
166 }
167
168 void DeformedCollisionBarrierForm::compute_partial_gradient(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const
169 {
171 const Eigen::MatrixXd displaced_surface = collision_mesh_.vertices(utils::unflatten(get_updated_mesh_nodes(x), state_.mesh->dimension()));
172 const Eigen::VectorXd grad = collision_mesh_.to_full_dof(barrier_potential_.gradient(collision_set, collision_mesh_, displaced_surface));
174 });
175 }
176
178 {
180
181 const Eigen::MatrixXd displaced_surface = collision_mesh_.vertices(utils::unflatten(get_updated_mesh_nodes(x), state_.mesh->dimension()));
182 build_collision_set(displaced_surface);
183 }
184
185 bool DeformedCollisionBarrierForm::is_step_collision_free(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const
186 {
187 // const Eigen::MatrixXd V0 = utils::unflatten(get_updated_mesh_nodes(x0), state_.mesh->dimension());
188 // const Eigen::MatrixXd V1 = utils::unflatten(get_updated_mesh_nodes(x1), state_.mesh->dimension());
189
190 // // Skip CCD if the displacement is zero.
191 // if ((V1 - V0).lpNorm<Eigen::Infinity>() == 0.0)
192 // return true;
193
194 // bool is_valid = ipc::is_step_collision_free(
195 // collision_mesh_,
196 // collision_mesh_.vertices(V0),
197 // collision_mesh_.vertices(V1),
198 // broad_phase_method_,
199 // 1e-6, 1e6);
200
201 return true; // is_valid;
202 }
203
204 double DeformedCollisionBarrierForm::max_step_size(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const
205 {
206 // const Eigen::MatrixXd V0 = utils::unflatten(get_updated_mesh_nodes(x0), state_.mesh->dimension());
207 // const Eigen::MatrixXd V1 = utils::unflatten(get_updated_mesh_nodes(x1), state_.mesh->dimension());
208
209 // double max_step = ipc::compute_collision_free_stepsize(
210 // collision_mesh_,
211 // collision_mesh_.vertices(V0),
212 // collision_mesh_.vertices(V1),
213 // broad_phase_method_, 1e-6, 1e6);
214
215 return 1; // max_step;
216 }
217
218 void DeformedCollisionBarrierForm::build_collision_set(const Eigen::MatrixXd &displaced_surface)
219 {
220 static Eigen::MatrixXd cached_displaced_surface;
221 if (cached_displaced_surface.size() == displaced_surface.size() && cached_displaced_surface == displaced_surface)
222 return;
223
224 collision_set.build(collision_mesh_, displaced_surface, dhat_, 0, broad_phase_method_);
225
226 cached_displaced_surface = displaced_surface;
227 }
228
235
237 const State &state,
238 const std::vector<int> &boundary_ids,
239 const double dhat,
240 const bool use_log_barrier,
241 const double dmin)
242 : CollisionBarrierForm(variable_to_simulations, state, dhat, dmin),
243 boundary_ids_(boundary_ids)
244 {
245 for (const auto &id : boundary_ids_)
246 boundary_ids_to_dof_[id] = std::set<int>();
247
249
250 if (!use_log_barrier)
251 barrier_potential_.set_barrier(std::make_shared<QuadraticBarrier>());
252 }
253
255 {
256 Eigen::MatrixXd node_positions;
257 Eigen::MatrixXi boundary_edges, boundary_triangles;
258 std::vector<Eigen::Triplet<double>> displacement_map_entries;
260 node_positions, boundary_edges, boundary_triangles, displacement_map_entries);
261
262 std::vector<bool> is_on_surface;
263 is_on_surface.resize(node_positions.rows(), false);
264
266 Eigen::MatrixXd points, uv, normals;
267 Eigen::VectorXd weights;
268 Eigen::VectorXi global_primitive_ids;
269 for (const auto &lb : state_.total_local_boundary)
270 {
271 const int e = lb.element_id();
272 bool has_samples = utils::BoundarySampler::boundary_quadrature(lb, state_.n_boundary_samples(), *state_.mesh, false, uv, points, normals, weights, global_primitive_ids);
273
274 if (!has_samples)
275 continue;
276
277 const basis::ElementBases &gbs = state_.geom_bases()[e];
278
279 vals.compute(e, state_.mesh->is_volume(), points, gbs, gbs);
280
281 for (int i = 0; i < lb.size(); ++i)
282 {
283 const int primitive_global_id = lb.global_primitive_id(i);
284 const auto nodes = gbs.local_nodes_for_primitive(primitive_global_id, *state_.mesh);
285 const int boundary_id = state_.mesh->get_boundary_id(primitive_global_id);
286
287 if (!std::count(boundary_ids_.begin(), boundary_ids_.end(), boundary_id))
288 continue;
289
290 for (long n = 0; n < nodes.size(); ++n)
291 {
292 const assembler::AssemblyValues &v = vals.basis_values[nodes(n)];
293 is_on_surface[v.global[0].index] = true;
294 assert(v.global[0].index < node_positions.rows());
295 boundary_ids_to_dof_[boundary_id].insert(v.global[0].index);
296 }
297 }
298 }
299
300 Eigen::SparseMatrix<double> displacement_map;
301 if (!displacement_map_entries.empty())
302 {
303 displacement_map.resize(node_positions.rows(), state_.n_geom_bases);
304 displacement_map.setFromTriplets(displacement_map_entries.begin(), displacement_map_entries.end());
305 }
306
307 // Fix boundary edges and boundary triangles to exclude vertices not on triangles
308 Eigen::MatrixXi boundary_edges_alt(0, 2), boundary_triangles_alt(0, 3);
309 {
310 for (int i = 0; i < boundary_edges.rows(); ++i)
311 {
312 bool on_surface = true;
313 for (int j = 0; j < boundary_edges.cols(); ++j)
314 on_surface &= is_on_surface[boundary_edges(i, j)];
315 if (on_surface)
316 {
317 boundary_edges_alt.conservativeResize(boundary_edges_alt.rows() + 1, 2);
318 boundary_edges_alt.row(boundary_edges_alt.rows() - 1) = boundary_edges.row(i);
319 }
320 }
321
322 if (state_.mesh->is_volume())
323 {
324 for (int i = 0; i < boundary_triangles.rows(); ++i)
325 {
326 bool on_surface = true;
327 for (int j = 0; j < boundary_triangles.cols(); ++j)
328 on_surface &= is_on_surface[boundary_triangles(i, j)];
329 if (on_surface)
330 {
331 boundary_triangles_alt.conservativeResize(boundary_triangles_alt.rows() + 1, 3);
332 boundary_triangles_alt.row(boundary_triangles_alt.rows() - 1) = boundary_triangles.row(i);
333 }
334 }
335 }
336 else
337 boundary_triangles_alt.resize(0, 0);
338 }
339
340 collision_mesh_ = ipc::CollisionMesh(is_on_surface,
341 node_positions,
342 boundary_edges_alt,
343 boundary_triangles_alt,
344 displacement_map);
345
346 can_collide_cache_.resize(collision_mesh_.num_vertices(), collision_mesh_.num_vertices());
347 for (int i = 0; i < can_collide_cache_.rows(); ++i)
348 {
349 int dof_idx_i = collision_mesh_.to_full_vertex_id(i);
350 if (!is_on_surface[dof_idx_i])
351 continue;
352 for (int j = 0; j < can_collide_cache_.cols(); ++j)
353 {
354 int dof_idx_j = collision_mesh_.to_full_vertex_id(j);
355 if (!is_on_surface[dof_idx_j])
356 continue;
357
358 bool collision_allowed = true;
359 for (const auto &id : boundary_ids_)
360 if (boundary_ids_to_dof_[id].count(dof_idx_i) && boundary_ids_to_dof_[id].count(dof_idx_j))
361 collision_allowed = false;
362 can_collide_cache_(i, j) = collision_allowed;
363 }
364 }
365
366 collision_mesh_.can_collide = [&](size_t vi, size_t vj) {
367 return (bool)can_collide_cache_(vi, vj);
368 };
369
370 collision_mesh_.init_area_jacobians();
371 }
372
373} // namespace polyfem::solver
int V
ElementAssemblyValues vals
Definition Assembler.cpp:22
const double weight_
int x
main class that contains the polyfem solver and all its state
Definition State.hpp:79
void get_vertices(Eigen::MatrixXd &vertices) const
Definition StateDiff.cpp:71
int n_bases
number of bases
Definition State.hpp:178
int n_boundary_samples() const
quadrature used for projecting boundary conditions
Definition State.hpp:263
const std::vector< basis::ElementBases > & geom_bases() const
Get a constant reference to the geometry mapping bases.
Definition State.hpp:223
Eigen::VectorXi in_node_to_node
Inpute nodes (including high-order) to polyfem nodes, only for isoparametric.
Definition State.hpp:455
mesh::Obstacle obstacle
Obstacles used in collisions.
Definition State.hpp:473
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
solver::DiffCache diff_cached
Definition State.hpp:656
int n_geom_bases
number of geometric bases
Definition State.hpp:182
std::vector< mesh::LocalBoundary > total_local_boundary
mapping from elements to nodes for all mesh
Definition State.hpp:436
void build_collision_mesh()
extracts the boundary mesh for collision, called in build_basis
Definition State.cpp:1357
stores per local bases evaluations
std::vector< basis::Local2Global > global
stores per element basis values at given quadrature points and geometric mapping
Stores the basis functions for a given element in a mesh (facet in 2d, cell in 3d).
Eigen::VectorXi local_nodes_for_primitive(const int local_index, const mesh::Mesh &mesh) const
static void extract_boundary_mesh(const mesh::Mesh &mesh, const int n_bases, const std::vector< basis::ElementBases > &bases, const std::vector< mesh::LocalBoundary > &total_local_boundary, Eigen::MatrixXd &node_positions, Eigen::MatrixXi &boundary_edges, Eigen::MatrixXi &boundary_triangles, std::vector< Eigen::Triplet< double > > &displacement_map_entries)
extracts the boundary mesh
Definition OutData.cpp:145
virtual void solution_changed(const Eigen::VectorXd &new_x) override
Update cached fields upon a change in the solution.
const VariableToSimulationGroup variable_to_simulations_
ipc::BarrierPotential barrier_potential_
std::string name() const override
void build_collision_set(const Eigen::MatrixXd &displaced_surface)
ipc::BroadPhaseMethod broad_phase_method_
bool is_step_collision_free(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const override
Checks if the step is collision free.
double max_step_size(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const override
Determine the maximum step size allowable between the current and next solution.
CollisionBarrierForm(const VariableToSimulationGroup &variable_to_simulation, const State &state, const double dhat, const double dmin=0)
void compute_partial_gradient(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const override
void solution_changed(const Eigen::VectorXd &x) override
Update cached fields upon a change in the solution.
double value_unweighted(const Eigen::VectorXd &x) const override
Compute the value of the form.
Eigen::VectorXd get_updated_mesh_nodes(const Eigen::VectorXd &x) const
bool is_step_collision_free(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const override
Checks if the step is collision free.
void compute_partial_gradient(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const override
double max_step_size(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) const override
Determine the maximum step size allowable between the current and next solution.
DeformedCollisionBarrierForm(const VariableToSimulationGroup &variable_to_simulation, const State &state, const double dhat)
double value_unweighted(const Eigen::VectorXd &x) const override
Compute the value of the form.
const ipc::BarrierPotential barrier_potential_
void build_collision_set(const Eigen::MatrixXd &displaced_surface)
Eigen::VectorXd get_updated_mesh_nodes(const Eigen::VectorXd &x) const
void solution_changed(const Eigen::VectorXd &x) override
Update cached fields upon a change in the solution.
Eigen::VectorXd u(int step) const
virtual double weight() const
Get the form's multiplicative constant weight.
Definition Form.hpp:128
LayerThicknessForm(const VariableToSimulationGroup &variable_to_simulations, const State &state, const std::vector< int > &boundary_ids, const double dhat, const bool use_log_barrier=false, const double dmin=0)
std::map< int, std::set< int > > boundary_ids_to_dof_
A collection of VariableToSimulation.
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.
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.
static bool boundary_quadrature(const mesh::LocalBoundary &local_boundary, const int order, const mesh::Mesh &mesh, const bool skip_computation, Eigen::MatrixXd &uv, Eigen::MatrixXd &points, Eigen::MatrixXd &normals, Eigen::VectorXd &weights, Eigen::VectorXi &global_primitive_ids)
Eigen::VectorXd map_primitive_to_node_order(const State &state, const Eigen::VectorXd &primitives)
Eigen::VectorXd map_node_to_primitive_order(const State &state, const Eigen::VectorXd &nodes)
std::tuple< bool, int, Tree > is_valid(const int dim, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const Eigen::VectorXd &u, const double threshold)
Definition Jacobian.cpp:241
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 & adjoint_logger()
Retrieves the current logger for adjoint.
Definition Logger.cpp:30
void log_and_throw_adjoint_error(const std::string &msg)
Definition Logger.cpp:79