Commit 8b466634 authored by Tomáš Oberhuber's avatar Tomáš Oberhuber
Browse files

Adding function generators and finite differences.

parent 3107e03f
Loading
Loading
Loading
Loading
+1 −1
Original line number Diff line number Diff line
@@ -33,7 +33,7 @@ void tnlFunctionDiscretizer< Mesh, Function, DiscreteFunction >::discretize( con
      typename Mesh::CoordinatesType c;
      mesh.getElementCoordinates( i, c );
      mesh.getElementCenter( c, v );
      discreteFunction[ i ] = function.template getF< XDiffOrder, YDiffOrder, XDiffOrder >( v );
      discreteFunction[ i ] = function.template getF< XDiffOrder, YDiffOrder, ZDiffOrder >( v );
      i++;
   }
}
+7 −3
Original line number Diff line number Diff line
@@ -43,8 +43,12 @@ template< typename Vertex, typename Device >
      typename Vertex::RealType tnlSinWaveFunction< 1, Vertex, Device >::getF( const Vertex& v ) const
{
   const RealType& x = v.x();
   if( XDiffOrder == 0 && YDiffOrder == 0 && ZDiffOrder == 0 )
      return this->amplitude * sin( this->phase + x / ( M_PI * this->waveLength) );
   if( YDiffOrder != 0 || ZDiffOrder != 0 )
      return 0.0;
   if( XDiffOrder == 0 )
      return this->amplitude * sin( this->phase + 2.0 * M_PI * x / this->waveLength );
   if( XDiffOrder == 1 )
      return 2.0 * M_PI / this->waveLength * this->amplitude * cos( 2.0 * M_PI * x / this->waveLength );
   return 0.0;
}

@@ -57,7 +61,7 @@ template< typename Vertex, typename Device >
   const RealType& y = v.y();
   if( XDiffOrder == 0 && YDiffOrder == 0 && ZDiffOrder == 0 )
   {
      return this->amplitude * sin( this->phase + sqrt( x * x + y * y ) / ( M_PI * this->waveLength) );
      return this->amplitude * sin( this->phase + 2.0 * M_PI * sqrt( x * x + y * y ) / this->waveLength );
   }
   return 0.0;
}
+136 −56
Original line number Diff line number Diff line
@@ -20,93 +20,173 @@

#include <schemes/tnlFiniteDifferences.h>


template< typename Real, typename Device, typename Index >
   template< typename GridFunction >
      Real tnlFiniteDifferences< tnlGrid< 1, Real, Device, Index, tnlIdenticalGridGeometry > >::
         getFirstForwardDifferenceX( const GridType& grid,
                                     const CoordinatesType& c,
                                     const GridFunction& function )
{
   tnlAssert( c.x() < grid.getDimensions().x() - 1,
             cerr << "c.x() = " << c.x() << " grid.getDimensions().x() - 1 = " << grid.getDimensions().x() - 1 );
   Real hx = grid.getParametricStep().x();
   Index i1 = grid.getElementIndex( c.x() );
   Index i2 = grid.getElementIndex( c.x() + 1 );
   return 1.0 / hx * ( function[ i2 ] - function[ i1 ] );
}

template< typename Real, typename Device, typename Index >
   template< typename GridFunction >
      void tnlFiniteDifferences< tnlGrid< 1, Real, Device, Index, tnlIdenticalGridGeometry > >::
         getFirstForwardDifferenceX( const GridType& grid,
   template< typename GridFunction,
             int XDifferenceOrder,
             int YDifferenceOrder,
             int ZDifferenceOrder,
             int XDifferenceDirection,
             int YDifferenceDirection,
             int ZDifferenceDirection >
Real tnlFiniteDifferences< tnlGrid< 1, Real, Device, Index, tnlIdenticalGridGeometry > >::getDifference( const GridType& grid,
                                                                                                         const GridFunction& inFunction,
                                                                                                         GridFunction& outFunction )
{
   for( Index i = 0; i < grid.getDimensions().x() - 1; i++ )
   {
   IndexType iBegin, iEnd;
   if( XDifferenceDirection == 0 || XDifferenceDirection == 1 )
      iBegin = 0;
   else
      iBegin = 1;
   if( XDifferenceDirection == 1 )
      iEnd = grid.getDimensions().x() - 1;
   else
      iEnd = grid.getDimensions().x();

   CoordinatesType c;
      c.x() = i;
      Index k = grid.getElementIndex( c );
      outFunction[ k ] = this->getFirstForwardDifferenceX( grid, c, inFunciton );
   for( c.x() = iBegin; c.x() < iEnd; c.x()++ )
   {
      outFunction[ grid.getElementIndex( c.x() ) ] =
               getDifference< GridFunction,
                              XDifferenceOrder,
                              YDifferenceOrder,
                              ZDifferenceOrder,
                              XDifferenceDirection,
                              YDifferenceDirection,
                              ZDifferenceDirection >( grid, c, inFunction );
   }
}

template< typename Real, typename Device, typename Index >
   template< typename GridFunction >
      Real tnlFiniteDifferences< tnlGrid< 1, Real, Device, Index, tnlIdenticalGridGeometry > >::
         getFirstBackwardDifferenceX( const GridType& grid,
   template< typename GridFunction,
             int XDifferenceOrder,
             int YDifferenceOrder,
             int ZDifferenceOrder,
             int XDifferenceDirection,
             int YDifferenceDirection,
             int ZDifferenceDirection >
Real tnlFiniteDifferences< tnlGrid< 1, Real, Device, Index, tnlIdenticalGridGeometry > >::getDifference( const GridType& grid,
                                                                                                         const CoordinatesType& c,
                                      GridFunction& function )
                                                                                                         const GridFunction& function )
{

   if( YDifferenceOrder > 0 || ZDifferenceOrder > 0 )
      return 0.0;
   const RealType hx = grid.getParametricStep().x();
   if( XDifferenceOrder == 1 )
   {
      if( XDifferenceDirection == 0 )
         return ( function[ grid.getElementIndex( c.x() + 1 ) ] -
                  function[ grid.getElementIndex( c.x() - 1 ) ] ) / ( 2.0 * hx );
      else
         return ( function[ grid.getElementIndex( c.x() + XDifferenceDirection ) ] -
                  function[ grid.getElementIndex( c.x() ) ] ) / ( XDifferenceDirection * hx );
   }
   if( XDifferenceOrder == 2 )
   {
      return ( function[ grid.getElementIndex( c.x() + 1 ) ] -
               2.0 * function[ grid.getElementIndex( c.x() ) ] +
               function[ grid.getElementIndex( c.x() - 1 ) ] ) / (  hx * hx );
   }

}

/****
 *  2D Grid
 */

template< typename Real, typename Device, typename Index >
   template< typename GridFunction >
      Real tnlFiniteDifferences< tnlGrid< 1, Real, Device, Index, tnlIdenticalGridGeometry > >::
         getFirstBackwardDifferenceX( const GridType& grid,
                                      GridFunction& function )
   template< typename GridFunction,
             int XDifferenceOrder,
             int YDifferenceOrder,
             int ZDifferenceOrder,
             int XDifferenceDirection,
             int YDifferenceDirection,
             int ZDifferenceDirection >
Real tnlFiniteDifferences< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry > >::getDifference( const GridType& grid,
                                                                                                         const GridFunction& inFunction,
                                                                                                         GridFunction& outFunction )
{

}

template< typename Real, typename Device, typename Index >
   template< typename GridFunction >
      Real tnlFiniteDifferences< tnlGrid< 1, Real, Device, Index, tnlIdenticalGridGeometry > >::
         getFirstCentralDifferenceX( const GridType& grid,
   template< typename GridFunction,
             int XDifferenceOrder,
             int YDifferenceOrder,
             int ZDifferenceOrder,
             int XDifferenceDirection,
             int YDifferenceDirection,
             int ZDifferenceDirection >
Real tnlFiniteDifferences< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry > >::getDifference( const GridType& grid,
                                                                                                         const CoordinatesType& c,
                                     GridFunction& function )
                                                                                                         const GridFunction& function )
{
   if( YDifferenceOrder > 0 || ZDifferenceOrder > 0 )
      return 0.0;
   const RealType hx = grid.getParametricSpaceStep();
   if( XDifferenceOrder == 1 )
   {
      return ( function[ grid.getElementIndex( c.x() + XDifferenceDirection ) ] -
               function[ grid.getElementIndex( c.x() ) ] ) / ( XDifferenceDirection * hx );
   }
   if( XDifferenceOrder == 2 )
   {
      return ( function[ grid.getElementIndex( c.x() + 1 ) ] -
               2.0 * function[ grid.getElementIndex( c.x() ) ] +
               function[ grid.getElementIndex( c.x() - 1 ) ] ) / (  hx * hx );
   }

}

/****
 *  3D Grid
 */

template< typename Real, typename Device, typename Index >
   template< typename GridFunction >
      Real tnlFiniteDifferences< tnlGrid< 1, Real, Device, Index, tnlIdenticalGridGeometry > >::
         getFirstCentralDifferenceX( const GridType& grid,
                                     GridFunction& function )
   template< typename GridFunction,
             int XDifferenceOrder,
             int YDifferenceOrder,
             int ZDifferenceOrder,
             int XDifferenceDirection,
             int YDifferenceDirection,
             int ZDifferenceDirection >
Real tnlFiniteDifferences< tnlGrid< 3, Real, Device, Index, tnlIdenticalGridGeometry > >::getDifference( const GridType& grid,
                                                                                                         const GridFunction& inFunction,
                                                                                                         GridFunction& outFunction )
{

}

template< typename Real, typename Device, typename Index >
   template< typename GridFunction >
      Real tnlFiniteDifferences< tnlGrid< 1, Real, Device, Index, tnlIdenticalGridGeometry > >::
         getSecondCentralDifferenceX( const GridType& grid,
   template< typename GridFunction,
             int XDifferenceOrder,
             int YDifferenceOrder,
             int ZDifferenceOrder,
             int XDifferenceDirection,
             int YDifferenceDirection,
             int ZDifferenceDirection >
Real tnlFiniteDifferences< tnlGrid< 3, Real, Device, Index, tnlIdenticalGridGeometry > >::getDifference( const GridType& grid,
                                                                                                         const CoordinatesType& c,
                                      GridFunction& function )
                                                                                                         const GridFunction& function )
{

   if( YDifferenceOrder > 0 || ZDifferenceOrder > 0 )
      return 0.0;
   const RealType hx = grid.getParametricSpaceStep();
   if( XDifferenceOrder == 1 )
   {
      return ( function[ grid.getElementIndex( c.x() + XDifferenceDirection ) ] -
               function[ grid.getElementIndex( c.x() ) ] ) / ( XDifferenceDirection * hx );
   }

template< typename Real, typename Device, typename Index >
   template< typename GridFunction >
      Real tnlFiniteDifferences< tnlGrid< 1, Real, Device, Index, tnlIdenticalGridGeometry > >::
         getSecondCentralDifferenceX( const GridType& grid,
                                      GridFunction& function )
   if( XDifferenceOrder == 2 )
   {
      return ( function[ grid.getElementIndex( c.x() + 1 ) ] -
               2.0 * function[ grid.getElementIndex( c.x() ) ] +
               function[ grid.getElementIndex( c.x() - 1 ) ] ) / (  hx * hx );
   }

}


#include <implementation/schemes/tnlFiniteDifferences_impl.h>


+72 −2
Original line number Diff line number Diff line
@@ -26,7 +26,7 @@ class tnlFiniteDifferences
{
};

template< typename Real = double, typename Device = tnlHost, typename Index = int >
template< typename Real, typename Device, typename Index >
class tnlFiniteDifferences< tnlGrid< 1, Real, Device, Index, tnlIdenticalGridGeometry > >
{
   public:
@@ -34,8 +34,8 @@ class tnlFiniteDifferences< tnlGrid< 1, Real, Device, Index, tnlIdenticalGridGeo
   typedef Real RealType;
   typedef Device DeviceType;
   typedef Index IndexType;
   typedef tnlTuple< 1, Index > CoordinatesType;
   typedef tnlGrid< 1, Real, Device, Index, tnlIdenticalGridGeometry > GridType;
   typedef typename GridType::CoordinatesType CoordinatesType;

   template< typename GridFunction,
             int XDifferenceOrder,
@@ -61,6 +61,76 @@ class tnlFiniteDifferences< tnlGrid< 1, Real, Device, Index, tnlIdenticalGridGeo

};

template< typename Real, typename Device, typename Index >
class tnlFiniteDifferences< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry > >
{
   public:

   typedef Real RealType;
   typedef Device DeviceType;
   typedef Index IndexType;
   typedef tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry > GridType;
   typedef typename GridType::CoordinatesType CoordinatesType;

   template< typename GridFunction,
             int XDifferenceOrder,
             int YDifferenceOrder,
             int ZDifferenceOrder,
             int XDifferenceDirection = XDifferenceOrder%2,
             int YDifferenceDirection = YDifferenceOrder%2,
             int ZDifferenceDirection = ZDifferenceOrder%2 >
   static RealType getDifference( const GridType& grid,
                                  const GridFunction& inFunction,
                                  GridFunction& outFunction );

   template< typename GridFunction,
             int XDifferenceOrder,
             int YDifferenceOrder,
             int ZDifferenceOrder,
             int XDifferenceDirection = 0,
             int YDifferenceDirection = 0,
             int ZDifferenceDirection = 0 >
   static RealType getDifference( const GridType& grid,
                                  const CoordinatesType& c,
                                  const GridFunction& function );

};

template< typename Real, typename Device, typename Index >
class tnlFiniteDifferences< tnlGrid< 3, Real, Device, Index, tnlIdenticalGridGeometry > >
{
   public:

   typedef Real RealType;
   typedef Device DeviceType;
   typedef Index IndexType;
   typedef tnlGrid< 3, Real, Device, Index, tnlIdenticalGridGeometry > GridType;
   typedef typename GridType::CoordinatesType CoordinatesType;

   template< typename GridFunction,
             int XDifferenceOrder,
             int YDifferenceOrder,
             int ZDifferenceOrder,
             int XDifferenceDirection = XDifferenceOrder%2,
             int YDifferenceDirection = YDifferenceOrder%2,
             int ZDifferenceDirection = ZDifferenceOrder%2 >
   static RealType getDifference( const GridType& grid,
                                  const GridFunction& inFunction,
                                  GridFunction& outFunction );

   template< typename GridFunction,
             int XDifferenceOrder,
             int YDifferenceOrder,
             int ZDifferenceOrder,
             int XDifferenceDirection = 0,
             int YDifferenceDirection = 0,
             int ZDifferenceDirection = 0 >
   static RealType getDifference( const GridType& grid,
                                  const CoordinatesType& c,
                                  const GridFunction& function );

};

#include <implementation/schemes/tnlFiniteDifferences_impl.h>


+2 −13
Original line number Diff line number Diff line
set( headers tnl-functions-test.h )
             
if( BUILD_CUDA )
    CUDA_ADD_EXECUTABLE( tnl-functions-test${mpiExt}${debugExt} ${headers} tnl-functions-test.cu )
    TARGET_LINK_LIBRARIES( tnl-functions-test${mpiExt}${debugExt} ${CPPUNIT_LIBRARIES}
                                                              tnl${mpiExt}${debugExt}-0.1 )
else()
    ADD_EXECUTABLE( tnl-functions-test${mpiExt}${debugExt} ${headers} tnl-functions-test.cpp )
    TARGET_LINK_LIBRARIES( tnl-functions-test${mpiExt}${debugExt} ${CPPUNIT_LIBRARIES}
                                                              tnl${mpiExt}${debugExt}-0.1 )
endif()

                                                             
INSTALL( FILES tnl-functions-test
         DESTINATION share/tnl-${tnlVersion} )                                                            
                                                                       
Loading