diff --git a/dune/gfe/sumcosseratenergy.hh b/dune/gfe/sumcosseratenergy.hh
deleted file mode 100644
index 51f6241a2090e51977d9fc79f0f426bf2e49911f..0000000000000000000000000000000000000000
--- a/dune/gfe/sumcosseratenergy.hh
+++ /dev/null
@@ -1,73 +0,0 @@
-#ifndef DUNE_GFE_SUMCOSSERATENERGY_HH
-#define DUNE_GFE_SUMCOSSERATENERGY_HH
-
-#include <dune/elasticity/assemblers/localfestiffness.hh>
-
-#include <dune/gfe/localenergy.hh>
-
-namespace Dune {
-
-namespace GFE {
-template<class Basis, class TargetSpace, class field_type=double>
-class SumCosseratEnergy
-: public GFE::LocalEnergy<Basis, TargetSpace>
-{
- // grid types
- typedef typename Basis::GridView GridView;
- typedef typename GridView::ctype ctype;
- typedef typename Basis::LocalView::Tree::FiniteElement LocalFiniteElement;
- typedef typename GridView::ctype DT;
- typedef typename TargetSpace::ctype RT;
-
- enum {dim=GridView::dimension};
-
-public:
-
- /** \brief Constructor with elastic energy and cosserat energy
- * \param elasticEnergy The elastic energy
- * \param cosseratEnergy The cosserat energy
- */
-#if DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 7)
- SumCosseratEnergy(std::shared_ptr<LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy,
-#elif DUNE_VERSION_GTE(DUNE_ELASTICITY, 2, 8)
- SumCosseratEnergy(std::shared_ptr<Dune::LocalEnergy<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy,
-#else
- SumCosseratEnergy(std::shared_ptr<Elasticity::LocalEnergy<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy,
-#endif
- std::shared_ptr<GFE::LocalEnergy<Basis, TargetSpace>> cosseratEnergy)
-
- : elasticEnergy_(elasticEnergy),
- cosseratEnergy_(cosseratEnergy)
- {}
-
- /** \brief Assemble the energy for a single element */
- RT energy (const typename Basis::LocalView& localView,
- const std::vector<TargetSpace>& localSolution) const
- {
- auto&& element = localView.element();
- auto&& localFiniteElement = localView.tree().finiteElement();
- std::vector<Dune::FieldVector<field_type,dim>> localConfiguration(localSolution.size());
- for (int i = 0; i < localSolution.size(); i++) {
- localConfiguration[i] = localSolution[i].r;
- }
- return elasticEnergy_->energy(element, localFiniteElement, localConfiguration) + cosseratEnergy_->energy(localView, localSolution);
- }
-
-private:
-
-#if DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 7)
- std::shared_ptr<LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy_;
-#elif DUNE_VERSION_GTE(DUNE_ELASTICITY, 2, 8)
- std::shared_ptr<Dune::LocalEnergy<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy_;
-#else
- std::shared_ptr<Elasticity::LocalEnergy<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy_;
-#endif
-
- std::shared_ptr<GFE::LocalEnergy<Basis, TargetSpace> > cosseratEnergy_;
-};
-
-} // namespace GFE
-
-} // namespace Dune
-
-#endif //#ifndef DUNE_GFE_SUMCOSSERATENERGY_HH
diff --git a/dune/gfe/surfacecosseratenergy.hh b/dune/gfe/surfacecosseratenergy.hh
index d6ddb215820fb89b096cc71ba5f096e2e675b82f..366e233747f0a98c2e986b40fb8c606be0963e4f 100644
--- a/dune/gfe/surfacecosseratenergy.hh
+++ b/dune/gfe/surfacecosseratenergy.hh
@@ -1,6 +1,7 @@
#ifndef DUNE_GFE_SURFACECOSSERATENERGY_HH
#define DUNE_GFE_SURFACECOSSERATENERGY_HH
+#include <dune/common/indices.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/fufem/functions/virtualgridfunction.hh>
@@ -16,14 +17,15 @@
#include <dune/gfe/tensor3.hh>
#include <dune/gfe/vertexnormals.hh>
+namespace Dune::GFE {
-template<class Basis, class TargetSpace, class field_type=double, class GradientRT=double>
+template<class Basis, class... TargetSpaces>
class SurfaceCosseratEnergy
-: public Dune::GFE::LocalEnergy<Basis,TargetSpace>
+: public Dune::GFE::LocalEnergy<Basis, TargetSpaces...>
{
using GridView = typename Basis::GridView;
using DT = typename GridView::ctype ;
- using RT = typename TargetSpace::ctype;
+ using RT = typename Dune::GFE::LocalEnergy<Basis, TargetSpaces...>::RT ;
using Entity = typename GridView::template Codim<0>::Entity ;
using RBM0 = RealTuple<RT,GridView::dimensionworld> ;
using RBM1 = Rotation<RT,GridView::dimensionworld> ;
@@ -94,13 +96,18 @@ public:
}
RT energy(const typename Basis::LocalView& localView,
- const std::vector<RigidBodyMotion<field_type,dimWorld> >& localSolution) const
-{
+ const std::vector<TargetSpaces>&... localSolutions) const
+{
+ static_assert(sizeof...(TargetSpaces) == 2, "SurfaceCosseratEnergy needs exactly two TargetSpace!");
+
+ using namespace Dune::Indices;
+ using TargetSpace0 = typename std::tuple_element<0, std::tuple<TargetSpaces...> >::type;
+ const std::vector<TargetSpace0>& localSolution0 = std::get<0>(std::forward_as_tuple(localSolutions...));
+ using TargetSpace1 = typename std::tuple_element<1, std::tuple<TargetSpaces...> >::type;
+ const std::vector<TargetSpace1>& localSolution1 = std::get<1>(std::forward_as_tuple(localSolutions...));
+
// The element geometry
auto element = localView.element();
-
- // The set of shape functions on this element
- const auto& localFiniteElement = localView.tree().finiteElement();
auto gridView = localView.globalBasis().gridView();
////////////////////////////////////////////////////////////////////////////////////
@@ -121,8 +128,34 @@ RT energy(const typename Basis::LocalView& localView,
////////////////////////////////////////////////////////////////////////////////////
// Set up the local nonlinear finite element function
////////////////////////////////////////////////////////////////////////////////////
- typedef LocalGeodesicFEFunction<gridDim, DT, typename Basis::LocalView::Tree::FiniteElement, TargetSpace> LocalGFEFunctionType;
- LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution);
+ using namespace Dune::TypeTree::Indices;
+
+ // The set of shape functions on this element
+ const auto& deformationLocalFiniteElement = localView.tree().child(_0,0).finiteElement();
+ const auto& orientationLocalFiniteElement = localView.tree().child(_1,0).finiteElement();
+
+ // to construt a local GFE function, in case they are the shape functions are the same, we can use use one GFE-Function
+#if MIXED_SPACE
+ std::vector<RBM0> localSolutionRBM0(localSolution0.size());
+ std::vector<RBM1> localSolutionRBM1(localSolution1.size());
+ for (int i = 0; i < localSolution0.size(); i++)
+ localSolutionRBM0[i] = localSolution0[i];
+ for (int i = 0; i< localSolution1.size(); i++)
+ localSolutionRBM1[i] = localSolution1[i];
+ typedef LocalGeodesicFEFunction<gridDim, DT, decltype(deformationLocalFiniteElement), RBM0> LocalGFEFunctionType0;
+ typedef LocalGeodesicFEFunction<gridDim, DT, decltype(orientationLocalFiniteElement), RBM1> LocalGFEFunctionType1;
+ LocalGFEFunctionType0 localGeodesicFEFunction0(deformationLocalFiniteElement,localSolutionRBM0);
+ LocalGFEFunctionType1 localGeodesicFEFunction1(orientationLocalFiniteElement,localSolutionRBM1);
+#else
+ std::vector<RBM> localSolutionRBM(localSolution0.size());
+ for (int i = 0; i < localSolution0.size(); i++) {
+ for (int j = 0; j < dimWorld; j++)
+ localSolutionRBM[i].r[j] = localSolution0[i][j];
+ localSolutionRBM[i].q = localSolution1[i];
+ }
+ typedef LocalGeodesicFEFunction<gridDim, DT, decltype(deformationLocalFiniteElement), RBM> LocalGFEFunctionType;
+ LocalGFEFunctionType localGeodesicFEFunction(deformationLocalFiniteElement,localSolutionRBM);
+#endif
RT energy = 0;
@@ -134,8 +167,8 @@ RT energy(const typename Basis::LocalView& localView,
auto id = idSet.subId(it.inside(), it.indexInInside(), 1);
auto boundaryGeometry = geometriesOnShellBoundary_.at(id);
- auto quadOrder = (it.type().isSimplex()) ? localFiniteElement.localBasis().order()
- : localFiniteElement.localBasis().order() * boundaryDim;
+ auto quadOrder = (it.type().isSimplex()) ? deformationLocalFiniteElement.localBasis().order()
+ : deformationLocalFiniteElement.localBasis().order() * boundaryDim;
const auto& quad = Dune::QuadratureRules<DT, boundaryDim>::rule(it.type(), quadOrder);
for (size_t pt=0; pt<quad.size(); pt++) {
@@ -152,18 +185,42 @@ RT energy(const typename Basis::LocalView& localView,
const DT integrationElement = boundaryGeometry.integrationElement(quad[pt].position());
- // The value of the local function
- RigidBodyMotion<field_type,dimWorld> value = localGeodesicFEFunction.evaluate(quadPos);
-
- // The derivative of the local function
- auto derivative = localGeodesicFEFunction.evaluateDerivative(quadPos,value);
+ RBM value;
+ Dune::FieldMatrix<RT, RBM::embeddedDim, dimWorld> derivative;
+ Dune::FieldMatrix<RT, RBM::embeddedDim, boundaryDim> derivative2D;
- Dune::FieldMatrix<RT,7,2> derivative2D;
- for (int i = 0; i < 7; i++) {
- for (int j = 0; j < 2; j++) {
- derivative2D[i][j] = derivative[i][j];
+#if MIXED_SPACE
+ // The value of the local function
+ RBM0 value0 = localGeodesicFEFunction0.evaluate(quadPos);
+ RBM1 value1 = localGeodesicFEFunction1.evaluate(quadPos);
+ // The derivatives of the local functions
+ auto derivative0 = localGeodesicFEFunction0.evaluateDerivative(quadPos,value0);
+ auto derivative1 = localGeodesicFEFunction1.evaluateDerivative(quadPos,value1);
+
+ // Put the value and the derivatives together from the separated values
+ for (int i = 0; i < RBM0::dim; i++)
+ value.r[i] = value0[i];
+ value.q = value1;
+ for (int j = 0; j < dimWorld; j++) {
+ for (int i = 0; i < RBM0::embeddedDim; i++) {
+ derivative[i][j] = derivative0[i][j];
+ if (j < boundaryDim)
+ derivative2D[i][j] = derivative0[i][j];
+ }
+ for (int i = 0; i < RBM1::embeddedDim; i++) {
+ derivative[RBM0::embeddedDim + i][j] = derivative1[i][j];
+ if (j < boundaryDim)
+ derivative2D[RBM0::embeddedDim + i][j] = derivative1[i][j];
}
}
+#else
+ value = localGeodesicFEFunction.evaluate(quadPos);
+ derivative = localGeodesicFEFunction.evaluateDerivative(quadPos,value);
+ for (int i = 0; i < RBM::embeddedDim; i++)
+ for (int j = 0; j < boundaryDim; j++)
+ derivative2D[i][j] = derivative[i][j];
+#endif
+
//////////////////////////////////////////////////////////
// The rotation and its derivative
// Note: we need it in matrix coordinates
@@ -343,5 +400,6 @@ private:
double b1_, b2_, b3_;
};
+} // namespace Dune::GFE
#endif //#ifndef DUNE_GFE_SURFACECOSSERATENERGY_HH
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index f0843d81f4e46dca10975d882c9b15ff7e45390a..94c64915ad8ecb0540c32699bf78fd3a569a0f98 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -1,3 +1,4 @@
+#define MIXED_SPACE 0
#include <iostream>
#include <fstream>
@@ -12,59 +13,52 @@
#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
#include <dune/common/bitsetvector.hh>
+#include <dune/common/indices.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/common/timer.hh>
+#include <dune/common/tuplevector.hh>
#include <dune/common/version.hh>
-#include <dune/grid/uggrid.hh>
-#include <dune/grid/utility/structuredgridfactory.hh>
-
-#include <dune/grid/io/file/gmshreader.hh>
-#include <dune/grid/io/file/vtk.hh>
-
-#if DUNE_VERSION_GTE(DUNE_ELASTICITY, 2, 8)
#include <dune/elasticity/materials/exphenckydensity.hh>
#include <dune/elasticity/materials/henckydensity.hh>
#include <dune/elasticity/materials/mooneyrivlindensity.hh>
#include <dune/elasticity/materials/neohookedensity.hh>
#include <dune/elasticity/materials/stvenantkirchhoffdensity.hh>
-#include <dune/elasticity/materials/sumenergy.hh>
-#include <dune/elasticity/materials/localintegralenergy.hh>
-#include <dune/elasticity/materials/neumannenergy.hh>
-#else
-#include <dune/elasticity/materials/exphenckyenergy.hh>
-#include <dune/elasticity/materials/henckyenergy.hh>
-#if !DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 7)
-#include <dune/elasticity/materials/mooneyrivlinenergy.hh>
-#endif
-#include <dune/elasticity/materials/neohookeenergy.hh>
-#include <dune/elasticity/materials/neumannenergy.hh>
-#include <dune/elasticity/materials/sumenergy.hh>
-#include <dune/elasticity/materials/stvenantkirchhoffenergy.hh>
-#endif
-
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/functions/functionspacebases/compositebasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
-#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functiontools/boundarydofs.hh>
-#include <dune/fufem/functiontools/basisinterpolator.hh>
-#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
#include <dune/fufem/dunepython.hh>
-#include <dune/gfe/rigidbodymotion.hh>
-#include <dune/gfe/localgeodesicfeadolcstiffness.hh>
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+#include <dune/grid/io/file/gmshreader.hh>
+#include <dune/grid/io/file/vtk.hh>
+
#include <dune/gfe/cosseratvtkwriter.hh>
-#include <dune/gfe/geodesicfeassembler.hh>
+#include <dune/gfe/localintegralenergy.hh>
+#include <dune/gfe/mixedgfeassembler.hh>
+#include <dune/gfe/mixedlocalgfeadolcstiffness.hh>
+#include <dune/gfe/neumannenergy.hh>
+#include <dune/gfe/surfacecosseratenergy.hh>
+#include <dune/gfe/sumenergy.hh>
+#include <dune/gfe/vertexnormals.hh>
+
+#if MIXED_SPACE
+#include <dune/gfe/mixedriemanniantrsolver.hh>
+#else
+#include <dune/gfe/geodesicfeassemblerwrapper.hh>
#include <dune/gfe/riemannianpnsolver.hh>
#include <dune/gfe/riemanniantrsolver.hh>
-#include <dune/gfe/vertexnormals.hh>
-#include <dune/gfe/surfacecosseratenergy.hh>
-#include <dune/gfe/sumcosseratenergy.hh>
+#include <dune/gfe/rigidbodymotion.hh>
+#endif
+
+#include <dune/istl/multitypeblockvector.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
#include <dune/solvers/norms/energynorm.hh>
@@ -74,7 +68,15 @@
# define WORLD_DIM 3
#endif
const int dim = WORLD_DIM;
-const int order = 2;
+
+const int targetDim = WORLD_DIM;
+
+const int displacementOrder = 2;
+const int rotationOrder = 2;
+
+#if !MIXED_SPACE
+static_assert(displacementOrder==rotationOrder, "displacement and rotation order do not match!");
+#endif
#if DUNE_VERSION_LT(DUNE_COMMON, 2, 7)
template<>
@@ -100,8 +102,14 @@ int main (int argc, char *argv[]) try
// initialize MPI, finalize is done automatically on exit
Dune::MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
- if (mpiHelper.rank()==0)
- std::cout << "ORDER = " << order << std::endl;
+ if (mpiHelper.rank()==0) {
+ std::cout << "DISPLACEMENTORDER = " << displacementOrder << ", ROTATIONORDER = " << rotationOrder << std::endl;
+#if MIXED_SPACE
+ std::cout << "MIXED_SPACE = 1" << std::endl;
+#else
+ std::cout << "MIXED_SPACE = 0" << std::endl;
+#endif
+ }
// Start Python interpreter
Python::start();
@@ -115,12 +123,6 @@ int main (int argc, char *argv[]) try
<< std::endl << "sys.path.append(os.getcwd() + '/../../problems/')"
<< std::endl;
- typedef BlockVector<FieldVector<double,dim> > SolutionType;
- typedef RigidBodyMotion<double,dim> TargetSpace;
- typedef std::vector<TargetSpace> SolutionTypeCosserat;
-
- const int blocksize = TargetSpace::TangentVector::dimension;
-
// parse data file
ParameterTree parameterSet;
@@ -146,15 +148,16 @@ int main (int argc, char *argv[]) try
const bool startFromFile = parameterSet.get<bool>("startFromFile");
const std::string resultPath = parameterSet.get("resultPath", "");
- // ///////////////////////////////////////
- // Create the grid
- // ///////////////////////////////////////
+ /////////////////////////////////////////////////////////////
+ // CREATE THE GRID
+ /////////////////////////////////////////////////////////////
typedef UGGrid<dim> GridType;
std::shared_ptr<GridType> grid;
FieldVector<double,dim> lower(0), upper(1);
+
if (parameterSet.get<bool>("structuredGrid")) {
lower = parameterSet.get<FieldVector<double,dim> >("lower");
@@ -206,17 +209,49 @@ int main (int argc, char *argv[]) try
typedef GridType::LeafGridView GridView;
GridView gridView = grid->leafGridView();
- // FE basis spanning the FE space that we are working in
- typedef Dune::Functions::LagrangeBasis<GridView, order> FEBasis;
- FEBasis feBasis(gridView);
+ /////////////////////////////////////////////////////////////
+ // DATA TYPES
+ /////////////////////////////////////////////////////////////
+
+ using namespace Dune::Indices;
- // dune-fufem-style FE basis for the transition from dune-fufem to dune-functions
- typedef DuneFunctionsBasis<FEBasis> FufemFEBasis;
- FufemFEBasis fufemFEBasis(feBasis);
+ typedef std::vector<RealTuple<double,dim> > DisplacementVector;
+ typedef std::vector<Rotation<double,dim>> RotationVector;
+ const int dimRotation = Rotation<double,dim>::embeddedDim;
+ typedef TupleVector<DisplacementVector, RotationVector> SolutionType;
- // /////////////////////////////////////////
- // Read Dirichlet values
- // /////////////////////////////////////////
+ /////////////////////////////////////////////////////////////
+ // FUNCTION SPACE
+ /////////////////////////////////////////////////////////////
+
+ using namespace Dune::Functions::BasisFactory;
+ auto compositeBasis = makeBasis(
+ gridView,
+ composite(
+ power<dim>(
+ lagrange<displacementOrder>()
+ ),
+ power<dimRotation>(
+ lagrange<rotationOrder>()
+ )
+ ));
+
+ auto deformationPowerBasis = makeBasis(
+ gridView,
+ power<dim>(
+ lagrange<displacementOrder>()
+ ));
+ typedef Dune::Functions::LagrangeBasis<GridView,displacementOrder> DeformationFEBasis;
+ typedef Dune::Functions::LagrangeBasis<GridView,rotationOrder> OrientationFEBasis;
+ DeformationFEBasis deformationFEBasis(gridView);
+ OrientationFEBasis orientationFEBasis(gridView);
+
+ using CompositeBasis = decltype(compositeBasis);
+ using LocalView = typename CompositeBasis::LocalView;
+
+ /////////////////////////////////////////////////////////////
+ // BOUNDARY DATA
+ /////////////////////////////////////////////////////////////
BitSetVector<1> dirichletVertices(gridView.size(dim), false);
BitSetVector<1> neumannVertices(gridView.size(dim), false);
@@ -240,97 +275,64 @@ int main (int argc, char *argv[]) try
auto neumannBoundary = std::make_shared<BoundaryPatch<GridView>>(gridView, neumannVertices);
BoundaryPatch<GridView> surfaceShellBoundary(gridView, surfaceShellVertices);
+ std::cout << "On rank " << mpiHelper.rank() << ": Dirichlet boundary has " << dirichletBoundary.numFaces() << " faces\n";
std::cout << "On rank " << mpiHelper.rank() << ": Neumann boundary has " << neumannBoundary->numFaces() << " faces\n";
std::cout << "On rank " << mpiHelper.rank() << ": Shell boundary has " << surfaceShellBoundary.numFaces() << " faces\n";
-
- BitSetVector<1> dirichletNodes(feBasis.size(), false);
-#if DUNE_VERSION_LT(DUNE_FUFEM, 2, 7)
- using FufemFEBasis = DuneFunctionsBasis<FEBasis>;
- FufemFEBasis fufemFeBasis(feBasis);
- constructBoundaryDofs(dirichletBoundary,fufemFeBasis,dirichletNodes);
-#else
- constructBoundaryDofs(dirichletBoundary,feBasis,dirichletNodes);
-#endif
-
- BitSetVector<1> neumannNodes(feBasis.size(), false);
-#if DUNE_VERSION_LT(DUNE_FUFEM, 2, 7)
- constructBoundaryDofs(*neumannBoundary,fufemFeBasis,neumannNodes);
-#else
- constructBoundaryDofs(*neumannBoundary,feBasis,neumannNodes);
-#endif
-
- BitSetVector<1> surfaceShellNodes(feBasis.size(), false);
-#if DUNE_VERSION_LT(DUNE_FUFEM, 2, 7)
- constructBoundaryDofs(surfaceShellBoundary,fufemFeBasis,surfaceShellNodes);
-#else
- constructBoundaryDofs(surfaceShellBoundary,feBasis,surfaceShellNodes);
-#endif
-
- BitSetVector<blocksize> dirichletDofs(feBasis.size(), false);
- for (size_t i=0; i<feBasis.size(); i++)
- for (int j=0; j<blocksize; j++) {
- if (dirichletNodes[i][0])
- dirichletDofs[i][j] = true;
+ BitSetVector<1> dirichletNodes(compositeBasis.size({0}),false);
+ BitSetVector<1> surfaceShellNodes(compositeBasis.size({1}),false);
+ constructBoundaryDofs(dirichletBoundary,deformationFEBasis,dirichletNodes);
+ constructBoundaryDofs(surfaceShellBoundary,orientationFEBasis,surfaceShellNodes);
+
+ //Create BitVector matching the tangential space
+ const int dimRotationTangent = Rotation<double,dim>::TangentVector::dimension;
+ typedef MultiTypeBlockVector<std::vector<FieldVector<double,dim> >, std::vector<FieldVector<double,dimRotationTangent>> > VectorForBit;
+ typedef Solvers::DefaultBitVector_t<VectorForBit> BitVector;
+
+ BitVector dirichletDofs;
+ dirichletDofs[_0].resize(compositeBasis.size({0}));
+ dirichletDofs[_1].resize(compositeBasis.size({1}));
+ for (size_t i = 0; i < compositeBasis.size({0}); i++){
+ for (int j = 0; j < dim; j++){
+ dirichletDofs[_0][i][j] = dirichletNodes[i][0];
}
- for (size_t i=0; i<feBasis.size(); i++)
- for (int j=3; j<blocksize; j++) {
- if (not surfaceShellNodes[i][0])
- dirichletDofs[i][j] = true;
+ }
+ for (size_t i = 0; i < compositeBasis.size({1}); i++) {
+ for (int j = 0; j < dimRotationTangent; j++){
+ dirichletDofs[_1][i][j] = (not surfaceShellNodes[i][0] or dirichletNodes[i][0]);
}
+ }
- // //////////////////////////
- // Initial iterate
- // //////////////////////////
+ /////////////////////////////////////////////////////////////
+ // INITIAL DATA
+ /////////////////////////////////////////////////////////////
- SolutionTypeCosserat x(feBasis.size());
+ SolutionType x;
+ x[_0].resize(compositeBasis.size({0}));
+ x[_1].resize(compositeBasis.size({1}));
- BlockVector<FieldVector<double,3> > v;
-
//Initial deformation of the underlying substrate
+ BlockVector<FieldVector<double,dim> > displacement(compositeBasis.size({0}));
lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
PythonFunction<FieldVector<double,dim>, FieldVector<double,dim> > pythonInitialDeformation(Python::evaluate(lambda));
- ::Functions::interpolate(fufemFEBasis, v, pythonInitialDeformation);
-
- //Copy over the initial deformation
- for (size_t i=0; i<x.size(); i++) {
- x[i].r = v[i];
- }
-
- ////////////////////////////////////////////////////////
- // Main homotopy loop
- ////////////////////////////////////////////////////////
-
- using namespace Dune::Functions::BasisFactory;
- auto powerBasis = makeBasis(
- gridView,
- power<dim>(
- lagrange<order>(),
- blockedInterleaved()
- ));
+ Dune::Functions::interpolate(deformationPowerBasis, displacement, pythonInitialDeformation);
- BlockVector<FieldVector<double,dim>> identity;
- Dune::Functions::interpolate(powerBasis, identity, [](FieldVector<double,dim> x){ return x; });
+ BlockVector<FieldVector<double,dim> > identity(compositeBasis.size({0}));
+ Dune::Functions::interpolate(deformationPowerBasis, identity, [](FieldVector<double,dim> x){ return x; });
- BlockVector<FieldVector<double,dim> > displacement(v.size());
- for (int i = 0; i < v.size(); i++) {
- displacement[i] = v[i];
+ for (int i = 0; i < displacement.size(); i++) {
+ x[_0][i] = displacement[i]; //Copy over the initial deformation to the deformation part of x
+ displacement[i] -= identity[i]; //Subtract identity to get the initial displacement as a function
}
- displacement -= identity;
-
- auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement);
- auto localDisplacementFunction = localFunction(displacementFunction);
-
- FieldVector<double,dim> neumannValues {0,0,0};
-
- if (parameterSet.hasKey("neumannValues"))
- neumannValues = parameterSet.get<FieldVector<double,dim> >("neumannValues");
- std::cout << "Neumann values: " << neumannValues << std::endl;
-
+ auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(deformationPowerBasis, displacement);
// We need to subsample, because VTK cannot natively display real second-order functions
- SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(order-1));
- vtkWriter.addVertexData(localDisplacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
+ SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(displacementOrder-1));
+ vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
vtkWriter.write(resultPath + "finite-strain_homotopy_" + parameterSet.get<std::string>("energy") + "_0");
+
+ /////////////////////////////////////////////////////////////
+ // INITIAL SURFACE SHELL DATA
+ /////////////////////////////////////////////////////////////
typedef MultiLinearGeometry<double, dim-1, dim> ML;
std::unordered_map<GridType::GlobalIdSet::IdType, ML> geometriesOnShellBoundary;
@@ -339,6 +341,7 @@ int main (int argc, char *argv[]) try
// Read in the grid deformation
if (startFromFile) {
+ // Create a basis of order 1 in order to deform the geometries on the surface shell boundary
const std::string pathToGridDeformationFile = parameterSet.get("pathToGridDeformationFile", "");
// for this, we create a basis of order 1 in order to deform the geometries on the surface shell boundary
auto feBasisOrder1 = makeBasis(
@@ -427,15 +430,30 @@ int main (int argc, char *argv[]) try
}
}
+ /////////////////////////////////////////////////////////////
+ // MAIN HOMOTOPY LOOP
+ /////////////////////////////////////////////////////////////
+ FieldVector<double,dim> neumannValues {0,0,0};
+ if (parameterSet.hasKey("neumannValues"))
+ neumannValues = parameterSet.get<FieldVector<double,dim> >("neumannValues");
+ std::cout << "Neumann values: " << neumannValues << std::endl;
+
+ // Vertex Normals for the 3D-Part
+ std::vector<UnitVector<double,dim> > vertexNormals(gridView.size(dim));
+ Dune::FieldVector<double,dim> vertexNormalRaw = {0,0,1};
+ for (int i = 0; i< vertexNormals.size(); i++) {
+ UnitVector vertexNormal(vertexNormalRaw);
+ vertexNormals[i] = vertexNormal;
+ }
+ //Function for the Cosserat material parameters
const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
Python::Reference surfaceShellClass = Python::import(materialParameters.get<std::string>("surfaceShellParameters"));
Python::Callable surfaceShellCallable = surfaceShellClass.get("SurfaceShellParameters");
-
Python::Reference pythonObject = surfaceShellCallable();
PythonFunction<Dune::FieldVector<double, dim>, double> fThickness(pythonObject.get("thickness"));
PythonFunction<Dune::FieldVector<double, dim>, Dune::FieldVector<double, 2>> fLame(pythonObject.get("lame"));
- for (int i=0; i<numHomotopySteps; i++)
+ for (int i = 0; i < numHomotopySteps; i++)
{
double homotopyParameter = (i+1)*(1.0/numHomotopySteps);
@@ -444,17 +462,13 @@ int main (int argc, char *argv[]) try
// ////////////////////////////////////////////////////////////
-#if DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 8)
- std::shared_ptr<NeumannFunction> neumannFunction;
- neumannFunction = std::make_shared<NeumannFunction>(neumannValues, homotopyParameter);
-#else
// A constant vector-valued function, for simple Neumann boundary values
std::shared_ptr<std::function<Dune::FieldVector<double,dim>(Dune::FieldVector<double,dim>)>> neumannFunctionPtr;
neumannFunctionPtr = std::make_shared<std::function<Dune::FieldVector<double,dim>(Dune::FieldVector<double,dim>)>>([&](FieldVector<double,dim> ) {
return neumannValues * (-homotopyParameter);
});
-#endif
+ const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
if (mpiHelper.rank() == 0)
{
std::cout << "Homotopy step: " << i << ", parameter: " << homotopyParameter << std::endl;
@@ -462,9 +476,12 @@ int main (int argc, char *argv[]) try
materialParameters.report();
}
- // Assembler using ADOL-C
std::cout << "Selected energy is: " << parameterSet.get<std::string>("energy") << std::endl;
-#if DUNE_VERSION_GTE(DUNE_ELASTICITY, 2,8)
+
+ // /////////////////////////////////////////////////
+ // Create the energy functional
+ // /////////////////////////////////////////////////
+
std::shared_ptr<Elasticity::LocalDensity<dim,ValueType>> elasticDensity;
if (parameterSet.get<std::string>("energy") == "stvenantkirchhoff")
elasticDensity = std::make_shared<Elasticity::StVenantKirchhoffDensity<dim,ValueType>>(materialParameters);
@@ -477,123 +494,96 @@ int main (int argc, char *argv[]) try
if (parameterSet.get<std::string>("energy") == "mooneyrivlin")
elasticDensity = std::make_shared<Elasticity::MooneyRivlinDensity<dim,ValueType>>(materialParameters);
- if(!elasticDensity)
- DUNE_THROW(Exception, "Error: Selected energy not available!");
-
- auto elasticEnergy = std::make_shared<LocalIntegralEnergy<GridView, FEBasis::LocalView::Tree::FiniteElement, ValueType>>(elasticDensity);
- auto neumannEnergy = std::make_shared<NeumannEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- ValueType> >(neumannBoundary,neumannFunctionPtr);
- auto elasticAndNeumann = std::make_shared<Dune::SumEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- ValueType>>(elasticEnergy, neumannEnergy);
-#else
-#if DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 7)
- std::shared_ptr<LocalFEStiffness<GridView,
-#else
- std::shared_ptr<Elasticity::LocalEnergy<GridView,
-#endif
- FEBasis::LocalView::Tree::FiniteElement,
- std::vector<Dune::FieldVector<ValueType, dim>> > > elasticEnergy;
-
- if (parameterSet.get<std::string>("energy") == "stvenantkirchhoff")
- elasticEnergy = std::make_shared<StVenantKirchhoffEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- ValueType> >(materialParameters);
-#if !DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 7)
- if (parameterSet.get<std::string>("energy") == "mooneyrivlin")
- elasticEnergy = std::make_shared<MooneyRivlinEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- ValueType> >(materialParameters);
-#endif
-
- if (parameterSet.get<std::string>("energy") == "neohooke")
- elasticEnergy = std::make_shared<NeoHookeEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- ValueType> >(materialParameters);
-
- if (parameterSet.get<std::string>("energy") == "hencky")
- elasticEnergy = std::make_shared<HenckyEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- ValueType> >(materialParameters);
-
- if (parameterSet.get<std::string>("energy") == "exphencky")
- elasticEnergy = std::make_shared<ExpHenckyEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- ValueType> >(materialParameters);
-
- if(!elasticEnergy)
+ if(!elasticDensity)
DUNE_THROW(Exception, "Error: Selected energy not available!");
- auto neumannEnergy = std::make_shared<NeumannEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- ValueType> >(neumannBoundary.get(),neumannFunction.get());
-
- auto elasticAndNeumann = std::make_shared<SumEnergy<GridView,
- FEBasis::LocalView::Tree::FiniteElement,
- ValueType>>(elasticEnergy, neumannEnergy);
-#endif
-
- using LocalEnergyBase = GFE::LocalEnergy<FEBasis,RigidBodyMotion<adouble, dim> >;
-
- std::shared_ptr<LocalEnergyBase> surfaceCosseratEnergy;
- std::vector<UnitVector<double,3> > vertexNormals(gridView.size(3));
- Dune::FieldVector<double,3> vertexNormalRaw = {0,0,1};
- for (int i = 0; i< vertexNormals.size(); i++) {
- UnitVector vertexNormal(vertexNormalRaw);
- vertexNormals[i] = vertexNormal;
- }
-
- surfaceCosseratEnergy = std::make_shared<SurfaceCosseratEnergy<FEBasis,RigidBodyMotion<adouble, dim>, adouble, adouble>>(materialParameters, std::move(vertexNormals), &surfaceShellBoundary, std::move(geometriesOnShellBoundary), fThickness, fLame);
-
- std::shared_ptr<LocalEnergyBase> totalEnergy;
- totalEnergy = std::make_shared<GFE::SumCosseratEnergy<FEBasis,RigidBodyMotion<adouble, dim>, adouble>> (elasticAndNeumann, surfaceCosseratEnergy);
-
-
- LocalGeodesicFEADOLCStiffness<FEBasis,
- TargetSpace> localGFEADOLCStiffness(totalEnergy.get());
+ auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,dim>>>(elasticDensity);
+ auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,dim>>>(neumannBoundary,*neumannFunctionPtr);
+ auto surfaceCosseratEnergy = std::make_shared<GFE::SurfaceCosseratEnergy<CompositeBasis, RealTuple<ValueType,dim>, Rotation<ValueType,dim> >>(materialParameters, std::move(vertexNormals), &surfaceShellBoundary, std::move(geometriesOnShellBoundary), fThickness, fLame);
- GeodesicFEAssembler<FEBasis,TargetSpace> assembler(gridView, &localGFEADOLCStiffness);
+ GFE::SumEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,targetDim>> sumEnergy;
+ sumEnergy.addLocalEnergy(neumannEnergy);
+ sumEnergy.addLocalEnergy(elasticEnergy);
+ sumEnergy.addLocalEnergy(surfaceCosseratEnergy);
+ MixedLocalGFEADOLCStiffness<CompositeBasis,
+ RealTuple<double,dim>,
+ Rotation<double,dim> > localGFEADOLCStiffness(&sumEnergy);
+ MixedGFEAssembler<CompositeBasis,
+ RealTuple<double,dim>,
+ Rotation<double,dim> > mixedAssembler(compositeBasis, &localGFEADOLCStiffness);
////////////////////////////////////////////////////////
// Set Dirichlet values
////////////////////////////////////////////////////////
- Python::Reference dirichletValuesClass = Python::import(parameterSet.get<std::string>("dirichletValues"));
+ BitSetVector<RealTuple<double,dim>::TangentVector::dimension> deformationDirichletDofs(dirichletDofs[_0].size(), false);
+ for(int i = 0; i < dirichletDofs[_0].size(); i++)
+ for (int j = 0; j < RealTuple<double,dim>::TangentVector::dimension; j++)
+ deformationDirichletDofs[i][j] = dirichletDofs[_0][i][j];
+ BitSetVector<Rotation<double,dim>::TangentVector::dimension> orientationDirichletDofs(dirichletDofs[_1].size(), false);
+ for(int i = 0; i < dirichletDofs[_1].size(); i++)
+ for (int j = 0; j < Rotation<double,dim>::TangentVector::dimension; j++)
+ orientationDirichletDofs[i][j] = dirichletDofs[_1][i][j];
+ Python::Reference dirichletValuesClass = Python::import(parameterSet.get<std::string>("dirichletValues"));
Python::Callable C = dirichletValuesClass.get("DirichletValues");
// Call a constructor.
Python::Reference dirichletValuesPythonObject = C(homotopyParameter);
// Extract object member functions as Dune functions
- PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > deformationDirichletValues(dirichletValuesPythonObject.get("deformation"));
- PythonFunction<FieldVector<double,dim>, FieldMatrix<double,3,3> > orientationDirichletValues(dirichletValuesPythonObject.get("orientation"));
-
- std::vector<FieldVector<double,3> > ddV;
- ::Functions::interpolate(fufemFEBasis, ddV, deformationDirichletValues, dirichletDofs);
-
- std::vector<FieldMatrix<double,3,3> > dOV;
- ::Functions::interpolate(fufemFEBasis, dOV, orientationDirichletValues, dirichletDofs);
-
- for (size_t j=0; j<x.size(); j++)
- if (dirichletNodes[j][0])
- {
- x[j].r = ddV[j];
- x[j].q.set(dOV[j]);
+ PythonFunction<FieldVector<double,dim>, FieldVector<double,targetDim> > deformationDirichletValues(dirichletValuesPythonObject.get("deformation"));
+ PythonFunction<FieldVector<double,dim>, FieldMatrix<double,targetDim,targetDim> > rotationalDirichletValues(dirichletValuesPythonObject.get("orientation"));
+
+ BlockVector<FieldVector<double,targetDim> > ddV;
+ Dune::Functions::interpolate(deformationPowerBasis, ddV, deformationDirichletValues, deformationDirichletDofs);
+
+ BlockVector<FieldMatrix<double,targetDim,targetDim> > dOV;
+ Dune::Functions::interpolate(orientationFEBasis, dOV, rotationalDirichletValues, orientationDirichletDofs);
+
+ for (int i = 0; i < compositeBasis.size({0}); i++)
+ if (dirichletDofs[_0][i][0])
+ x[_0][i] = ddV[i];
+ for (int i = 0; i < compositeBasis.size({1}); i++)
+ if (dirichletDofs[_1][i][0])
+ x[_1][i].set(dOV[i]);
+
+#if !MIXED_SPACE
+ //The MixedRiemannianTrustRegionSolver can treat the Displacement and Orientation Space as separate ones
+ //The RiemannianTrustRegionSolver can only treat the Displacement and Rotation together in a RigidBodyMotion
+ //Therefore, x and the dirichletDofs are converted to a RigidBodyMotion structure, as well as the Hessian and Gradient that are returned by the assembler
+ typedef RigidBodyMotion<double, dim> RBM;
+ std::vector<RBM> xRBM(compositeBasis.size({0}));
+ BitSetVector<RBM::TangentVector::dimension> dirichletDofsRBM(compositeBasis.size({0}), false);
+ for (int i = 0; i < compositeBasis.size({0}); i++) {
+ for (int j = 0; j < dim; j ++) { // Displacement part
+ xRBM[i].r[j] = x[_0][i][j];
+ dirichletDofsRBM[i][j] = dirichletDofs[_0][i][j];
}
- // /////////////////////////////////////////////////
- // Create the solver and solve
- // /////////////////////////////////////////////////
+ xRBM[i].q = x[_1][i]; // Rotation part
+ for (int j = dim; j < RBM::TangentVector::dimension; j ++)
+ dirichletDofsRBM[i][j] = dirichletDofs[_1][i][j-dim];
+ }
+ typedef Dune::GFE::GeodesicFEAssemblerWrapper<CompositeBasis, DeformationFEBasis, RBM, RealTuple<double, dim>, Rotation<double,dim>> GFEAssemblerWrapper;
+ GFEAssemblerWrapper assembler(&mixedAssembler, deformationFEBasis);
+#endif
+
+ ///////////////////////////////////////////////////
+ // Create the chosen solver and solve!
+ ///////////////////////////////////////////////////
if (parameterSet.get<std::string>("solvertype") == "multigrid") {
- RiemannianTrustRegionSolver<FEBasis,TargetSpace> solver;
+#if MIXED_SPACE
+ MixedRiemannianTrustRegionSolver<GridType, CompositeBasis, DeformationFEBasis, RealTuple<double,dim>, OrientationFEBasis, Rotation<double,dim>> solver;
solver.setup(*grid,
- &assembler,
+ &mixedAssembler,
+ deformationFEBasis,
+ orientationFEBasis,
x,
- dirichletDofs,
+ deformationDirichletDofs,
+ orientationDirichletDofs,
tolerance,
maxSolverSteps,
initialTrustRegionRadius,
@@ -607,48 +597,75 @@ int main (int argc, char *argv[]) try
solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
solver.setInitialIterate(x);
solver.solve();
-
x = solver.getSol();
+#else
+ RiemannianTrustRegionSolver<DeformationFEBasis, RBM, GFEAssemblerWrapper> solver;
+ solver.setup(*grid,
+ &assembler,
+ xRBM,
+ dirichletDofsRBM,
+ tolerance,
+ maxSolverSteps,
+ initialTrustRegionRadius,
+ multigridIterations,
+ mgTolerance,
+ mu, nu1, nu2,
+ baseIterations,
+ baseTolerance,
+ instrumented);
+
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+ solver.setInitialIterate(xRBM);
+ solver.solve();
+ xRBM = solver.getSol();
+ for (int i = 0; i < xRBM.size(); i++) {
+ x[_0][i] = xRBM[i].r;
+ x[_1][i] = xRBM[i].q;
+ }
+#endif
} else { //parameterSet.get<std::string>("solvertype") == "cholmod"
- RiemannianProximalNewtonSolver<FEBasis,TargetSpace> solver;
+
+#if MIXED_SPACE
+ DUNE_THROW(Exception, "Error: There is no MixedRiemannianProximalNewtonSolver!");
+#else
+ RiemannianProximalNewtonSolver<DeformationFEBasis, RBM, GFEAssemblerWrapper> solver;
solver.setup(*grid,
&assembler,
- x,
- dirichletDofs,
+ xRBM,
+ dirichletDofsRBM,
tolerance,
maxSolverSteps,
initialRegularization,
instrumented);
- solver.setInitialIterate(x);
+ solver.setInitialIterate(xRBM);
solver.solve();
-
- x = solver.getSol();
+ xRBM = solver.getSol();
+ for (int i = 0; i < xRBM.size(); i++) {
+ x[_0][i] = xRBM[i].r;
+ x[_1][i] = xRBM[i].q;
+ }
+#endif
}
- // /////////////////////////////////////////////////////
- // Solve!
- // /////////////////////////////////////////////////////
std::cout << "Overall calculation took " << overallTimer.elapsed() << " sec." << std::endl;
/////////////////////////////////
// Output result
/////////////////////////////////
- for (size_t i=0; i<x.size(); i++)
- v[i] = x[i].r;
// Compute the displacement
- BlockVector<FieldVector<double,dim> > displacement(v.size());
- for (int i = 0; i < v.size(); i++) {
- displacement[i] = v[i];
+ BlockVector<FieldVector<double,dim> > displacement(compositeBasis.size({0}));
+ for (int i = 0; i < compositeBasis.size({0}); i++) {
+ for (int j = 0; j < dim; j++) {
+ displacement[i][j] = x[_0][i][j];
+ }
+ displacement[i] -= identity[i];
}
- displacement -= identity;
-
- auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement);
- auto localDisplacementFunction = localFunction(displacementFunction);
+ auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(deformationPowerBasis, displacement);
// We need to subsample, because VTK cannot natively display real second-order functions
- SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(order-1));
- vtkWriter.addVertexData(localDisplacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
+ SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(displacementOrder-1));
+ vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
vtkWriter.write(resultPath + "finite-strain_homotopy_" + parameterSet.get<std::string>("energy") + "_" + std::to_string(neumannValues[0]) + "_" + std::to_string(i+1));
}
} catch (Exception& e) {