linear-elasticity.cc

Commented version of this example This example on gitlab

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#include <config.h>
 
#include <array>
#include <cmath>
#include <cstddef>
#include <iostream>
#include <string>
#include <vector>
 
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/common/stringutility.hh>
 
#include <dune/grid/utility/structuredgridfactory.hh>
#include <dune/grid/yaspgrid.hh>
 
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/cholmod.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/scaledidmatrix.hh>
#include <dune/istl/solvers.hh>
 
#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
 
#include <dune/vtk/datacollectors/discontinuouslagrangedatacollector.hh>
#include <dune/vtk/vtkwriter.hh>
 
#include <dune/assembler/backends/istlmatrixbackend.hh>
#include <dune/assembler/backends/istlvectorbackend.hh>
#include <dune/assembler/defaultglobalassembler.hh>
#include <dune/assembler/parallel/coloredrangeexecutor.hh>
#include <dune/assembler/parallel/gridcoloring.hh>
 
#include <dune/fufem/domains/boundary.hh>
#include <dune/fufem/domains/boundarypatch.hh>
#include <dune/fufem/constraints/affineconstraints.hh>
#include <dune/fufem/constraints/boundaryconstraints.hh>
#include <dune/fufem/forms/forms.hh>
#include <dune/fufem/logger.hh>
 
int main (int argc, char *argv[]) try
{
  auto log = Dune::Fufem::makeLogger(std::cout, "[% 8.3f] [% 8.3f] ");
 
  auto parameters = Dune::ParameterTree();
  Dune::ParameterTreeParser::readOptions(argc, argv, parameters);
 
  auto refinements  = parameters.get<std::size_t>("refinements", 1);
  auto threadCount  = parameters.get<std::size_t>("threads", std::thread::hardware_concurrency());
#if not(HAVE_SUITESPARSE_CHOLMOD)
  auto tol          = parameters.get<double>("tol", 1e-4);
  auto maxIter      = parameters.get<std::size_t>("maxIter", 1000);
  auto verbosity    = parameters.get<std::size_t>("verbosity", 2);
#endif
 
  log("Command line processed");
 
  Dune::MPIHelper::instance(argc, argv);
 
  log("MPI initialized");
  const int dim = 3;
 
  using LocalVector = Dune::FieldVector<double,dim>;
  using LocalMatrix = Dune::FieldMatrix<double,dim,dim>;
  double E_Young = 200e9; // in Pa = N/m^2 = 10^9 GPa
  double nu = 0.3;
 
  double lambda = E_Young*nu/((1.-2.*nu) * (1.+nu));
  double mu = E_Young/(2*(1+nu));
 
  double rho = 7.85e3;    // in kg/m^3 = 10^3 g/cm^3
  double g_n = 9.80665;   // in m/s^2 = N/kg
  auto g = LocalVector({0.0, -g_n*rho, 0.0});
  double thickness = .01;
  double length = 1.0;
 
  auto dirichletIndicator = [](auto x) { return x[0] <= 1e-10; };
  auto dirichletValues = [](const auto& x){ return LocalVector(0.0); };
  using Grid = Dune::YaspGrid<dim>;
 
  auto lowerCorner = LocalVector({0, 0, 0});
  auto upperCorner = LocalVector({length, thickness, thickness});
  auto elementCount = std::array{(unsigned int)(length/thickness), 1u, 1u};
 
  auto gridPtr = Dune::StructuredGridFactory<Grid>::createCubeGrid(lowerCorner, upperCorner, elementCount);
  auto& grid = *gridPtr;
 
  log("Grid created");
 
  grid.globalRefine(refinements);
 
  log(Dune::formatString("Grid refined %d times", refinements));
 
  auto gridView = grid.leafGridView();
  auto dirichletPatch = Dune::Fufem::BoundaryPatch(gridView);
  for (const auto& intersection : Dune::Fufem::Boundary(gridView))
    if (dirichletIndicator(intersection.geometry().center()))
      dirichletPatch.insertFace(intersection);
 
  log("Boundary patch set up");
  auto coloredRange = Dune::Assembler::Experimental::coloredElementRange(gridView);
 
  log("Grid coloring computed");
 
  auto executor = Dune::Assembler::Experimental::ColoredRangeExecutor(coloredRange, threadCount);
 
  log(Dune::formatString("Using parallelism with %d threads", threadCount));
  const int order = 2;
 
  using namespace Dune::Functions::BasisFactory;
  auto basis = makeBasis(gridView, power<dim>(lagrange<order>(), blockedInterleaved()));
 
  log(Dune::formatString("Basis created with dimension %d", basis.dimension()));
  using Vector = Dune::BlockVector<LocalVector>;
  using BitVector = std::vector<std::array<char,dim>>;
  using Matrix = Dune::BCRSMatrix<LocalMatrix>;
 
  auto sol = Vector();
  auto rhs = Vector();
  auto stiffnessMatrix = Matrix();
  auto constraints = Dune::Fufem::makeAffineConstraints<BitVector, Vector>(basis);
 
  Dune::Fufem::computeBoundaryConstraints(constraints, basis, dirichletValues, dirichletPatch);
 
  log("Boundary constraints computed");
  using namespace Dune::Fufem::Forms;
 
  auto u = trialFunction(basis);
  auto v = testFunction(basis);
 
  auto Id = Dune::ScaledIdentityMatrix<double,dim>{1.0};
 
  auto E = [&](const auto& v) {
    return symmetrize(grad(v));
  };
 
  auto sigma = 2*mu*E(u) + lambda*Id*trace(E(u));
 
  auto a = integrate(dot(sigma, E(v)));
  auto l = integrate(dot(g,v));
  auto matrixBackend = Dune::Assembler::ISTLMatrixBackend(stiffnessMatrix);
  auto rhsBackend = Dune::Assembler::ISTLVectorBackend(rhs);
  auto assembler = Dune::Assembler::Assembler(basis, basis, executor);
 
  log("Assembler set up");
  auto patternBuilder = matrixBackend.patternBuilder();
 
  patternBuilder.resize(basis, basis);
 
  assembler.assembleMatrixPattern(a, patternBuilder);
 
  log("Matrix pattern assembled");
  patternBuilder.setupMatrix();
 
  log("Matrix set up");
 
  assembler.assembleMatrixEntries(a, matrixBackend);
 
  log("Matrix assembled");
  rhsBackend.resize(basis);
  rhsBackend.setZero();
 
  assembler.assembleVectorEntries(l, rhsBackend);
 
  log("RHS assembled");
  constraints.constrainLinearSystem(stiffnessMatrix, rhs);
 
  log("Boundary condition incorporated into system");
  Dune::Assembler::ISTLVectorBackend(sol).resize(basis);
  Dune::Assembler::ISTLVectorBackend(sol).setZero();
#if HAVE_SUITESPARSE_CHOLMOD
  auto solver = Dune::Cholmod<Vector>();
  solver.setMatrix(stiffnessMatrix);
#else
  auto matrixOperator = Dune::MatrixAdapter<Matrix,Vector,Vector>(stiffnessMatrix);
  auto preconditioner = Dune::SeqILDL<Matrix,Vector,Vector>(stiffnessMatrix, 1.0);
 
  log("Preconditioner created");
 
  auto solver = Dune::CGSolver<Vector>(
                  matrixOperator,
                  preconditioner,
                  tol,
                  maxIter,
                  verbosity
                );
#endif
 
  log("Solver created");
 
  auto solverStatistics = Dune::InverseOperatorResult();
 
  solver.apply(sol, rhs, solverStatistics);
 
  log("Linear system solved");
  auto sigma_vol = 1./3. * trace(sigma);
  auto sigma_dev = sigma - sigma_vol*Id;
  auto sqrt = RangeOperator([](const auto& z) { return std::sqrt(z); });
  auto sigma_von_Mises_sol = sqrt(3./2.*Dune::Fufem::Forms::dot(sigma_dev|sol, sigma_dev|sol));
  auto vtkWriter = Dune::Vtk::UnstructuredGridWriter(Dune::Vtk::DiscontinuousLagrangeDataCollector(gridView, order));
  vtkWriter.addPointData(u|sol, Dune::VTK::FieldInfo("u", Dune::VTK::FieldInfo::Type::vector, dim));
  vtkWriter.addPointData(E(u)|sol, Dune::VTK::FieldInfo("Eu", Dune::VTK::FieldInfo::Type::tensor, dim*dim));
  vtkWriter.addPointData(sigma|sol, Dune::VTK::FieldInfo("sigma", Dune::VTK::FieldInfo::Type::tensor, dim*dim));
  vtkWriter.addPointData(sigma_dev|sol, Dune::VTK::FieldInfo("sigma_dev", Dune::VTK::FieldInfo::Type::tensor, dim*dim));
  vtkWriter.addPointData(sigma_vol|sol, Dune::VTK::FieldInfo("sigma_vol", Dune::VTK::FieldInfo::Type::scalar, 1));
  vtkWriter.addPointData(sigma_von_Mises_sol, Dune::VTK::FieldInfo("sigma_von_Mises", Dune::VTK::FieldInfo::Type::scalar, 1));
  vtkWriter.write("linear-elasticity");
 
  log("Solution written to vtk file");
}
catch (Dune::Exception& e) {
  std::cout << e.what() << std::endl;
  return 1;
}