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:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
 *   the Free Software Foundation, either version 3 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
#ifndef DUMUX_DISCRETIZATION_EXTRUSION_HH
#define DUMUX_DISCRETIZATION_EXTRUSION_HH
 
#include <dune/common/std/type_traits.hh>
#include <dumux/common/deprecated.hh>
 
namespace Dumux {
 
struct NoExtrusion
{
    template<class SCVF>
    [[deprecated("Will be removed after 3.6. Use area(fvGeometry, scvf)")]]
    static constexpr auto area(const SCVF& scvf)
    { return scvf.area(); }
 
    template<class SCV>
    [[deprecated("Will be removed after 3.6. Use volume(fvGeometry, scv)")]]
    static constexpr auto volume(const SCV& scv)
    { return scv.volume(); }
 
    template<class FVGeo, class SCVF>
    static constexpr auto area(const FVGeo&, const SCVF& scvf)
    { return scvf.area(); }
 
    template<class FVGeo, class SCV>
    static constexpr auto volume(const FVGeo&, const SCV& scv)
    { return scv.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 SCVF>
    [[deprecated("Will be removed after 3.6. Use area(fvGeometry, scvf)")]]
    static constexpr auto area(const SCVF& scvf)
    {
        static_assert(int(SCVF::Traits::Geometry::mydimension) == int(SCVF::Traits::Geometry::coorddimension-1), "Area element to be called with a codim-1-entity!");
        static_assert(SCVF::Traits::Geometry::coorddimension <= 2, "Axis rotation only makes sense for geometries up to 2D!");
        static_assert(radialAxis < int(SCVF::Traits::Geometry::coorddimension), "Illegal radial axis!");
 
        // Guldinus theorem
        return scvf.area()*2.0*M_PI*scvf.center()[radialAxis];
    }
 
    template<class FVGeo, class SCVF>
    static constexpr auto area(const FVGeo&, const SCVF& scvf)
    {
        static_assert(int(SCVF::Traits::Geometry::mydimension) == int(SCVF::Traits::Geometry::coorddimension-1), "Area element to be called with a codim-1-entity!");
        static_assert(SCVF::Traits::Geometry::coorddimension <= 2, "Axis rotation only makes sense for geometries up to 2D!");
        static_assert(radialAxis < int(SCVF::Traits::Geometry::coorddimension), "Illegal radial axis!");
 
        // Guldinus theorem
        return scvf.area()*2.0*M_PI*scvf.center()[radialAxis];
    }
 
    template<class SCV>
    [[deprecated("Will be removed after 3.6. Use volume(fvGeometry, scv)")]]
    static constexpr auto volume(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, class SCV>
    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 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 SCVF>
    [[deprecated("Will be removed after 3.6. Use area(fvGeometry, scvf)")]]
    static constexpr auto area(const SCVF& scvf)
    {
        static_assert(int(SCVF::Traits::Geometry::mydimension) == int(SCVF::Traits::Geometry::coorddimension-1), "Area element to be called with a codim-1-entity!");
        static_assert(SCVF::Traits::Geometry::coorddimension == 1, "Spherical rotation only makes sense for 1D geometries!");
 
        // sphere surface area
        const auto radius = scvf.center()[0];
        return 4.0*M_PI*radius*radius;
    }
 
    template<class FVGeo, class SCVF>
    static constexpr auto area(const FVGeo&, const SCVF& scvf)
    {
        static_assert(int(SCVF::Traits::Geometry::mydimension) == int(SCVF::Traits::Geometry::coorddimension-1), "Area element to be called with a codim-1-entity!");
        static_assert(SCVF::Traits::Geometry::coorddimension == 1, "Spherical rotation only makes sense for 1D geometries!");
 
        // sphere surface area
        const auto radius = scvf.center()[0];
        return 4.0*M_PI*radius*radius;
    }
 
    template<class SCV>
    [[deprecated("Will be removed after 3.6. Use volume(fvGeometry, scv)")]]
    static constexpr auto volume(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 radius0 = scv.corner(0)[0];
        const auto radius1 = scv.corner(1)[0];
        using std::abs;
        return 4.0/3.0*M_PI*abs(radius1*radius1*radius1 - radius0*radius0*radius0);
    }
 
    template<class FVGeo, class SCV>
    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
        if constexpr (Deprecated::hasSCVGeometryInterface<FVGeo>())
        {
            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);
        }
        else
        {
            // fall-back version for old interface (may trigger deprecation warning)
            const auto geo = scv.geometry();
            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