Use dune-parmg or use dune-functions-basis

Added the parallel solver form dune-parmg and added the use of dune-functions-basis For backwards compatibility, a flag USE_DUNE_FUNCTIONS_BASES and a deprecated message for the old implementations is placed. TODO for later: Currently, dune_parmg only works with the old scalar basis types.

dune/elasticity/assemblers/feassembler.hh

 #ifndef DUNE_ELASTICITY_ASSEMBLERS_FEASSEMBLER_HH
 #define DUNE_ELASTICITY_ASSEMBLERS_FEASSEMBLER_HH
 
-#include <dune/istl/bcrsmatrix.hh>
 #include <dune/common/fmatrix.hh>
+#include <dune/common/hybridutilities.hh>
+
+#include <dune/istl/bvector.hh>
+#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/matrixindexset.hh>
 #include <dune/istl/matrix.hh>
 
 #include "localfestiffness.hh"
 
+#if USE_DUNE_FUNCTIONS_BASES
+
+namespace Dune
+{
+
+namespace Elasticity
+{
+
+// detector for MultiTypeBlockVector
+template <class ABC> struct is_multi_type_block_vector : std::false_type {};
+template <class ...Args> struct is_multi_type_block_vector<Dune::MultiTypeBlockVector<Args... >> : std::true_type {};
+
+// detector for std::array
+template <class XYZ> struct is_std_array : std::false_type {};
+template <class T, size_t i> struct is_std_array<std::array<T,i>> : std::true_type {};
 
 /** \brief A global FE assembler for variational problems
  */
 template <class Basis, class VectorType>
 class FEAssembler {
 
+    using LocalView = typename Basis::LocalView;
+    using field_type = typename VectorType::value_type::field_type;
+
+    //! Dimension of the grid.
+    enum { gridDim = LocalView::GridView::dimension };
+
+public:
+
+    const Basis basis_;
+
+    //! Partitions to assemble the Gradient and Hessian, is interior here, because when doing parallel calculations, the assembled matrix and gradient are treated as if they do not overlap
+    static constexpr auto paritionsOverlap = Dune::Partitions::interiorBorder;
+
+protected:
+
+    std::shared_ptr<LocalFEStiffness<LocalView,field_type>> localStiffness_;
+
+public:
+
+    /** \brief Constructor for a given grid */
+    FEAssembler(const Basis& basis,
+                std::shared_ptr<LocalFEStiffness<LocalView, field_type>> localStiffness)
+        : basis_(basis),
+          localStiffness_(localStiffness)
+    {}
+
+    /** \brief Assemble the tangent stiffness matrix and the functional gradient together
+     *
+     * This may be more efficient than computing them separately
+     */
+    template<class MatrixType>
+    void assembleGradientAndHessian(const VectorType& sol,
+                                    VectorType& gradient,
+                                    MatrixType& hessian,
+                                    bool computeOccupationPattern=true) const;
+
+    /** \brief Compute the energy of a deformation state */
+    auto computeEnergy(const VectorType& sol) const;
+
+// protected:
+    template <class IndexSet>
+    void getNeighborsPerVertex(IndexSet& nb) const;
+
+
+
+}; // end class
+
+
+template <class Basis, class VectorType>
+template <class IndexSet>
+void FEAssembler<Basis,VectorType>::
+getNeighborsPerVertex(IndexSet& nb) const
+{
+    // do we have an array of index sets?
+    if constexpr ( is_std_array<IndexSet>::value )
+    {
+        const auto k = nb.size();
+        for( size_t i=0; i<k; i++ )
+            for( size_t j=0; j<k; j++ )
+            {
+                // BlockedInterleaved ordering of the powernodes is assumed!
+                nb[i][j].resize( basis_.size({i}), basis_.size({j}) );
+            }
+    }
+    else
+    {
+        const int n = basis_.size();
+        nb.resize(n, n);
+    }
+
+    auto localView = basis_.localView();
+
+    for (const auto& element : elements(basis_.gridView(), paritionsOverlap))
+    {
+        localView.bind(element);
+
+        for (size_t i=0; i<localView.size(); i++)
+        {
+            // multi index!
+            const auto iIdx = localView.index(i);
+
+            for (size_t j=0; j<localView.size(); j++)
+            {
+                // multi index!
+                const auto jIdx = localView.index(j);
+
+                if constexpr ( is_std_array<IndexSet>::value )
+                {
+                    // first entry: matrix block
+                    // second entry: index in the matrix block
+                    nb[iIdx[0]][jIdx[0]].add(iIdx[1], jIdx[1]);
+                }
+                else
+                {
+                    // we add only the first entry -> the second entry refers to the dense blocks
+                    nb.add(iIdx[0], jIdx[0]);
+                }
+            }
+        }
+    }
+}
+
+
+template <class Basis, class VectorType>
+template <class MultiTypeMatrix>
+void FEAssembler<Basis,VectorType>::
+assembleGradientAndHessian(const VectorType& sol,
+                           VectorType& gradient,
+                           MultiTypeMatrix& hessian,
+                           bool computeOccupationPattern) const
+{
+    if constexpr ( is_multi_type_block_vector<VectorType>::value )
+    {
+        constexpr int n = sol.size();
+
+        // currently only for n = 2 :(
+        if ( n != 2 )
+            DUNE_THROW(Exception,"assembleGradientAndHessian is currently only available for one or two block types");
+
+        if (computeOccupationPattern)
+        {
+            std::array<std::array<MatrixIndexSet,n>,n> neighborsPerVertex;
+            getNeighborsPerVertex(neighborsPerVertex);
+
+            // ready for every n > 1
+            Hybrid::forEach(Hybrid::integralRange(Hybrid::size(hessian)), [&](auto&& i) {
+                Hybrid::forEach(Hybrid::integralRange(Hybrid::size(hessian)), [&](auto&& j) {
+                    neighborsPerVertex[i][j].exportIdx(hessian[i][j]);
+                });
+            });
+        }
+
+        // ready for every n > 1
+        Hybrid::forEach(Hybrid::integralRange(Hybrid::size(gradient)), [&](auto&& i) {
+            gradient[i].resize(sol[i].size());
+        });
+
+    }
+    else // is_multi_type_block_vector = false
+    {
+        if (computeOccupationPattern)
+        {
+            Dune::MatrixIndexSet neighborsPerVertex;
+            getNeighborsPerVertex(neighborsPerVertex);
+            neighborsPerVertex.exportIdx(hessian);
+        }
+
+        gradient.resize(sol.size());
+    }
+
+    hessian = 0;
+    gradient = 0;
+
+    auto localView = basis_.localView();
+
+    for (const auto& element : elements(basis_.gridView(), paritionsOverlap))
+    {
+        using namespace Dune::Indices;
+
+        localView.bind(element);
+        const auto size = localView.size();
+
+
+        // Extract local solution
+        std::vector<field_type> localSolution(size);
+        std::vector<field_type> localGradient(size);
+        Matrix<FieldMatrix<field_type,1,1>> localHessian;
+        localHessian.setSize(size, size);
+
+        if constexpr ( is_multi_type_block_vector<VectorType>::value )
+        {
+            for (int i=0; i<size; i++)
+            {
+                auto idx = localView.index(i);
+                idx[0] == 0 ? localSolution[i] = sol[_0][ idx[1] ][ idx[2] ]
+                            : localSolution[i] = sol[_1][ idx[1] ][ idx[2] ];
+            }
+        }
+        else
+        {
+            for (int i=0; i<size; i++)
+            {
+                auto idx = localView.index(i);
+                localSolution[i]   = sol[idx[0]][idx[1]];
+
+
+            }
+        }
+
+        // setup local matrix and gradient
+        localStiffness_->assembleGradientAndHessian(localView, localSolution, localGradient, localHessian);
+
+        if constexpr ( is_multi_type_block_vector<VectorType>::value )
+        {
+            for(int i=0; i<size; i++)
+            {
+                const auto row = localView.index(i);
+
+                row[0] == 0 ? gradient[_0][row[1]][row[2]] += localGradient[i]
+                            : gradient[_1][row[1]][row[2]] += localGradient[i];
+
+                // hessian upper left
+                for (int j=0; j<size; j++ )
+                {
+                    auto col = localView.index(j);
+                    if ( row[0] == 0 and col[0] == 0)
+                        hessian[_0][_0][row[1]][col[1]][row[2]][col[2]] += localHessian[i][j];
+                    else if ( row[0] == 0 and col[0] == 1 )
+                        hessian[_0][_1][row[1]][col[1]][row[2]][col[2]] += localHessian[i][j];
+                    else if ( row[0] == 1 and col[0] == 0 )
+                        hessian[_1][_0][row[1]][col[1]][row[2]][col[2]] += localHessian[i][j];
+                    else
+                        hessian[_1][_1][row[1]][col[1]][row[2]][col[2]] += localHessian[i][j];
+                }
+            }
+        }
+        else
+        {
+            for(int i=0; i<size; i++)
+            {
+                auto row = localView.index(i);
+
+                gradient[row[0]][row[1]] += localGradient[i];
+
+                for (int j=0; j<size; j++ )
+                {
+                    auto col = localView.index(j);
+                    hessian[row[0]][col[0]][row[1]][col[1]] += localHessian[i][j];
+                }
+            }
+        }
+    }
+}
+
+template <class Basis, class VectorType>
+auto FEAssembler<Basis,VectorType>::
+computeEnergy(const VectorType& sol) const
+{
+    using namespace Dune::Indices;
+
+    if constexpr ( is_multi_type_block_vector<VectorType>::value )
+    {
+        constexpr int n = sol.size();
+
+        // currently only for n = 2 :(
+        if ( n != 2 )
+            DUNE_THROW(Exception,"assembleGradientAndHessian is currently only available for one or two block types");
+    }
+
+    double energy = 0;
+
+    if (sol.dim()!=basis_.dimension())
+        DUNE_THROW(Dune::Exception, "Solution vector doesn't match the basis!");
+
+    auto localView = basis_.localView();
+
+    // Loop over all elements
+    for (const auto& element : elements(basis_.gridView(), paritionsOverlap))
+    {
+        localView.bind(element);
+        const auto size = localView.size();
+
+        std::vector<field_type> localSolution(size);
+
+        for (size_t i=0; i<size; i++)
+        {
+            auto idx = localView.index(i);
+
+            if constexpr ( is_multi_type_block_vector<VectorType>::value )
+            {
+                idx[0] == 0 ? localSolution[i] = sol[_0][ idx[1] ][ idx[2] ]
+                            : localSolution[i] = sol[_1][ idx[1] ][ idx[2] ];
+            }
+            else
+            {
+                localSolution[i] = sol[idx[0]][idx[1]];
+            }
+        }
+
+        energy += localStiffness_->energy(localView, localSolution);
+    }
+
+    return energy;
+}
+
+} // namespace Elasticity
+
+} // namespace Dune
+
+
+#else
+
+namespace Dune
+{
+
+/** \brief A global FE assembler for variational problems
+ */
+template <class Basis, class VectorType>
+class DUNE_DEPRECATED_MSG("dune-elasticity with dune-fufem bases is now deprecated. Use dune-functions with LocalView concept! Enable the flag USE_DUNE_FUNCTIONS_BASES")
+FEAssembler
+{
+
     typedef typename Basis::GridView GridView;
 
     //! Dimension of the grid.
@@ -28,6 +348,9 @@ class FEAssembler {
 public:
     const Basis basis_;
 
+    //! Partitions to assemble the Gradient and Hessian, is interiorBorder here, because when doing parallel calculations, the assembled matrix and gradient are treated as if they do not overlap
+    static constexpr auto paritionsOverlap = Dune::Partitions::interiorBorder;
+
 protected:
 
     LocalFEStiffness<GridView,
@@ -70,7 +393,7 @@ getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
 
     nb.resize(n, n);
 
-    for (const auto& element : elements(basis_.getGridView(), Dune::Partitions::interior))
+    for (const auto& element : elements(basis_.getGridView(), paritionsOverlap))
     {
         const typename Basis::LocalFiniteElement& lfe = basis_.getLocalFiniteElement(element);
 
@@ -111,7 +434,7 @@ assembleGradientAndHessian(const VectorType& sol,
     gradient.resize(sol.size());
     gradient = 0;
 
-    for (const auto& element : elements(basis_.getGridView(), Dune::Partitions::interior))
+    for (const auto& element : elements(basis_.getGridView(), paritionsOverlap))
     {
         const int numOfBaseFct = basis_.getLocalFiniteElement(element).localBasis().size();
 
@@ -159,7 +482,7 @@ computeEnergy(const VectorType& sol) const
         DUNE_THROW(Dune::Exception, "Solution vector doesn't match the basis!");
 
     // Loop over all elements
-    for (const auto& element : elements(basis_.getGridView(), Dune::Partitions::interior))
+    for (const auto& element : elements(basis_.getGridView(), paritionsOverlap))
     {
         // Number of degrees of freedom on this element
         size_t nDofs = basis_.getLocalFiniteElement(element).localBasis().size();
@@ -176,4 +499,9 @@ computeEnergy(const VectorType& sol) const
     return energy;
 
 }
+
+} // namespace Dune
+
+#endif // USE_DUNE_FUNCTIONS_BASES
+
 #endif