Add exmaple with primal plasticity formulation

This shows how point-wise operators can be used to handle
matrix-valued unknowns.

src/CMakeLists.txt

@@ -18,3 +18,6 @@ dune_target_enable_all_packages("benchmark")
 
 add_executable("showcase" showcase.cc)
 dune_target_enable_all_packages("showcase")
+
+add_executable("primal-plasticity" primal-plasticity.cc)
+dune_target_enable_all_packages("primal-plasticity")

src/primal-plasticity.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#include <config.h>
+
+#include <array>
+#include <vector>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/indices.hh>
+
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
+
+#include <dune/istl/matrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/matrixindexset.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/multitypeblockmatrix.hh>
+#include <dune/istl/multitypeblockvector.hh>
+
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/backends/istlvectorbackend.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/functionspacebases/compositebasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/functions/gridfunctions/gridviewfunction.hh>
+
+#include <dune/fufem/parallel/elementcoloring.hh>
+#include <dune/fufem/assemblers/dunefunctionsoperatorassembler.hh>
+#include <dune/fufem/assemblers/dunefunctionsfunctionalassembler.hh>
+#include <dune/fufem/backends/istlmatrixbackend.hh>
+
+#include <dune/fufem-forms/forms.hh>
+
+
+
+int main (int argc, char *argv[]) try
+{
+  std::size_t threadCount = 4;
+  if (argc>1)
+    threadCount = std::stoul(std::string(argv[1]));
+
+  // Set up MPI, if available
+  Dune::MPIHelper::instance(argc, argv);
+
+  ///////////////////////////////////
+  //   Generate the grid
+  ///////////////////////////////////
+
+  const int dim = 2;
+  using GridType = Dune::YaspGrid<dim>;
+  Dune::FieldVector<double,dim> upperRight = {1, 1};
+  std::array<int,dim> elements = {{400, 400}};
+//  std::array<int,dim> elements = {{4, 4}};
+  GridType grid(upperRight,elements);
+
+
+  // Problem parameters
+  double lambda = 10;
+  double mu = 6.5;
+  double k1 = 3;
+  double k2 = 1; // ???
+
+  /////////////////////////////////////////////////////////
+  //   Choose a finite element space
+  /////////////////////////////////////////////////////////
+
+  using namespace Dune::Indices;
+  using namespace Dune::Functions::BasisBuilder;
+
+  // Dimension of symmetric trace-free matrix space
+  const int dp = dim*(dim+1)/2 - 1;
+
+  auto basis = makeBasis(
+    grid.leafGridView(),
+    composite(
+      power<dim>(lagrange<1>()),
+      power<dp>(lagrange<0>()),
+      lagrange<0>()
+    ));
+
+  /////////////////////////////////////////////////////////
+  //   Stiffness matrix and right hand side vector
+  /////////////////////////////////////////////////////////
+
+  using M00 = Dune::BCRSMatrix<Dune::FieldMatrix<double,dim,dim>>;
+  using M01 = Dune::BCRSMatrix<Dune::FieldMatrix<double,dim,dp>>;
+  using M02 = Dune::BCRSMatrix<Dune::Fufem::SingleColumnMatrix<Dune::FieldMatrix<double,dim,1>>>;
+
+  using M10 = Dune::BCRSMatrix<Dune::FieldMatrix<double,dp,dim>>;
+  using M11 = Dune::BCRSMatrix<Dune::FieldMatrix<double,dp,dp>>;
+  using M12 = Dune::BCRSMatrix<Dune::Fufem::SingleColumnMatrix<Dune::FieldMatrix<double,dp,1>>>;
+
+  using M20 = Dune::BCRSMatrix<Dune::Fufem::SingleRowMatrix<Dune::FieldMatrix<double,1,dim>>>;
+  using M21 = Dune::BCRSMatrix<Dune::Fufem::SingleRowMatrix<Dune::FieldMatrix<double,1,dp>>>;
+  using M22 = Dune::BCRSMatrix<double>;
+
+  using Matrix = Dune::MultiTypeBlockMatrix<
+                  Dune::MultiTypeBlockVector<M00, M01, M02>,
+                  Dune::MultiTypeBlockVector<M10, M11, M12>,
+                  Dune::MultiTypeBlockVector<M20, M21, M22>
+                  >;
+
+  using V0 = Dune::BlockVector<Dune::FieldVector<double,dim> >;
+  using V1 = Dune::BlockVector<Dune::FieldVector<double,dp> >;
+  using V2 = Dune::BlockVector<double>;
+
+  using Vector = Dune::MultiTypeBlockVector<V0, V1, V2>;
+
+  /////////////////////////////////////////////////////////
+  //  Assemble the system
+  /////////////////////////////////////////////////////////
+
+  Vector rhs;
+
+  auto rhsBackend = Dune::Functions::istlVectorBackend(rhs);
+
+  rhsBackend.resize(basis);
+  rhs = 0;
+
+  Matrix matrix;
+
+  // Compute colored partition of the grid view for parallel assembly
+  auto gridViewPartition = Dune::Fufem::coloredGridViewPartition(basis.gridView());
+
+  // Setup matrix pattern
+//  {
+//    Dune::Timer timer;
+//    auto matrixBackend = Dune::Fufem::istlMatrixBackend(matrix);
+//    auto operatorAssembler = Dune::Fufem::DuneFunctionsOperatorAssembler{basis, basis};
+//    auto patternBuilder = matrixBackend.patternBuilder();
+
+//    timer.reset();
+//    operatorAssembler.assembleBulkPattern(patternBuilder, gridViewPartition, threadCount);
+//    std::cout << "assembleBulkPattern took " << timer.elapsed() << "s." << std::endl;
+
+//    timer.reset();
+//    patternBuilder.setupMatrix();
+//    matrixBackend.assign(0);
+//    std::cout << "setupMatrix took " << timer.elapsed() << "s." << std::endl;
+//  }
+
+  // Assemble matrix
+  {
+    Dune::Timer timer;
+    using namespace Dune::Indices;
+    using namespace Dune::Fufem::Forms;
+
+    // Some point-wise linear operators used in the problem definition
+
+    // Symmetric part of a matrix
+    auto symmetricPart = RangeOperator([](const auto& M) {
+      return 0.5*(M + M.transposed());
+    });
+
+    // Trace of a matrix
+    auto trace = RangeOperator([&](const auto& M) {
+      double tr = 0;
+      for (auto k : Dune::range(dim))
+        tr += M[k][k];
+      return tr;
+    });
+
+    // Coordinate isometry for trace free symmetric matrices
+    // This can be used to convert coefficient vector into matrices
+    auto B = RangeOperator([](const auto& p) {
+      if constexpr (dim==2)
+        return 1.0/std::sqrt(2.0) * Dune::FieldMatrix<double,2,2>{{p[0], p[1]},{p[1], -p[0]}};
+    });
+
+    // Identity matrix (unfortunately, there's no operator+ for FieldMatrix and ScaledIdentityMatrix)
+    auto Id = Dune::FieldMatrix<double,dim,dim>(0);
+    for (auto k : Dune::range(dim))
+      Id[k][k] = 1;
+
+    // Elasticity tensor (for isotropic material)
+    auto C = RangeOperator([&](const auto& e) {
+      return 2*mu*e + lambda*Id*trace(e);
+    });
+
+    // Symmetric gradient
+    auto E = [&](const auto& v) {
+      return symmetricPart(grad(v));
+    };
+
+    // Ansatz functions
+    auto u = trialFunction(basis, _0);
+    auto p = trialFunction(basis, _1);
+    auto eta = trialFunction(basis, _2);
+    auto P = B(p);
+
+    // Test functions
+    auto v = testFunction(basis, _0);
+    auto q = testFunction(basis, _1);
+    auto nu = testFunction(basis, _2);
+    auto Q = B(q);
+
+    // Bilinear form for primal plasticity formulation with kinematic and isotropic hardening
+    auto A = integrate(dot(C(E(u)-P),E(v)-Q) + k1*dot(P,Q) + k2*dot(eta,nu));
+
+    // Alternatively we could e.g. use:
+//     auto sigma = C(E(u)-P);
+//     auto A = integrate(dot(sigma,E(v)-Q)  + k1*dot(P,Q) + k2*dot(eta,nu));
+
+    auto operatorAssembler = Dune::Fufem::DuneFunctionsOperatorAssembler{basis, basis};
+    auto matrixBackend = Dune::Fufem::istlMatrixBackend(matrix);
+    operatorAssembler.assembleBulk(matrixBackend, A, gridViewPartition, threadCount);
+
+    std::cout << "Assembling the problem took " << timer.elapsed() << "s" << std::endl;
+  }
+
+
+}
+catch (Dune::Exception& e) {
+  std::cout << e.what() << std::endl;
+}