diff --git a/dune/gfe/assemblers/localintegralenergy.hh b/dune/gfe/assemblers/localintegralenergy.hh
index 6316fe823fcf46fe411d7e3b3bf7db97575d28cf..b10d519db2b8e18660ed5431e16553d775223676 100644
--- a/dune/gfe/assemblers/localintegralenergy.hh
+++ b/dune/gfe/assemblers/localintegralenergy.hh
@@ -2,6 +2,7 @@
#define DUNE_GFE_ASSEMBLERS_LOCALINTEGRALENERGY_HH
#include <dune/common/fmatrix.hh>
+#include <dune/common/version.hh>
#include <dune/geometry/quadraturerules.hh>
@@ -11,6 +12,34 @@
namespace Dune::GFE {
+#if ! DUNE_VERSION_GTE(DUNE_LOCALFUNCTIONS, 2, 10)
+ namespace Impl
+ {
+ template <class Basis>
+ class LocalFiniteElementFactory
+ {
+ public:
+ static auto get(const typename Basis::LocalView& localView)
+ -> decltype(localView.tree().child(0).finiteElement())
+ {
+ return localView.tree().child(0).finiteElement();
+ }
+ };
+
+ /** \brief Specialize for scalar bases, here we cannot call tree().child() */
+ template <class GridView, int order>
+ class LocalFiniteElementFactory<Dune::Functions::LagrangeBasis<GridView,order> >
+ {
+ public:
+ static auto get(const typename Dune::Functions::LagrangeBasis<GridView,order>::LocalView& localView)
+ -> decltype(localView.tree().finiteElement())
+ {
+ return localView.tree().finiteElement();
+ }
+ };
+ }
+#endif
+
/** \brief An energy given as an integral over a density
*
* \tparam Basis The scalar finite element basis used to construct the interpolation rule
@@ -44,9 +73,27 @@ namespace Dune::GFE {
{
RT energy = 0;
- if constexpr (not Impl::LocalEnergyTypes<TargetSpace>::isProductManifold)
+ if constexpr (Basis::LocalView::Tree::isLeaf || Basis::LocalView::Tree::isPower)
{
- const auto& localFiniteElement = localView.tree().finiteElement();
+#if DUNE_VERSION_GTE(DUNE_LOCALFUNCTIONS, 2, 10)
+ // Get an appropriate scalar local finite element, to construct the interpolation rule with
+ // TODO: This is not a good design, for several reasons:
+ // * The interpolation rule could want to have state beyond what we know here.
+ // It should therefore be constructed outside of the LocalIntegralEnergy class
+ // * I don't really see why Basis should be allowed to be scalar-valued
+ // to begin with, but a lot of code still currently does that.
+ auto lfeGetter = [&localView]()
+ {
+ if constexpr (Basis::LocalView::Tree::isPower)
+ return localView.tree().child(0).finiteElement();
+ else
+ return localView.tree().finiteElement();
+ };
+
+ const auto& localFiniteElement = lfeGetter();
+#else
+ const auto& localFiniteElement = Impl::LocalFiniteElementFactory<Basis>::get(localView);
+#endif
LocalInterpolationRule localInterpolationRule(localFiniteElement,localConfiguration);
int quadOrder = (localFiniteElement.type().isSimplex())
@@ -97,6 +144,11 @@ namespace Dune::GFE {
}
}
}
+ else
+ {
+ // You need a scalar basis or a power basis when calling this method.
+ std::abort();
+ }
return energy;
}
@@ -106,7 +158,7 @@ namespace Dune::GFE {
{
RT energy = 0;
- if constexpr (Impl::LocalEnergyTypes<TargetSpace>::isProductManifold)
+ if constexpr (Impl::LocalEnergyTypes<TargetSpace>::isProductManifold && Basis::LocalView::Tree::isComposite)
{
static_assert(TargetSpace::size() == 2,
"LocalIntegralEnergy only implemented for product spaces with two factors!");
@@ -177,6 +229,8 @@ namespace Dune::GFE {
derivative);
}
}
+ else
+ DUNE_THROW(Dune::NotImplemented, "Non-product manifold or non-composite basis");
return energy;
}
diff --git a/test/localadolcstiffnesstest.cc b/test/localadolcstiffnesstest.cc
index 2cd55096f12f9d347cbd5e45b62939032f921b20..b8cefbd4d6d5ee9601d6b4133e528df4b818f1bd 100644
--- a/test/localadolcstiffnesstest.cc
+++ b/test/localadolcstiffnesstest.cc
@@ -37,11 +37,11 @@ typedef double FDType;
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
-#include <dune/gfe/assemblers/localgeodesicfestiffness.hh>
#include <dune/gfe/localgeodesicfefunction.hh>
-#include <dune/gfe/assemblers/cosseratenergystiffness.hh>
+#include <dune/gfe/assemblers/localintegralenergy.hh>
#include <dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh>
#include <dune/gfe/assemblers/localgeodesicfefdstiffness.hh>
+#include <dune/gfe/densities/planarcosseratshelldensity.hh>
#include <dune/gfe/spaces/productmanifold.hh>
#include <dune/gfe/spaces/realtuple.hh>
#include <dune/gfe/spaces/rotation.hh>
@@ -104,12 +104,25 @@ int main (int argc, char *argv[]) try
typedef GridType::LeafGridView GridView;
GridView gridView = grid.leafGridView();
+ ///////////////////////////////////////////////
+ // Construct the basis for the tangent spaces
+ ///////////////////////////////////////////////
typedef Functions::LagrangeBasis<GridView,1> FEBasis;
FEBasis feBasis(gridView);
- // /////////////////////////////////////////
- // Read Dirichlet values
- // /////////////////////////////////////////
+ using namespace Functions::BasisFactory;
+
+ const int dimRotation = Rotation<double,3>::TangentVector::dimension;
+
+ auto tangentBasis = makeBasis(
+ gridView,
+ power<3+dimRotation>(
+ lagrange<1>()
+ )
+ );
+
+ using TangentBasis = decltype(tangentBasis);
+
// //////////////////////////
// Initial iterate
@@ -139,7 +152,6 @@ int main (int argc, char *argv[]) try
};
std::vector<FieldVector<double,3> > v;
- using namespace Functions::BasisFactory;
auto powerBasis = makeBasis(
gridView,
@@ -165,29 +177,45 @@ int main (int argc, char *argv[]) try
materialParameters["L_c"] = "1";
materialParameters["q"] = "2";
materialParameters["kappa"] = "1";
- materialParameters["b1"] = "1";
- materialParameters["b2"] = "1";
- materialParameters["b3"] = "1";
- ///////////////////////////////////////////////////////////////////////
- // Assemblers for the Riemannian derivatives
- ///////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////
+ // Assembler using ADOL-C
+ //////////////////////////////////////////////////////
+
+ // The 'active' target spaces, i.e., the number type is replaced by adouble
+ using ATargetSpace = typename TargetSpace::template rebind<adouble>::other;
- // Assembler using ADOL-C
- auto cosseratLocalEnergy = std::make_shared<CosseratEnergyLocalStiffness<FEBasis, 3,adouble> >(materialParameters);
+ // Select geometric finite element interpolation method
+ using AInterpolationRule = LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
- LocalGeodesicFEADOLCStiffness<FEBasis,
- TargetSpace> localGFEADOLCStiffness(cosseratLocalEnergy);
+ auto activeDensity = std::make_shared<GFE::PlanarCosseratShellDensity<GridType::Codim<0>::Entity::Geometry::LocalCoordinate, adouble> >(materialParameters);
- // Assembler using finite differences
- CosseratEnergyLocalStiffness<FEBasis, 3,FDType>
- cosseratEnergyFDLocalStiffness(materialParameters);
- LocalGeodesicFEFDStiffness<FEBasis,
+ auto activeCosseratLocalEnergy = std::make_shared<GFE::LocalIntegralEnergy<TangentBasis,AInterpolationRule,ATargetSpace> >(activeDensity);
+
+ // The actual assembler
+ LocalGeodesicFEADOLCStiffness<TangentBasis,
+ TargetSpace> localGFEADOLCStiffness(activeCosseratLocalEnergy);
+
+ //////////////////////////////////////////////////////
+ // Assembler using finite differences
+ //////////////////////////////////////////////////////
+
+ // Select geometric finite element interpolation method
+ using InterpolationRule = LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
+
+ auto cosseratDensity = std::make_shared<GFE::PlanarCosseratShellDensity<GridType::Codim<0>::Entity::Geometry::LocalCoordinate, double> >(materialParameters);
+
+ auto cosseratLocalEnergy = std::make_shared<GFE::LocalIntegralEnergy<TangentBasis,InterpolationRule,TargetSpace> >(cosseratDensity);
+
+ // The actual assembler
+ LocalGeodesicFEFDStiffness<TangentBasis,
TargetSpace,
- FDType> localGFEFDStiffness(&cosseratEnergyFDLocalStiffness);
+ FDType> localGFEFDStiffness(cosseratLocalEnergy.get());
- // Compute and compare matrices
+ ////////////////////////////////////////////////////////
+ // Compute and compare tangent matrices
+ ////////////////////////////////////////////////////////
for (const auto& element : Dune::elements(gridView))
{
std::cout << " ++++ element " << gridView.indexSet().index(element) << " ++++" << std::endl;
@@ -195,6 +223,9 @@ int main (int argc, char *argv[]) try
auto localView = feBasis.localView();
localView.bind(element);
+ auto tangentLocalView = tangentBasis.localView();
+ tangentLocalView.bind(element);
+
const int numOfBaseFct = localView.size();
// Extract local configuration
@@ -210,12 +241,12 @@ int main (int argc, char *argv[]) try
Matrix<double> localRiemannianADHessian;
Matrix<double> localRiemannianFDHessian;
- localGFEADOLCStiffness.assembleGradientAndHessian(localView,
+ localGFEADOLCStiffness.assembleGradientAndHessian(tangentLocalView,
localSolution,
localRiemannianADGradient,
localRiemannianADHessian);
- localGFEFDStiffness.assembleGradientAndHessian(localView,
+ localGFEFDStiffness.assembleGradientAndHessian(tangentLocalView,
localSolution,
localRiemannianFDGradient,
localRiemannianFDHessian);