31 using namespace solver;
32 using namespace time_integrator;
34 using namespace utils;
48 std::unique_ptr<EnergyCSVWriter> energy_csv =
nullptr;
49 std::unique_ptr<RuntimeStatsCSVWriter> stats_csv =
nullptr;
51 if (
args[
"output"][
"stats"])
55 stats_csv = std::make_unique<RuntimeStatsCSVWriter>(
resolve_output_path(
"stats.csv"), *
this, t0, dt);
57 const bool remesh_enabled =
args[
"space"][
"remesh"][
"enabled"];
62 energy_csv->write(save_i, sol);
69 for (
int t = 1; t <= time_steps; ++t)
71 double forward_solve_time = 0, remeshing_time = 0, global_relaxation_time = 0;
81 energy_csv->write(save_i, sol);
88 remesh_success = this->
remesh(t0 + dt * t, dt, sol);
93 energy_csv->write(save_i, sol);
107 energy_csv->write(save_i, sol);
108 save_timestep(t0 + dt * t, t, t0, dt, sol, Eigen::MatrixXd());
127 logger().info(
"{}/{} t={}", t, time_steps, t0 + dt * t);
129 const std::string rest_mesh_path =
args[
"output"][
"data"][
"rest_mesh"].get<std::string>();
130 if (!rest_mesh_path.empty())
137 V,
F,
mesh->get_body_ids(),
mesh->is_volume(),
true);
141 if (!state_path.empty())
146 if (remesh_enabled && stats_csv)
147 stats_csv->write(t, forward_solve_time, remeshing_time, global_relaxation_time, sol);
153 assert(sol.cols() == 1);
172 const Eigen::MatrixXd displaced =
collision_mesh.displace_vertices(
175 if (ipc::has_intersections(
collision_mesh, displaced,
args[
"solver"][
"contact"][
"CCD"][
"broad_phase"]))
186 if (
problem->is_time_dependent())
188 if (init_time_integrator)
193 Eigen::MatrixXd solution, velocity, acceleration;
195 solution.col(0) = sol;
196 assert(solution.rows() == sol.size());
198 assert(velocity.rows() == sol.size());
200 assert(acceleration.rows() == sol.size());
210 const double dt =
args[
"time"][
"dt"];
246 args.value(
"/time/quasistatic"_json_pointer,
true),
mass,
249 args[
"solver"][
"advanced"][
"lagged_regularization_weight"],
250 args[
"solver"][
"advanced"][
"lagged_regularization_iterations"],
255 avg_mass,
args[
"contact"][
"use_convergent_formulation"] ? bool(
args[
"contact"][
"use_area_weighting"]) :
false,
256 args[
"contact"][
"use_convergent_formulation"] ? bool(
args[
"contact"][
"use_improved_max_operator"]) :
false,
257 args[
"contact"][
"use_convergent_formulation"] ? bool(
args[
"contact"][
"use_physical_barrier"]) :
false,
258 args[
"solver"][
"contact"][
"barrier_stiffness"],
259 args[
"solver"][
"contact"][
"CCD"][
"broad_phase"],
260 args[
"solver"][
"contact"][
"CCD"][
"tolerance"],
261 args[
"solver"][
"contact"][
"CCD"][
"max_iterations"],
264 args[
"contact"][
"adhesion"][
"adhesion_enabled"],
265 args[
"contact"][
"adhesion"][
"dhat_p"],
266 args[
"contact"][
"adhesion"][
"dhat_a"],
267 args[
"contact"][
"adhesion"][
"adhesion_strength"],
269 args[
"contact"][
"adhesion"][
"tangential_adhesion_coefficient"],
270 args[
"contact"][
"adhesion"][
"epsa"],
271 args[
"solver"][
"contact"][
"tangential_adhesion_iterations"],
277 args[
"contact"][
"friction_coefficient"],
278 args[
"contact"][
"epsv"],
279 args[
"solver"][
"contact"][
"friction_iterations"],
281 args[
"solver"][
"rayleigh_damping"]);
283 for (
const auto &form : forms)
295 polysolve::linear::Solver::create(
args[
"solver"][
"linear"],
logger()));
308 assert(sol.size() ==
rhs.size());
317 logger().info(
"Lagging iteration 1:");
323 int subsolve_count = 0;
328 std::shared_ptr<polysolve::nonlinear::Solver> nl_solver =
make_nl_solver(
true);
332 args[
"solver"][
"augmented_lagrangian"][
"initial_weight"],
333 args[
"solver"][
"augmented_lagrangian"][
"scaling"],
334 args[
"solver"][
"augmented_lagrangian"][
"max_weight"],
335 args[
"solver"][
"augmented_lagrangian"][
"eta"],
336 [&](
const Eigen::VectorXd &
x) {
337 this->solve_data.update_barrier_stiffness(sol);
342 {{
"type", al_weight > 0 ?
"al" :
"rc"},
344 {
"info", nl_solver->info()}});
350 Eigen::MatrixXd prev_sol = sol;
357 if (
args[
"space"][
"advanced"][
"count_flipped_els_continuous"])
360 logger().debug(
"Flipped elements (cnt {}) : {}", invalidList.size(), invalidList);
372 for (
int lag_i = 1; !lagging_converged; lag_i++)
380 Eigen::VectorXd grad;
382 const double delta_x_norm = (prev_sol - sol).lpNorm<Eigen::Infinity>();
383 logger().debug(
"Lagging convergence grad_norm={:g} tol={:g} (||Δx||={:g})", grad.norm(), lagging_tol, delta_x_norm);
384 if (grad.norm() <= lagging_tol)
387 "Lagging converged in {:d} iteration(s) (grad_norm={:g} tol={:g})",
388 lag_i, grad.norm(), lagging_tol);
389 lagging_converged =
true;
393 if (delta_x_norm <= 1e-12)
396 "Lagging produced tiny update between iterations {:d} and {:d} (grad_norm={:g} grad_tol={:g} ||Δx||={:g} Δx_tol={:g}); stopping early",
397 lag_i - 1, lag_i, grad.norm(), lagging_tol, delta_x_norm, 1e-6);
398 lagging_converged =
false;
406 "Lagging failed to converge with {:d} iteration(s) (grad_norm={:g} tol={:g})",
407 lag_i, grad.norm(), lagging_tol);
408 lagging_converged =
false;
413 logger().info(
"Lagging iteration {:d}:", lag_i + 1);
414 nl_problem.
init(sol);
417 nl_solver->minimize(nl_problem, tmp_sol);
427 {
"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
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::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
bool remesh(const double time, const double dt, Eigen::MatrixXd &sol)
Remesh the FE space and update solution(s).
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 have 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.
bool is_problem_linear() const
Returns whether the system is linear. Collisions and pressure add nonlinearity to the problem.
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
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 solve_reduced(NLProblem &nl_problem, Eigen::MatrixXd &sol, std::shared_ptr< polysolve::nonlinear::Solver > nl_solver)
std::function< void(const double)> post_subsolve
void solve_al(NLProblem &nl_problem, Eigen::MatrixXd &sol, std::shared_ptr< polysolve::nonlinear::Solver > nl_solver)
bool uses_lagging() const
int max_lagging_iterations() const
virtual void init(const TVector &x0) override
double normalize_forms() override
TVector full_to_reduced(const TVector &full) const
virtual void gradient(const TVector &x, TVector &gradv) override
void init_lagging(const TVector &x) override
TVector reduced_to_full(const TVector &reduced) const
void update_lagging(const TVector &x, const int iter_num) override
std::vector< std::shared_ptr< Form > > init_forms(const Units &units, const int dim, const double t, const Eigen::VectorXi &in_node_to_node, const int n_bases, std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &geom_bases, const assembler::Assembler &assembler, assembler::AssemblyValsCache &ass_vals_cache, const assembler::AssemblyValsCache &mass_ass_vals_cache, const double jacobian_threshold, const solver::ElementInversionCheck check_inversion, 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 std::vector< std::string > &hard_constraint_files, const std::vector< json > &soft_constraint_files, const bool contact_enabled, const ipc::CollisionMesh &collision_mesh, const double dhat, const double avg_mass, const bool use_area_weighting, const bool use_improved_max_operator, const bool use_physical_barrier, const json &barrier_stiffness, const ipc::BroadPhaseMethod broad_phase, const double ccd_tolerance, const long ccd_max_iterations, const bool enable_shape_derivatives, const bool adhesion_enabled, const double dhat_p, const double dhat_a, const double Y, const double tangential_adhesion_coefficient, const double epsa, const int tangential_adhesion_iterations, const assembler::MacroStrainValue ¯o_strain_constraint, const bool periodic_contact, const Eigen::VectorXi &tiled_to_single, const std::shared_ptr< utils::PeriodicBoundary > &periodic_bc, 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.
void update_dt()
updates the dt inside the different forms
std::shared_ptr< solver::NLProblem > nl_problem
std::shared_ptr< solver::ContactForm > contact_form
std::vector< std::shared_ptr< solver::AugmentedLagrangianForm > > al_form
std::shared_ptr< time_integrator::ImplicitTimeIntegrator > time_integrator
void update_barrier_stiffness(const Eigen::VectorXd &x)
update the barrier stiffness for the forms
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.
std::vector< int > count_invalid(const int dim, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const Eigen::VectorXd &u)
spdlog::logger & logger()
Retrieves the current logger.
void log_and_throw_error(const std::string &msg)