boundary.cc 3.06 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <iostream>

#if HAVE_DUNE_SPGRID
#include <dune/grid/spgrid.hh>
#endif

#include <amdis/AMDiS.hpp>
#include <amdis/Integrate.hpp>
#include <amdis/ProblemStat.hpp>
13
#include <amdis/LocalOperators.hpp>
14 15 16 17

using namespace AMDiS;

// 1 component with polynomial degree 2
18
using Param   = YaspGridBasis<GRIDDIM, 2>;
19 20 21 22 23 24 25 26 27 28
using ElliptProblem = ProblemStat<Param>;

template <class SetBoundary>
void run(SetBoundary setBoundary)
{
  ElliptProblem prob("ellipt");
  prob.initialize(INIT_ALL);
  setBoundary(prob.boundaryManager());

  auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
29
  prob.addMatrixOperator(opL);
30 31 32 33 34 35 36

  auto f = [](auto const& x) {
    double r2 = dot(x,x);
    double ux = std::exp(-10.0 * r2);
    return -(400.0 * r2 - 20.0 * x.size()) * ux;
  };
  auto opForce = makeOperator(tag::test{}, f, 6);
37
  prob.addVectorOperator(opForce);
38 39 40

  // set boundary condition
  auto g = [](auto const& x){ return std::exp(-10.0 * dot(x,x)); };
41
  prob.addDirichletBC(BoundaryType{1}, g);
42 43 44 45 46 47

  AdaptInfo adaptInfo("adapt");

  prob.assemble(adaptInfo);
  prob.solve(adaptInfo);

48
  double errorL2 = integrate(sqr(g - prob.solution()), prob.gridView(), 6);
49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69
  msg("error_L2 = {}", errorL2);
}

void run_periodic()
{
#if HAVE_DUNE_SPGRID
  Dune::SPCube<double,2> cube({0.0,0.0},{1.0,1.0});
  Dune::SPDomain<double,2> domain({cube}, Dune::SPTopology<2>(0b01));
  Dune::SPGrid<double,2> grid(domain, Dune::SPMultiIndex<2>({2,2}));
  using Grid = Dune::SPGrid<double,2>;
#else
  Dune::YaspGrid<2> grid({1.0,1.0},{2,2},std::bitset<2>("10"),0);
  using Grid = Dune::YaspGrid<2>;
#endif

  using Traits = LagrangeBasis<typename Grid::LeafGridView, 2>;
  ProblemStat<Traits> prob("ellipt", grid);
  prob.initialize(INIT_ALL);
  prob.boundaryManager().setBoxBoundary({-1,-1,1,1});

  auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
70
  prob.addMatrixOperator(opL);
71 72

  auto opForce = makeOperator(tag::test{}, 1.0, 6);
73
  prob.addVectorOperator(opForce);
74 75 76 77 78 79

  // set boundary condition
  auto g = [](auto const& x) {
    return std::sin(0.5*M_PI*x[1])*std::sin(2*M_PI * (0.25 + x[0]))
         + std::cos(0.5*M_PI*x[1])*std::sin(2*M_PI * (-0.25 + x[0]));
  };
80
  prob.addDirichletBC(BoundaryType{1}, g);
81 82 83 84 85 86 87 88 89

  typename ElliptProblem::WorldMatrix A{{1.0,0.0}, {0.0,1.0}};
  typename ElliptProblem::WorldVector b{1.0, 0.0};
  prob.addPeriodicBC(BoundaryType{-1}, A, b);

  AdaptInfo adaptInfo("adapt");

  prob.assemble(adaptInfo);
  prob.solve(adaptInfo);
Praetorius, Simon's avatar
Praetorius, Simon committed
90
  prob.writeFiles(adaptInfo);
91 92 93 94 95
}


int main(int argc, char** argv)
{
96
  Environment env(argc, argv);
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

  auto b = [](auto const& x){ return x[0] < 1.e-8 || x[1] < 1.e-8 || x[0] > 1.0-1.e-8 || x[1] > 1.0-1.e-8; };

  auto setBoundary1 = [](auto& boundaryManager)
  {
    boundaryManager.setBoxBoundary({1,1,1,1});
  };
  run(setBoundary1);

  auto setBoundary2 = [b](auto& boundaryManager)
  {
    boundaryManager.setIndicator([b](auto const& x) { return b(x) ? 1 : 0; });
  };
  run(setBoundary2);

  auto setBoundary3 = [b](auto& boundaryManager)
  {
    boundaryManager.setPredicate(b, 1);
  };
  run(setBoundary3);

  run_periodic();

  return 0;
}