Add a GeodesicFEAssemblerWrapper - it wraps a MixedGFEAssembler so it can be...

Add a GeodesicFEAssemblerWrapper - it wraps a MixedGFEAssembler so it can be used like a GeodesicFEAssembler

The GeodesicFEAssemblerWrapper assembles the Gradient and the Hessian using the MixedGFEAssembler and then
resturctures them so they can be used with the normal RiemannianTRSolver.
This only works, if the FE spaces have the same order.

dune/gfe/geodesicfeassemblerwrapper.hh

+#ifndef GLOBAL_GEODESIC_FE_ASSEMBLERWRAPPER_HH
+#define GLOBAL_GEODESIC_FE_ASSEMBLERWRAPPER_HH
+
+#include <dune/gfe/mixedgfeassembler.hh>
+#include <dune/common/tuplevector.hh>
+
+namespace Dune::GFE {
+
+/** \brief A wrapper that wraps a MixedGFEAssembler into an assembler that does not distinguish between the two finite element spaces
+
+    It reimplements the same methods as these two and transfers the structure of the gradient and the Hessian
+ */
+
+template <class Basis, class ScalarBasis, class TargetSpace, class MixedSpace0, class MixedSpace1>
+class 
+GeodesicFEAssemblerWrapper {
+
+    typedef typename Basis::GridView GridView;
+
+    //! Dimension of the grid.
+    enum { gridDim = GridView::dimension };
+
+    //! Dimension of the tangent space
+    enum { blocksize = TargetSpace::TangentVector::dimension };
+    enum { blocksize0 = MixedSpace0::TangentVector::dimension };
+    enum { blocksize1 = MixedSpace1::TangentVector::dimension };
+
+    //!
+    typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
+    typedef typename MixedGFEAssembler<Basis, MixedSpace0, MixedSpace1>::MatrixType MatrixType;
+
+protected:
+    MixedGFEAssembler<Basis, MixedSpace0, MixedSpace1>* mixedAssembler_;
+
+public:
+    const ScalarBasis& basis_;
+
+    /** \brief Constructor for a given grid */
+    GeodesicFEAssemblerWrapper(MixedGFEAssembler<Basis, MixedSpace0, MixedSpace1>* mixedAssembler, ScalarBasis& basis)
+        : mixedAssembler_(mixedAssembler),
+          basis_(basis)
+    {
+        hessianMixed_ = std::make_unique<MatrixType>();
+        // The two spaces from the mixed version need to have the same embeddedDim as the Target Space
+        assert(MixedSpace0::embeddedDim + MixedSpace1::embeddedDim == TargetSpace::embeddedDim);
+        assert(blocksize0 + blocksize1 == blocksize);
+    }
+
+    /** \brief Assemble the tangent stiffness matrix and the functional gradient together
+     *
+     * This is more efficient than computing them separately, because you need the gradient
+     * anyway to compute the Riemannian Hessian.
+     */
+    virtual void assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
+                                            Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient,
+                                            Dune::BCRSMatrix<MatrixBlock>& hessian,
+                                            bool computeOccupationPattern=true) const;
+
+    /** \brief Compute the energy of a deformation state */
+    virtual double computeEnergy(const std::vector<TargetSpace>& sol) const;
+
+    /** \brief Get the occupation structure of the Hessian */
+    virtual void getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const;
+
+private:
+    Dune::TupleVector<std::vector<MixedSpace0>,std::vector<MixedSpace1>> splitVector(const std::vector<TargetSpace>& sol) const;
+    std::unique_ptr<MatrixType> hessianMixed_;
+}; // end class
+} //end namespace
+
+
+template <class Basis, class ScalarBasis, class TargetSpace, class MixedSpace0, class MixedSpace1>
+Dune::TupleVector<std::vector<MixedSpace0>,std::vector<MixedSpace1>> Dune::GFE::GeodesicFEAssemblerWrapper<Basis, ScalarBasis, TargetSpace,MixedSpace0,MixedSpace1>::
+splitVector(const std::vector<TargetSpace>& sol) const {
+    using namespace Indices;
+    // Split the solution into the Deformation and the Rotational part
+    auto n = basis_.size();
+    assert (sol.size() == n);
+    Dune::TupleVector<std::vector<MixedSpace0>,std::vector<MixedSpace1>> solutionSplit;
+    solutionSplit[_0].resize(n);
+    solutionSplit[_1].resize(n);
+    for (int i = 0; i < n; i++) {
+        solutionSplit[_0][i] = sol[i].r; // Deformation part
+        solutionSplit[_1][i] = sol[i].q; // Rotational part
+    }
+    return solutionSplit;
+}
+
+template <class Basis, class ScalarBasis, class TargetSpace, class MixedSpace0, class MixedSpace1>
+void Dune::GFE::GeodesicFEAssemblerWrapper<Basis, ScalarBasis, TargetSpace,MixedSpace0,MixedSpace1>::
+getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
+{
+    auto n = basis_.size();
+    nb.resize(n, n);
+    //Retrieve occupation structure for the mixed space and convert it to the non-mixed space
+    //In the mixed space, each index set is for one part of the composite basis
+    //So: nb00 is for the displacement part, nb11 is for the rotation part and nb01 and nb10 (where nb01^T = nb10) is for the coupling
+    Dune::MatrixIndexSet nb00, nb01, nb10, nb11;
+    mixedAssembler_->getMatrixPattern(nb00, nb01, nb10, nb11);
+    auto size0 = nb00.rows();
+    auto size1 = nb11.rows();
+    assert(size0 == size1);
+    assert(size0 == n);
+    //After checking if the sizes match, we can just copy over the occupation pattern
+    //as all occupation patterns work on the same basis combination, so they are all equal
+    nb = nb00;
+}
+
+template <class Basis, class ScalarBasis, class TargetSpace, class MixedSpace0, class MixedSpace1>
+void Dune::GFE::GeodesicFEAssemblerWrapper<Basis, ScalarBasis, TargetSpace,MixedSpace0,MixedSpace1>::
+assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
+                           Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient,
+                           Dune::BCRSMatrix<MatrixBlock>& hessian,
+                           bool computeOccupationPattern) const
+{
+    using namespace Dune::TypeTree::Indices;
+    auto n = basis_.size();
+
+    // Get a split up version of the input
+    auto solutionSplit = splitVector(sol);
+
+    // Define the Matrix and the Gradient in Block Stucture
+    Dune::BlockVector<Dune::FieldVector<double, blocksize0> > gradient0(n);
+    Dune::BlockVector<Dune::FieldVector<double, blocksize1> > gradient1(n);
+    
+    if (computeOccupationPattern) {
+        Dune::MatrixIndexSet neighborsPerVertex;
+        getNeighborsPerVertex(neighborsPerVertex);
+        neighborsPerVertex.exportIdx((*hessianMixed_)[_0][_0]);
+        neighborsPerVertex.exportIdx((*hessianMixed_)[_0][_1]);
+        neighborsPerVertex.exportIdx((*hessianMixed_)[_1][_0]);
+        neighborsPerVertex.exportIdx((*hessianMixed_)[_1][_1]);
+        neighborsPerVertex.exportIdx(hessian);
+    }
+
+    *hessianMixed_ = 0;
+    mixedAssembler_->assembleGradientAndHessian(solutionSplit[_0], solutionSplit[_1], gradient0, gradient1, *hessianMixed_, false);
+
+    // Transform gradient and hessian back to the non-mixed structure
+    hessian = 0;
+    gradient.resize(n);
+    gradient = 0;
+    for (int i = 0; i < n; i++) {
+        for (int j = 0; j < blocksize0 + blocksize1; j++)
+            gradient[i][j] = j < blocksize0 ? gradient0[i][j] : gradient1[i][j - blocksize0];
+    }
+
+    // All 4 matrices are nxn;
+    // Each FieldMatrix in (*hessianMixed_)[_0][_0] is blocksize0 x blocksize0
+    // Each FieldMatrix in (*hessianMixed_)[_1][_0] is blocksize1 x blocksize0
+    // Each FieldMatrix in (*hessianMixed_)[_0][_1] is blocksize0 x blocksize1
+    // Each FieldMatrix in (*hessianMixed_)[_1][_1] is blocksize1 x blocksize1
+    // The hessian will then be a nxn matrix where each FieldMatrix is (blocksize0+blocksize1)x(blocksize0+blocksize1)
+    
+    for (size_t l = 0; l < blocksize0; l++) {
+        // Extract Upper Left Block of mixed matrix
+        for (auto rowIt = (*hessianMixed_)[_0][_0].begin(); rowIt != (*hessianMixed_)[_0][_0].end(); ++rowIt)
+            for (auto colIt = rowIt->begin(); colIt != rowIt->end(); ++colIt)
+                for(size_t k = 0; k < blocksize0; k++)
+                    hessian[rowIt.index()][colIt.index()][k][l] = (*colIt)[k][l];
+        // Extract Lower Left Block of mixed matrix
+        for (auto rowIt = (*hessianMixed_)[_1][_0].begin(); rowIt != (*hessianMixed_)[_1][_0].end(); ++rowIt)
+            for (auto colIt = rowIt->begin(); colIt != rowIt->end(); ++colIt)
+                for(size_t k = 0; k < blocksize1; k++)
+                    hessian[rowIt.index()][colIt.index()][k + blocksize0][l] = (*colIt)[k][l];
+    }
+    for (size_t l = 0; l < blocksize1; l++) {
+        // Extract Upper Right Block of mixed matrix
+        for (auto rowIt = (*hessianMixed_)[_0][_1].begin(); rowIt != (*hessianMixed_)[_0][_1].end(); ++rowIt)
+            for (auto colIt = rowIt->begin(); colIt != rowIt->end(); ++colIt)
+                for(size_t k = 0; k < blocksize0; k++)
+                        hessian[rowIt.index()][colIt.index()][k][l + blocksize0] = (*colIt)[k][l];
+        // Extract Lower Right Block of mixed matrix
+        for (auto rowIt = (*hessianMixed_)[_1][_1].begin(); rowIt != (*hessianMixed_)[_1][_1].end(); ++rowIt)
+            for (auto colIt = rowIt->begin(); colIt != rowIt->end(); ++colIt)
+                for(size_t k = 0; k < blocksize1; k++)
+                    hessian[rowIt.index()][colIt.index()][k + blocksize0][l + blocksize0] = (*colIt)[k][l];
+    }
+}
+
+template <class Basis, class ScalarBasis, class TargetSpace, class MixedSpace0, class MixedSpace1>
+double Dune::GFE::GeodesicFEAssemblerWrapper<Basis, ScalarBasis, TargetSpace,MixedSpace0,MixedSpace1>::
+computeEnergy(const std::vector<TargetSpace>& sol) const
+{
+    using namespace Dune::TypeTree::Indices;
+    auto solutionSplit = splitVector(sol);
+    return mixedAssembler_->computeEnergy(solutionSplit[_0], solutionSplit[_1]);
+}
+#endif //GLOBAL_GEODESIC_FE_ASSEMBLERWRAPPER_HH
\ No newline at end of file

test/CMakeLists.txt

@@ -26,5 +26,11 @@ dune_add_test(SOURCES harmonicmaptest.cc
               TIMEOUT 600
               CMAKE_GUARD MPI_FOUND)
 
+# Run distributed tests
+dune_add_test(SOURCES geodesicfeassemblerwrappertest.cc
+              MPI_RANKS 1 4
+              TIMEOUT 600
+              CMAKE_GUARD MPI_FOUND)
+
 # Copy the example grid used for testing into the build dir
 file(COPY grids/irregular-square.msh DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/grids)

test/geodesicfeassemblerwrappertest.cc

+#include <config.h>
+
+#include <array>
+
+// Includes for the ADOL-C automatic differentiation library
+// Need to come before (almost) all others.
+#include <adolc/adouble.h>
+#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
+
+#include <dune/common/typetraits.hh>
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/functions/functionspacebases/compositebasis.hh>
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include <dune/grid/utility/structuredgridfactory.hh>
+#include <dune/grid/uggrid.hh>
+
+#include <dune/gfe/cosseratenergystiffness.hh>
+#include <dune/gfe/geodesicfeassembler.hh>
+#include <dune/gfe/harmonicenergy.hh>
+#include <dune/gfe/localgeodesicfeadolcstiffness.hh>
+#include <dune/gfe/localprojectedfefunction.hh>
+#include <dune/gfe/mixedgfeassembler.hh>
+#include <dune/gfe/mixedlocalgfeadolcstiffness.hh>
+#include <dune/gfe/riemanniantrsolver.hh>
+
+#include <dune/gfe/geodesicfeassemblerwrapper.hh>
+
+#include <dune/istl/multitypeblockmatrix.hh>
+#include <dune/istl/multitypeblockvector.hh>
+
+// grid dimension
+const int gridDim = 2;
+
+// target dimension
+const int dim = 3;
+
+//order of the finite element space
+const int displacementOrder = 2;
+const int rotationOrder = 2;
+
+using namespace Dune;
+using namespace Indices;
+
+//differentiation method: ADOL-C
+using ValueType = adouble;
+
+//Types for the mixed space
+using DisplacementVector = std::vector<RealTuple<double,dim>>;
+using RotationVector =  std::vector<Rotation<double,dim>>;
+using Vector = MultiTypeBlockVector<DisplacementVector, RotationVector>;
+const int dimCR = Rotation<double,dim>::TangentVector::dimension; //dimCorrectionRotation = Dimension of the correction for rotations
+using CorrectionType = MultiTypeBlockVector<BlockVector<FieldVector<double,dim> >, BlockVector<FieldVector<double,dimCR>>>;
+
+using MatrixRow0 = MultiTypeBlockVector<BCRSMatrix<FieldMatrix<double,dim,dim>>,  BCRSMatrix<FieldMatrix<double,dim,dimCR>>>;
+using MatrixRow1 = MultiTypeBlockVector<BCRSMatrix<FieldMatrix<double,dimCR,dim>>, BCRSMatrix<FieldMatrix<double,dimCR,dimCR>>>;
+using MatrixType = MultiTypeBlockMatrix<MatrixRow0,MatrixRow1>;
+
+//Types for the Non-mixed space
+using RBM = RigidBodyMotion<double, dim>;
+const static int blocksize = RBM::TangentVector::dimension;
+
+struct NeumannFunction
+    : public Dune::VirtualFunction<FieldVector<double,gridDim>, FieldVector<double,dim> >
+{
+    NeumannFunction(){}
+
+    void evaluate(const FieldVector<double, gridDim>& x, FieldVector<double,dim>& out) const {
+        out = 0;
+        out.axpy(1.0, values_);
+    }
+    FieldVector<double,3> values_ = {3e4,2e4,1e4};
+};  
+using CorrectionTypeWrapped = BlockVector<FieldVector<double, blocksize> >;
+using MatrixTypeWrapped = BCRSMatrix<FieldMatrix<double, blocksize, blocksize> >;
+
+int main (int argc, char *argv[])
+{
+  MPIHelper::instance(argc, argv);
+
+  /////////////////////////////////////////////////////////////////////////
+  //    Create the grid
+  /////////////////////////////////////////////////////////////////////////
+  using GridType = UGGrid<gridDim>;
+  auto grid = StructuredGridFactory<GridType>::createCubeGrid({0,0}, {1,1}, {2,2});
+  grid->globalRefine(2);
+  grid->loadBalance();
+
+  using GridView = GridType::LeafGridView;
+  GridView gridView = grid->leafGridView();
+
+  std::function<bool(FieldVector<double,gridDim>)> isNeumann = [](FieldVector<double,gridDim> coordinate) {
+    return coordinate[0] > 0.99;
+  };
+
+  BitSetVector<1> neumannVertices(gridView.size(gridDim), false);
+  const GridView::IndexSet& indexSet = gridView.indexSet();
+  
+  for (auto&& vertex : vertices(gridView))
+      neumannVertices[indexSet.index(vertex)] = isNeumann(vertex.geometry().corner(0));
+
+  BoundaryPatch<GridView> neumannBoundary(gridView, neumannVertices);
+
+
+  /////////////////////////////////////////////////////////////////////////
+  //  Create a composite basis
+  /////////////////////////////////////////////////////////////////////////
+
+  using namespace Functions::BasisFactory;
+
+  auto compositeBasis = makeBasis(
+    gridView,
+    composite(
+      power<dim>(
+        lagrange<displacementOrder>()
+      ),
+      power<dim>(
+        lagrange<rotationOrder>()
+      )
+  ));
+
+  using CompositeBasis = decltype(compositeBasis);
+  using LocalView = typename CompositeBasis::LocalView;
+
+  /////////////////////////////////////////////////////////////////////////
+  //  Create the energy functions with their parameters
+  /////////////////////////////////////////////////////////////////////////
+
+  //Surface-Cosserat-Energy-Parameters
+  ParameterTree parameters;
+  parameters["thickness"] = "1";
+  parameters["mu"] = "2.7191e+4";
+  parameters["lambda"] = "4.4364e+4";
+  parameters["mu_c"] = "0";
+  parameters["L_c"] = "0.01";
+  parameters["q"] = "2";
+  parameters["kappa"] = "1";
+
+
+  auto neumannFunction = std::make_shared<NeumannFunction>();
+  CosseratEnergyLocalStiffness<decltype(compositeBasis), dim,adouble> cosseratEnergy(parameters,
+                                                                     &neumannBoundary,
+                                                                     neumannFunction,
+                                                                     nullptr);
+  MixedLocalGFEADOLCStiffness<CompositeBasis,
+                              RealTuple<double,dim>,
+                              Rotation<double,dim> > mixedLocalGFEADOLCStiffness(&cosseratEnergy);
+  MixedGFEAssembler<CompositeBasis,
+                    RealTuple<double,dim>,
+                    Rotation<double,dim> > mixedAssembler(compositeBasis, &mixedLocalGFEADOLCStiffness);
+
+  using RBM = RigidBodyMotion<double, dim>;
+  using DeformationFEBasis = Functions::LagrangeBasis<GridView,displacementOrder>;
+  DeformationFEBasis deformationFEBasis(gridView);
+  using GFEAssemblerWrapper = GFE::GeodesicFEAssemblerWrapper<CompositeBasis, DeformationFEBasis, RBM, RealTuple<double, dim>, Rotation<double,dim>>;
+  GFEAssemblerWrapper assembler(&mixedAssembler, deformationFEBasis);
+  
+  /////////////////////////////////////////////////////////////////////////
+  //  Prepare the iterate x where we want to assemble - identity in 2D with z = 0
+  /////////////////////////////////////////////////////////////////////////  
+  auto deformationPowerBasis = makeBasis(
+    gridView,
+    power<gridDim>(
+        lagrange<displacementOrder>()
+  ));
+  BlockVector<FieldVector<double,gridDim> > identity(compositeBasis.size({0}));
+  Functions::interpolate(deformationPowerBasis, identity, [](FieldVector<double,gridDim> x){ return x; });
+  BlockVector<FieldVector<double,dim> > initialDeformation(compositeBasis.size({0}));
+  initialDeformation = 0;
+
+  Vector x;
+  x[_0].resize(compositeBasis.size({0}));
+  x[_1].resize(compositeBasis.size({1}));
+  std::vector<RBM> xRBM(compositeBasis.size({0}));
+  for (int i = 0; i < compositeBasis.size({0}); i++) {
+    for (int j = 0; j < gridDim; j++)
+      initialDeformation[i][j] = identity[i][j];
+    x[_0][i] = initialDeformation[i];
+    for (int j = 0; j < dim; j ++) { // Displacement part
+      xRBM[i].r[j] = x[_0][i][j];
+    }
+    xRBM[i].q = x[_1][i]; // Rotation part
+  }
+
+  //////////////////////////////////////////////////////////////////////////////
+  //  Compute the energy, assemble the Gradient and Hessian using 
+  //  the GeodesicFEAssemblerWrapper and the MixedGFEAssembler and compare!
+  //////////////////////////////////////////////////////////////////////////////
+  CorrectionTypeWrapped gradient;
+  MatrixTypeWrapped hessianMatrix;
+  double energy = assembler.computeEnergy(xRBM);
+  assembler.assembleGradientAndHessian(xRBM, gradient, hessianMatrix, true);
+  double gradientTwoNorm = gradient.two_norm();
+  double gradientInfinityNorm = gradient.infinity_norm();
+  double matrixFrobeniusNorm = hessianMatrix.frobenius_norm();
+
+  CorrectionType gradientMixed;
+  gradientMixed[_0].resize(x[_0].size());
+  gradientMixed[_1].resize(x[_1].size());
+  MatrixType hessianMatrixMixed;
+  double energyMixed = mixedAssembler.computeEnergy(x[_0], x[_1]);
+  mixedAssembler.assembleGradientAndHessian(x[_0], x[_1], gradientMixed[_0], gradientMixed[_1], hessianMatrixMixed, true);
+  double gradientMixedTwoNorm = gradientMixed.two_norm();
+  double gradientMixedInfinityNorm = gradientMixed.infinity_norm();
+  double matrixMixedFrobeniusNorm = hessianMatrixMixed.frobenius_norm();
+
+  if (std::abs(energy - energyMixed)/energyMixed > 1e-8)
+  {
+      std::cerr << std::setprecision(9);
+      std::cerr << "The energy calculated by the GeodesicFEAssemblerWrapper is " << energy << " but "
+                << energyMixed << " (calculated by the MixedGFEAssembler) was expected!" << std::endl;
+      return 1;
+  }
+  if ( std::abs(gradientTwoNorm - gradientMixedTwoNorm)/gradientMixedTwoNorm > 1e-8 ||
+       std::abs(gradientInfinityNorm - gradientMixedInfinityNorm)/gradientMixedInfinityNorm > 1e-8)
+  {
+      std::cerr << std::setprecision(9);
+      std::cerr << "The gradient infinity norm calculated by the GeodesicFEAssemblerWrapper is " << gradientInfinityNorm << " but "
+                << gradientMixedInfinityNorm << " (calculated by the MixedGFEAssembler) was expected!" << std::endl;
+      std::cerr << "The gradient norm calculated by the GeodesicFEAssemblerWrapper is " << gradientTwoNorm << " but "
+                << gradientMixedTwoNorm << " (calculated by the MixedGFEAssembler) was expected!" << std::endl;
+      return 1;
+  }
+
+  if (std::abs(matrixFrobeniusNorm - matrixMixedFrobeniusNorm)/matrixMixedFrobeniusNorm > 1e-8)
+  {
+      std::cerr << std::setprecision(9);
+      std::cerr << "The matrix norm calculated by the GeodesicFEAssemblerWrapper is " << matrixFrobeniusNorm << " but "
+                << matrixMixedFrobeniusNorm << " (calculated by the MixedGFEAssembler) was expected!" << std::endl;
+      return 1;
+  }
+}