Merge branch 'issue/convection_diffusion_operator' into 'develop'

Corrected errors in ConvectionDiffusionOperator and added example

See merge request spraetor/dune-amdis!25

examples/CMakeLists.txt

 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)

examples/convection_diffusion.cc

+// -*- 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;
+}

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;