6#include <wmtk/ExecutionScheduler.hpp>
7#include <wmtk/utils/TriQualityUtils.hpp>
8#include <wmtk/utils/TetraQualityUtils.hpp>
9#include <wmtk/utils/TupleUtils.hpp>
10#include <wmtk/utils/AMIPS.h>
11#include <wmtk/utils/AMIPS2D.h>
15 template <
class WMTKMesh>
23 std::array<double, NDOF> T;
24 for (
auto i = 0; i < WildRemesher<WMTKMesh>::VERTICES_PER_ELEMENT; i++)
25 for (
auto j = 0; j < WildRemesher<WMTKMesh>::DIM; j++)
31 if constexpr (std::is_same_v<wmtk::TriMesh, WMTKMesh>)
32 energy = wmtk::AMIPS2D_energy(T);
34 energy = wmtk::AMIPS_energy(T);
37 if (!std::isfinite(energy))
38 return std::numeric_limits<double>::max();
43 template <
class WMTKMesh>
51 const std::vector<basis::ElementBases> from_bases =
53 const std::vector<basis::ElementBases> to_bases =
56 const Eigen::MatrixXd &from_projection_quantities = old_local_mesh.projection_quantities();
57 const Eigen::MatrixXd &to_projection_quantities = new_local_mesh.projection_quantities();
60 Eigen::SparseMatrix<double> M, A;
68 no_density, to_bases, to_bases,
cache, M);
69 assert(M.rows() == new_local_mesh.num_vertices() * m.
dim());
73 old_local_mesh.num_vertices(), from_bases, from_bases,
74 new_local_mesh.num_vertices(), to_bases, to_bases,
76 assert(A.rows() == new_local_mesh.num_vertices() * m.
dim());
77 assert(A.cols() == old_local_mesh.num_vertices() * m.
dim());
83 ipc::CollisionMesh collision_mesh;
97 std::vector<int> boundary_nodes;
98 for (
int d = 0; d < m.
dim(); ++d)
101 for (
const int fv : new_local_mesh.fixed_vertices())
102 boundary_nodes.push_back(fv * m.dim() + d);
104 for (
const int fv : new_local_mesh.boundary_facets().reshaped())
105 boundary_nodes.push_back(fv * m.dim() + d);
107 std::sort(boundary_nodes.begin(), boundary_nodes.end());
108 auto new_end = std::unique(boundary_nodes.begin(), boundary_nodes.end());
109 boundary_nodes.erase(new_end, boundary_nodes.end());
113 Eigen::MatrixXd projected_quantities = to_projection_quantities;
114 const int n_constrained_quantaties = projected_quantities.cols() / 3;
115 const int n_unconstrained_quantaties = projected_quantities.cols() - n_constrained_quantaties;
118 for (
int i = 0; i < n_constrained_quantaties; ++i)
120 const auto level_before =
logger().level();
121 logger().set_level(spdlog::level::warn);
125 M, A, from_projection_quantities.col(i),
127 new_local_mesh.rest_positions(), new_local_mesh.elements(), m.
dim(),
129 collision_mesh, m.
state.
args[
"contact"][
"dhat"],
133 m.state.args[
"contact"][
"use_convergent_formulation"],
134 m.state.args[
"solver"][
"contact"][
"CCD"][
"broad_phase"],
135 m.state.args[
"solver"][
"contact"][
"CCD"][
"tolerance"],
136 m.state.args[
"solver"][
"contact"][
"CCD"][
"max_iterations"],
138 boundary_nodes, 0, to_projection_quantities.col(i),
140 to_projection_quantities.col(i));
141 logger().set_level(level_before);
146 M, A, from_projection_quantities.rightCols(n_unconstrained_quantaties),
147 boundary_nodes, projected_quantities.rightCols(n_unconstrained_quantaties));
151 assert(projected_quantities.rows() == m.
dim() * new_local_mesh.num_vertices());
152 for (
int i = 0; i < new_local_mesh.num_vertices(); i++)
154 m.
vertex_attrs[new_local_mesh.local_to_global()[i]].projection_quantities =
155 projected_quantities.middleRows(m.
dim() * i, m.
dim());
161 template <
class WMTKMesh>
164 if (!WMTKMesh::smooth_before(v))
166 return !vertex_attrs[v.vid(*
this)].fixed;
169 template <
class WMTKMesh>
172 if (!Super::smooth_before(v))
175 if (this->op_cache ==
nullptr)
177 if constexpr (std::is_same_v<WMTKMesh, wmtk::TriMesh>)
178 this->op_cache = std::make_shared<TriOperationCache>();
180 this->op_cache = std::make_shared<TetOperationCache>();
183 this->op_cache->local_energy = local_mesh_energy(
184 vertex_attrs[v.vid(*
this)].rest_position);
191 template <
class WMTKMesh>
194 if (!WMTKMesh::smooth_before(v))
197 const size_t vid = v.vid(*
this);
199 const std::vector<Tuple> one_ring = this->get_one_ring_elements_for_vertex(v);
200 assert(one_ring.size() > 0);
207 double max_quality = 0;
208 for (
const Tuple &t : get_one_ring_elements_for_vertex(v))
209 max_quality = std::max(max_quality,
get_quality(*
this, t));
210 assert(max_quality > 0);
213 constexpr int NDOF = DIM * VERTICES_PER_ELEMENT;
214 std::vector<std::array<double, NDOF>> assembles(one_ring.size());
218 for (
int i = 0; i < one_ring.size(); i++)
220 const Tuple &t = one_ring[i];
221 auto t_id = t.fid(*
this);
223 assert(!is_inverted(t));
224 const auto local_verts = orient_preserve_element_reorder(element_vids(t), vid);
226 for (
auto j = 0; j < VERTICES_PER_ELEMENT; j++)
227 for (
auto d = 0; d < DIM; d++)
228 assembles[i][j * DIM + d] = vertex_attrs[local_verts[j]].rest_position[d];
233 const VectorNd old_rest_position = vertex_attrs[vid].rest_position;
236 const auto log_lvl = wmtk::logger().level();
237 wmtk::logger().set_level(spdlog::level::warn);
238 if constexpr (DIM == 2)
239 vertex_attrs[vid].rest_position = wmtk::newton_method_from_stack_2d(
240 assembles, wmtk::AMIPS2D_energy, wmtk::AMIPS2D_jacobian, wmtk::AMIPS2D_hessian);
242 vertex_attrs[vid].rest_position = wmtk::newton_method_from_stack(
243 assembles, wmtk::AMIPS_energy, wmtk::AMIPS_jacobian, wmtk::AMIPS_hessian);
244 wmtk::logger().set_level(log_lvl);
247 if (std::any_of(one_ring.begin(), one_ring.end(), [
this](
const Tuple &t) {
248 return this->is_rest_inverted(t);
259 vertex_attrs[vid].projection_quantities.leftCols(n_quantities() / 3) +=
260 old_rest_position - vertex_attrs[vid].rest_position;
267 if (std::any_of(one_ring.begin(), one_ring.end(), [
this](
const Tuple &t) {
268 return this->is_inverted(t);
278 template <
class WMTKMesh>
281 utils::Timer timer(this->timings[
"Smooth vertex after"]);
283 if (!Super::smooth_after(v))
291 const std::vector<Tuple> one_ring = this->get_one_ring_elements_for_vertex(v);
292 return local_relaxation(one_ring, args[
"smooth"][
"acceptance_tolerance"]);
295 template <
class WMTKMesh>
299 double max_quality = 0;
301 max_quality = std::max(max_quality,
get_quality(m, t));
302 assert(max_quality > 0);
305 executor.should_renew = [](
auto) {
return true; };
306 executor.renew_neighbor_tuples = [](
auto &,
auto,
auto &) ->
Operations {
return {}; };
308 const int max_iters = args[
"smooth"][
"max_iters"];
309 for (
int i = 0; i < max_iters; i++)
312 for (
auto &v : WMTKMesh::get_vertices())
313 smooths.emplace_back(
"vertex_smooth", v);
314 executor(*
this, smooths);
315 if (executor.cnt_success() == 0)
std::unique_ptr< MatrixCache > cache
#define POLYFEM_REMESHER_SCOPED_TIMER(name)
std::string formulation() const
return the formulation (checks if the problem is scalar or not and deals with multiphysics)
std::shared_ptr< polysolve::nonlinear::Solver > make_nl_solver(bool for_al) const
factory to create the nl solver depending on input
json args
main input arguments containing all defaults
solver::SolveData solve_data
timedependent stuff cached
Caches basis evaluation and geometric mapping at every element.
void assemble(const bool is_volume, const int size, const int n_basis, const Density &density, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const AssemblyValsCache &cache, StiffnessMatrix &mass) const
Assembles the mass matrix.
void assemble_cross(const bool is_volume, const int size, const int n_from_basis, const std::vector< basis::ElementBases > &from_bases, const std::vector< basis::ElementBases > &from_gbases, const int n_to_basis, const std::vector< basis::ElementBases > &to_bases, const std::vector< basis::ElementBases > &to_gbases, const AssemblyValsCache &cache, StiffnessMatrix &mass) const
Assembles the cross mass matrix between to function spaces.
bool smooth_before(const Tuple &t) override
bool smooth_after(const Tuple &t) override
typename Super::Tuple Tuple
virtual bool is_volume() const
Is the mesh a volumetric mesh.
const State & state
Reference to the simulation state.
std::array< size_t, VERTICES_PER_ELEMENT > element_vids(const Tuple &t) const
Get the vertex ids of an element.
std::vector< std::pair< std::string, Tuple > > Operations
virtual bool smooth_before(const Tuple &t) override
typename WMTKMesh::Tuple Tuple
wmtk::AttributeCollection< VertexAttributes > vertex_attrs
virtual void smooth_vertices()
std::vector< Tuple > get_one_ring_elements_for_vertex(const Tuple &t) const
Get the one ring of elements around a vertex.
int dim() const override
Dimension of the mesh.
Eigen::Matrix< double, DIM, 1 > VectorNd
virtual bool smooth_after(const Tuple &t) override
std::shared_ptr< solver::ContactForm > contact_form
Eigen::VectorXd constrained_L2_projection(std::shared_ptr< polysolve::nonlinear::Solver > nl_solver, const Eigen::SparseMatrix< double > &M, const Eigen::SparseMatrix< double > &A, const Eigen::VectorXd &y, const Eigen::MatrixXd &rest_positions, const Eigen::MatrixXi &elements, const int dim, const ipc::CollisionMesh &collision_mesh, const double dhat, const double barrier_stiffness, const bool use_convergent_formulation, const ipc::BroadPhaseMethod broad_phase_method, const double ccd_tolerance, const int ccd_max_iterations, const std::vector< int > &boundary_nodes, const size_t obstacle_ndof, const Eigen::VectorXd &target_x, const Eigen::VectorXd &x0)
void project_local_quantities(WildRemesher< WMTKMesh > &m, LocalMesh< WildRemesher< WMTKMesh > > &old_local_mesh, LocalMesh< WildRemesher< WMTKMesh > > &new_local_mesh)
void reduced_L2_projection(const Eigen::MatrixXd &M, const Eigen::MatrixXd &A, const Eigen::Ref< const Eigen::MatrixXd > &y, const std::vector< int > &boundary_nodes, Eigen::Ref< Eigen::MatrixXd > x)
double get_quality(const WildRemesher< WMTKMesh > &m, const typename WMTKMesh::Tuple &t)
spdlog::logger & logger()
Retrieves the current logger.