Commit 5308a150 authored by Ján Bobot's avatar Ján Bobot Committed by Jakub Klinkovský
Browse files

Added new getEntityMeasure specializations and function isPlanar for 3D polygons + refactoring 

- added getEntityMeasure specializations for 2D and 3D polygon, pyramid, wedge, polyhedron
- added unit tests for new getEntityMeasure specializations
- added function isPlanar for checking whether all points of a 3D polygon are on the same plane
- added unit test for function isPlanar
- refactored isDynamicTopology type trait to avoid warning
parent 9f5c2dda

src/TNL/Meshes/Geometry/getEntityCenter.h

+1 −1
Original line number Diff line number Diff line
@@ -50,7 +50,7 @@ getEntityCenter( const Mesh< MeshConfig, Device > & mesh, 
                 const MeshEntity< MeshConfig, Device, EntityTopology > & entity )

{

   using EntityType = MeshEntity< MeshConfig, Device, EntityTopology >;

   constexpr typename MeshConfig::LocalIndexType subvertices = EntityType::template SubentityTraits< 0 >::count;

   const typename MeshConfig::LocalIndexType subvertices = entity.template getSubentitiesCount< 0 >();

   typename MeshTraits< MeshConfig >::PointType c = 0;

   for( typename MeshConfig::LocalIndexType i = 0;

        i < subvertices;

src/TNL/Meshes/Geometry/getEntityMeasure.h

+177 −0
Original line number Diff line number Diff line
@@ -14,12 +14,17 @@ 
#include <TNL/Meshes/GridEntity.h>

#include <TNL/Meshes/Mesh.h>

#include <TNL/Meshes/MeshEntity.h>

#include <TNL/Meshes/Geometry/getOutwardNormalVector.h>

#include <TNL/Meshes/Topologies/Vertex.h>

#include <TNL/Meshes/Topologies/Edge.h>

#include <TNL/Meshes/Topologies/Triangle.h>

#include <TNL/Meshes/Topologies/Quadrangle.h>

#include <TNL/Meshes/Topologies/Tetrahedron.h>

#include <TNL/Meshes/Topologies/Hexahedron.h>

#include <TNL/Meshes/Topologies/Polygon.h>

#include <TNL/Meshes/Topologies/Wedge.h>

#include <TNL/Meshes/Topologies/Pyramid.h>

#include <TNL/Meshes/Topologies/Polyhedron.h>



namespace TNL {

namespace Meshes {
@@ -172,5 +177,177 @@ getEntityMeasure( const Mesh< MeshConfig, Device > & mesh, 
         + getTetrahedronVolume( v6 - v4, v2 - v4, v7 - v4 );

}



// Polygon

template< int Coord1,

          int Coord2,

          typename MeshConfig,

          typename Device >

__cuda_callable__

typename MeshConfig::RealType

getPolygon2DArea( const Mesh< MeshConfig, Device > & mesh,

                  const MeshEntity< MeshConfig, Device, Topologies::Polygon > & entity )

{

    // http://geomalgorithms.com/code.html (function area2D_Polygon)



    using Real = typename MeshConfig::RealType;

    using Index = typename MeshConfig::LocalIndexType;



    Real area{ 0.0 };

    const auto n = entity.template getSubentitiesCount< 0 >();

    for ( Index i = 1, j = 2, k = 0; j < n; i++, j++, k++ ) {

        const auto& v0 = mesh.getPoint( entity.template getSubentityIndex< 0 >( i ) );

        const auto& v1 = mesh.getPoint( entity.template getSubentityIndex< 0 >( j ) );

        const auto& v2 = mesh.getPoint( entity.template getSubentityIndex< 0 >( k ) );

        area += v0[Coord1] * ( v1[Coord2] - v2[Coord2] );

    }



    // 1. wrap around term

    {

        const auto& v0 = mesh.getPoint( entity.template getSubentityIndex< 0 >( n - 1 ) );

        const auto& v1 = mesh.getPoint( entity.template getSubentityIndex< 0 >( 0 ) );

        const auto& v2 = mesh.getPoint( entity.template getSubentityIndex< 0 >( n - 2 ) );

        area += v0[Coord1] * ( v1[Coord2] - v2[Coord2] );

    }



    // 2. wrap around term

    {

        const auto& v0 = mesh.getPoint( entity.template getSubentityIndex< 0 >( 0 ) );

        const auto& v1 = mesh.getPoint( entity.template getSubentityIndex< 0 >( 1 ) );

        const auto& v2 = mesh.getPoint( entity.template getSubentityIndex< 0 >( n - 1 ) );

        area += v0[Coord1] * ( v1[Coord2] - v2[Coord2] );

    }



    return Real( 0.5 ) * area;

}



template< typename MeshConfig,

          typename Device,

          std::enable_if_t< MeshConfig::spaceDimension == 2, bool > = true >

__cuda_callable__

typename MeshConfig::RealType

getEntityMeasure( const Mesh< MeshConfig, Device > & mesh,

                  const MeshEntity< MeshConfig, Device, Topologies::Polygon > & entity )

{

    const auto area = getPolygon2DArea< 0, 1 >( mesh, entity );

    return TNL::abs( area );

}



template< typename MeshConfig,

          typename Device,

          std::enable_if_t< MeshConfig::spaceDimension == 3, bool > = true >

__cuda_callable__

typename MeshConfig::RealType

getEntityMeasure( const Mesh< MeshConfig, Device > & mesh,

                  const MeshEntity< MeshConfig, Device, Topologies::Polygon > & entity )



{

    // http://geomalgorithms.com/code.html (function area3D_Polygon)



    using Real = typename MeshConfig::RealType;



    // select largest abs coordinate of normal vector to ignore for projection

    auto normal = getNormalVector( mesh, entity );

    normal = TNL::abs( normal );

    int coord = 2;  // ignore z-coord

    if ( normal.x() > normal.y() ) {

        if ( normal.x() > normal.z() ) coord = 0;  // ignore x-coord

    }

    else if ( normal.y() > normal.z() ) coord = 1; // ignore y-coord



    Real area;

    switch( coord ) {

        case 0: // ignored x-coord

            area = getPolygon2DArea< 1, 2 >( mesh, entity );

            area *= l2Norm( normal ) / normal.x();

            break;

        case 1: // ignored y-coord

            area = getPolygon2DArea< 0, 2 >( mesh, entity );

            area *= l2Norm( normal ) / normal.y();

            break;

        default: // ignored z-coord

            area = getPolygon2DArea< 0, 1 >( mesh, entity );

            area *= l2Norm( normal ) / normal.z();

            break;

    }

    return TNL::abs( area );

}



// Wedge

template< typename MeshConfig,

          typename Device >

__cuda_callable__

typename MeshConfig::RealType

getEntityMeasure( const Mesh< MeshConfig, Device > & mesh,

                  const MeshEntity< MeshConfig, Device, Topologies::Wedge > & entity )

{

    using Real = typename MeshConfig::RealType;



    const auto& v0 = mesh.getPoint( entity.template getSubentityIndex< 0 >( 0 ) );

    const auto& v1 = mesh.getPoint( entity.template getSubentityIndex< 0 >( 1 ) );

    const auto& v2 = mesh.getPoint( entity.template getSubentityIndex< 0 >( 2 ) );

    const auto& v3 = mesh.getPoint( entity.template getSubentityIndex< 0 >( 3 ) );

    const auto& v4 = mesh.getPoint( entity.template getSubentityIndex< 0 >( 4 ) );

    const auto& v5 = mesh.getPoint( entity.template getSubentityIndex< 0 >( 5 ) );

    // Partition wedge into three tetrahedrons.

    return getTetrahedronVolume( v2 - v3, v0 - v3, v1 - v3 )

         + getTetrahedronVolume( v2 - v3, v1 - v3, v4 - v3 )

         + getTetrahedronVolume( v2 - v3, v4 - v3, v5 - v3 );

}



// Pyramid

template< typename MeshConfig,

          typename Device >

__cuda_callable__

typename MeshConfig::RealType

getEntityMeasure( const Mesh< MeshConfig, Device > & mesh,

                  const MeshEntity< MeshConfig, Device, Topologies::Pyramid > & entity )

{

    using Real = typename MeshConfig::RealType;



    const auto& v0 = mesh.getPoint( entity.template getSubentityIndex< 0 >( 0 ) );

    const auto& v1 = mesh.getPoint( entity.template getSubentityIndex< 0 >( 1 ) );

    const auto& v2 = mesh.getPoint( entity.template getSubentityIndex< 0 >( 2 ) );

    const auto& v3 = mesh.getPoint( entity.template getSubentityIndex< 0 >( 3 ) );

    const auto& v4 = mesh.getPoint( entity.template getSubentityIndex< 0 >( 4 ) );

    // Partition pyramid into two tetrahedrons.

    return getTetrahedronVolume( v4 - v0, v3 - v0, v1 - v0 )

         + getTetrahedronVolume( v4 - v2, v1 - v2, v3 - v2 );

}



// Polyhedron

template< typename MeshConfig,

          typename Device >

__cuda_callable__

typename MeshConfig::RealType

getEntityMeasure( const Mesh< MeshConfig, Device > & mesh,

                  const MeshEntity< MeshConfig, Device, Topologies::Polyhedron > & entity )

{

    using Real = typename MeshConfig::RealType;

    using Index = typename MeshConfig::LocalIndexType;

    Real volume{ 0.0 };

    const Index facesCount = entity.template getSubentitiesCount< 2 >();

    for( Index faceIdx = 0; faceIdx < facesCount; faceIdx++ ) {

        const auto face = mesh.template getEntity< 2 >( entity.template getSubentityIndex< 2 >( faceIdx ) );

        const Index verticesCount = face.template getSubentitiesCount< 0 >();

        const auto& v0 = mesh.getPoint( face.template getSubentityIndex< 0 >( 0 ) );

        for( Index i = 1, j = 2; j < verticesCount; i++, j++ ) {

            const auto& v1 = mesh.getPoint( face.template getSubentityIndex< 0 >( i ) );

            const auto& v2 = mesh.getPoint( face.template getSubentityIndex< 0 >( j ) );

            // Partition polyhedron into tetrahedrons by triangulating faces and connecting each triangle to the origin point (0,0,0).

            // It is required that vertices of all faces are stored consistently in CW or CCW order as faces are viewed from the outside.

            // Otherwise signs of some tetrahedron volumes may be incorrect, resulting in overall incorrect volume.

            // https://stackoverflow.com/a/1849746



            // volume += dot(v0 x v1, v2)

            volume += Real {

                  ( v0.y() * v1.z() - v0.z() * v1.y() ) * v2.x()

                + ( v0.z() * v1.x() - v0.x() * v1.z() ) * v2.y()

                + ( v0.x() * v1.y() - v0.y() * v1.x() ) * v2.z()

            };

        }

    }

    return Real{ 1.0 / 6.0 } * TNL::abs( volume );

}



} // namespace Meshes

} // namespace TNL

src/TNL/Meshes/Geometry/getOutwardNormalVector.h

+22 −10
Original line number Diff line number Diff line
@@ -117,19 +117,14 @@ getOutwardNormalVector( const Mesh< MeshConfig, Device > & mesh, 
template< typename MeshConfig, typename Device, typename EntityTopology >

__cuda_callable__

typename MeshTraits< MeshConfig >::PointType

getOutwardNormalVector( const Mesh< MeshConfig, Device > & mesh,

                        const MeshEntity< MeshConfig, Device, EntityTopology > & face,

                        typename MeshTraits< MeshConfig >::PointType cellCenter )

getNormalVector( const Mesh< MeshConfig, Device > & mesh,

                 const MeshEntity< MeshConfig, Device, EntityTopology > & entity )

{

   using MeshType = Mesh< MeshConfig, Device >;

   using FaceType = MeshEntity< MeshConfig, Device, EntityTopology >;

   using PointType = typename MeshTraits< MeshConfig >::PointType;

   static_assert( std::is_same< typename MeshType::Face, FaceType >::value, "getOutwardNormalVector called for an entity which is not a face" );

   static_assert( MeshConfig::spaceDimension == 3, "general overload intended for 3D was called with the wrong space dimension" );

   

   const auto& v0 = mesh.getPoint( face.template getSubentityIndex< 0 >( 0 ) );

   const auto& v1 = mesh.getPoint( face.template getSubentityIndex< 0 >( 1 ) );

   const auto& v2 = mesh.getPoint( face.template getSubentityIndex< 0 >( 2 ) );

   const auto& v0 = mesh.getPoint( entity.template getSubentityIndex< 0 >( 0 ) );

   const auto& v1 = mesh.getPoint( entity.template getSubentityIndex< 0 >( 1 ) );

   const auto& v2 = mesh.getPoint( entity.template getSubentityIndex< 0 >( 2 ) );

   const PointType u1 = v0 - v1;

   const PointType u2 = v0 - v2;

   const PointType n {
@@ -137,6 +132,23 @@ getOutwardNormalVector( const Mesh< MeshConfig, Device > & mesh, 
      u1.z() * u2.x() - u1.x() * u2.z(),   // second component of the cross product

      u1.x() * u2.y() - u1.y() * u2.x()    // third component of the cross product

   };

   return n;

}



template< typename MeshConfig, typename Device, typename EntityTopology >

__cuda_callable__

typename MeshTraits< MeshConfig >::PointType

getOutwardNormalVector( const Mesh< MeshConfig, Device > & mesh,

                        const MeshEntity< MeshConfig, Device, EntityTopology > & face,

                        typename MeshTraits< MeshConfig >::PointType cellCenter )

{

   using MeshType = Mesh< MeshConfig, Device >;

   using FaceType = MeshEntity< MeshConfig, Device, EntityTopology >;

   using PointType = typename MeshTraits< MeshConfig >::PointType;

   static_assert( std::is_same< typename MeshType::Face, FaceType >::value, "getOutwardNormalVector called for an entity which is not a face" );

   static_assert( MeshConfig::spaceDimension == 3, "general overload intended for 3D was called with the wrong space dimension" );



   const PointType n = getNormalVector( mesh, face );



   // check on which side of the face is the reference cell center

   const PointType faceCenter = getEntityCenter( mesh, face );

src/TNL/Meshes/Geometry/isPlanar.h

0 → 100644
+34 −0
Original line number Diff line number Diff line 
#pragma once



#include <TNL/Meshes/Geometry/getEntityMeasure.h>



namespace TNL {

namespace Meshes {



// Polygon

template< typename MeshConfig,

          typename Device,

          std::enable_if_t< MeshConfig::spaceDimension == 3, bool > = true >

__cuda_callable__

bool

isPlanar( const Mesh< MeshConfig, Device > & mesh,

          const MeshEntity< MeshConfig, Device, Topologies::Polygon > & entity,

          const typename MeshConfig::RealType precision )

{

   using Real = typename MeshConfig::RealType;

   using Index = typename MeshConfig::LocalIndexType;

   const auto& v0 = mesh.getPoint( entity.template getSubentityIndex< 0 >( 0 ) );

   const auto& v1 = mesh.getPoint( entity.template getSubentityIndex< 0 >( 1 ) );

   const Index verticesCount = entity.template getSubentitiesCount< 0 >();

   for( Index i = 2, j = 3; j < verticesCount; i++, j++ ) {

      const auto& v2 = mesh.getPoint( entity.template getSubentityIndex< 0 >( i ) );

      const auto& v3 = mesh.getPoint( entity.template getSubentityIndex< 0 >( j ) );

      const Real volume{ getTetrahedronVolume( v0 - v1, v2 - v1, v3 - v1 ) };

      if( volume > precision )

         return false;

   }

   return true;

}



} // namespace Meshes

} // namespace TNL

src/TNL/Meshes/MeshDetails/initializer/EntitySeed.h

+2 −2
Original line number Diff line number Diff line
@@ -141,8 +141,8 @@ public: 
   {

      if( std::is_same< EntityTopology, Topologies::Polygon >::value )

         TNL_ASSERT_GE( cornersCount, 3, "polygons must have at least 3 corners" );

      else if( std::is_same< EntityTopology, Topologies::Polyhedron >::value )

         TNL_ASSERT_GE( cornersCount, 2, "polyhedron must have at least 2 faces" );

      /*else if( std::is_same< EntityTopology, Topologies::Polyhedron >::value )

         TNL_ASSERT_GE( cornersCount, 2, "polyhedron must have at least 2 faces" );*/



      this->cornerIds.setSize( cornersCount );

   }