From a66171205cd024d03c92c221fb85aec74a162d8b Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Thu, 26 Sep 2019 13:31:29 +0200
Subject: [PATCH] added ellipt demo corrections
---
examples/ellipt.cc | 47 +++++++++++++++++++++++-----------------------
1 file changed, 24 insertions(+), 23 deletions(-)
diff --git a/examples/ellipt.cc b/examples/ellipt.cc
index 411ef195..e032165b 100644
--- a/examples/ellipt.cc
+++ b/examples/ellipt.cc
@@ -1,26 +1,21 @@
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-#include <iostream>
-
#include <amdis/AMDiS.hpp>
#include <amdis/Integrate.hpp>
#include <amdis/LocalOperators.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/common/Literals.hpp>
+#include <dune/grid/yaspgrid.hh>
+using Grid = Dune::YaspGrid<GRIDDIM>;
+
using namespace AMDiS;
using namespace Dune::Indices;
-// 1 component with polynomial degree 1
-using Param = YaspGridBasis<GRIDDIM, 2>;
-using ElliptProblem = ProblemStat<Param>;
-
int main(int argc, char** argv)
{
Environment env(argc, argv);
+ Dune::Timer t;
- int numLevels = GRIDDIM == 2 ? 8 : 5;
+ int numLevels = GRIDDIM == 2 ? 6 : 4;
if (argc > 2)
numLevels = std::atoi(argv[2]);
@@ -34,7 +29,8 @@ int main(int argc, char** argv)
return {-20.0 * std::exp(-10.0 * dot(x,x)) * x};
};
- ElliptProblem prob("ellipt");
+ using Param = LagrangeBasis<Grid, 2>;
+ ProblemStat<Param> prob("ellipt");
prob.initialize(INIT_ALL);
auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
@@ -44,40 +40,44 @@ int main(int argc, char** argv)
prob.addVectorOperator(opForce, 0);
// set boundary condition
- prob.addDirichletBC(BoundaryType{1}, 0, 0, g);
+ prob.addDirichletBC(1, 0, 0, g);
AdaptInfo adaptInfo("adapt");
std::vector<double> errL2; errL2.reserve(numLevels);
std::vector<double> errH1; errH1.reserve(numLevels);
std::vector<double> widths; widths.reserve(numLevels);
- for (int i = 0; i < numLevels; ++i) {
- prob.grid()->globalRefine(1);
- auto gridView = prob.gridView();
+ for (int l = 0; l < numLevels; ++l) {
+ prob.globalRefine(1);
double h = 0;
- for (auto const& e : edges(gridView))
- h = std::max(h, e.geometry().volume());
+ for (auto const& e : elements(prob.gridView(), Dune::Partitions::interior)) {
+ auto refElem = Dune::referenceElement<double,GRIDDIM>(e.type());
+ auto geo = e.geometry();
+ for (int i = 0; i < refElem.size(GRIDDIM-1); ++i) { // edges
+ auto v0 = geo.global(refElem.position(refElem.subEntity(i,GRIDDIM-1,0,GRIDDIM), GRIDDIM));
+ auto v1 = geo.global(refElem.position(refElem.subEntity(i,GRIDDIM-1,1,GRIDDIM), GRIDDIM));
+ h = std::max(h, (v0 - v1).two_norm());
+ }
+ }
widths.push_back(h);
- prob.globalBasis()->update(gridView);
- prob.solutionVector()->compress();
prob.assemble(adaptInfo);
prob.solve(adaptInfo);
- double errorL2 = integrate(sqr(g - prob.solution(0)), gridView, 6);
+ double errorL2 = integrate(sqr(g - prob.solution(0)), prob.gridView(), 6);
errL2.push_back(std::sqrt(errorL2));
- double errorH1 = errorL2 + integrate(unary_dot(grad_g - gradientAtQP(prob.solution(0))), gridView, 6);
+ double errorH1 = errorL2 + integrate(unary_dot(grad_g - gradientAtQP(prob.solution(0))), prob.gridView(), 6);
errH1.push_back(std::sqrt(errorH1));
#if WRITE_FILES
Dune::VTKWriter<typename ElliptProblem::GridView> vtkWriter(gridView);
vtkWriter.addVertexData(prob.solution(0), Dune::VTK::FieldInfo("u", Dune::VTK::FieldInfo::Type::scalar, 1));
- vtkWriter.write("u_" + std::to_string(i));
+ vtkWriter.write("u_" + std::to_string(l));
#endif
}
- std::cout << "\n";
+ msg("");
msg("{:5} | {:12} | {:12} | {:12} | {:12} | {:12}",
"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");
@@ -92,5 +92,6 @@ int main(int argc, char** argv)
msg("{:<5} | {:<12} | {:<12} | {:<12} | {:<12} | {:<12}",
i+1, widths[i], errL2[i], errH1[i], eocL2[i], eocH1[i]);
+ msg("elapsed time: {} seconds", t.elapsed());
return 0;
}
--
GitLab