boundary.cc 3.18 KB
 Praetorius, Simon committed Jan 12, 2019 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 ``````#ifdef HAVE_CONFIG_H #include "config.h" #endif #include #if HAVE_DUNE_SPGRID #include #endif #include #include #include #include #include using namespace AMDiS; using namespace Dune::Indices; // 1 component with polynomial degree 2 using Param = YaspGridBasis; using ElliptProblem = ProblemStat; template void run(SetBoundary setBoundary) { ElliptProblem prob("ellipt"); prob.initialize(INIT_ALL); setBoundary(prob.boundaryManager()); auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0); prob.addMatrixOperator(opL, 0, 0); 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); prob.addVectorOperator(opForce, 0); // set boundary condition auto g = [](auto const& x){ return std::exp(-10.0 * dot(x,x)); }; prob.addDirichletBC(BoundaryType{1}, 0, 0, g); AdaptInfo adaptInfo("adapt"); prob.assemble(adaptInfo); prob.solve(adaptInfo); double errorL2 = integrate(sqr(g - prob.solution(0)), prob.gridView(), 6); msg("error_L2 = {}", errorL2); } void run_periodic() { #if HAVE_DUNE_SPGRID Dune::SPCube cube({0.0,0.0},{1.0,1.0}); Dune::SPDomain domain({cube}, Dune::SPTopology<2>(0b01)); Dune::SPGrid grid(domain, Dune::SPMultiIndex<2>({2,2})); using Grid = Dune::SPGrid; #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; ProblemStat prob("ellipt", grid); prob.initialize(INIT_ALL); prob.boundaryManager().setBoxBoundary({-1,-1,1,1}); auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0); prob.addMatrixOperator(opL, 0, 0); auto opForce = makeOperator(tag::test{}, 1.0, 6); prob.addVectorOperator(opForce, 0); // 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])); }; prob.addDirichletBC(BoundaryType{1}, 0, 0, g); 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); prob.writeFiles(adaptInfo, true); } int main(int argc, char** argv) { AMDiS::init(argc, argv); 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(); AMDiS::finalize(); return 0; }``````