9 class QuadraticBarrier :
public ipc::Barrier
12 QuadraticBarrier(
const double weight = 1) :
weight_(weight) {}
14 double operator()(
const double d,
const double dhat)
const override
19 return weight_ * (d - dhat) * (d - dhat);
21 double first_derivative(
const double d,
const double dhat)
const override
26 return 2 *
weight_ * (d - dhat);
28 double second_derivative(
const double d,
const double dhat)
const override
35 double units(
const double dhat)
const override
47 :
AdjointForm(variable_to_simulation), state_(state), dhat_(dhat), dmin_(dmin), barrier_potential_(dhat)
52 [
this](
const std::string &p) { return this->state_.resolve_input_path(p); },
91 if ((
V1 -
V0).lpNorm<Eigen::Infinity>() == 0.0)
94 const ipc::TightInclusionCCD tight_inclusion_ccd(1e-6, 1e6);
95 bool is_valid = ipc::is_step_collision_free(
101 tight_inclusion_ccd);
111 const ipc::TightInclusionCCD tight_inclusion_ccd(1e-6, 1e6);
112 double max_step = ipc::compute_collision_free_stepsize(
118 tight_inclusion_ccd);
127 static Eigen::MatrixXd cached_displaced_surface;
128 if (cached_displaced_surface.size() == displaced_surface.size() && cached_displaced_surface == displaced_surface)
133 cached_displaced_surface = displaced_surface;
138 Eigen::VectorXd X =
X_init;
145 const std::vector<int> &boundary_ids,
147 const bool use_log_barrier,
150 boundary_ids_(boundary_ids)
157 if (!use_log_barrier)
163 Eigen::MatrixXd node_positions;
164 Eigen::MatrixXi boundary_edges, boundary_triangles;
165 std::vector<Eigen::Triplet<double>> displacement_map_entries;
167 node_positions, boundary_edges, boundary_triangles, displacement_map_entries);
169 std::vector<bool> is_on_surface;
170 is_on_surface.resize(node_positions.rows(),
false);
173 Eigen::MatrixXd points, uv, normals;
174 Eigen::VectorXd weights;
175 Eigen::VectorXi global_primitive_ids;
178 const int e = lb.element_id();
188 for (
int i = 0; i < lb.size(); ++i)
190 const int primitive_global_id = lb.global_primitive_id(i);
192 const int boundary_id =
state_.
mesh->get_boundary_id(primitive_global_id);
197 for (
long n = 0; n < nodes.size(); ++n)
200 is_on_surface[v.
global[0].index] =
true;
201 assert(v.
global[0].index < node_positions.rows());
207 Eigen::SparseMatrix<double> displacement_map;
208 if (!displacement_map_entries.empty())
211 displacement_map.setFromTriplets(displacement_map_entries.begin(), displacement_map_entries.end());
215 Eigen::MatrixXi boundary_edges_alt(0, 2), boundary_triangles_alt(0, 3);
217 for (
int i = 0; i < boundary_edges.rows(); ++i)
219 bool on_surface =
true;
220 for (
int j = 0; j < boundary_edges.cols(); ++j)
221 on_surface &= is_on_surface[boundary_edges(i, j)];
224 boundary_edges_alt.conservativeResize(boundary_edges_alt.rows() + 1, 2);
225 boundary_edges_alt.row(boundary_edges_alt.rows() - 1) = boundary_edges.row(i);
231 for (
int i = 0; i < boundary_triangles.rows(); ++i)
233 bool on_surface =
true;
234 for (
int j = 0; j < boundary_triangles.cols(); ++j)
235 on_surface &= is_on_surface[boundary_triangles(i, j)];
238 boundary_triangles_alt.conservativeResize(boundary_triangles_alt.rows() + 1, 3);
239 boundary_triangles_alt.row(boundary_triangles_alt.rows() - 1) = boundary_triangles.row(i);
244 boundary_triangles_alt.resize(0, 0);
248 std::vector<bool>(is_on_surface.size(),
false),
251 boundary_triangles_alt,
258 if (!is_on_surface[dof_idx_i])
263 if (!is_on_surface[dof_idx_j])
266 bool collision_allowed =
true;
269 collision_allowed =
false;
282 :
AdjointForm(variable_to_simulation), state_(state), dhat_(dhat), barrier_potential_(dhat)
290 [
this](
const std::string &p) { return this->state_.resolve_input_path(p); },
359 static Eigen::MatrixXd cached_displaced_surface;
360 if (cached_displaced_surface.size() == displaced_surface.size() && cached_displaced_surface == displaced_surface)
365 cached_displaced_surface = displaced_surface;
370 Eigen::VectorXd X =
X_init;
381 params_(state.args[
"contact"][
"dhat"], state.args[
"contact"][
"alpha_t"], 0, state.args[
"contact"][
"alpha_n"], 0, state.mesh->is_volume() ? 2 : 1),
384 auto tmp_ids = args[
"surface_selection"].get<std::vector<int>>();
406 for (
const auto &bs : b.bases)
408 for (
const auto &g : bs.global())
410 is_obstacle[g.index] =
true;
415 collision_mesh_.can_collide = [
this, is_obstacle](
size_t vi,
size_t vj) {
422 ipc::SmoothCollisions collisions;
424 collisions.build(
collision_mesh_, displaced_surface, smooth_contact->get_params(), smooth_contact->using_adaptive_dhat(), smooth_contact->get_broad_phase());
438 Eigen::VectorXd coeff(forces.size());
441 return (coeff.array() * forces.array()).matrix().squaredNorm() / 2;
456 Eigen::VectorXd coeff(forces.size());
459 return weight() * (hessian * (coeff.array() * forces.array()).matrix());
477 Eigen::VectorXd coeff(forces.size());
480 Eigen::VectorXd grads = (hessian * (coeff.array() * forces.array()).matrix());
ElementAssemblyValues vals
main class that contains the polyfem solver and all its state
void get_vertices(Eigen::MatrixXd &vertices) const
int n_bases
number of bases
int n_boundary_samples() const
quadrature used for projecting boundary conditions
const std::vector< basis::ElementBases > & geom_bases() const
Get a constant reference to the geometry mapping bases.
Eigen::VectorXi in_node_to_node
Inpute nodes (including high-order) to polyfem nodes, only for isoparametric.
mesh::Obstacle obstacle
Obstacles used in collisions.
ipc::CollisionMesh collision_mesh
IPC collision mesh.
std::unique_ptr< mesh::Mesh > mesh
current mesh, it can be a Mesh2D or Mesh3D
json args
main input arguments containing all defaults
solver::DiffCache diff_cached
std::vector< basis::ElementBases > bases
FE bases, the size is #elements.
int n_geom_bases
number of geometric bases
std::vector< mesh::LocalBoundary > total_local_boundary
mapping from elements to nodes for all mesh
solver::SolveData solve_data
timedependent stuff cached
void build_collision_mesh()
extracts the boundary mesh for collision, called in build_basis
StiffnessMatrix basis_nodes_to_gbasis_nodes
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
Eigen::VectorXd u(int step) const
std::shared_ptr< solver::ContactForm > contact_form
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)
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)
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.
void log_and_throw_adjoint_error(const std::string &msg)
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix