Skip to content
Snippets Groups Projects
Commit 8ce8b3c0 authored by Jonathan Youett's avatar Jonathan Youett
Browse files

Cleanup

parent 17321126
No related branches found
No related tags found
No related merge requests found
......@@ -28,12 +28,10 @@ template<int mydim, int coorddim, class GridImp> class SubGridLocalGeometry;
because the Dune grid interface specifies the number and types of template
parameters of that class.
*/
template<int mydim, int coorddim, class GridImp> //, bool isLocal>
template<int mydim, int coorddim, class GridImp>
class SubGridGeometry
{
typedef typename GridImp::ctype ctype;
using ctype = typename GridImp::ctype;
public:
......@@ -41,12 +39,6 @@ class SubGridGeometry
enum {CodimInHostGrid = GridImp::HostGridType::dimension - mydim};
enum {DimensionWorld = GridImp::HostGridType::dimensionworld};
// Select the type we are wrapping depending on the isLocal flag
// - isLocal == true: LocalGeometry
// - isLocal == false: Geometry
// typedef typename std::conditional<isLocal,
//typename GridImp::HostGridType::Traits::template Codim<CodimInHostGrid>::LocalGeometry,
// typename GridImp::HostGridType::Traits::template Codim<CodimInHostGrid>::Geometry>::type HostGridGeometry;
using HostGridGeometry = typename GridImp::HostGridType::Traits::template Codim<CodimInHostGrid>::Geometry;
//! The type of the Jacobian matrix of the mapping from the reference element to this element
......@@ -126,24 +118,21 @@ class SubGridGeometry
return hostGeometry_.center();
}
//! The Jacobian matrix of the mapping from the reference element to this element
const JacobianTransposed jacobianTransposed (const FieldVector<ctype, mydim>& local) const {
return hostGeometry_.jacobianTransposed(local);
}
//! The inverse of the Jacobian matrix of the mapping from the reference element to this element
const JacobianInverseTransposed jacobianInverseTransposed (const FieldVector<ctype, mydim>& local) const {
return hostGeometry_.jacobianInverseTransposed(local);
}
private:
const HostGridGeometry hostGeometry_;
};
/** \brief Encapsulate a host grid LocalGeometry */
/** \brief Local SubGrid geometries. They can be either a wrapped hostgrid geometry or a multilineargeometry. */
template<int mydim, int coorddim, class GridImp>
class SubGridLocalGeometry
{
......@@ -160,7 +149,8 @@ class SubGridLocalGeometry
// TODO does this work for arbitratry co-dimensions?
using HostGridLocalGeometry = typename GridImp::HostGridType::Traits::template Codim<CodimInHostGrid>::LocalGeometry;
using MultiLinGeometry = MultiLinearGeometry<ctype, mydim, coorddim>;
// variant can be either of the template types
// use a variant to store either a HostGridLocalGeometry or a MultiLinearGeometry
using GeometryContainer = Std::variant<HostGridLocalGeometry, MultiLinGeometry>;
public:
......@@ -178,7 +168,7 @@ class SubGridLocalGeometry
using JacobianInverseTransposed = typename std::common_type<Dune::FieldMatrix<ctype,coorddim, mydim>,
typename HostGridLocalGeometry::JacobianInverseTransposed>::type;
/** Constructor for a given host grid geometry*/
/** \brief Construct from a given host grid geometry*/
SubGridLocalGeometry(HostGridLocalGeometry hostLocalGeometry)
: localGeometry_(std::move(hostLocalGeometry))
{}
......@@ -194,7 +184,6 @@ class SubGridLocalGeometry
return Std::visit([&](auto&& geom) {return geom.type();}, localGeometry_);
}
/** \brief Return true if the geometry mapping is affine and false otherwise */
bool affine() const {
return Std::visit([&](auto&& geom) {return geom.affine();}, localGeometry_);
......@@ -205,27 +194,23 @@ class SubGridLocalGeometry
return Std::visit([&](auto&& geom) {return geom.corners();}, localGeometry_);
}
//! access to coordinates of corners. Index is the number of the corner
FieldVector<ctype, coorddim> corner (int i) const {
return Std::visit([&](auto&& geom) {return geom.corner(i);}, localGeometry_);
}
/** \brief Maps a local coordinate within reference element to
* global coordinate in element */
FieldVector<ctype, coorddim> global (const FieldVector<ctype, mydim>& local) const{
return Std::visit([&](auto&& geom) {return geom.global(local);}, localGeometry_);
}
/** \brief Maps a global coordinate within the element to a
* local coordinate in its reference element */
FieldVector<ctype, mydim> local (const FieldVector<ctype, coorddim>& global) const {
return Std::visit([&](auto&& geom) {return geom.local(global);}, localGeometry_);
}
/**
*/
ctype integrationElement (const FieldVector<ctype, mydim>& local) const {
......@@ -250,19 +235,17 @@ class SubGridLocalGeometry
return Std::visit([&](auto&& geom) {return geom.center();}, localGeometry_);
}
//! The Jacobian matrix of the mapping from the reference element to this element
const JacobianTransposed jacobianTransposed (const FieldVector<ctype, mydim>& local) const {
return Std::visit([&](auto&& geom) {return geom.jacobianTransposed(local);}, localGeometry_);
}
//! The inverse of the Jacobian matrix of the mapping from the reference element to this element
const JacobianInverseTransposed jacobianInverseTransposed (const FieldVector<ctype, mydim>& local) const {
return Std::visit([&](auto&& geom) {return geom.jacobianInverseTransposed(local);}, localGeometry_);
}
private:
const GeometryContainer localGeometry_;
};
......
......@@ -92,7 +92,6 @@ class SubGridLeafIntersection :
using ctype = typename GridImp::ctype;
using Base::insideIntersect_;
using ctype = typename GridImp::ctype;
public:
typedef typename GridImp::template Codim<0>::Entity Entity;
......@@ -233,38 +232,39 @@ public:
return boundary() or ((outside_.level() == inside_.level()) and insideIntersect_.conforming());
}
//! intersection of codimension 1 of this neighbor with element where
//! iteration started.
//! Here returned element is in LOCAL coordinates of the element
//! where iteration started.
/** \brief Return the local geometry, mapping local coordinates of the intersection to
* local coordinates of the inside element.
* If the intersection is not part of the host grid, construct one using MultiLinearGeometry.
*/
LocalGeometry geometryInInside () const {
if (boundary())
return LocalGeometry(LocalGeometryImpl(insideIntersect_.geometryInInside()));
else {
// Case 0: inside_'s the smaller element, or both elements are on the same level.
if (inside_.level() >= outside_.level()) {
if (inside_.level() >= outside_.level())
return LocalGeometry(LocalGeometryImpl(insideIntersect_.geometryInInside()));
}
// Case 1: outsideIntersect_ is the intersection we want, embedded in inside_'s geometry
auto outsideGeo = outsideIntersect_.geometry();
auto outsideIsGeo = outsideIntersect_.geometry();
auto insideGeo = inside_.geometry();
using Coords = std::decay_t<decltype(outsideGeo.corner(0))>;
auto corners = std::vector<Coords>(outsideGeo.corners());
using Coords = std::decay_t<decltype(insideGeo.local(outsideIsGeo.corner(0)))>;
auto corners = std::vector<Coords>(outsideIsGeo.corners());
// Map the coordinates of the corners from the outside's local coordinates into inside's
// local coordinates
for (size_t i = 0; i < corners.size(); i++)
corners[i] = inside_.geometry().local(outsideGeo.corner(i));
corners[i] = insideGeo.local(outsideIsGeo.corner(i));
// create a new geometry with the coordinates we just computed
auto geo = Dune::MultiLinearGeometry<ctype, dim-1, dim>(outsideGeo.type(), corners);
auto geo = Dune::MultiLinearGeometry<ctype, dim-1, dim>(outsideIsGeo.type(), corners);
return LocalGeometry(LocalGeometryImpl(std::move(geo)));
}
}
//! intersection of codimension 1 of this neighbor with element where iteration started.
//! Here returned element is in LOCAL coordinates of neighbor
/** \brief Return the local geometry, mapping local coordinates of the intersection to
* local coordinates of the outside element.
*/
LocalGeometry geometryInOutside () const {
if(boundary())
DUNE_THROW(Dune::Exception, "No outside geometry available for boundary intersections");
......@@ -272,32 +272,33 @@ public:
// Case 0: outside_'s the smaller element, or both elements are on the same level.
if (inside_.level() < outside_.level()) {
auto outsideGeo = outsideIntersect_.geometry();
auto outsideIsGeo = outsideIntersect_.geometry();
auto outsideGeo = outside_.geometry();
using Coords = std::decay_t<decltype(outside_.geometry().local(outsideGeo.corner(0)))>;
auto corners = std::vector<Coords>(outsideGeo.corners());
using Coords = std::decay_t<decltype(outsideGeo.local(outsideIsGeo.corner(0)))>;
auto corners = std::vector<Coords>(outsideIsGeo.corners());
for (size_t i = 0; i < corners.size(); i++)
corners[i] = outside_.geometry().local(outsideGeo.corner(i));
corners[i] = outsideGeo.local(outsideIsGeo.corner(i));
// create a new geometry with the coordinates we just computed
auto geo = Dune::MultiLinearGeometry<ctype, dim-1, dim>(outsideGeo.type(), corners);
auto geo = Dune::MultiLinearGeometry<ctype, dim-1, dim>(outsideIsGeo.type(), corners);
return LocalGeometry(LocalGeometryImpl(geo));
}
// Case 1: outsideIntersect_ is the intersection we want, embedded in inside_'s geometry
auto insideIsGeo = insideIntersect_.geometry();
auto outsideGeo = outside_.geometry();
auto insideGeo = insideIntersect_.geometry();
using Coords = std::decay_t<decltype(outside_.geometry().local(insideGeo.corner(0)))>;
auto corners = std::vector<Coords>(insideGeo.corners());
using Coords = std::decay_t<decltype(outsideGeo.local(insideIsGeo.corner(0)))>;
auto corners = std::vector<Coords>(insideIsGeo.corners());
// Map the coordinates of the corners from the insides's local coordinates into outside's
// local coordinates
for (size_t i = 0; i < corners.size(); i++)
corners[i] = outside_.geometry().local(insideGeo.corner(i));
corners[i] = outsideGeo.local(insideIsGeo.corner(i));
// create a new geometry with the coordinates we just computed
auto geo = Dune::MultiLinearGeometry<ctype, dim-1, dim>(insideGeo.type(), corners);
auto geo = Dune::MultiLinearGeometry<ctype, dim-1, dim>(insideIsGeo.type(), corners);
return LocalGeometry(LocalGeometryImpl(std::move(geo)));
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment