Merge branch 'feature/constraints' into 'develop'

Boundary conditions and Constraints

See merge request spraetor/dune-amdis!71

dune.module

@@ -6,4 +6,4 @@ Module: amdis
 Version: 0.1
 Maintainer: simon.praetorius@tu-dresden.de
 Depends: dune-common (>= 2.6) dune-geometry (>= 2.6) dune-localfunctions (>= 2.6) dune-typetree (>= 2.6) dune-grid (>= 2.6) dune-functions (>= 2.6)
-Suggests: dune-uggrid dune-alugrid dune-foamgrid
+Suggests: dune-uggrid dune-alugrid dune-foamgrid dune-spgrid

examples/CMakeLists.txt

 add_custom_target(examples)
 
-set(projects2d "ellipt" "heat" "vecellipt" "stokes0" "stokes1" "stokes3" "navier_stokes" "convection_diffusion")
+set(projects2d
+  "ellipt"
+  "heat"
+  "vecellipt"
+  "stokes0"
+  "stokes1"
+  "stokes3"
+  "navier_stokes"
+  "convection_diffusion"
+  "boundary"
+  "periodic"
+  "neumann")
 
 foreach(project ${projects2d})
     add_executable(${project}.2d ${project}.cc)

examples/boundary.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+#include <iostream>
+
+#if HAVE_DUNE_SPGRID
+#include <dune/grid/spgrid.hh>
+#endif
+
+#include <amdis/AMDiS.hpp>
+#include <amdis/Integrate.hpp>
+#include <amdis/ProblemStat.hpp>
+#include <amdis/Operators.hpp>
+#include <amdis/common/Literals.hpp>
+
+using namespace AMDiS;
+using namespace Dune::Indices;
+
+// 1 component with polynomial degree 2
+using Param   = YaspGridBasis<AMDIS_DIM, 2>;
+using ElliptProblem = ProblemStat<Param>;
+
+template <class SetBoundary>
+void run(SetBoundary setBoundary)
+{
+  ElliptProblem prob("ellipt");
+  prob.initialize(INIT_ALL);
+  setBoundary(prob.boundaryManager());
+
+  auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
+  prob.addMatrixOperator(opL, 0, 0);
+
+  auto f = [](auto const& x) {
+    double r2 = dot(x,x);
+    double ux = std::exp(-10.0 * r2);
+    return -(400.0 * r2 - 20.0 * x.size()) * ux;
+  };
+  auto opForce = makeOperator(tag::test{}, f, 6);
+  prob.addVectorOperator(opForce, 0);
+
+  // set boundary condition
+  auto g = [](auto const& x){ return std::exp(-10.0 * dot(x,x)); };
+  prob.addDirichletBC(BoundaryType{1}, 0, 0, g);
+
+  AdaptInfo adaptInfo("adapt");
+
+  prob.assemble(adaptInfo);
+  prob.solve(adaptInfo);
+
+  double errorL2 = integrate(sqr(g - prob.solution(0)), prob.gridView(), 6);
+  msg("error_L2 = {}", errorL2);
+}
+
+void run_periodic()
+{
+#if HAVE_DUNE_SPGRID
+  Dune::SPCube<double,2> cube({0.0,0.0},{1.0,1.0});
+  Dune::SPDomain<double,2> domain({cube}, Dune::SPTopology<2>(0b01));
+  Dune::SPGrid<double,2> grid(domain, Dune::SPMultiIndex<2>({2,2}));
+  using Grid = Dune::SPGrid<double,2>;
+#else
+  Dune::YaspGrid<2> grid({1.0,1.0},{2,2},std::bitset<2>("10"),0);
+  using Grid = Dune::YaspGrid<2>;
+#endif
+
+  using Traits = LagrangeBasis<typename Grid::LeafGridView, 2>;
+  ProblemStat<Traits> prob("ellipt", grid);
+  prob.initialize(INIT_ALL);
+  prob.boundaryManager().setBoxBoundary({-1,-1,1,1});
+
+  auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
+  prob.addMatrixOperator(opL, 0, 0);
+
+  auto opForce = makeOperator(tag::test{}, 1.0, 6);
+  prob.addVectorOperator(opForce, 0);
+
+  // set boundary condition
+  auto g = [](auto const& x) {
+    return std::sin(0.5*M_PI*x[1])*std::sin(2*M_PI * (0.25 + x[0]))
+         + std::cos(0.5*M_PI*x[1])*std::sin(2*M_PI * (-0.25 + x[0]));
+  };
+  prob.addDirichletBC(BoundaryType{1}, 0, 0, g);
+
+  typename ElliptProblem::WorldMatrix A{{1.0,0.0}, {0.0,1.0}};
+  typename ElliptProblem::WorldVector b{1.0, 0.0};
+  prob.addPeriodicBC(BoundaryType{-1}, A, b);
+
+  AdaptInfo adaptInfo("adapt");
+
+  prob.assemble(adaptInfo);
+  prob.solve(adaptInfo);
+  prob.writeFiles(adaptInfo, true);
+}
+
+
+int main(int argc, char** argv)
+{
+  AMDiS::init(argc, argv);
+
+  auto b = [](auto const& x){ return x[0] < 1.e-8 || x[1] < 1.e-8 || x[0] > 1.0-1.e-8 || x[1] > 1.0-1.e-8; };
+
+  auto setBoundary1 = [](auto& boundaryManager)
+  {
+    boundaryManager.setBoxBoundary({1,1,1,1});
+  };
+  run(setBoundary1);
+
+  auto setBoundary2 = [b](auto& boundaryManager)
+  {
+    boundaryManager.setIndicator([b](auto const& x) { return b(x) ? 1 : 0; });
+  };
+  run(setBoundary2);
+
+  auto setBoundary3 = [b](auto& boundaryManager)
+  {
+    boundaryManager.setPredicate(b, 1);
+  };
+  run(setBoundary3);
+
+  run_periodic();
+
+  AMDiS::finalize();
+  return 0;
+}

examples/ellipt.cc

@@ -44,8 +44,7 @@ int main(int argc, char** argv)
   prob.addVectorOperator(opForce, 0);
 
   // set boundary condition
-  auto boundary = [](auto const& x){ return x[0] < 1.e-8 || x[1] < 1.e-8 || x[0] > 1.0-1.e-8 || x[1] > 1.0-1.e-8; };
-  prob.addDirichletBC(boundary, 0, 0, g);
+  prob.addDirichletBC(BoundaryType{1}, 0, 0, g);
 
   AdaptInfo adaptInfo("adapt");
 

examples/init/boundary.dat.2d

+elliptMesh->global refinements: 3
+
+ellipt->mesh:                   elliptMesh
+
+ellipt->solver->name:           bicgstab
+ellipt->solver->max iteration:  10000
+ellipt->solver->relative tolerance: 1.e-10
+ellipt->solver->info:           10
+ellipt->solver->left precon:    ilu
+ellipt->solver->right precon:   no
+ellipt->solver->precon->name:   ilu
+
+ellipt->output[0]->filename:    boundary.2d
+ellipt->output[0]->name:        u
+ellipt->output[0]->output directory: ./output

examples/neumann.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+#include <iostream>
+
+#if HAVE_DUNE_SPGRID
+#include <dune/grid/spgrid.hh>
+#endif
+#if HAVE_DUNE_ALUGRID
+#include <dune/alugrid/grid.hh>
+#endif
+#include <dune/grid/uggrid.hh>
+
+#include <amdis/AMDiS.hpp>
+#include <amdis/Integrate.hpp>
+#include <amdis/ProblemStat.hpp>
+#include <amdis/Operators.hpp>
+#include <amdis/common/Literals.hpp>
+
+using namespace AMDiS;
+using namespace Dune::Indices;
+
+template <class Grid>
+void run(Grid& grid)
+{
+  grid.globalRefine(Grid::dimension == 2 ? 4 : 2);
+
+  using Traits = LagrangeBasis<typename Grid::LeafGridView, 2>;
+  ProblemStat<Traits> prob("ellipt", grid);
+  prob.initialize(INIT_ALL);
+  prob.boundaryManager().setBoxBoundary({1,2,2,1});
+
+  auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
+  prob.addMatrixOperator(opL, 0, 0);
+
+  // set dirichlet boundary condition
+  auto g = [](auto const& x) { return 0.0; };
+  prob.addDirichletBC(BoundaryType{1}, 0, 0, g);
+
+  // set neumann boundary condition
+  auto opNeumann = makeOperator(tag::test{}, 1.0);
+  prob.addVectorOperator(BoundaryType{2}, opNeumann, 0);
+
+  AdaptInfo adaptInfo("adapt");
+
+  prob.assemble(adaptInfo);
+  prob.solve(adaptInfo);
+  prob.writeFiles(adaptInfo, true);
+}
+
+
+int main(int argc, char** argv)
+{
+  AMDiS::init(argc, argv);
+
+  // 2d grids
+
+  Dune::YaspGrid<2> grid0({1.0,1.0},{2,2});
+  run(grid0);
+
+#if HAVE_DUNE_SPGRID
+  Dune::SPDomain<double,2> domain({0.0,0.0}, {1.0,1.0});
+  Dune::SPGrid<double,2> grid1(domain, Dune::SPMultiIndex<2>({2,2}));
+  run(grid1);
+#endif
+
+#if HAVE_DUNE_ALUGRID
+  using Grid2 = Dune::ALUGrid<2,2,Dune::simplex,Dune::conforming>;
+  using Factory2 = Dune::StructuredGridFactory<Grid2>;
+  auto grid2 = Factory2::createSimplexGrid({0.0,0.0}, {1.0,1.0},
+                                            std::array<unsigned int,2>{2u,2u});
+  run(*grid2);
+#endif
+
+#if HAVE_DUNE_UGGRID
+  using Grid3 = Dune::UGGrid<2>;
+  using Factory3 = Dune::StructuredGridFactory<Grid3>;
+  auto grid3 = Factory3::createSimplexGrid({0.0,0.0}, {1.0,1.0},
+                                            std::array<unsigned int,2>{2u,2u});
+  run(*grid3);
+#endif
+
+  // 3d grids
+
+  Dune::YaspGrid<3> grid4({1.0,1.0,1.0},{2,2,2});
+  run(grid4);
+
+  AMDiS::finalize();
+  return 0;
+}

examples/periodic.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+#include <iostream>
+
+#if HAVE_DUNE_SPGRID
+#include <dune/grid/spgrid.hh>
+#endif
+
+#include <amdis/AMDiS.hpp>
+#include <amdis/Integrate.hpp>
+#include <amdis/ProblemStat.hpp>
+#include <amdis/Operators.hpp>
+#include <amdis/common/Literals.hpp>
+
+using namespace AMDiS;
+
+
+template <class Grid>
+void print(Grid const& grid)
+{
+  auto const& indexSet = grid.leafIndexSet();
+  std::cout << "vertices:\n";
+  for (auto const& v : vertices(grid.leafGridView())) {
+    auto coord = v.geometry().corner(0);
+    std::cout << indexSet.index(v) << ": [" << coord[0] << ", " << coord[1] << "]\n";
+  }
+  std::cout << "\n";
+
+  std::cout << "elements:\n";
+  for (auto const& e : elements(grid.leafGridView())) {
+    std::cout << indexSet.index(e) << ": {" << indexSet.subIndex(e,0,2);
+    for (unsigned int i = 1; i < e.subEntities(2); ++i) {
+      std::cout << ", " << indexSet.subIndex(e,i,2);
+     }
+     std::cout << "}\n";
+  }
+  std::cout << "\n";
+
+  std::cout << "boundary-intersections:\n";
+  for (auto const& e : elements(grid.leafGridView())) {
+    for (auto const& inter : intersections(grid.leafGridView(), e)) {
+      if (!inter.boundary())
+        continue;
+
+      std::cout << inter.boundarySegmentIndex() << ": [" << inter.geometry().center() << "] {" << indexSet.index(inter.inside());
+      if (inter.neighbor())
+        std::cout << ", " << indexSet.index(inter.outside());
+      std::cout << "}\n";
+    }
+  }
+}
+
+template <class Grid>
+void run(Grid& grid)
+{
+  using Traits = LagrangeBasis<typename Grid::LeafGridView, 1>;
+
+  ProblemStat<Traits> prob("ellipt", grid);
+  prob.initialize(INIT_ALL);
+  prob.boundaryManager().setBoxBoundary({-1,-1,1,1});
+
+  print(grid);
+
+  using BC = PeriodicBC<FieldVector<double,2>, typename Traits::GlobalBasis::MultiIndex>;
+  BC periodicBC(prob.boundaryManagerPtr(),-1,{{{1.0,0.0}, {0.0,1.0}}, {1.0, 0.0}});
+
+  periodicBC.init(prob.globalBasis(), prob.globalBasis());
+  std::cout << "periodicNodes:\n";
+  for (std::size_t i = 0; i < periodicBC.periodicNodes().size(); ++i)
+    std::cout << i << ": " << periodicBC.periodicNodes()[i] << "\n";
+  std::cout << "\n";
+
+  std::cout << "associations:\n";
+  for (auto const& a : periodicBC.associations())
+    std::cout << a.first << " -> " << a.second << "\n";
+  std::cout << "\n";
+}
+
+int main(int argc, char** argv)
+{
+  AMDiS::init(argc, argv);
+
+#if HAVE_DUNE_SPGRID
+  Dune::SPCube<double,2> cube({0.0,0.0},{1.0,1.0});
+  Dune::SPDomain<double,2> domain({cube}, Dune::SPTopology<2>(0b01));
+  Dune::SPGrid<double,2> grid1(domain, Dune::SPMultiIndex<2>({2,2}));
+  run(grid1);
+#endif
+
+  // NOTE: 'periodic' YaspGrid not supported yet, but a simple YaspGrid can be made periodic
+  Dune::YaspGrid<2> grid2({1.0,1.0}, {2,2});
+  run(grid2);
+
+  AMDiS::finalize();
+  return 0;
+}

src/amdis/Assembler.hpp

@@ -2,15 +2,30 @@
 
 namespace AMDiS
 {
-  template <class GridView, class Element, class Operators, class ElementAssembler>
+  template <class ElementAssembler, class IntersectionAssembler, class BoundaryAssembler>
+  struct AssemblerTriple
+  {
+    ElementAssembler elementAssembler;
+    IntersectionAssembler intersectionAssembler;
+    BoundaryAssembler boundaryAssembler;
+  };
+
+  template <class EA, class IA, class BA>
+  AssemblerTriple<EA, IA, BA> makeAssemblerTriple(EA const& ea, IA const& ia, BA const& ba)
+  {
+    return {ea, ia, ba};
+  }
+
+
+  template <class GridView, class Element, class Operators, class EA, class IA, class BA>
   void assembleOperators(
       GridView const& gridView,
       Element const& element,
       Operators& operators,
-      ElementAssembler const& localAssembler)
+      AssemblerTriple<EA,IA,BA> const& assemblerTriple)
   {
     // assemble element operators
-    localAssembler(element, operators.onElement());
+    assemblerTriple.elementAssembler(element, operators.onElement());
 
     // assemble intersection operators
     if (!operators.onIntersection().empty()
@@ -18,11 +33,53 @@ namespace AMDiS
     {
       for (auto const& intersection : intersections(gridView, element)) {
         if (intersection.boundary())
-          localAssembler(intersection, operators.onBoundary());
+          assemblerTriple.boundaryAssembler(intersection, operators.onBoundary());
         else
-          localAssembler(intersection, operators.onIntersection());
+          assemblerTriple.intersectionAssembler(intersection, operators.onIntersection());
+      }
+    }
   }
+
+
+  template <class Node, class ElementVector>
+  auto makeVectorAssembler(Node const& node, ElementVector& elementVector)
+  {
+    return makeAssemblerTriple(
+      [&](auto const& element, auto& operators) {
+        for (auto op : operators)
+          op->assemble(element, node, elementVector);
+      },
+      [&](auto const& is, auto& operators) {
+        for (auto op : operators)
+          op->assemble(is, node, elementVector);
+      },
+      [&](auto const& bis, auto& operators) {
+        for (auto data : operators) {
+          if (data.bc.onBoundary(bis))
+            data.op->assemble(bis, node, elementVector);
+        }
+      });
+  }
+
+
+  template <class RowNode, class ColNode, class ElementMatrix>
+  auto makeMatrixAssembler(RowNode const& rowNode, ColNode const& colNode, ElementMatrix& elementMatrix)
+  {
+    return makeAssemblerTriple(
+      [&](auto const& element, auto& operators) {
+        for (auto op : operators)
+          op->assemble(element, rowNode, colNode, elementMatrix);
+      },
+      [&](auto const& is, auto& operators) {
+        for (auto op : operators)
+          op->assemble(is, rowNode, colNode, elementMatrix);
+      },
+      [&](auto const& bis, auto& operators) {
+        for (auto data : operators) {
+          if (data.bc.onBoundary(bis))
+            data.op->assemble(bis, rowNode, colNode, elementMatrix);
         }
+      });
   }
 
 } // end namespace AMDiS

src/amdis/Boundary.hpp

@@ -2,6 +2,6 @@
 
 namespace AMDiS
 {
-  struct BoundaryType { int b; };
+  using BoundaryType = int;
 
 } // end namespace AMDiS

src/amdis/BoundaryCondition.hpp

+#pragma once
+
+#include <memory>
+
+#include <amdis/Boundary.hpp>
+#include <amdis/BoundaryManager.hpp>
+
+namespace AMDiS
+{
+  class BoundaryCondition
+  {
+  public:
+    BoundaryCondition() = default;
+    BoundaryCondition(std::shared_ptr<BoundaryManagerBase const> const& boundaryManager, BoundaryType id)
+      : boundaryManager_(boundaryManager)
+      , id_(id)
+    {}
+
+    /// Return true if intersection is on boundary with id
+    template <class Intersection>
+    bool onBoundary(Intersection const& is) const
+    {
+      return is.boundary() && (!boundaryManager_  || boundaryManager_->boundaryId(is) == id_);
+    }
+
+    template <class RowBasis, class ColBasis>
+    void init(RowBasis const& rowBasis, ColBasis const& colBasis) { /* do nothing */ }
+
+    template <class Matrix, class X, class B, class RN, class CN>
+    void fillBoundaryCondition(Matrix& A, X& x, B& b, RN const& rowNode, CN const& colNode) { /* do nothing */ }
+
+  protected:
+    std::shared_ptr<BoundaryManagerBase const> boundaryManager_{nullptr};
+    BoundaryType id_{0};
+  };
+
+} // end namespace AMDiS

src/amdis/BoundaryManager.hpp

+#pragma once
+
+#include <memory>
+#include <vector>
+
+#include <dune/common/hybridutilities.hh>
+#include <dune/common/std/type_traits.hh>
+
+#include <amdis/Boundary.hpp>
+#include <amdis/common/Concepts.hpp>
+
+namespace AMDiS
+{
+  class BoundaryManagerBase
+  {
+  public:
+    BoundaryManagerBase(std::size_t numBoundarySegments)
+      : boundaryIds_(numBoundarySegments, BoundaryType{1})
+    {}
+
+    /// Return the stored boundary id for the given intersection
+    template <class Intersection>
+    BoundaryType boundaryId(Intersection const& intersection) const
+    {
+      return boundaryIds_[intersection.boundarySegmentIndex()];
+    }
+
+  protected:
+    std::vector<BoundaryType> boundaryIds_; // maps a boundarySegementIndex to an ID
+  };
+
+
+  /// Manage boundary ids of boundary segments in a grid
+  /**
+   * Manager for bounary IDs, that can be initialized from different sources
+   * - cube domains can be assigned box boundaries, by specifying left-right,
+   *   front-back, and bottom-top ids
+   * - An indicator from global coodinates that returns an id. The indicator function
+   *   is evaluated in the center points of intersections.
+   * - A predicate that return true/false for a boundary part and sets the given id
+   * - Read ids from the grid. Those may be initialized by some grid readers. Note:
+   *   not all grids support `intersection.boundaryId()`, but e.g. AlbertaGrid and
+   *   ALUGrid do support this. Can be read by DGF parser and AlbertaGrid constructor
+   *   from macroFileName.
+   **/
+  template <class G>
+  class BoundaryManager
+      : public BoundaryManagerBase
+  {
+    using Super = BoundaryManagerBase;
+
+  public:
+    using Grid = G;
+    enum { dim = Grid::dimension };
+    enum { dow = Grid::dimensionworld };
+
+    using Segment = typename Grid::LevelGridView::Intersection;
+    using Domain = typename Grid::template Codim<0>::Geometry::GlobalCoordinate;
+
+  public:
+    /// Constructor. Stores a shared pointer to the grid and initializes all
+    /// boundary IDs to the value stored in the grid
+    BoundaryManager(std::shared_ptr<G> const& grid)
+      : Super(grid->numBoundarySegments())
+      , grid_(grid)
+    {
+      setBoundaryId();
+    }
+
+    /// Set boundary ids [left,right, front,back, bottom,top] for cube domains
+    void setBoxBoundary(std::array<BoundaryType, 2*dow> const& ids)
+    {
+      auto gv = grid_->levelGridView(0);
+      for (auto const& e : elements(gv))
+      {
+        for (auto const& segment : intersections(gv,e)) {
+          if (!segment.boundary())
+            continue;
+
+          auto n = segment.centerUnitOuterNormal();
+          auto index = segment.boundarySegmentIndex();
+
+          for (int i = 0; i < dow; ++i) {
+            if (n[i] < -0.5)
+              boundaryIds_[index] = ids[2*i];
+            else if (n[i] > 0.5)
+              boundaryIds_[index] = ids[2*i+1];
+          }
+        }
+      }
+    }
+
+
+    /// Set indicator(center) for all boundary intersections
+    template <class Indicator,
+      REQUIRES(Concepts::Functor<Indicator, int(Domain)>) >
+    void setIndicator(Indicator const& indicator)
+    {
+      auto gv = grid_->levelGridView(0);
+      for (auto const& e : elements(gv))
+      {
+        for (auto const& segment : intersections(gv,e)) {