43 if (state.
args[
"time"][
"integrator"].is_string())
45 if (state.
args[
"time"][
"integrator"][
"type"] ==
"ImplicitEuler")
47 if (state.
args[
"time"][
"integrator"][
"type"] ==
"BDF")
48 return state.
args[
"time"][
"integrator"][
"steps"].get<
int>();
54 double dot(
const Eigen::MatrixXd &A,
const Eigen::MatrixXd &B) {
return (A.array() * B.array()).sum(); }
56 class LocalThreadScalarStorage
63 LocalThreadScalarStorage()
69 class LocalThreadVecStorage
76 LocalThreadVecStorage(
const int size)
86 T triangle_area(
const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &
V)
88 Eigen::Matrix<T, Eigen::Dynamic, 1> l1 =
V.row(1) -
V.row(0);
89 Eigen::Matrix<T, Eigen::Dynamic, 1> l2 =
V.row(2) -
V.row(0);
90 T area = 0.5 * sqrt(pow(l1(1) * l2(2) - l1(2) * l2(1), 2) + pow(l1(0) * l2(2) - l1(2) * l2(0), 2) + pow(l1(1) * l2(0) - l1(0) * l2(1), 2));
94 Eigen::MatrixXd triangle_area_grad(
const Eigen::MatrixXd &
F)
97 Eigen::Matrix<Diff, Eigen::Dynamic, Eigen::Dynamic> full_diff(
F.rows(),
F.cols());
98 for (
int i = 0; i <
F.rows(); i++)
99 for (
int j = 0; j <
F.cols(); j++)
100 full_diff(i, j) =
Diff(i + j *
F.rows(),
F(i, j));
103 Eigen::MatrixXd
grad(
F.rows(),
F.cols());
104 for (
int i = 0; i <
F.rows(); ++i)
105 for (
int j = 0; j <
F.cols(); ++j)
106 grad(i, j) = reduced_diff.getGradient()(i + j *
F.rows());
111 template <
typename T>
112 T line_length(
const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &
V)
114 Eigen::Matrix<T, Eigen::Dynamic, 1> L =
V.row(1) -
V.row(0);
119 Eigen::MatrixXd line_length_grad(
const Eigen::MatrixXd &
F)
122 Eigen::Matrix<Diff, Eigen::Dynamic, Eigen::Dynamic> full_diff(
F.rows(),
F.cols());
123 for (
int i = 0; i <
F.rows(); i++)
124 for (
int j = 0; j <
F.cols(); j++)
125 full_diff(i, j) =
Diff(i + j *
F.rows(),
F(i, j));
126 auto reduced_diff = line_length(full_diff);
128 Eigen::MatrixXd
grad(
F.rows(),
F.cols());
129 for (
int i = 0; i <
F.rows(); ++i)
130 for (
int j = 0; j <
F.cols(); ++j)
131 grad(i, j) = reduced_diff.getGradient()(i + j *
F.rows());
136 template <
typename T>
137 Eigen::Matrix<T, 2, 1> edge_normal(
const Eigen::Matrix<T, 4, 1> &
V)
139 Eigen::Matrix<T, 2, 1> v1 =
V.segment(0, 2);
140 Eigen::Matrix<T, 2, 1> v2 =
V.segment(2, 2);
141 Eigen::Matrix<T, 2, 1> normal = v1 - v2;
143 normal = normal / normal.norm();
147 template <
typename T>
148 Eigen::Matrix<T, 3, 1> face_normal(
const Eigen::Matrix<T, 9, 1> &
V)
150 Eigen::Matrix<T, 3, 1> v1 =
V.segment(0, 3);
151 Eigen::Matrix<T, 3, 1> v2 =
V.segment(3, 3);
152 Eigen::Matrix<T, 3, 1> v3 =
V.segment(6, 3);
153 Eigen::Matrix<T, 3, 1> normal = (v2 - v1).
cross(v3 - v1);
154 normal = normal / normal.norm();
158 Eigen::MatrixXd extract_lame_params(
const std::map<std::string, Assembler::ParamFunc> &lame_params,
const int e,
const int t,
const Eigen::MatrixXd &local_pts,
const Eigen::MatrixXd &pts)
160 Eigen::MatrixXd params = Eigen::MatrixXd::Zero(local_pts.rows(), 2);
162 auto search_lambda = lame_params.find(
"lambda");
163 auto search_mu = lame_params.find(
"mu");
165 if (search_lambda == lame_params.end() || search_mu == lame_params.end())
168 for (
int p = 0; p < local_pts.rows(); p++)
170 params(p, 0) = search_lambda->second(local_pts.row(p), pts.row(p), t, e);
171 params(p, 1) = search_mu->second(local_pts.row(p), pts.row(p), t, e);
181 const Eigen::MatrixXd &solution,
182 const std::set<int> &interested_ids,
186 const auto &bases = state.
bases;
189 const int dim = state.
mesh->dimension();
190 const int actual_dim = state.
problem->is_scalar() ? 1 : dim;
191 const int n_elements = int(bases.size());
192 const double t0 = state.
problem->is_time_dependent() ? state.
args[
"time"][
"t0"].get<
double>() : 0.0;
193 const double dt = state.
problem->is_time_dependent() ? state.
args[
"time"][
"dt"].get<
double>() : 0.0;
203 params.
t = dt * cur_step + t0;
204 params.
step = cur_step;
206 Eigen::MatrixXd u, grad_u;
207 Eigen::MatrixXd result;
209 for (
int e = start; e < end; ++e)
211 if (interested_ids.size() != 0 && interested_ids.find(state.
mesh->get_body_id(e)) == interested_ids.end())
227 local_storage.val += dot(result, local_storage.da);
230 for (
const LocalThreadScalarStorage &local_storage : storage)
231 integral += local_storage.val;
237 LocalThreadScalarStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
240 Eigen::MatrixXd points, normal;
241 Eigen::VectorXd weights;
243 Eigen::MatrixXd u, grad_u;
244 Eigen::MatrixXd result;
245 IntegrableFunctional::ParameterType params;
246 params.t = dt * cur_step + t0;
247 params.step = cur_step;
249 for (int lb_id = start; lb_id < end; ++lb_id)
251 const auto &lb = state.total_local_boundary[lb_id];
252 const int e = lb.element_id();
254 for (int i = 0; i < lb.size(); i++)
256 const int global_primitive_id = lb.global_primitive_id(i);
257 if (interested_ids.size() != 0 && interested_ids.find(state.mesh->get_boundary_id(global_primitive_id)) == interested_ids.end())
260 utils::BoundarySampler::boundary_quadrature(lb, state.n_boundary_samples(), *state.mesh, i, false, uv, points, normal, weights);
262 assembler::ElementAssemblyValues &vals = local_storage.vals;
263 vals.compute(e, state.mesh->is_volume(), points, bases[e], gbases[e]);
264 io::Evaluator::interpolate_at_local_vals(e, dim, actual_dim, vals, solution, u, grad_u);
266 const Eigen::MatrixXd lame_params = extract_lame_params(state.assembler->parameters(), e, params.t, points, vals.val);
269 params.body_id = state.mesh->get_body_id(e);
270 params.boundary_id = state.mesh->get_boundary_id(global_primitive_id);
271 j.evaluate(lame_params, points, vals.val, u, grad_u, normal, vals, params, result);
273 local_storage.val += dot(result, weights);
277 for (
const LocalThreadScalarStorage &local_storage : storage)
278 integral += local_storage.val;
282 std::vector<bool> traversed(state.
n_bases,
false);
284 params.
t = dt * cur_step + t0;
285 params.
step = cur_step;
286 for (
int e = 0; e < bases.size(); e++)
288 const auto &bs = bases[e];
289 for (
int i = 0; i < bs.bases.size(); i++)
291 const auto &b = bs.bases[i];
292 assert(b.global().size() == 1);
293 const auto &g = b.global()[0];
294 if (traversed[g.index])
297 const Eigen::MatrixXd lame_params = extract_lame_params(state.
assembler->parameters(), e, params.
t, Eigen::MatrixXd::Zero(1, dim) , g.node);
299 params.
node = g.index;
303 j.evaluate(lame_params, Eigen::MatrixXd::Zero(1, dim) , g.node, solution.block(g.index * dim, 0, dim, 1).transpose(), Eigen::MatrixXd::Zero(1, dim * actual_dim) , Eigen::MatrixXd::Zero(0, 0) ,
assembler::ElementAssemblyValues(), params,
val);
305 traversed[g.index] =
true;
313 void AdjointTools::compute_shape_derivative_functional_term(
315 const Eigen::MatrixXd &solution,
317 const std::set<int> &interested_ids,
319 Eigen::VectorXd &term,
320 const int cur_time_step)
323 const auto &bases = state.
bases;
324 const int dim = state.
mesh->dimension();
325 const int actual_dim = state.
problem->is_scalar() ? 1 : dim;
326 const double t0 = state.
problem->is_time_dependent() ? state.
args[
"time"][
"t0"].get<
double>() : 0.0;
327 const double dt = state.
problem->is_time_dependent() ? state.
args[
"time"][
"dt"].get<
double>() : 0.0;
329 const int n_elements = int(bases.size());
332 auto storage = utils::create_thread_storage(LocalThreadVecStorage(term.size()));
334 if (spatial_integral_type == SpatialIntegralType::Volume)
336 utils::maybe_parallel_for(n_elements, [&](
int start,
int end,
int thread_id) {
337 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
339 Eigen::MatrixXd u, grad_u, j_val, dj_dgradu, dj_dx;
342 params.
t = cur_time_step * dt + t0;
343 params.
step = cur_time_step;
345 for (
int e = start; e < end; ++e)
347 if (interested_ids.size() != 0 && interested_ids.find(state.
mesh->get_body_id(e)) == interested_ids.end())
352 io::Evaluator::interpolate_at_local_vals(e, dim, actual_dim,
vals, solution, u, grad_u);
373 Eigen::MatrixXd tau_q, grad_u_q;
374 for (
auto &v :
gvals.basis_values)
376 for (
int q = 0; q < local_storage.da.size(); ++q)
378 local_storage.vec.block(v.global[0].index * dim, 0, dim, 1) += (j_val(q) * local_storage.da(q)) * v.grad_t_m.row(q).transpose();
381 local_storage.vec.block(v.global[0].index * dim, 0, dim, 1) += (v.val(q) * local_storage.da(q)) * dj_dx.row(q).transpose();
385 if (dim == actual_dim)
392 tau_q = dj_dgradu.row(q);
393 grad_u_q = grad_u.row(q);
395 for (
int d = 0; d < dim; d++)
396 local_storage.vec(v.global[0].index * dim + d) += -dot(tau_q, grad_u_q.col(d) * v.grad_t_m.row(q)) * local_storage.da(q);
403 else if (spatial_integral_type == SpatialIntegralType::Surface)
405 utils::maybe_parallel_for(state.
total_local_boundary.size(), [&](
int start,
int end,
int thread_id) {
406 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
408 Eigen::MatrixXd uv, points, normal;
409 Eigen::VectorXd &weights = local_storage.da;
411 Eigen::MatrixXd u, grad_u, x, grad_x, j_val, dj_dgradu, dj_dgradx, dj_dx;
413 IntegrableFunctional::ParameterType params;
414 params.t = cur_time_step * dt + t0;
415 params.step = cur_time_step;
417 for (int lb_id = start; lb_id < end; ++lb_id)
419 const auto &lb = state.total_local_boundary[lb_id];
420 const int e = lb.element_id();
422 for (int i = 0; i < lb.size(); i++)
424 const int global_primitive_id = lb.global_primitive_id(i);
425 if (interested_ids.size() != 0 && interested_ids.find(state.mesh->get_boundary_id(global_primitive_id)) == interested_ids.end())
428 utils::BoundarySampler::boundary_quadrature(lb, state.n_boundary_samples(), *state.mesh, i, false, uv, points, normal, weights);
430 assembler::ElementAssemblyValues &vals = local_storage.vals;
431 io::Evaluator::interpolate_at_local_vals(*state.mesh, state.problem->is_scalar(), bases, gbases, e, points, solution, u, grad_u);
434 vals.compute(e, state.mesh->is_volume(), points, gbases[e], gbases[e]);
438 const int n_loc_bases_ = int(vals.basis_values.size());
440 const Eigen::MatrixXd lame_params = extract_lame_params(state.assembler->parameters(), e, params.t, points, vals.val);
443 params.body_id = state.mesh->get_body_id(e);
444 params.boundary_id = state.mesh->get_boundary_id(global_primitive_id);
446 j.evaluate(lame_params, points, vals.val, u, grad_u, normal, vals, params, j_val);
447 j_val = j_val.array().colwise() * weights.array();
449 if (j.depend_on_gradu())
451 j.dj_dgradu(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_dgradu);
452 dj_dgradu = dj_dgradu.array().colwise() * weights.array();
455 if (j.depend_on_gradx())
457 j.dj_dgradx(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_dgradx);
458 dj_dgradx = dj_dgradx.array().colwise() * weights.array();
463 j.dj_dx(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_dx);
464 dj_dx = dj_dx.array().colwise() * weights.array();
467 const auto nodes = gbases[e].local_nodes_for_primitive(lb.global_primitive_id(i), *state.mesh);
469 if (nodes.size() != dim)
470 log_and_throw_adjoint_error(
"Only linear geometry is supported in differentiable surface integral functional!");
472 Eigen::MatrixXd velocity_div_mat;
473 if (state.mesh->is_volume())
476 for (int d = 0; d < 3; d++)
477 V.row(d) = gbases[e].bases[nodes(d)].global()[0].node;
478 velocity_div_mat = face_velocity_divergence(V);
483 for (int d = 0; d < 2; d++)
484 V.row(d) = gbases[e].bases[nodes(d)].global()[0].node;
485 velocity_div_mat = edge_velocity_divergence(V);
488 Eigen::MatrixXd grad_u_q, tau_q, grad_x_q;
489 for (long n = 0; n < nodes.size(); ++n)
491 const assembler::AssemblyValues &v = vals.basis_values[nodes(n)];
493 local_storage.vec.block(v.global[0].index * dim, 0, dim, 1) += j_val.sum() * velocity_div_mat.row(n).transpose();
496 for (long n = 0; n < n_loc_bases_; ++n)
498 const assembler::AssemblyValues &v = vals.basis_values[n];
501 local_storage.vec.block(v.global[0].index * dim, 0, dim, 1) += dj_dx.transpose() * v.val;
504 if (j.depend_on_gradu())
506 for (int q = 0; q < weights.size(); ++q)
508 if (dim == actual_dim)
510 vector2matrix(grad_u.row(q), grad_u_q);
511 vector2matrix(dj_dgradu.row(q), tau_q);
515 grad_u_q = grad_u.row(q);
516 tau_q = dj_dgradu.row(q);
519 for (int d = 0; d < dim; d++)
520 local_storage.vec(v.global[0].index * dim + d) += -dot(tau_q, grad_u_q.col(d) * v.grad_t_m.row(q));
524 if (j.depend_on_gradx())
526 for (int d = 0; d < dim; d++)
528 for (int q = 0; q < weights.size(); ++q)
529 local_storage.vec(v.global[0].index * dim + d) += dot(dj_dgradx.block(q, d * dim, 1, dim), v.grad.row(q));
537 else if (spatial_integral_type == SpatialIntegralType::VertexSum)
541 for (
const LocalThreadVecStorage &local_storage : storage)
542 term += local_storage.
vec;
544 term = utils::flatten(utils::unflatten(term, dim)(state.
primitive_to_node(), Eigen::all));
313 void AdjointTools::compute_shape_derivative_functional_term( {
…}
547 void AdjointTools::dJ_shape_static_adjoint_term(
549 const Eigen::MatrixXd &sol,
550 const Eigen::MatrixXd &adjoint,
551 Eigen::VectorXd &one_form)
553 Eigen::VectorXd elasticity_term, rhs_term, pressure_term, contact_term;
556 Eigen::MatrixXd adjoint_zeroed = adjoint;
565 rhs_term.setZero(one_form.size());
573 pressure_term.setZero(one_form.size());
581 contact_term.setZero(elasticity_term.size());
582 one_form -= elasticity_term + rhs_term + pressure_term + contact_term;
584 one_form = utils::flatten(utils::unflatten(one_form, state.
mesh->dimension())(state.
primitive_to_node(), Eigen::all));
547 void AdjointTools::dJ_shape_static_adjoint_term( {
…}
587 void AdjointTools::dJ_shape_homogenization_adjoint_term(
589 const Eigen::MatrixXd &sol,
590 const Eigen::MatrixXd &adjoint,
591 Eigen::VectorXd &one_form)
593 Eigen::VectorXd elasticity_term, contact_term;
595 std::shared_ptr<NLHomoProblem> homo_problem = std::dynamic_pointer_cast<NLHomoProblem>(state.
solve_data.
nl_problem);
596 assert(homo_problem);
598 const int dim = state.
mesh->dimension();
601 const Eigen::MatrixXd affine_adjoint = homo_problem->reduced_to_disp_grad(adjoint,
true);
602 const Eigen::VectorXd full_adjoint = homo_problem->NLProblem::reduced_to_full(adjoint.topRows(homo_problem->reduced_size())) + io::Evaluator::generate_linear_field(state.
n_bases, state.
mesh_nodes, affine_adjoint);
612 contact_term.setZero(elasticity_term.size());
614 one_form = -(elasticity_term + contact_term);
616 Eigen::VectorXd force;
617 homo_problem->FullNLProblem::gradient(sol, force);
620 one_form = utils::flatten(utils::unflatten(one_form, dim)(state.
primitive_to_node(), Eigen::all));
587 void AdjointTools::dJ_shape_homogenization_adjoint_term( {
…}
623 void AdjointTools::dJ_periodic_shape_adjoint_term(
626 const Eigen::VectorXd &periodic_mesh_representation,
627 const Eigen::MatrixXd &sol,
628 const Eigen::MatrixXd &adjoint,
629 Eigen::VectorXd &one_form)
631 std::shared_ptr<NLHomoProblem> homo_problem = std::dynamic_pointer_cast<NLHomoProblem>(state.
solve_data.
nl_problem);
632 assert(homo_problem);
634 const Eigen::MatrixXd reduced_sol = homo_problem->full_to_reduced(sol, state.
diff_cached.
disp_grad());
635 const Eigen::VectorXd extended_sol = homo_problem->reduced_to_extended(reduced_sol);
637 const Eigen::VectorXd extended_adjoint = homo_problem->reduced_to_extended(adjoint,
true);
638 const Eigen::MatrixXd affine_adjoint = homo_problem->reduced_to_disp_grad(adjoint,
true);
639 const Eigen::VectorXd full_adjoint = homo_problem->NLProblem::reduced_to_full(adjoint.topRows(homo_problem->reduced_size())) + io::Evaluator::generate_linear_field(state.
n_bases, state.
mesh_nodes, affine_adjoint);
641 const int dim = state.
mesh->dimension();
646 homo_problem->set_project_to_psd(
false);
647 homo_problem->FullNLProblem::hessian(sol, hessian);
648 Eigen::VectorXd partial_term = full_adjoint.transpose() * hessian;
650 one_form -= utils::flatten(utils::unflatten(partial_term, dim)(state.
primitive_to_node(), Eigen::all));
652 one_form = periodic_mesh_map.
apply_jacobian(one_form, periodic_mesh_representation);
656 Eigen::VectorXd contact_term;
659 one_form -= contact_term;
623 void AdjointTools::dJ_periodic_shape_adjoint_term( {
…}
663 void AdjointTools::dJ_shape_transient_adjoint_term(
665 const Eigen::MatrixXd &adjoint_nu,
666 const Eigen::MatrixXd &adjoint_p,
667 Eigen::VectorXd &one_form)
669 const double t0 = state.
args[
"time"][
"t0"];
670 const double dt = state.
args[
"time"][
"dt"];
671 const int time_steps = state.
args[
"time"][
"time_steps"];
672 const int bdf_order = get_bdf_order(state);
674 Eigen::VectorXd elasticity_term, rhs_term, pressure_term, damping_term, mass_term, contact_term, friction_term;
677 Eigen::VectorXd cur_p, cur_nu;
678 for (
int i = time_steps; i > 0; --i)
680 const int real_order = std::min(bdf_order, i);
681 double beta = time_integrator::BDF::betas(real_order - 1);
682 double beta_dt = beta * dt;
683 const double t = i * dt + t0;
687 cur_p = adjoint_p.col(i);
688 cur_nu = adjoint_nu.col(i);
702 damping_term.setZero(mass_term.size());
711 contact_term.setZero(mass_term.size());
720 friction_term.setZero(mass_term.size());
723 one_form += beta_dt * (elasticity_term + rhs_term + pressure_term + damping_term + contact_term + friction_term + mass_term);
727 Eigen::VectorXd sum_alpha_p;
729 sum_alpha_p.setZero(adjoint_p.rows());
730 int num = std::min(bdf_order, time_steps);
731 for (
int j = 0; j < num; ++j)
733 int order = std::min(bdf_order - 1, j);
734 sum_alpha_p -= time_integrator::BDF::alphas(order)[j] * adjoint_p.col(j + 1);
740 one_form += mass_term;
742 one_form = utils::flatten(utils::unflatten(one_form, state.
mesh->dimension())(state.
primitive_to_node(), Eigen::all));
663 void AdjointTools::dJ_shape_transient_adjoint_term( {
…}
745 void AdjointTools::dJ_material_static_adjoint_term(
747 const Eigen::MatrixXd &sol,
748 const Eigen::MatrixXd &adjoint,
749 Eigen::VectorXd &one_form)
751 Eigen::MatrixXd adjoint_zeroed = adjoint;
745 void AdjointTools::dJ_material_static_adjoint_term( {
…}
756 void AdjointTools::dJ_material_transient_adjoint_term(
758 const Eigen::MatrixXd &adjoint_nu,
759 const Eigen::MatrixXd &adjoint_p,
760 Eigen::VectorXd &one_form)
762 const double t0 = state.
args[
"time"][
"t0"];
763 const double dt = state.
args[
"time"][
"dt"];
764 const int time_steps = state.
args[
"time"][
"time_steps"];
765 const int bdf_order = get_bdf_order(state);
767 one_form.setZero(state.
bases.size() * 2);
769 auto storage = utils::create_thread_storage(LocalThreadVecStorage(one_form.size()));
771 utils::maybe_parallel_for(time_steps, [&](
int start,
int end,
int thread_id) {
772 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
773 Eigen::VectorXd elasticity_term;
774 for (
int i_aux = start; i_aux < end; ++i_aux)
776 const int i = time_steps - i_aux;
777 const int real_order = std::min(bdf_order, i);
778 double beta_dt = time_integrator::BDF::betas(real_order - 1) * dt;
780 Eigen::VectorXd cur_p = adjoint_p.col(i);
784 local_storage.vec += beta_dt * elasticity_term;
788 for (
const LocalThreadVecStorage &local_storage : storage)
789 one_form += local_storage.vec;
756 void AdjointTools::dJ_material_transient_adjoint_term( {
…}
792 void AdjointTools::dJ_friction_transient_adjoint_term(
794 const Eigen::MatrixXd &adjoint_nu,
795 const Eigen::MatrixXd &adjoint_p,
796 Eigen::VectorXd &one_form)
798 const double dt = state.
args[
"time"][
"dt"];
800 const int time_steps = state.
args[
"time"][
"time_steps"];
801 const int dim = state.
mesh->dimension();
802 const int bdf_order = get_bdf_order(state);
806 std::shared_ptr<time_integrator::ImplicitTimeIntegrator> time_integrator =
807 time_integrator::ImplicitTimeIntegrator::construct_time_integrator(state.
args[
"time"][
"integrator"]);
809 Eigen::MatrixXd solution, velocity, acceleration;
815 const double dt = state.
args[
"time"][
"dt"];
816 time_integrator->init(solution, velocity, acceleration, dt);
819 for (
int t = 1; t <= time_steps; ++t)
821 const int real_order = std::min(bdf_order, t);
822 double beta = time_integrator::BDF::betas(real_order - 1);
824 const Eigen::MatrixXd surface_solution_prev = state.
collision_mesh.vertices(utils::unflatten(state.
diff_cached.
u(t - 1), dim));
828 time_integrator->update_quantities(state.
diff_cached.
u(t));
835 surface_solution_prev,
841 Eigen::VectorXd cur_p = adjoint_p.col(t);
844 one_form(0) += dot(cur_p, force) * beta * dt;
792 void AdjointTools::dJ_friction_transient_adjoint_term( {
…}
848 void AdjointTools::dJ_damping_transient_adjoint_term(
850 const Eigen::MatrixXd &adjoint_nu,
851 const Eigen::MatrixXd &adjoint_p,
852 Eigen::VectorXd &one_form)
854 const double t0 = state.
args[
"time"][
"t0"];
855 const double dt = state.
args[
"time"][
"dt"];
856 const int time_steps = state.
args[
"time"][
"time_steps"];
857 const int bdf_order = get_bdf_order(state);
861 auto storage = utils::create_thread_storage(LocalThreadVecStorage(one_form.size()));
863 utils::maybe_parallel_for(time_steps, [&](
int start,
int end,
int thread_id) {
864 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
865 Eigen::VectorXd damping_term;
866 for (
int t_aux = start; t_aux < end; ++t_aux)
868 const int t = time_steps - t_aux;
869 const int real_order = std::min(bdf_order, t);
870 const double beta = time_integrator::BDF::betas(real_order - 1);
872 Eigen::VectorXd cur_p = adjoint_p.col(t);
876 local_storage.vec += (beta * dt) * damping_term;
880 for (
const LocalThreadVecStorage &local_storage : storage)
881 one_form += local_storage.vec;
848 void AdjointTools::dJ_damping_transient_adjoint_term( {
…}
884 void AdjointTools::dJ_initial_condition_adjoint_term(
886 const Eigen::MatrixXd &adjoint_nu,
887 const Eigen::MatrixXd &adjoint_p,
888 Eigen::VectorXd &one_form)
890 const int ndof = state.
ndof();
891 one_form.setZero(ndof * 2);
894 one_form.segment(0, ndof) = -adjoint_nu.col(0);
895 one_form.segment(ndof, ndof) = -adjoint_p.col(0);
900 one_form(ndof + b) = 0;
884 void AdjointTools::dJ_initial_condition_adjoint_term( {
…}
904 void AdjointTools::dJ_dirichlet_static_adjoint_term(
906 const Eigen::MatrixXd &adjoint,
907 Eigen::VectorXd &one_form)
912 gradd_h.prune([&boundary_nodes_set](
const Eigen::Index &row,
const Eigen::Index &col,
const FullNLProblem::Scalar &value) {
915 if (boundary_nodes_set.find(row) == boundary_nodes_set.end())
919 one_form.setZero(state.
ndof());
922 one_form = utils::flatten(utils::unflatten(one_form, state.
mesh->dimension())(state.
primitive_to_node(), Eigen::all));
904 void AdjointTools::dJ_dirichlet_static_adjoint_term( {
…}
925 void AdjointTools::dJ_dirichlet_transient_adjoint_term(
927 const Eigen::MatrixXd &adjoint_nu,
928 const Eigen::MatrixXd &adjoint_p,
929 Eigen::VectorXd &one_form)
931 const double dt = state.
args[
"time"][
"dt"];
932 const int time_steps = state.
args[
"time"][
"time_steps"];
933 const int bdf_order = get_bdf_order(state);
938 one_form.setZero(time_steps * n_dirichlet_dof);
939 for (
int i = 1; i <= time_steps; ++i)
941 const int real_order = std::min(bdf_order, i);
942 const double beta_dt = time_integrator::BDF::betas(real_order - 1) * dt;
944 one_form.segment((i - 1) * n_dirichlet_dof, n_dirichlet_dof) = -(1. / beta_dt) * adjoint_p(state.
boundary_nodes, i);
925 void AdjointTools::dJ_dirichlet_transient_adjoint_term( {
…}
948 void AdjointTools::dJ_pressure_static_adjoint_term(
950 const std::vector<int> &boundary_ids,
951 const Eigen::MatrixXd &sol,
952 const Eigen::MatrixXd &adjoint,
953 Eigen::VectorXd &one_form)
955 const int n_pressure_dof = boundary_ids.size();
957 one_form.setZero(n_pressure_dof);
959 for (
int i = 0; i < boundary_ids.size(); ++i)
967 one_form(i) = pressure_term;
948 void AdjointTools::dJ_pressure_static_adjoint_term( {
…}
971 void AdjointTools::dJ_pressure_transient_adjoint_term(
973 const std::vector<int> &boundary_ids,
974 const Eigen::MatrixXd &adjoint_nu,
975 const Eigen::MatrixXd &adjoint_p,
976 Eigen::VectorXd &one_form)
978 const double t0 = state.
args[
"time"][
"t0"];
979 const double dt = state.
args[
"time"][
"dt"];
980 const int time_steps = state.
args[
"time"][
"time_steps"];
981 const int bdf_order = get_bdf_order(state);
983 const int n_pressure_dof = boundary_ids.size();
985 one_form.setZero(time_steps * n_pressure_dof);
986 Eigen::VectorXd cur_p, cur_nu;
987 for (
int i = time_steps; i > 0; --i)
989 const int real_order = std::min(bdf_order, i);
990 double beta = time_integrator::BDF::betas(real_order - 1);
991 double beta_dt = beta * dt;
992 const double t = i * dt + t0;
994 cur_p = adjoint_p.col(i);
995 cur_nu = adjoint_nu.col(i);
999 for (
int b = 0; b < boundary_ids.size(); ++b)
1007 one_form((i - 1) * n_pressure_dof + b) = -beta_dt * pressure_term;
971 void AdjointTools::dJ_pressure_transient_adjoint_term( {
…}
1012 void AdjointTools::dJ_du_step(
1015 const Eigen::MatrixXd &solution,
1016 const std::set<int> &interested_ids,
1019 Eigen::VectorXd &term)
1021 const auto &bases = state.
bases;
1024 const int dim = state.
mesh->dimension();
1025 const int actual_dim = state.
problem->is_scalar() ? 1 : dim;
1026 const int n_elements = int(bases.size());
1027 const double t0 = state.
problem->is_time_dependent() ? state.
args[
"time"][
"t0"].get<
double>() : 0.0;
1028 const double dt = state.
problem->is_time_dependent() ? state.
args[
"time"][
"dt"].get<
double>() : 0.0;
1030 term = Eigen::MatrixXd::Zero(state.
n_bases * actual_dim, 1);
1035 if (spatial_integral_type == SpatialIntegralType::Volume)
1037 auto storage = utils::create_thread_storage(LocalThreadVecStorage(term.size()));
1038 utils::maybe_parallel_for(n_elements, [&](
int start,
int end,
int thread_id) {
1039 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
1041 Eigen::MatrixXd u, grad_u;
1042 Eigen::MatrixXd lambda, mu;
1043 Eigen::MatrixXd dj_du, dj_dgradu, dj_dgradx;
1046 params.
t = dt * cur_step + t0;
1047 params.
step = cur_step;
1049 for (
int e = start; e < end; ++e)
1051 if (interested_ids.size() != 0 && interested_ids.find(state.
mesh->get_body_id(e)) == interested_ids.end())
1062 const int n_loc_bases_ = int(
vals.basis_values.size());
1064 io::Evaluator::interpolate_at_local_vals(e, dim, actual_dim,
vals, solution, u, grad_u);
1069 dj_dgradu.resize(0, 0);
1073 for (
int q = 0; q < dj_dgradu.rows(); q++)
1074 dj_dgradu.row(q) *= local_storage.da(q);
1081 for (
int q = 0; q < dj_du.rows(); q++)
1082 dj_du.row(q) *= local_storage.da(q);
1085 for (
int i = 0; i < n_loc_bases_; ++i)
1088 assert(v.
global.size() == 1);
1089 for (
int d = 0; d < actual_dim; d++)
1096 for (
int q = 0; q < local_storage.da.size(); ++q)
1097 val += dot(dj_dgradu.block(q, d * dim, 1, dim), v.
grad_t_m.row(q));
1103 for (
int q = 0; q < local_storage.da.size(); ++q)
1104 val += dj_du(q, d) * v.
val(q);
1106 local_storage.vec(v.
global[0].index * actual_dim + d) +=
val;
1111 for (
const LocalThreadVecStorage &local_storage : storage)
1112 term += local_storage.vec;
1114 else if (spatial_integral_type == SpatialIntegralType::Surface)
1116 auto storage = utils::create_thread_storage(LocalThreadVecStorage(term.size()));
1117 utils::maybe_parallel_for(state.
total_local_boundary.size(), [&](
int start,
int end,
int thread_id) {
1118 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
1120 Eigen::MatrixXd uv, samples, gtmp;
1121 Eigen::MatrixXd points, normal;
1122 Eigen::VectorXd weights;
1124 Eigen::MatrixXd u, grad_u;
1125 Eigen::MatrixXd lambda, mu;
1126 Eigen::MatrixXd dj_du, dj_dgradu, dj_dgradu_local;
1128 IntegrableFunctional::ParameterType params;
1129 params.t = dt * cur_step + t0;
1130 params.step = cur_step;
1132 for (int lb_id = start; lb_id < end; ++lb_id)
1134 const auto &lb = state.total_local_boundary[lb_id];
1135 const int e = lb.element_id();
1137 for (int i = 0; i < lb.size(); i++)
1139 const int global_primitive_id = lb.global_primitive_id(i);
1140 if (interested_ids.size() != 0 && interested_ids.find(state.mesh->get_boundary_id(global_primitive_id)) == interested_ids.end())
1143 utils::BoundarySampler::boundary_quadrature(lb, state.n_boundary_samples(), *state.mesh, i, false, uv, points, normal, weights);
1145 assembler::ElementAssemblyValues &vals = local_storage.vals;
1146 vals.compute(e, state.mesh->is_volume(), points, bases[e], gbases[e]);
1147 io::Evaluator::interpolate_at_local_vals(e, dim, actual_dim, vals, solution, u, grad_u);
1149 const Eigen::MatrixXd lame_params = extract_lame_params(state.assembler->parameters(), e, params.t, points, vals.val);
1153 const int n_loc_bases_ = int(vals.basis_values.size());
1156 params.body_id = state.mesh->get_body_id(e);
1157 params.boundary_id = state.mesh->get_boundary_id(global_primitive_id);
1159 dj_dgradu.resize(0, 0);
1160 if (j.depend_on_gradu())
1162 j.dj_dgradu(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_dgradu);
1163 for (int q = 0; q < dj_dgradu.rows(); q++)
1164 dj_dgradu.row(q) *= weights(q);
1167 dj_dgradu_local.resize(0, 0);
1168 if (j.depend_on_gradu_local())
1170 j.dj_dgradu_local(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_dgradu_local);
1171 for (int q = 0; q < dj_dgradu_local.rows(); q++)
1172 dj_dgradu_local.row(q) *= weights(q);
1176 if (j.depend_on_u())
1178 j.dj_du(lame_params, points, vals.val, u, grad_u, normal, vals, params, dj_du);
1179 for (int q = 0; q < dj_du.rows(); q++)
1180 dj_du.row(q) *= weights(q);
1183 for (int l = 0; l < lb.size(); ++l)
1185 const auto nodes = bases[e].local_nodes_for_primitive(lb.global_primitive_id(l), *state.mesh);
1187 for (long n = 0; n < nodes.size(); ++n)
1189 const assembler::AssemblyValues &v = vals.basis_values[nodes(n)];
1190 assert(v.global.size() == 1);
1191 for (int d = 0; d < actual_dim; d++)
1196 if (j.depend_on_gradu())
1198 for (int q = 0; q < weights.size(); ++q)
1199 val += dot(dj_dgradu.block(q, d * dim, 1, dim), v.grad_t_m.row(q));
1202 if (j.depend_on_gradu_local())
1204 for (int q = 0; q < weights.size(); ++q)
1205 val += dot(dj_dgradu_local.block(q, d * dim, 1, dim), v.grad.row(q));
1208 if (j.depend_on_u())
1210 for (int q = 0; q < weights.size(); ++q)
1211 val += dj_du(q, d) * v.val(q);
1213 local_storage.vec(v.global[0].index * actual_dim + d) += val;
1220 for (
const LocalThreadVecStorage &local_storage : storage)
1221 term += local_storage.vec;
1223 else if (spatial_integral_type == SpatialIntegralType::VertexSum)
1225 std::vector<bool> traversed(state.
n_bases,
false);
1227 params.
t = dt * cur_step + t0;
1228 params.
step = cur_step;
1229 for (
int e = 0; e < bases.size(); e++)
1231 const auto &bs = bases[e];
1232 for (
int i = 0; i < bs.bases.size(); i++)
1234 const auto &b = bs.bases[i];
1235 assert(b.global().size() == 1);
1236 const auto &g = b.global()[0];
1237 if (traversed[g.index])
1240 const Eigen::MatrixXd lame_params = extract_lame_params(state.
assembler->parameters(), e, params.
t, Eigen::MatrixXd::Zero(1, dim) , g.node);
1242 params.
node = g.index;
1245 Eigen::MatrixXd
val;
1246 j.dj_du(lame_params, Eigen::MatrixXd::Zero(1, dim) , g.node, solution.block(g.index * dim, 0, dim, 1).transpose(), Eigen::MatrixXd::Zero(1, dim * actual_dim) , Eigen::MatrixXd::Zero(0, 0) ,
assembler::ElementAssemblyValues(), params,
val);
1247 term.block(g.index * actual_dim, 0, actual_dim, 1) +=
val.transpose();
1248 traversed[g.index] =
true;
1012 void AdjointTools::dJ_du_step( {
…}
1254 Eigen::VectorXd AdjointTools::map_primitive_to_node_order(
const State &state,
const Eigen::VectorXd &primitives)
1256 int dim = state.
mesh->dimension();
1257 assert(primitives.size() == (state.
n_geom_bases * dim));
1258 Eigen::VectorXd nodes(primitives.size());
1261 nodes.segment(map[v] * dim, dim) = primitives.segment(v * dim, dim);
1254 Eigen::VectorXd AdjointTools::map_primitive_to_node_order(
const State &state,
const Eigen::VectorXd &primitives) {
…}
1265 Eigen::VectorXd AdjointTools::map_node_to_primitive_order(
const State &state,
const Eigen::VectorXd &nodes)
1267 int dim = state.
mesh->dimension();
1269 Eigen::VectorXd primitives(nodes.size());
1272 primitives.segment(map[v] * dim, dim) = nodes.segment(v * dim, dim);
1265 Eigen::VectorXd AdjointTools::map_node_to_primitive_order(
const State &state,
const Eigen::VectorXd &nodes) {
…}
1276 Eigen::MatrixXd AdjointTools::edge_normal_gradient(
const Eigen::MatrixXd &
V)
1279 Eigen::Matrix<Diff, 4, 1> full_diff(4, 1);
1280 for (
int i = 0; i < 2; i++)
1281 for (
int j = 0; j < 2; j++)
1282 full_diff(i * 2 + j) =
Diff(i * 2 + j,
V(i, j));
1283 auto reduced_diff = edge_normal(full_diff);
1285 Eigen::MatrixXd grad(2, 4);
1286 for (
int i = 0; i < 2; ++i)
1287 grad.row(i) = reduced_diff[i].getGradient();
1276 Eigen::MatrixXd AdjointTools::edge_normal_gradient(
const Eigen::MatrixXd &
V) {
…}
1292 Eigen::MatrixXd AdjointTools::face_normal_gradient(
const Eigen::MatrixXd &
V)
1295 Eigen::Matrix<Diff, 9, 1> full_diff(9, 1);
1296 for (
int i = 0; i < 3; i++)
1297 for (
int j = 0; j < 3; j++)
1298 full_diff(i * 3 + j) =
Diff(i * 3 + j,
V(i, j));
1299 auto reduced_diff = face_normal(full_diff);
1301 Eigen::MatrixXd grad(3, 9);
1302 for (
int i = 0; i < 3; ++i)
1303 grad.row(i) = reduced_diff[i].getGradient();
1292 Eigen::MatrixXd AdjointTools::face_normal_gradient(
const Eigen::MatrixXd &
V) {
…}
1308 Eigen::MatrixXd AdjointTools::edge_velocity_divergence(
const Eigen::MatrixXd &
V)
1310 return line_length_grad(
V) / line_length<double>(
V);
1308 Eigen::MatrixXd AdjointTools::edge_velocity_divergence(
const Eigen::MatrixXd &
V) {
…}
1313 Eigen::MatrixXd AdjointTools::face_velocity_divergence(
const Eigen::MatrixXd &
V)
1315 return triangle_area_grad(
V) / triangle_area<double>(
V);
1313 Eigen::MatrixXd AdjointTools::face_velocity_divergence(
const Eigen::MatrixXd &
V) {
…}
1318 double AdjointTools::triangle_jacobian(
const Eigen::VectorXd &v1,
const Eigen::VectorXd &v2,
const Eigen::VectorXd &v3)
1320 Eigen::VectorXd a = v2 - v1, b = v3 - v1;
1321 return a(0) * b(1) - b(0) * a(1);
1318 double AdjointTools::triangle_jacobian(
const Eigen::VectorXd &v1,
const Eigen::VectorXd &v2,
const Eigen::VectorXd &v3) {
…}
1324 double AdjointTools::tet_determinant(
const Eigen::VectorXd &v1,
const Eigen::VectorXd &v2,
const Eigen::VectorXd &v3,
const Eigen::VectorXd &v4)
1326 Eigen::Matrix3d mat;
1327 mat.col(0) << v2 - v1;
1328 mat.col(1) << v3 - v1;
1329 mat.col(2) << v4 - v1;
1330 return mat.determinant();
1324 double AdjointTools::tet_determinant(
const Eigen::VectorXd &v1,
const Eigen::VectorXd &v2,
const Eigen::VectorXd &v3,
const Eigen::VectorXd &v4) {
…}
1333 void AdjointTools::scaled_jacobian(
const Eigen::MatrixXd &
V,
const Eigen::MatrixXi &F, Eigen::VectorXd &quality)
1335 const int dim = F.cols() - 1;
1337 quality.setZero(F.rows());
1340 for (
int i = 0; i < F.rows(); i++)
1342 Eigen::RowVector3d e0;
1344 e0.head(2) =
V.row(F(i, 2)) -
V.row(F(i, 1));
1345 Eigen::RowVector3d e1;
1347 e1.head(2) =
V.row(F(i, 0)) -
V.row(F(i, 2));
1348 Eigen::RowVector3d e2;
1350 e2.head(2) =
V.row(F(i, 1)) -
V.row(F(i, 0));
1352 double l0 = e0.norm();
1353 double l1 = e1.norm();
1354 double l2 = e2.norm();
1356 double A = 0.5 * (e0.cross(e1)).norm();
1357 double Lmax = std::max(l0 * l1, std::max(l1 * l2, l0 * l2));
1359 quality(i) = 2 * A * (2 / sqrt(3)) / Lmax;
1364 for (
int i = 0; i < F.rows(); i++)
1366 Eigen::RowVector3d e0 =
V.row(F(i, 1)) -
V.row(F(i, 0));
1367 Eigen::RowVector3d e1 =
V.row(F(i, 2)) -
V.row(F(i, 1));
1368 Eigen::RowVector3d e2 =
V.row(F(i, 0)) -
V.row(F(i, 2));
1369 Eigen::RowVector3d e3 =
V.row(F(i, 3)) -
V.row(F(i, 0));
1370 Eigen::RowVector3d e4 =
V.row(F(i, 3)) -
V.row(F(i, 1));
1371 Eigen::RowVector3d e5 =
V.row(F(i, 3)) -
V.row(F(i, 2));
1373 double l0 = e0.norm();
1374 double l1 = e1.norm();
1375 double l2 = e2.norm();
1376 double l3 = e3.norm();
1377 double l4 = e4.norm();
1378 double l5 = e5.norm();
1380 double J = std::abs((e0.cross(e3)).dot(e2));
1382 double a1 = l0 * l2 * l3;
1383 double a2 = l0 * l1 * l4;
1384 double a3 = l1 * l2 * l5;
1385 double a4 = l3 * l4 * l5;
1387 double a = std::max({a1, a2, a3, a4, J});
1388 quality(i) = J * sqrt(2) / a;
1333 void AdjointTools::scaled_jacobian(
const Eigen::MatrixXd &
V,
const Eigen::MatrixXi &F, Eigen::VectorXd &quality) {
…}
1393 bool AdjointTools::is_flipped(
const Eigen::MatrixXd &
V,
const Eigen::MatrixXi &F)
1397 for (
int i = 0; i < F.rows(); i++)
1401 else if (F.cols() == 4)
1403 for (
int i = 0; i < F.rows(); i++)
1393 bool AdjointTools::is_flipped(
const Eigen::MatrixXd &
V,
const Eigen::MatrixXi &F) {
…}
ElementAssemblyValues vals
assembler::ElementAssemblyValues gvals
bool depend_on_gradu_local() const
void dj_du(const Eigen::MatrixXd &elastic_params, const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &u, const Eigen::MatrixXd &grad_u, const Eigen::MatrixXd &reference_normals, const assembler::ElementAssemblyValues &vals, ParameterType ¶ms, Eigen::MatrixXd &val) const
bool depend_on_gradu() const
void dj_dgradu(const Eigen::MatrixXd &elastic_params, const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &u, const Eigen::MatrixXd &grad_u, const Eigen::MatrixXd &reference_normals, const assembler::ElementAssemblyValues &vals, ParameterType ¶ms, Eigen::MatrixXd &val) const
void evaluate(const Eigen::MatrixXd &elastic_params, const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &u, const Eigen::MatrixXd &grad_u, const Eigen::MatrixXd &reference_normals, const assembler::ElementAssemblyValues &vals, ParameterType ¶ms, Eigen::MatrixXd &val) const
void dj_dx(const Eigen::MatrixXd &elastic_params, const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &u, const Eigen::MatrixXd &grad_u, const Eigen::MatrixXd &reference_normals, const assembler::ElementAssemblyValues &vals, ParameterType ¶ms, Eigen::MatrixXd &val) const
main class that contains the polyfem solver and all its state
Eigen::MatrixXd initial_vel_update
int n_bases
number of bases
assembler::AssemblyValsCache ass_vals_cache
used to store assembly values for small problems
const std::vector< basis::ElementBases > & geom_bases() const
Get a constant reference to the geometry mapping bases.
std::shared_ptr< assembler::Assembler > assembler
assemblers
ipc::CollisionMesh collision_mesh
IPC collision mesh.
std::shared_ptr< assembler::Mass > mass_matrix_assembler
std::vector< int > primitive_to_node() const
std::unique_ptr< mesh::Mesh > mesh
current mesh, it can be a Mesh2D or Mesh3D
std::shared_ptr< polyfem::mesh::MeshNodes > mesh_nodes
Mapping from input nodes to FE nodes.
std::shared_ptr< assembler::Problem > problem
current problem, it contains rhs and bc
std::vector< int > node_to_primitive() const
json args
main input arguments containing all defaults
solver::DiffCache diff_cached
void initial_velocity(Eigen::MatrixXd &velocity) const
Load or compute the initial velocity.
void initial_acceleration(Eigen::MatrixXd &acceleration) const
Load or compute the initial acceleration.
std::vector< basis::ElementBases > bases
FE bases, the size is #elements.
int n_geom_bases
number of geometric bases
std::vector< mesh::LocalBoundary > total_local_boundary
mapping from elements to nodes for all mesh
assembler::AssemblyValsCache mass_ass_vals_cache
std::vector< int > boundary_nodes
list of boundary nodes
solver::SolveData solve_data
timedependent stuff cached
bool is_contact_enabled() const
does the simulation has contact
StiffnessMatrix basis_nodes_to_gbasis_nodes
void compute(const int el_index, const bool is_volume, const basis::ElementBases &basis, const basis::ElementBases &gbasis, ElementAssemblyValues &vals) const
retrieves cached basis evaluation and geometric for the given element if it doesn't exist,...
stores per local bases evaluations
std::vector< basis::Local2Global > global
stores per element basis values at given quadrature points and geometric mapping
void compute(const int el_index, const bool is_volume, const Eigen::MatrixXd &pts, const basis::ElementBases &basis, const basis::ElementBases &gbasis)
computes the per element values at the local (ref el) points (pts) sets basis_values,...
quadrature::Quadrature quadrature
static void interpolate_at_local_vals(const mesh::Mesh &mesh, const bool is_problem_scalar, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const int el_index, const Eigen::MatrixXd &local_pts, const Eigen::MatrixXd &fun, Eigen::MatrixXd &result, Eigen::MatrixXd &result_grad)
interpolate solution and gradient at element (calls interpolate_at_local_vals with sol)
const StiffnessMatrix & gradu_h(int step) const
Eigen::MatrixXd disp_grad(int step=0) const
const ipc::Collisions & collision_set(int step) const
Eigen::VectorXd v(int step) const
Eigen::VectorXd u(int step) const
const ipc::FrictionCollisions & friction_collision_set(int step) const
Eigen::VectorXd apply_jacobian(const Eigen::VectorXd &grad, const Eigen::VectorXd &x) const override
std::shared_ptr< solver::FrictionForm > friction_form
std::shared_ptr< solver::InertiaForm > inertia_form
std::shared_ptr< solver::PeriodicContactForm > periodic_contact_form
std::shared_ptr< solver::PressureForm > pressure_form
std::shared_ptr< solver::BodyForm > body_form
std::shared_ptr< solver::NLProblem > nl_problem
std::shared_ptr< solver::ContactForm > contact_form
std::shared_ptr< solver::ElasticForm > damping_form
std::shared_ptr< solver::ElasticForm > elastic_form
Eigen::Matrix< double, dim, 1 > cross(const Eigen::Matrix< double, dim, 1 > &x, const Eigen::Matrix< double, dim, 1 > &y)
DScalar1< double, Eigen::Matrix< double, Eigen::Dynamic, 1 > > Diff
void vector2matrix(const Eigen::VectorXd &vec, Eigen::MatrixXd &mat)
auto & get_local_thread_storage(Storages &storage, int thread_id)
auto create_thread_storage(const LocalStorage &initial_local_storage)
double triangle_area(const Eigen::MatrixXd V)
Compute the signed area of a triangle defined by three points.
void maybe_parallel_for(int size, const std::function< void(int, int, int)> &partial_for)
Eigen::Matrix< double, Eigen::Dynamic, 1, 0, MAX_QUAD_POINTS, 1 > QuadratureVector
void log_and_throw_adjoint_error(const std::string &msg)
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix
Automatic differentiation scalar with first-order derivatives.
static void setVariableCount(size_t value)
Set the independent variable count used by the automatic differentiation layer.
Parameters for the functional evaluation.