phasefield2.cc 4.2 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifdef HAVE_CONFIG_H
# include "config.h"
#endif

#if ! HAVE_DUNE_FUNCTIONS
#error "Require dune-functions!"
#endif

#if ! HAVE_DUNE_ALUGRID
#error "Require dune-alugrid!"
#endif

#if ! HAVE_ALBERTA
#error "Require Alberta!"
#endif

#include <iostream>
#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
#include <dune/common/exceptions.hh> // We use exceptions

#include <dune/alugrid/grid.hh>

#include <dune/multimesh/mmhierarchiciterator.hh>
#include <dune/multimesh/mmgridfactory.hh>

#include <dune/geometry/quadraturerules.hh>

#include <dune/grid/albertagrid/albertareader.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh>

#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/mtl/interface/umfpack_solve.hpp>

#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>

#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>

#include "phasefield2.hh"

using namespace Dune;
using namespace Dune::Indices;

struct MTLVector
{
  using value_type = double;

  MTLVector(mtl::dense_vector<double>& vec)
    : vec_(vec) {}

  void resize(std::size_t s) { vec_.change_dim(s); }

  std::size_t size() const { return mtl::size(vec_); }

  decltype(auto) operator[](std::size_t i) { return vec_[i]; }
  decltype(auto) operator[](std::size_t i) const { return vec_[i]; }

  decltype(auto) operator[](ReservedVector<long unsigned int, 1> i) { return vec_[i[0]]; }
  decltype(auto) operator[](ReservedVector<long unsigned int, 1> i) const { return vec_[i[0]]; }

  mtl::dense_vector<double>& vec_;
};

template <class VectorType, class Basis>
void writeFile(VectorType& x, Basis const& basis, std::string filename)
{
  auto gv = basis.gridView();
  auto xWrapper = MTLVector(x);
  auto xFct = Functions::makeDiscreteGlobalBasisFunction<double>(basis, TypeTree::hybridTreePath(), xWrapper);
  VTKWriter<decltype(gv)> vtkWriter(gv);
  vtkWriter.addVertexData(xFct, VTK::FieldInfo("u", VTK::FieldInfo::Type::scalar, 1));
  vtkWriter.write(filename);
}


int main(int argc, char** argv)
{
  MPIHelper::instance(argc, argv);

  assert(argc > 1);
  std::string filename = argv[1];

  using Grid = Dune::ALUGrid<2, 2, Dune::simplex, Dune::conforming>;
  AlbertaReader<Grid> albertaReader;
  GridFactory<Grid> factory;
  albertaReader.readGrid(filename, factory);

  std::unique_ptr<Grid> gridPtr(factory.createGrid());
  auto& grid = *gridPtr;

  double eps = 0.05;
  auto signedDistFct = [](auto const& x) { return x.two_norm() - 1.0; };

  grid.globalRefine(10);
  auto maxLevel = grid.maxLevel();

  // refine grid[1] along phase-field interface
  for (int i = 0; i < 4; ++i) {
    std::cout << "refine " << i << "...\n";
    for (auto const& elem : elements(grid.leafGridView())) {
      auto geo = elem.geometry();
      bool refine = false;
      for (int j = 0; j < geo.corners(); ++j)
        refine = refine || std::abs(signedDistFct(geo.corner(j))) < std::max(1,4-i)*eps;

      if (refine)
        grid.mark(1, elem);
    }

    grid.preAdapt();
    grid.adapt();
    grid.postAdapt();
  }

  using namespace Functions::BasisBuilder;
  auto fine_basis = makeBasis(grid.leafGridView(), lagrange<1>());

  using VectorType = mtl::dense_vector<double>;
  using MatrixType = mtl::compressed2D<double>;

  // interpolate phase-field to fine grid
  VectorType phase(fine_basis.dimension());
  auto phaseWrapper = MTLVector(phase);
  interpolate(fine_basis, phaseWrapper, [eps,d=signedDistFct](auto const& x)
  {
    return 0.5*(1.0 - std::tanh(3.0*d(x)/eps));
  });
  writeFile(phase, fine_basis, "phase");

  // assemble a sequence of systems for finer and finer grids
  for (int l = 2; l < maxLevel; ++l) {
    auto coarse_basis = makeBasis(grid.levelGridView(l), lagrange<1>());

    VectorType rhs(coarse_basis.dimension());
    MatrixType matrix(coarse_basis.dimension(), coarse_basis.dimension());
    assembleSystem(coarse_basis, fine_basis, phase, matrix, rhs, eps);

    VectorType u(coarse_basis.dimension());
    u = 0.0;
    umfpack_solve(matrix, u, rhs);

    writeFile(u, coarse_basis, "u_" + std::to_string(l));
  }
}