diff --git a/CHANGELOG.md b/CHANGELOG.md
index 663fa3503d2277a580c5b21362f6362585158746..556901d5fb71f9468ef36b402fb3c8ced40c84da 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,9 @@
# Master
+- The `SimoFoxEnergy` class (formerly known as `SimoFoxEnergyLocalStiffness`)
+ does not accept a surface load density anymore. Use the
+ `NeumannEnergy` class for assembling surface loads.
+
- The classes `LocalGeodesicFEFunction` and `LocalProjectedFEFunction`
now take a `dune-functions` basis as their first arguments, instead
of a `LocalFiniteElement`.
diff --git a/dune/gfe/assemblers/mixedgfeassembler.hh b/dune/gfe/assemblers/mixedgfeassembler.hh
index 426096d456df5da22e23df44c151535eeb5ce3a8..b38628c6299eed3d2f8f264a7b40e2c8b563280f 100644
--- a/dune/gfe/assemblers/mixedgfeassembler.hh
+++ b/dune/gfe/assemblers/mixedgfeassembler.hh
@@ -70,16 +70,16 @@ namespace Dune::GFE
* This is more efficient than computing them separately, because you need the gradient
* anyway to compute the Riemannian Hessian.
*/
- virtual void assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
- const std::vector<TargetSpace1>& configuration1,
- Dune::BlockVector<Dune::FieldVector<double, blocksize0> >& gradient0,
- Dune::BlockVector<Dune::FieldVector<double, blocksize1> >& gradient1,
- MatrixType& hessian,
- bool computeOccupationPattern=true) const;
+ void assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
+ const std::vector<TargetSpace1>& configuration1,
+ Dune::BlockVector<Dune::FieldVector<double, blocksize0> >& gradient0,
+ Dune::BlockVector<Dune::FieldVector<double, blocksize1> >& gradient1,
+ MatrixType& hessian,
+ bool computeOccupationPattern=true) const;
/** \brief Compute the energy of a deformation state */
- virtual double computeEnergy(const std::vector<TargetSpace0>& configuration0,
- const std::vector<TargetSpace1>& configuration1) const;
+ double computeEnergy(const std::vector<TargetSpace0>& configuration0,
+ const std::vector<TargetSpace1>& configuration1) const;
//protected:
void getMatrixPattern(Dune::MatrixIndexSet& nb00,
diff --git a/dune/gfe/assemblers/simofoxenergy.hh b/dune/gfe/assemblers/simofoxenergy.hh
index 424376beb1d0090b037dc20456f4fe5c6845b744..821dd1d198fefa6a4413ec6c89947754914cfc96 100644
--- a/dune/gfe/assemblers/simofoxenergy.hh
+++ b/dune/gfe/assemblers/simofoxenergy.hh
@@ -6,7 +6,6 @@
#include <dune/common/transpose.hh>
#include <dune/common/tuplevector.hh>
#include <dune/functions/functionspacebases/subspacebasis.hh>
-#include <dune/fufem/boundarypatch.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/gfe/linearalgebra.hh>
#include <dune/gfe/assemblers/localenergy.hh>
@@ -50,7 +49,7 @@ namespace Dune::GFE
* and director field.
*/
template <class Basis, template <typename, typename> typename LocalFEFunction, typename field_type = double>
- class SimoFoxEnergyLocalStiffness
+ class SimoFoxEnergy
: public LocalEnergy<Basis, ProductManifold<RealTuple<field_type, 3>,
UnitVector<field_type, 3> > >
{
@@ -71,13 +70,9 @@ namespace Dune::GFE
* \param parameters The material parameters
* \param x0 reference configuration
*/
- SimoFoxEnergyLocalStiffness(const Dune::ParameterTree ¶meters, const BoundaryPatch<GridView> *neumannBoundary,
- const std::function<Dune::FieldVector<double, 3>(Dune::FieldVector<double, 2>)> neumannFunction,
- const std::function<Dune::FieldVector<double, 3>(Dune::FieldVector<double, 2>)> volumeLoad,
- const Dune::TupleVector<std::vector<RealTuple<double, 3> >, std::vector<UnitVector<double, 3> > > &x0)
- : neumannBoundary_(neumannBoundary),
- neumannFunction_(neumannFunction),
- thickness_{parameters.template get<double>("thickness")}, // The sheeqll thickness
+ SimoFoxEnergy(const Dune::ParameterTree ¶meters,
+ const Dune::TupleVector<std::vector<RealTuple<double, 3> >, std::vector<UnitVector<double, 3> > > &x0)
+ : thickness_{parameters.template get<double>("thickness")}, // The sheeqll thickness
mu_{parameters.template get<double>("mu")}, // Lame constant 1
lambda_{parameters.template get<double>("lambda")}, // Lame constant 2
kappa_{parameters.template get<double>("kappa")}, // Shear correction factor
@@ -144,15 +139,6 @@ namespace Dune::GFE
/** \brief Calculate the Green-Lagrange strain components */
static Dune::FieldVector<RT, 8> calculateGreenLagrangianStrains(const KinematicVariables &kin);
- /** \brief Save all tangent base matrices for all nodes in one place*/
- Dune::BlockVector<Dune::FieldMatrix<field_type, 2, 3> > directorTangentSpaces;
-
- /** \brief The Neumann boundary */
- const BoundaryPatch<GridView> *neumannBoundary_;
-
- /** \brief The function implementing the Neumann data */
- const std::function<Dune::FieldVector<double, 3>(Dune::FieldVector<double, dimworld>)> neumannFunction_;
-
/** \brief The shell thickness */
double thickness_;
@@ -165,9 +151,6 @@ namespace Dune::GFE
/** \brief Material tangent matrix */
Dune::FieldMatrix<double, 8, 8> CMat_;
- /** \brief The function implementing a volume load */
- const std::function<Dune::FieldVector<double, 3>(Dune::FieldVector<double, dimworld>)> volumeLoad_;
-
/** \brief Stores the reference configuration of the midsurface and the director field */
const std::vector<RealTuple<double, 3> > &midSurfaceRefConfig;
const std::vector<UnitVector<double, 3> > &directorRefConfig;
@@ -181,7 +164,7 @@ namespace Dune::GFE
* gamma_1 and gamma_2 are the transverse shear in the two parametric directions
* Paper Equation 4.10 */
template <class Basis, template <typename, typename> typename LocalFEFunction, typename field_type>
- Dune::FieldVector<field_type, 8> SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::calculateGreenLagrangianStrains(
+ Dune::FieldVector<field_type, 8> SimoFoxEnergy<Basis, LocalFEFunction, field_type>::calculateGreenLagrangianStrains(
const KinematicVariables &kin)
{
Dune::FieldVector<RT, 8> egl;
@@ -210,28 +193,28 @@ namespace Dune::GFE
template <class Basis, template <typename, typename> typename LocalFEFunction, typename field_type>
template <typename Element, typename LocalDirectorFunction, typename LocalMidSurfaceFunction, typename LocalDirectorReferenceFunction,
typename LocalMidSurfaceReferenceFunction, typename IntegrationPointPosition>
- auto SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::kinematicVariablesFactory(
+ auto SimoFoxEnergy<Basis, LocalFEFunction, field_type>::kinematicVariablesFactory(
const Element &element, const LocalDirectorFunction &directorFunction, const LocalDirectorReferenceFunction &directorReferenceFunction,
const LocalMidSurfaceFunction &midSurfaceFunction, const LocalMidSurfaceReferenceFunction &midSurfaceReferenceFunction,
const LocalMidSurfaceFunction &midSurfaceDisplacementFunction, const IntegrationPointPosition &quadPos)
{
KinematicVariables kin{};
- const auto jInvT = element.geometry().jacobianInverseTransposed(quadPos);
+ const auto geometryJacobianInverse = element.geometry().jacobianInverse(quadPos);
kin.t = directorFunction.evaluate(quadPos).globalCoordinates();
kin.t0 = directorReferenceFunction.evaluate(quadPos).globalCoordinates();
- kin.a1anda2 = transpose(midSurfaceFunction.evaluateDerivative(quadPos) * transpose(jInvT));
- kin.A1andA2 = transpose(midSurfaceReferenceFunction.evaluateDerivative(quadPos) * transpose(jInvT));
- kin.ud1andud2 = transpose(midSurfaceDisplacementFunction.evaluateDerivative(quadPos) * transpose(jInvT));
- kin.t0d1Andt0d2 = transpose(directorReferenceFunction.evaluateDerivative(quadPos) * transpose(jInvT));
- kin.td1Andtd2 = transpose(directorFunction.evaluateDerivative(quadPos) * transpose(jInvT));
+ kin.a1anda2 = transpose(midSurfaceFunction.evaluateDerivative(quadPos) * geometryJacobianInverse);
+ kin.A1andA2 = transpose(midSurfaceReferenceFunction.evaluateDerivative(quadPos) * geometryJacobianInverse);
+ kin.ud1andud2 = transpose(midSurfaceDisplacementFunction.evaluateDerivative(quadPos) * geometryJacobianInverse);
+ kin.t0d1Andt0d2 = transpose(directorReferenceFunction.evaluateDerivative(quadPos) * geometryJacobianInverse);
+ kin.td1Andtd2 = transpose(directorFunction.evaluateDerivative(quadPos) * geometryJacobianInverse);
return kin;
}
template <class Basis, template <typename, typename> typename LocalFEFunction, typename field_type>
- auto SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::getReferenceLocalConfigurations(
+ auto SimoFoxEnergy<Basis, LocalFEFunction, field_type>::getReferenceLocalConfigurations(
const typename Basis::LocalView &localView) const {
using namespace Dune::Indices;
const int nDofs0 = localView.tree().child(_0, 0).finiteElement().size();
@@ -269,8 +252,8 @@ namespace Dune::GFE
* see for details Paper Equation 4.11,4.10 and 10.1
*/
template <class Basis, template <typename, typename> typename LocalFEFunction, typename field_type>
- typename SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::RT
- SimoFoxEnergyLocalStiffness<Basis, LocalFEFunction, field_type>::energy(
+ typename SimoFoxEnergy<Basis, LocalFEFunction, field_type>::RT
+ SimoFoxEnergy<Basis, LocalFEFunction, field_type>::energy(
const typename Basis::LocalView &localView,
const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients &localConfiguration) const
{
@@ -327,36 +310,6 @@ namespace Dune::GFE
energy += 0.5 * Egl * (CMat_ * Egl) * curQuad.weight() * integrationElement;
}
- //////////////////////////////////////////////////////////////////////////////
- // Assemble boundary contributions
- //////////////////////////////////////////////////////////////////////////////
-
- if (not neumannFunction_) return energy;
-
- for (auto &&it : intersections(neumannBoundary_->gridView(), element))
- {
- if (not neumannBoundary_ or not neumannBoundary_->contains(it)) continue;
-
- const auto &quadLine = QuadratureRules<DT, gridDim - 1>::rule(it.type(), quadOrder);
-
- for (const auto &curQuad : quadLine)
- {
- // Local position of the quadrature point
- const FieldVector<DT, gridDim> &quadPos = it.geometryInInside().global(curQuad.position());
-
- const DT integrationElement = it.geometry().integrationElement(curQuad.position());
-
- // The value of the local function
- RealTuple<field_type, 3> deformationValue = localMidSurfaceFunction.evaluate(quadPos);
-
- // Value of the Neumann data at the current position
- auto neumannValue = neumannFunction_(it.geometry().global(curQuad.position()));
-
- // Only translational dofs are affected by the Neumann force
- for (size_t i = 0; i < neumannValue.size(); i++)
- energy -= (neumannValue[i] * deformationValue.globalCoordinates()[i]) * curQuad.weight() * integrationElement;
- }
- }
return energy;
}
diff --git a/dune/gfe/densities/duneelasticitydensity.hh b/dune/gfe/densities/duneelasticitydensity.hh
index 4d104077e7ec9d6f90e63d57e429b6cefe99b8fa..b3e1de3a36c03f709134c4bc234c9fecda63343e 100644
--- a/dune/gfe/densities/duneelasticitydensity.hh
+++ b/dune/gfe/densities/duneelasticitydensity.hh
@@ -73,7 +73,7 @@ namespace Dune::GFE
}
// Construct a copy of this density but using 'adouble' as the number type
- virtual std::unique_ptr<LocalDensity<ElementOrIntersection,ATargetSpace> > makeActiveDensity() const
+ virtual std::unique_ptr<LocalDensity<ElementOrIntersection,ATargetSpace> > makeActiveDensity() const override
{
// The active dune-elasticity density
auto activeDensity = elasticityDensity_->makeActiveDensity();
diff --git a/dune/gfe/densities/localdensity.hh b/dune/gfe/densities/localdensity.hh
index 6bc6b4c220447d2f01c200f0c6f169b35f24155c..8f6e50abc12ac99d4af765979d226dbaf8be92b9 100644
--- a/dune/gfe/densities/localdensity.hh
+++ b/dune/gfe/densities/localdensity.hh
@@ -63,7 +63,7 @@ namespace Dune::GFE
}
- ~LocalDensity()
+ virtual ~LocalDensity()
{
myfree3(densityTangent_);
myfree3(Yppp_);
diff --git a/src/simofoxshell.cc b/src/simofoxshell.cc
index f437dab615561da64944b1d6e55aa1db5af0cf64..72e7291b60cdf8c28482907634e9edaeac960a73 100644
--- a/src/simofoxshell.cc
+++ b/src/simofoxshell.cc
@@ -28,10 +28,12 @@
#include <dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh>
#include <dune/gfe/assemblers/simofoxenergy.hh>
#include <dune/gfe/assemblers/mixedgfeassembler.hh>
+#include <dune/gfe/assemblers/sumenergy.hh>
#include <dune/gfe/functions/embeddedglobalgfefunction.hh>
#include <dune/gfe/functions/localgeodesicfefunction.hh>
#include <dune/gfe/functions/localprojectedfefunction.hh>
#include <dune/gfe/mixedriemanniantrsolver.hh>
+#include <dune/gfe/neumannenergy.hh>
#include <dune/gfe/spaces/unitvector.hh>
#if !MIXED_SPACE
@@ -195,10 +197,10 @@ int main(int argc, char *argv[]) try
}
BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
- BoundaryPatch<GridView> neumannBoundary(gridView, neumannVertices);
+ auto neumannBoundary = std::make_shared<BoundaryPatch<GridView> >(gridView, neumannVertices);
if (mpiHelper.rank() == 0)
- std::cout << "Neumann boundary has " << neumannBoundary.numFaces() << " faces\n";
+ std::cout << "Neumann boundary has " << neumannBoundary->numFaces() << " faces\n";
BitSetVector<1> deformationDirichletNodes(midsurfaceFEBasis.size(), false);
#if DUNE_VERSION_GTE(DUNE_FUFEM, 2, 10)
@@ -209,9 +211,9 @@ int main(int argc, char *argv[]) try
BitSetVector<1> neumannNodes(midsurfaceFEBasis.size(), false);
#if DUNE_VERSION_GTE(DUNE_FUFEM, 2, 10)
- Fufem::markBoundaryPatchDofs(neumannBoundary, directorFEBasis, neumannNodes);
+ Fufem::markBoundaryPatchDofs(*neumannBoundary, directorFEBasis, neumannNodes);
#else
- constructBoundaryDofs(neumannBoundary, directorFEBasis, neumannNodes);
+ constructBoundaryDofs(*neumannBoundary, directorFEBasis, neumannNodes);
#endif
BitSetVector<3> deformationDirichletDofs(midsurfaceFEBasis.size(), false);
@@ -316,19 +318,25 @@ int main(int argc, char *argv[]) try
materialParameters.report();
}
- // Assembler using ADOL-C
- auto simoFoxEnergyLocalStiffness
- = std::make_shared<GFE::SimoFoxEnergyLocalStiffness<decltype(compositeBasis),
+ // The total energy on one element
+ auto sumEnergy = std::make_shared<GFE::SumEnergy<decltype(compositeBasis), GFE::RealTuple<adouble,3>,GFE::UnitVector<adouble,3> > >();
+
+ // Internal energy of the shell
+ auto simoFoxEnergy
+ = std::make_shared<GFE::SimoFoxEnergy<decltype(compositeBasis),
LocalFEFunction,
- adouble> > (materialParameters,
- &neumannBoundary,
- neumannFunction,
- nullptr, x0);
+ adouble> > (materialParameters, x0);
+ sumEnergy->addLocalEnergy(simoFoxEnergy);
+
+ // The Neumann surface load term
+ auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<decltype(compositeBasis), GFE::RealTuple<adouble,3>, GFE::UnitVector<adouble,3> > >(neumannBoundary,neumannFunction);
+ sumEnergy->addLocalEnergy(neumannEnergy);
using TargetSpace = GFE::ProductManifold<GFE::RealTuple<double,3>,GFE::UnitVector<double,3> >;
+ // Use ADOL-C to assemble the element stiffness matrices
GFE::LocalGeodesicFEADOLCStiffness<decltype(compositeBasis),
- TargetSpace> localGFEADOLCStiffness(simoFoxEnergyLocalStiffness);
+ TargetSpace> localGFEADOLCStiffness(sumEnergy);
GFE::MixedGFEAssembler<decltype(compositeBasis),TargetSpace> assembler(compositeBasis, localGFEADOLCStiffness);
////////////////////////////////////////////////////////