diff --git a/dune.module b/dune.module
index 390ad16e23f809e35891886e8edf4c13e2b551ce..bd22ba1dd5e77e68e84713eeb8d7624f062f645a 100644
--- a/dune.module
+++ b/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
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 72c586a1570c6253d0392577eacafc381ad7a70c..14ca4fe498081040f84f9bc1fdffb6deeac8802a 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -1,6 +1,17 @@
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)
diff --git a/examples/boundary.cc b/examples/boundary.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c6f31b4b990a2e58321357bed1ebafedfb8e90f2
--- /dev/null
+++ b/examples/boundary.cc
@@ -0,0 +1,124 @@
+#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;
+}
diff --git a/examples/ellipt.cc b/examples/ellipt.cc
index 61e679262b84473cd3f5e1cfb0f1fcf19e3bec7f..3a19575a81544a55a9df0bf0a81d7baa6ea1a43b 100644
--- a/examples/ellipt.cc
+++ b/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");
diff --git a/examples/init/boundary.dat.2d b/examples/init/boundary.dat.2d
new file mode 100644
index 0000000000000000000000000000000000000000..1f284109be8700240d6856c65a04874f86fe2cf1
--- /dev/null
+++ b/examples/init/boundary.dat.2d
@@ -0,0 +1,15 @@
+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
diff --git a/examples/neumann.cc b/examples/neumann.cc
new file mode 100644
index 0000000000000000000000000000000000000000..65ab4016564cf489f33d0cdedc269b6049737c6f
--- /dev/null
+++ b/examples/neumann.cc
@@ -0,0 +1,90 @@
+#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;
+}
diff --git a/examples/periodic.cc b/examples/periodic.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f03e0a07637f8910c289b4c9f5f09883d6ff04cd
--- /dev/null
+++ b/examples/periodic.cc
@@ -0,0 +1,97 @@
+#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;
+}
diff --git a/src/amdis/Assembler.hpp b/src/amdis/Assembler.hpp
index 7452ea985a1cf682ab9c008f25d4eda85da22d7a..9da5ea595c4bade6807d37b712def39adc0c0836 100644
--- a/src/amdis/Assembler.hpp
+++ b/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
diff --git a/src/amdis/Boundary.hpp b/src/amdis/Boundary.hpp
index ee11eb4cb7781d425b87e6975675072752641037..881b5010cfabaed06810370387db14d18b855839 100644
--- a/src/amdis/Boundary.hpp
+++ b/src/amdis/Boundary.hpp
@@ -2,6 +2,6 @@
namespace AMDiS
{
- struct BoundaryType { int b; };
+ using BoundaryType = int;
} // end namespace AMDiS
diff --git a/src/amdis/BoundaryCondition.hpp b/src/amdis/BoundaryCondition.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e3a231aad23b6a1244480d7309019133baf73d3d
--- /dev/null
+++ b/src/amdis/BoundaryCondition.hpp
@@ -0,0 +1,37 @@
+#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
diff --git a/src/amdis/BoundaryManager.hpp b/src/amdis/BoundaryManager.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..375c334ae1cba89044f78358ca6b04560573a856
--- /dev/null
+++ b/src/amdis/BoundaryManager.hpp
@@ -0,0 +1,161 @@
+#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)) {
+ if (!segment.boundary())
+ continue;
+
+ auto index = segment.boundarySegmentIndex();
+ boundaryIds_[index] = indicator(segment.geometry().center());
+ }
+ }
+ }
+
+
+ /// Set id for all boundary intersections with pred(center) == true
+ template <class Predicate,
+ REQUIRES(Concepts::Functor<Predicate, bool(Domain)>) >
+ void setPredicate(Predicate const& pred, BoundaryType id)
+ {
+ auto gv = grid_->levelGridView(0);
+ for (auto const& e : elements(gv))
+ {
+ for (auto const& segment : intersections(gv,e)) {
+ if (!segment.boundary())
+ continue;
+
+ auto index = segment.boundarySegmentIndex();
+ if (pred(segment.geometry().center()))
+ boundaryIds_[index] = id;
+ }
+ }
+ }
+
+
+ template <class I>
+ using HasBoundaryId = decltype(std::declval<I>().boundaryId());
+
+ /// Set boundary ids as stored in the grid, e.g. read by grid reader
+ void setBoundaryId()
+ {
+ if (!Dune::Std::is_detected<HasBoundaryId, Segment>::value)
+ return;
+
+ auto gv = grid_->levelGridView(0);
+ for (auto const& e : elements(gv))
+ {
+ for (auto const& segment : intersections(gv,e)) {
+ if (!segment.boundary())
+ continue;
+
+ auto index = segment.boundarySegmentIndex();
+ Dune::Hybrid::ifElse(Dune::Std::is_detected<HasBoundaryId, Segment>{},
+ [&](auto id) { boundaryIds_[index] = id(segment).boundaryId(); });
+ }
+ }
+ }
+
+ private:
+ std::shared_ptr<Grid> grid_;
+ using Super::boundaryIds_;
+ };
+
+} // end namespace AMDiS
diff --git a/src/amdis/DirichletBC.hpp b/src/amdis/DirichletBC.hpp
index 82a12b18079d030f929963ac5dde39e2ef97a5ad..0d570e28f4852d9b6cba5b8e710223b10b86937f 100644
--- a/src/amdis/DirichletBC.hpp
+++ b/src/amdis/DirichletBC.hpp
@@ -2,18 +2,14 @@
#include <functional>
#include <list>
-#include <type_traits>
+#include <memory>
#include <vector>
-#include <dune/common/hybridutilities.hh>
-#include <dune/functions/functionspacebases/boundarydofs.hh>
-#include <dune/functions/functionspacebases/interpolate.hh>
-#include <dune/functions/functionspacebases/subspacebasis.hh>
+#include <dune/common/std/optional.hh>
#include <amdis/Boundary.hpp>
-#include <amdis/Output.hpp>
+#include <amdis/BoundaryCondition.hpp>
#include <amdis/common/Concepts.hpp>
-#include <amdis/common/ValueCategory.hpp>
#include <amdis/utility/RangeType.hpp>
#include <amdis/utility/TreeData.hpp>
#include <amdis/utility/Visitor.hpp>
@@ -37,82 +33,47 @@ namespace AMDiS
**/
template <class Domain, class Range = double>
class DirichletBC
+ : public BoundaryCondition
{
+ using Super = BoundaryCondition;
+
public:
+ template <class BM, class Values,
+ REQUIRES(Concepts::Functor<Values, Range(Domain)>) >
+ DirichletBC(BM&& boundaryManager, BoundaryType id, Values&& values)
+ : Super(std::forward<BM>(boundaryManager), id)
+ , values_(std::forward<Values>(values))
+ {}
+
template <class Predicate, class Values,
- REQUIRES(Concepts::Functor<Predicate, bool(Domain)> &&
- Concepts::Functor<Values, Range(Domain)>) >
+ REQUIRES(Concepts::Functor<Predicate, bool(Domain)>),
+ REQUIRES(Concepts::Functor<Values, Range(Domain)>)>
DirichletBC(Predicate&& predicate, Values&& values)
: predicate_(std::forward<Predicate>(predicate))
, values_(std::forward<Values>(values))
{}
-
- /// Prepare the matrix, rhs and solution block for assembling dirichlet
- /// boundary conditions, e.g. collect the corresponding indices of boundary
- /// DOFS for later elimination in \ref finish.
- template <class Matrix, class VectorX, class VectorB, class RowNode, class ColNode>
- void init(Matrix& /*matrix*/, VectorX& /*solution*/, VectorB& rhs,
- RowNode const& rowNode, ColNode const&)
+ template <class Intersection>
+ bool onBoundary(Intersection const& intersection) const
{
- if (!initialized_) {
- forEachLeafNode_(rowNode, [&](auto const& /*node*/, auto const& tp) {
- auto subBasis = Dune::Functions::subspaceBasis(rhs.basis(), cat(rowNode.treePath(),tp));
- this->initImpl(subBasis);
- });
- initialized_ = true;
- }
+ return predicate_ ? (*predicate_)(intersection.geometry().center())
+ : Super::onBoundary(intersection);
}
+ // fill \ref dirichletNodes_ with 1 or 0 whether DOF corresponds to the boundary or not.
+ template <class RowBasis, class ColBasis>
+ void init(RowBasis const& rowBasis, ColBasis const& colBasis);
/// \brief Apply dirichlet BC to matrix and vector, i.e., add a unit-row
/// to the matrix and optionally delete the corresponding matrix-column.
template <class Matrix, class VectorX, class VectorB, class RowNode, class ColNode>
- void finish(Matrix& matrix, VectorX& solution, VectorB& rhs,
- RowNode const& rowNode, ColNode const& colNode)
- {
- test_exit_dbg(initialized_, "Boundary condition not initialized!");
- auto columns = matrix.applyDirichletBC(dirichletNodes_);
-
- Dune::Hybrid::ifElse(std::is_same<RangeType_t<RowNode>, Range>{},
- [&](auto id) {
- auto subBasis = Dune::Functions::subspaceBasis(rhs.basis(), rowNode.treePath());
- Dune::Functions::interpolate(subBasis, rhs, values_, dirichletNodes_);
- });
-
- Dune::Hybrid::ifElse(std::is_same<RangeType_t<ColNode>, Range>{},
- [&](auto id) {
- auto subBasis = Dune::Functions::subspaceBasis(solution.basis(), colNode.treePath());
- Dune::Functions::interpolate(subBasis, solution, values_, dirichletNodes_);
- });
-
- // subtract columns of dirichlet nodes from rhs
- for (auto const& triplet : columns)
- rhs[triplet.row] -= triplet.value * solution[triplet.col];
-
- initialized_ = false;
- }
-
- protected:
- // fill \ref dirichletNodes_ with 1 or 0 whether DOF corresponds to the boundary or not.
- template <class SubBasis>
- void initImpl(SubBasis const& subBasis)
- {
- using Dune::Functions::forEachBoundaryDOF;
- using LV = typename SubBasis::LocalView;
- using I = typename SubBasis::GridView::Intersection;
- dirichletNodes_.resize(subBasis.dimension());
- forEachBoundaryDOF(subBasis, [&](int localIndex, LV const& localView, I const& intersection) {
- dirichletNodes_[localView.index(localIndex)] = predicate_(intersection.geometry().center());
- });
- }
+ void fillBoundaryCondition(Matrix& matrix, VectorX& solution, VectorB& rhs,
+ RowNode const& rowNode, ColNode const& colNode);
private:
- std::function<bool(Domain)> predicate_;
+ Dune::Std::optional<std::function<bool(Domain)>> predicate_;
std::function<Range(Domain)> values_;
-
- bool initialized_ = false;
- std::vector<char> dirichletNodes_;
+ std::vector<bool> dirichletNodes_;
};
@@ -126,6 +87,8 @@ namespace AMDiS
};
template <class RowBasis, class ColBasis>
- using Constraints = MatrixData<RowBasis, ColBasis, DirichletData<RowBasis>::template type>;
+ using DirichletBCs = MatrixData<RowBasis, ColBasis, DirichletData<RowBasis>::template type>;
} // end namespace AMDiS
+
+#include "DirichletBC.inc.hpp"
diff --git a/src/amdis/DirichletBC.inc.hpp b/src/amdis/DirichletBC.inc.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d3ab7385eacf67fa0676d8607cd02c5c0fad0efe
--- /dev/null
+++ b/src/amdis/DirichletBC.inc.hpp
@@ -0,0 +1,55 @@
+#pragma once
+
+#include <type_traits>
+
+#include <dune/common/hybridutilities.hh>
+#include <dune/functions/functionspacebases/boundarydofs.hh>
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/functionspacebases/subspacebasis.hh>
+
+#include <amdis/Output.hpp>
+#include <amdis/linear_algebra/Constraints.hpp>
+
+namespace AMDiS {
+
+template <class D, class R>
+ template <class RowBasis, class ColBasis>
+void DirichletBC<D,R>::
+init(RowBasis const& rowBasis, ColBasis const& colBasis)
+{
+ using Dune::Functions::forEachBoundaryDOF;
+ using LV = typename ColBasis::LocalView;
+ using I = typename ColBasis::GridView::Intersection;
+ dirichletNodes_.resize(colBasis.dimension());
+ forEachBoundaryDOF(colBasis, [&](int localIndex, LV const& localView, I const& intersection) {
+ dirichletNodes_[localView.index(localIndex)] = onBoundary(intersection);
+ });
+}
+
+
+template <class D, class R>
+ template <class Matrix, class VectorX, class VectorB, class RowNode, class ColNode>
+void DirichletBC<D,R>::
+fillBoundaryCondition(Matrix& matrix, VectorX& solution, VectorB& rhs,
+ RowNode const& rowNode, ColNode const& colNode)
+{
+ auto columns = Constraints<Matrix>::dirichletBC(matrix, dirichletNodes_);
+
+ Dune::Hybrid::ifElse(std::is_same<RangeType_t<RowNode>, R>{},
+ [&](auto id) {
+ auto subBasis = Dune::Functions::subspaceBasis(rhs.basis(), rowNode.treePath());
+ Dune::Functions::interpolate(subBasis, rhs, values_, dirichletNodes_);
+ });
+
+ // subtract columns of dirichlet nodes from rhs
+ for (auto const& triplet : columns)
+ rhs[triplet.row] -= triplet.value * solution[triplet.col];
+
+ Dune::Hybrid::ifElse(std::is_same<RangeType_t<ColNode>, R>{},
+ [&](auto id) {
+ auto subBasis = Dune::Functions::subspaceBasis(solution.basis(), colNode.treePath());
+ Dune::Functions::interpolate(subBasis, solution, values_, dirichletNodes_);
+ });
+}
+
+} // end namespace AMDiS
diff --git a/src/amdis/LinearAlgebra.hpp b/src/amdis/LinearAlgebra.hpp
index 6184b2042b8dc444e9f368bc56262c8f9fb83f2b..aa920278aeb6156d4224030cfc4f745352e4cea3 100644
--- a/src/amdis/LinearAlgebra.hpp
+++ b/src/amdis/LinearAlgebra.hpp
@@ -2,16 +2,19 @@
#if HAVE_MTL
+#include <amdis/linear_algebra/mtl/Constraints.hpp>
#include <amdis/linear_algebra/mtl/DOFMatrix.hpp>
#include <amdis/linear_algebra/mtl/DOFVector.hpp>
#elif HAVE_EIGEN
+#include <amdis/linear_algebra/eigen/Constraints.hpp>
#include <amdis/linear_algebra/eigen/DOFMatrix.hpp>
#include <amdis/linear_algebra/eigen/DOFVector.hpp>
#else // ISTL
+#include <amdis/linear_algebra/istl/Constraints.hpp>
#include <amdis/linear_algebra/istl/DOFMatrix.hpp>
#include <amdis/linear_algebra/istl/DOFVector.hpp>
diff --git a/src/amdis/LocalAssemblerList.hpp b/src/amdis/LocalAssemblerList.hpp
index b9e7fc694542bcdff421069de9050eedbd222729..a5e993d63d6f5e384b1f732555345c96909ae7b0 100644
--- a/src/amdis/LocalAssemblerList.hpp
+++ b/src/amdis/LocalAssemblerList.hpp
@@ -3,7 +3,7 @@
#include <list>
#include <memory>
-#include <amdis/Boundary.hpp>
+#include <amdis/BoundaryCondition.hpp>
#include <amdis/LocalAssemblerBase.hpp>
#include <amdis/utility/TreeData.hpp>
@@ -13,7 +13,11 @@ namespace AMDiS
{
template <class E> struct element_operator { using type = E; };
template <class I> struct intersection_operator { using type = I; };
- template <class I> struct boundary_operator { using type = I; BoundaryType id; };
+ template <class I> struct boundary_operator : public BoundaryCondition
+ {
+ using type = I;
+ using BoundaryCondition::BoundaryCondition;
+ };
}
template <class GridView>
@@ -26,7 +30,7 @@ namespace AMDiS
struct DataElement
{
std::shared_ptr<OperatorType> op;
- BoundaryType b = {0};
+ BoundaryCondition bc;
};
/// Lists of \ref DataElement on the Element, BoundaryIntersction, and InteriorIntersections
@@ -54,49 +58,49 @@ namespace AMDiS
template <class Geo>
void bind(Element const& elem, Geo const& geo)
{
- for (auto& scaled : element_) scaled.op->bind(elem,geo);
- for (auto& scaled : boundary_) scaled.op->bind(elem,geo);
- for (auto& scaled : intersection_) scaled.op->bind(elem,geo);
+ for (auto op : element_) op->bind(elem,geo);
+ for (auto op : intersection_) op->bind(elem,geo);
+ for (auto& data : boundary_) data.op->bind(elem,geo);
}
/// Unbind all operators from the element
void unbind()
{
- for (auto& scaled : element_) scaled.op->unbind();
- for (auto& scaled : boundary_) scaled.op->unbind();
- for (auto& scaled : intersection_) scaled.op->unbind();
+ for (auto op : element_) op->unbind();
+ for (auto op : intersection_) op->unbind();
+ for (auto& data : boundary_) data.op->unbind();
assembled_ = true;
}
- template <class... Args>
- void push(tag::element_operator<Element>, Args&&... args)
+ template <class Op>
+ void push(tag::element_operator<Element>, Op&& op)
{
- element_.push_back({std::forward<Args>(args)...});
+ element_.emplace_back(std::forward<Op>(op));
}
- template <class... Args>
- void push(tag::boundary_operator<Intersection> b, Args&&... args)
+ template <class Op>
+ void push(tag::intersection_operator<Intersection>, Op&& op)
{
- boundary_.push_back({std::forward<Args>(args)..., b.id});
+ intersection_.emplace_back(std::forward<Op>(op));
}
- template <class... Args>
- void push(tag::intersection_operator<Intersection>, Args&&... args)
+ template <class Op>
+ void push(tag::boundary_operator<Intersection> b, Op&& op)
{
- intersection_.push_back({std::forward<Args>(args)...});
+ boundary_.push_back({std::forward<Op>(op), b});
}
auto& onElement() { return element_; }
auto const& onElement() const { return element_; }
- auto& onBoundary() { return boundary_; }
- auto const& onBoundary() const { return boundary_; }
-
auto& onIntersection() { return intersection_; }
auto const& onIntersection() const { return intersection_; }
+ auto& onBoundary() { return boundary_; }
+ auto const& onBoundary() const { return boundary_; }
+
private:
/// The type of local operators associated with grid elements
using ElementOperator = LocalAssemblerBase<Element, Nodes...>;
@@ -105,14 +109,14 @@ namespace AMDiS
using IntersectionOperator = LocalAssemblerBase<Intersection, Nodes...>;
/// List of operators to be assembled on grid elements
- std::vector<DataElement<ElementOperator>> element_;
+ std::vector<std::shared_ptr<ElementOperator>> element_;
+
+ /// List of operators to be assembled on interior intersections
+ std::vector<std::shared_ptr<IntersectionOperator>> intersection_;
/// List of operators to be assembled on boundary intersections
std::vector<DataElement<IntersectionOperator>> boundary_;
- /// List of operators to be assembled on interior intersections
- std::vector<DataElement<IntersectionOperator>> intersection_;
-
/// if false, do reassemble of all operators
bool assembled_ = false;
@@ -133,9 +137,11 @@ namespace AMDiS
template <class RowBasis, class ColBasis>
- using MatrixOperators = MatrixData<RowBasis, ColBasis, OperatorLists<typename RowBasis::GridView>::template MatData>;
+ using MatrixOperators
+ = MatrixData<RowBasis, ColBasis, OperatorLists<typename RowBasis::GridView>::template MatData>;
template <class GlobalBasis>
- using VectorOperators = VectorData<GlobalBasis, OperatorLists<typename GlobalBasis::GridView>::template VecData>;
+ using VectorOperators
+ = VectorData<GlobalBasis, OperatorLists<typename GlobalBasis::GridView>::template VecData>;
} // end namespace AMDiS
diff --git a/src/amdis/LocalBasisCache.hpp b/src/amdis/LocalBasisCache.hpp
index 22c904fee9f46133c8db5e5eca8896e813087726..cab8d27c58129855977e1319f3f9af3a33ddf4bb 100644
--- a/src/amdis/LocalBasisCache.hpp
+++ b/src/amdis/LocalBasisCache.hpp
@@ -95,8 +95,8 @@ namespace AMDiS
{}
// evaluate local basis-function at all quadrature points
- template <class LocalContext>
- ShapeValuesContainer const& evaluateFunctionAtQp(LocalContext const& context, QuadratureRule const& quad) const
+ template <class LocalContext, int d>
+ ShapeValuesContainer const& evaluateFunctionAtQp(LocalContext const& context, Dune::QuadratureRule<RangeFieldType,d> const& quad) const
{
ValuesKey key{context.type().id(),quad.order(),quad.size()};
return ShapeValuesContainerCache::get(key, [&](ValuesKey const&)
@@ -109,8 +109,8 @@ namespace AMDiS
}
// evaluate local basis-gradients at all quadrature points
- template <class LocalContext>
- ShapeGradientsContainer const& evaluateJacobianAtQp(LocalContext const& context, QuadratureRule const& quad) const
+ template <class LocalContext, int d>
+ ShapeGradientsContainer const& evaluateJacobianAtQp(LocalContext const& context, Dune::QuadratureRule<RangeFieldType,d> const& quad) const
{
GradientsKey key{context.type().id(),quad.order(),quad.size()};
return ShapeGradientsContainerCache::get(key, [&](GradientsKey const&)
diff --git a/src/amdis/LocalOperator.hpp b/src/amdis/LocalOperator.hpp
index 5083e070b321bdc5555cf32753ea6229ebb93b54..ef6bfcf46b8b1b1d7d02e11b403ba7714d94bbd4 100644
--- a/src/amdis/LocalOperator.hpp
+++ b/src/amdis/LocalOperator.hpp
@@ -4,6 +4,7 @@
#include <type_traits>
#include <amdis/GridFunctions.hpp>
+#include <amdis/Output.hpp>
#include <amdis/common/Utility.hpp>
#include <amdis/utility/FiniteElementType.hpp>
diff --git a/src/amdis/PeriodicBC.hpp b/src/amdis/PeriodicBC.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5d10800a83c5343c4201ff1c1e6f672f9925ad98
--- /dev/null
+++ b/src/amdis/PeriodicBC.hpp
@@ -0,0 +1,110 @@
+#pragma once
+
+#include <functional>
+#include <list>
+#include <map>
+#include <type_traits>
+#include <vector>
+
+#include <amdis/common/ConceptsBase.hpp>
+
+#include <amdis/Output.hpp>
+
+namespace AMDiS
+{
+ template <class ct, int dow>
+ struct FaceTransformation
+ {
+ using WorldMatrix = FieldMatrix<ct, dow, dow>;
+ using WorldVector = FieldVector<ct, dow>;
+
+ WorldMatrix matrix_; // Note: We assume the matrix to be orthogonal, here
+ WorldVector shift_;
+
+ WorldVector evaluate(WorldVector const& x) const
+ {
+ WorldVector y = shift_;
+ matrix_.umv(x, y);
+ return y;
+ }
+
+ WorldVector evaluateInverse(WorldVector const& y) const
+ {
+ WorldVector ys = y - shift_;
+ WorldVector x;
+ matrix_.mtv(ys, x);
+ return x;
+ }
+ };
+
+
+ /// Implements a periodic boundary condition
+ template <class D, class MI>
+ class PeriodicBC
+ : public BoundaryCondition
+ {
+ using Super = BoundaryCondition;
+
+ public:
+ using Domain = D;
+ using MultiIndex = MI;
+ using FaceTrafo = FaceTransformation<typename D::field_type, D::dimension>;
+
+ public:
+ template <class BM>
+ PeriodicBC(BM&& boundaryManager, BoundaryType id, FaceTrafo faceTrafo)
+ : Super(std::forward<BM>(boundaryManager), id)
+ , faceTrafo_(std::move(faceTrafo))
+ {}
+
+ template <class RowBasis, class ColBasis>
+ void init(RowBasis const& rowBasis, ColBasis const& colBasis);
+
+ 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);
+
+ auto const& associations() const
+ {
+ return associations_;
+ }
+
+ auto const& periodicNodes() const
+ {
+ return periodicNodes_;
+ }
+
+ protected:
+ template <class Basis>
+ void initAssociations(Basis const& basis);
+
+ template <class Basis>
+ void initAssociations2(Basis const& basis);
+
+ // Get the coordinates of the DOFs
+ template <class Node>
+ std::vector<Domain> coords(Node const& node, std::vector<std::size_t> const& localIndices) const;
+
+ private:
+ FaceTrafo faceTrafo_;
+
+ std::vector<bool> periodicNodes_;
+ std::map<std::size_t, std::size_t> associations_;
+ Dune::Std::optional<bool> hasConnectedIntersections_;
+ };
+
+
+ template <class Basis>
+ struct PeriodicData
+ {
+ using WorldVector = typename Basis::GridView::template Codim<0>::Geometry::GlobalCoordinate;
+
+ template <class RowNode, class ColNode>
+ using type = std::list< std::shared_ptr<PeriodicBC<WorldVector, typename Basis::MultiIndex>> >;
+ };
+
+ template <class RowBasis, class ColBasis>
+ using PeriodicBCs = MatrixData<RowBasis, ColBasis, PeriodicData<RowBasis>::template type>;
+
+} // end namespace AMDiS
+
+#include "PeriodicBC.inc.hpp"
diff --git a/src/amdis/PeriodicBC.inc.hpp b/src/amdis/PeriodicBC.inc.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ea08e9ed52e0e1eb281fb5d62e6f3d573e19fbcf
--- /dev/null
+++ b/src/amdis/PeriodicBC.inc.hpp
@@ -0,0 +1,248 @@
+#pragma once
+
+#include <limits>
+
+#include <dune/common/reservedvector.hh>
+#include <dune/functions/common/functionfromcallable.hh>
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/functionspacebases/subentitydofs.hh>
+#include <amdis/LinearAlgebra.hpp>
+
+namespace AMDiS {
+
+template <class D, class MI>
+ template <class Basis, class ColBasis>
+void PeriodicBC<D,MI>::
+init(Basis const& basis, ColBasis const& /*colBasis*/)
+{
+ if (!bool(hasConnectedIntersections_)) {
+ hasConnectedIntersections_ = true;
+ const auto& gridView = basis.gridView().grid().levelGridView(0);
+ for (auto const& element : elements(gridView)) {
+ for (const auto& intersection : intersections(gridView, element)) {
+ if (!this->onBoundary(intersection))
+ continue;
+
+ if (!intersection.neighbor()) {
+ hasConnectedIntersections_ = false;
+ break;
+ }
+ }
+ }
+ }
+
+ if (*hasConnectedIntersections_)
+ initAssociations(basis);
+ else
+ initAssociations2(basis);
+}
+
+
+template <class D, class MI>
+ template <class Basis>
+void PeriodicBC<D,MI>::
+initAssociations(Basis const& basis)
+{
+ static const double tol = std::sqrt(std::numeric_limits<typename D::field_type>::epsilon());
+ periodicNodes_.clear();
+ periodicNodes_.resize(basis.dimension(), false);
+
+ auto localView = basis.localView();
+ auto seDOFs = subEntityDOFs(basis);
+ const auto& gridView = basis.gridView();
+ for (auto const& element : elements(gridView)) {
+ if (!element.hasBoundaryIntersections())
+ continue;
+
+ localView.bind(element);
+ for (const auto& intersection : intersections(gridView, element)) {
+ if (!this->onBoundary(intersection))
+ continue;
+
+ test_exit(intersection.neighbor(),
+ "Neighbors of periodic intersection not assigned. Use a periodic Grid instead!");
+
+ // collect DOFs of inside element on intersection
+ localView.bind(intersection.inside());
+ seDOFs.bind(localView, intersection.indexInInside(), 1);
+ std::vector<std::size_t> insideDOFs(seDOFs.begin(), seDOFs.end());
+ std::vector<MI> insideGlobalDOFs(insideDOFs.size());
+ std::transform(insideDOFs.begin(), insideDOFs.end(), insideGlobalDOFs.begin(),
+ [&](std::size_t localIndex) { return localView.index(localIndex); });
+ auto insideCoords = coords(localView.tree(), insideDOFs);
+
+ // collect DOFs of ouside element on intersection
+ localView.bind(intersection.outside());
+ seDOFs.bind(localView, intersection.indexInOutside(), 1);
+ std::vector<std::size_t> outsideDOFs(seDOFs.begin(), seDOFs.end());
+ auto outsideCoords = coords(localView.tree(), outsideDOFs);
+
+ // compare mapped coords of DOFs
+ assert(insideDOFs.size() == outsideDOFs.size());
+ for (std::size_t i = 0; i < insideCoords.size(); ++i) {
+ auto x = faceTrafo_.evaluate(insideCoords[i]);
+ for (std::size_t j = 0; j < outsideCoords.size(); ++j) {
+ auto const& y = outsideCoords[j];
+ if ((x-y).two_norm() < tol) {
+ periodicNodes_[insideGlobalDOFs[i]] = true;
+ associations_[insideGlobalDOFs[i]] = localView.index(outsideDOFs[j]);
+ }
+ }
+ }
+ }
+ }
+}
+
+namespace Impl
+{
+ template <class D>
+ class CompareTol
+ {
+ using ctype = typename D::field_type;
+
+ public:
+ CompareTol(ctype tol)
+ : tol_(tol)
+ {}
+
+ bool operator()(D const& lhs, D const& rhs) const
+ {
+ for (int i = 0; i < D::dimension; i++) {
+ if (std::abs(lhs[i] - rhs[i]) < tol_)
+ continue;
+ return lhs[i] < rhs[i];
+ }
+ return false;
+ }
+
+ private:
+ ctype tol_;
+ };
+}
+
+
+
+template <class D, class MI>
+ template <class Basis>
+void PeriodicBC<D,MI>::
+initAssociations2(Basis const& basis)
+{
+ static const double tol = std::sqrt(std::numeric_limits<typename D::field_type>::epsilon());
+ periodicNodes_.clear();
+ periodicNodes_.resize(basis.dimension(), false);
+
+ // Dune::ReservedVector<std::pair<EntitySeed,int>,2>
+ using EntitySeed = typename Basis::GridView::Grid::template Codim<0>::EntitySeed;
+ using CoordAssoc = std::map<D, std::vector<std::pair<EntitySeed,int>>, Impl::CompareTol<D>>;
+ CoordAssoc coordAssoc(Impl::CompareTol<D>{tol});
+
+ // generate periodic element pairs
+ const auto& gridView = basis.gridView();
+ for (auto const& element : elements(gridView)) {
+ for (const auto& intersection : intersections(gridView, element)) {
+ if (!this->onBoundary(intersection))
+ continue;
+
+ auto x = intersection.geometry().center();
+ auto seed = intersection.inside().seed();
+ coordAssoc[x].push_back({seed,intersection.indexInInside()});
+ coordAssoc[faceTrafo_.evaluate(x)].push_back({seed,intersection.indexInInside()});
+ }
+ }
+
+ auto localView = basis.localView();
+ auto seDOFs = subEntityDOFs(basis);
+ for (auto const& assoc : coordAssoc) {
+ auto const& seeds = assoc.second;
+ if (seeds.size() != 2)
+ continue;
+
+ // collect DOFs of inside element on intersection
+ auto el1 = gridView.grid().entity(seeds[0].first);
+ localView.bind(el1);
+ seDOFs.bind(localView, seeds[0].second, 1);
+ std::vector<std::size_t> insideDOFs(seDOFs.begin(), seDOFs.end());
+ std::vector<MI> insideGlobalDOFs(insideDOFs.size());
+ std::transform(insideDOFs.begin(), insideDOFs.end(), insideGlobalDOFs.begin(),
+ [&](std::size_t localIndex) { return localView.index(localIndex); });
+ auto insideCoords = coords(localView.tree(), insideDOFs);
+
+ // collect DOFs of ouside element on intersection
+ auto el2 = gridView.grid().entity(seeds[1].first);
+ localView.bind(el2);
+ seDOFs.bind(localView, seeds[1].second, 1);
+ std::vector<std::size_t> outsideDOFs(seDOFs.begin(), seDOFs.end());
+ auto outsideCoords = coords(localView.tree(), outsideDOFs);
+
+ // compare mapped coords of DOFs
+ assert(insideDOFs.size() == outsideDOFs.size());
+ for (std::size_t i = 0; i < insideCoords.size(); ++i) {
+ auto x = faceTrafo_.evaluate(insideCoords[i]);
+ for (std::size_t j = 0; j < outsideCoords.size(); ++j) {
+ auto const& y = outsideCoords[j];
+ if ((x-y).two_norm() < tol) {
+ periodicNodes_[insideGlobalDOFs[i]] = true;
+ associations_[insideGlobalDOFs[i]] = localView.index(outsideDOFs[j]);
+ }
+ }
+ }
+ }
+}
+
+
+template <class D, class MI>
+ template <class Node>
+std::vector<D> PeriodicBC<D,MI>::
+coords(Node const& tree, std::vector<std::size_t> const& localIndices) const
+{
+ std::vector<D> dofCoords(localIndices.size());
+ AMDiS::forEachLeafNode_(tree, [&](auto const& node, auto const& tp)
+ {
+ std::size_t size = node.finiteElement().size();
+ auto geometry = node.element().geometry();
+
+ auto const& localInterpol = node.finiteElement().localInterpolation();
+ using FiniteElement = std::decay_t<decltype(node.finiteElement())>;
+ using DomainType = typename FiniteElement::Traits::LocalBasisType::Traits::DomainType;
+ using RangeType = typename FiniteElement::Traits::LocalBasisType::Traits::RangeType;
+
+ std::array<std::vector<typename RangeType::field_type>, D::dimension> coeffs;
+ for (int d = 0; d < D::dimension; ++d) {
+ auto evalCoord = [&](DomainType const& local) -> RangeType { return geometry.global(local)[d]; };
+ using FFC = Dune::Functions::FunctionFromCallable<RangeType(DomainType), decltype(evalCoord)>;
+ localInterpol.interpolate(FFC(evalCoord), coeffs[d]);
+ }
+
+ for (std::size_t j = 0; j < localIndices.size(); ++j) {
+ D x;
+ for (std::size_t i = 0; i < size; ++i) {
+ if (node.localIndex(i) == localIndices[j]) {
+ for (int d = 0; d < D::dimension; ++d)
+ x[d] = coeffs[d][i];
+ break;
+ }
+ }
+ dofCoords[j] = std::move(x);
+ }
+ });
+
+ return dofCoords;
+}
+
+
+template <class D, class MI>
+ template <class Matrix, class VectorX, class VectorB, class RowNode, class ColNode>
+void PeriodicBC<D,MI>::
+fillBoundaryCondition(Matrix& matrix, VectorX& solution, VectorB& rhs,
+ RowNode const& rowNode, ColNode const& colNode)
+{
+ Constraints<Matrix>::periodicBC(matrix, periodicNodes_, associations_);
+
+ for (auto const& a : associations_) {
+ rhs[a.second] += rhs[a.first];
+ solution[a.second] = solution[a.first];
+ }
+ Dune::Functions::interpolate(rhs.basis(), rhs, [](auto const&) { return 0.0; }, periodicNodes_);
+}
+
+} // end namespace AMDiS
diff --git a/src/amdis/ProblemStat.hpp b/src/amdis/ProblemStat.hpp
index 2ae9e342ba9a079c929229506a60845ba0fb9482..7d028f8256e8c6ca1b20c7029b86fd9421e31bd9 100644
--- a/src/amdis/ProblemStat.hpp
+++ b/src/amdis/ProblemStat.hpp
@@ -24,6 +24,7 @@
#include <amdis/LocalAssemblerList.hpp>
#include <amdis/Marker.hpp>
#include <amdis/Mesh.hpp>
+#include <amdis/PeriodicBC.hpp>
#include <amdis/ProblemStatBase.hpp>
#include <amdis/ProblemStatTraits.hpp>
#include <amdis/StandardProblemIteration.hpp>
@@ -62,6 +63,7 @@ namespace AMDiS
using Grid = typename GridView::Grid;
using Element = typename GridView::template Codim<0>::Entity;
using WorldVector = typename Element::Geometry::GlobalCoordinate;
+ using WorldMatrix = FieldMatrix<typename WorldVector::field_type, WorldVector::dimension, WorldVector::dimension>;
/// Dimension of the grid
static constexpr int dim = Grid::dimension;
@@ -158,7 +160,7 @@ namespace AMDiS
void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath row = {}, ColTreePath col = {})
{
using I = typename GridView::Intersection;
- systemMatrix_->addOperator(tag::boundary_operator<I>{b}, op, row, col);
+ systemMatrix_->addOperator(tag::boundary_operator<I>{boundaryManager_,b}, op, row, col);
}
/** @} */
@@ -207,7 +209,7 @@ namespace AMDiS
void addVectorOperator(BoundaryType b, Operator const& op, TreePath path = {})
{
using I = typename GridView::Intersection;
- rhs_->addOperator(tag::boundary_operator<I>{b}, op, path);
+ rhs_->addOperator(tag::boundary_operator<I>{boundaryManager_,b}, op, path);
}
/** @} */
@@ -238,6 +240,15 @@ namespace AMDiS
void addDirichletBC(Predicate const& predicate,
RowTreePath row, ColTreePath col,
Values const& values);
+
+ template <class RowTreePath, class ColTreePath, class Values>
+ void addDirichletBC(BoundaryType id,
+ RowTreePath row, ColTreePath col,
+ Values const& values);
+
+ /// Add a periodic boundary conditions to the system, by specifying a face transformation
+ /// y = A*x + b of coordinates. We assume, that A is orthonormal.
+ void addPeriodicBC(BoundaryType id, WorldMatrix const& A, WorldVector const& b);
/** @} */
@@ -278,6 +289,12 @@ namespace AMDiS
/// Return the gridView of the leaf-level
GridView const& gridView() const { return globalBasis_->gridView(); }
+ /// Return the boundary manager to identify boundary segments
+ BoundaryManager<Grid>& boundaryManager() { return *boundaryManager_; }
+ BoundaryManager<Grid> const& boundaryManager() const { return *boundaryManager_; }
+
+ auto const& boundaryManagerPtr() { return boundaryManager_; }
+
/// Return the \ref globalBasis_
GlobalBasis& globalBasis() { return *globalBasis_; }
GlobalBasis const& globalBasis() const { return *globalBasis_; }
@@ -375,8 +392,15 @@ namespace AMDiS
}
void adoptGrid(std::shared_ptr<Grid> const& grid)
+ {
+ adoptGrid(grid, std::make_shared<BoundaryManager<Grid>>(grid));
+ }
+
+ void adoptGrid(std::shared_ptr<Grid> const& grid,
+ std::shared_ptr<BoundaryManager<Grid>> const& boundaryManager)
{
grid_ = grid;
+ boundaryManager_ = boundaryManager;
Parameters::get(name_ + "->mesh", gridName_);
}
@@ -416,6 +440,9 @@ namespace AMDiS
/// Name of the grid
std::string gridName_ = "mesh";
+ /// Management of boundary conditions
+ std::shared_ptr<BoundaryManager<Grid>> boundaryManager_;
+
/// FE spaces of this problem.
std::shared_ptr<GlobalBasis> globalBasis_;
@@ -448,7 +475,8 @@ namespace AMDiS
private: // some internal data-structures
- Constraints<GlobalBasis, GlobalBasis> constraints_;
+ DirichletBCs<GlobalBasis, GlobalBasis> dirichletBCs_;
+ PeriodicBCs<GlobalBasis, GlobalBasis> periodicBCs_;
};
diff --git a/src/amdis/ProblemStat.inc.hpp b/src/amdis/ProblemStat.inc.hpp
index f7d4cb1129ad53b5234babeb03293d15980b3ace..f1ff89b05ec1e286bf1b8e07ecdc36af386c3547 100644
--- a/src/amdis/ProblemStat.inc.hpp
+++ b/src/amdis/ProblemStat.inc.hpp
@@ -37,7 +37,7 @@ void ProblemStat<Traits>::initialize(
(adoptFlag.isSet(INIT_MESH) ||
adoptFlag.isSet(INIT_SYSTEM) ||
adoptFlag.isSet(INIT_FE_SPACE))) {
- grid_ = adoptProblem->grid_;
+ adoptGrid(adoptProblem->grid_, adoptProblem->boundaryManager_);
}
}
@@ -123,6 +123,7 @@ void ProblemStat<Traits>::createGrid()
{
Parameters::get(name_ + "->mesh", gridName_);
grid_ = MeshCreator<Grid>::create(gridName_);
+ boundaryManager_ = std::make_shared<BoundaryManager<Grid>>(grid_);
msg("Create grid:");
msg("#elements = {}" , grid_->size(0));
@@ -163,7 +164,8 @@ void ProblemStat<Traits>::createGlobalBasisImpl(std::false_type)
template <class Traits>
void ProblemStat<Traits>::initGlobalBasis(GlobalBasis const& globalBasis)
{
- constraints_.init(*globalBasis_, *globalBasis_);
+ dirichletBCs_.init(*globalBasis_, *globalBasis_);
+ periodicBCs_.init(*globalBasis_, *globalBasis_);
}
@@ -261,8 +263,38 @@ addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Val
using Range = RangeType_t<decltype(i)>;
using BcType = DirichletBC<WorldVector,Range>;
auto bc = std::make_unique<BcType>(predicate, valueGridFct);
- constraints_[i][j].push_back(std::move(bc));
- // TODO: make DirichletBC an abstract class and add specialization with gridfunction type
+ dirichletBCs_[i][j].push_back(std::move(bc));
+}
+
+
+// Adds a Dirichlet boundary condition
+template <class Traits>
+ template <class RowTreePath, class ColTreePath, class Values>
+void ProblemStat<Traits>::
+addDirichletBC(BoundaryType id, RowTreePath row, ColTreePath col, Values const& values)
+{
+ auto localView = globalBasis_->localView();
+ auto i = child(localView.tree(), makeTreePath(row));
+ auto j = child(localView.tree(), makeTreePath(col));
+
+ auto valueGridFct = makeGridFunction(values, globalBasis_->gridView());
+
+ using Range = RangeType_t<decltype(i)>;
+ using BcType = DirichletBC<WorldVector,Range>;
+ auto bc = std::make_unique<BcType>(boundaryManager_, id, valueGridFct);
+ dirichletBCs_[i][j].push_back(std::move(bc));
+}
+
+
+template <class Traits>
+void ProblemStat<Traits>::
+addPeriodicBC(BoundaryType id, WorldMatrix const& matrix, WorldVector const& vector)
+{
+ auto localView = globalBasis_->localView();
+ using BcType = PeriodicBC<WorldVector, typename GlobalBasis::MultiIndex>;
+ using FaceTrafo = FaceTransformation<typename WorldVector::field_type, WorldVector::dimension>;
+ auto bc = std::make_unique<BcType>(boundaryManager_, id, FaceTrafo{matrix, vector});
+ periodicBCs_[localView.tree()][localView.tree()].push_back(std::move(bc));
}
@@ -347,10 +379,14 @@ buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool as
rhs_->init(asmVector);
auto localView = globalBasis_->localView();
- forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto) {
- forEachNode_(localView.tree(), [&,this](auto const& colNode, auto) {
- for (auto bc : constraints_[rowNode][colNode])
- bc->init(*systemMatrix_, *solution_, *rhs_, rowNode, colNode);
+ forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTp) {
+ auto rowBasis = Dune::Functions::subspaceBasis(*globalBasis_, rowTp);
+ forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTp) {
+ auto colBasis = Dune::Functions::subspaceBasis(*globalBasis_, colTp);
+ for (auto bc : dirichletBCs_[rowNode][colNode])
+ bc->init(rowBasis, colBasis);
+ for (auto bc : periodicBCs_[rowNode][colNode])
+ bc->init(rowBasis, colBasis);
});
});
msg("{} total DOFs", globalBasis_->dimension());
@@ -374,8 +410,10 @@ buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool as
forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto) {
forEachNode_(localView.tree(), [&,this](auto const& colNode, auto) {
// finish boundary condition
- for (auto bc : constraints_[rowNode][colNode])
- bc->finish(*systemMatrix_, *solution_, *rhs_, rowNode, colNode);
+ for (auto bc : dirichletBCs_[rowNode][colNode])
+ bc->fillBoundaryCondition(*systemMatrix_, *solution_, *rhs_, rowNode, colNode);
+ for (auto bc : periodicBCs_[rowNode][colNode])
+ bc->fillBoundaryCondition(*systemMatrix_, *solution_, *rhs_, rowNode, colNode);
});
});
diff --git a/src/amdis/linear_algebra/CMakeLists.txt b/src/amdis/linear_algebra/CMakeLists.txt
index a436f7713b38c10fe2b334a1bc2dcd8541b5cc43..7027481fbf5e2d205a98cf50c6df0b359fb1f492 100644
--- a/src/amdis/linear_algebra/CMakeLists.txt
+++ b/src/amdis/linear_algebra/CMakeLists.txt
@@ -1,5 +1,6 @@
install(FILES
Common.hpp
+ Constraints.hpp
DOFMatrixBase.hpp
DOFMatrixBase.inc.hpp
DOFVectorBase.hpp
diff --git a/src/amdis/linear_algebra/Constraints.hpp b/src/amdis/linear_algebra/Constraints.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c774ff771eab520813a45402b5169d221dfc28eb
--- /dev/null
+++ b/src/amdis/linear_algebra/Constraints.hpp
@@ -0,0 +1,24 @@
+#pragma once
+
+#include <amdis/linear_algebra/DOFMatrixBase.hpp>
+
+namespace AMDiS
+{
+ template <class Matrix>
+ struct Constraints
+ {
+ template <class R, class C, class B, class BitVector>
+ static auto dirichletBC(DOFMatrixBase<R,C,B>& matrix, BitVector const& nodes, bool setDiagonal = true)
+ {
+ return Constraints<typename Matrix::BaseMatrix>::dirichletBC(matrix.matrix(), nodes, setDiagonal);
+ }
+
+ template <class R, class C, class B, class BitVector, class Association>
+ static auto periodicBC(DOFMatrixBase<R,C,B>& matrix, BitVector const& left, Association const& association,
+ bool setDiagonal = true)
+ {
+ return Constraints<typename Matrix::BaseMatrix>::periodicBC(matrix.matrix(), left, association, setDiagonal);
+ }
+ };
+
+} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/DOFMatrixBase.inc.hpp b/src/amdis/linear_algebra/DOFMatrixBase.inc.hpp
index 9f5e6971866d4684d61f5447f7070bc30396f544..adde13479ec5b705a466c4ef69fd1fccd20ee5f0 100644
--- a/src/amdis/linear_algebra/DOFMatrixBase.inc.hpp
+++ b/src/amdis/linear_algebra/DOFMatrixBase.inc.hpp
@@ -81,10 +81,7 @@ assemble(RowLocalView const& rowLocalView, ColLocalView const& colLocalView)
auto& matOp = operators_[rowNode][colNode];
if (matOp) {
matOp.bind(element, geometry);
- assembleOperators(gv, element, matOp, [&](auto const& context, auto& operators) {
- for (auto scaled : operators)
- scaled.op->assemble(context, rowNode, colNode, elementMatrix_);
- });
+ assembleOperators(gv, element, matOp, makeMatrixAssembler(rowNode, colNode, elementMatrix_));
matOp.unbind();
}
});
diff --git a/src/amdis/linear_algebra/DOFVectorBase.inc.hpp b/src/amdis/linear_algebra/DOFVectorBase.inc.hpp
index 884fad1d4e8dd668b17b66b350a0b8121989818a..7de4dcde61482e36d7991b920014dad31e3b3161 100644
--- a/src/amdis/linear_algebra/DOFVectorBase.inc.hpp
+++ b/src/amdis/linear_algebra/DOFVectorBase.inc.hpp
@@ -71,10 +71,7 @@ assemble(LocalView const& localView)
auto& rhsOp = operators_[node];
if (rhsOp) {
rhsOp.bind(element, geometry);
- assembleOperators(gv, element, rhsOp, [&](auto const& context, auto& operators) {
- for (auto scaled : operators)
- scaled.op->assemble(context, node, elementVector_);
- });
+ assembleOperators(gv, element, rhsOp, makeVectorAssembler(node, elementVector_));
rhsOp.unbind();
}
});
diff --git a/src/amdis/linear_algebra/eigen/CMakeLists.txt b/src/amdis/linear_algebra/eigen/CMakeLists.txt
index 41aacc5baefecd529dfada60665c42586a89aa27..8fbcf5e3f3627776a82ccb1800e3da71a85c1843 100644
--- a/src/amdis/linear_algebra/eigen/CMakeLists.txt
+++ b/src/amdis/linear_algebra/eigen/CMakeLists.txt
@@ -1,4 +1,5 @@
install(FILES
+ Constraints.hpp
DirectRunner.hpp
DOFMatrix.hpp
DOFVector.hpp
diff --git a/src/amdis/linear_algebra/eigen/Constraints.hpp b/src/amdis/linear_algebra/eigen/Constraints.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..fec59492f940bbbdb5d3c0d08b482bd089ba64c8
--- /dev/null
+++ b/src/amdis/linear_algebra/eigen/Constraints.hpp
@@ -0,0 +1,60 @@
+#pragma once
+
+#include <list>
+
+#include <Eigen/SparseCore>
+
+#include <amdis/common/Mpl.hpp>
+#include <amdis/linear_algebra/Common.hpp>
+#include <amdis/linear_algebra/Constraints.hpp>
+
+namespace AMDiS
+{
+ template <class T, int O>
+ struct Constraints<Eigen::SparseMatrix<T,O>>
+ {
+ using Matrix = Eigen::SparseMatrix<T,O>;
+
+ template <class BitVector>
+ static std::vector<Triplet<T>> dirichletBC(Matrix& mat, BitVector const& nodes, bool setDiagonal = true)
+ {
+ clearDirichletRow(mat, nodes, setDiagonal, int_<O>);
+ return {};
+ }
+
+ template <class BitVector, class Associations>
+ static std::vector<Triplet<T>> periodicBC(Matrix& mat, BitVector const& left, Associations const& left2right,
+ bool setDiagonal = true)
+ {
+ error_exit("Not implemented");
+ return {};
+ }
+
+ protected:
+ template <class BitVector>
+ static void clearDirichletRow(Matrix& mat, BitVector const& nodes, bool setDiagonal, int_t<Eigen::ColMajor>)
+ {
+ for (typename Matrix::Index c = 0; c < mat.outerSize(); ++c) {
+ for (typename Matrix::InnerIterator it(mat, c); it; ++it) {
+ if (nodes[it.row()]) {
+ it.valueRef() = (setDiagonal && it.row() == it.col() ? T(1) : T(0));
+ break;
+ }
+ }
+ }
+ }
+
+ template <class BitVector>
+ static void clearDirichletRow(Matrix& mat, BitVector const& nodes, bool setDiagonal, int_t<Eigen::RowMajor>)
+ {
+ for (typename Matrix::Index r = 0; r < mat.outerSize(); ++r) {
+ if (nodes[r]) {
+ for (typename Matrix::InnerIterator it(mat, r); it; ++it) {
+ it.valueRef() = (setDiagonal && it.row() == it.col() ? T(1) : T(0));
+ }
+ }
+ }
+ }
+ };
+
+} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/eigen/DOFMatrix.hpp b/src/amdis/linear_algebra/eigen/DOFMatrix.hpp
index 0ae8d5dfef7cad611824a2a8771ca50c4fda0b63..46920b20442500d9e650399ed402f8f47b0575a9 100644
--- a/src/amdis/linear_algebra/eigen/DOFMatrix.hpp
+++ b/src/amdis/linear_algebra/eigen/DOFMatrix.hpp
@@ -85,38 +85,6 @@ namespace AMDiS
return matrix_.nonZeros();
}
- /// \brief Deletes all rows with \p dirichletNodes != 0 and if \p apply is true, adds
- /// a one on the diagonal of the matrix.
- auto applyDirichletBC(std::vector<char> const& dirichletNodes)
- {
- for (std::size_t i = 0; i < dirichletNodes.size(); ++i) {
- if (dirichletNodes[i] == 1)
- clearDirichletRow(i, int_<Orientation>);
- }
- return std::list<Triplet<value_type>>{};
- }
-
-
- protected:
- void clearDirichletRow(size_type idx, int_t<Eigen::ColMajor>)
- {
- for (size_type k = 0; k < matrix_.outerSize(); ++k) {
- for (typename BaseMatrix::InnerIterator it(matrix_, k); it; ++it) {
- if (it.row() == idx) {
- it.valueRef() = (it.row() == it.col() ? value_type(1) : value_type(0));
- break;
- }
- }
- }
- }
-
- void clearDirichletRow(size_type idx, int_t<Eigen::RowMajor>)
- {
- Eigen::SparseVector<ValueType> unit_vector(matrix_.cols());
- unit_vector.coeffRef(idx) = value_type(1);
- matrix_.row(idx) = unit_vector;
- }
-
private:
/// The data-matrix
BaseMatrix matrix_;
diff --git a/src/amdis/linear_algebra/istl/CMakeLists.txt b/src/amdis/linear_algebra/istl/CMakeLists.txt
index d84972a0beae1be790e9e972845953fa061f363e..f342ea029c11494a5c06cf177fe6262b51e7cd26 100644
--- a/src/amdis/linear_algebra/istl/CMakeLists.txt
+++ b/src/amdis/linear_algebra/istl/CMakeLists.txt
@@ -1,4 +1,5 @@
install(FILES
+ Constraints.hpp
DirectRunner.hpp
DOFMatrix.hpp
DOFVector.hpp
diff --git a/src/amdis/linear_algebra/istl/Constraints.hpp b/src/amdis/linear_algebra/istl/Constraints.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5a11d6ff7d1e47cccdb79063e6bd446ca23a08a1
--- /dev/null
+++ b/src/amdis/linear_algebra/istl/Constraints.hpp
@@ -0,0 +1,43 @@
+#pragma once
+
+#include <list>
+
+#include <dune/istl/bcrsmatrix.hh>
+
+#include <amdis/linear_algebra/Common.hpp>
+#include <amdis/linear_algebra/Constraints.hpp>
+
+namespace AMDiS
+{
+ template <class T, class A>
+ struct Constraints<Dune::BCRSMatrix<T,A>>
+ {
+ using Matrix = Dune::BCRSMatrix<T,A>;
+
+ template <class BitVector>
+ static std::vector<Triplet<T>> dirichletBC(Matrix& mat, BitVector const& nodes, bool setDiagonal = true)
+ {
+ // loop over the matrix rows
+ for (std::size_t i = 0; i < mat.N(); ++i) {
+ if (nodes[i]) {
+ auto cIt = mat[i].begin();
+ auto cEndIt = mat[i].end();
+ // loop over nonzero matrix entries in current row
+ for (; cIt != cEndIt; ++cIt)
+ *cIt = (setDiagonal && i == cIt.index() ? T(1) : T(0));
+ }
+ }
+
+ return {};
+ }
+
+ template <class BitVector, class Associations>
+ static std::vector<Triplet<T>> periodicBC(Matrix& mat, BitVector const& left, Associations const& left2right,
+ bool setDiagonal = true)
+ {
+ error_exit("Not implemented");
+ return {};
+ }
+ };
+
+} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/istl/DOFMatrix.hpp b/src/amdis/linear_algebra/istl/DOFMatrix.hpp
index 678ef0b58b2d970db5e3ac9b33c6e079d8dfc8e1..572d2bd9bd01a262cc85ec143bf5c96541834700 100644
--- a/src/amdis/linear_algebra/istl/DOFMatrix.hpp
+++ b/src/amdis/linear_algebra/istl/DOFMatrix.hpp
@@ -102,21 +102,6 @@ namespace AMDiS
return matrix_.nonzeroes();
}
- auto applyDirichletBC(std::vector<char> const& dirichletNodes)
- {
- // loop over the matrix rows
- for (std::size_t i = 0; i < matrix_.N(); ++i) {
- if (dirichletNodes[i]) {
- auto cIt = matrix_[i].begin();
- auto cEndIt = matrix_[i].end();
- // loop over nonzero matrix entries in current row
- for (; cIt != cEndIt; ++cIt)
- *cIt = (i == cIt.index() ? 1 : 0);
- }
- }
- return std::list<Triplet<value_type>>{};
- }
-
private:
BaseMatrix matrix_;
diff --git a/src/amdis/linear_algebra/mtl/CMakeLists.txt b/src/amdis/linear_algebra/mtl/CMakeLists.txt
index e295f73f7767c0b95987825dd993f89720bc30f5..adcf1b591b2a56d86c1c884b5c1d03bc1e5af7a0 100644
--- a/src/amdis/linear_algebra/mtl/CMakeLists.txt
+++ b/src/amdis/linear_algebra/mtl/CMakeLists.txt
@@ -1,4 +1,5 @@
install(FILES
+ Constraints.hpp
DOFMatrix.hpp
DOFVector.hpp
ITL_Preconditioner.hpp
diff --git a/src/amdis/linear_algebra/mtl/Constraints.hpp b/src/amdis/linear_algebra/mtl/Constraints.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..bd20520609e12d897ee2b149ed48cdb3c4e3b53c
--- /dev/null
+++ b/src/amdis/linear_algebra/mtl/Constraints.hpp
@@ -0,0 +1,146 @@
+#pragma once
+
+#include <list>
+
+#include <boost/numeric/mtl/matrix/compressed2D.hpp>
+#include <boost/numeric/mtl/matrix/inserter.hpp>
+#include <boost/numeric/mtl/utility/property_map.hpp>
+#include <boost/numeric/mtl/utility/range_wrapper.hpp>
+
+#include <amdis/linear_algebra/Common.hpp>
+#include <amdis/linear_algebra/Constraints.hpp>
+
+namespace AMDiS
+{
+ template <class T, class P>
+ struct Constraints<mtl::compressed2D<T,P>>
+ {
+ using Matrix = mtl::compressed2D<T,P>;
+
+#ifdef SYMMETRIC_DIRICHLET_BC
+ template <class BitVector>
+ static std::vector<Triplet<T>> dirichletBC(Matrix& mat, BitVector const& nodes, bool setDiagonal = true)
+ {
+ std::vector<Triplet<T>> columns;
+ if (setDiagonal)
+ columns.reserve(std::size_t(mat.nnz()/(0.84*num_rows(mat))));
+
+ // Define the property maps
+ auto row = mtl::mat::row_map(mat);
+ auto col = mtl::mat::col_map(mat);
+ auto value = mtl::mat::value_map(mat);
+
+ // iterate over the matrix
+ std::size_t slotSize = 0;
+ for (auto r : mtl::rows_of(mat)) { // rows or columns
+ std::size_t rowSize = 0;
+ for (auto i : mtl::nz_of(r)) { // non-zeros within
+ ++rowSize;
+ if (nodes[row(i)]) {
+ // set identity row
+ value(i, T(0));
+ }
+ else if (setDiagonal && nodes[col(i)]) {
+ columns.push_back({row(i), col(i), value(i)});
+ value(i, T(0));
+ }
+ }
+ slotSize = std::max(slotSize, rowSize);
+ }
+
+ // set diagonal entry
+ if (setDiagonal) {
+ mtl::mat::inserter<Matrix, mtl::update_store<T> > ins(mat, slotSize);
+ for (std::size_t i = 0; i < nodes.size(); ++i) {
+ if (nodes[i])
+ ins[i][i] = T(1);
+ }
+ }
+
+ return columns;
+ }
+#else
+ template <class BitVector>
+ static std::vector<Triplet<T>> dirichletBC(Matrix& mat, BitVector const& nodes, bool setDiagonal = true)
+ {
+ // Define the property maps
+ auto row = mtl::mat::row_map(mat);
+ auto col = mtl::mat::col_map(mat);
+ auto value = mtl::mat::value_map(mat);
+
+ // iterate over the matrix
+ for (auto r : mtl::rows_of(mat)) { // rows of the matrix
+ if (nodes[r.value()]) {
+ for (auto i : mtl::nz_of(r)) { // non-zeros within
+ // set identity row
+ value(i, (setDiagonal && row(i) == col(i) ? T(1) : T(0)) );
+ }
+ }
+ }
+
+ return {};
+ }
+#endif
+
+
+ template <class Associations>
+ static std::size_t at(Associations const& m, std::size_t idx)
+ {
+ auto it = m.find(idx);
+ assert(it != m.end());
+ return it->second;
+ }
+
+#ifdef SYMMETRIC_PERIODIC_BC
+ template <class BitVector, class Associations>
+ static std::vector<Triplet<T>> periodicBC(Matrix& mat, BitVector const& left, Associations const& left2right,
+ bool setDiagonal = true)
+ {
+ error_exit("Not implemented");
+ }
+#else
+ template <class BitVector, class Associations>
+ static std::vector<Triplet<T>> periodicBC(Matrix& mat, BitVector const& left, Associations const& left2right,
+ bool setDiagonal = true)
+ {
+ std::vector<Triplet<T>> rowValues;
+ rowValues.reserve(left2right.size()*std::size_t(mat.nnz()/(0.9*num_rows(mat))));
+
+ // Define the property maps
+ // auto row = mtl::mat::row_map(mat);
+ auto col = mtl::mat::col_map(mat);
+ auto value = mtl::mat::value_map(mat);
+
+ // iterate over the matrix
+ std::size_t slotSize = 0;
+ for (auto r : mtl::rows_of(mat)) {
+ if (left[r.value()]) {
+ slotSize = std::max(slotSize, std::size_t(mat.nnz_local(r.value())));
+ std::size_t right = at(left2right,r.value());
+
+ for (auto i : mtl::nz_of(r)) {
+ rowValues.push_back({right,col(i),value(i)});
+ value(i, T(0));
+ }
+ }
+ }
+
+ mtl::mat::inserter<Matrix, mtl::update_plus<T> > ins(mat, 2*slotSize);
+ for (auto const& t : rowValues)
+ ins[t.row][t.col] += t.value;
+
+ if (setDiagonal) {
+ for (std::size_t i = 0; i < size(left); ++i) {
+ if (left[i]) {
+ ins[i][i] = T(1);
+ ins[i][at(left2right,i)] = T(-1);
+ }
+ }
+ }
+
+ return {};
+ }
+#endif
+ };
+
+} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/DOFMatrix.hpp b/src/amdis/linear_algebra/mtl/DOFMatrix.hpp
index f1213b7e06ce5f68de39b9bf2d4c68c740667514..26ecdf75afe01b658465bc5086a0932ae78924c4 100644
--- a/src/amdis/linear_algebra/mtl/DOFMatrix.hpp
+++ b/src/amdis/linear_algebra/mtl/DOFMatrix.hpp
@@ -94,28 +94,6 @@ namespace AMDiS
return matrix_.nnz();
}
- /// \brief Deletes all rows with \p dirichletNodes != 0 and if \p apply is true, adds
- /// a one on the diagonal of the matrix.
- auto applyDirichletBC(std::vector<char> const& dirichletNodes)
- {
- // Define the property maps
- auto row = mtl::mat::row_map(matrix_);
- auto col = mtl::mat::col_map(matrix_);
- auto value = mtl::mat::value_map(matrix_);
-
- // iterate over the matrix
- for (auto r : mtl::rows_of(matrix_)) { // rows of the matrix
- if (dirichletNodes[r.value()]) {
- for (auto i : mtl::nz_of(r)) { // non-zeros within
- // set identity row
- value(i, (row(i) == col(i) ? value_type(1) : value_type(0)) );
- }
- }
- }
-
- return std::list<Triplet<value_type>>{};
- }
-
protected:
// Estimates the slot-size used in the inserter to reserve enough space per row.
void calculateSlotSize()