19 using DefGradMatrix = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3>;
21 template <
typename Derived>
41 const std::function<Eigen::MatrixXd(
const Eigen::MatrixXd &)> &fun)
const override;
44 const Eigen::MatrixXd &mat,
45 Eigen::MatrixXd &stress,
46 Eigen::MatrixXd &result)
const override;
49 Eigen::MatrixXd &stress,
50 Eigen::MatrixXd &result)
const override;
53 const Eigen::MatrixXd &vect,
54 Eigen::MatrixXd &stress,
55 Eigen::MatrixXd &result)
const override;
58 Derived &
derived() {
return static_cast<Derived &
>(*this); }
60 const Derived &
derived()
const {
return static_cast<const Derived &
>(*this); }
72 virtual Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3>
gradient(
76 const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> &
F)
const
81 virtual Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 9, 9>
hessian(
85 const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> &
F)
const
95 typedef Eigen::Matrix<T, Eigen::Dynamic, 1>
AutoDiffVect;
96 typedef Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3>
AutoDiffGradMat;
97 typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> DoubleGradMat;
106 const int n_pts = data.
da.size();
107 for (
long p = 0; p < n_pts; ++p)
112 for (
int d = 0; d <
size(); ++d)
113 def_grad(d, d) += T(1);
117 DoubleGradMat tmp_jac_it = data.
vals.
jac_it[p];
118 tmp_jac_it = tmp_jac_it.inverse();
121 for (
long k = 0; k < jac_it.size(); ++k)
122 jac_it(k) = T(tmp_jac_it(k));
129 energy +=
val * data.
da(p);
143 template <
int n_basis,
int dim>
147 typedef Eigen::Matrix<Diff, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3>
AutoDiffGradMat;
148 typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> DoubleGradMat;
150 const int n_pts = data.
da.size();
153 local_disp.setZero();
157 for (
size_t ii = 0; ii < bs.global.size(); ++ii)
159 for (
int d = 0; d < dim; ++d)
161 local_disp(i, d) += bs.global[ii].val * data.
x(bs.global[ii].index *
size() + d);
166 Eigen::Matrix<double, dim, dim> def_grad(
size(),
size());
173 for (
long p = 0; p < n_pts; ++p)
180 Eigen::Matrix<double, dim, dim> jac_it = data.
vals.
jac_it[p];
181 Eigen::Matrix<double, n_basis, dim> G = grad * jac_it;
182 def_grad = local_disp.transpose() * G + Eigen::Matrix<double, dim, dim>::Identity(
size(),
size());
184 for (
int d1 = 0; d1 <
size(); ++d1)
185 for (
int d2 = 0; d2 <
size(); ++d2)
186 def_grad_ad(d1, d2) =
Diff(d1 *
size() + d2, def_grad(d1, d2));
190 DoubleGradMat tmp_jac_it = data.
vals.
jac_it[p];
191 tmp_jac_it = tmp_jac_it.inverse();
194 for (
long k = 0; k < jac_it.size(); ++k)
195 jac_it(k) =
Diff(tmp_jac_it(k));
197 def_grad_ad *= jac_it;
202 Eigen::Matrix<double, dim, dim> P;
203 for (
int i = 0; i <
size(); ++i)
204 for (
int j = 0; j <
size(); ++j)
205 P(i, j) =
val.getGradient()(i *
size() + j);
207 const Eigen::Matrix<double, n_basis, dim> Rloc = G * P.transpose() * data.
da(p);
210 for (
int d = 0; d <
size(); ++d)
211 res(a *
size() + d) += Rloc(a, d);
215 template <
int n_basis,
int dim>
218 typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> DoubleGradMat;
220 const int n_pts = data.
da.size();
223 local_disp.setZero();
227 for (
size_t ii = 0; ii < bs.global.size(); ++ii)
229 for (
int d = 0; d < dim; ++d)
231 local_disp(i, d) += bs.global[ii].val * data.
x(bs.global[ii].index *
size() + d);
236 Eigen::Matrix<double, dim, dim> def_grad(
size(),
size());
240 for (
long p = 0; p < n_pts; ++p)
247 Eigen::Matrix<double, dim, dim> jac_it = data.
vals.
jac_it[p];
248 Eigen::Matrix<double, n_basis, dim> G = grad * jac_it;
249 def_grad = local_disp.transpose() * G + Eigen::Matrix<double, dim, dim>::Identity(
size(),
size());
254 jac_it = jac_it.inverse();
260 const Eigen::Matrix<double, n_basis, dim> Bgrad =
derived().real_def_grad() ? G : grad;
261 const Eigen::Matrix<double, n_basis, dim> Rloc = Bgrad * P.transpose() * data.
da(p);
264 for (
int d = 0; d <
size(); ++d)
265 res(a *
size() + d) += Rloc(a, d);
269 template <
int n_basis,
int dim>
272 typedef DScalar2<double, Eigen::Matrix<double, Eigen::Dynamic, 1, 0, 9, 1>, Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 9, 9>> Diff2;
273 typedef Eigen::Matrix<Diff2, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3>
AutoDiffGradMat;
274 typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> DoubleGradMat;
276 const int d =
size();
278 const int n_pts = data.
da.size();
280 Eigen::Matrix<double, n_basis, dim> local_disp(nb, d);
281 local_disp.setZero();
285 for (
size_t ii = 0; ii < bs.global.size(); ++ii)
287 for (
int dd = 0; dd < d; ++dd)
288 local_disp(i, dd) += bs.global[ii].val * data.
x(bs.global[ii].index * d + dd);
294 Eigen::Matrix<double, dim, dim> def_grad(d, d);
296 H.resize(nb * d, nb * d);
299 for (
long p = 0; p < n_pts; ++p)
301 Eigen::Matrix<double, n_basis, dim> grad_ref(nb, d);
305 const Eigen::Matrix<double, dim, dim> jac_it = data.
vals.
jac_it[p];
306 const Eigen::Matrix<double, n_basis, dim> G = grad_ref * jac_it;
308 def_grad = local_disp.transpose() * G + Eigen::Matrix<double, dim, dim>::Identity();
310 for (
int i = 0; i < d; ++i)
311 for (
int j = 0; j < d; ++j)
312 def_grad_ad(i, j) = Diff2(i * d + j, def_grad(i, j));
316 DoubleGradMat tmp_jac_it = data.
vals.
jac_it[p];
317 tmp_jac_it = tmp_jac_it.inverse();
320 for (
long k = 0; k < jac_it.size(); ++k)
321 jac_it(k) = Diff2(tmp_jac_it(k));
323 def_grad_ad *= jac_it;
328 Eigen::Matrix<double, dim * dim, dim * dim> A(d * d, d * d);
333 for (
int r = 0; r < d * d; ++r)
334 for (
int c = 0; c < d * d; ++c)
335 A(r, c) =
val.getHessian()(r, c);
337 for (
int a = 0; a < nb; ++a)
339 const Eigen::Matrix<double, dim * dim, dim> Ba = compute_B_block<dim>(G.row(a));
341 for (
int b = 0; b < nb; ++b)
343 const Eigen::Matrix<double, dim * dim, dim> Bb = compute_B_block<dim>(G.row(b));
345 const Eigen::Matrix<double, dim, dim> Kab =
346 Ba.transpose() * A * Bb * data.
da(p);
348 H.template block<dim, dim>(a * d, b * d) += Kab;
354 template <
int n_basis,
int dim>
357 typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> DoubleGradMat;
359 const int d =
size();
361 const int n_pts = data.
da.size();
363 Eigen::Matrix<double, n_basis, dim> local_disp(nb, d);
364 local_disp.setZero();
368 for (
size_t ii = 0; ii < bs.global.size(); ++ii)
370 for (
int dd = 0; dd < d; ++dd)
371 local_disp(i, dd) += bs.global[ii].val * data.
x(bs.global[ii].index * d + dd);
375 Eigen::Matrix<double, dim, dim> def_grad(d, d);
377 H.resize(nb * d, nb * d);
380 for (
long p = 0; p < n_pts; ++p)
382 Eigen::Matrix<double, n_basis, dim> grad_ref(nb, d);
386 const Eigen::Matrix<double, dim, dim> jac_it = data.
vals.
jac_it[p];
387 const Eigen::Matrix<double, n_basis, dim> G = grad_ref * jac_it;
389 def_grad = local_disp.transpose() * G + Eigen::Matrix<double, dim, dim>::Identity();
394 jac_it = jac_it.inverse();
403 const Eigen::Matrix<double, n_basis, dim> Bgrad =
derived().real_def_grad() ? G : grad_ref;
405 for (
int a = 0; a < nb; ++a)
407 const Eigen::Matrix<double, dim * dim, dim> Ba = compute_B_block<dim>(Bgrad.row(a));
409 for (
int b = 0; b < nb; ++b)
411 const Eigen::Matrix<double, dim * dim, dim> Bb = compute_B_block<dim>(Bgrad.row(b));
413 const Eigen::Matrix<double, dim, dim> Kab =
414 Ba.transpose() * A * Bb * data.
da(p);
416 H.template block<dim, dim>(a * d, b * d) += Kab;
423 Eigen::Matrix<double, dim * dim, dim>
compute_B_block(
const Eigen::Matrix<double, 1, dim> &g)
const;
virtual void assemble_gradient(const bool is_volume, const int n_basis, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const AssemblyValsCache &cache, const double t, const double dt, const Eigen::MatrixXd &displacement, const Eigen::MatrixXd &displacement_prev, Eigen::MatrixXd &rhs) const
virtual void assemble_hessian(const bool is_volume, const int n_basis, const bool project_to_psd, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const AssemblyValsCache &cache, const double t, const double dt, const Eigen::MatrixXd &displacement, const Eigen::MatrixXd &displacement_prev, utils::MatrixCache &mat_cache, StiffnessMatrix &grad) const
virtual double assemble_energy(const bool is_volume, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const AssemblyValsCache &cache, const double t, const double dt, const Eigen::MatrixXd &displacement, const Eigen::MatrixXd &displacement_prev) const
std::vector< AssemblyValues > basis_values
std::vector< Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3 > > jac_it
void compute_stress_grad_multiply_stress(const OptAssemblerData &data, Eigen::MatrixXd &stress, Eigen::MatrixXd &result) const override
void compute_hessian_from_stress_noad(const NonLinearAssemblerData &data, Eigen::MatrixXd &H) 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
Eigen::Matrix< double, dim *dim, dim > compute_B_block(const Eigen::Matrix< double, 1, dim > &g) const
Eigen::VectorXd assemble_gradient_full_ad(const NonLinearAssemblerData &data) const
virtual Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, 0, 9, 9 > hessian(const RowVectorNd &p, const double t, const int el_id, const Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3 > &F) const
Derived & derived()
Returns this as a reference to derived class.
Eigen::MatrixXd assemble_hessian_full_ad(const NonLinearAssemblerData &data) const
const Derived & derived() const
Returns this as a const reference to derived class.
Eigen::MatrixXd assemble_hessian_stress_noad(const NonLinearAssemblerData &data) const
virtual ~GenericElastic()=default
Eigen::VectorXd assemble_gradient_stress_noad(const NonLinearAssemblerData &data) const
Eigen::MatrixXd assemble_hessian_stress_ad(const NonLinearAssemblerData &data) const
Eigen::MatrixXd assemble_hessian(const NonLinearAssemblerData &data) const override
virtual void add_multimaterial(const int index, const json ¶ms, const Units &units) override=0
void compute_stress_grad_multiply_vect(const OptAssemblerData &data, const Eigen::MatrixXd &vect, Eigen::MatrixXd &stress, Eigen::MatrixXd &result) const override
Eigen::VectorXd assemble_gradient(const NonLinearAssemblerData &data) const override
bool allow_inversion() const override
Eigen::VectorXd assemble_gradient_stress_ad(const NonLinearAssemblerData &data) const
void compute_hessian_from_stress(const NonLinearAssemblerData &data, Eigen::MatrixXd &H) const
void compute_gradient_from_stress_noad(const NonLinearAssemblerData &data, Eigen::VectorXd &res) const
virtual Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3 > gradient(const RowVectorNd &p, const double t, const int el_id, const Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3 > &F) const
AutodiffType autodiff_type_
void compute_gradient_from_stress(const NonLinearAssemblerData &data, Eigen::VectorXd &res) const
double compute_energy(const NonLinearAssemblerData &data) const override
T compute_energy_aux(const NonLinearAssemblerData &data) const
void compute_stress_grad_multiply_mat(const OptAssemblerData &data, const Eigen::MatrixXd &mat, Eigen::MatrixXd &stress, Eigen::MatrixXd &result) const override
virtual bool real_def_grad() const
const ElementAssemblyValues & vals
const Eigen::MatrixXd & x
const QuadratureVector & da
Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3 > DefGradMatrix
Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3 > AutoDiffGradMat
Eigen::Matrix< T, Eigen::Dynamic, 1 > AutoDiffVect
DScalar1< double, Eigen::Matrix< double, Eigen::Dynamic, 1 > > Diff
void compute_disp_grad_at_quad(const assembler::NonLinearAssemblerData &data, const AutoDiffVect &local_disp, const int p, const int size, AutoDiffGradMat &def_grad)
void get_local_disp(const assembler::NonLinearAssemblerData &data, const int size, AutoDiffVect &local_disp)
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > RowVectorNd
void log_and_throw_error(const std::string &msg)
Automatic differentiation scalar with first-order derivatives.
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.