Merge branch 'feature/move-to-dune-functions-bases' into 'master'

Migrate from dune-fufem to dune-functions bases

See merge request agnumpde/dune-elasticity!18

CHANGELOG.md

@@ -2,7 +2,9 @@
 
 - Introduce class `LocalDensity` for material-specific implementations
 - Introduce class `LocalIntegralEnergy` to work with the densities
+- Local energies and FE assemblers use now `dune-functions` power bases instead of scalar `dune-fufem` bases; the key element is the `LocalView` which contains the information for each element
 
 ## Deprecations and removals
 
 - Redundant implementations of `LocalEnergy` classes are now replaced by combining `LocalDensity` and `LocalIntegralEnergy`
+- Local energies and FE assemblers with `dune-fufem` scalar bases are now deprecated

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/fufem/assemblers/istlbackend.hh>
+
+#include <dune/istl/bvector.hh>
+#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/matrixindexset.hh>
 #include <dune/istl/matrix.hh>
 
+#include <dune/functions/functionspacebases/concepts.hh>
+#include <dune/functions/backends/istlvectorbackend.hh>
+
 #include "localfestiffness.hh"
 
 
-/** \brief A global FE assembler for variational problems
+namespace Dune::Elasticity {
+
+/** \brief A global FE assembler for variational problems for dune-functions global bases
  */
 template <class Basis, class VectorType>
 class FEAssembler {
 
+    using LocalView = typename Basis::LocalView;
+    using field_type = typename VectorType::field_type;
+
+    //! Dimension of the grid.
+    enum { gridDim = LocalView::GridView::dimension };
+
+public:
+
+    const Basis basis_;
+
+    /** \brief Partition type on which to assemble
+     *
+     * We want a fully nonoverlapping partition, and therefore need to skip ghost elements.
+     */
+    static constexpr auto partitionType = Dune::Partitions::interiorBorder;
+
+protected:
+
+    std::shared_ptr<LocalFEStiffness<LocalView,field_type>> localStiffness_;
+
+public:
+
+    /** \brief Constructor for a given basis and local assembler */
+    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& configuration,
+                                    VectorType& gradient,
+                                    MatrixType& hessian,
+                                    bool computeOccupationPattern=true) const;
+
+    /** \brief Compute the energy of a deformation state */
+    auto computeEnergy(const VectorType& configuration) const;
+
+
+}; // end class
+
+
+template <class Basis, class VectorType>
+template <class MatrixType>
+void FEAssembler<Basis,VectorType>::
+assembleGradientAndHessian(const VectorType& configuration,
+                           VectorType& gradient,
+                           MatrixType& hessian,
+                           bool computeOccupationPattern) const
+{
+  // create backends for multi index access
+  auto hessianBackend  = Dune::Fufem::ISTLMatrixBackend(hessian);
+  auto configurationBackend      = Dune::Functions::istlVectorBackend(configuration);
+  auto gradientBackend = Dune::Functions::istlVectorBackend(gradient);
+
+  if (computeOccupationPattern)
+  {
+      auto patternBuilder = hessianBackend.patternBuilder();
+      patternBuilder.resize(basis_,basis_);
+
+      auto localView = basis_.localView();
+
+      for (const auto& element : elements(basis_.gridView(), partitionType))
+      {
+        localView.bind(element);
+        for (size_t i=0; i<localView.size(); i++)
+          for (size_t j=0; j<localView.size(); j++)
+            patternBuilder.insertEntry( localView.index(i), localView.index(j) );
+      }
+      patternBuilder.setupMatrix();
+  }
+
+  gradientBackend.resize(basis_);
+
+  hessian = 0;
+  gradient = 0;
+
+  auto localView = basis_.localView();
+
+  for (const auto& element : elements(basis_.gridView(), partitionType))
+  {
+    localView.bind(element);
+    const auto size = localView.size();
+
+    // Extract local configuration
+    std::vector<field_type> localConfiguration(size);
+    std::vector<field_type> localGradient(size);
+    Matrix<field_type> localHessian;
+    localHessian.setSize(size, size);
+
+    for (int i=0; i<size; i++)
+      localConfiguration[i] = configurationBackend[ localView.index(i) ];
+
+    // setup local matrix and gradient
+    localStiffness_->assembleGradientAndHessian(localView, localConfiguration, localGradient, localHessian);
+
+    for(int i=0; i<size; i++)
+    {
+      auto row = localView.index(i);
+      gradientBackend[row] += localGradient[i];
+
+      for (int j=0; j<size; j++ )
+      {
+        auto col = localView.index(j);
+        // fufem ISTLBackend uses operator() for entry access
+        hessianBackend(row,col) += localHessian[i][j];
+      }
+    }
+  }
+}
+
+template <class Basis, class VectorType>
+auto FEAssembler<Basis,VectorType>::
+computeEnergy(const VectorType& configuration) const
+{
+  double energy = 0;
+
+  if (configuration.dim()!=basis_.dimension())
+      DUNE_THROW(Dune::Exception, "Solution vector doesn't match the basis!");
+
+  auto localView = basis_.localView();
+  // fufem multi index access
+  auto configurationBackend = Functions::istlVectorBackend(configuration);
+
+  // Loop over all elements
+  for (const auto& element : elements(basis_.gridView(), partitionType))
+  {
+      localView.bind(element);
+      const auto size = localView.size();
+
+      std::vector<field_type> localConfiguration(size);
+
+      for (size_t i=0; i<size; i++)
+        localConfiguration[i] = configurationBackend[ localView.index(i) ];
+
+      energy += localStiffness_->energy(localView, localConfiguration);
+  }
+
+  return energy;
+}
+
+} // namespace Dune::Elasticity
+
+
+
+namespace Dune
+{
+
+/** \brief A global FE assembler for variational problems (old fufem bases version)
+ */
+template <class Basis, class VectorType >
+class DUNE_DEPRECATED_MSG("dune-elasticity with dune-fufem bases is now deprecated. Use Elasticity::FEAssembler with dune-functions bases.")
+FEAssembler
+{
+
     typedef typename Basis::GridView GridView;
 
     //! Dimension of the grid.
@@ -182,4 +351,7 @@ computeEnergy(const VectorType& sol) const
     return energy;
 
 }
+
+} // namespace Dune
+
 #endif

dune/elasticity/assemblers/localadolcstiffness.hh

@@ -14,12 +14,130 @@
 
 #include <dune/elasticity/assemblers/localfestiffness.hh>
 
-//#define ADOLC_VECTOR_MODE
+namespace Dune::Elasticity {
 
 /** \brief Assembles energy gradient and Hessian with ADOL-C (automatic differentiation)
  */
-template<class GridView, class LocalFiniteElement, class VectorType>
+template<class LocalView>
 class LocalADOLCStiffness
+    : public LocalFEStiffness<LocalView,double>
+{
+    enum {gridDim=LocalView::GridView::dimension};
+
+public:
+
+    // accepts ADOL_C's "adouble" only
+    LocalADOLCStiffness(std::shared_ptr<Dune::Elasticity::LocalEnergy<LocalView, adouble>> energy)
+    : localEnergy_(energy)
+    {}
+
+    /** \brief Compute the energy at the current configuration */
+    virtual double energy (const LocalView& localView,
+                           const std::vector<double>& localConfiguration) const;
+
+    /** \brief Assemble the local stiffness matrix at the current position
+     *
+     *     This uses the automatic differentiation toolbox ADOL_C.
+     */
+    virtual void assembleGradientAndHessian(const LocalView& localView,
+                                            const std::vector<double>& localConfiguration,
+                                            std::vector<double>& localGradient,
+                                            Dune::Matrix<double>& localHessian);
+
+    std::shared_ptr<Dune::Elasticity::LocalEnergy<LocalView, adouble>> localEnergy_;
+
+};
+
+
+template<class LocalView>
+double LocalADOLCStiffness<LocalView>::
+energy(const LocalView& localView,
+       const std::vector<double>& localConfiguration) const
+{
+    int rank = Dune::MPIHelper::getCollectiveCommunication().rank();
+    double pureEnergy;
+
+    std::vector<adouble> localAConfiguration(localConfiguration.size());
+
+    trace_on(rank);
+
+    adouble energy = 0;
+
+    for (size_t i=0; i<localConfiguration.size(); i++)
+        localAConfiguration[i] <<= localConfiguration[i];
+
+    energy = localEnergy_->energy(localView,localAConfiguration);
+
+    energy >>= pureEnergy;
+
+    trace_off(rank);
+
+    return pureEnergy;
+}
+
+
+template<class LocalView>
+void LocalADOLCStiffness<LocalView>::
+assembleGradientAndHessian(const LocalView& localView,
+                           const std::vector<double>& localConfiguration,
+                           std::vector<double>& localGradient,
+                           Dune::Matrix<double>& localHessian
+                          )
+{
+    int rank = Dune::MPIHelper::getCollectiveCommunication().rank();
+    // Tape energy computation.  We may not have to do this every time, but it's comparatively cheap.
+    energy(localView, localConfiguration);
+
+    /////////////////////////////////////////////////////////////////
+    // Compute the energy gradient
+    /////////////////////////////////////////////////////////////////
+
+    // Compute the actual gradient
+    size_t nDoubles = localConfiguration.size();
+
+    // Compute gradient
+    localGradient.resize(nDoubles);
+    gradient(rank,nDoubles,localConfiguration.data(),localGradient.data());
+
+    /////////////////////////////////////////////////////////////////
+    // Compute Hessian
+    /////////////////////////////////////////////////////////////////
+
+    localHessian.setSize(nDoubles,nDoubles);
+
+    std::vector<double> v(nDoubles);
+    std::vector<double> w(nDoubles);
+
+    std::fill(v.begin(), v.end(), 0.0);
+
+    for (size_t i=0; i<nDoubles; i++)
+    {
+      // Evaluate Hessian multiplied with the i-th canonical basis vector (i.e.: the i-th column of the Hessian)
+      v[i] = 1.;
+
+      int rc= 3;
+      MINDEC(rc, hess_vec(rank, nDoubles, const_cast<double*>(localConfiguration.data()), v.data(), w.data())); // ADOL-C won't accept pointers with const values
+      if (rc < 0)
+        DUNE_THROW(Dune::Exception, "ADOL-C has returned with error code " << rc << "!");
+
+      for (size_t j=0; j<nDoubles; j++)
+        localHessian[i][j] = w[j];
+
+      // Make v the null vector again
+      v[i] = 0.;
+    }
+
+    // TODO: implement ADOLC_VECTOR_MODE variant
+}
+
+} // namespace Dune::Elasticity
+
+
+/** \brief Assembles energy gradient and Hessian with ADOL-C (automatic differentiation)
+ */
+template<class GridView, class LocalFiniteElement, class VectorType>
+class DUNE_DEPRECATED_MSG("dune-elasticity with dune-fufem bases is now deprecated. Use Elasticity::LocalADOLCStiffness with LocalView concept!")
+LocalADOLCStiffness
     : public LocalFEStiffness<GridView,LocalFiniteElement,VectorType>
 {
     // grid types
@@ -38,7 +156,7 @@ public:
     //! Dimension of a tangent space
     enum { blocksize = VectorType::value_type::dimension };
 
-    LocalADOLCStiffness(const Dune::Elasticity::LocalEnergy<GridView, LocalFiniteElement, AVectorType>* energy)
+    LocalADOLCStiffness(const Dune::LocalEnergy<GridView, LocalFiniteElement, AVectorType>* energy)
     : localEnergy_(energy)
     {}
 
@@ -56,7 +174,7 @@ public:
                          const VectorType& localConfiguration,
                          VectorType& localGradient);
 
-    const Dune::Elasticity::LocalEnergy<GridView, LocalFiniteElement, AVectorType>* localEnergy_;
+    const Dune::LocalEnergy<GridView, LocalFiniteElement, AVectorType>* localEnergy_;
 
 };
 

dune/elasticity/assemblers/localenergy.hh

@@ -3,13 +3,33 @@
 
 #include <vector>
 
-namespace Dune {
+namespace Dune::Elasticity {
 
-namespace Elasticity {
+/** \brief Base class for energies defined by integrating over one grid element */
+template<class LocalView, class field_type = double>
+class LocalEnergy
+{
+
+public:
+
+  /** \brief Compute the energy at the current configuration
+   *
+   * \param localView Container with current element and local function basis
+   * \param localConfiguration The coefficients of a FE function on the current element
+   */
+  virtual field_type energy (const LocalView& localView,
+                             const std::vector<field_type>& localConfiguration) const = 0;
+
+};
+
+} // namespace Dune::Elasticity
+
+namespace Dune {
 
 /** \brief Base class for energies defined by integrating over one grid element */
 template<class GridView, class LocalFiniteElement, class VectorType>
-class LocalEnergy
+class DUNE_DEPRECATED_MSG("dune-elasticity with dune-fufem bases is now deprecated. Use Elasticity::LocalEnergy with LocalView concept!")
+LocalEnergy
 {
     typedef typename VectorType::value_type::field_type RT;
     typedef typename GridView::template Codim<0>::Entity Element;
@@ -29,8 +49,6 @@ public:
 
 };
 
-}  // namespace Elasticity
-
 }  // namespace Dune
 
 #endif  // DUNE_ELASTICITY_ASSEMBLERS_LOCALENERGY_HH

dune/elasticity/assemblers/localfestiffness.hh

@@ -6,9 +6,33 @@
 
 #include <dune/elasticity/assemblers/localenergy.hh>
 
-template<class GridView, class LocalFiniteElement, class VectorType>
+namespace Dune::Elasticity {
+
+template<class LocalView, class RT = double>
 class LocalFEStiffness
-: public Dune::Elasticity::LocalEnergy<GridView, LocalFiniteElement, VectorType>
+: public Dune::Elasticity::LocalEnergy<LocalView, RT>
+{
+    enum {gridDim=LocalView::GridView::dimension};
+
+public:
+
+    /** \brief Assemble the local gradient and tangent matrix at the current position
+    */
+    virtual void assembleGradientAndHessian(const LocalView& localView,
+                                            const std::vector<RT>& localConfiguration,
+                                            std::vector<RT>& localGradient,
+                                            Dune::Matrix<RT>& localHessian
+                                           ) = 0;
+
+};
+
+} // namespace Dune::Elasticity
+
+
+template<class GridView, class LocalFiniteElement, class VectorType>
+class DUNE_DEPRECATED_MSG("dune-elasticity with dune-fufem bases is now deprecated. Use Elasticity::LocalFEStiffness with LocalView concept!")
+LocalFEStiffness
+: public Dune::LocalEnergy<GridView, LocalFiniteElement, VectorType>
 {
     // grid types
     typedef typename GridView::Grid::ctype DT;

dune/elasticity/common/trustregionsolver.cc

@@ -8,9 +8,11 @@
 #include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
 #include <dune/fufem/functionspacebases/p1nodalbasis.hh>
 #include <dune/fufem/assemblers/operatorassembler.hh>
+#include <dune/fufem/assemblers/dunefunctionsoperatorassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/massassembler.hh>
 #include <dune/fufem/assemblers/basisinterpolationmatrixassembler.hh>
+#include <dune/fufem/assemblers/istlbackend.hh>
 
 // Using a monotone multigrid as the inner solver
 #include <dune/solvers/iterationsteps/trustregiongsstep.hh>
@@ -27,7 +29,10 @@
 template <class BasisType, class VectorType>
 void TrustRegionSolver<BasisType,VectorType>::
 setup(const typename BasisType::GridView::Grid& grid,
-      const FEAssembler<DuneFunctionsBasis<BasisType>, VectorType>* assembler,
+         std::conditional_t< // do we have a dune-functions basis? -> choose right assembler type
+           Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
+             const Dune::Elasticity::FEAssembler<BasisType,VectorType>*,
+             const Dune::FEAssembler<DuneFunctionsBasis<BasisType>,VectorType>* > assembler,
          const SolutionType& x,
          const Dune::BitSetVector<blocksize>& dirichletNodes,
          double tolerance,
@@ -65,7 +70,11 @@ setup(const typename BasisType::GridView::Grid& grid,
     mgSetup_->ignore(std::const_pointer_cast<Dune::BitSetVector<blocksize> >(ignoreNodes_));
 
     BasisType feBasis(grid.levelGridView(grid.maxLevel()));
-    DuneFunctionsBasis<BasisType> basis(feBasis);
+
+    std::conditional_t< // do we have a dune-functions basis? -> choose right assembler type
+      Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
+        BasisType,
+        DuneFunctionsBasis<BasisType> > basis(feBasis);
 
     innerMultigridStep_ = std::make_unique<Dune::ParMG::Multigrid<VectorType> >();
     innerMultigridStep_->mu_ = mu;
@@ -75,8 +84,11 @@ setup(const typename BasisType::GridView::Grid& grid,
     if (grid.comm().rank()==0)
         std::cout << "Parallel multigrid setup took " << setupTimer.stop() << " seconds." << std::endl;
 #else
-    DuneFunctionsBasis<BasisType> basis(grid.leafGridView());
 
+    std::conditional_t< // do we have a dune-functions basis?
+      Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
+        BasisType,
+        DuneFunctionsBasis<BasisType> > basis(grid.leafGridView());
     // ////////////////////////////////
     //   Create a multigrid solver
     // ////////////////////////////////
@@ -117,15 +129,46 @@ setup(const typename BasisType::GridView::Grid& grid,
     // //////////////////////////////////////////////////////////////////////////////////////
     //   Assemble a Laplace matrix to create a norm that's equivalent to the H1-norm
     // //////////////////////////////////////////////////////////////////////////////////////
-    OperatorAssembler<DuneFunctionsBasis<BasisType>,DuneFunctionsBasis<BasisType>,Dune::Partitions::Interior> operatorAssembler(basis, basis);
 
-    LaplaceAssembler<GridType,
-                     typename DuneFunctionsBasis<BasisType>::LocalFiniteElement,
-                     typename DuneFunctionsBasis<BasisType>::LocalFiniteElement> laplaceStiffness;
-    typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > ScalarMatrixType;
+
+    using ScalarMatrixType = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >;
     ScalarMatrixType localA;
 
-    operatorAssembler.assemble(laplaceStiffness, localA);
+    std::conditional_t< // do we have a dune-functions basis?
+      Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
+        Dune::Fufem::DuneFunctionsOperatorAssembler<BasisType,BasisType>,
+        OperatorAssembler<DuneFunctionsBasis<BasisType>,DuneFunctionsBasis<BasisType>,Dune::Partitions::Interior> > operatorAssembler(basis,basis);
+
+    if constexpr ( Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>() )
+    {
+      using FiniteElement = std::decay_t<decltype(basis.localView().tree().child(0).finiteElement())>;
+      // construct a Fufem Basis Assembler
+      auto laplaceStiffness = LaplaceAssembler<GridType, FiniteElement, FiniteElement>();
+      // transform it to a dune-functions assembler
+      auto localAssembler = [&](const auto& element, auto& localMatrix, auto&& trialLocalView, auto&& ansatzLocalView){
+          laplaceStiffness.assemble(element,
+                                    localMatrix,
+                                    trialLocalView.tree().child(0).finiteElement(),
+                                    ansatzLocalView.tree().child(0).finiteElement());
+      };
+
+
+      auto matrixBackend = Dune::Fufem::istlMatrixBackend(localA);
+      auto patternBuilder = matrixBackend.patternBuilder();
+
+      operatorAssembler.assembleBulkPattern(patternBuilder);
+      patternBuilder.setupMatrix();
+
+      operatorAssembler.assembleBulkEntries(matrixBackend, localAssembler);
+    }
+    else
+    {
+      LaplaceAssembler<GridType,
+                      typename DuneFunctionsBasis<BasisType>::LocalFiniteElement,
+                      typename DuneFunctionsBasis<BasisType>::LocalFiniteElement> laplaceStiffness;
+
+      operatorAssembler.assemble(laplaceStiffness, localA);
+    }
 
     ScalarMatrixType* A = new ScalarMatrixType(localA);
 
@@ -136,12 +179,37 @@ setup(const typename BasisType::GridView::Grid& grid,
     //   This will be used to monitor the gradient
     // //////////////////////////////////////////////////////////////////////////////////////
 
-    MassAssembler<GridType,
-                  typename DuneFunctionsBasis<BasisType>::LocalFiniteElement,
-                  typename DuneFunctionsBasis<BasisType>::LocalFiniteElement> massStiffness;
     ScalarMatrixType localMassMatrix;
 
-    operatorAssembler.assemble(massStiffness, localMassMatrix);
+    if constexpr ( Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>() )
+    {
+      // construct a Fufem Basis Assembler
+      using FiniteElement = std::decay_t<decltype(basis.localView().tree().child(0).finiteElement())>;
+      auto massStiffness = MassAssembler<GridType, FiniteElement, FiniteElement>();
+      // transform it to a dune-functions assembler
+      auto localMassAssembler = [&](const auto& element, auto& localMatrix, auto&& trialLocalView, auto&& ansatzLocalView){
+        massStiffness.assemble(element,
+                               localMatrix,
+                               trialLocalView.tree().child(0).finiteElement(),
+                               ansatzLocalView.tree().child(0).finiteElement());
+      };
+
+      auto massMatrixBackend = Dune::Fufem::istlMatrixBackend(localMassMatrix);
+      auto massPatternBuilder = massMatrixBackend.patternBuilder();
+
+      operatorAssembler.assembleBulkPattern(massPatternBuilder);
+      massPatternBuilder.setupMatrix();
+
+      operatorAssembler.assembleBulkEntries(massMatrixBackend, localMassAssembler);
+    }
+    else
+    {
+      MassAssembler<GridType,
+                    typename DuneFunctionsBasis<BasisType>::LocalFiniteElement,
+                    typename DuneFunctionsBasis<BasisType>::LocalFiniteElement> massStiffness;
+
+      operatorAssembler.assemble(massStiffness, localMassMatrix);
+    }
 
     ScalarMatrixType* massMatrix = new ScalarMatrixType(localMassMatrix);
     l2Norm_ = std::make_shared<H1SemiNorm<CorrectionType> >(*massMatrix);
@@ -151,10 +219,19 @@ setup(const typename BasisType::GridView::Grid& grid,
     // ////////////////////////////////////////////////////////////
 
     hessianMatrix_ = std::make_shared<MatrixType>();
-    Dune::MatrixIndexSet indices(grid_->size(1), grid_->size(1));
-    assembler_->getNeighborsPerVertex(indices);
-    indices.exportIdx(*hessianMatrix_);
-
+    if constexpr ( Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>() )
+    {
+      auto hessianBackend = Dune::Fufem::istlMatrixBackend(*hessianMatrix_);
+      auto hessianPatternBuilder = hessianBackend.patternBuilder();