one-body-sample.cc

/* -*- mode:c++; mode:semantic -*- */

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <exception>
#include <iostream>

#include <dune/common/exceptions.hh>
#include <dune/common/stdstreams.hh>

#include <dune/common/bitsetvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>

#include <dune/grid/yaspgrid.hh>
#include <dune/grid/common/mcmgmapper.hh>

#include <dune/fufem/assemblers/boundaryfunctionalassembler.hh>
#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
#include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
#include <dune/fufem/assemblers/operatorassembler.hh>
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functions/constantfunction.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>

#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>

#include <dune/solvers/common/numproc.hh> // Solver::FULL
#include <dune/solvers/iterationsteps/blockgsstep.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>

#include <dune/tectonic/globalnonlinearity.hh>
#include <dune/tectonic/myconvexproblem.hh>
#include <dune/tectonic/myblockproblem.hh>

int const dim = 2;

template <class GridView>
void setup_boundary(GridView const &gridView,
                    Dune::FieldVector<double, dim> const &end_points,
                    Dune::BitSetVector<dim> &ignoreNodes,
                    Dune::BitSetVector<1> &neumannNodes,
                    Dune::BitSetVector<1> &frictionalNodes) {
  typedef typename GridView::template Codim<dim>::Iterator VertexLeafIterator;

  typedef Dune::MultipleCodimMultipleGeomTypeMapper<
      GridView, Dune::MCMGVertexLayout> VertexMapper;
  VertexMapper const myVertexMapper(gridView);

  size_t dirichlet_nodes = 0;
  size_t neumann_nodes = 0;
  size_t frictional_nodes = 0;
  for (VertexLeafIterator it = gridView.template begin<dim>();
       it != gridView.template end<dim>(); ++it) {
    Dune::GeometryType const gt = it->type();
    assert(it->geometry().corners() == 1);
    Dune::FieldVector<double, dim> const coordinates = it->geometry().corner(0);
    if (coordinates[0] == end_points[0]) {
      ++dirichlet_nodes;
      size_t const id = myVertexMapper.map(*it);
      ignoreNodes[id] = true;
    } else if (coordinates[1] == 0) {
      ++neumann_nodes;
      size_t const id = myVertexMapper.map(*it);
      neumannNodes[id] = true;
    } else if (coordinates[0] == 0) {
      ++frictional_nodes;
      size_t const id = myVertexMapper.map(*it);
      frictionalNodes[id] = true;
    }
  }
  std::cout << "Number of Neumann nodes: " << neumann_nodes << std::endl;
  std::cout << "Number of Dirichlet nodes: " << dirichlet_nodes << std::endl;
  std::cout << "Number of frictional nodes: " << dirichlet_nodes << std::endl;
}

// Assembles Neumann boundary term in f
template <class GridType, class GridView, class LocalVectorType, class FEBasis>
void assemble_neumann(GridView const &gridView, FEBasis const &feBasis,
                      Dune::BitSetVector<1> const &neumannNodes,
                      Dune::BlockVector<LocalVectorType> &
                          f) { // constant sample function on neumann boundary
  BoundaryPatch<GridView> neumannBoundary(gridView, neumannNodes);
  LocalVectorType SampleVector;
  // FIXME: random values
  SampleVector[0] = 1;
  SampleVector[1] = 2;
  ConstantFunction<LocalVectorType, LocalVectorType> fNeumann(SampleVector);
  NeumannBoundaryAssembler<GridType, LocalVectorType> neumannBoundaryAssembler(
      fNeumann);

  BoundaryFunctionalAssembler<FEBasis> boundaryFunctionalAssembler(
      feBasis, neumannBoundary);
  boundaryFunctionalAssembler.assemble(
      neumannBoundaryAssembler, f,
      true); // resize the output vector and zero all of its entries
}

// Assembles constant 1-function on frictional boundary in nodalIntegrals
template <class GridType, class GridView, class LocalVectorType, class FEBasis>
void assemble_frictional(
    GridView const &gridView, FEBasis const &feBasis,
    Dune::BitSetVector<1> const &frictionalNodes,
    std::vector<Dune::FieldVector<double, 1>> &nodalIntegrals) {
  BoundaryPatch<GridView> frictionalBoundary(gridView, frictionalNodes);
  ConstantFunction<LocalVectorType, Dune::FieldVector<double, 1>>
  constantOneFunction(1);
  NeumannBoundaryAssembler<GridType, Dune::FieldVector<double, 1>>
  frictionalBoundaryAssembler(constantOneFunction);

  BoundaryFunctionalAssembler<FEBasis> boundaryFunctionalAssembler(
      feBasis, frictionalBoundary);
  boundaryFunctionalAssembler.assemble(
      frictionalBoundaryAssembler, nodalIntegrals,
      true); // resize the output vector and zero all of its entries
}

int main() {
  try {
    typedef Dune::FieldVector<double, dim> SmallVector;
    typedef Dune::FieldMatrix<double, dim, dim> SmallMatrix;
    typedef Dune::BCRSMatrix<SmallMatrix> OperatorType;
    typedef Dune::BlockVector<SmallVector> VectorType;

    // FIXME: Random values
    size_t const runs = 3;
    double const E = 1;
    double const nu = 0.3;
    int const refinements = 5;
    size_t const solver_maxIterations = 100;
    double const solver_tolerance = 1e-4;

    // {{{ Set up grid
    typedef Dune::YaspGrid<dim> GridType;
    Dune::FieldVector<double, dim> const end_points(
        1); // nth dimension (zero-indexed) goes from 0 to end_points[n]
    GridType grid(
        end_points,
        Dune::FieldVector<int, dim>(2), // number of elements in each direction
        Dune::FieldVector<bool, dim>(false), // non-periodic in each direction
        0);                                  // zero overlap (whatever that is)
    grid.globalRefine(refinements);

    typedef GridType::LeafGridView GridView;
    GridView const leafView = grid.leafView();
    // }}}

    // Set up nodal basis
    typedef P1NodalBasis<GridType::LeafGridView, double> P1Basis;
    P1Basis const p1Basis(leafView);

    // Assemble elastic force on the body
    StVenantKirchhoffAssembler<GridType, P1Basis::LocalFiniteElement,
                               P1Basis::LocalFiniteElement> const
    localStiffness(E, nu);
    OperatorType stiffnessMatrix;
    OperatorAssembler<P1Basis, P1Basis>(p1Basis, p1Basis)
        .assemble(localStiffness, stiffnessMatrix);
    EnergyNorm<OperatorType, VectorType> energyNorm(stiffnessMatrix);

    Dune::BitSetVector<dim> ignoreNodes(grid.size(grid.maxLevel(), dim), false);
    Dune::BitSetVector<1> neumannNodes(grid.size(grid.maxLevel(), dim), false);
    Dune::BitSetVector<1> frictionalNodes(grid.size(grid.maxLevel(), dim),
                                          false);
    setup_boundary<GridType::LeafGridView>(leafView, end_points, ignoreNodes,
                                           neumannNodes, frictionalNodes);

    typedef MyConvexProblem<OperatorType, VectorType, Dune::LinearFunction>
    MyConvexProblemType;
    typedef MyBlockProblem<MyConvexProblemType> MyBlockProblemType;
    GenericNonlinearGS<MyBlockProblemType> nonlinearGSStep;
    nonlinearGSStep.ignoreNodes_ = &ignoreNodes;

    VectorType u1(grid.size(grid.maxLevel(), dim));
    u1 = 0;
    VectorType u2 = u1;

    for (size_t run = 1; run <= runs; ++run) {
      VectorType f(grid.size(grid.maxLevel(), dim));
      f = 0;
      VectorType neumannTerm(grid.size(grid.maxLevel(), dim));
      assemble_neumann<GridType, GridType::LeafGridView, SmallVector, P1Basis>(
          leafView, p1Basis, neumannNodes, neumannTerm);
      f += neumannTerm;

      // {{{ Assemble terms for the nonlinearity
      std::vector<Dune::FieldVector<double, 1>> nodalIntegrals;
      assemble_frictional<GridType, GridType::LeafGridView, SmallVector,
                          P1Basis>(leafView, p1Basis, frictionalNodes,
                                   nodalIntegrals);

      // TODO: populate on S_F
      std::vector<double> normalStress;
      normalStress.resize(grid.size(grid.maxLevel(), dim));
      std::fill(normalStress.begin(), normalStress.end(), 0.0);

      // TODO: populate on S_F
      std::vector<double> coefficientOfFriction;
      coefficientOfFriction.resize(grid.size(grid.maxLevel(), dim));
      std::fill(coefficientOfFriction.begin(), coefficientOfFriction.end(),
                0.0);

      Dune::GlobalNonlinearity<dim, Dune::LinearFunction> myGlobalNonlinearity(
          coefficientOfFriction, normalStress, nodalIntegrals);
      // }}}

      {
        MyConvexProblemType myConvexProblem(stiffnessMatrix,
                                            myGlobalNonlinearity, f, u1);
        MyBlockProblemType myBlockProblem(myConvexProblem);
        nonlinearGSStep.setProblem(u1, myBlockProblem);

        LoopSolver<VectorType> solver(&nonlinearGSStep, solver_maxIterations,
                                      solver_tolerance, &energyNorm,
                                      Solver::FULL); // Solver::QUIET);
        solver.solve();
      }

      // Use a linear solver for comparison; should return roughly the
      // same results if phi vanishes (e.g. because the normalstress is zero)
      {
        // TODO: Why does blockGSStep even provide a default constructor?
        BlockGSStep<OperatorType, VectorType> blockGSStep(stiffnessMatrix, u2,
                                                          f);
        blockGSStep.ignoreNodes_ = &ignoreNodes;

        LoopSolver<VectorType> solver(&blockGSStep, solver_maxIterations,
                                      solver_tolerance, &energyNorm,
                                      Solver::FULL); // Solver::QUIET);
        solver.solve();
      }
    }

    VectorType diff = u2;
    diff -= u1;
    std::cout << "Infinity norm of the difference of the two solutions: "
              << diff.infinity_norm() << std::endl;
  }
  catch (Dune::Exception &e) {
    Dune::derr << "Dune reported error: " << e << std::endl;
  }
  catch (std::exception &e) {
    std::cout << "Standard exception: " << e.what() << std::endl;
  }
}