From 8a62780d9e836ed46a6c5c9192c38e4e7ad32ba7 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Sat, 27 Feb 2021 17:43:01 +0100
Subject: [PATCH] Modernize LocalIntegralEnergy
---
.../materials/localintegralenergy.hh | 18 ++++++++++--------
1 file changed, 10 insertions(+), 8 deletions(-)
diff --git a/dune/elasticity/materials/localintegralenergy.hh b/dune/elasticity/materials/localintegralenergy.hh
index e19d9c8..6ef3815 100644
--- a/dune/elasticity/materials/localintegralenergy.hh
+++ b/dune/elasticity/materials/localintegralenergy.hh
@@ -2,6 +2,7 @@
#define DUNE_ELASTICITY_MATERIALS_LOCALINTEGRALENERGY_HH
#include <dune/common/fmatrix.hh>
+#include <dune/common/transpose.hh>
#include <dune/geometry/quadraturerules.hh>
@@ -53,8 +54,9 @@ energy(const LocalView& localView,
field_type energy = 0;
// store gradients of shape functions and base functions
- std::vector<FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size());
- std::vector<FieldVector<DT,gridDim> > gradients(localFiniteElement.size());
+ using FiniteElement = typename LocalView::Tree::ChildType::FiniteElement;
+ using Jacobian = typename FiniteElement::Traits::LocalBasisType::Traits::JacobianType;
+ std::vector<Jacobian> jacobians(localFiniteElement.size());
int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
: localFiniteElement.localBasis().order() * gridDim;
@@ -65,23 +67,23 @@ energy(const LocalView& localView,
{
const DT integrationElement = element.geometry().integrationElement(qp.position());
- const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(qp.position());
+ const auto geometryJacobianIT = element.geometry().jacobianInverseTransposed(qp.position());
// Global position
auto x = element.geometry().global(qp.position());
// Get gradients of shape functions
- localFiniteElement.localBasis().evaluateJacobian(qp.position(), referenceGradients);
+ localFiniteElement.localBasis().evaluateJacobian(qp.position(), jacobians);
// compute gradients of base functions
- for (size_t i=0; i<gradients.size(); ++i)
- jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+ for (size_t i=0; i<jacobians.size(); ++i)
+ jacobians[i] = jacobians[i] * transpose(geometryJacobianIT);
// Deformation gradient
FieldMatrix<field_type,gridDim,gridDim> deformationGradient(0);
- for (size_t i=0; i<gradients.size(); i++)
+ for (size_t i=0; i<jacobians.size(); i++)
for (int j=0; j<gridDim; j++)
- deformationGradient[j].axpy(localConfiguration[ localView.tree().child(j).localIndex(i) ], gradients[i]);
+ deformationGradient[j].axpy(localConfiguration[ localView.tree().child(j).localIndex(i) ], jacobians[i][0]);
// Integrate the energy density
energy += qp.weight() * integrationElement * (*localDensity_)(x, deformationGradient);
}
--
GitLab