extrusion.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
//
// SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_DISCRETIZATION_EXTRUSION_HH
#define DUMUX_DISCRETIZATION_EXTRUSION_HH
 
#include <concepts>
#include <type_traits>
 
#include <dune/common/std/type_traits.hh>
 
namespace Dumux {
 
namespace Detail {

template<class T, class FVGeo, int codim>
concept EntityGeometry = std::same_as<
    std::decay_t<T>,
    typename FVGeo::GridGeometry::GridView::template Codim<codim>::Entity::Geometry
>;
 
} // end namespace Detail

struct NoExtrusion
{
    template<class FVGeo, class Face>
        requires (!Detail::EntityGeometry<Face, FVGeo, 1>)
    static constexpr auto area(const FVGeo&, const Face& face)
    { return face.area(); }
 
    template<class FVGeo>
    static constexpr auto area(const FVGeo&, const typename FVGeo::GridGeometry::GridView::template Codim<1>::Entity::Geometry& geo)
    { return geo.volume(); }
 
    template<class FVGeo, class SCV>
        requires (!Detail::EntityGeometry<SCV, FVGeo, 0>)
    static constexpr auto volume(const FVGeo&, const SCV& scv)
    { return scv.volume(); }
 
    template<class FVGeo>
    static constexpr auto volume(const FVGeo&, const typename FVGeo::GridGeometry::GridView::template Codim<0>::Entity::Geometry& geo)
    { return geo.volume(); }
 
    template<class Geometry>
    static constexpr auto integrationElement(const Geometry& geo, const typename Geometry::LocalCoordinate& x)
    { return geo.integrationElement(x); }
};

template<int radAx = 0>
struct RotationalExtrusion
{
    static constexpr int radialAxis = radAx;

    template<class FVGeo, class Face>
        requires (!Detail::EntityGeometry<Face, FVGeo, 1>)
    static constexpr auto area(const FVGeo&, const Face& face)
    {
        static_assert(int(Face::Traits::Geometry::mydimension) == int(Face::Traits::Geometry::coorddimension-1), "Area element to be called with a codim-1-entity!");
        static_assert(Face::Traits::Geometry::coorddimension <= 2, "Axis rotation only makes sense for geometries up to 2D!");
        static_assert(radialAxis < int(Face::Traits::Geometry::coorddimension), "Illegal radial axis!");
 
        // Guldinus theorem
        return face.area()*2.0*M_PI*face.center()[radialAxis];
    }

    template<class FVGeo>
    static constexpr auto area(const FVGeo&, const typename FVGeo::GridGeometry::GridView::template Codim<1>::Entity::Geometry& geo)
    {
        using Geometry = typename FVGeo::GridGeometry::GridView::template Codim<1>::Entity::Geometry;
        static_assert(Geometry::coorddimension <= 2, "Axis rotation only makes sense for geometries up to 2D!");
        static_assert(radialAxis < int(Geometry::coorddimension), "Illegal radial axis!");
 
        // Guldinus theorem
        return geo.volume()*2.0*M_PI*geo.center()[radialAxis];
    }

    template<class FVGeo, class SCV>
        requires (!Detail::EntityGeometry<SCV, FVGeo, 0>)
    static constexpr auto volume(const FVGeo&, const SCV& scv)
    {
        static_assert(int(SCV::Traits::Geometry::mydimension) == int(SCV::Traits::Geometry::coorddimension), "Volume element to be called with a codim-0-entity!");
        static_assert(SCV::Traits::Geometry::coorddimension <= 2, "Axis rotation only makes sense for geometries up to 2D!");
        static_assert(radialAxis < int(SCV::Traits::Geometry::coorddimension), "Illegal radial axis!");
 
        // Guldinus theorem
        return scv.volume()*2.0*M_PI*scv.center()[radialAxis];
    }

    template<class FVGeo>
    static constexpr auto volume(const FVGeo&, const typename FVGeo::GridGeometry::GridView::template Codim<0>::Entity::Geometry& geo)
    {
        using Geometry = typename FVGeo::GridGeometry::GridView::template Codim<0>::Entity::Geometry;
        static_assert(Geometry::coorddimension <= 2, "Axis rotation only makes sense for geometries up to 2D!");
        static_assert(radialAxis < int(Geometry::coorddimension), "Illegal radial axis!");
 
        // Guldinus theorem
        return geo.volume()*2.0*M_PI*geo.center()[radialAxis];
    }

    template<class Geometry>
    static constexpr auto integrationElement(const Geometry& geo, const typename Geometry::LocalCoordinate& x)
    {
        static_assert(Geometry::coorddimension <= 2, "Axis rotation only makes sense for geometries up to 2D!");
        static_assert(radialAxis < int(Geometry::coorddimension), "Illegal radial axis!");
 
        // Multiply with the polar extrusion factor (2*pi) and the determinant of the transformation Jacobian (radius)
        return geo.integrationElement(x)*2.0*M_PI*geo.global(x)[radialAxis];
    }
};

struct SphericalExtrusion
{
    template<class FVGeo, class Face>
        requires (!Detail::EntityGeometry<Face, FVGeo, 1>)
    static constexpr auto area(const FVGeo&, const Face& face)
    {
        static_assert(int(Face::Traits::Geometry::mydimension) == int(Face::Traits::Geometry::coorddimension-1), "Area element to be called with a codim-1-entity!");
        static_assert(Face::Traits::Geometry::coorddimension == 1, "Spherical rotation only makes sense for 1D geometries!");
 
        // sphere surface area
        const auto radius = face.center()[0];
        return 4.0*M_PI*radius*radius;
    }

    template<class FVGeo>
    static constexpr auto area(const FVGeo&, const typename FVGeo::GridGeometry::GridView::template Codim<1>::Entity::Geometry& geo)
    {
        using Geometry = typename FVGeo::GridGeometry::GridView::template Codim<1>::Entity::Geometry;
        static_assert(Geometry::coorddimension == 1, "Spherical rotation only makes sense for 1D geometries!");
 
        // sphere surface area
        const auto radius = geo.center()[0];
        return 4.0*M_PI*radius*radius;
    }

    template<class FVGeo, class SCV>
        requires (!Detail::EntityGeometry<SCV, FVGeo, 0>)
    static constexpr auto volume(const FVGeo& fvGeometry, const SCV& scv)
    {
        static_assert(int(SCV::Traits::Geometry::mydimension) == int(SCV::Traits::Geometry::coorddimension), "Volume element to be called with a codim-0-entity!");
        static_assert(SCV::Traits::Geometry::coorddimension == 1, "Spherical rotation only makes sense for 1D geometries!");
 
        // subtract two balls
        const auto geo = fvGeometry.geometry(scv);
        const auto radius0 = geo.corner(0)[0];
        const auto radius1 = geo.corner(1)[0];
        using std::abs;
        return 4.0/3.0*M_PI*abs(radius1*radius1*radius1 - radius0*radius0*radius0);
    }

    template<class FVGeo>
    static constexpr auto volume(const FVGeo&, const typename FVGeo::GridGeometry::GridView::template Codim<0>::Entity::Geometry& geo)
    {
        using Geometry = typename FVGeo::GridGeometry::GridView::template Codim<0>::Entity::Geometry;
        static_assert(Geometry::coorddimension == 1, "Spherical rotation only makes sense for 1D geometries!");
 
        // subtract two balls
        const auto radius0 = geo.corner(0)[0];
        const auto radius1 = geo.corner(1)[0];
        using std::abs;
        return 4.0/3.0*M_PI*abs(radius1*radius1*radius1 - radius0*radius0*radius0);
    }

    template<class Geometry>
    static constexpr auto integrationElement(const Geometry& geo, const typename Geometry::LocalCoordinate& x)
    {
        static_assert(Geometry::coorddimension == 1, "Spherical rotation only makes sense for 1D geometries!");
 
        // Multiply with the constant spherical extrusion factor (int_0^2pi int_0^pi sin(phi) dphi dtheta = 4pi)
        // and the remaining (radius-dependent) part of determinant of the transformation Jacobian (radius*radius)
        const auto radius = geo.global(x)[0];
        return geo.integrationElement(x)*4.0*M_PI*radius*radius;
    }
};

template<class T>
class Extrusion
{
    template<class G>
    using E = typename G::Extrusion;
public:
    using type = typename Dune::Std::detected_or<NoExtrusion, E, T>::type;
};

template<class T>
using Extrusion_t = typename Extrusion<T>::type;

template<class T>
inline constexpr bool isRotationalExtrusion = false;

template<int radialAxis>
inline constexpr bool isRotationalExtrusion<RotationalExtrusion<radialAxis>> = true;
 
} // end namespace Dumux
 
#endif