Loading [MathJax]/extensions/tex2jax.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
State.hpp
Go to the documentation of this file.
1#pragma once
2
3#include <polyfem/Common.hpp>
4
5#include <polyfem/Units.hpp>
6
9
18
19#include <polyfem/mesh/Mesh.hpp>
23
26
32
34
35#include <polysolve/linear/Solver.hpp>
36
37#include <Eigen/Dense>
38#include <Eigen/Sparse>
39
40#include <memory>
41#include <string>
42#include <unordered_map>
43
44#include <spdlog/sinks/basic_file_sink.h>
45
46#include <ipc/collision_mesh.hpp>
47#include <ipc/utils/logger.hpp>
48
49// Forward declaration
51{
52 class Solver;
53}
54
55namespace polyfem::assembler
56{
57 class Mass;
58 class ViscousDamping;
59 class ViscousDampingPrev;
60} // namespace polyfem::assembler
61
62namespace polyfem
63{
64 namespace mesh
65 {
66 class Mesh2D;
67 class Mesh3D;
68 } // namespace mesh
69
70 enum class CacheLevel
71 {
72 None,
75 };
76
78 class State
79 {
80 public:
81 //---------------------------------------------------
82 //-----------------initialization--------------------
83 //---------------------------------------------------
84
85 ~State() = default;
87 State();
88
90 void set_max_threads(const int max_threads = std::numeric_limits<int>::max());
91
95 void init(const json &args, const bool strict_validation);
96
98 void init_time();
99
102
103 //---------------------------------------------------
104 //-----------------logger----------------------------
105 //---------------------------------------------------
106
112 void init_logger(
113 const std::string &log_file,
114 const spdlog::level::level_enum log_level,
115 const spdlog::level::level_enum file_log_level,
116 const bool is_quiet);
117
121 void init_logger(std::ostream &os, const spdlog::level::level_enum log_level);
122
125 void set_log_level(const spdlog::level::level_enum log_level);
126
130 std::string get_log(const Eigen::MatrixXd &sol)
131 {
132 std::stringstream ss;
133 save_json(sol, ss);
134 return ss.str();
135 }
136
137 private:
139 void init_logger(const std::vector<spdlog::sink_ptr> &sinks, const spdlog::level::level_enum log_level);
140
142 spdlog::sink_ptr console_sink_ = nullptr;
143 spdlog::sink_ptr file_sink_ = nullptr;
144
145 public:
146 //---------------------------------------------------
147 //-----------------assembly--------------------------
148 //---------------------------------------------------
149
151
153
155 std::shared_ptr<assembler::Assembler> assembler = nullptr;
156
157 std::shared_ptr<assembler::Mass> mass_matrix_assembler = nullptr;
158
159 std::shared_ptr<assembler::MixedAssembler> mixed_assembler = nullptr;
160 std::shared_ptr<assembler::Assembler> pressure_assembler = nullptr;
161
162 std::shared_ptr<assembler::PressureAssembler> elasticity_pressure_assembler = nullptr;
163
164 std::shared_ptr<assembler::ViscousDamping> damping_assembler = nullptr;
165 std::shared_ptr<assembler::ViscousDampingPrev> damping_prev_assembler = nullptr;
166
168 std::shared_ptr<assembler::Problem> problem;
169
171 std::vector<basis::ElementBases> bases;
173 std::vector<basis::ElementBases> pressure_bases;
175 std::vector<basis::ElementBases> geom_bases_;
176
183
185 std::map<int, Eigen::MatrixXd> polys;
187 std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> polys_3d;
188
190 Eigen::VectorXi disc_orders;
191
193 std::shared_ptr<polyfem::mesh::MeshNodes> mesh_nodes, geom_mesh_nodes, pressure_mesh_nodes;
194
200
204 double avg_mass;
205
207 Eigen::MatrixXd rhs;
208
212
215 std::string formulation() const;
216
219 bool iso_parametric() const;
220
223 const std::vector<basis::ElementBases> &geom_bases() const
224 {
225 return iso_parametric() ? bases : geom_bases_;
226 }
227
232 void build_basis();
236 void assemble_rhs();
240 void assemble_mass_mat();
241
243 std::shared_ptr<assembler::RhsAssembler> build_rhs_assembler(
244 const int n_bases,
245 const std::vector<basis::ElementBases> &bases,
248 std::shared_ptr<assembler::RhsAssembler> build_rhs_assembler() const
249 {
251 }
252
253 std::shared_ptr<assembler::PressureAssembler> build_pressure_assembler(
254 const int n_bases_,
255 const std::vector<basis::ElementBases> &bases_) const;
256 std::shared_ptr<assembler::PressureAssembler> build_pressure_assembler() const
257 {
259 }
260
264 {
266 const int n_b_samples_j = args["space"]["advanced"]["n_boundary_samples"];
267 const int gdiscr_order = mesh->orders().size() <= 0 ? 1 : mesh->orders().maxCoeff();
268 const int discr_order = std::max(disc_orders.maxCoeff(), gdiscr_order);
269
270 const int n_b_samples = std::max(n_b_samples_j, AssemblerUtils::quadrature_order("Mass", discr_order, AssemblerUtils::BasisType::POLY, mesh->dimension()));
271
272 return n_b_samples;
273 }
274
275 private:
279 void sol_to_pressure(Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure);
282
283 public:
286 void set_materials(std::vector<std::shared_ptr<assembler::Assembler>> &assemblers) const;
290
291 //---------------------------------------------------
292 //-----------------solver----------------------------
293 //---------------------------------------------------
294
295 public:
299 void solve_problem(Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure);
303 void solve(Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
304 {
305 if (!mesh)
306 {
307 logger().error("Load the mesh first!");
308 return;
309 }
311
312 build_basis();
313
314 assemble_rhs();
316
317 solve_problem(sol, pressure);
318 }
319
325 void init_solve(Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure);
331 void solve_transient_navier_stokes_split(const int time_steps, const double dt, Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure);
338 void solve_transient_navier_stokes(const int time_steps, const double t0, const double dt, Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure);
345 void solve_transient_linear(const int time_steps, const double t0, const double dt, Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure);
351 void solve_transient_tensor_nonlinear(const int time_steps, const double t0, const double dt, Eigen::MatrixXd &sol);
355 void init_nonlinear_tensor_solve(Eigen::MatrixXd &sol, const double t = 1.0, const bool init_time_integrator = true);
358 void init_linear_solve(Eigen::MatrixXd &sol, const double t = 1.0);
361 void initial_solution(Eigen::MatrixXd &solution) const;
364 void initial_velocity(Eigen::MatrixXd &velocity) const;
367 void initial_acceleration(Eigen::MatrixXd &acceleration) const;
371 void solve_linear(Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure);
375 void solve_navier_stokes(Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure);
379 void solve_tensor_nonlinear(Eigen::MatrixXd &sol, const int t = 0, const bool init_lagging = true);
380
383 std::shared_ptr<polysolve::nonlinear::Solver> make_nl_solver(bool for_al) const;
384
386 std::shared_ptr<utils::PeriodicBoundary> periodic_bc;
387 bool has_periodic_bc() const
388 {
389 return args["boundary_conditions"]["periodic_boundary"]["enabled"].get<bool>();
390 }
391
399 void solve_linear(
400 const std::unique_ptr<polysolve::linear::Solver> &solver,
402 Eigen::VectorXd &b,
403 const bool compute_spectrum,
404 Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure);
405
407 bool is_problem_linear() const { return assembler->is_linear() && !is_contact_enabled() && !is_pressure_enabled() && !has_constraints(); }
408
409 bool has_constraints() const
410 {
411 return has_constraints_;
412 }
413
414 private:
416
417 public:
420 void build_stiffness_mat(StiffnessMatrix &stiffness);
421
422 //---------------------------------------------------
423 //-----------------nodes flags-----------------------
424 //---------------------------------------------------
425
426 public:
428 std::unordered_map<int, std::array<bool, 3>>
429 boundary_conditions_ids(const std::string &bc_type) const;
430
432 std::vector<int> boundary_nodes;
434 std::vector<int> pressure_boundary_nodes;
436 std::vector<mesh::LocalBoundary> total_local_boundary;
438 std::vector<mesh::LocalBoundary> local_boundary;
440 std::vector<mesh::LocalBoundary> local_neumann_boundary;
442 std::vector<mesh::LocalBoundary> local_pressure_boundary;
444 std::unordered_map<int, std::vector<mesh::LocalBoundary>> local_pressure_cavity;
446 std::map<int, basis::InterfaceData> poly_edge_to_data;
448 std::vector<int> dirichlet_nodes;
449 std::vector<RowVectorNd> dirichlet_nodes_position;
451 std::vector<int> neumann_nodes;
452 std::vector<RowVectorNd> neumann_nodes_position;
453
455 Eigen::VectorXi in_node_to_node;
458
459 std::vector<int> primitive_to_node() const;
460 std::vector<int> node_to_primitive() const;
461
462 private:
464 void build_node_mapping();
465
466 //---------------------------------------------------
467 //-----------------Geometry--------------------------
468 //---------------------------------------------------
469 public:
471 std::unique_ptr<mesh::Mesh> mesh;
474
480 void load_mesh(bool non_conforming = false,
481 const std::vector<std::string> &names = std::vector<std::string>(),
482 const std::vector<Eigen::MatrixXi> &cells = std::vector<Eigen::MatrixXi>(),
483 const std::vector<Eigen::MatrixXd> &vertices = std::vector<Eigen::MatrixXd>());
484
490 void load_mesh(GEO::Mesh &meshin, const std::function<int(const size_t, const std::vector<int> &, const RowVectorNd &, bool)> &boundary_marker, bool non_conforming = false, bool skip_boundary_sideset = false);
491
496 void load_mesh(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, bool non_conforming = false)
497 {
498 mesh = mesh::Mesh::create(V, F, non_conforming);
499 load_mesh(non_conforming);
500 }
501
503 void reset_mesh();
504
506 void build_mesh_matrices(Eigen::MatrixXd &V, Eigen::MatrixXi &F);
507
513 bool remesh(const double time, const double dt, Eigen::MatrixXd &sol);
514
515 //---------------------------------------------------
516 //-----------------IPC-------------------------------
517 //---------------------------------------------------
518
520 ipc::CollisionMesh collision_mesh;
521
523 ipc::CollisionMesh periodic_collision_mesh;
526
528 static void build_collision_mesh(
529 const mesh::Mesh &mesh,
530 const int n_bases,
531 const std::vector<basis::ElementBases> &bases,
532 const std::vector<basis::ElementBases> &geom_bases,
533 const std::vector<mesh::LocalBoundary> &total_local_boundary,
535 const json &args,
536 const std::function<std::string(const std::string &)> &resolve_input_path,
537 const Eigen::VectorXi &in_node_to_node,
538 ipc::CollisionMesh &collision_mesh);
539
543
547 bool is_obstacle_vertex(const size_t vi) const
548 {
549 // The obstalce vertices are at the bottom of the collision mesh vertices
550 return vi >= collision_mesh.full_num_vertices() - obstacle.n_vertices();
551 }
552
557 {
558 return args["contact"]["enabled"];
559 }
560
565 {
566 return args["contact"]["adhesion"]["adhesion_enabled"];
567 }
568
573 {
574 return (args["boundary_conditions"]["pressure_boundary"].size() > 0)
575 || (args["boundary_conditions"]["pressure_cavity"].size() > 0);
576 }
577
579 bool has_dhat = false;
580
581 //---------------------------------------------------
582 //-----------------OUTPUT----------------------------
583 //---------------------------------------------------
584 public:
586 std::string output_dir;
596
600 void export_data(const Eigen::MatrixXd &sol, const Eigen::MatrixXd &pressure);
601
609 void save_timestep(const double time, const int t, const double t0, const double dt, const Eigen::MatrixXd &sol, const Eigen::MatrixXd &pressure);
610
616 void save_subsolve(const int i, const int t, const Eigen::MatrixXd &sol, const Eigen::MatrixXd &pressure);
617
621 void save_json(const Eigen::MatrixXd &sol, std::ostream &out);
622
625 void save_json(const Eigen::MatrixXd &sol);
626
628 void compute_errors(const Eigen::MatrixXd &sol);
629
632 void save_restart_json(const double t0, const double dt, const int t) const;
633
634 //-----------PATH management
637 std::string root_path() const;
638
643 std::string resolve_input_path(const std::string &path, const bool only_if_exists = false) const;
644
648 std::string resolve_output_path(const std::string &path) const;
649
650 //---------------------------------------------------
651 //-----------------differentiable--------------------
652 //---------------------------------------------------
653 public:
655 void cache_transient_adjoint_quantities(const int current_step, const Eigen::MatrixXd &sol, const Eigen::MatrixXd &disp_grad);
657
658 std::unique_ptr<polysolve::linear::Solver> lin_solver_cached; // matrix factorization of last linear solve
659
660 int ndof() const
661 {
662 const int actual_dim = problem->is_scalar() ? 1 : mesh->dimension();
663 if (mixed_assembler == nullptr)
664 return actual_dim * n_bases;
665 else
666 return actual_dim * n_bases + n_pressure_bases;
667 }
668
669 // Aux functions for setting up adjoint equations
670 void compute_force_jacobian(const Eigen::MatrixXd &sol, const Eigen::MatrixXd &disp_grad, StiffnessMatrix &hessian);
671 void compute_force_jacobian_prev(const int force_step, const int sol_step, StiffnessMatrix &hessian_prev) const;
672 // Solves the adjoint PDE for derivatives and caches
673 void solve_adjoint_cached(const Eigen::MatrixXd &rhs);
674 Eigen::MatrixXd solve_adjoint(const Eigen::MatrixXd &rhs) const;
675 // Returns cached adjoint solve
676 Eigen::MatrixXd get_adjoint_mat(int type) const
677 {
678 assert(diff_cached.adjoint_mat().size() > 0);
679 if (problem->is_time_dependent())
680 {
681 if (type == 0)
682 return diff_cached.adjoint_mat().leftCols(diff_cached.adjoint_mat().cols() / 2);
683 else if (type == 1)
684 return diff_cached.adjoint_mat().middleCols(diff_cached.adjoint_mat().cols() / 2, diff_cached.adjoint_mat().cols() / 2);
685 else
686 log_and_throw_adjoint_error("Invalid adjoint type!");
687 }
688
689 return diff_cached.adjoint_mat();
690 }
691 Eigen::MatrixXd solve_static_adjoint(const Eigen::MatrixXd &adjoint_rhs) const;
692 Eigen::MatrixXd solve_transient_adjoint(const Eigen::MatrixXd &adjoint_rhs) const;
693 // Change geometric node positions
694 void set_mesh_vertex(int v_id, const Eigen::VectorXd &vertex);
695 void get_vertices(Eigen::MatrixXd &vertices) const;
696 void get_elements(Eigen::MatrixXi &elements) const;
697
698 // Get geometric node indices for surface/volume
699 void compute_surface_node_ids(const int surface_selection, std::vector<int> &node_ids) const;
700 void compute_total_surface_node_ids(std::vector<int> &node_ids) const;
701 void compute_volume_node_ids(const int volume_selection, std::vector<int> &node_ids) const;
702
703 // to replace the initial condition in json during initial condition optimization
705 // mapping from positions of FE basis nodes to positions of geometry nodes
707
708 //---------------------------------------------------
709 //-----------------homogenization--------------------
710 //---------------------------------------------------
711 public:
713
715 void solve_homogenization_step(Eigen::MatrixXd &sol, const int t = 0, bool adaptive_initial_weight = false); // sol is the extended solution, i.e. [periodic fluctuation, macro strain]
716 void init_homogenization_solve(const double t);
717 void solve_homogenization(const int time_steps, const double t0, const double dt, Eigen::MatrixXd &sol);
718 bool is_homogenization() const
719 {
720 return args["boundary_conditions"]["periodic_boundary"]["linear_displacement_offset"].size() > 0;
721 }
722 };
723
724} // namespace polyfem
int V
main class that contains the polyfem solver and all its state
Definition State.hpp:79
void get_vertices(Eigen::MatrixXd &vertices) const
Definition StateDiff.cpp:71
Eigen::MatrixXd initial_vel_update
Definition State.hpp:704
io::OutRuntimeData timings
runtime statistics
Definition State.hpp:590
assembler::MacroStrainValue macro_strain_constraint
Definition State.hpp:712
std::shared_ptr< utils::PeriodicBoundary > periodic_bc
periodic BC and periodic mesh utils
Definition State.hpp:386
int n_bases
number of bases
Definition State.hpp:178
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
Definition State.hpp:180
std::string formulation() const
return the formulation (checks if the problem is scalar or not and deals with multiphysics)
Definition State.cpp:328
assembler::AssemblyValsCache ass_vals_cache
used to store assembly values for small problems
Definition State.hpp:196
int n_boundary_samples() const
quadrature used for projecting boundary conditions
Definition State.hpp:263
void set_mesh_vertex(int v_id, const Eigen::VectorXd &vertex)
Definition StateDiff.cpp:96
bool is_adhesion_enabled() const
does the simulation have adhesion
Definition State.hpp:564
const std::vector< basis::ElementBases > & geom_bases() const
Get a constant reference to the geometry mapping bases.
Definition State.hpp:223
double starting_max_edge_length
Definition State.hpp:594
std::vector< mesh::LocalBoundary > local_pressure_boundary
mapping from elements to nodes for pressure boundary conditions
Definition State.hpp:442
void sol_to_pressure(Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
splits the solution in solution and pressure for mixed problems
Definition State.cpp:371
std::shared_ptr< assembler::ViscousDamping > damping_assembler
Definition State.hpp:164
void init(const json &args, const bool strict_validation)
initialize the polyfem solver with a json settings
std::string root_path() const
Get the root path for the state (e.g., args["root_path"] or ".")
bool is_pressure_enabled() const
does the simulation has pressure
Definition State.hpp:572
void assemble_mass_mat()
assemble mass, step 4 of solve build mass matrix based on defined basis modifies mass (and maybe more...
Definition State.cpp:1480
Eigen::VectorXi in_node_to_node
Inpute nodes (including high-order) to polyfem nodes, only for isoparametric.
Definition State.hpp:455
void build_polygonal_basis()
builds bases for polygons, called inside build_basis
Definition State.cpp:1061
void solve(Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
solves the problem, call other methods
Definition State.hpp:303
mesh::Obstacle obstacle
Obstacles used in collisions.
Definition State.hpp:473
std::shared_ptr< assembler::Assembler > assembler
assemblers
Definition State.hpp:155
void compute_errors(const Eigen::MatrixXd &sol)
computes all errors
ipc::CollisionMesh periodic_collision_mesh
IPC collision mesh under periodic BC.
Definition State.hpp:523
void build_periodic_collision_mesh()
Definition State.cpp:1224
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
io::OutGeometryData out_geom
visualization stuff
Definition State.hpp:588
bool use_avg_pressure
use average pressure for stokes problem to fix the additional dofs, true by default if false,...
Definition State.hpp:211
std::string resolve_output_path(const std::string &path) const
Resolve output path relative to output_dir if the path is not absolute.
std::vector< int > neumann_nodes
per node neumann
Definition State.hpp:451
ipc::CollisionMesh collision_mesh
IPC collision mesh.
Definition State.hpp:520
std::string output_dir
Directory for output files.
Definition State.hpp:586
void load_mesh(bool non_conforming=false, const std::vector< std::string > &names=std::vector< std::string >(), const std::vector< Eigen::MatrixXi > &cells=std::vector< Eigen::MatrixXi >(), const std::vector< Eigen::MatrixXd > &vertices=std::vector< Eigen::MatrixXd >())
loads the mesh from the json arguments
Definition StateLoad.cpp:91
std::vector< basis::ElementBases > geom_bases_
Geometric mapping bases, if the elements are isoparametric, this list is empty.
Definition State.hpp:175
solver::CacheLevel optimization_enabled
Definition State.hpp:654
void solve_homogenization_step(Eigen::MatrixXd &sol, const int t=0, bool adaptive_initial_weight=false)
In Elasticity PDE, solve for "min W(disp_grad + \grad u)" instead of "min W(\grad u)".
Eigen::VectorXi in_primitive_to_primitive
maps in vertices/edges/faces/cells to polyfem vertices/edges/faces/cells
Definition State.hpp:457
void set_materials(std::vector< std::shared_ptr< assembler::Assembler > > &assemblers) const
set the material and the problem dimension
std::string get_log(const Eigen::MatrixXd &sol)
gets the output log as json this is not what gets printed but more informative information,...
Definition State.hpp:130
std::map< int, basis::InterfaceData > poly_edge_to_data
nodes on the boundary of polygonal elements, used for harmonic bases
Definition State.hpp:446
std::vector< int > pressure_boundary_nodes
list of neumann boundary nodes
Definition State.hpp:434
void set_max_threads(const int max_threads=std::numeric_limits< int >::max())
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
std::vector< basis::ElementBases > pressure_bases
FE pressure bases for mixed elements, the size is #elements.
Definition State.hpp:173
void save_json(const Eigen::MatrixXd &sol, std::ostream &out)
saves the output statistic to a stream
bool is_homogenization() const
Definition State.hpp:718
Eigen::MatrixXd get_adjoint_mat(int type) const
Definition State.hpp:676
void load_mesh(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, bool non_conforming=false)
loads the mesh from V and F,
Definition State.hpp:496
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.
bool has_constraints() const
Definition State.hpp:409
std::shared_ptr< assembler::Assembler > pressure_assembler
Definition State.hpp:160
std::shared_ptr< polysolve::nonlinear::Solver > make_nl_solver(bool for_al) const
factory to create the nl solver depending on input
double starting_min_edge_length
Definition State.hpp:593
std::shared_ptr< assembler::Mass > mass_matrix_assembler
Definition State.hpp:157
void build_node_mapping()
build the mapping from input nodes to polyfem nodes
Definition State.cpp:245
bool remesh(const double time, const double dt, Eigen::MatrixXd &sol)
Remesh the FE space and update solution(s).
std::vector< int > primitive_to_node() const
Definition State.cpp:228
bool has_dhat
stores if input json contains dhat
Definition State.hpp:579
Eigen::MatrixXd solve_static_adjoint(const Eigen::MatrixXd &adjoint_rhs) const
void compute_force_jacobian(const Eigen::MatrixXd &sol, const Eigen::MatrixXd &disp_grad, StiffnessMatrix &hessian)
std::shared_ptr< polyfem::mesh::MeshNodes > pressure_mesh_nodes
Definition State.hpp:193
std::unique_ptr< mesh::Mesh > mesh
current mesh, it can be a Mesh2D or Mesh3D
Definition State.hpp:471
void reset_mesh()
Resets the mesh.
Definition StateLoad.cpp:22
std::unordered_map< int, std::array< bool, 3 > > boundary_conditions_ids(const std::string &bc_type) const
Construct a vector of boundary conditions ids with their dimension flags.
std::shared_ptr< polyfem::mesh::MeshNodes > mesh_nodes
Mapping from input nodes to FE nodes.
Definition State.hpp:193
StiffnessMatrix mass
Mass matrix, it is computed only for time dependent problems.
Definition State.hpp:202
void solve_homogenization(const int time_steps, const double t0, const double dt, Eigen::MatrixXd &sol)
std::vector< int > dirichlet_nodes
per node dirichlet
Definition State.hpp:448
bool has_periodic_bc() const
Definition State.hpp:387
void solve_problem(Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
solves the problems
Definition State.cpp:1676
void build_basis()
builds the bases step 2 of solve modifies bases, pressure_bases, geom_bases_, boundary_nodes,...
Definition State.cpp:506
void solve_navier_stokes(Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
solves a navier stokes
std::shared_ptr< assembler::Problem > problem
current problem, it contains rhs and bc
Definition State.hpp:168
bool has_constraints_
Definition State.hpp:415
std::vector< RowVectorNd > dirichlet_nodes_position
Definition State.hpp:449
spdlog::sink_ptr file_sink_
Definition State.hpp:143
spdlog::sink_ptr console_sink_
logger sink to stdout
Definition State.hpp:142
State()
Constructor.
Definition StateInit.cpp:56
std::map< int, std::pair< Eigen::MatrixXd, Eigen::MatrixXi > > polys_3d
polyhedra, used since poly have no geom mapping
Definition State.hpp:187
std::vector< int > node_to_primitive() const
Definition State.cpp:235
~State()=default
json args
main input arguments containing all defaults
Definition State.hpp:101
void set_log_level(const spdlog::level::level_enum log_level)
change log level
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.
solver::DiffCache diff_cached
Definition State.hpp:656
void initial_velocity(Eigen::MatrixXd &velocity) const
Load or compute the initial velocity.
io::OutStatsData stats
Other statistics.
Definition State.hpp:592
std::vector< RowVectorNd > neumann_nodes_position
Definition State.hpp:452
void assemble_rhs()
compute rhs, step 3 of solve build rhs vector based on defined basis and given rhs of the problem mod...
Definition State.cpp:1604
double min_boundary_edge_length
Definition State.hpp:595
Eigen::VectorXi periodic_collision_mesh_to_basis
index mapping from periodic 2x2 collision mesh to FE periodic mesh
Definition State.hpp:525
void init_time()
initialize time settings if args contains "time"
void initial_acceleration(Eigen::MatrixXd &acceleration) const
Load or compute the initial acceleration.
std::vector< basis::ElementBases > bases
FE bases, the size is #elements.
Definition State.hpp:171
void solve_adjoint_cached(const Eigen::MatrixXd &rhs)
void solve_transient_navier_stokes_split(const int time_steps, const double dt, Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
solves transient navier stokes with operator splitting
Eigen::MatrixXd solve_adjoint(const Eigen::MatrixXd &rhs) const
void get_elements(Eigen::MatrixXi &elements) const
Definition StateDiff.cpp:79
int n_geom_bases
number of geometric bases
Definition State.hpp:182
assembler::AssemblyValsCache pressure_ass_vals_cache
used to store assembly values for pressure for small problems
Definition State.hpp:199
void build_stiffness_mat(StiffnessMatrix &stiffness)
utility that builds the stiffness matrix and collects stats, used only for linear problems
std::vector< mesh::LocalBoundary > local_boundary
mapping from elements to nodes for dirichlet boundary conditions
Definition State.hpp:438
double avg_mass
average system mass, used for contact with IPC
Definition State.hpp:204
bool iso_parametric() const
check if using iso parametric bases
Definition State.cpp:462
std::shared_ptr< assembler::RhsAssembler > build_rhs_assembler() const
build a RhsAssembler for the problem
Definition State.hpp:248
void compute_volume_node_ids(const int volume_selection, std::vector< int > &node_ids) const
std::vector< mesh::LocalBoundary > total_local_boundary
mapping from elements to nodes for all mesh
Definition State.hpp:436
void solve_transient_navier_stokes(const int time_steps, const double t0, const double dt, Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
solves transient navier stokes with FEM
void init_solve(Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
initialize solver
void init_linear_solve(Eigen::MatrixXd &sol, const double t=1.0)
initialize the linear solve
int ndof() const
Definition State.hpp:660
void export_data(const Eigen::MatrixXd &sol, const Eigen::MatrixXd &pressure)
saves all data on the disk according to the input params
assembler::AssemblyValsCache mass_ass_vals_cache
Definition State.hpp:197
bool is_obstacle_vertex(const size_t vi) const
checks if vertex is obstacle
Definition State.hpp:547
void solve_tensor_nonlinear(Eigen::MatrixXd &sol, const int t=0, const bool init_lagging=true)
solves nonlinear problems
std::string resolve_input_path(const std::string &path, const bool only_if_exists=false) const
Resolve input path relative to root_path() if the path is not absolute.
std::vector< int > boundary_nodes
list of boundary nodes
Definition State.hpp:432
solver::SolveData solve_data
timedependent stuff cached
Definition State.hpp:321
std::unique_ptr< polysolve::linear::Solver > lin_solver_cached
Definition State.hpp:658
std::shared_ptr< polyfem::mesh::MeshNodes > geom_mesh_nodes
Definition State.hpp:193
void solve_linear(Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
solves a linear problem
void init_nonlinear_tensor_solve(Eigen::MatrixXd &sol, const double t=1.0, const bool init_time_integrator=true)
initialize the nonlinear solver
void init_logger(const std::string &log_file, const spdlog::level::level_enum log_level, const spdlog::level::level_enum file_log_level, const bool is_quiet)
initializing the logger
Definition StateInit.cpp:65
std::shared_ptr< assembler::ViscousDampingPrev > damping_prev_assembler
Definition State.hpp:165
bool is_contact_enabled() const
does the simulation have contact
Definition State.hpp:556
Eigen::MatrixXd solve_transient_adjoint(const Eigen::MatrixXd &adjoint_rhs) const
void build_collision_mesh()
extracts the boundary mesh for collision, called in build_basis
Definition State.cpp:1357
Eigen::VectorXi disc_orders
vector of discretization orders, used when not all elements have the same degree, one per element
Definition State.hpp:190
void compute_force_jacobian_prev(const int force_step, const int sol_step, StiffnessMatrix &hessian_prev) const
std::vector< mesh::LocalBoundary > local_neumann_boundary
mapping from elements to nodes for neumann boundary conditions
Definition State.hpp:440
std::shared_ptr< assembler::PressureAssembler > build_pressure_assembler() const
Definition State.hpp:256
void solve_transient_linear(const int time_steps, const double t0, const double dt, Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
solves transient linear problem
void compute_total_surface_node_ids(std::vector< int > &node_ids) const
std::shared_ptr< assembler::PressureAssembler > elasticity_pressure_assembler
Definition State.hpp:162
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.
void compute_surface_node_ids(const int surface_selection, std::vector< int > &node_ids) const
bool is_problem_linear() const
Returns whether the system is linear. Collisions and pressure add nonlinearity to the problem.
Definition State.hpp:407
std::shared_ptr< assembler::MixedAssembler > mixed_assembler
Definition State.hpp:159
Eigen::MatrixXd initial_sol_update
Definition State.hpp:704
std::map< int, Eigen::MatrixXd > polys
polygons, used since poly have no geom mapping
Definition State.hpp:185
Eigen::MatrixXd rhs
System right-hand side.
Definition State.hpp:207
std::unordered_map< int, std::vector< mesh::LocalBoundary > > local_pressure_cavity
mapping from elements to nodes for pressure boundary conditions
Definition State.hpp:444
void init_homogenization_solve(const double t)
StiffnessMatrix basis_nodes_to_gbasis_nodes
Definition State.hpp:706
static int quadrature_order(const std::string &assembler, const int basis_degree, const BasisType &b_type, const int dim)
utility for retrieving the needed quadrature order to precisely integrate the given form on the given...
Caches basis evaluation and geometric mapping at every element.
Utilies related to export of geometry.
Definition OutData.hpp:33
stores all runtime data
Definition OutData.hpp:341
all stats from polyfem
Definition OutData.hpp:368
void compute_mesh_stats(const polyfem::mesh::Mesh &mesh)
compute stats (counts els type, mesh lenght, etc), step 1 of solve
Definition OutData.cpp:2743
Abstract mesh class to capture 2d/3d conforming and non-conforming meshes.
Definition Mesh.hpp:39
static std::unique_ptr< Mesh > create(const std::string &path, const bool non_conforming=false)
factory to build the proper mesh
Definition Mesh.cpp:173
const Eigen::MatrixXd & adjoint_mat() const
class to store time stepping data
Definition SolveData.hpp:59
Used for test only.
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:44
nlohmann::json json
Definition Common.hpp:9
void log_and_throw_adjoint_error(const std::string &msg)
Definition Logger.cpp:79
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > RowVectorNd
Definition Types.hpp:13
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix
Definition Types.hpp:22