Merge branch 'feature/geodesicfeassemblerwrapper' into 'master'

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

See merge request !55

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;
+  }
+}