From 4db9083a87a5f89efe46a042d46a375f10f8d4a1 Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Mon, 2 Jul 2018 16:06:12 +0200
Subject: [PATCH] Corrected errors in ConvectionDiffusionOperator and added
example
---
examples/CMakeLists.txt | 2 +-
examples/convection_diffusion.cc | 44 ++++++++++++++++
.../assembler/ConvectionDiffusionOperator.hpp | 51 ++++++++++---------
3 files changed, 73 insertions(+), 24 deletions(-)
create mode 100644 examples/convection_diffusion.cc
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 5a8d6233..ab96c269 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -1,6 +1,6 @@
add_custom_target(examples)
-set(projects2d "ellipt" "heat" "vecellipt" "stokes0" "stokes1" "stokes3" "navier_stokes")
+set(projects2d "ellipt" "heat" "vecellipt" "stokes0" "stokes1" "stokes3" "navier_stokes" "convection_diffusion")
foreach(project ${projects2d})
add_executable(${project}.2d ${project}.cc)
diff --git a/examples/convection_diffusion.cc b/examples/convection_diffusion.cc
new file mode 100644
index 00000000..40871822
--- /dev/null
+++ b/examples/convection_diffusion.cc
@@ -0,0 +1,44 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#include <iostream>
+
+#include <amdis/AMDiS.hpp>
+#include <amdis/ProblemStat.hpp>
+#include <amdis/assembler/ConvectionDiffusionOperator.hpp>
+#include <amdis/common/Literals.hpp>
+
+using namespace AMDiS;
+
+// 1 component with polynomial degree 1
+//using Grid = Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>;
+using ElliptParam = YaspGridBasis<AMDIS_DIM, 1>;
+using ElliptProblem = ProblemStat<ElliptParam>;
+
+int main(int argc, char** argv)
+{
+ AMDiS::init(argc, argv);
+
+ using namespace Dune::Indices;
+
+ ElliptProblem prob("ellipt");
+ prob.initialize(INIT_ALL);
+
+ // -div(A*grad(u)) + div(b*u) + c*u = f
+ auto opCD = convectionDiffusion(/*A=*/1.0, /*b=*/0.0, /*c=*/1.0, /*f=*/1.0);
+ prob.addMatrixOperator(opCD, _0, _0);
+ prob.addVectorOperator(opCD, _0);
+
+ // set boundary condition
+ auto predicate = [](auto const& x){ return x[0] < 1.e-8 || x[1] < 1.e-8; }; // define boundary
+ auto dbcValues = [](auto const& x){ return 0.0; }; // set value
+ prob.addDirichletBC(predicate, _0, _0, dbcValues);
+
+ AdaptInfo adaptInfo("adapt");
+ prob.buildAfterCoarsen(adaptInfo, Flag(0));
+ prob.solve(adaptInfo);
+ prob.writeFiles(adaptInfo, true);
+
+ AMDiS::finalize();
+ return 0;
+}
diff --git a/src/amdis/assembler/ConvectionDiffusionOperator.hpp b/src/amdis/assembler/ConvectionDiffusionOperator.hpp
index a2f938e7..388370ab 100644
--- a/src/amdis/assembler/ConvectionDiffusionOperator.hpp
+++ b/src/amdis/assembler/ConvectionDiffusionOperator.hpp
@@ -14,7 +14,7 @@ namespace AMDiS
* @{
**/
- /// convection-diffusion operator
+ /// convection-diffusion operator, see \ref convectionDiffusion
/// <A*grad(u),grad(v)> - <b*u, grad(v)> + <c*u, v> = <f, v> (conserving) or
/// <A*grad(u),grad(v)> + <b*grad(u), v> + <c*u, v> = <f, v> (non conserving)
template <class LocalContext, class GridFctA, class GridFctB, class GridFctC, class GridFctF, bool conserving = true>
@@ -71,7 +71,7 @@ namespace AMDiS
this->getDegree(1,coeffOrder(localFctB),rowNode,colNode),
this->getDegree(0,coeffOrder(localFctC),rowNode,colNode)});
- using QuadratureRules = Dune::QuadratureRules<typename Context::Geometry::ctype, Context::dimension>;
+ using QuadratureRules = Dune::QuadratureRules<typename Context::Geometry::ctype, Context::LocalContext::mydimension>;
auto const& quad = QuadratureRules::rule(context.type(), quadDeg);
for (auto const& qp : quad) {
@@ -93,30 +93,33 @@ namespace AMDiS
localFE.localBasis().evaluateJacobian(local, shapeGradients);
// Compute the shape function gradients on the real element
- thread_local std::vector<Dune::FieldVector<RangeFieldType,Context::dow> > gradients(shapeGradients.size());
+ using WorldVector = FieldVector<RangeFieldType,Context::dow>;
+ thread_local std::vector<WorldVector> gradients(shapeGradients.size());
for (std::size_t i = 0; i < gradients.size(); ++i)
jacobian.mv(shapeGradients[i][0], gradients[i]);
const auto A = localFctA(local);
- const auto b = localFctB(local);
+ WorldVector b; b = localFctB(local);
const auto c = localFctC(local);
IF_CONSTEXPR(conserving) {
+ WorldVector gradAi, gradBi;
for (std::size_t i = 0; i < localFE.size(); ++i) {
const int local_i = rowNode.localIndex(i);
- const auto gradAi = A * gradients[i];
- const auto gradBi = b * gradients[i];
+ gradAi = A * gradients[i];
+ gradBi = b * gradients[i];
for (std::size_t j = 0; j < localFE.size(); ++j) {
const int local_j = colNode.localIndex(j);
elementMatrix[local_i][local_j] += (dot(gradAi, gradients[j]) + (c * shapeValues[i] - gradBi) * shapeValues[j]) * factor;
}
}
} else {
+ WorldVector grad_i;
for (std::size_t i = 0; i < localFE.size(); ++i) {
const int local_i = rowNode.localIndex(i);
- const auto grad_i = A * gradients[i];
- grad_i += b * shapeValues[i];
+ grad_i = A * gradients[i];
+ grad_i+= b * shapeValues[i];
for (std::size_t j = 0; j < localFE.size(); ++j) {
const int local_j = colNode.localIndex(j);
elementMatrix[local_i][local_j] += (dot(grad_i, gradients[j]) + c * shapeValues[i] * shapeValues[j]) * factor;
@@ -136,10 +139,10 @@ namespace AMDiS
Node const& node,
ElementVector& elementVector)
{
- static_assert(Node::isLeaf
+ static_assert(Node::isLeaf,
"Operator can be applied to Leaf-Nodes only." );
- using RangeType = typename FiniteElementType_t<RowNode>::Traits::LocalBasisType::Traits::RangeType;
+ using RangeType = typename FiniteElementType_t<Node>::Traits::LocalBasisType::Traits::RangeType;
auto localFctF = localFunction(gridFctF_); localFctF.bind(context.element());
@@ -147,7 +150,7 @@ namespace AMDiS
int quad_order = this->getDegree(0,coeffOrder(localFctF),node);
- using QuadratureRules = Dune::QuadratureRules<typename Context::Geometry::ctype, Context::dimension>;
+ using QuadratureRules = Dune::QuadratureRules<typename Context::Geometry::ctype, Context::LocalContext::dimension>;
auto const& quad = QuadratureRules::rule(context.type(), quad_order);
for (auto const& qp : quad) {
@@ -164,8 +167,8 @@ namespace AMDiS
const auto f = localFctF(local);
for (std::size_t i = 0; i < localFE.size(); ++i) {
- const int local_i = rowNode.localIndex(i);
- elementVector[local_i][local_j] += f * shapeValues[i] * factor;
+ const int local_i = node.localIndex(i);
+ elementVector[local_i] += f * shapeValues[i] * factor;
}
}
@@ -180,22 +183,23 @@ namespace AMDiS
template <class LocalFct>
int coeffOrder(LocalFct const& localFct)
{
- using Concept = Dune::Std::is_detected<HasLocalFunctionOrder, LocalFunction>;
+ using Concept = Dune::Std::is_detected<HasLocalFunctionOrder, LocalFct>;
return Dune::Hybrid::ifElse(Concept{},
[&](auto id) { return order(id(localFct)); },
[] (auto) { return 0; });
}
private:
- GridFctA gridFctA;
- GridFctB gridFctB;
- GridFctC gridFctC;
- GridFctF gridFctF;
+ GridFctA gridFctA_;
+ GridFctB gridFctB_;
+ GridFctC gridFctC_;
+ GridFctF gridFctF_;
};
- template <class PreGridFctA, class PreGridFctB, class PreGridFctC, class PreGridFctF, bool conserving>
+ template <class PreGridFctA, class PreGridFctB, class PreGridFctC, class PreGridFctF, class c>
struct PreConvectionDiffusionOperator
{
+ static constexpr bool conserving = c::value;
PreGridFctA gridFctA;
PreGridFctB gridFctB;
PreGridFctC gridFctC;
@@ -207,20 +211,21 @@ namespace AMDiS
PreGridFctC const& gridFctC, PreGridFctF const& gridFctF,
bool_t<conserving> = {})
{
- using Pre = PreConvectionDiffusionOperator<PreGridFctA, PreGridFctB, PreGridFctC, PreGridFctF, conserving>;
+ using Pre = PreConvectionDiffusionOperator<PreGridFctA, PreGridFctB, PreGridFctC, PreGridFctF, bool_t<conserving>>;
return Pre{gridFctA, gridFctB, gridFctC, gridFctF};
}
- template <class LocalContext, class... GrdFcts, bool conserving, class GridView>
- auto makeLocalOperator(PreConvectionDiffusionOperator<GridFcts..., conserving> const& pre, GridView const& gridView)
+ template <class LocalContext, class... T, class GridView>
+ auto makeLocalOperator(PreConvectionDiffusionOperator<T...> const& pre, GridView const& gridView)
{
+ using Pre = PreConvectionDiffusionOperator<T...>;
auto gridFctA = makeGridFunction(pre.gridFctA, gridView);
auto gridFctB = makeGridFunction(pre.gridFctB, gridView);
auto gridFctC = makeGridFunction(pre.gridFctC, gridView);
auto gridFctF = makeGridFunction(pre.gridFctF, gridView);
using GridFctOp = ConvectionDiffusionOperator<LocalContext,
- decltype(gridFctA), decltype(gridFctB), decltype(gridFctC), decltype(gridFctF), conserving>;
+ decltype(gridFctA), decltype(gridFctB), decltype(gridFctC), decltype(gridFctF), Pre::conserving>;
GridFctOp localOperator{gridFctA, gridFctB, gridFctC, gridFctF};
return localOperator;
--
GitLab