8 class QuadraticBarrier :
public ipc::Barrier
11 QuadraticBarrier(
const double weight = 1) :
weight_(weight) {}
13 double operator()(
const double d,
const double dhat)
const override
18 return weight_ * (d - dhat) * (d - dhat);
20 double first_derivative(
const double d,
const double dhat)
const override
25 return 2 *
weight_ * (d - dhat);
27 double second_derivative(
const double d,
const double dhat)
const override
34 double units(
const double dhat)
const override
46 :
AdjointForm(variable_to_simulation), state_(state), dhat_(dhat), dmin_(dmin), barrier_potential_(dhat)
51 [
this](
const std::string &p) { return this->state_.resolve_input_path(p); },
90 if ((
V1 -
V0).lpNorm<Eigen::Infinity>() == 0.0)
93 const ipc::TightInclusionCCD tight_inclusion_ccd(1e-6, 1e6);
94 bool is_valid = ipc::is_step_collision_free(
100 tight_inclusion_ccd);
110 const ipc::TightInclusionCCD tight_inclusion_ccd(1e-6, 1e6);
111 double max_step = ipc::compute_collision_free_stepsize(
117 tight_inclusion_ccd);
126 static Eigen::MatrixXd cached_displaced_surface;
127 if (cached_displaced_surface.size() == displaced_surface.size() && cached_displaced_surface == displaced_surface)
132 cached_displaced_surface = displaced_surface;
137 Eigen::VectorXd X =
X_init;
143 :
AdjointForm(variable_to_simulation), state_(state), dhat_(dhat), barrier_potential_(dhat)
151 [
this](
const std::string &p) { return this->state_.resolve_input_path(p); },
220 static Eigen::MatrixXd cached_displaced_surface;
221 if (cached_displaced_surface.size() == displaced_surface.size() && cached_displaced_surface == displaced_surface)
226 cached_displaced_surface = displaced_surface;
231 Eigen::VectorXd X =
X_init;
238 const std::vector<int> &boundary_ids,
240 const bool use_log_barrier,
243 boundary_ids_(boundary_ids)
250 if (!use_log_barrier)
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);
262 std::vector<bool> is_on_surface;
263 is_on_surface.resize(node_positions.rows(),
false);
266 Eigen::MatrixXd points, uv, normals;
267 Eigen::VectorXd weights;
268 Eigen::VectorXi global_primitive_ids;
271 const int e = lb.element_id();
281 for (
int i = 0; i < lb.size(); ++i)
283 const int primitive_global_id = lb.global_primitive_id(i);
285 const int boundary_id =
state_.
mesh->get_boundary_id(primitive_global_id);
290 for (
long n = 0; n < nodes.size(); ++n)
293 is_on_surface[v.
global[0].index] =
true;
294 assert(v.
global[0].index < node_positions.rows());
300 Eigen::SparseMatrix<double> displacement_map;
301 if (!displacement_map_entries.empty())
304 displacement_map.setFromTriplets(displacement_map_entries.begin(), displacement_map_entries.end());
308 Eigen::MatrixXi boundary_edges_alt(0, 2), boundary_triangles_alt(0, 3);
310 for (
int i = 0; i < boundary_edges.rows(); ++i)
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)];
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);
324 for (
int i = 0; i < boundary_triangles.rows(); ++i)
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)];
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);
337 boundary_triangles_alt.resize(0, 0);
343 boundary_triangles_alt,
350 if (!is_on_surface[dof_idx_i])
355 if (!is_on_surface[dof_idx_j])
358 bool collision_allowed =
true;
361 collision_allowed =
false;
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.
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
int n_geom_bases
number of geometric bases
std::vector< mesh::LocalBoundary > total_local_boundary
mapping from elements to nodes for all mesh
void build_collision_mesh()
extracts the boundary mesh for collision, called in build_basis
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
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)