Loading [MathJax]/extensions/tex2jax.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
StateHomogenization.cpp
Go to the documentation of this file.
1#include <polyfem/State.hpp>
7#include <polysolve/linear/FEMSolver.hpp>
8#include <polysolve/nonlinear/Solver.hpp>
9
13
14#include <unsupported/Eigen/SparseExtra>
15
18
19#include <ipc/ipc.hpp>
20
21namespace polyfem
22{
23
24 using namespace assembler;
25 using namespace mesh;
26 using namespace solver;
27 using namespace utils;
28 using namespace quadrature;
29
31 {
32 const int dim = mesh->dimension();
33 const int ndof = n_bases * dim;
34
35 const std::vector<std::shared_ptr<Form>> forms = solve_data.init_forms(
36 // General
37 units,
38 mesh->dimension(), t, in_node_to_node,
39 // Elastic form
41 args["solver"]["advanced"]["jacobian_threshold"], args["solver"]["advanced"]["check_inversion"],
42 // Body form
44 n_boundary_samples(), rhs, Eigen::VectorXd::Zero(ndof) /* only to set neumann BC, not used*/, mass_matrix_assembler->density(),
45 // Pressure form
47 // Inertia form
48 args.value("/time/quasistatic"_json_pointer, true), mass,
49 nullptr,
50 // Lagged regularization form
51 args["solver"]["advanced"]["lagged_regularization_weight"],
52 args["solver"]["advanced"]["lagged_regularization_iterations"],
53 // Augmented lagrangian form
54 obstacle.ndof(), args["constraints"]["hard"], args["constraints"]["soft"],
55 // Contact form
56 args["contact"]["enabled"], args["contact"]["periodic"].get<bool>() ? periodic_collision_mesh : collision_mesh, args["contact"]["dhat"],
57 avg_mass, args["contact"]["use_convergent_formulation"] ? bool(args["contact"]["use_area_weighting"]) : false,
58 args["contact"]["use_convergent_formulation"] ? bool(args["contact"]["use_improved_max_operator"]) : false,
59 args["contact"]["use_convergent_formulation"] ? bool(args["contact"]["use_physical_barrier"]) : false,
60 args["solver"]["contact"]["barrier_stiffness"],
61 args["solver"]["contact"]["CCD"]["broad_phase"],
62 args["solver"]["contact"]["CCD"]["tolerance"],
63 args["solver"]["contact"]["CCD"]["max_iterations"],
65 // Normal Adhesion Form
66 args["contact"]["adhesion"]["adhesion_enabled"],
67 args["contact"]["adhesion"]["dhat_p"],
68 args["contact"]["adhesion"]["dhat_a"],
69 args["contact"]["adhesion"]["adhesion_strength"],
70 // Tangential Adhesion Form
71 args["contact"]["adhesion"]["tangential_adhesion_coefficient"],
72 args["contact"]["adhesion"]["epsa"],
73 args["solver"]["contact"]["tangential_adhesion_iterations"],
74 // Homogenization
76 // Periodic contact
77 args["contact"]["periodic"], periodic_collision_mesh_to_basis, periodic_bc,
78 // Friction form
79 args["contact"]["friction_coefficient"],
80 args["contact"]["epsv"],
81 args["solver"]["contact"]["friction_iterations"],
82 // Rayleigh damping form
83 args["solver"]["rayleigh_damping"]);
84
85 for (const auto &[name, form] : solve_data.named_forms())
86 {
87 if (name == "augmented_lagrangian")
88 {
89 form->set_weight(0);
90 form->disable();
91 }
92 }
93
94 bool solve_symmetric_flag = false;
95 {
96 const auto &fixed_entry = macro_strain_constraint.get_fixed_entry();
97 for (int i = 0; i < dim; i++)
98 {
99 for (int j = 0; j < i; j++)
100 {
101 if (std::find(fixed_entry.data(), fixed_entry.data() + fixed_entry.size(), i + j * dim) == fixed_entry.data() + fixed_entry.size() && std::find(fixed_entry.data(), fixed_entry.data() + fixed_entry.size(), j + i * dim) == fixed_entry.data() + fixed_entry.size())
102 {
103 logger().info("Strain entry [{},{}] and [{},{}] are not fixed, solve for symmetric strain...", i, j, j, i);
104 solve_symmetric_flag = true;
105 break;
106 }
107 }
108 if (solve_symmetric_flag)
109 break;
110 }
111 }
112
113 std::shared_ptr<NLHomoProblem> homo_problem = std::make_shared<NLHomoProblem>(
114 ndof,
116 *this, t, forms, solve_data.al_form, solve_symmetric_flag, polysolve::linear::Solver::create(args["solver"]["linear"], logger()));
118 homo_problem->add_form(solve_data.periodic_contact_form);
120 homo_problem->add_form(solve_data.strain_al_lagr_form);
121
122 solve_data.nl_problem = homo_problem;
123 solve_data.nl_problem->init(Eigen::VectorXd::Zero(homo_problem->reduced_size() + homo_problem->macro_reduced_size()));
124 solve_data.nl_problem->update_quantities(t, Eigen::VectorXd::Zero(homo_problem->reduced_size() + homo_problem->macro_reduced_size()));
125 }
126
127 void State::solve_homogenization_step(Eigen::MatrixXd &sol, const int t, bool adaptive_initial_weight)
128 {
129 const int dim = mesh->dimension();
130 const int ndof = n_bases * dim;
131
132 auto homo_problem = std::dynamic_pointer_cast<NLHomoProblem>(solve_data.nl_problem);
133
134 Eigen::VectorXd extended_sol;
135 extended_sol.setZero(ndof + dim * dim);
136
137 if (sol.size() == extended_sol.size())
138 extended_sol = sol;
139
140 const auto &fixed_entry = macro_strain_constraint.get_fixed_entry();
141 homo_problem->set_fixed_entry({});
142 {
143 std::shared_ptr<polysolve::nonlinear::Solver> nl_solver = make_nl_solver(true);
144
145 Eigen::VectorXi al_indices = fixed_entry.array() + homo_problem->full_size();
146 Eigen::VectorXd al_values = utils::flatten(macro_strain_constraint.eval(t))(fixed_entry);
147
148 std::shared_ptr<MacroStrainLagrangianForm> lagr_form = solve_data.strain_al_lagr_form;
149 lagr_form->enable();
150
151 const double initial_weight = args["solver"]["augmented_lagrangian"]["initial_weight"];
152 const double max_weight = args["solver"]["augmented_lagrangian"]["max_weight"];
153 const double eta_tol = args["solver"]["augmented_lagrangian"]["eta"];
154 const double scaling = args["solver"]["augmented_lagrangian"]["scaling"];
155 double al_weight = initial_weight;
156
157 Eigen::VectorXd tmp_sol = homo_problem->extended_to_reduced(extended_sol);
158 const Eigen::VectorXd initial_sol = tmp_sol;
159 const double initial_error = lagr_form->compute_error(extended_sol);
160 double current_error = initial_error;
161
162 // try to enforce fixed values on macro strain
163 extended_sol(al_indices) = al_values;
164 Eigen::VectorXd reduced_sol = homo_problem->extended_to_reduced(extended_sol);
165
166 homo_problem->line_search_begin(tmp_sol, reduced_sol);
167 int al_steps = 0;
168 bool force_al = true;
169
170 lagr_form->set_initial_weight(al_weight);
171
172 while (force_al
173 || !std::isfinite(homo_problem->value(reduced_sol))
174 || !homo_problem->is_step_valid(tmp_sol, reduced_sol)
175 || !homo_problem->is_step_collision_free(tmp_sol, reduced_sol))
176 {
177 force_al = false;
178 homo_problem->line_search_end();
179
180 logger().info("Solving AL Problem with weight {}", al_weight);
181
182 homo_problem->init(tmp_sol);
183 try
184 {
185 homo_problem->normalize_forms();
186 nl_solver->minimize(*homo_problem, tmp_sol);
187 }
188 catch (const std::runtime_error &e)
189 {
190 logger().error("AL solve failed!");
191 }
192
193 extended_sol = homo_problem->reduced_to_extended(tmp_sol);
194 logger().debug("Current macro strain: {}", extended_sol.tail(dim * dim));
195
196 current_error = lagr_form->compute_error(extended_sol);
197 const double eta = 1 - sqrt(current_error / initial_error);
198
199 logger().info("Current eta = {}, current error = {}, initial error = {}", eta, current_error, initial_error);
200
201 if (eta < eta_tol && al_weight < max_weight)
202 al_weight *= scaling;
203 else
204 lagr_form->update_lagrangian(extended_sol, al_weight);
205
206 if (eta <= 0)
207 {
208 if (adaptive_initial_weight)
209 {
210 args["solver"]["augmented_lagrangian"]["initial_weight"] = args["solver"]["augmented_lagrangian"]["initial_weight"].get<double>() * scaling;
211 {
212 json tmp = json::object();
213 tmp["/solver/augmented_lagrangian/initial_weight"_json_pointer] = args["solver"]["augmented_lagrangian"]["initial_weight"];
214 }
215 logger().warn("AL weight too small, increase weight and revert solution, new initial weight is {}", args["solver"]["augmented_lagrangian"]["initial_weight"].get<double>());
216 }
217 tmp_sol = initial_sol;
218 }
219
220 // try to enforce fixed values on macro strain
221 extended_sol(al_indices) = al_values;
222 reduced_sol = homo_problem->extended_to_reduced(extended_sol);
223
224 homo_problem->line_search_begin(tmp_sol, reduced_sol);
225 }
226 homo_problem->line_search_end();
227 lagr_form->disable();
228 }
229
230 homo_problem->set_fixed_entry(fixed_entry);
231
232 Eigen::VectorXd reduced_sol = homo_problem->extended_to_reduced(extended_sol);
233
234 homo_problem->init(reduced_sol);
235 std::shared_ptr<polysolve::nonlinear::Solver> nl_solver = make_nl_solver(false);
236 homo_problem->normalize_forms();
237 nl_solver->minimize(*homo_problem, reduced_sol);
238
239 logger().info("Macro Strain: {}", extended_sol.tail(dim * dim).transpose());
240
241 // check saddle point
242 {
243 json linear_args = args["solver"]["linear"];
244 std::string solver_name = linear_args["solver"];
245 if (solver_name.find("Pardiso") != std::string::npos)
246 {
247 linear_args["solver"] = "Eigen::PardisoLLT";
248 std::unique_ptr<polysolve::linear::Solver> solver =
249 polysolve::linear::Solver::create(linear_args, logger());
250
252 homo_problem->hessian(reduced_sol, A);
253 Eigen::VectorXd x, b = Eigen::VectorXd::Zero(A.rows());
254 try
255 {
256 dirichlet_solve(
257 *solver, A, b, {}, x, A.rows(), args["output"]["data"]["stiffness_mat"], false, false, false);
258 }
259 catch (const std::runtime_error &error)
260 {
261 logger().error("The solution is a saddle point!");
262 }
263 }
264 }
265
266 sol = homo_problem->reduced_to_extended(reduced_sol);
267 if (args["/boundary_conditions/periodic_boundary/force_zero_mean"_json_pointer].get<bool>())
268 {
269 Eigen::VectorXd integral = io::Evaluator::integrate_function(bases, geom_bases(), sol, dim, dim);
270 double area = io::Evaluator::integrate_function(bases, geom_bases(), Eigen::VectorXd::Ones(n_bases), dim, 1)(0);
271 for (int d = 0; d < dim; d++)
272 sol(Eigen::seqN(d, n_bases, dim), 0).array() -= integral(d) / area;
273
274 reduced_sol = homo_problem->extended_to_reduced(sol);
275 }
276
278 cache_transient_adjoint_quantities(t, homo_problem->reduced_to_full(reduced_sol), utils::unflatten(sol.bottomRows(dim * dim), dim));
279 }
280
281 void State::solve_homogenization(const int time_steps, const double t0, const double dt, Eigen::MatrixXd &sol)
282 {
283 bool is_static = !is_param_valid(args, "time");
284 if (!is_static && !args["time"]["quasistatic"])
285 log_and_throw_error("Transient homogenization can only do quasi-static!");
286
288
289 const int dim = mesh->dimension();
290 Eigen::MatrixXd extended_sol;
291 for (int t = 0; t <= time_steps; ++t)
292 {
293 double forward_solve_time = 0, remeshing_time = 0, global_relaxation_time = 0;
294
295 {
296 POLYFEM_SCOPED_TIMER(forward_solve_time);
297 solve_homogenization_step(extended_sol, t, false);
298 }
299 sol = extended_sol.topRows(extended_sol.size() - dim * dim) + io::Evaluator::generate_linear_field(n_bases, mesh_nodes, utils::unflatten(extended_sol.bottomRows(dim * dim), dim));
300
301 if (is_static)
302 return;
303
304 // Always save the solution for consistency
305 save_timestep(t0 + dt * t, t, t0, dt, sol, Eigen::MatrixXd()); // no pressure
306
307 {
308 POLYFEM_SCOPED_TIMER("Update quantities");
309
310 // solve_data.time_integrator->update_quantities(sol);
311
312 solve_data.nl_problem->update_quantities(t0 + (t + 1) * dt, sol);
313
316 }
317
318 logger().info("{}/{} t={}", t, time_steps, t0 + dt * t);
319
320 // const std::string rest_mesh_path = args["output"]["data"]["rest_mesh"].get<std::string>();
321 // if (!rest_mesh_path.empty())
322 // {
323 // Eigen::MatrixXd V;
324 // Eigen::MatrixXi F;
325 // build_mesh_matrices(V, F);
326 // io::MshWriter::write(
327 // resolve_output_path(fmt::format(args["output"]["data"]["rest_mesh"], t)),
328 // V, F, mesh->get_body_ids(), mesh->is_volume(), /*binary=*/true);
329 // }
330
331 // const std::string &state_path = resolve_output_path(fmt::format(args["output"]["data"]["state"], t));
332 // if (!state_path.empty())
333 // solve_data.time_integrator->save_state(state_path);
334
335 // save restart file
336 save_restart_json(t0, dt, t);
337 // stats_csv.write(t, forward_solve_time, remeshing_time, global_relaxation_time, sol);
338 }
339 }
340
341} // namespace polyfem
Quadrature quadrature
int x
#define POLYFEM_SCOPED_TIMER(...)
Definition Timer.hpp:10
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
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
const std::vector< basis::ElementBases > & geom_bases() const
Get a constant reference to the geometry mapping bases.
Definition State.hpp:223
std::vector< mesh::LocalBoundary > local_pressure_boundary
mapping from elements to nodes for pressure boundary conditions
Definition State.hpp:442
Eigen::VectorXi in_node_to_node
Inpute nodes (including high-order) to polyfem nodes, only for isoparametric.
Definition State.hpp:455
mesh::Obstacle obstacle
Obstacles used in collisions.
Definition State.hpp:473
std::shared_ptr< assembler::Assembler > assembler
assemblers
Definition State.hpp:155
ipc::CollisionMesh periodic_collision_mesh
IPC collision mesh under periodic BC.
Definition State.hpp:523
ipc::CollisionMesh collision_mesh
IPC collision mesh.
Definition State.hpp:520
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)".
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::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
Definition State.hpp:157
std::unique_ptr< mesh::Mesh > mesh
current mesh, it can be a Mesh2D or Mesh3D
Definition State.hpp:471
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)
json args
main input arguments containing all defaults
Definition State.hpp:101
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.
Eigen::VectorXi periodic_collision_mesh_to_basis
index mapping from periodic 2x2 collision mesh to FE periodic mesh
Definition State.hpp:525
std::vector< basis::ElementBases > bases
FE bases, the size is #elements.
Definition State.hpp:171
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
int ndof() const
Definition State.hpp:660
assembler::AssemblyValsCache mass_ass_vals_cache
Definition State.hpp:197
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::vector< mesh::LocalBoundary > local_neumann_boundary
mapping from elements to nodes for neumann boundary conditions
Definition State.hpp:440
std::shared_ptr< assembler::PressureAssembler > elasticity_pressure_assembler
Definition State.hpp:162
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)
Eigen::MatrixXd eval(const double t) const
const Eigen::VectorXi & get_fixed_entry() const
static Eigen::MatrixXd generate_linear_field(const int n_bases, const std::shared_ptr< mesh::MeshNodes > mesh_nodes, const Eigen::MatrixXd &grad)
static Eigen::VectorXd integrate_function(const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const Eigen::MatrixXd &fun, const int dim, const int actual_dim)
std::shared_ptr< solver::PeriodicContactForm > periodic_contact_form
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 &macro_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.
Definition SolveData.cpp:34
void update_dt()
updates the dt inside the different forms
std::shared_ptr< solver::NLProblem > nl_problem
std::shared_ptr< solver::MacroStrainLagrangianForm > strain_al_lagr_form
std::vector< std::pair< std::string, std::shared_ptr< solver::Form > > > named_forms() const
std::vector< std::shared_ptr< solver::AugmentedLagrangianForm > > al_form
void update_barrier_stiffness(const Eigen::VectorXd &x)
update the barrier stiffness for the forms
bool is_param_valid(const json &params, const std::string &key)
Determine if a key exists and is non-null in a json object.
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 & logger()
Retrieves the current logger.
Definition Logger.cpp:44
nlohmann::json json
Definition Common.hpp:9
void log_and_throw_error(const std::string &msg)
Definition Logger.cpp:73
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix
Definition Types.hpp:22