Commit 6d723a8b authored by Tomáš Oberhuber's avatar Tomáš Oberhuber
Browse files

Refactoring mean curvature flow.

parent ac16dd07

src/operators/diffusion/nonlinear-diffusion-operators/tnlFiniteVolumeNonlinearOperator.h

+9 −9
Original line number Diff line number Diff line
@@ -37,7 +37,7 @@ class tnlFiniteVolumeNonlinearOperator< tnlGrid< 1,MeshReal, Device, MeshIndex > 
   template< typename MeshEntity,

             typename Vector >

   __cuda_callable__

   Real getValue( const MeshEntity& entity,

   Real operator()( const MeshEntity& entity,

                    const Vector& u,

                    const RealType& time) const;
@@ -90,7 +90,7 @@ class tnlFiniteVolumeNonlinearOperator< tnlGrid< 2, MeshReal, Device, MeshIndex 
   template< typename MeshEntity,

             typename Vector >

   __cuda_callable__

   Real getValue( const MeshEntity& entity,

   Real operator()( const MeshEntity& entity,

                    const Vector& u,

                    const RealType& time) const;
@@ -143,7 +143,7 @@ class tnlFiniteVolumeNonlinearOperator< tnlGrid< 3, MeshReal, Device, MeshIndex 
   template< typename MeshEntity, 

             typename Vector >

   __cuda_callable__

   Real getValue( const MeshEntity& entity,

   Real operator()( const MeshEntity& entity,

                    const Vector& u,

                    const RealType& time) const;

src/operators/diffusion/nonlinear-diffusion-operators/tnlFiniteVolumeNonlinearOperator_impl.h

+51 −51
Original line number Diff line number Diff line
@@ -34,7 +34,7 @@ template< typename MeshEntity, 
__cuda_callable__

Real

tnlFiniteVolumeNonlinearOperator< tnlGrid< 1, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::

getValue( const MeshEntity& entity,

operator()( const MeshEntity& entity,

          const Vector& u,

          const Real& time ) const

{
@@ -111,18 +111,18 @@ template< typename MeshEntity, 
__cuda_callable__

Real

tnlFiniteVolumeNonlinearOperator< tnlGrid< 2, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::

getValue( const MeshEntity& entity,

operator()( const MeshEntity& entity,

            const Vector& u,

            const Real& time ) const

{

   const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities();      

   const typename MeshEntity::MeshType& mesh = entity.getMesh();

   const IndexType& cellIndex = entity.getIndex();

   return operatorQ.getValue( entity, u, time ) * 

      ( (  u[ neighbourEntities.template getEntityIndex<  1, 0 >() ] - u[ cellIndex ] ) * mesh.template getSpaceStepsProducts< -2, 0 >() / operatorQ.getValue( entity, u, time, 1 )

      + (  u[ neighbourEntities.template getEntityIndex<  0, 1 >() ] - u[ cellIndex ] ) * mesh.template getSpaceStepsProducts< 0, -2 >() / operatorQ.getValue( entity, u, time, 0, 1 ) 

      - ( -u[ neighbourEntities.template getEntityIndex< -1, 0 >() ] + u[ cellIndex ] ) * mesh.template getSpaceStepsProducts< -2, 0 >() / operatorQ.getValue( entity, u, time, -1)

      - ( -u[ neighbourEntities.template getEntityIndex<  0,-1 >() ] + u[ cellIndex ] ) * mesh.template getSpaceStepsProducts< 0, -2 >() / operatorQ.getValue( entity, u, time, 0, -1) );

   return operatorQ( entity, u, time ) * 

      ( (  u[ neighbourEntities.template getEntityIndex<  1, 0 >() ] - u[ cellIndex ] ) * mesh.template getSpaceStepsProducts< -2, 0 >() / operatorQ.operator()( entity, u, time, 1 )

      + (  u[ neighbourEntities.template getEntityIndex<  0, 1 >() ] - u[ cellIndex ] ) * mesh.template getSpaceStepsProducts< 0, -2 >() / operatorQ.operator()( entity, u, time, 0, 1 ) 

      - ( -u[ neighbourEntities.template getEntityIndex< -1, 0 >() ] + u[ cellIndex ] ) * mesh.template getSpaceStepsProducts< -2, 0 >() / operatorQ.operator()( entity, u, time, -1)

      - ( -u[ neighbourEntities.template getEntityIndex<  0,-1 >() ] + u[ cellIndex ] ) * mesh.template getSpaceStepsProducts< 0, -2 >() / operatorQ.operator()( entity, u, time, 0, -1) );

}



template< typename MeshReal,
@@ -166,19 +166,19 @@ updateLinearSystem( const RealType& time, 
{

   typename Matrix::MatrixRow matrixRow = matrix.getRow( index );

   const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities();

   const RealType aCoef = - tau * operatorQ.getValue( entity, u, time ) * mesh.template getSpaceStepsProducts< 0, -2 >() / 

                       operatorQ.getValue( entity, u, time, 0, -1 );

   const RealType bCoef = - tau * operatorQ.getValue( entity, u, time ) * mesh.template getSpaceStepsProducts< -2, 0 >() / 

                       operatorQ.getValue( entity, u, time, -1 );

   const RealType cCoef = tau * operatorQ.getValue( entity, u, time ) * ( mesh.template getSpaceStepsProducts< -2, 0 >() / 

                       operatorQ.getValue( entity, u, time, 1 ) + mesh.template getSpaceStepsProducts< 0, -2 >() / 

                       operatorQ.getValue( entity, u, time, 0, 1 )

                       + mesh.template getSpaceStepsProducts< -2, 0 >() / operatorQ.getValue( entity, u, time, -1 ) + 

                       mesh.template getSpaceStepsProducts< 0, -2 >() / operatorQ.getValue( entity, u, time, 0, -1 ) );

   const RealType dCoef = - tau * operatorQ.getValue( entity, u, time ) * mesh.template getSpaceStepsProducts< -2, 0 >() / 

                       operatorQ.getValue( entity, u, time, 1 );

   const RealType eCoef = - tau * operatorQ.getValue( entity, u, time ) * mesh.template getSpaceStepsProducts< 0, -2 >() / 

                       operatorQ.getValue(  entity, u, time, 0, 1 );

   const RealType aCoef = - tau * operatorQ.operator()( entity, u, time ) * mesh.template getSpaceStepsProducts< 0, -2 >() / 

                       operatorQ.operator()( entity, u, time, 0, -1 );

   const RealType bCoef = - tau * operatorQ.operator()( entity, u, time ) * mesh.template getSpaceStepsProducts< -2, 0 >() / 

                       operatorQ.operator()( entity, u, time, -1 );

   const RealType cCoef = tau * operatorQ.operator()( entity, u, time ) * ( mesh.template getSpaceStepsProducts< -2, 0 >() / 

                       operatorQ.operator()( entity, u, time, 1 ) + mesh.template getSpaceStepsProducts< 0, -2 >() / 

                       operatorQ.operator()( entity, u, time, 0, 1 )

                       + mesh.template getSpaceStepsProducts< -2, 0 >() / operatorQ.operator()( entity, u, time, -1 ) + 

                       mesh.template getSpaceStepsProducts< 0, -2 >() / operatorQ.operator()( entity, u, time, 0, -1 ) );

   const RealType dCoef = - tau * operatorQ.operator()( entity, u, time ) * mesh.template getSpaceStepsProducts< -2, 0 >() / 

                       operatorQ.operator()( entity, u, time, 1 );

   const RealType eCoef = - tau * operatorQ.operator()( entity, u, time ) * mesh.template getSpaceStepsProducts< 0, -2 >() / 

                       operatorQ.operator()(  entity, u, time, 0, 1 );

   matrixRow.setElement( 0, neighbourEntities.template getEntityIndex<  0, -1 >(), aCoef );

   matrixRow.setElement( 1, neighbourEntities.template getEntityIndex< -1,  0 >(), bCoef );

   matrixRow.setElement( 2, entity.getIndex(),                                     cCoef );
@@ -214,26 +214,26 @@ template< typename MeshEntity, 
__cuda_callable__

Real

tnlFiniteVolumeNonlinearOperator< tnlGrid< 3, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::

getValue( const MeshEntity& entity,

operator()( const MeshEntity& entity,

            const Vector& u,

            const Real& time ) const

{

   const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities();

   const typename MeshEntity::MeshType& mesh = entity.getMesh();

   const IndexType& cellIndex = entity.getIndex();

   return operatorQ.getValue( entity, u, time ) * 

   return operatorQ( entity, u, time ) * 

      ( (u[ neighbourEntities.template getEntityIndex< 1,0,0 >() ] - u[ cellIndex ]) 

          * mesh.template getSpaceStepsProducts< -2, 0, 0 >() / operatorQ.getValue( entity, u, time, 1 )

          * mesh.template getSpaceStepsProducts< -2, 0, 0 >() / operatorQ( entity, u, time, 1 )

          + ( u[ neighbourEntities.template getEntityIndex< 0,1,0 >() ] - u[ cellIndex ]) * mesh.template getSpaceStepsProducts< 0, -2, 0 >()/

          operatorQ.getValue( entity, u, time, 0, 1 ) 

          operatorQ( entity, u, time, 0, 1 ) 

          + ( u[ neighbourEntities.template getEntityIndex< 0,0,1 >() ] - u[ cellIndex ]) * mesh.template getSpaceStepsProducts< 0, 0, -2 >()/

          operatorQ.getValue( entity, u, time, 0, 0, 1 ) 

          operatorQ( entity, u, time, 0, 0, 1 ) 

          - ( - u[ neighbourEntities.template getEntityIndex< -1,0,0 >() ]  + u[ cellIndex ]) 

          * mesh.template getSpaceStepsProducts< -2, 0, 0 >() / operatorQ.getValue( entity, u, time, -1)

          * mesh.template getSpaceStepsProducts< -2, 0, 0 >() / operatorQ( entity, u, time, -1)

          -( - u[ neighbourEntities.template getEntityIndex< 0,-1,0 >() ] + u[ cellIndex ]) * mesh.template getSpaceStepsProducts< 0, -2, 0 >()

          /operatorQ.getValue( entity, u, time, 0, -1) 

          /operatorQ( entity, u, time, 0, -1) 

          -( - u[ neighbourEntities.template getEntityIndex< 0,0,-1 >() ] + u[ cellIndex ]) * mesh.template getSpaceStepsProducts< 0, 0, -2 >()

          /operatorQ.getValue( entity, u, time, 0, 0, -1) );

          /operatorQ( entity, u, time, 0, 0, -1) );

}



template< typename MeshReal,
@@ -279,25 +279,25 @@ updateLinearSystem( const RealType& time, 
{

   typename Matrix::MatrixRow matrixRow = matrix.getRow( index );

   const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities();

   const RealType aCoef = - tau * operatorQ.getValue( entity, u, time ) *

                       mesh.template getSpaceStepsProducts< 0, 0, -2 >() / operatorQ.getValue( entity, u, time, 0, 0, -1 );

   const RealType bCoef = - tau * operatorQ.getValue( entity, u, time ) * 

                       mesh.template getSpaceStepsProducts< 0, -2, 0 >() / operatorQ.getValue( entity, u, time, 0, -1, 0 );

   const RealType cCoef = - tau * operatorQ.getValue( entity, u, time ) * 

                       mesh.template getSpaceStepsProducts< -2, 0, 0 >() / operatorQ.getValue( entity, u, time, -1, 0, 0 );

   const RealType dCoef = tau * operatorQ.getValue( entity, u, time ) * ( 

                       mesh.template getSpaceStepsProducts< -2, 0, 0 >() / operatorQ.getValue( entity, u, time, 1, 0, 0 ) + 

                       mesh.template getSpaceStepsProducts< 0, -2, 0 >() / operatorQ.getValue( entity, u, time, 0, 1, 0 ) +

                       mesh.template getSpaceStepsProducts< 0, 0, -2 >() / operatorQ.getValue( entity, u, time, 0, 0, 1 ) + 

                       mesh.template getSpaceStepsProducts< -2, 0, 0 >() / operatorQ.getValue( entity, u, time, -1, 0, 0 ) +

                       mesh.template getSpaceStepsProducts< 0, -2, 0 >() / operatorQ.getValue( entity, u, time, 0, -1, 0 ) + 

                       mesh.template getSpaceStepsProducts< 0, 0, -2 >() / operatorQ.getValue( entity, u, time, 0, 0, -1 ) );

   const RealType eCoef = - tau * operatorQ.getValue( entity, u, time ) *

                       mesh.template getSpaceStepsProducts< -2, 0, 0 >() / operatorQ.getValue( entity, u, time, 1, 0, 0 );

   const RealType fCoef = - tau * operatorQ.getValue( entity, u, time ) *

                       mesh.template getSpaceStepsProducts< 0, -2, 0 >() / operatorQ.getValue( entity, u, time, 0, 1, 0 );

   const RealType gCoef = - tau * operatorQ.getValue( entity, u, time ) * 

                       mesh.template getSpaceStepsProducts< 0, 0, -2 >() / operatorQ.getValue( entity, u, time, 0, 0, 1 );

   const RealType aCoef = - tau * operatorQ( entity, u, time ) *

                       mesh.template getSpaceStepsProducts< 0, 0, -2 >() / operatorQ.operator()( entity, u, time, 0, 0, -1 );

   const RealType bCoef = - tau * operatorQ( entity, u, time ) * 

                       mesh.template getSpaceStepsProducts< 0, -2, 0 >() / operatorQ.operator()( entity, u, time, 0, -1, 0 );

   const RealType cCoef = - tau * operatorQ( entity, u, time ) * 

                       mesh.template getSpaceStepsProducts< -2, 0, 0 >() / operatorQ.operator()( entity, u, time, -1, 0, 0 );

   const RealType dCoef = tau * operatorQ( entity, u, time ) * ( 

                       mesh.template getSpaceStepsProducts< -2, 0, 0 >() / operatorQ.operator()( entity, u, time, 1, 0, 0 ) + 

                       mesh.template getSpaceStepsProducts< 0, -2, 0 >() / operatorQ.operator()( entity, u, time, 0, 1, 0 ) +

                       mesh.template getSpaceStepsProducts< 0, 0, -2 >() / operatorQ.operator()( entity, u, time, 0, 0, 1 ) + 

                       mesh.template getSpaceStepsProducts< -2, 0, 0 >() / operatorQ.operator()( entity, u, time, -1, 0, 0 ) +

                       mesh.template getSpaceStepsProducts< 0, -2, 0 >() / operatorQ.operator()( entity, u, time, 0, -1, 0 ) + 

                       mesh.template getSpaceStepsProducts< 0, 0, -2 >() / operatorQ.operator()( entity, u, time, 0, 0, -1 ) );

   const RealType eCoef = - tau * operatorQ.operator()( entity, u, time ) *

                       mesh.template getSpaceStepsProducts< -2, 0, 0 >() / operatorQ.operator()( entity, u, time, 1, 0, 0 );

   const RealType fCoef = - tau * operatorQ.operator()( entity, u, time ) *

                       mesh.template getSpaceStepsProducts< 0, -2, 0 >() / operatorQ.operator()( entity, u, time, 0, 1, 0 );

   const RealType gCoef = - tau * operatorQ.operator()( entity, u, time ) * 

                       mesh.template getSpaceStepsProducts< 0, 0, -2 >() / operatorQ.operator()( entity, u, time, 0, 0, 1 );

   matrixRow.setElement( 0, neighbourEntities.template getEntityIndex< 0,0,-1 >(), aCoef );

   matrixRow.setElement( 1, neighbourEntities.template getEntityIndex< 0,-1,0 >(), bCoef );

   matrixRow.setElement( 2, neighbourEntities.template getEntityIndex< -1,0,0 >(), cCoef );

src/operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator.h

+9 −9
Original line number Diff line number Diff line
@@ -37,7 +37,7 @@ class tnlOneSideDiffNonlinearOperator< tnlGrid< 1,MeshReal, Device, MeshIndex >, 
   template< typename Vector,

             typename MeshEntity >

   __cuda_callable__

   Real getValue( const MeshEntity& entity,

   Real operator()( const MeshEntity& entity,

                    const Vector& u,

                    const RealType& time) const;
@@ -90,7 +90,7 @@ class tnlOneSideDiffNonlinearOperator< tnlGrid< 2, MeshReal, Device, MeshIndex > 
   template< typename Vector,

             typename MeshEntity >

   __cuda_callable__

   Real getValue( const MeshEntity& entity,

   Real operator()( const MeshEntity& entity,

                    const Vector& u,

                    const RealType& time) const;
@@ -142,7 +142,7 @@ class tnlOneSideDiffNonlinearOperator< tnlGrid< 3, MeshReal, Device, MeshIndex > 
   template< typename Vector,

             typename MeshEntity >

   __cuda_callable__

   Real getValue( const MeshEntity& entity,

   Real operator()( const MeshEntity& entity,

                    const Vector& u,

                    const RealType& time) const;

src/operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator_impl.h

+63 −63

File changed.

Preview size limit exceeded, changes collapsed.

src/operators/diffusion/tnlNonlinearDiffusion.h

+9 −9
Original line number Diff line number Diff line
@@ -37,7 +37,7 @@ class tnlNonlinearDiffusion< tnlGrid< 1,MeshReal, Device, MeshIndex >, Nonlinear 
      template< typename MeshEntity,

                typename Vector >

      __cuda_callable__

      Real getValue( const MeshEntity& entity,

      Real operator()( const MeshEntity& entity,

                       const Vector& u,

                       const RealType& time) const;
@@ -90,7 +90,7 @@ class tnlNonlinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Nonlinea 
      template< typename MeshEntity,

                typename Vector >

      __cuda_callable__

      Real getValue( const MeshEntity& entity,

      Real operator()( const MeshEntity& entity,

                       const Vector& u,

                       const RealType& time) const;
@@ -142,7 +142,7 @@ class tnlNonlinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Nonlinea 
      template< typename MeshEntity,

                typename MeshFunction >

      __cuda_callable__

      Real getValue( const MeshEntity& entity,

      Real operator()( const MeshEntity& entity,

                       const MeshFunction& u,

                       const RealType& time) const;