30 using namespace solver;
31 using namespace time_integrator;
33 using namespace utils;
48 const bool remesh_enabled =
args[
"space"][
"remesh"][
"enabled"];
49#ifndef POLYFEM_WITH_REMESHING
51 log_and_throw_error(
"Remeshing is not enabled in this build! Set POLYFEM_WITH_REMESHING=ON in CMake to enable it.");
56 energy_csv.
write(save_i, sol);
63 for (
int t = 1; t <= time_steps; ++t)
65 double forward_solve_time = 0, remeshing_time = 0, global_relaxation_time = 0;
72#ifdef POLYFEM_WITH_REMESHING
75 energy_csv.
write(save_i, sol);
82 remesh_success = this->remesh(t0 + dt * t, dt, sol);
86 energy_csv.
write(save_i, sol);
99 energy_csv.
write(save_i, sol);
100 save_timestep(t0 + dt * t, t, t0, dt, sol, Eigen::MatrixXd());
119 logger().info(
"{}/{} t={}", t, time_steps, t0 + dt * t);
121 const std::string rest_mesh_path =
args[
"output"][
"data"][
"rest_mesh"].get<std::string>();
122 if (!rest_mesh_path.empty())
129 V,
F,
mesh->get_body_ids(),
mesh->is_volume(),
true);
133 if (!state_path.empty())
138 stats_csv.
write(t, forward_solve_time, remeshing_time, global_relaxation_time, sol);
144 assert(sol.cols() == 1);
163 const Eigen::MatrixXd displaced =
collision_mesh.displace_vertices(
166 if (ipc::has_intersections(
collision_mesh, displaced,
args[
"solver"][
"contact"][
"CCD"][
"broad_phase"]))
177 if (
problem->is_time_dependent())
179 if (init_time_integrator)
184 Eigen::MatrixXd solution, velocity, acceleration;
186 solution.col(0) = sol;
187 assert(solution.rows() == sol.size());
189 assert(velocity.rows() == sol.size());
191 assert(acceleration.rows() == sol.size());
201 const double dt =
args[
"time"][
"dt"];
226 mesh->dimension(), t,
236 args.value(
"/time/quasistatic"_json_pointer,
true),
mass,
239 args[
"solver"][
"advanced"][
"lagged_regularization_weight"],
240 args[
"solver"][
"advanced"][
"lagged_regularization_iterations"],
245 avg_mass,
args[
"contact"][
"use_convergent_formulation"],
246 args[
"solver"][
"contact"][
"barrier_stiffness"],
247 args[
"solver"][
"contact"][
"CCD"][
"broad_phase"],
248 args[
"solver"][
"contact"][
"CCD"][
"tolerance"],
249 args[
"solver"][
"contact"][
"CCD"][
"max_iterations"],
256 args[
"contact"][
"friction_coefficient"],
257 args[
"contact"][
"epsv"],
258 args[
"solver"][
"contact"][
"friction_iterations"],
260 args[
"solver"][
"rayleigh_damping"]);
262 for (
const auto &form : forms)
287 assert(sol.size() ==
rhs.size());
296 logger().info(
"Lagging iteration 1:");
302 int subsolve_count = 0;
307 std::shared_ptr<polysolve::nonlinear::Solver> nl_solver =
make_nl_solver(
true);
311 args[
"solver"][
"augmented_lagrangian"][
"initial_weight"],
312 args[
"solver"][
"augmented_lagrangian"][
"scaling"],
313 args[
"solver"][
"augmented_lagrangian"][
"max_weight"],
314 args[
"solver"][
"augmented_lagrangian"][
"eta"],
315 [&](
const Eigen::VectorXd &
x) {
316 this->solve_data.update_barrier_stiffness(sol);
321 {{
"type", al_weight > 0 ?
"al" :
"rc"},
323 {
"info", nl_solver->info()}});
329 Eigen::MatrixXd prev_sol = sol;
330 al_solver.
solve_al(nl_solver, nl_problem, sol);
344 for (
int lag_i = 1; !lagging_converged; lag_i++)
352 Eigen::VectorXd grad;
354 const double delta_x_norm = (prev_sol - sol).lpNorm<Eigen::Infinity>();
355 logger().debug(
"Lagging convergence grad_norm={:g} tol={:g} (||Δx||={:g})", grad.norm(), lagging_tol, delta_x_norm);
356 if (grad.norm() <= lagging_tol)
359 "Lagging converged in {:d} iteration(s) (grad_norm={:g} tol={:g})",
360 lag_i, grad.norm(), lagging_tol);
361 lagging_converged =
true;
365 if (delta_x_norm <= 1e-12)
368 "Lagging produced tiny update between iterations {:d} and {:d} (grad_norm={:g} grad_tol={:g} ||Δx||={:g} Δx_tol={:g}); stopping early",
369 lag_i - 1, lag_i, grad.norm(), lagging_tol, delta_x_norm, 1e-6);
370 lagging_converged =
false;
378 "Lagging failed to converge with {:d} iteration(s) (grad_norm={:g} tol={:g})",
379 lag_i, grad.norm(), lagging_tol);
380 lagging_converged =
false;
385 logger().info(
"Lagging iteration {:d}:", lag_i + 1);
386 nl_problem.
init(sol);
388 nl_solver->minimize(nl_problem, tmp_sol);
397 {
"info", nl_solver->info()}});
#define POLYFEM_SCOPED_TIMER(...)
Eigen::MatrixXd initial_vel_update
assembler::MacroStrainValue macro_strain_constraint
std::shared_ptr< utils::PeriodicBoundary > periodic_bc
periodic BC and periodic mesh utils
int n_bases
number of bases
void cache_transient_adjoint_quantities(const int current_step, const Eigen::MatrixXd &sol, const Eigen::MatrixXd &disp_grad)
int n_pressure_bases
number of pressure bases
assembler::AssemblyValsCache ass_vals_cache
used to store assembly values for small problems
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.
std::vector< mesh::LocalBoundary > local_pressure_boundary
mapping from elements to nodes for pressure boundary conditions
std::shared_ptr< assembler::ViscousDamping > damping_assembler
mesh::Obstacle obstacle
Obstacles used in collisions.
std::shared_ptr< assembler::Assembler > assembler
assemblers
ipc::CollisionMesh periodic_collision_mesh
IPC collision mesh under periodic BC.
void save_subsolve(const int i, const int t, const Eigen::MatrixXd &sol, const Eigen::MatrixXd &pressure)
saves a subsolve when save_solve_sequence_debug is true
std::string resolve_output_path(const std::string &path) const
Resolve output path relative to output_dir if the path is not absolute.
ipc::CollisionMesh collision_mesh
IPC collision mesh.
std::string output_dir
Directory for output files.
solver::CacheLevel optimization_enabled
void set_materials(std::vector< std::shared_ptr< assembler::Assembler > > &assemblers) const
set the material and the problem dimension
void save_timestep(const double time, const int t, const double t0, const double dt, const Eigen::MatrixXd &sol, const Eigen::MatrixXd &pressure)
saves a timestep
void solve_transient_tensor_nonlinear(const int time_steps, const double t0, const double dt, Eigen::MatrixXd &sol)
solves transient tensor nonlinear problem
void initial_solution(Eigen::MatrixXd &solution) const
Load or compute the initial solution.
std::shared_ptr< polysolve::nonlinear::Solver > make_nl_solver(bool for_al) const
factory to create the nl solver depending on input
std::shared_ptr< assembler::Mass > mass_matrix_assembler
std::unique_ptr< mesh::Mesh > mesh
current mesh, it can be a Mesh2D or Mesh3D
StiffnessMatrix mass
Mass matrix, it is computed only for time dependent problems.
std::shared_ptr< assembler::Problem > problem
current problem, it contains rhs and bc
json args
main input arguments containing all defaults
void save_restart_json(const double t0, const double dt, const int t) const
Save a JSON sim file for restarting the simulation at time t.
void initial_velocity(Eigen::MatrixXd &velocity) const
Load or compute the initial velocity.
io::OutStatsData stats
Other statistics.
Eigen::VectorXi periodic_collision_mesh_to_basis
index mapping from periodic 2x2 collision mesh to FE periodic mesh
void initial_acceleration(Eigen::MatrixXd &acceleration) const
Load or compute the initial acceleration.
std::vector< basis::ElementBases > bases
FE bases, the size is #elements.
std::vector< mesh::LocalBoundary > local_boundary
mapping from elements to nodes for dirichlet boundary conditions
double avg_mass
average system mass, used for contact with IPC
assembler::AssemblyValsCache mass_ass_vals_cache
void solve_tensor_nonlinear(Eigen::MatrixXd &sol, const int t=0, const bool init_lagging=true)
solves nonlinear problems
std::vector< int > boundary_nodes
list of boundary nodes
solver::SolveData solve_data
timedependent stuff cached
void init_nonlinear_tensor_solve(Eigen::MatrixXd &sol, const double t=1.0, const bool init_time_integrator=true)
initialize the nonlinear solver
std::shared_ptr< assembler::ViscousDampingPrev > damping_prev_assembler
bool is_contact_enabled() const
does the simulation has contact
std::vector< mesh::LocalBoundary > local_neumann_boundary
mapping from elements to nodes for neumann boundary conditions
std::shared_ptr< assembler::PressureAssembler > build_pressure_assembler() const
std::shared_ptr< assembler::PressureAssembler > elasticity_pressure_assembler
void build_mesh_matrices(Eigen::MatrixXd &V, Eigen::MatrixXi &F)
Build the mesh matrices (vertices and elements) from the mesh using the bases node ordering.
std::shared_ptr< assembler::MixedAssembler > mixed_assembler
Eigen::MatrixXd initial_sol_update
Eigen::MatrixXd rhs
System right-hand side.
std::unordered_map< int, std::vector< mesh::LocalBoundary > > local_pressure_cavity
mapping from elements to nodes for pressure boundary conditions
double characteristic_length() const
void write(const int i, const Eigen::MatrixXd &sol)
static void write(const std::string &path, const mesh::Mesh &mesh, const bool binary)
saves the mesh
static bool write(const std::string &path, const Eigen::MatrixXd &v, const Eigen::MatrixXi &e, const Eigen::MatrixXi &f)
json solver_info
information of the solver, eg num iteration, time, errors, etc the informations varies depending on t...
void write(const int t, const double forward, const double remeshing, const double global_relaxation, const Eigen::MatrixXd &sol)
std::function< void(const double)> post_subsolve
void solve_reduced(std::shared_ptr< NLSolver > nl_solver, NLProblem &nl_problem, Eigen::MatrixXd &sol)
void solve_al(std::shared_ptr< NLSolver > nl_solver, NLProblem &nl_problem, Eigen::MatrixXd &sol)
bool uses_lagging() const
int max_lagging_iterations() const
virtual void init(const TVector &x0) override
virtual TVector full_to_reduced(const TVector &full) const
virtual void gradient(const TVector &x, TVector &gradv) override
void init_lagging(const TVector &x) override
virtual TVector reduced_to_full(const TVector &reduced) const
void update_lagging(const TVector &x, const int iter_num) override
void update_dt()
updates the dt inside the different forms
std::shared_ptr< solver::BCPenaltyForm > al_pen_form
std::shared_ptr< solver::NLProblem > nl_problem
std::shared_ptr< solver::ContactForm > contact_form
std::vector< std::shared_ptr< Form > > init_forms(const Units &units, const int dim, const double t, const int n_bases, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &geom_bases, const assembler::Assembler &assembler, const assembler::AssemblyValsCache &ass_vals_cache, const assembler::AssemblyValsCache &mass_ass_vals_cache, const int n_pressure_bases, const std::vector< int > &boundary_nodes, const std::vector< mesh::LocalBoundary > &local_boundary, const std::vector< mesh::LocalBoundary > &local_neumann_boundary, const int n_boundary_samples, const Eigen::MatrixXd &rhs, const Eigen::MatrixXd &sol, const assembler::Density &density, const std::vector< mesh::LocalBoundary > &local_pressure_boundary, const std::unordered_map< int, std::vector< mesh::LocalBoundary > > &local_pressure_cavity, const std::shared_ptr< assembler::PressureAssembler > pressure_assembler, const bool ignore_inertia, const StiffnessMatrix &mass, const std::shared_ptr< assembler::ViscousDamping > damping_assembler, const double lagged_regularization_weight, const int lagged_regularization_iterations, const size_t obstacle_ndof, const bool contact_enabled, const ipc::CollisionMesh &collision_mesh, const double dhat, const double avg_mass, const bool use_convergent_contact_formulation, const json &barrier_stiffness, const ipc::BroadPhaseMethod broad_phase, const double ccd_tolerance, const long ccd_max_iterations, const bool enable_shape_derivatives, const assembler::MacroStrainValue ¯o_strain_constraint, const bool periodic_contact, const Eigen::VectorXi &tiled_to_single, const double friction_coefficient, const double epsv, const int friction_iterations, const json &rayleigh_damping)
Initialize the forms and return a vector of pointers to them.
std::shared_ptr< time_integrator::ImplicitTimeIntegrator > time_integrator
void update_barrier_stiffness(const Eigen::VectorXd &x)
update the barrier stiffness for the forms
std::shared_ptr< solver::BCLagrangianForm > al_lagr_form
std::shared_ptr< assembler::RhsAssembler > rhs_assembler
static std::shared_ptr< ImplicitTimeIntegrator > construct_time_integrator(const json ¶ms)
Factory method for constructing implicit time integrators from the name of the integrator.
Eigen::MatrixXd unflatten(const Eigen::VectorXd &x, int dim)
Unflatten rowwises, so every dim elements in x become a row.
spdlog::logger & logger()
Retrieves the current logger.
void log_and_throw_error(const std::string &msg)