ellipt.cc 3.08 KB
Newer Older
 Praetorius, Simon committed Sep 17, 2018 1 2 3 ``````#ifdef HAVE_CONFIG_H #include "config.h" #endif `````` Praetorius, Simon committed May 04, 2016 4 5 ``````#include `````` Praetorius, Simon committed Jan 28, 2018 6 7 8 9 ``````#include #include #include #include `````` Praetorius, Simon committed Jul 03, 2018 10 ``````#include `````` Praetorius, Simon committed May 04, 2016 11 `````` `````` Praetorius, Simon committed May 04, 2016 12 ``````using namespace AMDiS; `````` Praetorius, Simon committed Jul 03, 2018 13 ``````using namespace Dune::Indices; `````` Praetorius, Simon committed May 04, 2016 14 `````` `````` Praetorius, Simon committed May 04, 2016 15 ``````// 1 component with polynomial degree 1 `````` Praetorius, Simon committed Sep 23, 2018 16 ``````using Param = YaspGridBasis; `````` Praetorius, Simon committed Jul 03, 2018 17 ``````using ElliptProblem = ProblemStat; `````` Praetorius, Simon committed May 04, 2016 18 19 20 21 `````` int main(int argc, char** argv) { AMDiS::init(argc, argv); `````` Praetorius, Simon committed Sep 19, 2016 22 `````` `````` Praetorius, Simon committed Jul 03, 2018 23 24 25 26 27 28 29 30 31 32 33 34 35 `````` int numLevels = AMDIS_DIM == 2 ? 8 : 5; if (argc > 2) numLevels = std::atoi(argv[2]); 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 g = [](auto const& x){ return std::exp(-10.0 * dot(x,x)); }; auto grad_g = [](auto const& x) -> FieldMatrix { return {-20.0 * std::exp(-10.0 * dot(x,x)) * x}; }; `````` Praetorius, Simon committed Jan 22, 2018 36 `````` `````` Praetorius, Simon committed Jul 04, 2018 37 `````` ElliptProblem prob("ellipt"); `````` Praetorius, Simon committed May 04, 2016 38 `````` prob.initialize(INIT_ALL); `````` Praetorius, Simon committed Sep 19, 2016 39 `````` `````` Praetorius, Simon committed Dec 18, 2017 40 `````` auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0); `````` Praetorius, Simon committed Jul 04, 2018 41 `````` prob.addMatrixOperator(opL, _0, _0); `````` Praetorius, Simon committed Sep 19, 2016 42 `````` `````` Praetorius, Simon committed Jul 03, 2018 43 `````` auto opForce = makeOperator(tag::test{}, f, 6); `````` Praetorius, Simon committed Dec 08, 2017 44 `````` prob.addVectorOperator(opForce, _0); `````` Praetorius, Simon committed Sep 19, 2016 45 `````` `````` Praetorius, Simon committed May 04, 2016 46 `````` // set boundary condition `````` Praetorius, Simon committed Jul 03, 2018 47 48 49 50 `````` auto boundary = [](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; }; prob.addDirichletBC(boundary, _0, _0, g); AdaptInfo adaptInfo("adapt"); `````` Praetorius, Simon committed Sep 19, 2016 51 `````` `````` Praetorius, Simon committed Jul 03, 2018 52 53 54 55 `````` std::vector errL2; errL2.reserve(numLevels); std::vector errH1; errH1.reserve(numLevels); std::vector widths; widths.reserve(numLevels); for (int i = 0; i < numLevels; ++i) { `````` Praetorius, Simon committed Jul 04, 2018 56 `````` prob.grid().globalRefine(1); `````` Praetorius, Simon committed Jul 04, 2018 57 `````` auto gridView = prob.gridView(); `````` Praetorius, Simon committed Jul 03, 2018 58 59 60 61 62 63 `````` double h = 0; for (auto const& e : edges(gridView)) h = std::max(h, e.geometry().volume()); widths.push_back(h); `````` Praetorius, Simon committed Oct 02, 2018 64 `````` prob.globalBasis().update(gridView); `````` Praetorius, Simon committed Oct 02, 2018 65 `````` prob.buildAfterAdapt(adaptInfo, Flag(0)); `````` Praetorius, Simon committed Jul 03, 2018 66 67 68 69 70 71 72 73 `````` prob.solve(adaptInfo); double errorL2 = integrate(sqr(g - prob.getSolution(_0)), gridView, 6); errL2.push_back(std::sqrt(errorL2)); double errorH1 = errorL2 + integrate(unary_dot(grad_g - gradientAtQP(prob.getSolution(_0))), gridView, 6); errH1.push_back(std::sqrt(errorH1)); #if WRITE_FILES `````` Praetorius, Simon committed Jul 04, 2018 74 `````` Dune::VTKWriter vtkWriter(gridView); `````` Praetorius, Simon committed Jul 03, 2018 75 76 77 78 `````` vtkWriter.addVertexData(prob.getSolution(_0), Dune::VTK::FieldInfo("u", Dune::VTK::FieldInfo::Type::scalar, 1)); vtkWriter.write("u_" + std::to_string(i)); #endif } `````` Praetorius, Simon committed Sep 19, 2016 79 `````` `````` Praetorius, Simon committed Jul 03, 2018 80 `````` std::cout << "\n"; `````` Praetorius, Simon committed Sep 30, 2018 81 82 83 `````` msg("{:5} | {:12} | {:12} | {:12} | {:12} | {:12}\n", "level", "h", "|u - u*|_L2","|u - u*|_H1","eoc_L2","eoc_H1"); msg("{0:->7}{0:->15}{0:->15}{0:->15}{0:->15}{1:->14}","+","\n"); `````` Praetorius, Simon committed Nov 10, 2017 84 `````` `````` Praetorius, Simon committed Jul 03, 2018 85 86 87 88 `````` std::vector eocL2(numLevels, 0.0), eocH1(numLevels, 0.0); for (int i = 1; i < numLevels; ++i) { eocL2[i] = std::log(errL2[i]/errL2[i-1]) / std::log(widths[i]/widths[i-1]); eocH1[i] = std::log(errH1[i]/errH1[i-1]) / std::log(widths[i]/widths[i-1]); `````` Praetorius, Simon committed May 04, 2016 89 `````` } `````` Praetorius, Simon committed Sep 19, 2016 90 `````` `````` Praetorius, Simon committed Jul 03, 2018 91 `````` for (int i = 0; i < numLevels; ++i) `````` Praetorius, Simon committed Sep 30, 2018 92 93 `````` msg("{:<5} | {:<12} | {:<12} | {:<12} | {:<12} | {:<12}\n", i+1, widths[i], errL2[i], errH1[i], eocL2[i], eocH1[i]); `````` Praetorius, Simon committed Sep 19, 2016 94 `````` `````` Praetorius, Simon committed May 04, 2016 95 96 `````` AMDiS::finalize(); return 0; `````` Praetorius, Simon committed Sep 19, 2016 97 ``}``