Skip to content
Snippets Groups Projects
Commit a6617120 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

added ellipt demo corrections

parent f7cc7c36
No related branches found
No related tags found
1 merge request!107Restructuring of Matrix and Vector and linear-algebra backends
This commit is part of merge request !107. Comments created here will be created in the context of that merge request.
#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;
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment