diff --git a/src/operators/hamilton-jacobi/upwind-eikonal/upwindEikonal.h b/src/operators/hamilton-jacobi/upwind-eikonal/upwindEikonal.h index 9a44fb41609394f6cb926a3573fe8db53034b083..10cb70ecde2b8aacc072641c4823e33a7054a82c 100644 --- a/src/operators/hamilton-jacobi/upwind-eikonal/upwindEikonal.h +++ b/src/operators/hamilton-jacobi/upwind-eikonal/upwindEikonal.h @@ -44,48 +44,47 @@ template< typename MeshReal, typename Index > class upwindEikonalScheme< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index > { - -public: - typedef Real RealType; - typedef Device DeviceType; - typedef Index IndexType; - typedef tnlGrid< 1, Real, Device, Index > MeshType; - typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType; - typedef typename MeshType::CoordinatesType CoordinatesType; - - - static tnlString getType(); - - template< typename PreimageFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const PreimageFunction& u, - const MeshEntity& entity, - const RealType& signU ) const - { - const RealType& hx_inv = entity.getMesh().template getSpaceStepsProducts< -1 >(); - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - - const RealType& u_c = u[ entity.getIndex() ]; - const RealType& u_e = u[ neighbourEntities.template getEntityIndex< 1 >() ]; - const RealType& u_w = u[ neighbourEntities.template getEntityIndex< -1 >() ]; - - if( signU > 0.0 ) - { - const RealType& xf = negativePart( ( u_e - u_c ) * hx_inv ); - const RealType& xb = positivePart( ( u_c - u_w ) * hx_inv ); - return sqrt( xf * xf + xb * xb ); - } - else if( signU < 0.0 ) - { - const RealType& xf = negativePart( ( u_c - u_w ) * hx_inv ); - const RealType& xb = positivePart( ( u_e - u_c ) * hx_inv ); - return sqrt( xf * xf + xb * xb ); - } - else + public: + typedef Real RealType; + typedef Device DeviceType; + typedef Index IndexType; + typedef tnlGrid< 1, Real, Device, Index > MeshType; + typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType; + typedef typename MeshType::CoordinatesType CoordinatesType; + + + static tnlString getType(); + + template< typename PreimageFunction, typename MeshEntity > + __cuda_callable__ + Real operator()( const PreimageFunction& u, + const MeshEntity& entity, + const RealType& signU ) const { - return 0.0; - } - }; + const RealType& hx_inv = entity.getMesh().template getSpaceStepsProducts< -1 >(); + const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); + + const RealType& u_c = u[ entity.getIndex() ]; + const RealType& u_e = u[ neighbourEntities.template getEntityIndex< 1 >() ]; + const RealType& u_w = u[ neighbourEntities.template getEntityIndex< -1 >() ]; + + if( signU > 0.0 ) + { + const RealType& xf = negativePart( ( u_e - u_c ) * hx_inv ); + const RealType& xb = positivePart( ( u_c - u_w ) * hx_inv ); + return sqrt( xf * xf + xb * xb ); + } + else if( signU < 0.0 ) + { + const RealType& xf = negativePart( ( u_c - u_w ) * hx_inv ); + const RealType& xb = positivePart( ( u_e - u_c ) * hx_inv ); + return sqrt( xf * xf + xb * xb ); + } + else + { + return 0.0; + } + }; }; /**** @@ -98,51 +97,50 @@ template< typename MeshReal, typename Index > class upwindEikonalScheme< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index > { - -public: - typedef Real RealType; - typedef Device DeviceType; - typedef Index IndexType; - typedef tnlGrid< 2, Real, Device, Index > MeshType; - typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType; - typedef typename MeshType::CoordinatesType CoordinatesType; - - static tnlString getType(); - - template< typename PreimageFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const PreimageFunction& u, - const MeshEntity& entity, - const RealType& signU ) const - { - const RealType& hx_inv = entity.getMesh().template getSpaceStepsProducts< -1, 0 >(); - const RealType& hy_inv = entity.getMesh().template getSpaceStepsProducts< 0, -1 >(); - - const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& u_c = u[ entity.getIndex() ]; - const RealType& u_e = u[ neighbourEntities.template getEntityIndex< 1, 0 >() ]; - const RealType& u_w = u[ neighbourEntities.template getEntityIndex< -1, 0 >() ]; - const RealType& u_n = u[ neighbourEntities.template getEntityIndex< 0, 1 >() ]; - const RealType& u_s = u[ neighbourEntities.template getEntityIndex< 0, -1 >() ]; - if( signU > 0.0 ) - { - const RealType xf = negativePart( ( u_e - u_c ) * hx_inv ); - const RealType xb = positivePart( ( u_c - u_w ) * hx_inv ); - const RealType yf = negativePart( ( u_n - u_c ) * hy_inv ); - const RealType yb = positivePart( ( u_c - u_s ) * hy_inv ); - return sqrt( xf * xf + xb * xb + yf * yf + yb * yb ); + public: + typedef Real RealType; + typedef Device DeviceType; + typedef Index IndexType; + typedef tnlGrid< 2, Real, Device, Index > MeshType; + typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType; + typedef typename MeshType::CoordinatesType CoordinatesType; + + static tnlString getType(); + + template< typename PreimageFunction, typename MeshEntity > + __cuda_callable__ + Real operator()( const PreimageFunction& u, + const MeshEntity& entity, + const RealType& signU ) const + { + const RealType& hx_inv = entity.getMesh().template getSpaceStepsProducts< -1, 0 >(); + const RealType& hy_inv = entity.getMesh().template getSpaceStepsProducts< 0, -1 >(); + + const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); + const RealType& u_c = u[ entity.getIndex() ]; + const RealType& u_e = u[ neighbourEntities.template getEntityIndex< 1, 0 >() ]; + const RealType& u_w = u[ neighbourEntities.template getEntityIndex< -1, 0 >() ]; + const RealType& u_n = u[ neighbourEntities.template getEntityIndex< 0, 1 >() ]; + const RealType& u_s = u[ neighbourEntities.template getEntityIndex< 0, -1 >() ]; + if( signU > 0.0 ) + { + const RealType xf = negativePart( ( u_e - u_c ) * hx_inv ); + const RealType xb = positivePart( ( u_c - u_w ) * hx_inv ); + const RealType yf = negativePart( ( u_n - u_c ) * hy_inv ); + const RealType yb = positivePart( ( u_c - u_s ) * hy_inv ); + return sqrt( xf * xf + xb * xb + yf * yf + yb * yb ); + } + else if( signU < 0.0 ) + { + const RealType xf = positivePart( ( u_e - u_c ) * hx_inv ); + const RealType xb = negativePart( ( u_c - u_w ) * hx_inv ); + const RealType yf = positivePart( ( u_n - u_c ) * hy_inv ); + const RealType yb = negativePart( ( u_c - u_s ) * hy_inv ); + return sqrt( xf * xf + xb * xb + yf * yf + yb * yb ); + } + else + return 0.0; } - else if( signU < 0.0 ) - { - const RealType xf = positivePart( ( u_e - u_c ) * hx_inv ); - const RealType xb = negativePart( ( u_c - u_w ) * hx_inv ); - const RealType yf = positivePart( ( u_n - u_c ) * hy_inv ); - const RealType yb = negativePart( ( u_c - u_s ) * hy_inv ); - return sqrt( xf * xf + xb * xb + yf * yf + yb * yb ); - } - else - return 0.0; - } }; /**** @@ -156,59 +154,59 @@ template< typename MeshReal, class upwindEikonalScheme< tnlGrid< 3,MeshReal, Device, MeshIndex >, Real, Index > { -public: - typedef Real RealType; - typedef Device DeviceType; - typedef Index IndexType; - typedef tnlGrid< 3, Real, Device, Index > MeshType; - typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType; - typedef typename MeshType::CoordinatesType CoordinatesType; - - - static tnlString getType(); - - template< typename PreimageFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const PreimageFunction& u, - const MeshEntity& entity, - const RealType& signU ) const - { - const RealType& hx_inv = entity.getMesh().template getSpaceStepsProducts< -1, 0, 0 >(); - const RealType& hy_inv = entity.getMesh().template getSpaceStepsProducts< 0, -1, 0 >(); - const RealType& hz_inv = entity.getMesh().template getSpaceStepsProducts< 0, 0, -1 >(); - - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& u_c = u[ entity.getIndex() ]; - const RealType& u_e = u[ neighbourEntities.template getEntityIndex< 1, 0, 0 >() ]; - const RealType& u_w = u[ neighbourEntities.template getEntityIndex< -1, 0, 0 >() ]; - const RealType& u_n = u[ neighbourEntities.template getEntityIndex< 0, 1, 0 >() ]; - const RealType& u_s = u[ neighbourEntities.template getEntityIndex< 0, -1, 0 >() ]; - const RealType& u_t = u[ neighbourEntities.template getEntityIndex< 0, 0, 1 >() ]; - const RealType& u_b = u[ neighbourEntities.template getEntityIndex< 0, 0, -1 >() ]; - - if( signU > 0.0 ) - { - const RealType xf = negativePart( ( u_e - u_c ) * hx_inv ); - const RealType xb = positivePart( ( u_c - u_w ) * hx_inv ); - const RealType yf = negativePart( ( u_n - u_c ) * hy_inv ); - const RealType yb = positivePart( ( u_c - u_s ) * hy_inv ); - const RealType zf = negativePart( ( u_t - u_c ) * hz_inv ); - const RealType zb = positivePart( ( u_c - u_b ) * hz_inv ); - return sqrt( xf * xf + xb * xb + yf * yf + yb * yb + zf * zf + zb * zb ); - } - else if( signU < 0.0 ) + public: + typedef Real RealType; + typedef Device DeviceType; + typedef Index IndexType; + typedef tnlGrid< 3, Real, Device, Index > MeshType; + typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType; + typedef typename MeshType::CoordinatesType CoordinatesType; + + + static tnlString getType(); + + template< typename PreimageFunction, typename MeshEntity > + __cuda_callable__ + Real operator()( const PreimageFunction& u, + const MeshEntity& entity, + const RealType& signU ) const { - const RealType xf = positivePart( ( u_e - u_c ) * hx_inv ); - const RealType xb = negativePart( ( u_c - u_w ) * hx_inv ); - const RealType yf = positivePart( ( u_n - u_c ) * hy_inv ); - const RealType yb = negativePart( ( u_c - u_s ) * hy_inv ); - const RealType zf = positivePart( ( u_t - u_c ) * hz_inv ); - const RealType zb = negativePart( ( u_c - u_b ) * hz_inv ); - return sqrt( xf * xf + xb * xb + yf * yf + yb * yb + zf * zf + zb * zb ); - } - else - return 0.0; - }; + const RealType& hx_inv = entity.getMesh().template getSpaceStepsProducts< -1, 0, 0 >(); + const RealType& hy_inv = entity.getMesh().template getSpaceStepsProducts< 0, -1, 0 >(); + const RealType& hz_inv = entity.getMesh().template getSpaceStepsProducts< 0, 0, -1 >(); + + const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); + const RealType& u_c = u[ entity.getIndex() ]; + const RealType& u_e = u[ neighbourEntities.template getEntityIndex< 1, 0, 0 >() ]; + const RealType& u_w = u[ neighbourEntities.template getEntityIndex< -1, 0, 0 >() ]; + const RealType& u_n = u[ neighbourEntities.template getEntityIndex< 0, 1, 0 >() ]; + const RealType& u_s = u[ neighbourEntities.template getEntityIndex< 0, -1, 0 >() ]; + const RealType& u_t = u[ neighbourEntities.template getEntityIndex< 0, 0, 1 >() ]; + const RealType& u_b = u[ neighbourEntities.template getEntityIndex< 0, 0, -1 >() ]; + + if( signU > 0.0 ) + { + const RealType xf = negativePart( ( u_e - u_c ) * hx_inv ); + const RealType xb = positivePart( ( u_c - u_w ) * hx_inv ); + const RealType yf = negativePart( ( u_n - u_c ) * hy_inv ); + const RealType yb = positivePart( ( u_c - u_s ) * hy_inv ); + const RealType zf = negativePart( ( u_t - u_c ) * hz_inv ); + const RealType zb = positivePart( ( u_c - u_b ) * hz_inv ); + return sqrt( xf * xf + xb * xb + yf * yf + yb * yb + zf * zf + zb * zb ); + } + else if( signU < 0.0 ) + { + const RealType xf = positivePart( ( u_e - u_c ) * hx_inv ); + const RealType xb = negativePart( ( u_c - u_w ) * hx_inv ); + const RealType yf = positivePart( ( u_n - u_c ) * hy_inv ); + const RealType yb = negativePart( ( u_c - u_s ) * hy_inv ); + const RealType zf = positivePart( ( u_t - u_c ) * hz_inv ); + const RealType zb = negativePart( ( u_c - u_b ) * hz_inv ); + return sqrt( xf * xf + xb * xb + yf * yf + yb * yb + zf * zf + zb * zb ); + } + else + return 0.0; + }; };