boundaryterms nearly finished

dune/amdis/Assembler.hpp

+#pragma once
+
+#include <memory>
+#include <dune/amdis/QuadratureRules.hpp>
+#include <dune/amdis/common/TypeDefs.hpp>
+
+namespace AMDiS
+{
+  template <class GridView, class Element>
+  class Assembler
+  {
+    static constexpr int dim = GridView::dimension;
+    using ctype = typename GridView::ctype;
+
+    using QuadratureRule = QuadratureRuleFactory_t<Element, ctype, dim>;
+    using Geometry = typename Impl::Get<Element>::Geometry;
+
+  public:
+    template <class Operator>
+    Assembler(Operator& op, Element const& element, int degree, FirstOrderType type = GRD_PHI)
+      : geometry(get_geometry(element))
+    {
+      int order = op.getQuadratureDegree(geometry.type(), geometry, degree, type);
+      auto const& quad = Dune::QuadratureRules<ctype, Element::mydimension>::rule(geometry.type(), order);
+
+      quadrature.reset(new QuadratureRule(element, quad));
+    }
+
+    auto const& getQuadrature() const
+    {
+      return quadrature->getRule();
+    }
+
+    auto const& getTransformedQuadrature() const
+    {
+      return *quadrature;
+    }
+
+    Geometry const& getGeometry() const
+    {
+      return geometry;
+    }
+
+  private:
+    std::shared_ptr<QuadratureRule const> quadrature;
+    Geometry geometry; // TODO: create geometry just once for each element, e.g. in ProblemStat when traversing the grid
+  };
+
+} // end namespace AMDiS

dune/amdis/CMakeLists.txt

@@ -15,6 +15,7 @@ dune_add_library("duneamdis" NO_EXPORT
     linear_algebra/istl/SystemVector.cpp
     linear_algebra/mtl/SystemMatrix.cpp
     linear_algebra/mtl/SystemVector.cpp
+    utility/Filesystem.cpp
     )
 add_dune_alberta_flags("duneamdis" OBJECT USE_GENERIC)
 target_compile_definitions("duneamdis" PUBLIC AMDIS_BACKEND_MTL=1)

dune/amdis/Concepts.hpp

+#pragma once
+
+#include <dune/common/concept.hh>
+
+namespace AMDiS
+{
+  namespace Concepts
+  {
+    namespace Definition
+    {
+      // Concept for a Dune::Entity
+      struct Entity
+      {
+        template<class E>
+        auto require(const E& e) -> decltype(
+          e.level(),
+          e.partitionType(),
+          e.geometry(),
+          e.type(),
+          e.subEntities((unsigned int)(0)),
+          e.seed()
+        );
+      };
+
+      // Concept for a Dune::Intersection
+      struct Intersection
+      {
+        template<class I>
+        auto require(const I& i) -> decltype(
+          i.boundary(),
+          i.neighbor(),
+          i.inside(),
+          i.outside(),
+          i.geometryInInside(),
+          i.geometryInOutside(),
+          i.geometry(),
+          i.type(),
+          i.indexInInside(),
+          i.indexInOutside()
+        );
+      };
+
+    } // end namespace Definition
+
+    template <class E>
+    constexpr bool Entity = Dune::models<Definition::Entity, E>();
+
+    template <class I>
+    constexpr bool Intersection = Dune::models<Definition::Intersection, I>();
+
+  } // end namespace Concepts
+
+} // end namespace AMDiS

dune/amdis/FileWriter.hpp

@@ -11,6 +11,7 @@
 
 #include <dune/amdis/Initfile.hpp>
 #include <dune/amdis/common/Loops.hpp>
+#include <dune/amdis/utility/Filesystem.hpp>
 
 namespace AMDiS
 {
@@ -45,8 +46,11 @@ namespace AMDiS
       dir = "output";
       Parameters::get(base + "->output directory", dir);
 
-      vtkWriter = make_shared<Dune::VTKWriter<MeshView>>(meshView);
-      vtkSeqWriter = make_shared<Dune::VTKSequenceWriter<MeshView>>(vtkWriter, filename, dir, "");
+      if (!filesystem::exists(dir))
+        error_exit("Output directory '",dir,"' does not exist!");
+
+      vtkWriter = std::make_shared<Dune::VTKWriter<MeshView>>(meshView);
+      vtkSeqWriter = std::make_shared<Dune::VTKSequenceWriter<MeshView>>(vtkWriter, filename, dir, "");
     }
 
 

dune/amdis/FirstOrderAssembler.hpp

+#pragma once
+
+#include <vector>
+
+#include <dune/amdis/Assembler.hpp>
+#include <dune/amdis/common/Mpl.hpp>
+#include <dune/amdis/common/TypeDefs.hpp>
+#include <dune/amdis/utility/GetEntity.hpp>
+
+namespace AMDiS
+{
+  template <class GridView, class Element, FirstOrderType type>
+  class FirstOrderAssembler
+      : public Assembler<GridView, Element>
+  {
+    using Super = Assembler<GridView, Element>;
+
+    static constexpr int dim = GridView::dimension;
+    static constexpr int dow = GridView::dimensionworld;
+
+    using ElementMatrix = Impl::ElementMatrix;
+    using ElementVector = Impl::ElementVector;
+
+  public:
+    template <class Operator>
+    FirstOrderAssembler(Operator& op, Element const& element)
+      : Super(op, element, 1, type)
+    {
+      if (type == GRD_PHI)
+        op.initFot1(element, Super::getTransformedQuadrature());
+      else
+        op.initFot2(element, Super::getTransformedQuadrature());
+    }
+
+
+    template <class Operator, class RowView, class ColView>
+    void calculateElementMatrix(Operator& op, RowView const& rowView, ColView const& colView,
+                                ElementMatrix& elementMatrix)
+    {
+      calculateElementMatrix(op, rowView, colView, elementMatrix, std::integral_constant<FirstOrderType, type>{});
+    }
+
+    template <class Operator, class RowView>
+    void calculateElementVector(Operator& op, RowView const& rowView, ElementVector& elementVector)
+    {
+      calculateElementVector(op, rowView, elementVector, std::integral_constant<FirstOrderType, type>{});
+    }
+
+
+    template <class Operator, class RowView, class ColView>
+    void calculateElementMatrix(Operator& op,
+                                RowView const& rowView,
+                                ColView const& colView,
+                                ElementMatrix& elementMatrix,
+                                std::integral_constant<FirstOrderType, GRD_PHI>)
+    {
+      auto geometry = rowView.element().geometry();
+      auto const& quad_geometry = Super::getGeometry();
+
+      auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
+      auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
+      auto const& quad = Super::getQuadrature();
+
+      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+        // Position of the current quadrature point in the reference element
+        decltype(auto) quadPos = get_position<Element>(quad_geometry, quad[iq].position());
+
+        // The transposed inverse Jacobian of the map from the reference element to the element
+        const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
+
+        // The multiplicative factor in the integral transformation formula
+        const double factor = quad_geometry.integrationElement(quad[iq].position()) * quad[iq].weight();
+
+        std::vector<Dune::FieldVector<double,1> > rowShapeValues;
+        rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
+
+        // The gradients of the shape functions on the reference element
+        std::vector<Dune::FieldMatrix<double,1,dim> > colReferenceGradients;
+        colLocalBasis.evaluateJacobian(quadPos, colReferenceGradients);
+
+        // Compute the shape function gradients on the real element
+        std::vector<Dune::FieldVector<double,dow> > colGradients(colReferenceGradients.size());
+
+        for (std::size_t i = 0; i < colGradients.size(); ++i)
+          jacobian.mv(colReferenceGradients[i][0], colGradients[i]);
+
+        for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) {
+          for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) {
+            const int local_i = rowView.tree().localIndex(i);
+            const int local_j = colView.tree().localIndex(j);
+            op.evalFot1(elementMatrix[local_i][local_j], iq, factor, rowShapeValues[i], colGradients[j]);
+          }
+        }
+      }
+    }
+
+
+    template <class Operator, class RowView, class ColView>
+    void calculateElementMatrix(Operator& op,
+                                RowView const& rowView,
+                                ColView const& colView,
+                                ElementMatrix& elementMatrix,
+                                std::integral_constant<FirstOrderType, GRD_PSI>)
+    {
+      auto geometry = rowView.element().geometry();
+      auto const& quad_geometry = Super::getGeometry();
+
+      auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
+      auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
+      auto const& quad = Super::getQuadrature();
+
+      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+        // Position of the current quadrature point in the reference element
+        decltype(auto) quadPos = get_position<Element>(quad_geometry, quad[iq].position());
+
+        // The transposed inverse Jacobian of the map from the reference element to the element
+        const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
+
+        // The multiplicative factor in the integral transformation formula
+        const double factor = quad_geometry.integrationElement(quad[iq].position()) * quad[iq].weight();
+
+        // The gradients of the shape functions on the reference element
+        std::vector<Dune::FieldMatrix<double,1,dim> > rowReferenceGradients;
+        rowLocalBasis.evaluateJacobian(quadPos, rowReferenceGradients);
+
+        // Compute the shape function gradients on the real element
+        std::vector<Dune::FieldVector<double,dow> > rowGradients(rowReferenceGradients.size());
+
+        for (std::size_t i = 0; i < rowGradients.size(); ++i)
+          jacobian.mv(rowReferenceGradients[i][0], rowGradients[i]);
+
+        std::vector<Dune::FieldVector<double,1> > colShapeValues;
+        colLocalBasis.evaluateFunction(quadPos, colShapeValues);
+
+        for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) {
+          for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) {
+            const int local_i = rowView.tree().localIndex(i);
+            const int local_j = colView.tree().localIndex(j);
+            op.evalFot2(elementMatrix[local_i][local_j], iq, factor, rowGradients[i], colShapeValues[j]);
+          }
+        }
+      }
+    }
+
+
+    template <class Operator, class RowView>
+    void calculateElementVector(Operator& op,
+                                RowView const& rowView,
+                                ElementVector& elementVector,
+                                std::integral_constant<FirstOrderType, GRD_PSI>)
+    {
+      auto geometry = rowView.element().geometry();
+      auto const& quad_geometry = Super::getGeometry();
+
+      auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
+      auto const& quad = Super::getQuadrature();
+
+      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+        // Position of the current quadrature point in the reference element
+        decltype(auto) quadPos = get_position<Element>(quad_geometry, quad[iq].position());
+
+        // The transposed inverse Jacobian of the map from the reference element to the element
+        const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
+
+        // The multiplicative factor in the integral transformation formula
+        const double factor = quad_geometry.integrationElement(quad[iq].position()) * quad[iq].weight();
+
+        // The gradients of the shape functions on the reference element
+        std::vector<Dune::FieldMatrix<double,1,dim> > rowReferenceGradients;
+        rowLocalBasis.evaluateJacobian(quadPos, rowReferenceGradients);
+
+        // Compute the shape function gradients on the real element
+        std::vector<Dune::FieldVector<double,dow> > rowGradients(rowReferenceGradients.size());
+
+        for (std::size_t i = 0; i < rowGradients.size(); ++i)
+          jacobian.mv(rowReferenceGradients[i][0], rowGradients[i]);
+
+        for (std::size_t i = 0; i < size(elementVector); ++i) {
+          const int local_i = rowView.tree().localIndex(i);
+          op.evalFot2(elementVector[local_i], iq, factor, rowGradients[i]);
+        }
+      }
+    }
+  };
+
+
+
+} // end namespace AMDiS

dune/amdis/Mesh.hpp

@@ -19,7 +19,7 @@
 
 #include <dune/amdis/Initfile.hpp>
 #include <dune/amdis/Output.hpp>
-#include <dune/amdis/common/Filesystem.hpp>
+#include <dune/amdis/utility/Filesystem.hpp>
 
 namespace AMDiS
 {
@@ -112,7 +112,7 @@ namespace AMDiS
       Parameters::get(meshName + "->macro file name", filename);
 
       if (!filename.empty()) {
-        Path fn(filename);
+        filesystem::path fn(filename);
         auto ext = fn.extension();
 
 #if HAVE_ALBERTA

dune/amdis/Operator.hpp

@@ -6,6 +6,9 @@
 #include <dune/geometry/type.hh>
 
 #include <dune/amdis/OperatorTermBase.hpp>
+#include <dune/amdis/ZeroOrderAssembler.hpp>
+#include <dune/amdis/FirstOrderAssembler.hpp>
+#include <dune/amdis/SecondOrderAssembler.hpp>
 #include <dune/amdis/terms/TermGenerator.hpp>
 
 #include <dune/amdis/common/TypeDefs.hpp>
@@ -13,24 +16,21 @@
 
 namespace AMDiS
 {
-  enum FirstOrderType
-  {
-    GRD_PSI,
-    GRD_PHI
-  };
-
-
-  template <class MeshView>
+  template <class MeshView, class Element>
   class Operator
   {
     using Self = Operator;
-    using OperatorTermType = OperatorTerm<MeshView>;
+    using OperatorTermType = OperatorTerm<Element>;
 
     using ElementVector = Impl::ElementVector;
     using ElementMatrix = Impl::ElementMatrix;
 
     using IdType = typename MeshView::Grid::LocalIdSet::IdType;
 
+    template <class, class> friend class ZeroOrderAssembler;
+    template <class, class, FirstOrderType> friend class FirstOrderAssembler;
+    template <class, class> friend class SecondOrderAssembler;
+
   public:
     /// Add coefficients for zero-order operator < c * u, v >.
     /// The coefficient \p c must be a scalar expression.
@@ -124,62 +124,23 @@ namespace AMDiS
     /// This takes into account the polynomial degree of trial and test functions
     /// and the polynomial degree of the coefficient functions
     template <class Geometry>
-    int getQuadratureDegree(Dune::GeometryType t, Geometry const& geometry, int order, FirstOrderType firstOrderType = GRD_PHI);
+    int getQuadratureDegree(Dune::GeometryType t, Geometry const& geometry, int order,
+                            FirstOrderType firstOrderType = GRD_PHI);
 
 
     template <class RowView, class ColView>
-    bool getElementMatrix(RowView const& rowView,
+    bool getElementMatrix(Element const& element,
+                          RowView const& rowView,
                           ColView const& colView,
                           ElementMatrix& elementMatrix,
                           double* factor = NULL);
 
-    template <class Intersection, class RowView, class ColView>
-    bool getBoundaryElementMatrix(Intersection const& intersection,
-                                  RowView const& rowView,
-                                  ColView const& colView,
-                                  ElementMatrix& elementMatrix,
-                                  double* factor = NULL);
-
     template <class RowView>
-    bool getElementVector(RowView const& rowView,
+    bool getElementVector(Element const& element,
+                          RowView const& rowView,
                           ElementVector& elementVector,
                           double* factor = NULL);
 
-    template <class Intersection, class RowView>
-    bool getBoundaryElementVector(Intersection const& intersection,
-                                  RowView const& rowView,
-                                  ElementVector& elementVector,
-                                  double* factor = NULL);
-
-  protected:
-    template <class RowView, class ColView>
-    void assembleZeroOrderTerms(RowView const& rowView,
-				ColView const& colView,
-				ElementMatrix& elementMatrix);
-
-    template <class RowView>
-    void assembleZeroOrderTerms(RowView const& rowView,
-				ElementVector& elementVector);
-
-    template <class RowView, class ColView>
-    void assembleFirstOrderTermsGrdPhi(RowView const& rowView,
-				       ColView const& colView,
-				       ElementMatrix& elementMatrix);
-
-    template <class RowView, class ColView>
-    void assembleFirstOrderTermsGrdPsi(RowView const& rowView,
-				       ColView const& colView,
-				       ElementMatrix& elementMatrix);
-
-    template <class RowView>
-    void assembleFirstOrderTermsGrdPsi(RowView const& rowView,
-                                       ElementVector& elementVector);
-
-    template <class RowView, class ColView>
-    void assembleSecondOrderTerms(RowView const& rowView,
-				  ColView const& colView,
-				  ElementMatrix& elementMatrix);
-
   private:
     template <class Term>
     Self& addZOTImpl(Term const& term);
@@ -215,6 +176,65 @@ namespace AMDiS
       return sot(std::forward<Args>(args)...);
     }
 
+  protected:
+
+    template <class... Args>
+    void evalZot(double& result, int iq, double factor, Args const&... args)
+    {
+      for (auto* operatorTerm : zeroOrder)
+        result += operatorTerm->evalZot(iq, args...) * factor;
+    }
+
+    template <class... Args>
+    void evalFot1(double& result, int iq, double factor, Args const&... args)
+    {
+      for (auto* operatorTerm : firstOrderGrdPhi)
+        result += operatorTerm->evalFot1(iq, args...) * factor;
+    }
+
+    template <class... Args>
+    void evalFot2(double& result, int iq, double factor, Args const&... args)
+    {
+      for (auto* operatorTerm : firstOrderGrdPsi)
+        result += operatorTerm->evalFot2(iq, args...) * factor;
+    }
+
+    template <class... Args>
+    void evalSot(double& result, int iq, double factor, Args const&... args)
+    {
+      for (auto* operatorTerm : secondOrder)
+        result += operatorTerm->evalSot(iq, args...) * factor;
+    }
+
+
+    template <class E, class Quadrature>
+    void initZot(E const& element, Quadrature const& quad)
+    {
+      for (auto* operatorTerm : zeroOrder)
+        operatorTerm->init(element, quad);
+    }
+
+    template <class E, class Quadrature>
+    void initFot1(E const& element, Quadrature const& quad)
+    {
+      for (auto* operatorTerm : firstOrderGrdPhi)
+        operatorTerm->init(element, quad);
+    }
+
+    template <class E, class Quadrature>
+    void initFot2(E const& element, Quadrature const& quad)
+    {
+      for (auto* operatorTerm : firstOrderGrdPsi)
+        operatorTerm->init(element, quad);
+    }
+
+    template <class E, class Quadrature>
+    void initSot(E const& element, Quadrature const& quad)
+    {
+      for (auto* operatorTerm : secondOrder)
+        operatorTerm->init(element, quad);
+    }
+
   private:
     template <class LocalView>
     IdType getElementId(LocalView const& localView);

dune/amdis/OperatorTermBase.hpp

@@ -7,30 +7,27 @@
 #include <vector>
 #include <type_traits>
 
-#include <dune/geometry/quadraturerules.hh>
+#include <dune/amdis/QuadratureRules.hpp>
 
 #include <dune/amdis/OperatorEvaluation.hpp>
 #include <dune/amdis/common/ValueCategory.hpp>
+#include <dune/amdis/utility/GetEntity.hpp>
 
 namespace AMDiS
 {
-
   /// Abstract base class of all operator terms
-  template <class MeshView>
+  template <class Element>
   class OperatorTerm
   {
   protected:
-    using Codim0  = typename MeshView::template Codim<0>;
-    using Element = typename Codim0::Entity;
-
-    static constexpr int dim = MeshView::dimension;
-    static constexpr int dow = MeshView::dimensionworld;
+    static constexpr int dim = Element::dimension;
+    static constexpr int dow = Element::Geometry::coorddimension;
 
-    using QuadratureRule = Dune::QuadratureRule<double, dim>;
-    using PointList = std::vector<Dune::QuadraturePoint<double, dim>>;
+    using QuadratureRule = QuadratureRuleFactory_t<Element,double,dim>;
+    // using PointList = std::vector<Dune::QuadraturePoint<double, dim>>;
 
   public:
-    virtual void init(Element const& element, PointList const& points) = 0;
+    virtual void init(Element const& element, QuadratureRule const& points) = 0;
 
     virtual double evalZot(std::size_t iq,
         Dune::FieldVector<double,1> const& test,
@@ -54,17 +51,16 @@ namespace AMDiS
 
 
   /// Base class for all operators based on expressions
-  template <class MeshView, class Term, class Traits = tag::none>
+  template <class Element, class Term, class Traits = tag::none>
   class GenericOperatorTerm
-      : public OperatorTerm<MeshView>
+      : public OperatorTerm<Element>
       , public OperatorEvaluation // NOTE: better: use static functions
   {
-    using Super   = OperatorTerm<MeshView>;
-    using Element = typename Super::Element;
-    using PointList = typename Super::PointList;
+    using Super   = OperatorTerm<Element>;
+    using QuadratureRule = typename Super::QuadratureRule;
 
-    static constexpr int dim = MeshView::dimension;
-    static constexpr int dow = MeshView::dimensionworld;
+    static constexpr int dim = Super::dim;
+    static constexpr int dow = Super::dow;
 
   public:
     GenericOperatorTerm(Term const& term, Traits traits = {})
@@ -72,9 +68,9 @@ namespace AMDiS
       , traits(traits)
     {}
 
-    virtual void init(Element const& element, PointList const& points) override
+    virtual void init(Element const& element, QuadratureRule const& points) override
     {
-      term.init(element, points);
+      term.init(get_entity(element), points);
 
       // cache term evaluation
       values.resize(points.size());

dune/amdis/ProblemStat.hpp

@@ -73,7 +73,9 @@ namespace AMDiS
     using SystemVectorType = SystemVector<FeSpaces>;
     using SystemMatrixType = SystemMatrix<FeSpaces>;
 
-    using OperatorType = Operator<MeshView>;
+    using OperatorType = Operator<MeshView, Element>;
+    using BoundaryOperatorType = Operator<MeshView, typename MeshView::Intersection>;
+
     using LinearSolverType = LinearSolverInterface<typename SystemMatrixType::MultiMatrix,