Loading [MathJax]/extensions/tex2jax.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
StateInit.cpp
Go to the documentation of this file.
1#include <polyfem/State.hpp>
2
6
9
14
16
17#include <jse/jse.h>
18
19#include <spdlog/sinks/stdout_color_sinks.h>
20#include <spdlog/sinks/basic_file_sink.h>
21#include <spdlog/sinks/ostream_sink.h>
22
23#include <ipc/utils/logger.hpp>
24#include <wmtk/utils/Logger.hpp>
25
28
29#include <sstream>
30
31namespace spdlog::level
32{
34 spdlog::level::level_enum,
35 {{spdlog::level::level_enum::trace, "trace"},
36 {spdlog::level::level_enum::debug, "debug"},
37 {spdlog::level::level_enum::info, "info"},
38 {spdlog::level::level_enum::warn, "warning"},
39 {spdlog::level::level_enum::err, "error"},
40 {spdlog::level::level_enum::critical, "critical"},
41 {spdlog::level::level_enum::off, "off"},
42 {spdlog::level::level_enum::trace, 0},
43 {spdlog::level::level_enum::debug, 1},
44 {spdlog::level::level_enum::info, 2},
45 {spdlog::level::level_enum::warn, 3},
46 {spdlog::level::level_enum::err, 3},
47 {spdlog::level::level_enum::critical, 4},
48 {spdlog::level::level_enum::off, 5}})
49}
50
51namespace polyfem
52{
53 using namespace problem;
54 using namespace utils;
55
57 {
58 using namespace polysolve;
59
61
63 }
64
66 const std::string &log_file,
67 const spdlog::level::level_enum log_level,
68 const spdlog::level::level_enum file_log_level,
69 const bool is_quiet)
70 {
71 std::vector<spdlog::sink_ptr> sinks;
72
73 if (!is_quiet)
74 {
75 console_sink_ = std::make_shared<spdlog::sinks::stdout_color_sink_mt>();
76 sinks.emplace_back(console_sink_);
77 }
78
79 if (!log_file.empty())
80 {
81 file_sink_ = std::make_shared<spdlog::sinks::basic_file_sink_mt>(log_file, /*truncate=*/true);
82 // Set the file sink separately from the console so it can save all messages
83 file_sink_->set_level(file_log_level);
84 sinks.push_back(file_sink_);
85 }
86
87 init_logger(sinks, log_level);
88 spdlog::flush_every(std::chrono::seconds(3));
89 }
90
91 void State::init_logger(std::ostream &os, const spdlog::level::level_enum log_level)
92 {
93 std::vector<spdlog::sink_ptr> sinks;
94 sinks.emplace_back(std::make_shared<spdlog::sinks::ostream_sink_mt>(os, false));
95 init_logger(sinks, log_level);
96 }
97
99 const std::vector<spdlog::sink_ptr> &sinks,
100 const spdlog::level::level_enum log_level)
101 {
102 set_logger(std::make_shared<spdlog::logger>("polyfem", sinks.begin(), sinks.end()));
104
105 ipc::set_logger(std::make_shared<spdlog::logger>("ipctk", sinks.begin(), sinks.end()));
106
107 wmtk::set_logger(std::make_shared<spdlog::logger>("wmtk", sinks.begin(), sinks.end()));
108
109 // Set the logger at the lowest level, so all messages are passed to the sinks
110 logger().set_level(spdlog::level::trace);
111 ipc::logger().set_level(spdlog::level::trace);
112 wmtk::logger().set_level(spdlog::level::trace);
113
114 set_log_level(log_level);
115 }
116
117 void State::set_log_level(const spdlog::level::level_enum log_level)
118 {
119 spdlog::set_level(log_level);
120 if (console_sink_)
121 {
122 // Set only the level of the console
123 console_sink_->set_level(log_level); // Shared by all loggers
124 }
125 else
126 {
127 // Set the level of all sinks
128 logger().set_level(log_level);
129 ipc::logger().set_level(log_level);
130 }
131 }
132
133 void State::init(const json &p_args_in, const bool strict_validation)
134 {
135 json args_in = p_args_in; // mutable copy
136
137 has_constraints_ = p_args_in.contains("constraints");
138
139 apply_common_params(args_in);
140
141 // CHECK validity json
142 json rules;
143 jse::JSE jse;
144 {
145 jse.strict = strict_validation;
146 const std::string polyfem_input_spec = POLYFEM_INPUT_SPEC;
147 std::ifstream file(polyfem_input_spec);
148
149 if (file.is_open())
150 file >> rules;
151 else
152 {
153 logger().error("unable to open {} rules", polyfem_input_spec);
154 throw std::runtime_error("Invald spec file");
155 }
156
157 jse.include_directories.push_back(POLYFEM_JSON_SPEC_DIR);
158 jse.include_directories.push_back(POLYSOLVE_JSON_SPEC_DIR);
159 rules = jse.inject_include(rules);
160
161 polysolve::linear::Solver::apply_default_solver(rules, "/solver/linear");
162 polysolve::linear::Solver::apply_default_solver(rules, "/solver/adjoint_linear");
163 }
164
165 polysolve::linear::Solver::select_valid_solver(args_in["solver"]["linear"], logger());
166 if (args_in["solver"]["adjoint_linear"].is_null())
167 args_in["solver"]["adjoint_linear"] = args_in["solver"]["linear"];
168 else
169 polysolve::linear::Solver::select_valid_solver(args_in["solver"]["adjoint_linear"], logger());
170
171 // Use the /solver/nonlinear settings as the default for /solver/augmented_lagrangian/nonlinear
172 if (args_in.contains("/solver/nonlinear"_json_pointer))
173 {
174 if (args_in.contains("/solver/augmented_lagrangian/nonlinear"_json_pointer))
175 {
176 assert(args_in["solver"]["augmented_lagrangian"]["nonlinear"].is_object());
177 // Merge the augmented lagrangian settings into the nonlinear settings,
178 // and then replace the augmented lagrangian settings with the merged settings.
179 json nonlinear = args_in["solver"]["nonlinear"]; // copy
180 nonlinear.merge_patch(args_in["solver"]["augmented_lagrangian"]["nonlinear"]);
181 args_in["solver"]["augmented_lagrangian"]["nonlinear"] = nonlinear;
182 }
183 else
184 {
185 // Copy the nonlinear settings to the augmented_lagrangian settings
186 args_in["solver"]["augmented_lagrangian"]["nonlinear"] = args_in["solver"]["nonlinear"];
187 }
188 }
189 const bool valid_input = jse.verify_json(args_in, rules);
190
191 if (!valid_input)
192 {
193 logger().error("invalid input json:\n{}", jse.log2str());
194 throw std::runtime_error("Invald input json file");
195 }
196 // end of check
197
198 this->args = jse.inject_defaults(args_in, rules);
199 units.init(this->args["units"]);
200
201 if (!args_in.contains("/space/advanced/bc_method"_json_pointer) && this->args["space"]["basis_type"] != "Lagrange")
202 {
203 logger().warn("Setting bc method to lsq for non-Lagrange basis");
204 this->args["space"]["advanced"]["bc_method"] = "lsq";
205 }
206
207 // Save output directory and resolve output paths dynamically
208 const std::string output_dir = resolve_input_path(this->args["output"]["directory"]);
209 if (!output_dir.empty())
210 {
211 std::filesystem::create_directories(output_dir);
212 }
213 this->output_dir = output_dir;
214
215 std::string out_path_log = this->args["output"]["log"]["path"];
216 if (!out_path_log.empty())
217 {
218 out_path_log = resolve_output_path(out_path_log);
219 }
220
221 for (auto &path : this->args["constraints"]["hard"])
222 {
223 path = resolve_input_path(path);
224 }
225 for (auto &path : this->args["constraints"]["soft"])
226 {
227 path["data"] = resolve_input_path(path["data"]);
228 }
229
231 out_path_log,
232 this->args["output"]["log"]["level"],
233 this->args["output"]["log"]["file_level"],
234 this->args["output"]["log"]["quiet"]);
235
236 logger().info("Saving output to {}", output_dir);
237
238 const unsigned int thread_in = this->args["solver"]["max_threads"];
239 set_max_threads(thread_in);
240
241 has_dhat = args_in["contact"].contains("dhat");
242
243 init_time();
244
245 if (is_contact_enabled())
246 {
247 if (args["solver"]["contact"]["friction_iterations"] == 0)
248 {
249 logger().info("specified friction_iterations is 0; disabling friction");
250 args["contact"]["friction_coefficient"] = 0.0;
251 }
252 else if (args["solver"]["contact"]["friction_iterations"] < 0)
253 {
254 args["solver"]["contact"]["friction_iterations"] = std::numeric_limits<int>::max();
255 }
256 if (args["contact"]["friction_coefficient"] == 0.0)
257 {
258 args["solver"]["contact"]["friction_iterations"] = 0;
259 }
260 }
261 else
262 {
263 args["solver"]["contact"]["friction_iterations"] = 0;
264 args["contact"]["friction_coefficient"] = 0;
265 args["contact"]["periodic"] = false;
266 }
267
268 const std::string formulation = this->formulation();
270 assert(assembler->name() == formulation);
271 mass_matrix_assembler = std::make_shared<assembler::Mass>();
273
274 if (!other_name.empty())
275 {
278 }
279
280 if (args["solver"]["advanced"]["check_inversion"] == "Conservative")
281 {
282 if (auto elastic_assembler = std::dynamic_pointer_cast<assembler::ElasticityAssembler>(assembler))
283 elastic_assembler->set_use_robust_jacobian();
284 }
285
286 if (!args.contains("preset_problem"))
287 {
288 if (!assembler->is_tensor())
289 problem = std::make_shared<assembler::GenericScalarProblem>("GenericScalar");
290 else
291 problem = std::make_shared<assembler::GenericTensorProblem>("GenericTensor");
292
293 problem->clear();
294 if (!args["time"].is_null())
295 {
296 const auto tmp = R"({"is_time_dependent": true})"_json;
297 problem->set_parameters(tmp);
298 }
299 // important for the BC
300
301 auto bc = args["boundary_conditions"];
302 bc["root_path"] = root_path();
303 problem->set_parameters(bc);
304 problem->set_parameters(args["initial_conditions"]);
305
306 problem->set_parameters(args["output"]);
307 }
308 else
309 {
310 if (args["preset_problem"]["type"] == "Kernel")
311 {
312 problem = std::make_shared<KernelProblem>("Kernel", *assembler);
313 problem->clear();
314 KernelProblem &kprob = *dynamic_cast<KernelProblem *>(problem.get());
315 }
316 else
317 {
318 problem = ProblemFactory::factory().get_problem(args["preset_problem"]["type"]);
319 problem->clear();
320 }
321 // important for the BC
322 problem->set_parameters(args["preset_problem"]);
323 }
324
325 problem->set_units(*assembler, units);
326
328 {
329 if (is_contact_enabled())
330 {
331 if (!args["contact"]["use_convergent_formulation"])
332 {
333 args["contact"]["use_convergent_formulation"] = true;
334 logger().info("Use convergent formulation for differentiable contact...");
335 }
336 if (args["/solver/contact/barrier_stiffness"_json_pointer].is_string())
337 {
338 logger().error("Only constant barrier stiffness is supported in differentiable contact!");
339 }
340 }
341
342 if (args.contains("boundary_conditions") && args["boundary_conditions"].contains("rhs"))
343 {
344 json rhs = args["boundary_conditions"]["rhs"];
345 if ((rhs.is_array() && rhs.size() > 0 && rhs[0].is_string()) || rhs.is_string())
346 logger().error("Only constant rhs over space is supported in differentiable code!");
347 }
348 }
349 }
350
351 void State::set_max_threads(const int max_threads)
352 {
353 NThread::get().set_num_threads(max_threads);
354 }
355
357 {
358 if (!is_param_valid(args, "time"))
359 return;
360
361 const double t0 = Units::convert(args["time"]["t0"], units.time());
362 double tend, dt;
363 int time_steps;
364
365 // from "tend", "dt", "time_steps" only two can be used at a time
366 const int num_valid = is_param_valid(args["time"], "tend")
367 + is_param_valid(args["time"], "dt")
368 + is_param_valid(args["time"], "time_steps");
369 if (num_valid < 2)
370 {
371 log_and_throw_error("Exactly two of (tend, dt, time_steps) must be specified");
372 }
373 else if (num_valid == 2)
374 {
375 if (is_param_valid(args["time"], "tend"))
376 {
377 tend = Units::convert(args["time"]["tend"], units.time());
378 assert(tend > t0);
379 if (is_param_valid(args["time"], "dt"))
380 {
381 dt = Units::convert(args["time"]["dt"], units.time());
382 assert(dt > 0);
383 time_steps = int(ceil((tend - t0) / dt));
384 assert(time_steps > 0);
385 }
386 else if (is_param_valid(args["time"], "time_steps"))
387 {
388 time_steps = args["time"]["time_steps"];
389 assert(time_steps > 0);
390 dt = (tend - t0) / time_steps;
391 assert(dt > 0);
392 }
393 else
394 {
395 throw std::runtime_error("This code should be unreachable!");
396 }
397 }
398 else if (is_param_valid(args["time"], "dt"))
399 {
400 // tend is already confirmed to be invalid, so time_steps must be valid
401 assert(is_param_valid(args["time"], "time_steps"));
402
403 dt = Units::convert(args["time"]["dt"], units.time());
404 assert(dt > 0);
405
406 time_steps = args["time"]["time_steps"];
407 assert(time_steps > 0);
408
409 tend = t0 + time_steps * dt;
410 }
411 else
412 {
413 // tend and dt are already confirmed to be invalid
414 throw std::runtime_error("This code should be unreachable!");
415 }
416 }
417 else if (num_valid == 3)
418 {
419 tend = Units::convert(args["time"]["tend"], units.time());
420 dt = Units::convert(args["time"]["dt"], units.time());
421 time_steps = args["time"]["time_steps"];
422
423 // Check that all parameters agree
424 if (abs(t0 + dt * time_steps - tend) > 1e-12)
425 {
426 log_and_throw_error("Exactly two of (tend, dt, time_steps) must be specified");
427 }
428 }
429
430 // Store these for use later
431 args["time"]["tend"] = tend;
432 args["time"]["dt"] = dt;
433 args["time"]["time_steps"] = time_steps;
434
436
437 logger().info("t0={}, dt={}, tend={}", t0, dt, tend);
438 }
439
440 void State::set_materials(std::vector<std::shared_ptr<assembler::Assembler>> &assemblers) const
441 {
442 const int size = (assembler->is_tensor() || assembler->is_fluid()) ? mesh->dimension() : 1;
443 for (auto &a : assemblers)
444 a->set_size(size);
445
446 if (!utils::is_param_valid(args, "materials"))
447 return;
448
449 if (!args["materials"].is_array() && args["materials"]["type"] == "AMIPSAutodiff")
450 {
451 json transform_params = {};
452 transform_params["canonical_transformation"] = json::array();
453 if (!mesh->is_volume())
454 {
455 Eigen::MatrixXd regular_tri(3, 3);
456 regular_tri << 0, 0, 1,
457 1, 0, 1,
458 1. / 2., std::sqrt(3) / 2., 1;
459 regular_tri.transposeInPlace();
460 Eigen::MatrixXd regular_tri_inv = regular_tri.inverse();
461
462 const auto &mesh2d = *dynamic_cast<mesh::Mesh2D *>(mesh.get());
463 for (int e = 0; e < mesh->n_elements(); e++)
464 {
465 Eigen::MatrixXd transform;
466 mesh2d.compute_face_jacobian(e, regular_tri_inv, transform);
467 transform_params["canonical_transformation"].push_back(json({
468 {
469 transform(0, 0),
470 transform(0, 1),
471 },
472 {
473 transform(1, 0),
474 transform(1, 1),
475 },
476 }));
477 }
478 }
479 else
480 {
481 Eigen::MatrixXd regular_tet(4, 4);
482 regular_tet << 0, 0, 0, 1,
483 1, 0, 0, 1,
484 1. / 2., std::sqrt(3) / 2., 0, 1,
485 1. / 2., 1. / 2. / std::sqrt(3), std::sqrt(3) / 2., 1;
486 regular_tet.transposeInPlace();
487 Eigen::MatrixXd regular_tet_inv = regular_tet.inverse();
488
489 const auto &mesh3d = *dynamic_cast<mesh::Mesh3D *>(mesh.get());
490 for (int e = 0; e < mesh->n_elements(); e++)
491 {
492 Eigen::MatrixXd transform;
493 mesh3d.compute_cell_jacobian(e, regular_tet_inv, transform);
494 transform_params["canonical_transformation"].push_back(json({
495 {
496 transform(0, 0),
497 transform(0, 1),
498 transform(0, 2),
499 },
500 {
501 transform(1, 0),
502 transform(1, 1),
503 transform(1, 2),
504 },
505 {
506 transform(2, 0),
507 transform(2, 1),
508 transform(2, 2),
509 },
510 }));
511 }
512 }
513 transform_params["solve_displacement"] = true;
514 assembler->set_materials({}, transform_params, units);
515
516 return;
517 }
518
519 std::vector<int> body_ids(mesh->n_elements());
520 for (int i = 0; i < mesh->n_elements(); ++i)
521 body_ids[i] = mesh->get_body_id(i);
522
523 for (auto &a : assemblers)
524 a->set_materials(body_ids, args["materials"], units);
525 }
526
528 {
529 const int size = (this->assembler->is_tensor() || this->assembler->is_fluid()) ? this->mesh->dimension() : 1;
530 assembler.set_size(size);
531
532 if (!utils::is_param_valid(args, "materials"))
533 return;
534
535 std::vector<int> body_ids(mesh->n_elements());
536 for (int i = 0; i < mesh->n_elements(); ++i)
537 body_ids[i] = mesh->get_body_id(i);
538
539 assembler.set_materials(body_ids, args["materials"], units);
540 }
541
542} // namespace polyfem
std::string formulation() const
return the formulation (checks if the problem is scalar or not and deals with multiphysics)
Definition State.cpp:328
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 ".")
std::shared_ptr< assembler::Assembler > assembler
assemblers
Definition State.hpp:155
std::string resolve_output_path(const std::string &path) const
Resolve output path relative to output_dir if the path is not absolute.
std::string output_dir
Directory for output files.
Definition State.hpp:586
solver::CacheLevel optimization_enabled
Definition State.hpp:654
void set_materials(std::vector< std::shared_ptr< assembler::Assembler > > &assemblers) const
set the material and the problem dimension
void set_max_threads(const int max_threads=std::numeric_limits< int >::max())
std::shared_ptr< assembler::Assembler > pressure_assembler
Definition State.hpp:160
std::shared_ptr< assembler::Mass > mass_matrix_assembler
Definition State.hpp:157
bool has_dhat
stores if input json contains dhat
Definition State.hpp:579
std::unique_ptr< mesh::Mesh > mesh
current mesh, it can be a Mesh2D or Mesh3D
Definition State.hpp:471
std::shared_ptr< assembler::Problem > problem
current problem, it contains rhs and bc
Definition State.hpp:168
bool has_constraints_
Definition State.hpp:415
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
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 init_time()
initialize time settings if args contains "time"
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.
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
bool is_contact_enabled() const
does the simulation have contact
Definition State.hpp:556
std::shared_ptr< assembler::MixedAssembler > mixed_assembler
Definition State.hpp:159
Eigen::MatrixXd rhs
System right-hand side.
Definition State.hpp:207
void init(const json &json)
Definition Units.cpp:9
double characteristic_length() const
Definition Units.hpp:22
static double convert(const json &val, const std::string &unit_type)
Definition Units.cpp:31
const std::string & time() const
Definition Units.hpp:21
virtual bool is_tensor() const
virtual void set_size(const int size)
Definition Assembler.hpp:64
virtual bool is_fluid() const
static std::shared_ptr< MixedAssembler > make_mixed_assembler(const std::string &formulation)
static std::string other_assembler_name(const std::string &formulation)
static std::shared_ptr< Assembler > make_assembler(const std::string &formulation)
void compute_face_jacobian(const int el_id, const Eigen::MatrixXd &reference_map, Eigen::MatrixXd &jacobian) const
Definition Mesh2D.cpp:101
void compute_cell_jacobian(const int el_id, const Eigen::MatrixXd &reference_map, Eigen::MatrixXd &jacobian) const
Definition Mesh3D.cpp:399
static const ProblemFactory & factory()
std::shared_ptr< assembler::Problem > get_problem(const std::string &problem) const
void set_logger(spdlog::logger &logger)
static GeogramUtils & instance()
static NThread & get()
Definition par_for.hpp:19
void set_num_threads(const int max_threads)
Definition par_for.hpp:27
NLOHMANN_JSON_SERIALIZE_ENUM(CollisionProxyTessellation, {{CollisionProxyTessellation::REGULAR, "regular"}, {CollisionProxyTessellation::IRREGULAR, "irregular"}})
bool is_param_valid(const json &params, const std::string &key)
Determine if a key exists and is non-null in a json object.
void apply_common_params(json &args)
Definition JSONUtils.cpp:14
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:44
nlohmann::json json
Definition Common.hpp:9
void set_logger(std::shared_ptr< spdlog::logger > p_logger)
Setup a logger object to be used by Polyfem.
Definition Logger.cpp:62
void log_and_throw_error(const std::string &msg)
Definition Logger.cpp:73