bulk-surface.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 <vector>
#include <cmath>
#include <numbers>
#include <thread>
 
#include <dune/common/exceptions.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/common/rangeutilities.hh>
#include <dune/common/stringutility.hh>
 
#include <dune/geometry/type.hh>
 
#include <dune/grid/io/file/gmshreader.hh>
#include <dune/grid/uggrid.hh>
#include <dune/grid/common/gridfactory.hh>
 
#include <dune/istl/bvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/solvers.hh>
 
#include <dune/functions/backends/istlvectorbackend.hh>
#include <dune/functions/common/subdomain.hh>
#include <dune/functions/functionspacebases/compositebasis.hh>
#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/restrictedbasis.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/functions/gridfunctions/gridviewfunction.hh>
 
#include <dune/vtk/vtkwriter.hh>
#include <dune/vtk/datacollectors/discontinuouslagrangedatacollector.hh>
 
#include <dune/assembler/backends/istlvectorbackend.hh>
#include <dune/assembler/backends/istlmatrixbackend.hh>
#include <dune/assembler/defaultglobalassembler.hh>
#include <dune/assembler/parallel/gridcoloring.hh>
#include <dune/assembler/parallel/coloredrangeexecutor.hh>
 
#include <dune/fufem/grid/elementdata.hh>
#include <dune/fufem/logger.hh>
#include <dune/fufem/forms/forms.hh>
#include <dune/fufem/constraints/affineconstraints.hh>
#include <dune/fufem/constraints/boundaryconstraints.hh>
 
#include "utilities.hh"
 
 
 
template<class Constraints, class Basis, class IntersectionSet>
void computeTraceConstraints(Constraints& constraints, const Basis& basis, const IntersectionSet& intersectionSet)
{
  using Dune::Functions::istlVectorBackend;
  auto isIntersectionDOF = constraints.isConstrained();
  auto isIntersectionDOFBackend = istlVectorBackend(isIntersectionDOF);
  Dune::Fufem::markIntersectionDofs(intersectionSet, basis, isIntersectionDOFBackend);
 
  auto&& isConstrainedBackend = istlVectorBackend(constraints.isConstrained());
  auto&& constraintValuesBackend = istlVectorBackend(constraints.constantTerm());
 
  auto localView = basis.localView();
  const auto& gridView = basis.gridView();
  for(auto&& element : elements(gridView))
  {
    localView.bind(element);
    for(auto i : Dune::range(localView.tree().size()))
    {
      auto i_local = localView.tree().localIndex(i);
      auto i_global = localView.index(i_local);
      if (not isIntersectionDOFBackend[i_global])
      {
        isConstrainedBackend[i_global] = true;
        constraintValuesBackend[i_global] = 0;
      }
    }
  }
  constraints.resolveDependencies();
}
int main (int argc, char *argv[]) try
{
 
  auto log = Dune::Fufem::makeLogger(std::cout, "[% 8.3f] [% 8.3f] ");
 
  Dune::MPIHelper::instance(argc, argv);
 
  log("MPI initialized");
  auto parameters = Dune::ParameterTree();
  Dune::ParameterTreeParser::readOptions(argc, argv, parameters);
 
  auto meshFile     = parameters.get<std::string>("meshFile", "bulk-surface.msh");
  auto refinements  = parameters.get<std::size_t>("refinements", 0);
  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);
  auto alpha        = parameters.get<double>("alpha", 10.0);
  auto delta        = parameters.get<double>("delta", 2.0);
 
  log("Command line processed");
 
  auto rhs_A = [] (const auto& x) {
    if (x[0]>1)
      return 10;
    else
      return 0;
    };
  auto rhs_B = [] (const auto& x) { return 0; };
  auto rhs_I = [] (const auto& x) { return 0; };
 
  auto dirichletValues = [] (const auto& x) { return 0; };
  using Grid = Dune::UGGrid<2>;
 
  auto factory = Dune::GridFactory<Grid>();
  auto reader = Dune::GmshReader<Grid>(meshFile, factory);
  auto gridPtr = factory.createGrid();
  auto& grid = *gridPtr;
  auto elementData = Dune::Fufem::ElementData(reader, factory, grid);
 
  log("Grid created");
 
  grid.globalRefine(refinements);
  auto gridView = grid.leafGridView();
  elementData.update(gridView);
 
  log("Grid refined");
  using Dune::Functions::Experimental::SubDomain;
  using Dune::Functions::Experimental::SubDomainInterface;
 
  auto subDomain_A = SubDomain(gridView);
  for(auto&& element : elements(gridView))
    if (elementData(element) == 1)
      subDomain_A.insertElement(element);
 
  log(Dune::formatString("SubDomain A contains %d elements", subDomain_A.indexSet().size(0)));
 
  auto subDomain_B = SubDomain(gridView);
  for(auto&& element : elements(gridView))
    if (elementData(element) == 2)
      subDomain_B.insertElement(element);
 
  log(Dune::formatString("SubDomain B contains %d elements", subDomain_B.indexSet().size(0)));
 
  auto interface = SubDomainInterface(subDomain_A, subDomain_B);
  using namespace Dune::Functions::BasisFactory;
  using namespace Dune::Functions::Experimental::BasisFactory;
 
  auto basis = makeBasis(gridView,
    composite(
      restrict(lagrange<1>(), subDomain_A),
      restrict(lagrange<1>(), subDomain_B),
      restrict(lagrange<1>(), subDomain_A),
      flatLexicographic()));
 
  using namespace Dune::Indices;
  auto subDomainBasis_A = Dune::Functions::subspaceBasis(basis, _0);
  auto subDomainBasis_B = Dune::Functions::subspaceBasis(basis, _1);
  auto interfaceBasis = Dune::Functions::subspaceBasis(basis, _2);
 
  log(Dune::formatString("Basis created with dimension %d", basis.dimension()));
  using Vector = Dune::BlockVector<double>;
  using Matrix = Dune::BCRSMatrix<double>;
 
  auto sol = Vector();
  auto rhs = Vector();
  auto stiffnessMatrix = Matrix();
  auto constraints = Dune::Fufem::makeAffineConstraints(basis);
 
  computeBoundaryConstraints(constraints, subDomainBasis_B, dirichletValues);
 
  log("Boundary constraints computed");
 
  computeTraceConstraints(constraints, interfaceBasis, interface);
 
  log("Trace constraints computed");
  using namespace Dune::Fufem::Forms;
 
  auto normal = faceNormal(gridView);
  auto P = compose([&](auto n) {
    constexpr auto dim = decltype(n)::size();
    Dune::FieldMatrix<double, dim, dim> T;
    for(auto i : Dune::range(dim))
      for(auto j : Dune::range(dim))
        T[i][j] = (i==j) - n[i]*n[j];
    return T;
  }, normal);
 
  auto u_A = trialFunction(subDomainBasis_A);
  auto v_A = testFunction(subDomainBasis_A);
 
  auto u_B = trialFunction(subDomainBasis_B);
  auto v_B = testFunction(subDomainBasis_B);
 
  auto u_I = trialFunction(interfaceBasis);
  auto v_I = testFunction(interfaceBasis);
 
  auto f_A = Coefficient(Dune::Functions::makeGridViewFunction(rhs_A, gridView));
  auto f_B = Coefficient(Dune::Functions::makeGridViewFunction(rhs_B, gridView));
  auto f_I = Coefficient(Dune::Functions::makeGridViewFunction(rhs_I, gridView));
 
  auto a_A = integrate(dot(grad(u_A), grad(v_A)), subDomain_A);
  auto a_B = integrate(dot(grad(u_B), grad(v_B)), subDomain_B);
  auto a_I = integrate(delta*dot(P * grad(u_I), grad(v_I)), interface);
 
  auto a_coupling = integrate(alpha*(u_A - u_I)*(v_A - v_I) + alpha*(outside(u_B) - u_I)*(outside(v_B) - v_I), interface);
 
  auto a = a_A + a_B + a_I + a_coupling;
 
  auto l_A = integrate(f_A*v_A, subDomain_A);
  auto l_B = integrate(f_B*v_B, subDomain_B);
  auto l_I = integrate(f_I*v_I, interface);
 
  auto l = l_A + l_B + l_I;
  auto matrixBackend = Dune::Assembler::ISTLMatrixBackend(stiffnessMatrix);
  auto rhsBackend = Dune::Assembler::ISTLVectorBackend(rhs);
  auto assembler = Dune::Assembler::Assembler(basis, basis);
 
  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();
 
  auto matrixOperator = Dune::MatrixAdapter<Matrix,Vector,Vector>(stiffnessMatrix);
 
  auto preconditioner = makeAMGPreconditioner<Vector>(stiffnessMatrix);
 
  log("Preconditioner created");
 
  auto solver = Dune::CGSolver<Vector>(
                  matrixOperator,
                  preconditioner,
                  tol,
                  maxIter,
                  verbosity
                );
 
  log("Solver created");
 
  auto solverStatistics = Dune::InverseOperatorResult();
 
  solver.apply(sol, rhs, solverStatistics);
 
  log("Linear system solved");
 
  constraints.interpolate(sol);
  auto vtkWriter = Dune::Vtk::UnstructuredGridWriter(Dune::Vtk::DiscontinuousLagrangeDataCollector(basis.gridView(), 1));
  vtkWriter.addPointData(u_A|sol, Dune::VTK::FieldInfo("sol_A", Dune::VTK::FieldInfo::Type::scalar, 1));
  vtkWriter.addPointData(u_B|sol, Dune::VTK::FieldInfo("sol_B", Dune::VTK::FieldInfo::Type::scalar, 1));
  vtkWriter.addPointData((u_A + u_B)|sol, Dune::VTK::FieldInfo("sol_bulk", Dune::VTK::FieldInfo::Type::scalar, 1));
  vtkWriter.addPointData(u_I|sol, Dune::VTK::FieldInfo("sol_I", Dune::VTK::FieldInfo::Type::scalar, 1));
  vtkWriter.write("bulk-surface");
 
  log("Solution written to vtk file");
}
catch (Dune::Exception& e)
{
  std::cout << e.what() << std::endl;
  return 1;
}