Commit 51438329 authored by Praetorius, Simon's avatar Praetorius, Simon

intersection and boundary operators added

parent a9f97a0f
Pipeline #907 passed with stage
in 13 minutes and 4 seconds
......@@ -16,6 +16,8 @@ namespace AMDiS
using Geometry = typename Impl::Get<Element>::Geometry;
public:
/// Constructor, initializes the geometry of the element and a
/// quadrature-rule wrapper
template <class Operator>
Assembler(Operator& op, Element const& element, int degree, FirstOrderType type = GRD_PHI)
: geometry(get_geometry(element))
......@@ -26,24 +28,45 @@ namespace AMDiS
quadrature.reset(new QuadratureRule(element, quad));
}
/// \brief Returns reference to the transformed quadrature rule
/**
* For intersections this corresponds to the points in the intersection
* geometry, but coordinates extended to the whole inside-entity. For
* elements this is a wrapper around the classical element quadrature rule.
* Access the underlying dune-quadrature rule, with `quadrature->getRule()`.
**/
auto const& getQuadrature() const
{
return quadrature->getRule();
}
auto const& getTransformedQuadrature() const
{
return *quadrature;
}
/// Return the geometry of the Object
/**
* The geometry is either the element geometry, or the geometry of the
* inside-element for intersections.
**/
Geometry const& getGeometry() const
{
return geometry;
}
protected:
// transform coords from intersection coords to element coords
template <class Coordinate>
decltype(auto) map(Coordinate const& p) const
{
return get_position<Element>(geometry, p);
}
private:
/// transformed quadrature rule
std::shared_ptr<QuadratureRule const> quadrature;
Geometry geometry; // TODO: create geometry just once for each element, e.g. in ProblemStat when traversing the grid
/// geometry() for entities or geometryInInside() for intersections
Geometry geometry;
// TODO: create geometry just once for each element, e.g. in ProblemStat when traversing the grid
};
} // end namespace AMDiS
#pragma once
namespace AMDiS
{
template <class MeshView>
class BoundaryFacetIterator;
// An Iterator over all elements and when element hasBoundaryIntersections
template <class MeshView>
class BoundaryElementIterator
{
friend class BoundaryFacetIterator<MeshView>;
using Self = BoundaryElementIterator;
using Element = typename MeshView::template Codim<0>::Entity;
using ElementIterator = typename MeshView::template Codim<0>::Iterator;
class Iterator
{
public:
Iterator(ElementIterator elementIt, ElementIterator endIt)
: elementIt(elementIt)
, endIt(endIt)
{}
Iterator(Iterator const&) = default;
Iterator& operator=(Iterator const&) = default;
// iterate to next element that has boundary intersections
Iterator& operator++()
{
++elementIt;
while (elementIt != endIt && !elementIt->hasBoundaryIntersections())
++elementIt;
return *this;
}
Iterator operator++(int)
{
auto tmp = *this;
++(*this);
return tmp;
}
Element const& operator*() const
{
return *elementIt;
}
Element const* operator->() const
{
return &(*elementIt);
}
bool operator==(Iterator const& that) const
{
return elementIt == that.elementIt;
}
bool operator!=(Iterator const& that) const
{
return !(*this == that);
}
private:
ElementIterator elementIt;
ElementIterator endIt;
};
public:
/// Constructor.
BoundaryElementIterator(MeshView const& meshView)
: meshView(meshView)
{}
Iterator begin() {
auto elementIt = elements(meshView).begin();
auto endIt = elements(meshView).end();
while (elementIt != endIt && !elementIt->hasBoundaryIntersections())
++elementIt;
return {elementIt, endIt};
}
Iterator end() {
return {elements(meshView).end(), elements(meshView).end()};
}
private:
MeshView const& meshView;
};
/// Generator function for the boundary-element iterator
template <class MeshView>
BoundaryElementIterator<MeshView> boundary_elements(MeshView const& meshView)
{
return {meshView};
}
// An Iterator over all elements and when element hasBoundaryIntersections, then
// iterate over all boundary-intersections with given Index, oder for thos with
// predicate returns true
template <class MeshView>
class BoundaryFacetIterator
{
using Element = typename MeshView::template Codim<0>::Entity;
using Facet = typename MeshView::template Codim<1>::Entity;
using ElementIterator = typename BoundaryElementIterator<MeshView>::Iterator;
using FacetIterator = typename MeshView::IntersectionIterator;
class Iterator
{
public:
Iterator(MeshView const& meshView,
ElementIterator elementIt,
FacetIterator facetIt,
ElementIterator endIt)
: meshView(meshView)
, elementIt(elementIt)
, facetIt(facetIt)
, endIt(endIt)
{}
Iterator(Iterator const&) = default;
Iterator& operator=(Iterator const&) = default;
// iterate to next boundary face within current element or to the first
// boundary face in the next boundary element
Iterator& operator++()
{
++facetIt;
do {
auto facetEndIt = intersections(meshView, *elementIt).end();
while (facetIt != facetEndIt && !facetIt->boundary())
++facetIt;
if (facetIt != facetEndIt)
break;
++elementIt;
facetIt = intersections(meshView, *elementIt).begin();
} while (elementIt != endIt);
return *this;
}
Iterator operator++(int)
{
auto tmp = *this;
++(*this);
return tmp;
}
// returns Dune::Intersection
auto const& operator*() const
{
return *facetIt;
}
auto const* operator->() const
{
return &(*facetIt);
}
bool operator==(Iterator const& that) const
{
return elementIt == that.elementIt && (elementIt == endIt || facetIt == that.facetIt);
}
bool operator!=(Iterator const& that) const
{
return !(*this == that);
}
private:
MeshView const& meshView;
ElementIterator elementIt;
FacetIterator facetIt;
ElementIterator endIt;
};
public:
/// Constructor.
BoundaryFacetIterator(MeshView const& meshView)
: meshView(meshView)
{}
Iterator begin() {
auto elementIt = boundary_elements(meshView).begin();
auto endElementIt = boundary_elements(meshView).end();
auto facetIt = intersections(meshView, *elementIt).begin();
while (!facetIt->boundary())
++facetIt;
return {meshView, elementIt, facetIt, endElementIt};
}
Iterator end() {
auto elementIt = boundary_elements(meshView).begin();
auto endElementIt = boundary_elements(meshView).end();
auto facetIt = intersections(meshView, *elementIt).begin();
return {meshView, endElementIt, facetIt, endElementIt};
}
private:
MeshView const& meshView;
};
/// Generator function for the boundary-element iterator
template <class MeshView>
BoundaryFacetIterator<MeshView> boundary_facets(MeshView const& meshView)
{
return {meshView};
}
} // end namespace AMDiS
......@@ -8,9 +8,9 @@ namespace AMDiS
template <class WorldVector>
template <class Matrix, class VectorX, class VectorB>
void DirichletBC<WorldVector>::init(bool apply,
Matrix& matrix,
Matrix& matrix,
VectorX& solution,
VectorB& rhs)
VectorB& rhs)
{
using Dune::Functions::interpolate;
......@@ -24,13 +24,13 @@ namespace AMDiS
template <class WorldVector>
template <class Matrix, class VectorX, class VectorB>
void DirichletBC<WorldVector>::finish(bool apply,
Matrix& matrix,
Matrix& matrix,
VectorX& solution,
VectorB& rhs)
VectorB& rhs)
{
using Dune::Functions::interpolate;
test_exit(initialized, "Boundary condition not initialized!");
test_exit_dbg(initialized, "Boundary condition not initialized!");
matrix.clearDirichletRows(dirichletNodes, apply);
if (apply) {
......
......@@ -37,9 +37,9 @@ namespace AMDiS
public:
/// Constructor.
FileWriter(std::string base, MeshView const& meshView, std::vector<std::string> const& names_)
FileWriter(std::string base, MeshView const& meshView, std::vector<std::string> const& names)
: meshView(meshView)
, names(names_)
, names(names)
{
filename = "solution";
Parameters::get(base + "->filename", filename);
......
......@@ -27,23 +27,23 @@ namespace AMDiS
: Super(op, element, 1, type)
{
if (type == GRD_PHI)
op.initFot1(element, Super::getTransformedQuadrature());
op.initFot1(element, Super::getQuadrature());
else
op.initFot2(element, Super::getTransformedQuadrature());
op.initFot2(element, Super::getQuadrature());
}
// tag dispatching for FirstOrderType...
template <class Operator, class RowView, class ColView>
void calculateElementMatrix(Operator& op, RowView const& rowView, ColView const& colView,
ElementMatrix& elementMatrix)
void calculateElementMatrix(Operator& op, RowView const& rowView, ColView const& colView, ElementMatrix& elementMatrix)
{
calculateElementMatrix(op, rowView, colView, elementMatrix, std::integral_constant<FirstOrderType, type>{});
calculateElementMatrix(op, rowView, colView, elementMatrix, FirstOrderType_<type>);
}
template <class Operator, class RowView>
void calculateElementVector(Operator& op, RowView const& rowView, ElementVector& elementVector)
{
calculateElementVector(op, rowView, elementVector, std::integral_constant<FirstOrderType, type>{});
calculateElementVector(op, rowView, elementVector, FirstOrderType_<type>);
}
......@@ -52,18 +52,18 @@ namespace AMDiS
RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix,
std::integral_constant<FirstOrderType, GRD_PHI>)
FirstOrderType_t<GRD_PHI>)
{
auto geometry = rowView.element().geometry();
auto const& quad_geometry = Super::getGeometry();
auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
auto const& quad = Super::getQuadrature();
auto const& quad = Super::getQuadrature().getRule();
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
decltype(auto) quadPos = get_position<Element>(quad_geometry, quad[iq].position());
decltype(auto) quadPos = Super::map(quad[iq].position());
// The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
......@@ -100,18 +100,18 @@ namespace AMDiS
RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix,
std::integral_constant<FirstOrderType, GRD_PSI>)
FirstOrderType_t<GRD_PSI>)
{
auto geometry = rowView.element().geometry();
auto const& quad_geometry = Super::getGeometry();
auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
auto const& quad = Super::getQuadrature();
auto const& quad = Super::getQuadrature().getRule();
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
decltype(auto) quadPos = get_position<Element>(quad_geometry, quad[iq].position());
decltype(auto) quadPos = Super::map(quad[iq].position());
// The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
......@@ -147,17 +147,17 @@ namespace AMDiS
void calculateElementVector(Operator& op,
RowView const& rowView,
ElementVector& elementVector,
std::integral_constant<FirstOrderType, GRD_PSI>)
FirstOrderType_t<GRD_PSI>)
{
auto geometry = rowView.element().geometry();
auto const& quad_geometry = Super::getGeometry();
auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
auto const& quad = Super::getQuadrature();
auto const& quad = Super::getQuadrature().getRule();
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
decltype(auto) quadPos = get_position<Element>(quad_geometry, quad[iq].position());
decltype(auto) quadPos = Super::map(quad[iq].position());
// The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
......@@ -183,6 +183,4 @@ namespace AMDiS
}
};
} // end namespace AMDiS
......@@ -12,28 +12,33 @@ namespace AMDiS
{
public:
/// Constructs a unset Flag
Flag() : flags(0) {}
constexpr Flag() = default;
/// Constructs a Flag initialized by f
Flag(const unsigned long f) : flags(f) {}
constexpr Flag(const unsigned long f)
: flags(f)
{}
/// Copy constructor
Flag(Flag const& f) : flags(f.flags) {}
constexpr Flag(Flag const&) = default;
/// Move constructor
constexpr Flag(Flag&&) = default;
/// Compares two Flags
bool operator==(Flag const& f) const
constexpr bool operator==(Flag const& f) const
{
return (flags == f.flags);
}
/// Compares two Flags
bool operator!=(Flag const& f) const
constexpr bool operator!=(Flag const& f) const
{
return !(f == *this);
}
/// Assignment operator
Flag& operator=(Flag const& f)
constexpr Flag& operator=(Flag const& f)
{
if (this != &f)
flags = f.flags;
......@@ -41,101 +46,96 @@ namespace AMDiS
}
/// Typecast
operator bool() const
constexpr operator bool() const
{
return isAnySet();
}
/// Set \ref flags
void setFlags(const unsigned long f)
constexpr void setFlags(const unsigned long f)
{
flags = f;
}
/// Set \ref flags
void setFlags(Flag const& f)
constexpr void setFlags(Flag const& f)
{
flags = f.flags;
}
/// Sets \ref flags to \ref flags | f
void setFlag(const unsigned long f)
constexpr void setFlag(const unsigned long f)
{
flags |= f;
}
/// Sets \ref flags to \ref flags | f.flags
void setFlag(Flag const& f)
constexpr void setFlag(Flag const& f)
{
flags |= f.flags;
}
/// Sets \ref flags to \ref flags & ~f
void unsetFlag(const unsigned long f)
constexpr void unsetFlag(const unsigned long f)
{
flags &= ~f;
}
/// Sets \ref flags to \ref flags & ~f.flags
void unsetFlag(Flag const& f)
constexpr void unsetFlag(Flag const& f)
{
flags &= ~f.flags;
}
unsigned long getFlags() const
constexpr unsigned long getFlags() const
{
return flags;
}
/// Returns \ref flags | f.flags
Flag operator+(Flag const& f) const
friend Flag operator+(Flag r, Flag const& f)
{
Flag r(flags);
r.setFlag(f);
return r;
}
/// Returns \ref flags & ~f.flags
Flag operator-(Flag const& f) const
friend Flag operator-(Flag r, Flag const& f)
{
Flag r(flags);
r.unsetFlag(f);
return r;
}
/// Returns \ref flags | f.flags
Flag operator|(Flag const& f) const
friend Flag operator|(Flag r, Flag const& f)
{
Flag r(flags);
r.setFlag(f);
return r;
}
/// Returns \ref flags & f.flags
Flag operator&(Flag const& f) const
friend Flag operator&(Flag r, Flag const& f)
{
Flag r(flags);
r.flags &= f.flags;
return r;
}
/// Sets \ref flags to \ref flags &= f.flags
Flag operator&=(Flag const& f)
constexpr Flag operator&=(Flag const& f)
{
flags &= f.flags;
return *this;
}
/// Returns \ref flags ^ f.flags
Flag operator^(Flag const& f) const
friend Flag operator^(Flag r, Flag const& f)
{
Flag r(flags);
r.flags ^= f.flags;
return r;
}
/// Sets \ref flags to \ref flags & f.flags
Flag& operator|=(Flag const& f)
constexpr Flag& operator|=(Flag const& f)
{
if (this != &f)
flags |= f.flags;
......@@ -143,7 +143,7 @@ namespace AMDiS
}
/// Returns ~\ref flags
Flag operator~() const
constexpr Flag operator~() const
{
Flag r;
r.flags = ~flags;
......@@ -151,26 +151,26 @@ namespace AMDiS
}
/// Checks whether all set bits of f.flags are set in \ref flags too.
bool isSet(Flag const& f) const
constexpr bool isSet(Flag const& f) const
{
return ((flags&f.flags) == f.flags);
}
/// Returns !\ref isSet(f)
bool isUnset(Flag const& f) const
constexpr bool isUnset(Flag const& f) const
{
return !isSet(f);
}
/// Returns true if \ref flags != 0
bool isAnySet() const
constexpr bool isAnySet() const
{
return (flags != 0);
}
protected:
/// Internal flag representation
unsigned long flags;
unsigned long flags = 0;
};
} // end namespace AMDiS
......@@ -51,22 +51,22 @@ namespace AMDiS
/// Add coefficients for first-order operator < psi, term * grad(phi) > or
/// < grad(psi), term * phi >, where the first corresponds to
/// \p firstOrderType == GRD_PHI and the second one to \p firstOrderType == GRD_PSI.
/// \p type == GRD_PHI and the second one to \p type == GRD_PSI.
/// The coefficient \p b must be a vector expression.
template <class Term>
Self& addFOT(Term const& b, FirstOrderType firstOrderType)
Self& addFOT(Term const& b, FirstOrderType type)
{
return addFOTImpl(toTerm(b), firstOrderType);
return addFOTImpl(toTerm(b), type);
}
/// Add coefficients for first-order operator < psi, b * d_i(phi) > or
/// < d_i(psi), b * phi >, where the first corresponds to
/// \p firstOrderType == GRD_PHI and the second one to \p firstOrderType == GRD_PSI.
/// \p type == GRD_PHI and the second one to \p type == GRD_PSI.
/// The coefficient \p b must be a scalar expression.
template <class Term>
Self& addFOT(Term const& b, std::size_t i, FirstOrderType firstOrderType)
Self& addFOT(Term const& b, std::size_t i, FirstOrderType type)
{
return addFOTImpl(toTerm(b), i, firstOrderType);
return addFOTImpl(toTerm(b), i, type);
}
template <class... Args>
......@@ -105,14 +105,7 @@ namespace AMDiS
}
/// Calls \ref zot(), \ref for() or \ref sot(), depending on template
/// parameter \p Order.
template <std::size_t Order, class... Args>
static shared_ptr<Self> create(Args&&... args)
{
return create(index_<Order>, std::forward<Args>(args)...);
}
public: // initialize and assemble operator on element
/// Extract the polynomial degree from \p rowFeSpace and \p colFeSpace.
template <class RowFeSpace, class ColFeSpace>
......@@ -141,7 +134,9 @@ namespace AMDiS
ElementVector& elementVector,
double* factor = NULL);
private:
private: // implementation details
template <class Term>
Self& addZOTImpl(Term const& term);
......@@ -158,25 +153,7 @@ namespace AMDiS
Self& addSOTImpl(Term const& term, std::size_t i, std::size_t j);
template <class... Args>
static shared_ptr<Self> create(index_t<0>, Args&&... args)
{
return zot(std::forward<Args>(args)...);
}
template <class... Args>
static shared_ptr<Self> create(index_t<1>, Args&&... args)
{
return fot(std::forward<Args>(args)...);
}
template <class... Args>
static shared_ptr<Self> create(index_t<2>, Args&&... args)
{
return sot(std::forward<Args>(args)...);
}
protected:
protected: // sum up constribution from all operator-terms
template <class... Args>
void evalZot(double& result, int iq, double factor, Args const&... args)
......@@ -207,6 +184,8 @@ namespace AMDiS
}
protected: // initialize opoerator-terms
template <class E, class Quadrature>
void initZot(E const& element, Quadrature const& quad)
{
......@@ -235,11 +214,15 @@ namespace AMDiS
operatorTerm->init(element, quad);
}