diff --git a/dumux/common/pointsource.hh b/dumux/common/pointsource.hh index 2ca28690edaf53502abd10fbb1fbd0abb7f84ffe..eaf3b3f6bd901c2b07a52964bdb3fcf9b63ed353 100644 --- a/dumux/common/pointsource.hh +++ b/dumux/common/pointsource.hh @@ -316,7 +316,7 @@ public: constexpr int dim = GridGeometry::GridView::dimension; Dune::ReservedVector<std::size_t, 1<<dim> scvIndices; for (const auto& scv : scvs(fvGeometry)) - if (intersectsPointGeometry(globalPos, fvGeometry.geometry(scv))) + if (intersectsPointGeometry(globalPos, scv.geometry())) scvIndices.push_back(scv.indexInElement()); // for all scvs that tested positive add the point sources diff --git a/dumux/discretization/box/fvelementgeometry.hh b/dumux/discretization/box/fvelementgeometry.hh index aeaf3c984e01a2db731d3d4370334e4b355ee3a9..4424e4c64127ffc9c74c60327cd572a8ff50bacc 100644 --- a/dumux/discretization/box/fvelementgeometry.hh +++ b/dumux/discretization/box/fvelementgeometry.hh @@ -28,10 +28,6 @@ #include <optional> #include <utility> -#include <unordered_map> -#include <array> -#include <vector> - #include <dune/geometry/type.hh> #include <dune/localfunctions/lagrange/pqkfactory.hh> @@ -63,10 +59,6 @@ class BoxFVElementGeometry<GG, true> using LocalIndexType = typename IndexTraits<GridView>::LocalIndex; using CoordScalar = typename GridView::ctype; using FeLocalBasis = typename GG::FeCache::FiniteElementType::Traits::LocalBasisType; - - using GeometryHelper = BoxGeometryHelper<GridView, dim, - typename GG::SubControlVolume, - typename GG::SubControlVolumeFace>; public: //! export the element type using Element = typename GridView::template Codim<0>::Entity; @@ -193,20 +185,6 @@ public: bool hasBoundaryScvf() const { return gridGeometry().hasBoundaryScvf(eIdx_); } - //! Geometry of a sub control volume - typename SubControlVolume::Traits::Geometry geometry(const SubControlVolume& scv) const - { - assert(isBound()); - return gridGeometryPtr_->geometry(element(), scv); - } - - //! Geometry of a sub control volume face - typename SubControlVolumeFace::Traits::Geometry geometry(const SubControlVolumeFace& scvf) const - { - assert(isBound()); - return gridGeometryPtr_->geometry(element(), scvf); - } - private: const GridGeometry* gridGeometryPtr_; GridIndexType eIdx_; @@ -353,32 +331,7 @@ public: bool hasBoundaryScvf() const { return hasBoundaryScvf_; } - //! Geometry of a sub control volume - typename SubControlVolume::Traits::Geometry geometry(const SubControlVolume& scv) const - { - assert(isBound()); - const auto geo = element().geometry(); - return { Dune::GeometryTypes::cube(dim), GeometryHelper(geo).getScvCorners(scv.indexInElement()) }; - } - - //! Geometry of a sub control volume face - typename SubControlVolumeFace::Traits::Geometry geometry(const SubControlVolumeFace& scvf) const - { - assert(isBound()); - const auto geo = element().geometry(); - if (scvf.boundary()) - { - const auto localBoundaryIndex = scvf.index() - numInnerScvf_(element()); - const auto& key = scvfBoundaryGeometryKeys_[localBoundaryIndex]; - return { Dune::GeometryTypes::cube(dim-1), GeometryHelper(geo).getBoundaryScvfCorners(key[0], key[1]) }; - } - else - return { Dune::GeometryTypes::cube(dim-1), GeometryHelper(geo).getScvfCorners(scvf.index()) }; - } - private: - unsigned int numInnerScvf_(const Element& element) const - { return (dim==1) ? 1 : element.subEntities(dim-1); } void makeElementGeometries_() { @@ -396,7 +349,7 @@ private: scvs_.resize(elementGeometry.corners()); for (LocalIndexType scvLocalIdx = 0; scvLocalIdx < elementGeometry.corners(); ++scvLocalIdx) { - // get associated dof index + // get asssociated dof index const auto dofIdxGlobal = gridGeometry().vertexMapper().subIndex(element, scvLocalIdx, dim); // add scv to the local container @@ -407,9 +360,8 @@ private: } // construct the sub control volume faces - const auto numInnerScvf = numInnerScvf_(element); + const auto numInnerScvf = (dim==1) ? 1 : element.subEntities(dim-1); scvfs_.resize(numInnerScvf); - scvfBoundaryGeometryKeys_.clear(); LocalIndexType scvfLocalIdx = 0; for (; scvfLocalIdx < numInnerScvf; ++scvfLocalIdx) @@ -448,11 +400,6 @@ private: std::move(localScvIndices), true); - scvfBoundaryGeometryKeys_.emplace_back(std::array<LocalIndexType, 2>{{ - static_cast<LocalIndexType>(intersection.indexInInside()), - static_cast<LocalIndexType>(isScvfLocalIdx) - }}); - // increment local counter scvfLocalIdx++; } @@ -470,7 +417,6 @@ private: //! vectors to store the geometries locally after binding an element std::vector<SubControlVolume> scvs_; std::vector<SubControlVolumeFace> scvfs_; - std::vector<std::array<LocalIndexType, 2>> scvfBoundaryGeometryKeys_; bool hasBoundaryScvf_ = false; }; diff --git a/dumux/discretization/box/fvgridgeometry.hh b/dumux/discretization/box/fvgridgeometry.hh index 3e74e7bf25084382c1fac22c92866584cdd8cf26..b26b54d3c28b628e491361870960b93b90776968 100644 --- a/dumux/discretization/box/fvgridgeometry.hh +++ b/dumux/discretization/box/fvgridgeometry.hh @@ -28,8 +28,6 @@ #include <utility> #include <unordered_map> -#include <array> -#include <vector> #include <dune/localfunctions/lagrange/lagrangelfecache.hh> @@ -201,36 +199,11 @@ public: bool hasBoundaryScvf(GridIndexType eIdx) const { return hasBoundaryScvf_[eIdx]; } - //! Geometry of a sub control volume - typename SubControlVolume::Traits::Geometry geometry(const Element& element, const SubControlVolume& scv) const - { - assert(isBound()); - const auto geo = element.geometry(); - return { Dune::GeometryTypes::cube(dim), GeometryHelper(geo).getScvCorners(scv.indexInElement()) }; - } - - //! Geometry of a sub control volume face - typename SubControlVolumeFace::Traits::Geometry geometry(const Element& element, const SubControlVolumeFace& scvf) const - { - assert(isBound()); - const auto eIdx = this->elementMapper().index(element); - const auto geo = element.geometry(); - if (scvf.boundary()) - { - const auto localBoundaryIndex = scvf.index() - numInnerScvf_(element); - const auto& key = scvfBoundaryGeometryKeys_.at(eIdx)[localBoundaryIndex]; - return { Dune::GeometryTypes::cube(dim-1), GeometryHelper(geo).getBoundaryScvfCorners(key[0], key[1]) }; - } - else - return { Dune::GeometryTypes::cube(dim-1), GeometryHelper(geo).getScvfCorners(scvf.index()) }; - } - private: void update_() { scvs_.clear(); scvfs_.clear(); - scvfBoundaryGeometryKeys_.clear(); auto numElements = this->gridView().size(0); scvs_.resize(numElements); @@ -314,11 +287,6 @@ private: std::move(localScvIndices), true); - scvfBoundaryGeometryKeys_[eIdx].emplace_back(std::array<LocalIndexType, 2>{{ - static_cast<LocalIndexType>(intersection.indexInInside()), - static_cast<LocalIndexType>(isScvfLocalIdx) - }}); - // increment local counter scvfLocalIdx++; } @@ -379,15 +347,10 @@ private: DUNE_THROW(Dune::NotImplemented, "Periodic boundaries for box method for parallel simulations!"); } - unsigned int numInnerScvf_(const Element& element) const - { return (dim==1) ? 1 : element.subEntities(dim-1); } - const FeCache feCache_; std::vector<std::vector<SubControlVolume>> scvs_; std::vector<std::vector<SubControlVolumeFace>> scvfs_; - std::unordered_map<GridIndexType, std::vector<std::array<LocalIndexType, 2>>> scvfBoundaryGeometryKeys_; - // TODO do we need those? std::size_t numScv_; std::size_t numScvf_; diff --git a/test/discretization/box/test_boxfvgeometry.cc b/test/discretization/box/test_boxfvgeometry.cc index 7550d624eda1eb8de9837e3776c915b3f43d47d6..98ae68f315d698a7415b334424143bbe38a233fa 100644 --- a/test/discretization/box/test_boxfvgeometry.cc +++ b/test/discretization/box/test_boxfvgeometry.cc @@ -100,20 +100,6 @@ int main (int argc, char *argv[]) for (auto&& scv : scvs(fvGeometry)) { std::cout << "-- scv " << scv.localDofIndex() << " center at: " << scv.center() << " , volume: " << scv.volume() << std::endl; - - const auto scvGeo = fvGeometry.geometry(scv); - - if ((scv.center() - scvGeo.center()).two_norm2() > 1e-15) - DUNE_THROW(Dune::Exception, - "scv.center() and fvGeometry.geometry(scv).center() do not match! " - << scv.center() << " :: " << scvGeo.center() - ); - - if (std::abs(scv.volume() - scvGeo.volume()) > 1e-7) - DUNE_THROW(Dune::Exception, - "scv.volume() and fvGeometry.geometry(scv).volume() do not match! " - << scv.volume() << " :: " << scvGeo.volume() - ); } auto range2 = scvfs(fvGeometry); @@ -125,21 +111,6 @@ int main (int argc, char *argv[]) for (auto&& scvf : scvfs(fvGeometry)) { std::cout << "-- scvf " << scvf.index() << " ip at: " << scvf.ipGlobal() << " normal: " << scvf.unitOuterNormal(); - - const auto scvfGeo = fvGeometry.geometry(scvf); - - if ((scvf.ipGlobal() - scvfGeo.center()).two_norm2() > 1e-15) - DUNE_THROW(Dune::Exception, - "scvf.ipGlobal() and fvGeometry.geometry(scvf).center() do not match! " - << scvf.ipGlobal() << " :: " << scvfGeo.center() - ); - - if (std::abs(scvf.area() - scvfGeo.volume()) > 1e-7) - DUNE_THROW(Dune::Exception, - "scvf.area() and fvGeometry.geometry(scvf).volume() do not match! " - << scvf.area() << " :: " << scvfGeo.volume() - ); - if (scvf.boundary()) { ++boundaryCount;