From d984ed41f20f99737f7ec46894586d9ef6ce0648 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Wed, 18 Sep 2024 07:25:58 +0200
Subject: [PATCH] Move external loads out of NonplanarCosseratShellEnergy
That is one step towards splitting off the energy density
into a separate class.
Also, it helps to unify the code paths for the different grid dimensions
in cosserat-continuum.cc
---
.../nonplanarcosseratshellenergy.hh | 115 +-----------------
src/cosserat-continuum.cc | 98 ++++++++-------
test/nonplanarcosseratenergytest.cc | 6 +-
3 files changed, 56 insertions(+), 163 deletions(-)
diff --git a/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh b/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
index 9e5a8765..6fd7055a 100644
--- a/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
+++ b/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
@@ -95,14 +95,8 @@ public:
* \param stressFreeStateGridFunction Pointer to a parametrization representing the Cosserat shell in a stress-free state
*/
NonplanarCosseratShellEnergy(const Dune::ParameterTree& parameters,
- const StressFreeStateGridFunction* stressFreeStateGridFunction,
- const BoundaryPatch<GridView>* neumannBoundary,
- const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> neumannFunction,
- const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> volumeLoad)
- : stressFreeStateGridFunction_(stressFreeStateGridFunction),
- neumannBoundary_(neumannBoundary),
- neumannFunction_(neumannFunction),
- volumeLoad_(volumeLoad)
+ const StressFreeStateGridFunction* stressFreeStateGridFunction)
+ : stressFreeStateGridFunction_(stressFreeStateGridFunction)
{
// The shell thickness
thickness_ = parameters.template get<double>("thickness");
@@ -211,15 +205,6 @@ public:
/** \brief The geometry of the reference deformation used for assembling */
const StressFreeStateGridFunction* stressFreeStateGridFunction_;
-
- /** \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 function implementing a volume load */
- const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> volumeLoad_;
};
template <class Basis, int dim, class field_type, class StressFreeStateGridFunction>
@@ -416,54 +401,6 @@ energy(const typename Basis::LocalView& localView,
// Add energy density
energy += quad[pt].weight() * integrationElement * energyDensity;
-
- ///////////////////////////////////////////////////////////
- // Volume load contribution
- ///////////////////////////////////////////////////////////
-
- if (not volumeLoad_)
- continue;
-
- // Value of the volume load density at the current position
- Dune::FieldVector<double,3> volumeLoadDensity = volumeLoad_(geometry.global(quad[pt].position()));
-
- // Only translational dofs are affected by the volume load
- for (size_t i=0; i<volumeLoadDensity.size(); i++)
- energy += (volumeLoadDensity[i] * value[_0].globalCoordinates()[i]) * quad[pt].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& quad = Dune::QuadratureRules<DT, gridDim-1>::rule(it.type(), quadOrder);
-
- for (size_t pt=0; pt<quad.size(); pt++)
- {
- // Local position of the quadrature point
- const Dune::FieldVector<DT,gridDim>& quadPos = it.geometryInInside().global(quad[pt].position());
-
- const DT integrationElement = it.geometry().integrationElement(quad[pt].position());
-
- // The value of the local function
- Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > value = localGeodesicFEFunction.evaluate(quadPos);
-
- // Value of the Neumann data at the current position
- Dune::FieldVector<double,3> neumannValue = neumannFunction_(it.geometry().global(quad[pt].position()));
-
- // Only translational dofs are affected by the Neumann force
- for (size_t i=0; i<neumannValue.size(); i++)
- energy += (neumannValue[i] * value[_0].globalCoordinates()[i]) * quad[pt].weight() * integrationElement;
- }
}
return energy;
@@ -659,54 +596,6 @@ energy(const typename Basis::LocalView& localView,
// Add energy density
energy += quad[pt].weight() * integrationElement * energyDensity;
-
- ///////////////////////////////////////////////////////////
- // Volume load contribution
- ///////////////////////////////////////////////////////////
-
- if (not volumeLoad_)
- continue;
-
- // Value of the volume load density at the current position
- Dune::FieldVector<double,3> volumeLoadDensity = volumeLoad_(geometry.global(quad[pt].position()));
-
- // Only translational dofs are affected by the volume load
- for (size_t i=0; i<volumeLoadDensity.size(); i++)
- energy += (volumeLoadDensity[i] * deformationValue[i]) * quad[pt].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& quad = Dune::QuadratureRules<DT, gridDim-1>::rule(it.type(), quadOrder);
-
- for (size_t pt=0; pt<quad.size(); pt++)
- {
- // Local position of the quadrature point
- const Dune::FieldVector<DT,gridDim>& quadPos = it.geometryInInside().global(quad[pt].position());
-
- const DT integrationElement = it.geometry().integrationElement(quad[pt].position());
-
- // The value of the local function
- RealTuple<field_type,dim> deformationValue = localDeformationGFEFunction.evaluate(quadPos);
-
- // Value of the Neumann data at the current position
- Dune::FieldVector<double,3> neumannValue = neumannFunction_(it.geometry().global(quad[pt].position()));
-
- // Only translational dofs are affected by the Neumann force
- for (size_t i=0; i<neumannValue.size(); i++)
- energy += (neumannValue[i] * deformationValue[i]) * quad[pt].weight() * integrationElement;
- }
}
return energy;
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index 5c1c451d..f4250ceb 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -546,48 +546,69 @@ int main (int argc, char *argv[]) try
if (orientationDirichletDofs[i][0])
x[_1][i].set(dOV[i]);
- if (dim==dimworld) {
- auto sumEnergy = std::make_shared<GFE::SumEnergy<CompositeBasis, RealTuple<adouble,3>,Rotation<adouble,3> > >();
+ ////////////////////////////////////////////////////////////////
+ // Build the assembler for the tangent problems
+ ////////////////////////////////////////////////////////////////
+
+ // Construct the interpolation rule, i.e., the geometric finite element
+ using ScalarDeformationLocalFiniteElement = decltype(compositeBasis.localView().tree().child(_0,0).finiteElement());
+ using ScalarRotationLocalFiniteElement = decltype(compositeBasis.localView().tree().child(_1,0).finiteElement());
- // The Cosserat shell energy
- using ScalarDeformationLocalFiniteElement = decltype(compositeBasis.localView().tree().child(_0,0).finiteElement());
- using ScalarRotationLocalFiniteElement = decltype(compositeBasis.localView().tree().child(_1,0).finiteElement());
+ using AInterpolationRule = std::tuple<LocalGeodesicFEFunction<dim, double, ScalarDeformationLocalFiniteElement, RealTuple<adouble,3> >,
+ LocalGeodesicFEFunction<dim, double, ScalarRotationLocalFiniteElement, Rotation<adouble,3> > >;
- using AInterpolationRule = std::tuple<LocalGeodesicFEFunction<dim, double, ScalarDeformationLocalFiniteElement, RealTuple<adouble,3> >,
- LocalGeodesicFEFunction<dim, double, ScalarRotationLocalFiniteElement, Rotation<adouble,3> > >;
+ using ATargetSpace = typename TargetSpace::rebind<adouble>::other;
- using ATargetSpace = typename TargetSpace::rebind<adouble>::other;
+ using LocalCoordinate = typename GridType::Codim<0>::Entity::Geometry::LocalCoordinate;
- using LocalCoordinate = typename GridType::Codim<0>::Entity::Geometry::LocalCoordinate;
+ // The energy on one element
+ auto sumEnergy = std::make_shared<GFE::SumEnergy<CompositeBasis, RealTuple<adouble,3>,Rotation<adouble,3> > >();
+
+ if (dim==dimworld) {
auto cosseratDensity = createDensity<LocalCoordinate>(materialParameters);
auto localCosseratEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis,AInterpolationRule,ATargetSpace> >(cosseratDensity);
sumEnergy->addLocalEnergy(localCosseratEnergy);
+ }
+ else
+ {
+ static_assert((dim==dimworld) || (dim==2 && dimworld==3),
+ "For dim!=dimworld, only the case dim==2 and dimworld==3 is supported!");
+
+ auto localCosseratEnergy
+ = std::make_shared<NonplanarCosseratShellEnergy<CompositeBasis, 3, adouble, decltype(creator)> >(materialParameters,
+ &creator);
- // The Neumann surface load term
- auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<adouble,3>, Rotation<adouble,3> > >(neumannBoundary,neumannFunction);
- sumEnergy->addLocalEnergy(neumannEnergy);
-
- // The volume load term
- // TODO: There is a bug here: The volumeLoad function is currently defined in global coordinates,
- // but the density expects it to be in local coordinates with respect to the element being
- // integrated over. I have to think about where the binding should happen.
- // Practically, the bug does not really show, because all our volume loads
- // are constant in space anyway.
- auto volumeLoadDensity = std::make_shared<GFE::CosseratVolumeLoadDensity<LocalCoordinate,adouble> >(volumeLoad);
- auto volumeLoadEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, AInterpolationRule, ATargetSpace> >(volumeLoadDensity);
- sumEnergy->addLocalEnergy(volumeLoadEnergy);
-
- // The local assembler
- LocalGeodesicFEADOLCStiffness<CompositeBasis,TargetSpace> localGFEADOLCStiffness(sumEnergy,
- adolcScalarMode);
-
- MixedGFEAssembler<CompositeBasis,TargetSpace> mixedAssembler(compositeBasis, localGFEADOLCStiffness);
-
- ////////////////////////////////////////////
- // Set up the solver
- ////////////////////////////////////////////
+ sumEnergy->addLocalEnergy(localCosseratEnergy);
+ }
+
+ // The Neumann surface load term
+ auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<adouble,3>, Rotation<adouble,3> > >(neumannBoundary,neumannFunction);
+ sumEnergy->addLocalEnergy(neumannEnergy);
+
+ // The volume load term
+ // TODO: There is a bug here: The volumeLoad function is currently defined in global coordinates,
+ // but the density expects it to be in local coordinates with respect to the element being
+ // integrated over. I have to think about where the binding should happen.
+ // Practically, the bug does not really show, because all our volume loads
+ // are constant in space anyway.
+ auto volumeLoadDensity = std::make_shared<GFE::CosseratVolumeLoadDensity<LocalCoordinate,adouble> >(volumeLoad);
+ auto volumeLoadEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, AInterpolationRule, ATargetSpace> >(volumeLoadDensity);
+ sumEnergy->addLocalEnergy(volumeLoadEnergy);
+
+ // The local assembler
+ LocalGeodesicFEADOLCStiffness<CompositeBasis,TargetSpace> localGFEADOLCStiffness(sumEnergy,
+ adolcScalarMode);
+
+ MixedGFEAssembler<CompositeBasis,TargetSpace> mixedAssembler(compositeBasis, localGFEADOLCStiffness);
+
+ ////////////////////////////////////////////
+ // Set up the solver
+ ////////////////////////////////////////////
+ // TODO: There is no need why the solver setup code should depend on the grid dimension
+ if (dim==dimworld)
+ {
#if MIXED_SPACE
if (parameterSet.get<std::string>("solvertype", "trustRegion") == "trustRegion")
{
@@ -707,19 +728,6 @@ int main (int argc, char *argv[]) try
}
#endif
} else { //dim != dimworld
- using StiffnessType = LocalGeodesicFEADOLCStiffness<CompositeBasis, TargetSpace>;
- std::shared_ptr<StiffnessType> localGFEStiffness;
-
-#if HAVE_DUNE_CURVEDGEOMETRY && WORLD_DIM == 3 && GRID_DIM == 2
- auto localCosseratEnergy = std::make_shared<NonplanarCosseratShellEnergy<CompositeBasis, 3, adouble, decltype(creator)> >(materialParameters,
- &creator,
- neumannBoundary.get(),
- neumannFunction,
- volumeLoad);
-
- localGFEStiffness = std::make_shared<StiffnessType>(localCosseratEnergy, adolcScalarMode);
-#endif
- MixedGFEAssembler<CompositeBasis,TargetSpace> mixedAssembler(compositeBasis, localGFEStiffness);
#if MIXED_SPACE
MixedRiemannianTrustRegionSolver<GridType,
CompositeBasis,
diff --git a/test/nonplanarcosseratenergytest.cc b/test/nonplanarcosseratenergytest.cc
index a23461d4..44a5d254 100644
--- a/test/nonplanarcosseratenergytest.cc
+++ b/test/nonplanarcosseratenergytest.cc
@@ -124,11 +124,7 @@ double calculateEnergy(const FlatGridView& flatGridView,
double,
GridGeometry>;
- ShellEnergy nonplanarCosseratShellEnergy(materialParameters,
- &curvedGridGeometry,
- nullptr,
- nullptr,
- nullptr);
+ ShellEnergy nonplanarCosseratShellEnergy(materialParameters, &curvedGridGeometry);
///////////////////////////////////////////////////
// Compute the energy
--
GitLab