Commit 7800e1b5 authored by Praetorius, Simon's avatar Praetorius, Simon

Boundary conditions and Constraints

parent 4fddde92
......@@ -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
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)
......
#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;
}
......@@ -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");
......
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
#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;
}
#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;
}
......@@ -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
......@@ -2,6 +2,6 @@
namespace AMDiS
{
struct BoundaryType { int b; };
using BoundaryType = int;
} // end namespace AMDiS
#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
#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 :