Liebe Gitlab-Nutzer, lieber Gitlab-Nutzer, es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Ein Anmelden über dieses erzeugt ein neues Konto. Das alte Konto ist über den Reiter "Standard" erreichbar. Die Administratoren

Dear Gitlab user, it is now possible to log in to our service using the ZIH login/LDAP. Logging in via this will create a new account. The old account can be accessed via the "Standard" tab. The administrators

boundary.cc 3.02 KB
Newer Older
Praetorius, Simon's avatar
Praetorius, Simon committed
1 2
#include <config.h>

3 4 5 6 7 8 9 10 11
#include <iostream>

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

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

using namespace AMDiS;

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

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

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

  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);
36
  prob.addVectorOperator(opForce);
37 38 39

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

  AdaptInfo adaptInfo("adapt");

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

47
  double errorL2 = integrate(sqr(g - prob.solution()), prob.gridView(), 6);
48 49 50 51 52 53 54 55 56 57 58 59 60 61 62
  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

63
  using Traits = LagrangeBasis<Grid, 2>;
64 65
  ProblemStat<Traits> prob("ellipt", grid);
  prob.initialize(INIT_ALL);
66
  prob.boundaryManager()->setBoxBoundary({-1,-1,1,1});
67 68

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

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

  // 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]));
  };
79
  prob.addDirichletBC(BoundaryType{1}, g);
80 81 82 83 84 85 86 87 88

  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
89
  prob.writeFiles(adaptInfo);
90 91 92 93 94
}


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

  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;
}