diff --git a/dune/gfe/assemblers/CMakeLists.txt b/dune/gfe/assemblers/CMakeLists.txt
index 7ebb988f4f42a1cbbf3b74b920c4d1d514b5222a..935f1c8180271cc826faa03ca3ceca575218fe38 100644
--- a/dune/gfe/assemblers/CMakeLists.txt
+++ b/dune/gfe/assemblers/CMakeLists.txt
@@ -1,6 +1,5 @@
#install headers
install(FILES
- cosseratenergystiffness.hh
cosseratrodenergy.hh
geodesicfeassembler.hh
geodesicfeassemblerwrapper.hh
diff --git a/dune/gfe/assemblers/cosseratenergystiffness.hh b/dune/gfe/assemblers/cosseratenergystiffness.hh
deleted file mode 100644
index 2cfcfdb626c066ce1f5c7084d3e66bed0b8bb276..0000000000000000000000000000000000000000
--- a/dune/gfe/assemblers/cosseratenergystiffness.hh
+++ /dev/null
@@ -1,435 +0,0 @@
-#ifndef COSSERAT_ENERGY_LOCAL_STIFFNESS_HH
-#define COSSERAT_ENERGY_LOCAL_STIFFNESS_HH
-
-#define DONT_USE_CURL
-
-//#define QUADRATIC_2006 //only relevant for gridDim == 3
-// Use the formula for the quadratic membrane energy (first equation of (4.4)) in Neff's paper from 2006: A geometrically exact planar Cosserat shell model with microstructure: Existence of minimizers for zero Cosserat couple modulus
-// instead of the equation (2.27) of 2019: Reļ¬ned dimensional reduction for isotropic elastic Cosserat shells with initial curvature
-
-//#define QUADRATIC_MEMBRANE_ENERGY //only relevant for gridDim == 2
-
-#include <dune/common/fmatrix.hh>
-#include <dune/common/parametertree.hh>
-#include <dune/common/version.hh>
-
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/fufem/boundarypatch.hh>
-
-#include <dune/gfe/assemblers/localenergy.hh>
-#include <dune/gfe/densities/bulkcosseratdensity.hh>
-#include <dune/gfe/densities/planarcosseratshelldensity.hh>
-#ifdef PROJECTED_INTERPOLATION
-#include <dune/gfe/localprojectedfefunction.hh>
-#else
-#include <dune/gfe/localgeodesicfefunction.hh>
-#endif
-#include <dune/gfe/tensor3.hh>
-#include <dune/gfe/cosseratstrain.hh>
-#include <dune/gfe/spaces/orthogonalmatrix.hh>
-#include <dune/gfe/spaces/productmanifold.hh>
-#include <dune/gfe/spaces/realtuple.hh>
-#include <dune/gfe/spaces/rotation.hh>
-
-/** \brief Get LocalFiniteElements from a localView, for different tree depths of the local view
- *
- * We instantiate the CosseratEnergyLocalStiffness class with two different kinds of Basis:
- * A scalar one and a composite one that combines two scalar ones. But code for accessing the
- * finite elements in the basis tree only work for one kind of basis, not for the other.
- * To allow both kinds of basis in a single class we need this trickery below.
- */
-template <class Basis, std::size_t i>
-class LocalFiniteElementFactory
-{
-public:
- static auto get(const typename Basis::LocalView& localView,
- std::integral_constant<std::size_t, i> iType)
- -> decltype(localView.tree().child(iType,0).finiteElement())
- {
- return localView.tree().child(iType,0).finiteElement();
- }
-};
-
-/** \brief Specialize for scalar bases, here we cannot call tree().child() */
-template <class GridView, int order, std::size_t i>
-class LocalFiniteElementFactory<Dune::Functions::LagrangeBasis<GridView,order>,i>
-{
-public:
- static auto get(const typename Dune::Functions::LagrangeBasis<GridView,order>::LocalView& localView,
- std::integral_constant<std::size_t, i> iType)
- -> decltype(localView.tree().finiteElement())
- {
- return localView.tree().finiteElement();
- }
-};
-
-template<class Basis, int dim, class field_type=double>
-class CosseratEnergyLocalStiffness
- : public Dune::GFE::LocalEnergy<Basis,Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > >
-{
- // grid types
- typedef typename Basis::GridView GridView;
- typedef typename GridView::ctype DT;
- typedef Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > TargetSpace;
- typedef typename TargetSpace::ctype RT;
- typedef typename GridView::template Codim<0>::Entity Entity;
-
- // some other sizes
- constexpr static int gridDim = GridView::dimension;
- constexpr static int dimworld = GridView::dimensionworld;
-
-public:
-
- /** \brief Constructor with a set of material parameters
- * \param parameters The material parameters
- */
- CosseratEnergyLocalStiffness(const Dune::ParameterTree& parameters)
- : bulkDensity_(parameters),
- planarCosseratShellDensity_(parameters)
- {}
-
- /** \brief Constructor with material parameters and external loads
- * \param parameters The material parameters
- * \param neumannBoundary The part of the boundary where a surface load is applied
- * \param neumannFunction Surface load density
- * \param volumeLoad Volume load density
- */
- CosseratEnergyLocalStiffness(const Dune::ParameterTree& parameters,
- 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)
- : bulkDensity_(parameters),
- planarCosseratShellDensity_(parameters),
- neumannBoundary_(neumannBoundary),
- neumannFunction_(neumannFunction),
- volumeLoad_(volumeLoad)
- {}
-
- /** \brief Assemble the energy for a single element */
- RT energy (const typename Basis::LocalView& localView,
- const std::vector<TargetSpace>& localSolution) const override;
-
- RT energy (const typename Basis::LocalView& localView,
- const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override;
-
- /** The energy density if the grid is three-dimensional */
- Dune::GFE::BulkCosseratDensity<Dune::FieldVector<double,3>,field_type> bulkDensity_;
-
- /** The energy density if the grid is two-dimensional */
- Dune::GFE::PlanarCosseratShellDensity<Dune::FieldVector<double,2>,field_type> planarCosseratShellDensity_;
-
- /** \brief The Neumann boundary */
- const BoundaryPatch<GridView>* neumannBoundary_ = nullptr;
-
- /** \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>
-[[deprecated("Use an std::vector<RealTuple<field_type,dim>> and an std::vector<Rotation<field_type,dim>> together with the MixedGFEAssembler and the GFEAssemblerWrapper instead of std::vector<Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> >>.")]]
-typename CosseratEnergyLocalStiffness<Basis,dim,field_type>::RT
-CosseratEnergyLocalStiffness<Basis,dim,field_type>::
-energy(const typename Basis::LocalView& localView,
- const std::vector<Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > >& localSolution) const
-{
- using namespace Dune::Indices;
-
- RT energy = 0;
-
- auto element = localView.element();
-
- using namespace Dune::Indices;
- const auto& localFiniteElement = LocalFiniteElementFactory<Basis,0>::get(localView,_0);
-#ifdef PROJECTED_INTERPOLATION
- typedef Dune::GFE::LocalProjectedFEFunction<gridDim, DT, decltype(localFiniteElement), TargetSpace> LocalGFEFunctionType;
-#else
- typedef LocalGeodesicFEFunction<gridDim, DT, decltype(localFiniteElement), TargetSpace> LocalGFEFunctionType;
-#endif
- LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution);
-
- int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
- : localFiniteElement.localBasis().order() * gridDim;
-
- const auto& quad = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
-
- for (size_t pt=0; pt<quad.size(); pt++) {
-
- // Local position of the quadrature point
- const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
-
- const DT integrationElement = element.geometry().integrationElement(quadPos);
-
- const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
-
- DT weight = quad[pt].weight() * integrationElement;
-
- // The value of the local function
- Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > value = localGeodesicFEFunction.evaluate(quadPos);
-
- // The derivative of the local function defined on the reference element
- typename LocalGFEFunctionType::DerivativeType referenceDerivative = localGeodesicFEFunction.evaluateDerivative(quadPos,value);
-
- // The derivative of the function defined on the actual element
- typename LocalGFEFunctionType::DerivativeType derivative(0);
-
- for (size_t comp=0; comp<referenceDerivative.N(); comp++)
- jacobianInverseTransposed.umv(referenceDerivative[comp], derivative[comp]);
-
- /////////////////////////////////////////////////////////
- // compute U, the Cosserat strain
- /////////////////////////////////////////////////////////
- static_assert(dim>=gridDim, "Codim of the grid must be nonnegative");
-
- //
- Dune::FieldMatrix<field_type,dim,dim> R;
- value[_1].matrix(R);
-
- Dune::GFE::CosseratStrain<field_type,dim,gridDim> U(derivative,R);
-
- //////////////////////////////////////////////////////////
- // Compute the derivative of the rotation
- // Note: we need it in matrix coordinates
- //////////////////////////////////////////////////////////
-
- // Extract the orientation derivative
- Dune::FieldMatrix<field_type,4,gridDim> orientationDerivative;
- for (size_t i=0; i<4; ++i)
- for (size_t j=0; j<gridDim; ++j)
- orientationDerivative[i][j] = derivative[3+i][j];
-
- Tensor3<field_type,3,3,gridDim> DR = value[_1].quaternionTangentToMatrixTangent(orientationDerivative);
-
- // Add the local energy density
- if constexpr (gridDim==2) {
- energy += weight * planarCosseratShellDensity_(quadPos,
- value.globalCoordinates(),
- derivative);
- } else if constexpr (gridDim==3) {
- energy += weight * bulkDensity_(quadPos,
- value.globalCoordinates(),
- derivative);
- } else
- DUNE_THROW(Dune::NotImplemented, "CosseratEnergyStiffness for 1d grids");
-
- ///////////////////////////////////////////////////////////
- // Volume load contribution
- ///////////////////////////////////////////////////////////
-
- if (not volumeLoad_)
- continue;
-
- // Value of the volume load density at the current position
- auto volumeLoadDensity = volumeLoad_(element.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 Dune::QuadratureRule<DT, gridDim-1>& 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
- auto 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;
-}
-
-template <class Basis, int dim, class field_type>
-typename CosseratEnergyLocalStiffness<Basis,dim,field_type>::RT
-CosseratEnergyLocalStiffness<Basis,dim,field_type>::
-energy(const typename Basis::LocalView& localView,
- const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& localConfiguration) const
-{
- auto element = localView.element();
-
- RT energy = 0;
-
- using namespace Dune::Indices;
- const auto& deformationLocalFiniteElement = LocalFiniteElementFactory<Basis,0>::get(localView,_0);
- const auto& orientationLocalFiniteElement = LocalFiniteElementFactory<Basis,1>::get(localView,_1);
-
-#ifdef PROJECTED_INTERPOLATION
- typedef Dune::GFE::LocalProjectedFEFunction<gridDim, DT, decltype(deformationLocalFiniteElement), RealTuple<field_type,dim> >
- LocalDeformationGFEFunctionType;
-#else
- typedef LocalGeodesicFEFunction<gridDim, DT, decltype(deformationLocalFiniteElement), RealTuple<field_type,dim> >
- LocalDeformationGFEFunctionType;
-#endif
- LocalDeformationGFEFunctionType localDeformationGFEFunction(deformationLocalFiniteElement,localConfiguration[_0]);
-
-#ifdef PROJECTED_INTERPOLATION
- typedef Dune::GFE::LocalProjectedFEFunction<gridDim, DT, decltype(orientationLocalFiniteElement), Rotation<field_type,dim> > LocalOrientationGFEFunctionType;
-#else
- typedef LocalGeodesicFEFunction<gridDim, DT, decltype(orientationLocalFiniteElement), Rotation<field_type,dim> > LocalOrientationGFEFunctionType;
-#endif
- LocalOrientationGFEFunctionType localOrientationGFEFunction(orientationLocalFiniteElement,localConfiguration[_1]);
-
- // \todo Implement smarter quadrature rule selection for more efficiency, i.e., less evaluations of the Rotation GFE function
- int quadOrder = deformationLocalFiniteElement.localBasis().order() * ((element.type().isSimplex()) ? 1 : gridDim);
-
- const auto& quad = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
-
- for (size_t pt=0; pt<quad.size(); pt++)
- {
- // Local position of the quadrature point
- const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
-
- const DT integrationElement = element.geometry().integrationElement(quadPos);
-
- const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
-
- DT weight = quad[pt].weight() * integrationElement;
-
- // The value of the local deformation
- RealTuple<field_type,dim> deformationValue = localDeformationGFEFunction.evaluate(quadPos);
- Rotation<field_type,dim> orientationValue = localOrientationGFEFunction.evaluate(quadPos);
-
- TargetSpace value;
- value[_0] = deformationValue;
- value[_1] = orientationValue;
-
- // The derivative of the local function defined on the reference element
- typename LocalDeformationGFEFunctionType::DerivativeType deformationReferenceDerivative = localDeformationGFEFunction.evaluateDerivative(quadPos,deformationValue);
- typename LocalOrientationGFEFunctionType::DerivativeType orientationReferenceDerivative = localOrientationGFEFunction.evaluateDerivative(quadPos,orientationValue);
-
- // The derivative of the function defined on the actual element
- typename LocalDeformationGFEFunctionType::DerivativeType deformationDerivative;
- typename LocalOrientationGFEFunctionType::DerivativeType orientationDerivative;
-
- for (size_t comp=0; comp<deformationReferenceDerivative.N(); comp++)
- jacobianInverseTransposed.mv(deformationReferenceDerivative[comp], deformationDerivative[comp]);
-
- for (size_t comp=0; comp<orientationReferenceDerivative.N(); comp++)
- jacobianInverseTransposed.mv(orientationReferenceDerivative[comp], orientationDerivative[comp]);
-
- // The derivative of both factors together
- typename Dune::FieldMatrix<field_type,deformationDerivative.N() + orientationDerivative.N(), gridDim> derivative(0);
-
- size_t rowCount = 0;
- for (auto&& row : deformationDerivative)
- derivative[rowCount++] = row;
-
- for (auto&& row : orientationDerivative)
- derivative[rowCount++] = row;
-
- /////////////////////////////////////////////////////////
- // compute U, the Cosserat strain
- /////////////////////////////////////////////////////////
- static_assert(dim>=gridDim, "Codim of the grid must be nonnegative");
-
- //
- Dune::FieldMatrix<field_type,dim,dim> R;
- orientationValue.matrix(R);
-
- Dune::GFE::CosseratStrain<field_type,dim,gridDim> U(deformationDerivative,R);
-
- //////////////////////////////////////////////////////////
- // Compute the derivative of the rotation
- // Note: we need it in matrix coordinates
- //////////////////////////////////////////////////////////
-
- Tensor3<field_type,3,3,gridDim> DR = orientationValue.quaternionTangentToMatrixTangent(orientationDerivative);
-
- // Add the local energy density
- if constexpr (gridDim==2) {
- energy += weight * planarCosseratShellDensity_(quadPos,
- value.globalCoordinates(),
- derivative);
- } else if constexpr (gridDim==3) {
- energy += weight * bulkDensity_(quadPos,
- value.globalCoordinates(),
- derivative);
- } else
- DUNE_THROW(Dune::NotImplemented, "CosseratEnergyStiffness for 1d grids");
-
- ///////////////////////////////////////////////////////////
- // Volume load contribution
- ///////////////////////////////////////////////////////////
- if (not volumeLoad_)
- continue;
-
- // Value of the volume load density at the current position
- auto volumeLoadDensity = volumeLoad_(element.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
- auto 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.globalCoordinates()[i]) * quad[pt].weight() * integrationElement;
-
- }
-
- }
-
- return energy;
-}
-
-#endif //#ifndef COSSERAT_ENERGY_LOCAL_STIFFNESS_HH
diff --git a/dune/gfe/assemblers/localintegralenergy.hh b/dune/gfe/assemblers/localintegralenergy.hh
index b10d519db2b8e18660ed5431e16553d775223676..2dd2d5ed06264e4021d9ea9cfdbc87a4145fddb0 100644
--- a/dune/gfe/assemblers/localintegralenergy.hh
+++ b/dune/gfe/assemblers/localintegralenergy.hh
@@ -158,7 +158,9 @@ namespace Dune::GFE {
{
RT energy = 0;
- if constexpr (Impl::LocalEnergyTypes<TargetSpace>::isProductManifold && Basis::LocalView::Tree::isComposite)
+ if constexpr (Impl::LocalEnergyTypes<TargetSpace>::isProductManifold
+ && Basis::LocalView::Tree::isComposite
+ && gridDim==GridView::dimensionworld) // TODO: Implement the case gridDim!=dimworld
{
static_assert(TargetSpace::size() == 2,
"LocalIntegralEnergy only implemented for product spaces with two factors!");
@@ -230,7 +232,7 @@ namespace Dune::GFE {
}
}
else
- DUNE_THROW(Dune::NotImplemented, "Non-product manifold or non-composite basis");
+ DUNE_THROW(Dune::NotImplemented, "Non-product manifold or non-composite basis or gridDim!=dimworld");
return energy;
}
diff --git a/dune/gfe/neumannenergy.hh b/dune/gfe/neumannenergy.hh
index 7f408836cb027e8ab0c3c82fe10dd4a3f71b6780..de6c78a6b8f6adc01ef2676fadadfb38ba7dcfb7 100644
--- a/dune/gfe/neumannenergy.hh
+++ b/dune/gfe/neumannenergy.hh
@@ -27,6 +27,7 @@ namespace Dune::GFE {
using RT = typename Dune::GFE::LocalEnergy<Basis,TargetSpace>::RT;
constexpr static int dim = GridView::dimension;
+ constexpr static int dimworld = GridView::dimensionworld;
// TODO: Remove the hard-coded first factor space!
using WorldVector = typename std::tuple_element_t<0, std::tuple<TargetSpaces...> >::EmbeddedTangentVector;
@@ -37,7 +38,7 @@ namespace Dune::GFE {
* \param parameters The material parameters
*/
NeumannEnergy(const std::shared_ptr<BoundaryPatch<GridView> >& neumannBoundary,
- std::function<WorldVector(Dune::FieldVector<DT,dim>)> neumannFunction)
+ std::function<WorldVector(Dune::FieldVector<DT,dimworld>)> neumannFunction)
: neumannBoundary_(neumannBoundary),
neumannFunction_(neumannFunction)
{}
@@ -107,7 +108,7 @@ namespace Dune::GFE {
const std::shared_ptr<BoundaryPatch<GridView> > neumannBoundary_;
/** \brief The function implementing the Neumann data */
- std::function<WorldVector(Dune::FieldVector<DT,dim>)> neumannFunction_;
+ std::function<WorldVector(Dune::FieldVector<DT,dimworld>)> neumannFunction_;
};
} // namespace Dune::GFE
diff --git a/problems/cosserat-continuum-cantilever.py b/problems/cosserat-continuum-cantilever.py
index 0dfec2fed041e962900575ddb87ceda2bb36b943..47f68a9dee8df54b2f0bc5a889687827014b51f7 100644
--- a/problems/cosserat-continuum-cantilever.py
+++ b/problems/cosserat-continuum-cantilever.py
@@ -95,12 +95,6 @@ parameterSet.materialParameters.q = 2
# Shear correction factor
parameterSet.materialParameters.kappa = 1
-# TODO: These three parameters are not actually used,
-# but the current implementation requires them nevertheless.
-parameterSet.materialParameters.b1 = 1
-parameterSet.materialParameters.b2 = 1
-parameterSet.materialParameters.b3 = 1
-
#############################################
# Boundary values
diff --git a/problems/cosserat-continuum-twisted-strip.py b/problems/cosserat-continuum-twisted-strip.py
index 8439f111db0bd471dd25b25ddcb9f4fe6bf09ba1..4795e50f41fd698ad41c278d27058aaeba9b539b 100644
--- a/problems/cosserat-continuum-twisted-strip.py
+++ b/problems/cosserat-continuum-twisted-strip.py
@@ -92,12 +92,6 @@ parameterSet.materialParameters.q = 2
# Shear correction factor
parameterSet.materialParameters.kappa = 1
-# TODO: These three parameters are not actually used,
-# but the current implementation requires them nevertheless.
-parameterSet.materialParameters.b1 = 1
-parameterSet.materialParameters.b2 = 1
-parameterSet.materialParameters.b3 = 1
-
#############################################
# Boundary values
diff --git a/problems/cosserat-continuum-wong-pellegrino-wrinkling.py b/problems/cosserat-continuum-wong-pellegrino-wrinkling.py
index 8bce680b32d804e01c458fc767beccdaa7f377bb..83f1d60b4f5c4c6b3ef935661593b5cf2cff54bf 100644
--- a/problems/cosserat-continuum-wong-pellegrino-wrinkling.py
+++ b/problems/cosserat-continuum-wong-pellegrino-wrinkling.py
@@ -91,12 +91,6 @@ parameterSet.materialParameters.q = 2
# Shear correction factor
parameterSet.materialParameters.kappa = 1
-# TODO: These three parameters are not actually used,
-# but the current implementation requires them nevertheless.
-parameterSet.materialParameters.b1 = 1
-parameterSet.materialParameters.b2 = 1
-parameterSet.materialParameters.b3 = 1
-
#############################################
# Boundary values
diff --git a/problems/cosserat-continuum-wriggers-l-shape.py b/problems/cosserat-continuum-wriggers-l-shape.py
index edd8a7065a846dc6a5239882763e7e062c32cbe9..b91a935213e8fd7bd3e5c104626c87d08ea43161 100644
--- a/problems/cosserat-continuum-wriggers-l-shape.py
+++ b/problems/cosserat-continuum-wriggers-l-shape.py
@@ -87,12 +87,6 @@ parameterSet.materialParameters.q = 2
# Shear correction factor
parameterSet.materialParameters.kappa = 1
-# TODO: These three parameters are not actually used,
-# but the current implementation requires them nevertheless.
-parameterSet.materialParameters.b1 = 1
-parameterSet.materialParameters.b2 = 1
-parameterSet.materialParameters.b3 = 1
-
#############################################
# Boundary values
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index cdfe2623f36ea932d8699be6eb0dd6b9227f5163..51d8ce7be8dd48aec2cd27ba31efaca1a80c11e4 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -58,14 +58,18 @@
#include <dune/gfe/localgeodesicfefunction.hh>
#include <dune/gfe/localprojectedfefunction.hh>
+#include <dune/gfe/neumannenergy.hh>
#include <dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh>
-#include <dune/gfe/assemblers/cosseratenergystiffness.hh>
+#include <dune/gfe/assemblers/localintegralenergy.hh>
#include <dune/gfe/assemblers/nonplanarcosseratshellenergy.hh>
#include <dune/gfe/cosseratvtkwriter.hh>
#include <dune/gfe/cosseratvtkreader.hh>
#include <dune/gfe/assemblers/geodesicfeassembler.hh>
#include <dune/gfe/embeddedglobalgfefunction.hh>
#include <dune/gfe/assemblers/mixedgfeassembler.hh>
+#include <dune/gfe/assemblers/sumenergy.hh>
+#include <dune/gfe/densities/bulkcosseratdensity.hh>
+#include <dune/gfe/densities/planarcosseratshelldensity.hh>
#if MIXED_SPACE
#include <dune/gfe/mixedriemannianpnsolver.hh>
@@ -105,6 +109,21 @@ static_assert(displacementOrder==rotationOrder, "displacement and rotation order
// Image space of the geodesic fe functions
using TargetSpace = GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> >;
+// Method to construct a density that matches the grid dimension.
+// This cannot be done inside the 'main' method, because 'constexpr if' only works
+// when its argument depends on a template parameter.
+template <typename LocalCoordinate>
+auto createDensity(const ParameterTree& materialParameters)
+{
+ if constexpr (LocalCoordinate::size()==2)
+ {
+ return std::make_shared<GFE::PlanarCosseratShellDensity<LocalCoordinate, adouble> >(materialParameters);
+ }
+ else
+ {
+ return std::make_shared<GFE::BulkCosseratDensity<LocalCoordinate, adouble> >(materialParameters);
+ }
+}
int main (int argc, char *argv[]) try
{
@@ -322,16 +341,16 @@ int main (int argc, char *argv[]) try
neumannVertices[indexSet.index(vertex)] = isNeumann;
}
- BoundaryPatch<GridView> neumannBoundary(gridView, neumannVertices);
+ auto neumannBoundary = std::make_shared<BoundaryPatch<GridView> >(gridView, neumannVertices);
- std::cout << "On rank " << mpiHelper.rank() << ": Neumann boundary has " << neumannBoundary.numFaces()
+ std::cout << "On rank " << mpiHelper.rank() << ": Neumann boundary has " << neumannBoundary->numFaces()
<< " faces and " << neumannVertices.count() << " degrees of freedom.\n";
BitSetVector<1> neumannNodes(deformationFEBasis.size(), false);
#if DUNE_VERSION_GTE(DUNE_FUFEM, 2, 10)
- Fufem::markBoundaryPatchDofs(neumannBoundary,deformationFEBasis,neumannNodes);
+ Fufem::markBoundaryPatchDofs(*neumannBoundary,deformationFEBasis,neumannNodes);
#else
- constructBoundaryDofs(neumannBoundary,deformationFEBasis,neumannNodes);
+ constructBoundaryDofs(*neumannBoundary,deformationFEBasis,neumannNodes);
#endif
for (size_t i=0; i<deformationFEBasis.size(); i++) {
@@ -485,7 +504,14 @@ int main (int argc, char *argv[]) try
return nV;
};
+ if (parameterSet.get<std::string>("volumeLoadPythonFunction", "zero-volume-load") != "zero-volume-load")
+ {
+ std::cerr << "cosserat-continuum.cc: Volume loads are not fully implemented yet." << std::endl;
+ std::abort();
+ }
+
Python::Reference volumeLoadClass = Python::import(parameterSet.get<std::string>("volumeLoadPythonFunction", "zero-volume-load"));
+
Python::Callable volumeLoadCallable = volumeLoadClass.get("VolumeLoad");
// Call a constructor
@@ -532,14 +558,37 @@ int main (int argc, char *argv[]) try
x[_1][i].set(dOV[i]);
if (dim==dimworld) {
- auto localCosseratEnergy = std::make_shared<CosseratEnergyLocalStiffness<CompositeBasis, 3,adouble> > (materialParameters,
- &neumannBoundary,
- neumannFunction,
- volumeLoad);
+ auto sumEnergy = std::make_shared<GFE::SumEnergy<CompositeBasis, RealTuple<adouble,3>,Rotation<adouble,3> > >();
- LocalGeodesicFEADOLCStiffness<CompositeBasis,TargetSpace> localGFEADOLCStiffness(localCosseratEnergy,
+ // 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 ATargetSpace = typename TargetSpace::rebind<adouble>::other;
+
+ using LocalCoordinate = typename GridType::Codim<0>::Entity::Geometry::LocalCoordinate;
+ auto cosseratDensity = createDensity<LocalCoordinate>(materialParameters);
+
+ auto localCosseratEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis,AInterpolationRule,ATargetSpace> >(cosseratDensity);
+
+ 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 local assembler
+ LocalGeodesicFEADOLCStiffness<CompositeBasis,TargetSpace> localGFEADOLCStiffness(sumEnergy,
adolcScalarMode);
+
MixedGFEAssembler<CompositeBasis,TargetSpace> mixedAssembler(compositeBasis, localGFEADOLCStiffness);
+
+ ////////////////////////////////////////////
+ // Set up the solver
+ ////////////////////////////////////////////
#if MIXED_SPACE
if (parameterSet.get<std::string>("solvertype", "trustRegion") == "trustRegion")
{
@@ -665,7 +714,7 @@ int main (int argc, char *argv[]) try
#if HAVE_DUNE_CURVEDGEOMETRY && WORLD_DIM == 3 && GRID_DIM == 2
auto localCosseratEnergy = std::make_shared<NonplanarCosseratShellEnergy<CompositeBasis, 3, adouble, decltype(creator)> >(materialParameters,
&creator,
- &neumannBoundary,
+ neumannBoundary.get(),
neumannFunction,
volumeLoad);