one-body-sample.cc

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

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

#include <dune/common/exceptions.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/fufem/functionspacebases/p1nodalbasis.hh>
#include <dune/fufem/assemblers/operatorassembler.hh>
#include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>

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

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

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

#include <exception>

int const dim = 2;

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
    double const E = 1e8;
    double const nu = 0.3;
    int const refinements = 5;
    size_t const solver_maxIterations = 1000;
    double const solver_tolerance = 1e-5;

    typedef Dune::YaspGrid<dim> GridType;

    Dune::FieldVector<double, dim> const end_points(
        1); // nth dimension (zero-indexed) goes from 0 to end_points[n]
    Dune::FieldVector<int, dim> const elements(
        2); // number of elements in each direction
    Dune::FieldVector<bool, dim> const periodic(false);

    GridType grid(end_points, elements, periodic, 0);

    grid.globalRefine(refinements);

    typedef P1NodalBasis<GridType::LeafGridView, double> P1Basis;
    P1Basis const p1Basis(grid.leafView());

    OperatorAssembler<P1Basis, P1Basis> const globalAssembler(p1Basis, p1Basis);

    StVenantKirchhoffAssembler<GridType, P1Basis::LocalFiniteElement,
                               P1Basis::LocalFiniteElement> const
    localStiffness(E, nu);

    OperatorType stiffnessMatrix;
    globalAssembler.assemble(localStiffness, stiffnessMatrix);

    VectorType f(grid.size(grid.maxLevel(), dim));
    // Fill right-hand side with semi-random numbers
    for (size_t i = 0; i < f.size(); ++i)
      for (size_t j = 0; j < dim; ++j)
        f[i][j] = i + j;

    VectorType u(grid.size(grid.maxLevel(), dim));
    // Fill initial guess with other semi-random numbers
    for (size_t i = 0; i < f.size(); ++i)
      for (size_t j = 0; j < dim; ++j)
        f[i][j] = i * j;

    // TODO: Why does blockGSStep even provide a default constructor?
    BlockGSStep<OperatorType, VectorType> blockGSStep(stiffnessMatrix, u, f);
    EnergyNorm<OperatorType, VectorType> energyNorm(blockGSStep);

    Dune::BitSetVector<VectorType::block_type::dimension> ignoreNodes(
        grid.size(grid.maxLevel(), dim), false);

    { // Play around with the boundary
      typedef GridType::LeafGridView GridView;
      GridView leafView = grid.leafView();
      typedef GridView::Codim<dim>::Iterator VertexLeafIterator;

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

      size_t bounding_nodes = 0;
      for (VertexLeafIterator it = leafView.begin<dim>();
           it != leafView.end<dim>(); ++it) {
        Dune::GeometryType gt = it->type();
        assert(it->geometry().corners() == 1);
        SmallVector coordinates = it->geometry().corner(0);
        bool bounding(false);
        for (int i = 0; i < dim; ++i) {
          if (coordinates[i] == end_points[i] || coordinates[i] == 0) {
            bounding = true;
            break;
          }
        }
        if (bounding) {
          ++bounding_nodes;
          size_t id = myVertexMapper.map(*it);
          std::cout << "Ignoring id #" << id << std::endl;
          ignoreNodes[id] = true;
        }
      }
      std::cout << "Number of bounding nodes: " << bounding_nodes << std::endl;
    }

    blockGSStep.ignoreNodes_ = &ignoreNodes;

    LoopSolver<VectorType> solver(&blockGSStep, solver_maxIterations,
                                  solver_tolerance, &energyNorm, Solver::FULL);

    solver.solve();
    return 0;
  }
  catch (Dune::Exception& e) {
    Dune::derr << "Dune reported error: " << e << std::endl;
  }
  catch (std::exception& e) {
    std::cout << "Standard exception: " << e.what() << std::endl;
  }
}