Commit 14974267 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Fix a wrong behavior in the ContextGeometry implementation for Intersections

parent b0f55468
......@@ -10,11 +10,67 @@ namespace AMDiS
{
namespace Impl
{
template <class ct, int dim>
class IdentityGeometry
{
public:
using ctype = ct;
using Volume = ct;
static const int mydimension = dim;
static const int coorddimension = dim;
using Coordinate = Dune::FieldVector<ctype, dim>;
using LocalCoordinate = Coordinate;
using GlobalCoordinate = Coordinate;
using IdentityMatrix = Dune::DiagonalMatrix<ctype, dim>;
using JacobianTransposed = IdentityMatrix;
using JacobianInverseTransposed = IdentityMatrix;
private:
using ReferenceElement = Dune::Geo::ReferenceElement<Dune::Geo::ReferenceElementImplementation<ct, dim>>;
using ReferenceElements = Dune::Geo::ReferenceElements<ct, dim>;
public:
IdentityGeometry (const ReferenceElement& refElement)
: refElement_(refElement)
, identityMatrix_(ctype(1))
{}
IdentityGeometry (Dune::GeometryType gt)
: IdentityGeometry(ReferenceElements::general(gt))
{}
bool affine () const { return true; }
Dune::GeometryType type () const { return refElement_.type(); }
int corners () const { return refElement_.size(dim); }
Coordinate corner (int i) const { return refElement_.position(i, dim); }
Coordinate center () const { return refElement_.position(0, 0); }
const Coordinate& global (const Coordinate& local) const { return local; }
const Coordinate& local (const Coordinate& global) const { return global; }
ctype integrationElement (const Coordinate& /*local*/) const { return ctype(1); }
Volume volume () const { return refElement_.volume(); }
const IdentityMatrix& jacobianTransposed (const Coordinate& /*local*/) const { return identityMatrix_; }
const IdentityMatrix& jacobianInverseTransposed (const Coordinate& /*local*/) const { return identityMatrix_; }
friend ReferenceElement referenceElement (const IdentityGeometry& geometry) { return geometry.refElement_; }
private:
ReferenceElement refElement_;
IdentityMatrix identityMatrix_;
};
template <class E, class = std::void_t<>>
struct ContextTypes
{
using Entity = E;
using LocalGeometry = typename E::Geometry;
using Geometry = typename E::Geometry;
using LocalGeometry = IdentityGeometry<typename Geometry::ctype, Geometry::mydimension>;
};
// specialization for intersections
......@@ -47,7 +103,8 @@ namespace AMDiS
public:
using LocalContext = LC;
using Element = typename ContextType::Entity;
using Geometry = typename Element::Geometry;
using ElementGeometry = typename Element::Geometry;
using Geometry = typename ContextType::Geometry;
using LocalGeometry = typename ContextType::LocalGeometry;
using IsEntity = std::is_same<Element, LocalContext>;
......@@ -61,7 +118,7 @@ namespace AMDiS
ContextGeometry(LC const& localContext, Element const& element, Geometry const& geometry)
: localContext_(&localContext)
, element_(&element)
, geometry_(&geometry)
, elementGeometry_(&geometry)
{}
// Prevent from storing pointer to rvalue-reference
......@@ -80,33 +137,51 @@ namespace AMDiS
return *localContext_;
}
/// Return the geometry of the \ref Element
/// Geometry of the local context
Geometry const& geometry() const
{
return *geometry_;
if constexpr (IsEntity{})
return elementGeometry_;
else {
if (!geometry_)
geometry_.emplace(localContext_->geometry());
return *geometry_;
}
}
/// Return the geometry of the element, or geometryInInside of the intersection
LocalGeometry const& localGeometry() const
/// Return the geometry transform of the context in the element
decltype(auto) geometryInElement() const
{
return localGeometry_impl(IsEntity{});
if constexpr(IsEntity{})
return LocalGeometry{elementGeometry_.type()};
else {
if (!localGeometry_)
localGeometry_.emplace(localContext_->geometryInInside());
return *localGeometry_;
}
}
public:
/// Return the geometry of the \ref Element
ElementGeometry const& elementGeometry() const
{
return *elementGeometry_;
}
/// Coordinate `p` given in `localGeometry`, transformed to coordinate in geometry of the LocalContext.
template <class Coordinate>
decltype(auto) local(Coordinate const& p) const
decltype(auto) coordInElement(Coordinate const& p) const
{
return local_impl(p, IsEntity{});
if constexpr(IsEntity{})
return p;
else
return geometryInElement().global(p);
}
/// Transformation of coordinate `p` given in `localGeometry` to world space coordinates.
template <class Coordinate>
decltype(auto) global(Coordinate const& p) const
{
return geometry_->global(p);
return elementGeometry_->global(coordInElement(p));
}
/// Return the geometry-type of the localContext
......@@ -115,52 +190,21 @@ namespace AMDiS
return localContext_->type();
}
/// The integration element from the `localGeometry`, the quadrature points are
/// defined in.
/// The integration element geometry of the localContext
template <class Coordinate>
auto integrationElement(Coordinate const& p) const
{
return localGeometry().integrationElement(p);
}
private: // implementation detail
// position for elements
template <class Coordinate>
Coordinate const& local_impl(Coordinate const& p, std::true_type) const
{
return p;
}
// position for intersection
template <class Coordinate>
auto local_impl(Coordinate const& p, std::false_type) const
{
return localGeometry().global(p);
}
// local-geometry is the same as geometry
Geometry const& localGeometry_impl(std::true_type) const
{
return *geometry_;
}
// local-geometry of intersection in inside element
LocalGeometry const& localGeometry_impl(std::false_type) const
{
if (!localGeometry_)
localGeometry_.emplace(localContext_->geometryInInside());
return *localGeometry_;
return geometry().integrationElement(p);
}
private:
LocalContext const* localContext_;
Element const* element_;
Geometry const* geometry_;
ElementGeometry const* elementGeometry_;
// The localGeometry may be constructed only if needed
mutable std::optional<LocalGeometry> localGeometry_;
mutable std::optional<Geometry> geometry_;
};
......@@ -182,7 +226,7 @@ namespace AMDiS
GlobalContext(GV const& gridView, E&& element, G&& geometry)
: gridView_(gridView)
, element_(Dune::wrap_or_move(FWD(element)))
, geometry_(Dune::wrap_or_move(FWD(geometry)))
, elementGeometry_(Dune::wrap_or_move(FWD(geometry)))
{}
public:
......@@ -201,13 +245,13 @@ namespace AMDiS
/// Return the geometry of the \ref Element
Geometry const& geometry() const
{
return *geometry_;
return *elementGeometry_;
}
private:
GridView gridView_;
std::shared_ptr<Element const> element_;
std::shared_ptr<Geometry const> geometry_;
std::shared_ptr<Geometry const> elementGeometry_;
};
} // end namespace AMDiS
......@@ -147,7 +147,7 @@ namespace AMDiS
void assemble(CG const& contextGeo, RN const& rowNode, CN const& colNode,
Mat& elementMatrix) const
{
auto const& quad = getQuadratureRule(contextGeo.localContext().geometry(),
auto const& quad = getQuadratureRule(contextGeo.coordInElementContext().geometry(),
derivDeg_, localFctOrder(), rowNode, colNode);
impl().assemble(contextGeo, rowNode, colNode, quad, localFct_, elementMatrix);
}
......@@ -161,7 +161,7 @@ namespace AMDiS
void assemble(CG const& contextGeo, Node const& node,
Vec& elementVector) const
{
auto const& quad = getQuadratureRule(contextGeo.localContext().geometry(),
auto const& quad = getQuadratureRule(contextGeo.coordInElementContext().geometry(),
derivDeg_, localFctOrder(), node);
impl().assemble(contextGeo, node, quad, localFct_, elementVector);
}
......
......@@ -117,7 +117,7 @@ namespace AMDiS
std::size_t numLocalFe = rowNode.size();
auto contextGeometry = contextGeo.localContext().geometry();
auto contextGeometry = contextGeo.coordInElementContext().geometry();
int quadDeg = std::max({
getQuadratureDegree(contextGeometry,2,coeffOrder(localFctA_),rowNode,colNode),
getQuadratureDegree(contextGeometry,1,coeffOrder(localFctB_),rowNode,colNode),
......@@ -128,7 +128,7 @@ namespace AMDiS
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
auto&& local = contextGeo.local(quad[iq].position());
auto&& local = contextGeo.coordInElement(quad[iq].position());
// The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
......@@ -190,11 +190,11 @@ namespace AMDiS
std::size_t numLocalFe = node.size();
auto const& quad = getQuadratureRule(contextGeo.localContext().geometry(), 0, coeffOrder(localFctF_), node);
auto const& quad = getQuadratureRule(contextGeo.coordInElementContext().geometry(), 0, coeffOrder(localFctF_), node);
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
auto&& local = contextGeo.local(quad[iq].position());
auto&& local = contextGeo.coordInElement(quad[iq].position());
// the values of the shape functions on the reference element at the quadrature point
auto const& shapeValues = node.localBasisValuesAt(local);
......
......@@ -41,7 +41,7 @@ namespace AMDiS
for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
auto&& local = context.local(qp.position());
auto&& local = context.coordInElement(qp.position());
// The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = context.geometry().jacobianInverseTransposed(local);
......
......@@ -39,7 +39,7 @@ namespace AMDiS
for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
auto&& local = contextGeo.local(qp.position());
auto&& local = contextGeo.coordInElement(qp.position());
// The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
......
......@@ -42,7 +42,7 @@ namespace AMDiS
for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
decltype(auto) local = contextGeo.local(qp.position());
decltype(auto) local = contextGeo.coordInElement(qp.position());
// The multiplicative factor in the integral transformation formula
auto dx = contextGeo.integrationElement(qp.position()) * qp.weight();
......
......@@ -44,7 +44,7 @@ namespace AMDiS
for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
auto&& local = contextGeo.local(qp.position());
auto&& local = contextGeo.coordInElement(qp.position());
// The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
......
......@@ -42,7 +42,7 @@ namespace AMDiS
for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
auto&& local = contextGeo.local(qp.position());
auto&& local = contextGeo.coordInElement(qp.position());
// The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
......
......@@ -51,7 +51,7 @@ namespace AMDiS
LocalToGlobalBasisAdapter colLocalBasis(colNode, contextGeo.geometry());
for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
decltype(auto) local = contextGeo.local(qp.position());
decltype(auto) local = contextGeo.coordInElement(qp.position());
// The multiplicative factor in the integral transformation formula
auto factor = contextGeo.integrationElement(qp.position()) * qp.weight();
......
......@@ -45,7 +45,7 @@ namespace AMDiS
for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
decltype(auto) local = contextGeo.local(qp.position());
decltype(auto) local = contextGeo.coordInElement(qp.position());
// The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
......
......@@ -64,7 +64,7 @@ namespace AMDiS
for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
decltype(auto) local = contextGeo.local(qp.position());
decltype(auto) local = contextGeo.coordInElement(qp.position());
// The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
......@@ -115,7 +115,7 @@ namespace AMDiS
for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
auto&& local = contextGeo.local(qp.position());
auto&& local = contextGeo.coordInElement(qp.position());
// The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
......
......@@ -65,7 +65,7 @@ namespace AMDiS
for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
auto&& local = contextGeo.local(qp.position());
auto&& local = contextGeo.coordInElement(qp.position());
// The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
......@@ -106,7 +106,7 @@ namespace AMDiS
for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
auto&& local = contextGeo.local(qp.position());
auto&& local = contextGeo.coordInElement(qp.position());
// The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
......@@ -151,7 +151,7 @@ namespace AMDiS
for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
auto&& local = contextGeo.local(qp.position());
auto&& local = contextGeo.coordInElement(qp.position());
// The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
......
......@@ -48,7 +48,7 @@ namespace AMDiS
for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
decltype(auto) local = contextGeo.local(qp.position());
decltype(auto) local = contextGeo.coordInElement(qp.position());
// The multiplicative factor in the integral transformation formula
auto factor = contextGeo.integrationElement(qp.position()) * qp.weight();
......
......@@ -56,7 +56,7 @@ namespace AMDiS
for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
auto&& local = contextGeo.local(qp.position());
auto&& local = contextGeo.coordInElement(qp.position());
// The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
......
......@@ -36,7 +36,7 @@ namespace AMDiS
for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
auto&& local = contextGeo.local(qp.position());
auto&& local = contextGeo.coordInElement(qp.position());
// The multiplicative factor in the integral transformation formula
auto factor = contextGeo.integrationElement(qp.position()) * qp.weight();
......
......@@ -53,7 +53,7 @@ namespace AMDiS
for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
auto&& local = contextGeo.local(qp.position());
auto&& local = contextGeo.coordInElement(qp.position());
// The multiplicative factor in the integral transformation formula
const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
......@@ -81,7 +81,7 @@ namespace AMDiS
for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
auto&& local = contextGeo.local(qp.position());
auto&& local = contextGeo.coordInElement(qp.position());
// The multiplicative factor in the integral transformation formula
const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
......@@ -114,7 +114,7 @@ namespace AMDiS
for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
auto&& local = contextGeo.local(qp.position());
auto&& local = contextGeo.coordInElement(qp.position());
// The multiplicative factor in the integral transformation formula
const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
......
......@@ -41,7 +41,7 @@ namespace AMDiS
for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
auto&& local = contextGeo.local(qp.position());
auto&& local = contextGeo.coordInElement(qp.position());
// The multiplicative factor in the integral transformation formula
const auto factor = contextGeo.integrationElement(qp.position()) * qp.weight();
......
......@@ -40,7 +40,7 @@ namespace AMDiS
for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
auto&& local = contextGeo.local(qp.position());
auto&& local = contextGeo.coordInElement(qp.position());
// The multiplicative factor in the integral transformation formula
const auto factor = contextGeo.integrationElement(qp.position()) * qp.weight();
......
......@@ -63,7 +63,7 @@ namespace AMDiS
for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
auto&& local = contextGeo.local(qp.position());
auto&& local = contextGeo.coordInElement(qp.position());
// The multiplicative factor in the integral transformation formula
const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
......@@ -99,7 +99,7 @@ namespace AMDiS
for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
auto&& local = contextGeo.local(qp.position());
auto&& local = contextGeo.coordInElement(qp.position());
// The multiplicative factor in the integral transformation formula
const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
......@@ -144,7 +144,7 @@ namespace AMDiS
for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
auto&& local = contextGeo.local(qp.position());
auto&& local = contextGeo.coordInElement(qp.position());
// The multiplicative factor in the integral transformation formula
const auto factor = contextGeo.integrationElement(qp.position()) * qp.weight();
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment