From d3a696ab0820e394cdd6f4f07e678ce99a5a0e0f Mon Sep 17 00:00:00 2001 From: Andreas Longva Date: Fri, 12 May 2023 17:41:59 +0200 Subject: [PATCH] Stopgap design for having Hyperelastic materials depend on displacement gradient By depending on nabla u instead of the deformation gradient F, we can avoid the loss of accuracy in forming F, which causes a severe loss of accuracy for very stiff materials subject to small displacements. --- fenris-solid/src/lib.rs | 262 +++++++++++++++++++++++- fenris-solid/src/logdet.rs | 86 ++++++++ fenris-solid/src/materials.rs | 31 +-- fenris-solid/tests/unit_tests/logdet.rs | 50 +++++ fenris-solid/tests/unit_tests/mod.rs | 1 + src/assembly/local/elliptic.rs | 4 +- src/assembly/operators.rs | 16 +- 7 files changed, 422 insertions(+), 28 deletions(-) create mode 100644 fenris-solid/src/logdet.rs create mode 100644 fenris-solid/tests/unit_tests/logdet.rs diff --git a/fenris-solid/src/lib.rs b/fenris-solid/src/lib.rs index 81d17211..832fae68 100644 --- a/fenris-solid/src/lib.rs +++ b/fenris-solid/src/lib.rs @@ -1,15 +1,61 @@ //! Solid mechanics functionality for `fenris`. use fenris::allocators::DimAllocator; use fenris::assembly::operators::{EllipticContraction, EllipticEnergy, EllipticOperator, Operator}; -use fenris::nalgebra::{DMatrixViewMut, DVectorView, DefaultAllocator, DimName, OMatrix, OVector}; +use fenris::nalgebra::{ + DMatrixViewMut, DVectorView, DefaultAllocator, DimName, OMatrix, OVector, RealField, U1, U2, U3, +}; use fenris::{Real, SmallDim, Symmetry}; use std::cmp::min; pub mod materials; +mod logdet; +pub use logdet::log_det_F; + mod gravity_source; pub use gravity_source::GravitySource; +/// Compute the deformation gradient $\vec F$ given the displacement gradient $\nabla \vec u$. +#[allow(non_snake_case)] +pub fn deformation_gradient(u_grad: &OMatrix) -> OMatrix +where + T: RealField, + D: DimName, + DefaultAllocator: DimAllocator, +{ + let I = OMatrix::::identity(); + let F = I + u_grad.transpose(); + F +} + +/// Compute the displacement gradient $\nabla \vec u$ given the deformation gradient $\vec F$. +#[allow(non_snake_case)] +pub fn u_grad_from_F(F: &OMatrix) -> OMatrix +where + T: RealField, + D: DimName, + DefaultAllocator: DimAllocator, +{ + let I = OMatrix::::identity(); + (F - I).transpose() +} + +/// A hyperelastic material defined by its energy density $\psi(\vec F)$. +/// +/// Although hyperelastic materials in literature are defined in terms of the deformation gradient +/// $\vec F$, computing $\vec F = \vec I + \nabla \vec u^T$ incurs a significant loss in accuracy +/// if $\nabla \vec u$ is very small (such as for stiff materials). Therefore, in addition +/// to the required methods that work with $\vec F$, such as +/// [`compute_energy_density`](Self::compute_energy_density), there are additional methods with the +/// `_du` suffix which refer to a variant of the method where the input is the displacement +/// gradient $\nabla \vec u$ instead of $\vec F$. These methods have default impls which rely +/// on computing $\vec F$, but can be overridden by individual materials for more accurate results. +/// +/// See the [libCEED documentation](https://libceed.org/en/latest/examples/solids) for a discussion +/// of the numerical problems associated with forming $\vec F$ and derived quantities such as +/// $\log \det \vec F$. +/// +/// In the future, the trait may be changed to always work with $\nabla \vec u$ instead of $\vec F$. pub trait HyperelasticMaterial where T: Real, @@ -25,6 +71,18 @@ where parameters: &Self::Parameters, ) -> T; + /// Compute the energy density $\psi = \psi(\vec F)$ associated with the material. + /// + /// Takes as input the displacement gradient $\nabla \vec u$ instead of the deformation + /// gradient $\vec F$. See [the trait docs](Self) for more information. + fn compute_energy_density_du( + &self, + u_grad: &OMatrix, + parameters: &Self::Parameters, + ) -> T { + self.compute_energy_density(&deformation_gradient(u_grad), parameters) + } + /// Compute the First Piola-Kirchhoff stress tensor $\vec P = \vec P(\vec F)$. fn compute_stress_tensor( &self, @@ -32,6 +90,19 @@ where parameters: &Self::Parameters, ) -> OMatrix; + /// Compute the First Piola-Kirchhoff stress tensor $\vec P = \vec P(\vec F)$. + /// + /// Takes as input the displacement gradient $\nabla \vec u$ instead of the deformation + /// gradient $\vec F$. See [the trait docs](Self) for more information. + #[allow(non_snake_case)] + fn compute_stress_tensor_du( + &self, + u_grad: &OMatrix, + parameters: &Self::Parameters, + ) -> OMatrix { + self.compute_stress_tensor(&deformation_gradient(u_grad), parameters) + } + /// Compute the stress contraction operator $\\mathcal{C}\_{\vec P}(\vec F, \vec a, \vec b)$ with the given /// material parameters. /// @@ -48,6 +119,21 @@ where parameters: &Self::Parameters, ) -> OMatrix; + /// Compute the stress contraction operator $\\mathcal{C}\_{\vec P}(\vec F, \vec a, \vec b)$ with the given + /// material parameters. + /// + /// Takes as input the displacement gradient $\nabla \vec u$ instead of the deformation + /// gradient $\vec F$. See [the trait docs](Self) for more information. + fn compute_stress_contraction_du( + &self, + u_grad: &OMatrix, + a: &OVector, + b: &OVector, + parameters: &Self::Parameters, + ) -> OMatrix { + self.compute_stress_contraction(&deformation_gradient(u_grad), a, b, parameters) + } + /// Compute the contraction for a number of vectors at the same time, with the given /// parameters. /// @@ -128,8 +214,132 @@ where self.compute_stress_contraction(deformation_gradient, a_I, b_J, parameters) }) } + + /// Compute the contraction for a number of vectors at the same time, with the given + /// parameters. + /// + /// Takes as input the displacement gradient $\nabla \vec u$ instead of the deformation + /// gradient $\vec F$. See [the trait docs](Self) for more information. + #[allow(non_snake_case)] + fn accumulate_stress_contractions_du_into( + &self, + output: DMatrixViewMut, + alpha: T, + u_grad: &OMatrix, + a: DVectorView, + b: DVectorView, + parameters: &Self::Parameters, + ) { + compute_batch_contraction(output, alpha, a, b, |a_I, b_J| { + self.compute_stress_contraction_du(&u_grad, a_I, b_J, parameters) + }) + } } +// TODO: Remove this or develop it further. The idea is to be able to +// easily implement materials only in terms of du/dX and have some mechanism for automatically +// producing the whole HyperelasticMaterial trait. An alternative, simpler direction, +// would be to simply have the HyperelasticMaterial trait only work directly with u_grad. + +// pub trait HyperelasticMaterialDu +// where +// T: Real, +// GeometryDim: DimName, +// DefaultAllocator: DimAllocator, +// { +// type Parameters: Clone + Default + 'static; +// +// /// TODO: Docs +// fn compute_energy_density_du( +// &self, +// u_grad: &OMatrix, +// parameters: &Self::Parameters, +// ) -> T; +// +// /// TODO: Docs +// fn compute_stress_tensor_du( +// &self, +// u_grad: &OMatrix, +// parameters: &Self::Parameters, +// ) -> OMatrix; +// +// /// TODO: Docs +// fn compute_stress_contraction_du( +// &self, +// u_grad: &OMatrix, +// a: &OVector, +// b: &OVector, +// parameters: &Self::Parameters, +// ) -> OMatrix; +// +// /// TODO: Docs +// #[allow(non_snake_case)] +// fn accumulate_stress_contractions_du_into( +// &self, +// output: DMatrixViewMut, +// alpha: T, +// u_grad: &OMatrix, +// a: DVectorView, +// b: DVectorView, +// parameters: &Self::Parameters, +// ) { +// compute_batch_contraction(output, alpha, a, b, |a_I, b_J| { +// self.compute_stress_contraction_du(&u_grad, a_I, b_J, parameters) +// }) +// } +// } + +// #[macro_export] +// macro_rules! impl_hyperelastic_from_du_impl { +// ($material:ty) => { +// impl HyperelasticMaterial for $material +// where +// T: Real, +// D: DimName, +// $material: HyperelasticMaterialDu, +// DefaultAllocator: DimAllocator, +// { +// type Parameters = <$material as HyperelasticMaterialDu>::Parameters; +// +// fn compute_energy_density(&self, deformation_gradient: &OMatrix, parameters: &Self::Parameters) -> T { +// let u_grad = crate::u_grad_from_F(deformation_gradient); +// >::compute_energy_density_du(self, &u_grad, parameters) +// } +// +// fn compute_energy_density_du(&self, u_grad: &OMatrix, parameters: &Self::Parameters) -> T { +// >::compute_energy_density_du(self, u_grad, parameters) +// } +// +// fn compute_stress_tensor(&self, deformation_gradient: &OMatrix, parameters: &Self::Parameters) -> OMatrix { +// let u_grad = crate::u_grad_from_F(deformation_gradient); +// >::compute_stress_tensor_du(self, &u_grad, parameters) +// } +// +// fn compute_stress_tensor_du(&self, u_grad: &OMatrix, parameters: &Self::Parameters) -> OMatrix { +// >::compute_stress_tensor_du(self, u_grad, parameters) +// } +// +// fn compute_stress_contraction(&self, deformation_gradient: &OMatrix, a: &OVector, b: &OVector, parameters: &Self::Parameters) -> OMatrix { +// let u_grad = crate::u_grad_from_F(deformation_gradient); +// >::compute_stress_contraction_du(self, &u_grad, a, b, parameters) +// } +// +// fn compute_stress_contraction_du(&self, u_grad: &OMatrix, a: &OVector, b: &OVector, parameters: &Self::Parameters) -> OMatrix { +// >::compute_stress_contraction_du(self, u_grad, a, b, parameters) +// } +// +// fn accumulate_stress_contractions_into(&self, output: DMatrixViewMut, alpha: T, deformation_gradient: &OMatrix, a: DVectorView, b: DVectorView, parameters: &Self::Parameters) { +// let u_grad = crate::u_grad_from_F(deformation_gradient); +// >::accumulate_stress_contractions_du_into(self, output, alpha, &u_grad, a, b, parameters) +// } +// +// fn accumulate_stress_contractions_du_into(&self, output: DMatrixViewMut, alpha: T, u_grad: &OMatrix, a: DVectorView, b: DVectorView, parameters: &Self::Parameters) { +// >::accumulate_stress_contractions_du_into(self, output, alpha, u_grad, a, b, parameters) +// } +// } +// } +// } + /// Helper function to ease implementation of [`HyperelasticMaterial::accumulate_stress_contractions_into`]. /// /// Often implementations of this method will tend to look very similar. This method is provided as a means @@ -226,8 +436,7 @@ where DefaultAllocator: DimAllocator, { fn compute_energy(&self, u_grad: &OMatrix, parameters: &Self::Parameters) -> T { - let f = u_grad.transpose() + OMatrix::::identity(); - self.0.compute_energy_density(&f, parameters) + self.0.compute_energy_density_du(u_grad, parameters) } } @@ -243,9 +452,21 @@ where u_grad: &OMatrix, parameters: &Self::Parameters, ) -> OMatrix { - let f = u_grad.transpose() + OMatrix::::identity(); - let p = self.0.compute_stress_tensor(&f, parameters); - p.transpose() + // The material operator g is related to the stress tensor P by g = P^T, + // so it's more convenient to implement the transpose opereation g^T = P, + // which is anyway what is called by the assembler + self.compute_elliptic_operator_transpose(u_grad, parameters) + .transpose() + } + + fn compute_elliptic_operator_transpose( + &self, + u_grad: &OMatrix, + parameters: &Self::Parameters, + ) -> OMatrix { + // We avoid forming the deformation gradient here so that we can avoid + // the loss of accuracy implied by forming `F = I + grad u` for small `grad u` + self.0.compute_stress_tensor_du(u_grad, parameters) } } @@ -263,8 +484,8 @@ where b: &OVector, parameters: &Self::Parameters, ) -> OMatrix { - let f = u_grad.transpose() + OMatrix::::identity(); - self.0.compute_stress_contraction(&f, a, b, parameters) + self.0 + .compute_stress_contraction_du(u_grad, a, b, parameters) } fn symmetry(&self) -> Symmetry { @@ -281,8 +502,29 @@ where b: DVectorView, parameters: &Self::Parameters, ) { - let f = u_grad.transpose() + OMatrix::::identity(); self.0 - .accumulate_stress_contractions_into(output, alpha, &f, a, b, parameters) + .accumulate_stress_contractions_du_into(output, alpha, u_grad, a, b, parameters) } } + +mod internal { + use fenris::nalgebra::{U1, U2, U3}; + + pub trait Sealed {} + + impl Sealed for U1 {} + impl Sealed for U2 {} + impl Sealed for U3 {} +} + +/// A fixed-size dimension corresponding to physical space. +/// +/// Physical dimensions are comprised of the dimensions $1$, $2$ and $3$. The primary utility +/// of this trait is to support writing generic code that needs to evaluate functions +/// whose implementation differs from dimension to dimension, and that might not have an easily +/// accessible n-dimensional variant. +pub trait PhysicalDim: internal::Sealed + SmallDim {} + +impl PhysicalDim for U1 {} +impl PhysicalDim for U2 {} +impl PhysicalDim for U3 {} diff --git a/fenris-solid/src/logdet.rs b/fenris-solid/src/logdet.rs new file mode 100644 index 00000000..2d8723c1 --- /dev/null +++ b/fenris-solid/src/logdet.rs @@ -0,0 +1,86 @@ +use crate::PhysicalDim; +use fenris::allocators::DimAllocator; +use fenris::nalgebra::{DefaultAllocator, Matrix1, Matrix2, Matrix3, OMatrix}; +use fenris::util::try_transmute_ref; +use fenris::Real; +use numeric_literals::replace_float_literals; + +/// Compute $\log(\det \vec F)$ for the deformation gradient $\vec F$ given $\pd{\vec u}{\vec X}$. +/// +/// The implementation may be far more accurate than first forming $\vec F$, which discards +/// valuable information for small deformations (small $\pd{\vec u}{\vec X}$). +/// +/// The implementation is based on the technique described in the +/// [libCEED documentation](https://libceed.org/en/latest/examples/solids). +#[allow(non_snake_case)] +#[replace_float_literals(T::from_f64(literal).unwrap())] +pub fn log_det_F(du_dX: &OMatrix) -> Option +where + T: Real, + D: PhysicalDim, + DefaultAllocator: DimAllocator, +{ + match D::USIZE { + 1 => { + let du_dX: &Matrix1 = try_transmute_ref(du_dX).unwrap(); + let det = du_dX[(0, 0)]; + (det > 0.0).then(|| det.ln()) + } + 2 => log_det_F_2d(try_transmute_ref(du_dX).unwrap()), + 3 => log_det_F_3d(try_transmute_ref(du_dX).unwrap()), + _ => unreachable!("Physical dimensions do not extend past 3 dimensions"), + } +} + +#[allow(non_snake_case)] +#[replace_float_literals(T::from_f64(literal).unwrap())] +fn log_det_F_2d(du_dX: &Matrix2) -> Option { + // See comments in 3D impl for more elaborate explanation + // Given a matrix A = [a, b; c, d] the determinant is + // det(A) = ad - bc + // Let U = du_dX be the displacement Jacobian and uij its entries. Then we can write our determinant as + // det(F) = det(I + U) = a*d - b*c = (1+u11)*(1+u22) - b*c + // = 1 + u11*u22 + u11 + u22 - b*c + // = 1 + γ + // and so + // log(det(F)) = log(1 + γ) = log1p(γ) + let u11 = du_dX[(0, 0)]; + let u22 = du_dX[(1, 1)]; + let b = du_dX[(0, 1)]; + let c = du_dX[(1, 0)]; + let gamma = u11 * u22 + u11 + u22 - b * c; + (gamma > -1.0).then(|| T::ln_1p(gamma)) +} + +#[allow(non_snake_case)] +#[replace_float_literals(T::from_f64(literal).unwrap())] +fn log_det_F_3d(du_dX: &Matrix3) -> Option { + // Given a matrix A = [a, b, c; d, e, f; g, h, i] the determinant is + // det(A) = aei + bfg + cdh - ceg - bdi - afh + // The first term is the product of the diagonals. Since F = I + du_dX, + // the diagonal of F will be close to 1 for small du_dX, + // and so the determinant is dominated by this first term. + // Let U = du_dX be the displacement Jacobian and uij its entries. Then we can write our determinant as + // det(F) = det(I + U) = (1+u11)*(1+u22)*(1+u33) + ... + // = 1 + u11*u22*u33 + u11*u22 + u11*u33 + u22*u33 + u11 + u22 + u33 + ... + // = 1 + γ + // and so + // log(det(F)) = log(1 + γ) = log1p(γ) + let u11 = du_dX[(0, 0)]; + let u22 = du_dX[(1, 1)]; + let u33 = du_dX[(2, 2)]; + let a = 1.0 + u11; + let e = 1.0 + u22; + let i = 1.0 + u33; + let b = du_dX[(0, 1)]; + let c = du_dX[(0, 2)]; + let d = du_dX[(1, 0)]; + let f = du_dX[(1, 2)]; + let g = du_dX[(2, 0)]; + let h = du_dX[(2, 1)]; + let gamma = u11 * u22 * u33 + u11 * u22 + u11 * u33 + u22 * u33 + u11 + u22 + u33 + b * f * g + c * d * h + - c * e * g + - b * d * i + - a * f * h; + (gamma > -1.0).then(|| T::ln_1p(gamma)) +} diff --git a/fenris-solid/src/materials.rs b/fenris-solid/src/materials.rs index 257bd103..b0e41557 100644 --- a/fenris-solid/src/materials.rs +++ b/fenris-solid/src/materials.rs @@ -1,7 +1,7 @@ -use crate::{compute_batch_contraction, HyperelasticMaterial}; +use crate::{compute_batch_contraction, log_det_F, u_grad_from_F, HyperelasticMaterial, PhysicalDim}; use fenris::allocators::DimAllocator; use fenris::nalgebra::{DMatrixViewMut, DVectorView, DefaultAllocator, DimName, OMatrix, OVector}; -use fenris::{Real, SmallDim}; +use fenris::Real; use numeric_literals::replace_float_literals; use serde::{Deserialize, Serialize}; @@ -236,25 +236,28 @@ pub struct NeoHookeanMaterial; impl HyperelasticMaterial for NeoHookeanMaterial where T: Real, - D: SmallDim, + D: PhysicalDim, DefaultAllocator: DimAllocator, { type Parameters = LameParameters; fn compute_energy_density(&self, deformation_gradient: &OMatrix, parameters: &Self::Parameters) -> T { - let LameParameters { mu, lambda } = parameters.clone(); - let F = deformation_gradient; - let J = F.determinant(); + let u_grad = u_grad_from_F(deformation_gradient); + self.compute_energy_density_du(&u_grad, parameters) + } - if J <= T::zero() { - T::from_f64(f64::INFINITY).expect("T must be able to represent infinity") + fn compute_energy_density_du(&self, u_grad: &OMatrix, parameters: &Self::Parameters) -> T { + let LameParameters { mu, lambda } = parameters.clone(); + let du_dX = u_grad.transpose(); + if let Some(logJ) = log_det_F(&du_dX) { + // Original expression + // 0.5 * mu * (I_C - d) - mu * logJ + (0.5 * lambda) * (logJ.powi(2)) + // We have that + // I_C - d = tr(C) - d = tr(C - I) = 2 tr(E) + let tr_E = du_dX.trace() + 0.5 * du_dX.norm_squared(); + mu * tr_E - mu * logJ + (0.5 * lambda) * (logJ.powi(2)) } else { - let C = F.transpose() * F; - let I_C = C.trace(); - let logJ = J.ln(); - // Note: 2D/3D need different constants to ensure rest state has zero energy - let d = T::from_usize(D::dim()).unwrap(); - mu / 2.0 * (I_C - d) - mu * logJ + (lambda / 2.0) * (logJ.powi(2)) + T::from_f64(f64::INFINITY).expect("T must be able to represent infinity") } } diff --git a/fenris-solid/tests/unit_tests/logdet.rs b/fenris-solid/tests/unit_tests/logdet.rs new file mode 100644 index 00000000..59aa1404 --- /dev/null +++ b/fenris-solid/tests/unit_tests/logdet.rs @@ -0,0 +1,50 @@ +use fenris::nalgebra; +use fenris::nalgebra::{matrix, vector, Matrix3, Rotation3}; +use fenris_solid::{log_det_F, u_grad_from_F}; +use matrixcompare::assert_scalar_eq; + +#[allow(non_snake_case)] +fn arbitrary_F() -> Matrix3 { + // We construct F by an ad-hoc SVD construction. This way we can control the singular values + // and determinant (since rotations have determinant 1) + let rot1 = Rotation3::from_scaled_axis(vector![1.0, 2.0, -3.0]) + .matrix() + .clone_owned(); + let rot2 = Rotation3::from_scaled_axis(vector![4.0, -5.0, 6.0]) + .matrix() + .clone_owned(); + let sigma = matrix![1.0 + 0.1, 0.0, 0.0; + 0.0, 1.0 - 0.2, 0.0; + 0.0, 0.0, 1.0 + 0.3]; + rot1 * sigma * rot2 +} + +#[test] +#[allow(non_snake_case)] +fn test_log_det_F() { + let F = arbitrary_F(); + let du_dX = F - Matrix3::identity(); + assert_scalar_eq!( + log_det_F(&du_dX).unwrap(), + F.determinant().ln(), + comp = abs, + tol = 1e-12 + ); +} + +#[test] +#[allow(non_snake_case)] +fn log_det_negative_determinant() { + let rot1 = Rotation3::from_scaled_axis(vector![1.0, 2.0, -3.0]) + .matrix() + .clone_owned(); + let rot2 = Rotation3::from_scaled_axis(vector![4.0, -5.0, 6.0]) + .matrix() + .clone_owned(); + let sigma = matrix![-1e-6, 0.0, 0.0; + 0.0, 1.0, 0.0; + 0.0, 0.0, 1.0 + 0.3]; + let F = rot1 * sigma * rot2; + let du_dX = u_grad_from_F(&F).transpose(); + assert!(log_det_F(&du_dX).is_none()); +} diff --git a/fenris-solid/tests/unit_tests/mod.rs b/fenris-solid/tests/unit_tests/mod.rs index ff7a7636..9b5ceca8 100644 --- a/fenris-solid/tests/unit_tests/mod.rs +++ b/fenris-solid/tests/unit_tests/mod.rs @@ -4,6 +4,7 @@ use fenris::nalgebra::{matrix, Matrix2, Matrix3, Point3}; use fenris_solid::materials::LameParameters; mod gravity_source; +mod logdet; mod material_elliptic_operator; mod materials; diff --git a/src/assembly/local/elliptic.rs b/src/assembly/local/elliptic.rs index 014ef8ca..f98f4c0e 100644 --- a/src/assembly/local/elliptic.rs +++ b/src/assembly/local/elliptic.rs @@ -522,8 +522,8 @@ where let mut output = MatrixViewMut::from_slice_generic(output.as_mut_slice(), Operator::SolutionDim::name(), Dyn(n)); - let g = operator.compute_elliptic_operator(&u_grad, data); - let g_t_j_inv_t = g.transpose() * j_inv_t; + let g_t = operator.compute_elliptic_operator_transpose(&u_grad, data); + let g_t_j_inv_t = g_t * j_inv_t; output.gemm(weight * j_det.abs(), &g_t_j_inv_t, &phi_grad_ref, T::one()); } diff --git a/src/assembly/operators.rs b/src/assembly/operators.rs index bd211423..ddba2679 100644 --- a/src/assembly/operators.rs +++ b/src/assembly/operators.rs @@ -1,5 +1,4 @@ use crate::allocators::BiDimAllocator; -use crate::nalgebra::allocator::Allocator; use crate::nalgebra::{DMatrixViewMut, DVectorView, DefaultAllocator, DimName, OMatrix, OVector, Scalar}; use crate::{Real, SmallDim, Symmetry}; @@ -22,7 +21,7 @@ pub trait EllipticOperator: Operator where T: Scalar, GeometryDim: SmallDim, - DefaultAllocator: Allocator, + DefaultAllocator: BiDimAllocator, { /// Compute the elliptic operator $g = g(\nabla u)$ with the provided /// [operator parameters](Operator::Parameters). @@ -31,6 +30,19 @@ where gradient: &OMatrix, parameters: &Self::Parameters, ) -> OMatrix; + + /// Compute the transpose $g^T$ of the elliptic operator $g$. + /// + /// See [`compute_elliptic_operator`](Self::compute_elliptic_operator). Implementing + /// this function can often avoid unnecessary transposition in consumers. + fn compute_elliptic_operator_transpose( + &self, + gradient: &OMatrix, + parameters: &Self::Parameters, + ) -> OMatrix { + self.compute_elliptic_operator(gradient, parameters) + .transpose() + } } /// A contraction operator encoding derivative information for an elliptic operator.