Implementing the gradient operator for the deformed grids. (30c8e380) · Commits · TNL / tnl-dev

src/core/mfuncs.h

+2 −2

Original line number	Diff line number	Diff line
		@@ -68,13 +68,13 @@ inline double tnlAbs( const double& d )
		return fabs( d );
		};

		template< typename T >
		/*template< typename T >
		void swap( T& a, T& b)
		{
		T aux;
		aux = a;
		a = b;
		b = aux;
		}
		}*/

		#endif

src/implementation/mesh/tnlGrid2D_impl.h

+5 −5

Original line number	Diff line number	Diff line
		@@ -247,17 +247,17 @@ Real tnlGrid< 2, Real, Device, Index, Geometry > :: getElementMeasure( const Coo
		return geometry. getElementMeasure( coordinates );
		}

		/*template< typename Real,
		template< typename Real,
		typename Device,
		typename Index,
		template< int, typename, typename, typename > class Geometry >
		template< Index dy, Index dx >
		Real tnlGrid< 2, Real, Device, Index, Geometry > :: getElementsDistance( const Index j,
		const Index i ) const
		Real tnlGrid< 2, Real, Device, Index, Geometry > :: getElementsDistance( const CoordinatesType& c1,
		const CoordinatesType& c2 ) const
		{
		return geometry. getElementsDistance< dy, dx >( j, i );
		return geometry. getElementsDistance( c1, c2 );
		}

		/*
		template< typename Real,
		typename Device,
		typename Index,

src/implementation/mesh/tnlIdenticalGridGeometry_impl.h

+11 −11

Original line number	Diff line number	Diff line
		@@ -163,22 +163,22 @@ Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementMeasure( co
		return elementMeasure;
		}

		/*template< typename Real,
		template< typename Real,
		typename Device,
		typename Index >
		template< Index dx, Index dy >
		Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementsDistance( const Index i,
		const Index j ) const
		Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementsDistance( const CoordinatesType& c1,
		const CoordinatesType& c2 ) const
		{
		if( dy == 0 && dx == 1 )
		return parametricStep. x();
		if( dy == 1 && dx == 0 )
		return parametricStep. y();
		const Real x = dx * parametricStep. x();
		const Real y = dy * parametricStep. y();
		return sqrt( dx * dx + dy * dy );
		CoordinatesType c( c2 );
		c -= c1;
		if( c. y() == 0 )
		return parametricStep. x() * fabs( c. x() );
		if( c. x() == 0 )
		return parametricStep. y() * fabs( c. y() );
		return sqrt( c. x() * c. x() + c. y() * c. y() );
		}

		/*
		template< typename Real,
		typename Device,
		typename Index >

src/implementation/schemes/euler/fvm/tnlLaxFridrichs_impl.h

+54 −6

Original line number	Diff line number	Diff line
		@@ -139,8 +139,8 @@ void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient > ::
		/****
		* Compute the pressure gradient
		*/
		RealType p_x( 0.0 ), p_y( 0.0 );
		pressureGradient -> getGradient( c, p_x, p_y );
		VertexType grad_p;
		pressureGradient -> getGradient( c, grad_p );

		/****
		* ( rho * u1 )_t + ( rho * u1 * u1 )_x + ( rho * u1 * u2 )_y - p_x = 0
		@@ -148,14 +148,14 @@ void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient > ::
		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 )
		- p_x;
		- grad_p. x();
		/****
		* ( rho * u1 )_t + ( rho * u1 * u1 )_x + ( rho * u1 * u2 )_y - p_y = 0
		* ( 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 )
		- p_y;
		- grad_p. y();


		/****
		@@ -210,16 +210,64 @@ void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient > ::
		const RealType rho_g_n = 0.5 * ( rho[ c ] * u2_c + rho[ n ] * u2_n );
		const RealType rho_g_s = 0.5 * ( rho[ c ] * u2_c + rho[ s ] * u2_s );

		const RealType rho_u1_f_e = 0.5 * ( rho_u1[ c ] * u1_c + rho_u1[ e ] * u1_e );
		const RealType rho_u1_f_w = 0.5 * ( rho_u1[ c ] * u1_c + rho_u1[ w ] * u1_w );
		const RealType rho_u1_f_n = 0.5 * ( rho_u1[ c ] * u1_c + rho_u1[ n ] * u1_n );
		const RealType rho_u1_f_s = 0.5 * ( rho_u1[ c ] * u1_c + rho_u1[ s ] * u1_s );
		const RealType rho_u1_g_e = 0.5 * ( rho_u1[ c ] * u2_c + rho_u1[ e ] * u2_e );
		const RealType rho_u1_g_w = 0.5 * ( rho_u1[ c ] * u2_c + rho_u1[ w ] * u2_w );
		const RealType rho_u1_g_n = 0.5 * ( rho_u1[ c ] * u2_c + rho_u1[ n ] * u2_n );
		const RealType rho_u1_g_s = 0.5 * ( rho_u1[ c ] * u2_c + rho_u1[ s ] * u2_s );

		const RealType rho_u2_f_e = 0.5 * ( rho_u2[ c ] * u1_c + rho_u2[ e ] * u1_e );
		const RealType rho_u2_f_w = 0.5 * ( rho_u2[ c ] * u1_c + rho_u2[ w ] * u1_w );
		const RealType rho_u2_f_n = 0.5 * ( rho_u2[ c ] * u1_c + rho_u2[ n ] * u1_n );
		const RealType rho_u2_f_s = 0.5 * ( rho_u2[ c ] * u1_c + rho_u2[ s ] * u1_s );
		const RealType rho_u2_g_e = 0.5 * ( rho_u2[ c ] * u2_c + rho_u2[ e ] * u2_e );
		const RealType rho_u2_g_w = 0.5 * ( rho_u2[ c ] * u2_c + rho_u2[ w ] * u2_w );
		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 );

		/****
		* rho_t + ( rho u_1 )_x + ( rho u_2 )_y = 0
		*/
		rho_t = - 1.0 / mu_D_c * ( rho_f_e * e_normal. x() + rho_g_e * e_normal. y() +
		rho_f_n * n_normal. x() + rho_g_n * n_normal. y() +
		rho_f_w * w_normal. x() + rho_g_w * w_normal. y() +
		rho_f_s * s_normal. x() + rho_g_s * s_normal. y() )
		+ 1.0 / ( 8.0 * mu_D_c ) *
		+ this -> viscosityCoefficient * 1.0 / 8.0 * mu_D_c *
		( ( mu_D_c + mu_D_e ) * ( rho[ e ] - rho[ c ] ) +
		( mu_D_c + mu_D_n ) * ( rho[ n ] - rho[ c ] ) +
		( mu_D_c + mu_D_w ) * ( rho[ w ] - rho[ c ] ) +
		( mu_D_c + mu_D_s ) * ( rho[ s ] - rho[ c ] ) );

		/****
		* ( rho * u1 )_t + ( rho * u1 * u1 )_x + ( rho * u1 * u2 )_y - p_x = 0
		*/
		rho_u1_t = - 1.0 / mu_D_c * ( rho_u1_f_e * e_normal. x() + rho_u1_g_e * e_normal. y() +
		rho_u1_f_n * n_normal. x() + rho_u1_g_n * n_normal. y() +
		rho_u1_f_w * w_normal. x() + rho_u1_g_w * w_normal. y() +
		rho_u1_f_s * s_normal. x() + rho_u1_g_s * s_normal. y() )
		+ this -> viscosityCoefficient * 1.0 / 8.0 * mu_D_c *
		( ( mu_D_c + mu_D_e ) * ( rho_u1[ e ] - rho_u1[ c ] ) +
		( mu_D_c + mu_D_n ) * ( rho_u1[ n ] - rho_u1[ c ] ) +
		( mu_D_c + mu_D_w ) * ( rho_u1[ w ] - rho_u1[ c ] ) +
		( mu_D_c + mu_D_s ) * ( rho_u1[ s ] - rho_u1[ c ] ) )
		- grad_p. x();

		/****
		* ( rho * u2 )_t + ( rho * u2 * u1 )_x + ( rho * u2 * u2 )_y - p_y = 0
		*/
		rho_u2_t = - 1.0 / mu_D_c * ( rho_u2_f_e * e_normal. x() + rho_u2_g_e * e_normal. y() +
		rho_u2_f_n * n_normal. x() + rho_u2_g_n * n_normal. y() +
		rho_u2_f_w * w_normal. x() + rho_u2_g_w * w_normal. y() +
		rho_u2_f_s * s_normal. x() + rho_u2_g_s * s_normal. y() )
		+ this -> viscosityCoefficient * 1.0 / 8.0 * mu_D_c *
		( ( mu_D_c + mu_D_e ) * ( rho_u2[ e ] - rho_u2[ c ] ) +
		( mu_D_c + mu_D_n ) * ( rho_u2[ n ] - rho_u2[ c ] ) +
		( mu_D_c + mu_D_w ) * ( rho_u2[ w ] - rho_u2[ c ] ) +
		( mu_D_c + mu_D_s ) * ( rho_u2[ s ] - rho_u2[ c ] ) )
		- grad_p. y();
		}

		template< typename Real,

src/implementation/schemes/gradient/tnlCentralFDMGradient_impl.h

+17 −4

Original line number	Diff line number	Diff line
		@@ -40,8 +40,7 @@ void tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index > > :: setFunction(

		template< typename Real, typename Device, typename Index >
		void tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index > > :: getGradient( const Index& i,
		RealType& f_x,
		RealType& f_y ) const
		VertexType& grad_f ) const
		{
		tnlAssert( this -> mesh, cerr << "No mesh was set in tnlCentralFDMGradient. Use the bindMesh method." );

		@@ -49,9 +48,23 @@ void tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index > > :: getGradient(
		const Index w = mesh -> getElementNeighbour( i, -1, 0 );
		const Index n = mesh -> getElementNeighbour( i, 0, 1 );
		const Index s = mesh -> getElementNeighbour( i, 0, -1 );
		CoordinatesType cCoordinates;
		mesh -> getElementCoordinates( i, cCoordinates );
		CoordinatesType eCoordinates( cCoordinates ),
		wCoordinates( cCoordinates ),
		nCoordinates( cCoordinates ),
		sCoordinates( cCoordinates );
		eCoordinates. x() ++;
		wCoordinates. x() --;
		nCoordinates. y() ++;
		sCoordinates. y() --;

		f_x = ( f[ e ] - f[ w ] ) / ( 2.0 * mesh -> getParametricStep(). x() );
		f_y = ( f[ n ] - f[ s ] ) / ( 2.0 * mesh -> getParametricStep(). 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 ) );
		}

		#endif