1 #include <amdis/AMDiS.hpp> 3 #include <amdis/Integrate.hpp> 4 #include <amdis/LocalOperators.hpp> 5 #include <amdis/ProblemStat.hpp> 6 #include <amdis/common/Literals.hpp> 10 #include <dune/grid/yaspgrid.hh> 11 using Grid = Dune::YaspGrid<GRIDDIM>;
15 using namespace AMDiS;
19 int main(
int argc,
char** argv)
26 int numLevels = GRIDDIM == 2 ? 6 : 4;
28 numLevels = std::atoi(argv[2]);
30 auto f = [](
auto const& x) {
32 double ux = std::exp(-10.0 * r2);
33 return -(400.0 * r2 - 20.0 * x.size()) * ux;
35 auto g = [](
auto const& x){
return std::exp(-10.0 * dot(x,x)); };
36 auto grad_g = [](
auto const& x) -> FieldMatrix<double,1,GRIDDIM> {
37 return {-20.0 * std::exp(-10.0 * dot(x,x)) * x};
44 prob.initialize(INIT_ALL);
50 prob.addMatrixOperator(opL);
54 prob.addVectorOperator(opForce);
59 prob.addDirichletBC(1, g);
65 std::vector<double> errL2; errL2.reserve(numLevels);
66 std::vector<double> errH1; errH1.reserve(numLevels);
67 std::vector<double> widths; widths.reserve(numLevels);
71 for (
int l = 0; l < numLevels; ++l) {
77 for (
auto const& e : elements(prob.gridView(), Dune::Partitions::interior)) {
78 auto refElem = Dune::referenceElement<double,GRIDDIM>(e.type());
79 auto geo = e.geometry();
80 for (
int i = 0; i < refElem.size(GRIDDIM-1); ++i) {
81 auto v0 = geo.global(refElem.position(refElem.subEntity(i,GRIDDIM-1,0,GRIDDIM), GRIDDIM));
82 auto v1 = geo.global(refElem.position(refElem.subEntity(i,GRIDDIM-1,1,GRIDDIM), GRIDDIM));
83 h = std::max(h, (v0 - v1).two_norm());
90 prob.assemble(adaptInfo);
91 prob.solve(adaptInfo);
95 double errorL2 =
integrate(sqr(g - prob.solution()), prob.gridView(), 6);
96 errL2.push_back(std::sqrt(errorL2));
97 double errorH1 = errorL2 +
integrate(unary_dot(grad_g -
gradientOf(prob.solution())), prob.gridView(), 6);
98 errH1.push_back(std::sqrt(errorH1));
104 prob.writeFiles(adaptInfo);
107 msg(
"{:5} | {:12} | {:12} | {:12} | {:12} | {:12}",
108 "level",
"h",
"|u - u*|_L2",
"|u - u*|_H1",
"eoc_L2",
"eoc_H1");
109 msg_(
"{0:->7}{0:->15}{0:->15}{0:->15}{0:->15}{1:->14}",
"+",
"\n");
111 std::vector<double> eocL2(numLevels, 0.0), eocH1(numLevels, 0.0);
112 for (
int i = 1; i < numLevels; ++i) {
113 eocL2[i] = std::log(errL2[i]/errL2[i-1]) / std::log(widths[i]/widths[i-1]);
114 eocH1[i] = std::log(errH1[i]/errH1[i-1]) / std::log(widths[i]/widths[i-1]);
117 for (
int i = 0; i < numLevels; ++i)
118 msg(
"{:<5} | {:<12} | {:<12} | {:<12} | {:<12} | {:<12}",
119 i+1, widths[i], errL2[i], errH1[i], eocL2[i], eocH1[i]);
121 msg(
"elapsed time: {} seconds", t.elapsed());
Definition: SecondOrderGradTestGradTrial.hpp:19
auto makeOperator(Tag tag, Expr &&expr, QuadratureArgs &&... args)
Store tag and expression into a PreGridFunctionOperator to create a GridFunctionOperator.
Definition: GridFunctionOperator.hpp:220
Contains all classes needed for solving linear and non linear equation systems.
Definition: AdaptBase.hpp:6
Establishes an environment for sequential and parallel AMDiS programs.
Definition: Environment.hpp:19
ProblemStatTraits for a (composite) basis composed of lagrange bases of different degree...
Definition: ProblemStatTraits.hpp:111
Definition: ProblemStat.hpp:52
void msg(std::string const &str, Args &&... args)
print a message
Definition: Output.hpp:98
auto gradientOf(Expr const &expr)
Definition: DerivativeGridFunction.hpp:179
void msg_(std::string const &str, Args &&... args)
print a message (without appended newline)
Definition: Output.hpp:120
Definition: ZeroOrderTest.hpp:17
Holds adapt parameters and infos about the problem.
Definition: AdaptInfo.hpp:25
auto integrate(Expr &&expr, GridView const &gridView)
Integrate expression with quadrature rule determined by polynomial order of expression.
Definition: Integrate.hpp:49