10 bool delta(
int i,
int j)
12 return (i == j) ? true :
false;
15 template <
class T,
int p = 3>
18 constexpr static double C = 1e6;
21 static_assert(p % 2 == 1);
26 const T tmp1 = 2 * J - 1;
27 const T tmp2 = pow(tmp1, p);
28 return C * (1 / (tmp2 + 1) - 1);
31 static T first_derivatives(T J)
35 const T tmp1 = 2 * J - 1;
36 const T tmp2 = pow(tmp1, p);
37 return C * -2 * p * tmp2 / tmp1 / pow(1 + tmp2, 2);
40 static T second_derivatives(T J)
44 const T tmp1 = 2 * J - 1;
45 const T tmp2 = pow(tmp1, p);
46 return C * 4 * p * tmp2 / pow(tmp1, 2) * ((1 - p) + (1 + p) * tmp2) / pow(1 + tmp2, 3);
51 Eigen::Matrix<double, dim, dim>
hat(
const Eigen::Matrix<double, dim, 1> &
x)
54 Eigen::Matrix<double, dim, dim> prod;
68 Eigen::Matrix<double, dim, 1>
cross(
const Eigen::Matrix<double, dim, 1> &
x,
const Eigen::Matrix<double, dim, 1> &
y)
71 Eigen::Matrix<double, dim, 1>
z;
74 z(0) =
x(1) *
y(2) -
x(2) *
y(1);
75 z(1) =
x(2) *
y(0) -
x(0) *
y(2);
76 z(2) =
x(0) *
y(1) -
x(1) *
y(0);
81 template <
int dimt,
class T>
82 Eigen::Matrix<T, dimt, dimt> get_standard(
const int dim,
const bool use_rest_pose)
84 Eigen::Matrix<double, dimt, dimt> standard(dim, dim);
87 standard.setIdentity();
93 0.5, std::sqrt(3) / 2;
96 0.5, std::sqrt(3) / 2., 0,
97 0.5, 0.5 / std::sqrt(3), std::sqrt(3) / 2.;
98 standard = standard.inverse().transpose().eval();
101 Eigen::Matrix<T, dimt, dimt> res(dim, dim);
102 for (
int i = 0; i < dim; ++i)
104 for (
int j = 0; j < dim; ++j)
106 res(i, j) = T(standard(i, j));
115 return compute_energy_aux<double>(data);
122 size(), n_bases, data,
123 [&](
const NonLinearAssemblerData &data) {
return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 6, 1>>>(data); },
124 [&](
const NonLinearAssemblerData &data) {
return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 8, 1>>>(data); },
125 [&](
const NonLinearAssemblerData &data) {
return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 12, 1>>>(data); },
126 [&](
const NonLinearAssemblerData &data) {
return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 18, 1>>>(data); },
127 [&](
const NonLinearAssemblerData &data) {
return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 24, 1>>>(data); },
128 [&](
const NonLinearAssemblerData &data) {
return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 30, 1>>>(data); },
129 [&](
const NonLinearAssemblerData &data) {
return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 60, 1>>>(data); },
130 [&](
const NonLinearAssemblerData &data) {
return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 81, 1>>>(data); },
131 [&](
const NonLinearAssemblerData &data) {
return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, Eigen::Dynamic, 1, 0, SMALL_N, 1>>>(data); },
132 [&](
const NonLinearAssemblerData &data) {
return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, Eigen::Dynamic, 1, 0, BIG_N, 1>>>(data); },
133 [&](
const NonLinearAssemblerData &data) {
return compute_energy_aux<DScalar1<double, Eigen::VectorXd>>(data); });
199 template <
typename T>
202 typedef Eigen::Matrix<T, Eigen::Dynamic, 1>
AutoDiffVect;
203 typedef Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3>
AutoDiffGradMat;
204 typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> DoubleGradMat;
208 power =
size() == 2 ? 1. : (2. / 3.);
210 power =
size() == 2 ? 2. : 5. / 3.;
226 const int n_pts = data.
da.size();
227 for (
long p = 0; p < n_pts; ++p)
232 for (
int d = 0; d <
size(); ++d)
233 def_grad(d, d) += T(1);
237 DoubleGradMat tmp_jac_it = data.
vals.
jac_it[p];
238 tmp_jac_it = tmp_jac_it.inverse();
241 for (
long k = 0; k < jac_it.size(); ++k)
242 jac_it(k) = T(tmp_jac_it(k));
245 def_grad = def_grad * standard;
251 energy = std::nan(
"");
255 const T powJ = pow(det, power);
256 const T
val = (def_grad.transpose() * def_grad).trace() / powJ;
258 energy +=
val * data.
da(p);
263 template <
int n_basis,
int dim>
266 assert(data.
x.cols() == 1);
268 const int n_pts = data.
da.size();
271 local_disp.setZero();
275 for (
size_t ii = 0; ii < bs.global.size(); ++ii)
277 for (
int d = 0; d <
size(); ++d)
279 local_disp(i, d) += bs.global[ii].val * data.
x(bs.global[ii].index *
size() + d);
284 Eigen::Matrix<double, dim, dim> def_grad(
size(),
size());
289 const Eigen::Matrix<double, dim, dim> standard = get_standard<dim, double>(
size(),
use_rest_pose_);
291 for (
long p = 0; p < n_pts; ++p)
300 Eigen::Matrix<double, dim, dim> jac_it = data.
vals.
jac_it[p];
303 def_grad = (local_disp.transpose() * grad + jac_it.inverse()) * standard;
305 double J = def_grad.determinant();
309 Eigen::Matrix<double, dim, dim> delJ_delF(
size(),
size());
315 delJ_delF(0, 0) = def_grad(1, 1);
316 delJ_delF(0, 1) = -def_grad(1, 0);
317 delJ_delF(1, 0) = -def_grad(0, 1);
318 delJ_delF(1, 1) = def_grad(0, 0);
323 Eigen::Matrix<double, dim, 1> u(def_grad.rows());
324 Eigen::Matrix<double, dim, 1> v(def_grad.rows());
325 Eigen::Matrix<double, dim, 1> w(def_grad.rows());
331 delJ_delF.col(0) = cross<dim>(v, w);
332 delJ_delF.col(1) = cross<dim>(w, u);
333 delJ_delF.col(2) = cross<dim>(u, v);
336 const double m = (dim == 2) ? 2. : 5. / 3.;
337 const double powJ = pow(J, m);
338 Eigen::Matrix<double, dim, dim> gradient_temp = (2 * def_grad - (def_grad.squaredNorm() * m) / J * delJ_delF) / powJ;
339 Eigen::Matrix<double, n_basis, dim> gradient = grad * standard * gradient_temp.transpose();
341 G.noalias() += gradient * data.
da(p);
344 Eigen::Matrix<double, dim, n_basis> G_T = G.transpose();
346 constexpr int N = (n_basis == Eigen::Dynamic) ? Eigen::Dynamic : n_basis * dim;
347 Eigen::Matrix<double, N, 1> temp(Eigen::Map<Eigen::Matrix<double, N, 1>>(G_T.data(), G_T.size()));
355 size(), n_bases, data,
356 [&](
const NonLinearAssemblerData &data) {
return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 6, 1>, Eigen::Matrix<double, 6, 6>>>(data); },
357 [&](
const NonLinearAssemblerData &data) {
return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 8, 1>, Eigen::Matrix<double, 8, 8>>>(data); },
358 [&](
const NonLinearAssemblerData &data) {
return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 12, 1>, Eigen::Matrix<double, 12, 12>>>(data); },
359 [&](
const NonLinearAssemblerData &data) {
return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 18, 1>, Eigen::Matrix<double, 18, 18>>>(data); },
360 [&](
const NonLinearAssemblerData &data) {
return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 24, 1>, Eigen::Matrix<double, 24, 24>>>(data); },
361 [&](
const NonLinearAssemblerData &data) {
return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 30, 1>, Eigen::Matrix<double, 30, 30>>>(data); },
362 [&](
const NonLinearAssemblerData &data) {
return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 60, 1>, Eigen::Matrix<double, 60, 60>>>(data); },
363 [&](
const NonLinearAssemblerData &data) {
return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 81, 1>, Eigen::Matrix<double, 81, 81>>>(data); },
364 [&](
const NonLinearAssemblerData &data) {
return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, Eigen::Dynamic, 1, 0, SMALL_N, 1>, Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, SMALL_N, SMALL_N>>>(data); },
365 [&](
const NonLinearAssemblerData &data) {
return compute_energy_aux<DScalar2<double, Eigen::VectorXd, Eigen::MatrixXd>>(data); });
442 template <
int n_basis,
int dim>
445 assert(data.
x.cols() == 1);
447 constexpr int N = (n_basis == Eigen::Dynamic) ? Eigen::Dynamic : n_basis * dim;
448 const int n_pts = data.
da.size();
451 local_disp.setZero();
455 for (
size_t ii = 0; ii < bs.global.size(); ++ii)
457 for (
int d = 0; d <
size(); ++d)
459 local_disp(i, d) += bs.global[ii].val * data.
x(bs.global[ii].index *
size() + d);
464 Eigen::Matrix<double, dim, dim> def_grad(
size(),
size());
466 const Eigen::Matrix<double, dim, dim> standard = get_standard<dim, double>(
size(),
use_rest_pose_);
468 for (
long p = 0; p < n_pts; ++p)
477 Eigen::Matrix<double, dim, dim> jac_it = data.
vals.
jac_it[p];
482 def_grad = local_disp.transpose() * grad + Eigen::Matrix<double, dim, dim>::Identity(
size(),
size());
486 def_grad = (local_disp.transpose() * grad + jac_it.inverse()) * standard;
489 Eigen::Matrix<double, dim * dim, dim * dim> hessian_temp;
494 Eigen::Matrix<Diff, dim, dim> def_grad_ad(dim, dim);
495 for (
int i = 0; i < def_grad.size(); i++)
496 def_grad_ad(i) =
Diff(i, def_grad(i));
500 hessian_temp =
val.getHessian();
555 Eigen::Matrix<double, dim * dim, N> delF_delU_tensor(jac_it.size(), grad.size());
557 for (
size_t i = 0; i < local_disp.rows(); ++i)
559 for (
size_t j = 0; j < local_disp.cols(); ++j)
561 Eigen::Matrix<double, dim, dim> temp(
size(),
size());
563 temp.row(j) = grad.row(i);
564 temp = temp * standard;
565 Eigen::Matrix<double, dim * dim, 1> temp_flattened(Eigen::Map<Eigen::Matrix<double, dim * dim, 1>>(temp.data(), temp.size()));
566 delF_delU_tensor.col(i *
size() + j) = temp_flattened;
570 Eigen::Matrix<double, N, N> hessian = delF_delU_tensor.transpose() * hessian_temp * delF_delU_tensor;
572 H += hessian * data.
da(p);
579 Eigen::MatrixXd &all,
580 const std::function<Eigen::MatrixXd(
const Eigen::MatrixXd &)> &fun)
const
582 const auto &displacement = data.
fun;
584 const auto &bs = data.
bs;
585 const auto &gbs = data.
gbs;
586 const auto el_id = data.
el_id;
587 const auto t = data.
t;
589 Eigen::MatrixXd displacement_grad(
size(),
size());
591 assert(displacement.cols() == 1);
593 all.setZero(local_pts.rows(), all_size);
600 if (params.contains(
"use_rest_pose"))
613 return std::map<std::string, ParamFunc>();
618 if (params.contains(
"canonical_transformation"))
621 for (
int i = 0; i < params[
"canonical_transformation"].size(); ++i)
623 Eigen::MatrixXd transform_matrix(
size(),
size());
624 for (
int j = 0; j <
size(); ++j)
625 for (
int k = 0; k <
size(); ++k)
626 transform_matrix(j, k) = params[
"canonical_transformation"][i][j][k];
std::map< std::string, ParamFunc > parameters() const override
std::vector< Eigen::MatrixXd > canonical_transformation_
void add_multimaterial(const int index, const json ¶ms, const Units &units) override
double compute_energy(const NonLinearAssemblerData &data) const override
T compute_energy_aux(const NonLinearAssemblerData &data) const
void assign_stress_tensor(const OutputData &data, const int all_size, const ElasticityTensorType &type, Eigen::MatrixXd &all, const std::function< Eigen::MatrixXd(const Eigen::MatrixXd &)> &fun) const override
void compute_energy_hessian_aux_fast(const NonLinearAssemblerData &data, Eigen::MatrixXd &H) const
void add_multimaterial(const int index, const json ¶ms, const Units &units) override
Eigen::MatrixXd assemble_hessian(const NonLinearAssemblerData &data) const override
void compute_energy_aux_gradient_fast(const NonLinearAssemblerData &data, Eigen::VectorXd &G_flattened) const
Eigen::VectorXd assemble_gradient(const NonLinearAssemblerData &data) const override
std::vector< AssemblyValues > basis_values
std::vector< Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3 > > jac_it
const ElementAssemblyValues & vals
const Eigen::MatrixXd & x
const QuadratureVector & da
const basis::ElementBases & bs
const Eigen::MatrixXd & fun
const Eigen::MatrixXd & local_pts
const basis::ElementBases & gbs
Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3 > AutoDiffGradMat
Eigen::Matrix< double, dim, dim > hat(const Eigen::Matrix< double, dim, 1 > &x)
Eigen::Matrix< double, dim, 1 > cross(const Eigen::Matrix< double, dim, 1 > &x, const Eigen::Matrix< double, dim, 1 > &y)
Eigen::Matrix< T, Eigen::Dynamic, 1 > AutoDiffVect
DScalar1< double, Eigen::Matrix< double, Eigen::Dynamic, 1 > > Diff
T determinant(const Eigen::Matrix< T, rows, cols, option, maxRow, maxCol > &mat)
void compute_disp_grad_at_quad(const assembler::NonLinearAssemblerData &data, const AutoDiffVect &local_disp, const int p, const int size, AutoDiffGradMat &def_grad)
Eigen::MatrixXd hessian_from_energy(const int size, const int n_bases, const assembler::NonLinearAssemblerData &data, const std::function< DScalar2< double, Eigen::Matrix< double, 6, 1 >, Eigen::Matrix< double, 6, 6 > >(const assembler::NonLinearAssemblerData &)> &fun6, const std::function< DScalar2< double, Eigen::Matrix< double, 8, 1 >, Eigen::Matrix< double, 8, 8 > >(const assembler::NonLinearAssemblerData &)> &fun8, const std::function< DScalar2< double, Eigen::Matrix< double, 12, 1 >, Eigen::Matrix< double, 12, 12 > >(const assembler::NonLinearAssemblerData &)> &fun12, const std::function< DScalar2< double, Eigen::Matrix< double, 18, 1 >, Eigen::Matrix< double, 18, 18 > >(const assembler::NonLinearAssemblerData &)> &fun18, const std::function< DScalar2< double, Eigen::Matrix< double, 24, 1 >, Eigen::Matrix< double, 24, 24 > >(const assembler::NonLinearAssemblerData &)> &fun24, const std::function< DScalar2< double, Eigen::Matrix< double, 30, 1 >, Eigen::Matrix< double, 30, 30 > >(const assembler::NonLinearAssemblerData &)> &fun30, const std::function< DScalar2< double, Eigen::Matrix< double, 60, 1 >, Eigen::Matrix< double, 60, 60 > >(const assembler::NonLinearAssemblerData &)> &fun60, const std::function< DScalar2< double, Eigen::Matrix< double, 81, 1 >, Eigen::Matrix< double, 81, 81 > >(const assembler::NonLinearAssemblerData &)> &fun81, const std::function< DScalar2< double, Eigen::Matrix< double, Eigen::Dynamic, 1, 0, SMALL_N, 1 >, Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, 0, SMALL_N, SMALL_N > >(const assembler::NonLinearAssemblerData &)> &funN, const std::function< DScalar2< double, Eigen::VectorXd, Eigen::MatrixXd >(const assembler::NonLinearAssemblerData &)> &funn)
void get_local_disp(const assembler::NonLinearAssemblerData &data, const int size, AutoDiffVect &local_disp)
Eigen::VectorXd gradient_from_energy(const int size, const int n_bases, const assembler::NonLinearAssemblerData &data, const std::function< DScalar1< double, Eigen::Matrix< double, 6, 1 > >(const assembler::NonLinearAssemblerData &)> &fun6, const std::function< DScalar1< double, Eigen::Matrix< double, 8, 1 > >(const assembler::NonLinearAssemblerData &)> &fun8, const std::function< DScalar1< double, Eigen::Matrix< double, 12, 1 > >(const assembler::NonLinearAssemblerData &)> &fun12, const std::function< DScalar1< double, Eigen::Matrix< double, 18, 1 > >(const assembler::NonLinearAssemblerData &)> &fun18, const std::function< DScalar1< double, Eigen::Matrix< double, 24, 1 > >(const assembler::NonLinearAssemblerData &)> &fun24, const std::function< DScalar1< double, Eigen::Matrix< double, 30, 1 > >(const assembler::NonLinearAssemblerData &)> &fun30, const std::function< DScalar1< double, Eigen::Matrix< double, 60, 1 > >(const assembler::NonLinearAssemblerData &)> &fun60, const std::function< DScalar1< double, Eigen::Matrix< double, 81, 1 > >(const assembler::NonLinearAssemblerData &)> &fun81, const std::function< DScalar1< double, Eigen::Matrix< double, Eigen::Dynamic, 1, 0, SMALL_N, 1 > >(const assembler::NonLinearAssemblerData &)> &funN, const std::function< DScalar1< double, Eigen::Matrix< double, Eigen::Dynamic, 1, 0, BIG_N, 1 > >(const assembler::NonLinearAssemblerData &)> &funBigN, const std::function< DScalar1< double, Eigen::VectorXd >(const assembler::NonLinearAssemblerData &)> &funn)
Automatic differentiation scalar with first- and second-order derivatives.
static void setVariableCount(size_t value)
Set the independent variable count used by the automatic differentiation layer.