refinedsimplexgeometry.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_FUFEM_GEOMETRY_REFINEDSIMPLEXGEOMETRY_HH
#define DUNE_FUFEM_GEOMETRY_REFINEDSIMPLEXGEOMETRY_HH
 
#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
 
#include <dune/geometry/type.hh>
#include <dune/geometry/affinegeometry.hh>
 
template <class ct, int localdim, int worlddim>
class RefinedSimplexGeometry {
    public:
    typedef ct ctype;
    typedef Dune::FieldVector<ctype,worlddim> GlobalCoordinate;
 
    template <class VertexContainerType>
    RefinedSimplexGeometry(const VertexContainerType& vertices):
        vertices_(vertices)
    {
        DUNE_THROW(Dune::NotImplemented, "RefinedSimplexGeometry is implemented only for dim=2.");
    }
 
    private:
    const std::array<GlobalCoordinate,2*(worlddim+1)> vertices_;
};
 
template <class ct>
class RefinedSimplexGeometry<ct,2,3> {
    public:
 
    typedef ct ctype;
    typedef Dune::FieldVector<ctype,2> LocalCoordinate;
    typedef Dune::FieldVector<ctype,3> GlobalCoordinate;
 
    typedef Dune::FieldMatrix<ctype,2,3> JacobianTransposed;
    class JacobianInverseTransposed;
 
    template <class VertexContainerType>
    RefinedSimplexGeometry(const VertexContainerType& vertices):
        vertices_(vertices)
    {
        indexMap_[0][0] = 0;
        indexMap_[0][1] = 3;
        indexMap_[0][2] = 4;
 
        indexMap_[1][0] = 3;
        indexMap_[1][1] = 1;
        indexMap_[1][2] = 5;
 
        indexMap_[2][0] = 4;
        indexMap_[2][1] = 5;
        indexMap_[2][2] = 2;
 
        indexMap_[3][0] = 5;
        indexMap_[3][1] = 4;
        indexMap_[3][2] = 3;
 
        for (std::size_t i = 0; i<4; ++i)
        {
            translation_[i] = vertices_[indexMap_[i][0]];
            jacobianTransposed_[i][0] = linearMapSubElementLocalToWorldTransposed_[i][0] = vertices_[indexMap_[i][1]] - translation_[i];
            jacobianTransposed_[i][1] = linearMapSubElementLocalToWorldTransposed_[i][1] = vertices_[indexMap_[i][2]] - translation_[i];
 
            jacobianTransposed_[i] *= std::pow(-1,int(i==3))*2;
 
            jacobianInverseTransposed_[i].setup(jacobianTransposed_[i]);
 
            integrationElement_[i] = jacobianInverseTransposed_[i].detInv();
        }
    }
 
    GlobalCoordinate global(const LocalCoordinate& local) const
    {
        LocalCoordinate subElementLocal;
        int subElement=-1;
 
        getSubElement(local,subElement,subElementLocal);
 
        //std::cout << "local: " << local << ", subElement= " << subElement << ", subElementlocal: " << subElementLocal << std::endl;
 
        GlobalCoordinate returnValue = translation_[subElement];
        //std::cout << "translation: " << returnValue << ", jacobianTransposed: " << jacobianTransposed_[subElement] << std::endl;
        linearMapSubElementLocalToWorldTransposed_[subElement].umtv(subElementLocal,returnValue);
 
        //std::cout << "global= " << returnValue << std::endl;
 
        return returnValue;
    }
 
    bool affine() const
    {
        return false;
    }
 
    Dune::GeometryType type() const
    {
        DUNE_THROW(Dune::NotImplemented,"Method type() not implemented for RefinedSimplexGeometry");
    }
 
    int corners() const
    {
        return 6;
    }
 
    GlobalCoordinate corner(int i) const
    {
        return vertices_[i];
    }
 
    LocalCoordinate local(const GlobalCoordinate&) const
    {
        DUNE_THROW(Dune::NotImplemented,"Method local() not implemented for RefinedSimplexGeometry<ctype,2,3>");
    }
 
    GlobalCoordinate center() const
    {
        DUNE_THROW(Dune::NotImplemented,"Method center() not implemented for RefinedSimplexGeometry<ctype,2,3>");
    }
 
    ctype volume() const
    {
        DUNE_THROW(Dune::NotImplemented,"Method volume() not implemented for RefinedSimplexGeometry<ctype,2,3>");
    }
 
    ctype integrationElement(const LocalCoordinate& local) const
    {
        return integrationElement_[getSubElement(local)];
    }
 
    JacobianTransposed jacobianTransposed(const LocalCoordinate& local) const
    {
        return jacobianTransposed_[getSubElement(local)];
    }
 
    JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate& local) const
    {
        return jacobianInverseTransposed_[getSubElement(local)];
    }
 
    class JacobianInverseTransposed
        : public Dune::FieldMatrix< ctype, 3, 2 >
    {
        public:
        typedef Dune::FieldMatrix< ctype, 3, 2 > Base;
 
        void setup ( const JacobianTransposed &jt )
        {
            detInv_ = Dune::Impl::FieldMatrixHelper< ctype >::template rightInvA< 2, 3 >( jt, static_cast< Base & >( *this ) );
        }
 
        void setupDeterminant ( const JacobianTransposed &jt )
        {
            detInv_ = Dune::Impl::FieldMatrixHelper< ctype >::template sqrtDetAAT< 2, 3 >( jt );
        }
 
        ctype det () const { return ctype( 1 ) / detInv_; }
        ctype detInv () const { return detInv_; }
 
        private:
        ctype detInv_;
    };
 
    private:
    const std::array<GlobalCoordinate,6>    vertices_;
    std::array<JacobianTransposed,4>        linearMapSubElementLocalToWorldTransposed_;
    std::array<JacobianTransposed,4>        jacobianTransposed_;
    std::array<JacobianInverseTransposed,4> jacobianInverseTransposed_;
    std::array<GlobalCoordinate,4>          translation_;
    std::array<ctype,4>                     integrationElement_;
 
    Dune::FieldMatrix<std::size_t,4,3>      indexMap_;
    static int getSubElement(const LocalCoordinate& local)
    {
      if (local[0] + local[1] <= 0.5)
        return 0;
      else if (local[0] >= 0.5)
        return 1;
      else if (local[1] >= 0.5)
        return 2;
 
      return 3;
    }
 
    static void getSubElement(const LocalCoordinate& local,
                              int& subElement,
                              LocalCoordinate& subElementLocal)
    {
      if (local[0] + local[1] <= 0.5) {
        subElement = 0;
        subElementLocal[0] = 2*local[0];
        subElementLocal[1] = 2*local[1];
        return;
      } else if (local[0] >= 0.5) {
        subElement = 1;
        subElementLocal[0] = 2*local[0]-1;
        subElementLocal[1] = 2*local[1];
        return;
      } else if (local[1] >= 0.5) {
        subElement = 2;
        subElementLocal[0] = 2*local[0];
        subElementLocal[1] = 2*local[1]-1;
        return;
      }
 
      subElement = 3;
      subElementLocal[0] = -2 * local[0] + 1;
      subElementLocal[1] = -2 * local[1] + 1;
 
    }
};
 
template <class Element, class Transformation>
RefinedSimplexGeometry<typename Element::Geometry::ctype,2,3> makeRefinedSimplexGeometry(const Element& element, const Transformation& transformation)
{
    if(not element.type().isTriangle())
        DUNE_THROW(Dune::NotImplemented, "Free method makeRefinedSimplexGeometry(.,.) is only implemented for 2D simplex elements");
 
    std::array<typename Element::Geometry::GlobalCoordinate,6> movedVirtualVertices;
 
    typename Element::Geometry eltGeometry = element.geometry();
 
    for (std::size_t i=0; i<(std::size_t)eltGeometry.corners(); ++i)
        movedVirtualVertices[i] = eltGeometry.corner(i);
 
    for (std::size_t i=0; i<(std::size_t)element.subEntities(Element::dimension - 1); ++i)
    {
        auto y = transformation(element.template subEntity<Element::dimension -1>(i).geometry().center());
        movedVirtualVertices[eltGeometry.corners() + i] = y;
    }
 
    return RefinedSimplexGeometry<typename Element::Geometry::ctype,2,3>(movedVirtualVertices);
}
 
#endif