From 89afbd5a1cf80b558f810aeb01d045da937bc163 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Carsten=20Gr=C3=A4ser?= <graeser@mi.fu-berlin.de>
Date: Wed, 24 May 2023 10:58:46 +0200
Subject: [PATCH] Add exmaple with primal plasticity formulation
This shows how point-wise operators can be used to handle
matrix-valued unknowns.
---
src/CMakeLists.txt | 3 +
src/primal-plasticity.cc | 219 +++++++++++++++++++++++++++++++++++++++
2 files changed, 222 insertions(+)
create mode 100644 src/primal-plasticity.cc
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 789fb5d..31dbcd6 100644
--- a/src/CMakeLists.txt
+++ b/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")
diff --git a/src/primal-plasticity.cc b/src/primal-plasticity.cc
new file mode 100644
index 0000000..1bdfae1
--- /dev/null
+++ b/src/primal-plasticity.cc
@@ -0,0 +1,219 @@
+// -*- 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;
+}
--
GitLab