Merge branch 'feature/first_order_operators' into 'develop'

added some missing operators and a test for all operators

See merge request spraetor/dune-amdis!67

src/amdis/Operators.hpp

@@ -46,6 +46,10 @@
 #include <amdis/assembler/ZeroOrderTestvecTrialvec.hpp>  // <Psi, A * Phi>, <Psi, c * Phi>
 
 // first-order operators
+#include <amdis/assembler/FirstOrderPartialTest.hpp>      // <d_i(psi), c>
+#include <amdis/assembler/FirstOrderGradTest.hpp>         // <grad(psi), b>
+#include <amdis/assembler/FirstOrderDivTestvec.hpp>       // <div(Psi), c>
+
 #include <amdis/assembler/FirstOrderDivTestvecTrial.hpp>   // <div(Psi), c * phi>
 #include <amdis/assembler/FirstOrderGradTestTrial.hpp>     // <grad(psi), b * phi>
 #include <amdis/assembler/FirstOrderGradTestTrialvec.hpp>  // <grad(psi), c * Phi>

src/amdis/assembler/ConvectionDiffusionOperator.hpp

@@ -108,11 +108,11 @@ namespace AMDiS
         const auto c = localFctC(local);
 
         IF_CONSTEXPR(conserving) {
-          WorldVector gradAi, gradBi;
+          WorldVector gradAi;
           for (std::size_t i = 0; i < numLocalFe; ++i) {
             const int local_i = rowNode.localIndex(i);
             gradAi = A * gradients[i];
-            gradBi = b * gradients[i];
+            auto gradBi = b * gradients[i];
             for (std::size_t j = 0; j < numLocalFe; ++j) {
               const int local_j = colNode.localIndex(j);
               elementMatrix[local_i][local_j] += (dot(gradAi, gradients[j]) + (c * shapeValues[i] - gradBi) * shapeValues[j]) * factor;

src/amdis/assembler/FirstOrderDivTestvec.hpp

+#pragma once
+
+#include <type_traits>
+
+#include <amdis/GridFunctionOperator.hpp>
+#include <amdis/LocalBasisCache.hpp>
+#include <amdis/Output.hpp>
+#include <amdis/common/Size.hpp>
+
+namespace AMDiS
+{
+  /**
+   * \addtogroup operators
+   * @{
+   **/
+
+  namespace tag
+  {
+    struct divtestvec {};
+  }
+
+
+  /// first-order operator \f$ \langle\nabla\cdot\Psi, c\rangle \f$
+  template <class LocalContext, class GridFct>
+  class GridFunctionOperator<tag::divtestvec, LocalContext, GridFct>
+      : public GridFunctionOperatorBase<GridFunctionOperator<tag::divtestvec, LocalContext, GridFct>,
+                                        LocalContext, GridFct>
+  {
+    using Self = GridFunctionOperator;
+    using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
+
+    static_assert( Size_v<typename GridFct::Range> == 1, "Expression must be of scalar type." );
+
+  public:
+    GridFunctionOperator(tag::divtestvec, GridFct const& expr)
+      : Super(expr, 1)
+    {}
+
+    template <class Context, class Node, class ElementVector>
+    void getElementVector(Context const& context,
+                          Node const& node,
+                          ElementVector& elementVector)
+    {
+      static_assert(Node::isPower, "Node must be a Power-Node.");
+
+      static const std::size_t CHILDREN = Node::CHILDREN;
+      static_assert( CHILDREN == Context::dow, "" );
+
+      auto const& quad = this->getQuadratureRule(context.type(), node);
+      auto const& localFE = node.child(0).finiteElement();
+      std::size_t feSize = localFE.size();
+
+      using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
+      using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;
+
+      LocalBasisCache<LocalBasisType> cache(localFE.localBasis());
+
+      auto const& shapeGradientsCache = cache.evaluateJacobianAtQp(context, quad);
+      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+        // Position of the current quadrature point in the reference element
+        decltype(auto) local = context.local(quad[iq].position());
+
+        // The transposed inverse Jacobian of the map from the reference element to the element
+        const auto jacobian = context.geometry().jacobianInverseTransposed(local);
+
+        // The multiplicative factor in the integral transformation formula
+        const auto factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
+
+        // The gradients of the shape functions on the reference element
+        auto const& shapeGradients = shapeGradientsCache[iq];
+
+        // Compute the shape function gradients on the real element
+        using WorldVector = FieldVector<RangeFieldType,Context::dow>;
+        thread_local std::vector<WorldVector> gradients;
+        gradients.resize(shapeGradients.size());
+
+        for (std::size_t i = 0; i < gradients.size(); ++i)
+          jacobian.mv(shapeGradients[i][0], gradients[i]);
+
+        for (std::size_t j = 0; j < feSize; ++j) {
+          for (std::size_t k = 0; k < CHILDREN; ++k) {
+            const auto local_kj = node.child(k).localIndex(j);
+            elementVector[local_kj] += factor * gradients[j][k];
+          }
+        }
+      }
+
+    }
+  };
+
+  /** @} **/
+
+} // end namespace AMDiS

src/amdis/assembler/FirstOrderGradTest.hpp

+#pragma once
+
+#include <type_traits>
+
+#include <amdis/GridFunctionOperator.hpp>
+#include <amdis/LocalBasisCache.hpp>
+#include <amdis/Output.hpp>
+#include <amdis/common/Size.hpp>
+
+namespace AMDiS
+{
+  /**
+   * \addtogroup operators
+   * @{
+   **/
+
+  namespace tag
+  {
+    struct gradtest {};
+  }
+
+
+  /// first-order operator \f$ \langle\nabla\psi, b\rangle \f$
+  template <class LocalContext, class GridFct>
+  class GridFunctionOperator<tag::gradtest, LocalContext, GridFct>
+      : public GridFunctionOperatorBase<GridFunctionOperator<tag::gradtest, LocalContext, GridFct>,
+                                        LocalContext, GridFct>
+  {
+    static const int dow = ContextGeometry<LocalContext>::dow;
+    using Self = GridFunctionOperator;
+    using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
+
+    static_assert( Size_v<typename GridFct::Range> == dow, "Expression must be of vector type." );
+
+  public:
+    GridFunctionOperator(tag::gradtest, GridFct const& expr)
+      : Super(expr, 1)
+    {}
+
+    template <class Context, class Node, class ElementVector>
+    void getElementVector(Context const& context,
+                          Node const& node,
+                          ElementVector& elementVector)
+    {
+      static_assert(Node::isLeaf, "Node must be Leaf-Node.");
+
+      auto const& quad = this->getQuadratureRule(context.type(), node);
+      auto const& localFE = node.finiteElement();
+      std::size_t feSize = localFE.size();
+
+      using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
+      using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;
+
+      LocalBasisCache<LocalBasisType> cache(localFE.localBasis());
+
+      auto const& shapeGradientsCache = cache.evaluateJacobianAtQp(context, quad);
+      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+        // Position of the current quadrature point in the reference element
+        decltype(auto) local = context.local(quad[iq].position());
+
+        // The transposed inverse Jacobian of the map from the reference element to the element
+        const auto jacobian = context.geometry().jacobianInverseTransposed(local);
+
+        // The multiplicative factor in the integral transformation formula
+        const auto factor = Super::coefficient(local);
+        const auto dx = context.integrationElement(quad[iq].position()) * quad[iq].weight();
+
+        // The gradients of the shape functions on the reference element
+        auto const& shapeGradients = shapeGradientsCache[iq];
+
+        // Compute the shape function gradients on the real element
+        using WorldVector = FieldVector<RangeFieldType,Context::dow>;
+        thread_local std::vector<WorldVector> gradients;
+        gradients.resize(shapeGradients.size());
+
+        for (std::size_t i = 0; i < gradients.size(); ++i)
+          jacobian.mv(shapeGradients[i][0], gradients[i]);
+
+        for (std::size_t i = 0; i < feSize; ++i) {
+          const auto local_i = node.localIndex(i);
+          elementVector[local_i] += dx * (factor * gradients[i]);
+        }
+      }
+
+    }
+  };
+
+  /** @} **/
+
+} // end namespace AMDiS

src/amdis/assembler/FirstOrderPartialTest.hpp

+#pragma once
+
+#include <type_traits>
+
+#include <amdis/assembler/FirstOrderTestPartialTrial.hpp>
+
+namespace AMDiS
+{
+  /**
+   * \addtogroup operators
+   * @{
+   **/
+
+  namespace tag
+  {
+    struct partialtest
+    {
+      std::size_t comp;
+    };
+  }
+
+
+  /// first-order operator \f$ \langle\partial_i\psi, c\rangle \f$
+
+  template <class LocalContext, class GridFct>
+  class GridFunctionOperator<tag::partialtest, LocalContext, GridFct>
+      : public GridFunctionOperatorBase<GridFunctionOperator<tag::partialtest, LocalContext, GridFct>,
+                                        LocalContext, GridFct>
+  {
+    using Self = GridFunctionOperator;
+    using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
+
+    static_assert( Size_v<typename GridFct::Range> == 1, "Expression must be of scalar type." );
+
+  public:
+    GridFunctionOperator(tag::partialtest tag, GridFct const& expr)
+      : Super(expr, 1)
+      , comp_(tag.comp)
+    {}
+
+    template <class Context, class Node, class ElementVector>
+    void getElementVector(Context const& context,
+                          Node const& node,
+                          ElementVector& elementVector)
+    {
+      static_assert(Node::isLeaf, "Operator can be applied to Leaf-Nodes only.");
+
+      auto const& quad = this->getQuadratureRule(context.type(), node);
+      auto const& localFE = node.finiteElement();
+      std::size_t feSize = localFE.size();
+
+      using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
+      using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;
+
+      LocalBasisCache<LocalBasisType> cache(localFE.localBasis());
+
+      auto const& shapeGradientsCache = cache.evaluateJacobianAtQp(context, quad);
+      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+        // Position of the current quadrature point in the reference element
+        decltype(auto) local = context.local(quad[iq].position());
+
+        // The transposed inverse Jacobian of the map from the reference element to the element
+        const auto jacobian = context.geometry().jacobianInverseTransposed(local);
+        assert(jacobian.M() == Context::dim);
+
+        // The multiplicative factor in the integral transformation formula
+        const auto dx = context.integrationElement(quad[iq].position()) * quad[iq].weight();
+        const auto c = Super::coefficient(local);
+
+        // The gradients of the shape functions on the reference element
+        auto const& shapeGradients = shapeGradientsCache[iq];
+
+        // Compute the shape function gradients on the real element
+        thread_local std::vector<RangeFieldType> partial;
+        partial.resize(shapeGradients.size());
+        for (std::size_t i = 0; i < partial.size(); ++i) {
+          partial[i] = jacobian[comp_][0] * shapeGradients[i][0][0];
+          for (std::size_t j = 1; j < jacobian.M(); ++j)
+            partial[i] += jacobian[comp_][j] * shapeGradients[i][0][j];
+        }
+
+        for (std::size_t j = 0; j < feSize; ++j) {
+          const auto local_j = node.localIndex(j);
+          elementVector[local_j] += dx * (c * partial[j]);
+        }
+      }
+    }
+
+  private:
+    std::size_t comp_;
+  };
+
+  /** @} **/
+
+} // end namespace AMDiS

src/amdis/assembler/SecondOrderPartialTestPartialTrial.hpp

@@ -88,7 +88,7 @@ namespace AMDiS
         rowPartial.resize(rowShapeGradients.size());
         for (std::size_t i = 0; i < rowPartial.size(); ++i) {
           rowPartial[i] = jacobian[compTest_][0] * rowShapeGradients[i][0][0];
-          for (std::size_t j = 1; j < jacobian.cols(); ++j)
+          for (std::size_t j = 1; j < jacobian.cols; ++j)
             rowPartial[i] += jacobian[compTest_][j] * rowShapeGradients[i][0][j];
         }
 
@@ -96,7 +96,7 @@ namespace AMDiS
         colPartial.resize(colShapeGradients.size());
         for (std::size_t i = 0; i < colPartial.size(); ++i) {
           colPartial[i] = jacobian[compTrial_][0] * colShapeGradients[i][0][0];
-          for (std::size_t j = 1; j < jacobian.cols(); ++j)
+          for (std::size_t j = 1; j < jacobian.cols; ++j)
             colPartial[i] += jacobian[compTrial_][j] * colShapeGradients[i][0][j];
         }
 

test/CMakeLists.txt

@@ -40,6 +40,9 @@ dune_add_test(SOURCES MultiTypeVectorTest.cpp
 dune_add_test(SOURCES MultiTypeMatrixTest.cpp
   LINK_LIBRARIES amdis)
 
+dune_add_test(SOURCES OperatorsTest.cpp
+  LINK_LIBRARIES amdis)
+
 dune_add_test(SOURCES RangeTypeTest.cpp
   LINK_LIBRARIES amdis)
 

test/OperatorsTest.cpp

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#include "config.h"
+#include <amdis/AMDiS.hpp>
+#include <amdis/ProblemStat.hpp>
+#include <amdis/Operators.hpp>
+#include <amdis/assembler/ConvectionDiffusionOperator.hpp>
+#include <amdis/assembler/StokesOperator.hpp>
+
+
+using namespace AMDiS;
+
+
+using Grid = Dune::YaspGrid<2>;
+using Param   = TaylorHoodBasis<typename Grid::LeafGridView>;
+using Problem = ProblemStat<Param>;
+
+int main(int argc, char** argv)
+{
+  AMDiS::init(argc, argv);
+
+  using namespace Dune::Indices;
+
+  Problem prob("ellipt");
+  prob.initialize(INIT_ALL);
+
+  auto _v = _0;
+  auto _p = _1;
+
+  FieldVector<double,1> scalar = 1.0;
+  FieldVector<double,2> vec; vec = 1.0;
+  FieldMatrix<double,2,2> mat; mat = 1.0;
+
+  auto op1 = makeOperator(tag::divtestvec{}, scalar);
+  prob.addVectorOperator(op1, _v);
+
+  auto op2 = makeOperator(tag::divtestvec_trial{}, scalar);
+  prob.addMatrixOperator(op2, _v, _p);
+
+  auto op3 = makeOperator(tag::gradtest{}, vec);
+  prob.addVectorOperator(op3, _p);
+
+  auto op4 = makeOperator(tag::gradtest_trial{}, vec);
+  prob.addMatrixOperator(op4, _p, _p);
+
+  auto op5a = makeOperator(tag::gradtest_trialvec{}, scalar);
+  // auto op5b = makeOperator(tag::gradtest_trialvec{}, mat); // TODO: needs to be implemented
+  prob.addMatrixOperator(op5a, _p, _v);
+
+  auto op6 = makeOperator(tag::partialtest{0}, scalar);
+  prob.addVectorOperator(op6, _p);
+
+  auto op7 = makeOperator(tag::partialtest_trial{0}, scalar);
+  prob.addMatrixOperator(op7, _p, _p);
+
+  auto op8 = makeOperator(tag::test_divtrialvec{}, scalar);
+  prob.addMatrixOperator(op8, _p, _v);
+
+  auto op9 = makeOperator(tag::test_gradtrial{}, vec);
+  prob.addMatrixOperator(op9, _p, _p);
+
+  auto op10 = makeOperator(tag::test_partialtrial{0}, scalar);
+  prob.addMatrixOperator(op10, _p, _p);
+
+  auto op11a = makeOperator(tag::testvec_gradtrial{}, scalar);
+  // auto op11b = makeOperator(tag::testvec_gradtrial{}, mat); // TODO: needs to be implemented
+  prob.addMatrixOperator(op11a, _v, _p);
+
+  auto op12 = makeOperator(tag::divtestvec_divtrialvec{}, scalar);
+  prob.addMatrixOperator(op12, _v, _v);
+
+  auto op13a = makeOperator(tag::gradtest_gradtrial{}, scalar);
+  auto op13b = makeOperator(tag::gradtest_gradtrial{}, mat);
+  prob.addMatrixOperator(op13a, _p, _p);
+  prob.addMatrixOperator(op13b, _p, _p);
+
+  auto op14 = makeOperator(tag::partialtest_partialtrial{0,0}, scalar);
+  prob.addMatrixOperator(op14, _p, _p);
+
+  auto op15 = makeOperator(tag::test{}, scalar);
+  prob.addVectorOperator(op15, _p);
+
+  auto op16 = makeOperator(tag::test_trial{}, scalar);
+  prob.addMatrixOperator(op16, _p, _p);
+
+  auto op17 = makeOperator(tag::test_trialvec{}, vec);
+  prob.addMatrixOperator(op17, _p, _v);
+
+  auto op18 = makeOperator(tag::testvec{}, vec);
+  prob.addVectorOperator(op18, _v);
+
+  auto op19 = makeOperator(tag::testvec_trial{}, vec);
+  prob.addMatrixOperator(op19, _v, _p);
+
+  auto op20a = makeOperator(tag::testvec_trialvec{}, scalar);
+  auto op20b = makeOperator(tag::testvec_trialvec{}, mat);
+  prob.addMatrixOperator(op20a, _v, _v);
+  prob.addMatrixOperator(op20b, _v, _v);
+
+  // some special operators
+
+  auto opStokes = makeOperator(tag::stokes{}, scalar);
+  prob.addMatrixOperator(opStokes);
+  prob.addMatrixOperator(opStokes, {}, {});
+
+  auto opCDa = convectionDiffusion(scalar, scalar, scalar, scalar);
+  prob.addMatrixOperator(opCDa, _p, _p);
+  prob.addVectorOperator(opCDa, _p);
+
+  auto opCDb = convectionDiffusion(mat, vec, scalar, scalar, std::false_type{});
+  prob.addMatrixOperator(opCDb, _p, _p);
+  prob.addVectorOperator(opCDb, _p);
+
+  AMDiS::finalize();
+  return 0;
+}