Liebe Gitlab-Nutzerin, lieber Gitlab-Nutzer,
es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Die Konten der externen Nutzer:innen sind ü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. The accounts of external users 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;
}