From 8b4280d31c30d0dc4b04cc00622a6787ffd23e41 Mon Sep 17 00:00:00 2001
From: Tomas Oberhuber <tomas.oberhuber@fjfi.cvut.cz>
Date: Fri, 10 May 2013 15:09:36 +0200
Subject: [PATCH] Finished 2D diffusion on deformed grids. Template
specialization of numerical schemes for the grids with or without
deformations.
---
examples/navier-stokes/Makefile | 2 +-
src/implementation/mesh/tnlGrid2D_impl.h | 8 +-
.../mesh/tnlIdenticalGridGeometry_impl.h | 2 +-
.../diffusion/tnlLinearDiffusion_impl.h | 182 ++++++++++---
.../schemes/euler/fvm/tnlLaxFridrichs_impl.h | 245 ++++++++++++++----
.../gradient/tnlCentralFDMGradient_impl.h | 74 +++++-
src/schemes/diffusion/tnlLinearDiffusion.h | 43 ++-
src/schemes/euler/fvm/tnlLaxFridrichs.h | 59 ++++-
src/schemes/gradient/tnlCentralFDMGradient.h | 49 +++-
9 files changed, 537 insertions(+), 127 deletions(-)
diff --git a/examples/navier-stokes/Makefile b/examples/navier-stokes/Makefile
index d25b1b9779..de16a7dbc3 100644
--- a/examples/navier-stokes/Makefile
+++ b/examples/navier-stokes/Makefile
@@ -8,7 +8,7 @@ CONFIG_FILE = $(TARGET).cfg.desc
INSTALL_DIR = ${HOME}/local
CXX = g++
CUDA_CXX = nvcc
-OMP_FLAGS = -DHAVE_OPENMP -fopenmp
+#OMP_FLAGS = -DHAVE_OPENMP -fopenmp
#CXX_FLAGS = -std=gnu++0x -I$(TNL_INCLUDE_DIR) -O0 -g -DDEBUG $(OMP_FLAGS) -DTEMPLATE_EXPLICIT_INSTANTIATION
CXX_FLAGS = -std=gnu++0x -I$(TNL_INCLUDE_DIR) -O3 $(OMP_FLAGS) -DTEMPLATE_EXPLICIT_INSTANTIATION
#CXX_FLAGS = -std=gnu++0x -I$(TNL_INCLUDE_DIR) -O3 $(OMP_FLAGS) -pg -DTEMPLATE_EXPLICIT_INSTANTIATION
diff --git a/src/implementation/mesh/tnlGrid2D_impl.h b/src/implementation/mesh/tnlGrid2D_impl.h
index 54f32cc510..fc8ce76308 100644
--- a/src/implementation/mesh/tnlGrid2D_impl.h
+++ b/src/implementation/mesh/tnlGrid2D_impl.h
@@ -391,16 +391,14 @@ bool tnlGrid< 2, Real, Device, Index, Geometry > :: write( const MeshFunction& f
cerr << "I am not able to open the file " << fileName << "." << endl;
return false;
}
- const RealType hx = getParametricStep(). x();
- const RealType hy = getParametricStep(). y();
if( format == "gnuplot" )
for( IndexType j = 0; j < getDimensions(). y(); j++ )
{
for( IndexType i = 0; i < getDimensions(). x(); i++ )
{
- const RealType x = this -> getOrigin(). x() + i * hx;
- const RealType y = this -> getOrigin(). y() + j * hy;
- file << x << " " << " " << y << " " << function[ this -> getElementIndex( i, j ) ] << endl;
+ VertexType v;
+ this -> getElementCenter( CoordinatesType( i, j ), v );
+ file << v. x() << " " << " " << v. y() << " " << function[ this -> getElementIndex( i, j ) ] << endl;
}
file << endl;
}
diff --git a/src/implementation/mesh/tnlIdenticalGridGeometry_impl.h b/src/implementation/mesh/tnlIdenticalGridGeometry_impl.h
index 20cc9c8645..f288d9edbe 100644
--- a/src/implementation/mesh/tnlIdenticalGridGeometry_impl.h
+++ b/src/implementation/mesh/tnlIdenticalGridGeometry_impl.h
@@ -169,7 +169,7 @@ template< typename Real,
template< int dx, int dy >
Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementCoVolumeMeasure( const CoordinatesType& coordinates ) const
{
- return elementMeasure;
+ return 0.5 * elementMeasure;
}
diff --git a/src/implementation/schemes/diffusion/tnlLinearDiffusion_impl.h b/src/implementation/schemes/diffusion/tnlLinearDiffusion_impl.h
index 4c4e5abaec..06bff8b96d 100644
--- a/src/implementation/schemes/diffusion/tnlLinearDiffusion_impl.h
+++ b/src/implementation/schemes/diffusion/tnlLinearDiffusion_impl.h
@@ -18,52 +18,55 @@
#ifndef TNLLINEARDIFFUSION_IMPL_H_
#define TNLLINEARDIFFUSION_IMPL_H_
-template< typename Real, typename Device, typename Index >
-tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index > > :: tnlLinearDiffusion()
+template< typename Real,
+ typename Device,
+ typename Index,
+ template< int, typename, typename, typename > class GridGeometry >
+tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index, GridGeometry > > :: tnlLinearDiffusion()
: mesh( 0 )
{
}
-template< typename Real, typename Device, typename Index >
-void tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index > > :: bindMesh( const tnlGrid< 2, RealType, DeviceType, IndexType >& mesh )
+template< typename Real,
+ typename Device,
+ typename Index,
+ template< int, typename, typename, typename > class GridGeometry >
+void tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index, GridGeometry > > :: bindMesh( const MeshType& mesh )
{
this -> mesh = &mesh;
}
-template< typename Real, typename Device, typename Index >
+template< typename Real,
+ typename Device,
+ typename Index,
+ template< int, typename, typename, typename > class GridGeometry >
template< typename Vector >
-void tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index > > :: setFunction( Vector& f )
+void tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index, GridGeometry > > :: setFunction( Vector& f )
{
this -> f. bind( f );
this -> f. setName( tnlString( "bind Of " ) + f. getName() );
}
-template< typename Real, typename Device, typename Index >
-Real tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index > > :: getDiffusion( const Index& c ) const
+template< typename Real,
+ typename Device,
+ typename Index,
+ template< int, typename, typename, typename > class GridGeometry >
+Real tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index, GridGeometry > > :: getDiffusion( const Index& c ) const
{
tnlAssert( this -> mesh, cerr << "No mesh was set in tnlLinearDiffusion. Use the bindMesh method." );
- const Real hx = mesh -> getParametricStep(). x();
- const Real hy = mesh -> getParametricStep(). y();
-
const Index e = mesh -> getElementNeighbour( c, 1, 0 );
const Index w = mesh -> getElementNeighbour( c, -1, 0 );
const Index n = mesh -> getElementNeighbour( c, 0, 1 );
const Index s = mesh -> getElementNeighbour( c, 0, -1 );
- //return ( f[ e ] - 2.0 * f[ c ] + f[ w ] ) / ( hx * hx ) +
- // ( f[ n ] - 2.0 * f[ c ] + f[ s ] ) / ( hy * hy );
-
- /****
- * Diffusion for the deformed grids
- */
const Index en = mesh -> getElementNeighbour( c, 1, 1 );
const Index es = mesh -> getElementNeighbour( c, 1, -1 );
const Index wn = mesh -> getElementNeighbour( c, -1, 1 );
const Index ws = mesh -> getElementNeighbour( c, -1, -1 );
CoordinatesType cCoordinates, eCoordinates, nCoordinates, wCoordinates, sCoordinates,
- enCoordinates, esCoordinates, wnCoordinates, wsCoordinates;
+ enCoordinates, esCoordinates, wnCoordinates, wsCoordinates;
this -> mesh -> getElementCoordinates( c, cCoordinates );
eCoordinates = nCoordinates = wCoordinates = sCoordinates = cCoordinates;
eCoordinates. x() ++;
@@ -119,54 +122,159 @@ Real tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index > > :: getDiffusion( co
this -> mesh -> template getVertex< -1, -1 >( cCoordinates, wsVertex );
const RealType f_x_e = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< 1, 0 >( cCoordinates ) *
- ( f_cen * ( enVertex. y() - cCenter. y() ) +
+ ( f_cen * ( cCenter. y() - enVertex. y() ) +
f_ces * ( esVertex. y() - cCenter. y() ) +
- 0.5 * ( f_es + f[ e ] ) * ( esVertex. y() - eCenter. y() ) +
+ 0.5 * ( f_es + f[ e ] ) * ( eCenter. y() - esVertex. y() ) +
0.5 * ( f_en + f[ e ] ) * ( enVertex. y() - eCenter. y() ) );
const RealType f_y_e = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< 1, 0 >( cCoordinates ) *
( f_cen * ( enVertex. x() - cCenter. x() ) +
- f_ces * ( esVertex. x() - cCenter. x() ) +
+ f_ces * ( cCenter. x() - esVertex. x() ) +
0.5 * ( f_es + f[ e ] ) * ( esVertex. x() - eCenter. x() ) +
- 0.5 * ( f_en + f[ e ] ) * ( enVertex. x() - eCenter. x() ) );
+ 0.5 * ( f_en + f[ e ] ) * ( eCenter. x() ) - enVertex. x() );
+
const RealType f_x_w = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< -1, 0 >( cCoordinates ) *
( f_cwn * ( wnVertex. y() - cCenter. y() ) +
- f_cws * ( wsVertex. y() - cCenter. y() ) +
+ f_cws * ( cCenter. y() - wsVertex. y() ) +
0.5 * ( f_ws + f[ w ] ) * ( wsVertex. y() - wCenter. y() ) +
- 0.5 * ( f_wn + f[ w ] ) * ( wnVertex. y() - wCenter. y() ) );
+ 0.5 * ( f_wn + f[ w ] ) * ( wCenter. y() - wnVertex. y() ) );
const RealType f_y_w = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< -1, 0 >( cCoordinates ) *
- ( f_cwn * ( wnVertex. x() - cCenter. x() ) +
+ ( f_cwn * ( cCenter. x() - wnVertex. x() ) +
f_cws * ( wsVertex. x() - cCenter. x() ) +
- 0.5 * ( f_ws + f[ w ] ) * ( wsVertex. x() - wCenter. x() ) +
+ 0.5 * ( f_ws + f[ w ] ) * ( wCenter. x() - wsVertex. x() ) +
0.5 * ( f_wn + f[ w ] ) * ( wnVertex. x() - wCenter. x() ) );
+
const RealType f_x_n = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< 0, 1 >( cCoordinates ) *
( f_cen * ( enVertex. y() - cCenter. y() ) +
- f_cwn * ( wnVertex. y() - cCenter. y() ) +
- 0.5 * ( f_en + f[ n ] ) * ( enVertex. y() - nCenter. y() ) +
+ f_cwn * ( cCenter. y() - wnVertex. y() ) +
+ 0.5 * ( f_en + f[ n ] ) * ( nCenter. y() - enVertex. y() ) +
0.5 * ( f_wn + f[ n ] ) * ( wnVertex. y() - nCenter. y() ) );
const RealType f_y_n = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< 0, 1 >( cCoordinates ) *
- ( f_cen * ( enVertex. x() - cCenter. x() ) +
+ ( f_cen * ( cCenter. x() - enVertex. x() ) +
f_cwn * ( wnVertex. x() - cCenter. x() ) +
0.5 * ( f_en + f[ n ] ) * ( enVertex. x() - nCenter. x() ) +
- 0.5 * ( f_wn + f[ n ] ) * ( wnVertex. x() - nCenter. x() ) );
+ 0.5 * ( f_wn + f[ n ] ) * ( nCenter. x() - wnVertex. x() ) );
+
const RealType f_x_s = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< 0, -1 >( cCoordinates ) *
- ( f_ces * ( esVertex. y() - cCenter. y() ) +
+ ( f_ces * ( cCenter. y() - esVertex. y() ) +
f_cws * ( wsVertex. y() - cCenter. y() ) +
0.5 * ( f_es + f[ s ] ) * ( esVertex. y() - sCenter. y() ) +
- 0.5 * ( f_ws + f[ s ] ) * ( wsVertex. y() - sCenter. y() ) );
+ 0.5 * ( f_ws + f[ s ] ) * ( sCenter. y() - wsVertex. y() ) );
const RealType f_y_s = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< 0, -1 >( cCoordinates ) *
( f_ces * ( esVertex. x() - cCenter. x() ) +
- f_cws * ( wsVertex. x() - cCenter. x() ) +
- 0.5 * ( f_es + f[ s ] ) * ( enVertex. x() - sCenter. x() ) +
- 0.5 * ( f_ws + f[ s ] ) * ( wnVertex. x() - sCenter. x() ) );
+ f_cws * ( cCenter. x() - wsVertex. x() ) +
+ 0.5 * ( f_es + f[ s ] ) * ( sCenter. x() - esVertex. x() ) +
+ 0.5 * ( f_ws + f[ s ] ) * ( wsVertex. x() - sCenter. x() ) );
VertexType eNormal, nNormal, wNormal, sNormal;
this -> mesh -> template getEdgeNormal< 1, 0 >( cCoordinates, eNormal );
this -> mesh -> template getEdgeNormal< -1, 0 >( cCoordinates, wNormal );
this -> mesh -> template getEdgeNormal< 0, 1 >( cCoordinates, nNormal );
this -> mesh -> template getEdgeNormal< 0, -1 >( cCoordinates, sNormal );
+
+ /*const RealType eps = 0.000001;
+ if( fabs( f_x_e - ( f[ e ] - f[ c ] ) / hx ) > eps ||
+ fabs( f_x_w - ( f[ c ] - f[ w ] ) / hx ) > eps ||
+ fabs( f_y_n - ( f[ n ] - f[ c ] ) / hy ) > eps ||
+ fabs( f_y_s - ( f[ c ] - f[ s ] ) / hy ) > eps )
+ {
+
+ cout << "cCoordinates = " << cCoordinates << endl;
+ cout << "this -> mesh -> template getElementCoVolumeMeasure< 1, 0 >( cCoordinates ) = " <<
+ this -> mesh -> template getElementCoVolumeMeasure< 1, 0 >( cCoordinates ) << endl;
+ cout << "enVertex. y() - cCenter. y() = " << enVertex. y() - cCenter. y() << endl;
+ cout << "esVertex. y() - cCenter. y() = " << esVertex. y() - cCenter. y() << endl;
+ cout << "enVertex. y() - eCenter. y() = " << enVertex. y() - eCenter. y() << endl;
+ cout << "esVertex. y() - eCenter. y() = " << esVertex. y() - eCenter. y() << endl;
+ cout << "cMeasure = " << cMeasure << endl;
+ cout << "eMeasure = " << eMeasure << endl;
+ cout << "nMeasure = " << nMeasure << endl;
+ cout << "wMeasure = " << wMeasure << endl;
+ cout << "sMeasure = " << sMeasure << endl;
+ cout << "f[ c ] = " << f[ c ] << endl;
+ cout << "f[ e ] = " << f[ e ] << endl;
+ cout << "f[ n ] = " << f[ n ] << endl;
+ cout << "f[ w ] = " << f[ w ] << endl;
+ cout << "f[ s ] = " << f[ s ] << endl;
+ cout << "f[ en ] = " << f[ en ] << endl;
+ cout << "f[ es ] = " << f[ es ] << endl;
+ cout << "f[ wn ] = " << f[ wn ] << endl;
+ cout << "f[ ws ] = " << f[ ws ] << endl;
+ cout << "f_en = " << f_en << endl;
+ cout << "f_es = " << f_es << endl;
+ cout << "f_wn = " << f_wn << endl;
+ cout << "f_ws = " << f_ws << endl;
+ cout << "f_cen = " << f_cen << endl;
+ cout << "f_ces = " << f_ces << endl;
+ cout << "f_cwn = " << f_cwn << endl;
+ cout << "f_cws = " << f_cws << endl;
+ cout << "0.5 * ( f_es + f[ e ] ) = " << 0.5 * ( f_es + f[ e ] ) << endl;
+ cout << "0.5 * ( f_en + f[ e ] ) = " << 0.5 * ( f_en + f[ e ] ) << endl;
+
+ cout << "( f[ e ] - f[ c ] ) / hx = " << ( f[ e ] - f[ c ] ) / hx << endl;
+ cout << "( f[ c ] - f[ w ] ) / hx = " << ( f[ c ] - f[ w ] ) / hx << endl;
+ cout << "( f[ n ] - f[ c ] ) / hy = " << ( f[ n ] - f[ c ] ) / hy << endl;
+ cout << "( f[ c ] - f[ s ] ) / hy = " << ( f[ c ] - f[ s ] ) / hy << endl;
+ cout << " f_x_e = " << f_x_e << endl;
+ cout << " f_x_w = " << f_x_w << endl;
+ cout << " f_y_n = " << f_y_n << endl;
+ cout << " f_y_s = " << f_y_s << endl;
+
+ getchar();
+ }*/
+
+ //return ( ( f[ e ] - f[ c ] ) / hx - ( f[ c ] - f[ w ] ) / hx ) / hx +
+ // ( ( f[ n ] - f[ c ] ) / hy - ( f[ c ] - f[ s ] ) / hy ) / hy;
return 1.0 / this -> mesh -> getElementMeasure( cCoordinates ) *
- ( f_x_e * eNormal. x() + f_x_w * wNormal. x() +
- f_y_n * nNormal. y() + f_y_s * sNormal. y() );
+ ( f_x_e * eNormal. x() +
+ f_x_n * nNormal. x() +
+ f_x_w * wNormal. x() +
+ f_x_s * sNormal. x() +
+ f_y_e * eNormal. y() +
+ f_y_n * nNormal. y() +
+ f_y_w * wNormal. y() +
+ f_y_s * sNormal. y() );
}
+/****
+ * Specialization for the grids with no deformations (Identical grid geometry)
+ */
+
+template< typename Real, typename Device, typename Index >
+tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry > > :: tnlLinearDiffusion()
+: mesh( 0 )
+{
+}
+
+template< typename Real, typename Device, typename Index >
+void tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry > > :: bindMesh( const tnlGrid< 2, RealType, DeviceType, IndexType, tnlIdenticalGridGeometry >& mesh )
+{
+ this -> mesh = &mesh;
+}
+
+template< typename Real, typename Device, typename Index >
+ template< typename Vector >
+void tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry > > :: setFunction( Vector& f )
+{
+ this -> f. bind( f );
+ this -> f. setName( tnlString( "bind Of " ) + f. getName() );
+}
+
+template< typename Real, typename Device, typename Index >
+Real tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry > > :: getDiffusion( const Index& c ) const
+{
+ tnlAssert( this -> mesh, cerr << "No mesh was set in tnlLinearDiffusion. Use the bindMesh method." );
+
+ const Real hx = mesh -> getParametricStep(). x();
+ const Real hy = mesh -> getParametricStep(). y();
+
+ const Index e = mesh -> getElementNeighbour( c, 1, 0 );
+ const Index w = mesh -> getElementNeighbour( c, -1, 0 );
+ const Index n = mesh -> getElementNeighbour( c, 0, 1 );
+ const Index s = mesh -> getElementNeighbour( c, 0, -1 );
+
+ return ( f[ e ] - 2.0 * f[ c ] + f[ w ] ) / ( hx * hx ) +
+ ( f[ n ] - 2.0 * f[ c ] + f[ s ] ) / ( hy * hy );
+}
+
+
#endif
diff --git a/src/implementation/schemes/euler/fvm/tnlLaxFridrichs_impl.h b/src/implementation/schemes/euler/fvm/tnlLaxFridrichs_impl.h
index 4352d01236..dcbd6dad28 100644
--- a/src/implementation/schemes/euler/fvm/tnlLaxFridrichs_impl.h
+++ b/src/implementation/schemes/euler/fvm/tnlLaxFridrichs_impl.h
@@ -22,8 +22,9 @@
template< typename Real,
typename Device,
typename Index,
- typename PressureGradient >
-tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >,
+ typename PressureGradient,
+ template< int, typename, typename, typename > class GridGeometry >
+tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index, GridGeometry >,
PressureGradient > :: tnlLaxFridrichs()
: regularizeEps( 0.0 ),
viscosityCoefficient( 1.0 ),
@@ -35,8 +36,9 @@ tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >,
template< typename Real,
typename Device,
typename Index,
- typename PressureGradient >
-void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient > :: bindMesh( const MeshType& mesh )
+ typename PressureGradient,
+ template< int, typename, typename, typename > class GridGeometry >
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index, GridGeometry >, PressureGradient > :: bindMesh( const MeshType& mesh )
{
this -> mesh = &mesh;
}
@@ -44,8 +46,9 @@ void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient > ::
template< typename Real,
typename Device,
typename Index,
- typename PressureGradient >
-void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient > :: setRegularization( const RealType& epsilon )
+ typename PressureGradient,
+ template< int, typename, typename, typename > class GridGeometry >
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index, GridGeometry >, PressureGradient > :: setRegularization( const RealType& epsilon )
{
this -> regularizeEps = epsilon;
}
@@ -53,8 +56,9 @@ void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient > ::
template< typename Real,
typename Device,
typename Index,
- typename PressureGradient >
-void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient > :: setViscosityCoefficient( const RealType& v )
+ typename PressureGradient,
+ template< int, typename, typename, typename > class GridGeometry >
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index, GridGeometry >, PressureGradient > :: setViscosityCoefficient( const RealType& v )
{
this -> viscosityCoefficient = v;
}
@@ -62,9 +66,10 @@ void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient > ::
template< typename Real,
typename Device,
typename Index,
- typename PressureGradient >
+ typename PressureGradient,
+ template< int, typename, typename, typename > class GridGeometry >
template< typename Vector >
-void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient > :: setRho( Vector& rho )
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index, GridGeometry >, PressureGradient > :: setRho( Vector& rho )
{
this -> rho. bind( rho );
this -> rho. setName( tnlString( "bind Of " ) + rho. getName() );
@@ -73,9 +78,10 @@ void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient > ::
template< typename Real,
typename Device,
typename Index,
- typename PressureGradient >
+ typename PressureGradient,
+ template< int, typename, typename, typename > class GridGeometry >
template< typename Vector >
-void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient > :: setRhoU1( Vector& rho_u1 )
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index, GridGeometry >, PressureGradient > :: setRhoU1( Vector& rho_u1 )
{
this -> rho_u1. bind( rho_u1 );
this -> rho_u1. setName( tnlString( "bind Of " ) + rho_u1. getName() );
@@ -84,9 +90,10 @@ void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient > ::
template< typename Real,
typename Device,
typename Index,
- typename PressureGradient >
+ typename PressureGradient,
+ template< int, typename, typename, typename > class GridGeometry >
template< typename Vector >
-void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient > :: setRhoU2( Vector& rho_u2 )
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index, GridGeometry >, PressureGradient > :: setRhoU2( Vector& rho_u2 )
{
this -> rho_u2. bind( rho_u2 );
this -> rho_u2. setName( tnlString( "bind Of " ) + rho_u2. getName() );
@@ -95,9 +102,10 @@ void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient > ::
template< typename Real,
typename Device,
typename Index,
- typename PressureGradient >
+ typename PressureGradient,
+ template< int, typename, typename, typename > class GridGeometry >
template< typename Vector >
-void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient > :: setPressureGradient( Vector& grad_p )
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index, GridGeometry >, PressureGradient > :: setPressureGradient( Vector& grad_p )
{
this -> pressureGradient = &grad_p;
}
@@ -105,8 +113,9 @@ void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient > ::
template< typename Real,
typename Device,
typename Index,
- typename PressureGradient >
-void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient > :: getExplicitRhs( const IndexType centralVolume,
+ typename PressureGradient,
+ template< int, typename, typename, typename > class GridGeometry >
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index, GridGeometry >, PressureGradient > :: getExplicitRhs( const IndexType centralVolume,
RealType& rho_t,
RealType& rho_u1_t,
RealType& rho_u2_t ) const
@@ -114,53 +123,16 @@ void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient > ::
tnlAssert( mesh, cerr << "No mesh has been binded with the Lax-Fridrichs scheme." );
tnlAssert( pressureGradient, cerr << "No pressure gradient was set in the the Lax-Fridrichs scheme." )
- const IndexType& xSize = this -> mesh -> getDimensions(). x();
- const IndexType& ySize = this -> mesh -> getDimensions(). y();
- const RealType hx = this -> mesh -> getParametricStep(). x();
- const RealType hy = this -> mesh -> getParametricStep(). y();
-
const IndexType& c = centralVolume;
const IndexType e = this -> mesh -> getElementNeighbour( centralVolume, 1, 0 );
const IndexType w = this -> mesh -> getElementNeighbour( centralVolume, -1, 0 );
const IndexType n = this -> mesh -> getElementNeighbour( centralVolume, 0, 1 );
const IndexType s = this -> mesh -> getElementNeighbour( centralVolume, 0, -1 );
- /****
- * rho_t + ( rho u_1 )_x + ( rho u_2 )_y = 0
- */
const RealType u1_e = rho_u1[ e ] / regularize( rho[ e ] );
const RealType u1_w = rho_u1[ w ] / regularize( rho[ w ] );
const RealType u2_n = rho_u2[ n ] / regularize( rho[ n ] );
const RealType u2_s = rho_u2[ s ] / regularize( rho[ s ] );
- rho_t = this -> viscosityCoefficient * 0.25 * ( rho[ e ] + rho[ w ] + rho[ s ] + rho[ n ] - 4.0 * rho[ c ] )
- - ( rho[ e ] * u1_e - rho[ w ] * u1_w ) / ( 2.0 * hx )
- - ( rho[ n ] * u2_n - rho[ s ] * u2_s ) / ( 2.0 * hy );
-
- /****
- * Compute the pressure gradient
- */
- VertexType grad_p;
- pressureGradient -> getGradient( c, grad_p );
-
- /****
- * ( rho * u1 )_t + ( rho * u1 * u1 )_x + ( rho * u1 * u2 )_y - p_x = 0
- */
- rho_u1_t = this -> viscosityCoefficient * 0.25 * ( rho_u1[ e ] + rho_u1[ w ] + rho_u1[ s ] + rho_u1[ n ] - 4.0 * rho_u1[ c ] )
- - ( rho_u1[ e ] * u1_e - rho_u1[ w ] * u1_w ) / ( 2.0 * hx )
- - ( rho_u1[ n ] * u2_n - rho_u1[ s ] * u2_s ) / ( 2.0 * hy )
- - grad_p. x();
- /****
- * ( rho * u2 )_t + ( rho * u2 * u1 )_x + ( rho * u2 * u2 )_y - p_y = 0
- */
- rho_u2_t = this -> viscosityCoefficient * 0.25 * ( rho_u2[ e ] + rho_u2[ w ] + rho_u2[ s ] + rho_u2[ n ] - 4.0 * rho_u2[ c ] )
- - ( rho_u2[ e ] * u1_e - rho_u2[ w ] * u1_w ) / ( 2.0 * hx )
- - ( rho_u2[ n ] * u2_n - rho_u2[ s ] * u2_s ) / ( 2.0 * hy )
- - grad_p. y();
-
-
- /****
- * Scheme for deformed grids
- */
const RealType u1_c = rho_u1[ c ] / regularize( rho[ c ] );
const RealType u1_n = rho_u1[ n ] / regularize( rho[ n ] );
const RealType u1_s = rho_u1[ s ] / regularize( rho[ s ] );
@@ -168,7 +140,6 @@ void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient > ::
const RealType u2_e = rho_u2[ e ] / regularize( rho[ e ] );
const RealType u2_w = rho_u2[ w ] / regularize( rho[ w ] );
-
/****
* Get the central volume and its neighbours (east, north, west, south) coordinates
*/
@@ -228,6 +199,12 @@ void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient > ::
const RealType rho_u2_g_n = 0.5 * ( rho_u2[ c ] * u2_c + rho_u2[ n ] * u2_n );
const RealType rho_u2_g_s = 0.5 * ( rho_u2[ c ] * u2_c + rho_u2[ s ] * u2_s );
+ /****
+ * Compute the pressure gradient
+ */
+ VertexType grad_p;
+ pressureGradient -> getGradient( c, grad_p );
+
/****
* rho_t + ( rho u_1 )_x + ( rho u_2 )_y = 0
*/
@@ -270,13 +247,167 @@ void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient > ::
- grad_p. y();
}
+template< typename Real,
+ typename Device,
+ typename Index,
+ typename PressureGradient,
+ template< int, typename, typename, typename > class GridGeometry >
+Real tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index, GridGeometry >, PressureGradient > :: regularize( const Real& r ) const
+{
+ return r + ( ( r >= 0 ) - ( r < 0 ) ) * this -> regularizeEps;
+}
+
+/****
+ * Specialization for the grids with no deformations (Identical grid geometry)
+ */
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ typename PressureGradient >
+tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry >,
+ PressureGradient > :: tnlLaxFridrichs()
+: regularizeEps( 0.0 ),
+ viscosityCoefficient( 1.0 ),
+ mesh( 0 ),
+ pressureGradient( 0 )
+{
+}
+
template< typename Real,
typename Device,
typename Index,
typename PressureGradient >
-Real tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient > :: regularize( const Real& r ) const
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry >, PressureGradient > :: bindMesh( const MeshType& mesh )
+{
+ this -> mesh = &mesh;
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ typename PressureGradient >
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry >, PressureGradient > :: setRegularization( const RealType& epsilon )
+{
+ this -> regularizeEps = epsilon;
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ typename PressureGradient >
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry >, PressureGradient > :: setViscosityCoefficient( const RealType& v )
+{
+ this -> viscosityCoefficient = v;
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ typename PressureGradient >
+ template< typename Vector >
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry >, PressureGradient > :: setRho( Vector& rho )
+{
+ this -> rho. bind( rho );
+ this -> rho. setName( tnlString( "bind Of " ) + rho. getName() );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ typename PressureGradient >
+ template< typename Vector >
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry >, PressureGradient > :: setRhoU1( Vector& rho_u1 )
+{
+ this -> rho_u1. bind( rho_u1 );
+ this -> rho_u1. setName( tnlString( "bind Of " ) + rho_u1. getName() );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ typename PressureGradient >
+ template< typename Vector >
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry >, PressureGradient > :: setRhoU2( Vector& rho_u2 )
+{
+ this -> rho_u2. bind( rho_u2 );
+ this -> rho_u2. setName( tnlString( "bind Of " ) + rho_u2. getName() );
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ typename PressureGradient >
+ template< typename Vector >
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry >, PressureGradient > :: setPressureGradient( Vector& grad_p )
+{
+ this -> pressureGradient = &grad_p;
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ typename PressureGradient >
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry >, PressureGradient > :: getExplicitRhs( const IndexType centralVolume,
+ RealType& rho_t,
+ RealType& rho_u1_t,
+ RealType& rho_u2_t ) const
+{
+ tnlAssert( mesh, cerr << "No mesh has been binded with the Lax-Fridrichs scheme." );
+ tnlAssert( pressureGradient, cerr << "No pressure gradient was set in the the Lax-Fridrichs scheme." )
+
+ const IndexType& xSize = this -> mesh -> getDimensions(). x();
+ const IndexType& ySize = this -> mesh -> getDimensions(). y();
+ const RealType hx = this -> mesh -> getParametricStep(). x();
+ const RealType hy = this -> mesh -> getParametricStep(). y();
+
+ const IndexType& c = centralVolume;
+ const IndexType e = this -> mesh -> getElementNeighbour( centralVolume, 1, 0 );
+ const IndexType w = this -> mesh -> getElementNeighbour( centralVolume, -1, 0 );
+ const IndexType n = this -> mesh -> getElementNeighbour( centralVolume, 0, 1 );
+ const IndexType s = this -> mesh -> getElementNeighbour( centralVolume, 0, -1 );
+
+ /****
+ * rho_t + ( rho u_1 )_x + ( rho u_2 )_y = 0
+ */
+ const RealType u1_e = rho_u1[ e ] / regularize( rho[ e ] );
+ const RealType u1_w = rho_u1[ w ] / regularize( rho[ w ] );
+ const RealType u2_n = rho_u2[ n ] / regularize( rho[ n ] );
+ const RealType u2_s = rho_u2[ s ] / regularize( rho[ s ] );
+ rho_t = this -> viscosityCoefficient * 0.25 * ( rho[ e ] + rho[ w ] + rho[ s ] + rho[ n ] - 4.0 * rho[ c ] )
+ - ( rho[ e ] * u1_e - rho[ w ] * u1_w ) / ( 2.0 * hx )
+ - ( rho[ n ] * u2_n - rho[ s ] * u2_s ) / ( 2.0 * hy );
+
+ /****
+ * Compute the pressure gradient
+ */
+ VertexType grad_p;
+ pressureGradient -> getGradient( c, grad_p );
+
+ /****
+ * ( rho * u1 )_t + ( rho * u1 * u1 )_x + ( rho * u1 * u2 )_y - p_x = 0
+ */
+ rho_u1_t = this -> viscosityCoefficient * 0.25 * ( rho_u1[ e ] + rho_u1[ w ] + rho_u1[ s ] + rho_u1[ n ] - 4.0 * rho_u1[ c ] )
+ - ( rho_u1[ e ] * u1_e - rho_u1[ w ] * u1_w ) / ( 2.0 * hx )
+ - ( rho_u1[ n ] * u2_n - rho_u1[ s ] * u2_s ) / ( 2.0 * hy )
+ - grad_p. x();
+ /****
+ * ( rho * u2 )_t + ( rho * u2 * u1 )_x + ( rho * u2 * u2 )_y - p_y = 0
+ */
+ rho_u2_t = this -> viscosityCoefficient * 0.25 * ( rho_u2[ e ] + rho_u2[ w ] + rho_u2[ s ] + rho_u2[ n ] - 4.0 * rho_u2[ c ] )
+ - ( rho_u2[ e ] * u1_e - rho_u2[ w ] * u1_w ) / ( 2.0 * hx )
+ - ( rho_u2[ n ] * u2_n - rho_u2[ s ] * u2_s ) / ( 2.0 * hy )
+ - grad_p. y();
+}
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ typename PressureGradient >
+Real tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry >, PressureGradient > :: regularize( const Real& r ) const
{
return r + ( ( r >= 0 ) - ( r < 0 ) ) * this -> regularizeEps;
}
+
#endif
diff --git a/src/implementation/schemes/gradient/tnlCentralFDMGradient_impl.h b/src/implementation/schemes/gradient/tnlCentralFDMGradient_impl.h
index 558c1f24db..cc68ec6a2e 100644
--- a/src/implementation/schemes/gradient/tnlCentralFDMGradient_impl.h
+++ b/src/implementation/schemes/gradient/tnlCentralFDMGradient_impl.h
@@ -18,29 +18,41 @@
#ifndef TNLCENTRALFDMGRADIENT_IMPL_H_
#define TNLCENTRALFDMGRADIENT_IMPL_H_
-template< typename Real, typename Device, typename Index >
-tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index > > :: tnlCentralFDMGradient()
+template< typename Real,
+ typename Device,
+ typename Index,
+ template< int, typename, typename, typename > class GridGeometry >
+tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index, GridGeometry > > :: tnlCentralFDMGradient()
: mesh( 0 )
{
}
-template< typename Real, typename Device, typename Index >
-void tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index > > :: bindMesh( const tnlGrid< 2, RealType, DeviceType, IndexType >& mesh )
+template< typename Real,
+ typename Device,
+ typename Index,
+ template< int, typename, typename, typename > class GridGeometry >
+void tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index, GridGeometry > > :: bindMesh( const MeshType& mesh )
{
this -> mesh = &mesh;
}
-template< typename Real, typename Device, typename Index >
+template< typename Real,
+ typename Device,
+ typename Index,
+ template< int, typename, typename, typename > class GridGeometry >
template< typename Vector >
-void tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index > > :: setFunction( Vector& f )
+void tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index, GridGeometry > > :: setFunction( Vector& f )
{
this -> f. bind( f );
this -> f. setName( tnlString( "bind Of " ) + f. getName() );
}
-template< typename Real, typename Device, typename Index >
-void tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index > > :: getGradient( const Index& i,
- VertexType& grad_f ) const
+template< typename Real,
+ typename Device,
+ typename Index,
+ template< int, typename, typename, typename > class GridGeometry >
+void tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index, GridGeometry > > :: getGradient( const Index& i,
+ VertexType& grad_f ) const
{
tnlAssert( this -> mesh, cerr << "No mesh was set in tnlCentralFDMGradient. Use the bindMesh method." );
@@ -59,12 +71,48 @@ void tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index > > :: getGradient(
nCoordinates. y() ++;
sCoordinates. y() --;
-
- //grad_f. x() = ( f[ e ] - f[ w ] ) / ( 2.0 * mesh -> getParametricStep(). x() );
- //grad_f. y() = ( f[ n ] - f[ s ] ) / ( 2.0 * mesh -> getParametricStep(). y() );
-
grad_f. x() = ( f[ e ] - f[ w ] ) / ( mesh -> getElementsDistance( eCoordinates, wCoordinates ) );
grad_f. y() = ( f[ n ] - f[ s ] ) / ( mesh -> getElementsDistance( nCoordinates, sCoordinates ) );
}
+/****
+ * Specialization for the grids with no deformations (Identical grid geometry)
+ */
+
+template< typename Real, typename Device, typename Index >
+tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry > > :: tnlCentralFDMGradient()
+: mesh( 0 )
+{
+}
+
+template< typename Real, typename Device, typename Index >
+void tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry > > :: bindMesh( const MeshType& mesh )
+{
+ this -> mesh = &mesh;
+}
+
+template< typename Real, typename Device, typename Index >
+ template< typename Vector >
+void tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry > > :: setFunction( Vector& f )
+{
+ this -> f. bind( f );
+ this -> f. setName( tnlString( "bind Of " ) + f. getName() );
+}
+
+template< typename Real, typename Device, typename Index >
+void tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry > > :: getGradient( const Index& i,
+ VertexType& grad_f ) const
+{
+ tnlAssert( this -> mesh, cerr << "No mesh was set in tnlCentralFDMGradient. Use the bindMesh method." );
+
+ const Index e = mesh -> getElementNeighbour( i, 1, 0 );
+ const Index w = mesh -> getElementNeighbour( i, -1, 0 );
+ const Index n = mesh -> getElementNeighbour( i, 0, 1 );
+ const Index s = mesh -> getElementNeighbour( i, 0, -1 );
+
+ grad_f. x() = ( f[ e ] - f[ w ] ) / ( 2.0 * mesh -> getParametricStep(). x() );
+ grad_f. y() = ( f[ n ] - f[ s ] ) / ( 2.0 * mesh -> getParametricStep(). y() );
+}
+
+
#endif
diff --git a/src/schemes/diffusion/tnlLinearDiffusion.h b/src/schemes/diffusion/tnlLinearDiffusion.h
index f4589984c6..0c1a0d00d0 100644
--- a/src/schemes/diffusion/tnlLinearDiffusion.h
+++ b/src/schemes/diffusion/tnlLinearDiffusion.h
@@ -19,6 +19,7 @@
#define TNLLINEARDIFFUSION_H_
#include <mesh/tnlGrid.h>
+#include <mesh/tnlIdenticalGridGeometry.h>
#include <core/tnlHost.h>
template< typename Mesh >
@@ -26,20 +27,51 @@ class tnlLinearDiffusion
{
};
+template< typename Real,
+ typename Device,
+ typename Index,
+ template< int, typename, typename, typename > class GridGeometry >
+class tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index, GridGeometry > >
+{
+ public:
+
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef tnlGrid< 2, Real, Device, Index, GridGeometry > MeshType;
+ typedef typename MeshType :: CoordinatesType CoordinatesType;
+ typedef typename MeshType :: VertexType VertexType;
+
+ tnlLinearDiffusion();
+
+ void bindMesh( const MeshType& mesh );
+
+ template< typename Vector >
+ void setFunction( Vector& f ); // TODO: add const
+
+ RealType getDiffusion( const Index& i ) const;
+ protected:
+
+ // TODO: change to ConstSharedVector
+ tnlSharedVector< RealType, DeviceType, IndexType > f;
+
+ const MeshType* mesh;
+};
+
template< typename Real, typename Device, typename Index >
-class tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index > >
+class tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry > >
{
public:
typedef Real RealType;
typedef Device DeviceType;
typedef Index IndexType;
- typedef typename tnlGrid< 2, Real, Device, Index > :: CoordinatesType CoordinatesType;
- typedef typename tnlGrid< 2, Real, Device, Index > :: VertexType VertexType;
+ typedef typename tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry > :: CoordinatesType CoordinatesType;
+ typedef typename tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry > :: VertexType VertexType;
tnlLinearDiffusion();
- void bindMesh( const tnlGrid< 2, RealType, DeviceType, IndexType >& mesh );
+ void bindMesh( const tnlGrid< 2, RealType, DeviceType, IndexType, tnlIdenticalGridGeometry >& mesh );
template< typename Vector >
void setFunction( Vector& f ); // TODO: add const
@@ -50,9 +82,10 @@ class tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index > >
// TODO: change to ConstSharedVector
tnlSharedVector< RealType, DeviceType, IndexType > f;
- const tnlGrid< 2, RealType, DeviceType, IndexType >* mesh;
+ const tnlGrid< 2, RealType, DeviceType, IndexType, tnlIdenticalGridGeometry >* mesh;
};
+
#include <implementation/schemes/diffusion/tnlLinearDiffusion_impl.h>
#endif
diff --git a/src/schemes/euler/fvm/tnlLaxFridrichs.h b/src/schemes/euler/fvm/tnlLaxFridrichs.h
index 76a496c368..8ccd7232fc 100644
--- a/src/schemes/euler/fvm/tnlLaxFridrichs.h
+++ b/src/schemes/euler/fvm/tnlLaxFridrichs.h
@@ -20,6 +20,7 @@
#include <core/tnlSharedVector.h>
#include <mesh/tnlGrid.h>
+#include <mesh/tnlIdenticalGridGeometry.h>
#include <schemes/gradient/tnlCentralFDMGradient.h>
template< typename MeshType,
@@ -28,15 +29,69 @@ class tnlLaxFridrichs
{
};
+template< typename Real,
+ typename Device,
+ typename Index,
+ typename PressureGradient,
+ template< int, typename, typename, typename > class GridGeometry >
+class tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index, GridGeometry >, PressureGradient >
+{
+ public:
+
+ typedef tnlGrid< 2, Real, Device, Index, GridGeometry > MeshType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef typename MeshType :: VertexType VertexType;
+ typedef typename MeshType :: CoordinatesType CoordinatesType;
+
+ tnlLaxFridrichs();
+
+ void getExplicitRhs( const IndexType centralVolume,
+ RealType& rho_t,
+ RealType& rho_u1_t,
+ RealType& rho_u2_t ) const;
+
+ void setRegularization( const RealType& epsilon );
+
+ void setViscosityCoefficient( const RealType& v );
+
+ void bindMesh( const MeshType& mesh );
+
+ template< typename Vector >
+ void setRho( Vector& rho ); // TODO: add const
+
+ template< typename Vector >
+ void setRhoU1( Vector& rho_u1 ); // TODO: add const
+
+ template< typename Vector >
+ void setRhoU2( Vector& rho_u2 ); // TODO: add const
+
+ template< typename Vector >
+ void setPressureGradient( Vector& grad_p ); // TODO: add const
+
+ protected:
+
+ RealType regularize( const RealType& r ) const;
+
+ RealType regularizeEps, viscosityCoefficient;
+
+ const MeshType* mesh;
+
+ const PressureGradient* pressureGradient;
+
+ tnlSharedVector< RealType, DeviceType, IndexType > rho, rho_u1, rho_u2;
+};
+
template< typename Real,
typename Device,
typename Index,
typename PressureGradient >
-class tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient >
+class tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry >, PressureGradient >
{
public:
- typedef tnlGrid< 2, Real, Device, Index > MeshType;
+ typedef tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry > MeshType;
typedef Real RealType;
typedef Device DeviceType;
typedef Index IndexType;
diff --git a/src/schemes/gradient/tnlCentralFDMGradient.h b/src/schemes/gradient/tnlCentralFDMGradient.h
index 4e07c69196..b88d455ecf 100644
--- a/src/schemes/gradient/tnlCentralFDMGradient.h
+++ b/src/schemes/gradient/tnlCentralFDMGradient.h
@@ -19,6 +19,7 @@
#define TNLCENTRALFDMGRADIENT_H_
#include <mesh/tnlGrid.h>
+#include <mesh/tnlIdenticalGridGeometry.h>
#include <core/tnlHost.h>
template< typename Mesh >
@@ -26,20 +27,24 @@ class tnlCentralFDMGradient
{
};
-template< typename Real, typename Device, typename Index >
-class tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index > >
+template< typename Real,
+ typename Device,
+ typename Index,
+ template< int, typename, typename, typename > class GridGeometry >
+class tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index, GridGeometry > >
{
public:
typedef Real RealType;
typedef Device DeviceType;
typedef Index IndexType;
- typedef typename tnlGrid< 2, Real, Device, Index > :: CoordinatesType CoordinatesType;
- typedef typename tnlGrid< 2, Real, Device, Index > :: VertexType VertexType;
+ typedef tnlGrid< 2, Real, Device, Index, GridGeometry > MeshType;
+ typedef typename MeshType :: CoordinatesType CoordinatesType;
+ typedef typename MeshType :: VertexType VertexType;
tnlCentralFDMGradient();
- void bindMesh( const tnlGrid< 2, RealType, DeviceType, IndexType >& mesh );
+ void bindMesh( const MeshType& mesh );
template< typename Vector >
void setFunction( Vector& f ); // TODO: add const
@@ -51,9 +56,41 @@ class tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index > >
// TODO: change to ConstSharedVector
tnlSharedVector< RealType, DeviceType, IndexType > f;
- const tnlGrid< 2, RealType, DeviceType, IndexType >* mesh;
+ const MeshType* mesh;
};
+template< typename Real,
+ typename Device,
+ typename Index >
+class tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry > >
+{
+ public:
+
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry > MeshType;
+ typedef typename MeshType :: CoordinatesType CoordinatesType;
+ typedef typename MeshType :: VertexType VertexType;
+
+ tnlCentralFDMGradient();
+
+ void bindMesh( const MeshType& mesh );
+
+ template< typename Vector >
+ void setFunction( Vector& f ); // TODO: add const
+
+ void getGradient( const Index& i,
+ VertexType& grad_f ) const;
+ protected:
+
+ // TODO: change to ConstSharedVector
+ tnlSharedVector< RealType, DeviceType, IndexType > f;
+
+ const MeshType* mesh;
+};
+
+
#include <implementation/schemes/gradient/tnlCentralFDMGradient_impl.h>
#endif
--
GitLab