From 4591fed4c7e6bf7e9fc37705a509b874aca17811 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jan=20Sch=C3=A4fer?= <schafjan@fjfi.cvut.cz> Date: Mon, 28 Mar 2016 16:40:40 +0200 Subject: [PATCH] dodelano podle pressureGetter --- .../inviscid-flow/1d/EulerPressureGetter.h | 3 +- .../inviscid-flow/1d/EulerPressureGetter.h~ | 215 ++--------- examples/inviscid-flow/1d/EulerVelGetter.h | 173 ++------- examples/inviscid-flow/1d/EulerVelGetter.h~ | 173 ++------- .../inviscid-flow/1d/EulerVelGetter_impl.h | 342 ------------------ examples/inviscid-flow/1d/eulerProblem.h | 2 +- examples/inviscid-flow/1d/eulerProblem.h~ | 23 +- examples/inviscid-flow/1d/eulerProblem_impl.h | 92 ++--- .../inviscid-flow/1d/eulerProblem_impl.h~ | 130 +++---- 9 files changed, 166 insertions(+), 987 deletions(-) delete mode 100644 examples/inviscid-flow/1d/EulerVelGetter_impl.h diff --git a/examples/inviscid-flow/1d/EulerPressureGetter.h b/examples/inviscid-flow/1d/EulerPressureGetter.h index 7e8becf84b..7743b96876 100644 --- a/examples/inviscid-flow/1d/EulerPressureGetter.h +++ b/examples/inviscid-flow/1d/EulerPressureGetter.h @@ -9,13 +9,12 @@ template< typename Mesh, typename Real = typename Mesh::RealType, typename Index = typename Mesh::IndexType > class EulerPressureGetter -: public tnlDomain< Mesh::geMeshDimensions(), MeshDomain > +: public tnlDomain< Mesh::getMeshDimensions(), MeshDomain > { public: typedef Mesh MeshType; typedef Real RealType; - typedef Device DeviceType; typedef Index IndexType; typedef tnlMeshFunction< MeshType > MeshFunctionType; enum { Dimensions = MeshType::getMeshDimensions() }; diff --git a/examples/inviscid-flow/1d/EulerPressureGetter.h~ b/examples/inviscid-flow/1d/EulerPressureGetter.h~ index daa62c4879..0200ac5666 100644 --- a/examples/inviscid-flow/1d/EulerPressureGetter.h~ +++ b/examples/inviscid-flow/1d/EulerPressureGetter.h~ @@ -3,213 +3,58 @@ #include <core/vectors/tnlVector.h> #include <mesh/tnlGrid.h> +#include <functions/tnlDomain.h> template< typename Mesh, typename Real = typename Mesh::RealType, typename Index = typename Mesh::IndexType > class EulerPressureGetter -{ -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class EulerPressureGetter< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index > +: public tnlDomain< Mesh::getMeshDimensions(), MeshDomain > { public: - typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType; - typedef typename MeshType::CoordinatesType CoordinatesType; + + typedef Mesh MeshType; typedef Real RealType; - typedef Device DeviceType; +// typedef Device DeviceType; typedef Index IndexType; typedef tnlMeshFunction< MeshType > MeshFunctionType; enum { Dimensions = MeshType::getMeshDimensions() }; static tnlString getType(); - Real gamma; - MeshFunctionType& velocity; - MeshFunctionType& rhoVel; - MeshFunctionType& energy; + + EulerPressureGetter( const MeshFunctionType& velocity, + const MeshFunctionType& rhoVel, + const MeshFunctionType& energy, + const RealType& gamma ) + : velocity( velocity ), rhoVel( rhoVel ), energy( energy ), gamma( gamma ) + {} - void setGamma(const Real& gamma) - { - this->gamma = gamma; - }; - - void setVelocity(const MeshFunctionType& velocity) - { - this->velocity = velocity; - }; - - void setRhoVel(const MeshFunctionType& rhoVel) - { - this->rhoVel = rhoVel; - }; - - void setEnergy(const MeshFunctionType& energy) - { - this->energy = energy; - }; - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - __cuda_callable__ template< typename MeshEntity > - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class EulerPressureGetter< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index > -{ - public: - typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType; - typedef typename MeshType::CoordinatesType CoordinatesType; - typedef Real RealType; - typedef Device DeviceType; - typedef Index IndexType; - typedef tnlMeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static tnlString getType(); - Real gamma; - MeshFunctionType& velocity; - MeshFunctionType& rhoVel; - MeshFunctionType& energy; - - void setGamma(const Real& gamma) - { - this->gamma = gamma; - }; - - void setVelocity(const MeshFunctionType& velocity) - { - this->velocity = velocity; - }; - - void setRhoVel(const MeshFunctionType& rhoVel) + Real operator()( const MeshEntity& entity, + const RealType& time = 0.0 ) const { - this->rhoVel = rhoVel; - }; - - void setEnergy(const MeshFunctionType& energy) - { - this->energy = energy; - }; - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - __cuda_callable__ + return this->operator[]( entity.getIndex() ); + } + template< typename MeshEntity > - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class EulerPressureGetter< tnlGrid< 3,MeshReal, Device, MeshIndex >, Real, Index > -{ - public: - typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType; - typedef typename MeshType::CoordinatesType CoordinatesType; - typedef Real RealType; - typedef Device DeviceType; - typedef Index IndexType; - typedef tnlMeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static tnlString getType(); - Real gamma; - MeshFunctionType& velocity; - MeshFunctionType& rhoVel; - MeshFunctionType& energy; - - void setGamma(const Real& gamma) - { - this->gamma = gamma; - }; - - void setVelocity(const MeshFunctionType& velocity) - { - this->velocity = velocity; - }; - - void setRhoVel(const MeshFunctionType& rhoVel) - { - this->rhoVel = rhoVel; - }; - - void setEnergy(const MeshFunctionType& energy) + Real operator[]( const IndexType& idx ) const { - this->energy = energy; - }; + return ( this->gamma - 1.0 ) * ( this->energy[ idx ] - 0.5 * this->rhoVel[ idx ] * this->velocity[ idx ]); + } - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - __cuda_callable__ - template< typename MeshEntity > - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; + + protected: + + Real gamma; + + const MeshFunctionType& velocity; + + const MeshFunctionType& rhoVel; + + const MeshFunctionType& energy; - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; }; - -#include "EulerPressureGetter_impl.h" - #endif /* EulerPressureGetter_H */ diff --git a/examples/inviscid-flow/1d/EulerVelGetter.h b/examples/inviscid-flow/1d/EulerVelGetter.h index 890a2e6696..7abb577da7 100644 --- a/examples/inviscid-flow/1d/EulerVelGetter.h +++ b/examples/inviscid-flow/1d/EulerVelGetter.h @@ -8,174 +8,45 @@ template< typename Mesh, typename Real = typename Mesh::RealType, typename Index = typename Mesh::IndexType > class EulerVelGetter -{ -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class EulerVelGetter< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index > +: public tnlDomain< Mesh::getMeshDimensions(), MeshDomain > { public: - typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType; - typedef typename MeshType::CoordinatesType CoordinatesType; + + typedef Mesh MeshType; typedef Real RealType; - typedef Device DeviceType; typedef Index IndexType; typedef tnlMeshFunction< MeshType > MeshFunctionType; enum { Dimensions = MeshType::getMeshDimensions() }; static tnlString getType(); - MeshFunctionType rhoVel; - MeshFunctionType rho; - - void setRhoVel(const MeshFunctionType& rhoVel) - { - this->rhoVel = rhoVel; - }; - - void setRho(const MeshFunctionType& rho) - { - this->rho = rho; - }; + + EulerVelGetter( const MeshFunctionType& rho, + const MeshFunctionType& rhoVel) + : rho( rho ), rhoVel( rhoVel ) + {} - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - __cuda_callable__ template< typename MeshEntity > - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class EulerVelGetter< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index > -{ - public: - typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType; - typedef typename MeshType::CoordinatesType CoordinatesType; - typedef Real RealType; - typedef Device DeviceType; - typedef Index IndexType; - typedef tnlMeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static tnlString getType(); - MeshFunctionType rhoVel; - MeshFunctionType rho; - - void setRhoVel(const MeshFunctionType& rhoVel) - { - this->rhoVel = rhoVel; - }; - - void setRho(const MeshFunctionType& rho) + Real operator()( const MeshEntity& entity, + const RealType& time = 0.0 ) const { - this->rho = rho; - }; - - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - __cuda_callable__ + return this->operator[]( entity.getIndex() ); + } + template< typename MeshEntity > - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class EulerVelGetter< tnlGrid< 3,MeshReal, Device, MeshIndex >, Real, Index > -{ - public: - typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType; - typedef typename MeshType::CoordinatesType CoordinatesType; - typedef Real RealType; - typedef Device DeviceType; - typedef Index IndexType; - typedef tnlMeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static tnlString getType(); - MeshFunctionType rhoVel; - MeshFunctionType rho; - - void setRhoVel(const MeshFunctionType& rhoVel) - { - this->rhoVel = rhoVel; - }; - - void setRho(const MeshFunctionType& rho) + Real operator[]( const IndexType& idx ) const { - this->rho = rho; - }; - - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; + return this->rho[ idx ] / this->rhoVel[ idx ]; + } - __cuda_callable__ - template< typename MeshEntity > - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; + + protected: + + const MeshFunctionType& rho; + + const MeshFunctionType& rhoVel; - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; }; - -#include "EulerVelGetter_impl.h" - #endif /* EulerVelGetter_H */ diff --git a/examples/inviscid-flow/1d/EulerVelGetter.h~ b/examples/inviscid-flow/1d/EulerVelGetter.h~ index 890a2e6696..7abb577da7 100644 --- a/examples/inviscid-flow/1d/EulerVelGetter.h~ +++ b/examples/inviscid-flow/1d/EulerVelGetter.h~ @@ -8,174 +8,45 @@ template< typename Mesh, typename Real = typename Mesh::RealType, typename Index = typename Mesh::IndexType > class EulerVelGetter -{ -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class EulerVelGetter< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index > +: public tnlDomain< Mesh::getMeshDimensions(), MeshDomain > { public: - typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType; - typedef typename MeshType::CoordinatesType CoordinatesType; + + typedef Mesh MeshType; typedef Real RealType; - typedef Device DeviceType; typedef Index IndexType; typedef tnlMeshFunction< MeshType > MeshFunctionType; enum { Dimensions = MeshType::getMeshDimensions() }; static tnlString getType(); - MeshFunctionType rhoVel; - MeshFunctionType rho; - - void setRhoVel(const MeshFunctionType& rhoVel) - { - this->rhoVel = rhoVel; - }; - - void setRho(const MeshFunctionType& rho) - { - this->rho = rho; - }; + + EulerVelGetter( const MeshFunctionType& rho, + const MeshFunctionType& rhoVel) + : rho( rho ), rhoVel( rhoVel ) + {} - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - __cuda_callable__ template< typename MeshEntity > - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class EulerVelGetter< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index > -{ - public: - typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType; - typedef typename MeshType::CoordinatesType CoordinatesType; - typedef Real RealType; - typedef Device DeviceType; - typedef Index IndexType; - typedef tnlMeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static tnlString getType(); - MeshFunctionType rhoVel; - MeshFunctionType rho; - - void setRhoVel(const MeshFunctionType& rhoVel) - { - this->rhoVel = rhoVel; - }; - - void setRho(const MeshFunctionType& rho) + Real operator()( const MeshEntity& entity, + const RealType& time = 0.0 ) const { - this->rho = rho; - }; - - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - __cuda_callable__ + return this->operator[]( entity.getIndex() ); + } + template< typename MeshEntity > - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class EulerVelGetter< tnlGrid< 3,MeshReal, Device, MeshIndex >, Real, Index > -{ - public: - typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType; - typedef typename MeshType::CoordinatesType CoordinatesType; - typedef Real RealType; - typedef Device DeviceType; - typedef Index IndexType; - typedef tnlMeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static tnlString getType(); - MeshFunctionType rhoVel; - MeshFunctionType rho; - - void setRhoVel(const MeshFunctionType& rhoVel) - { - this->rhoVel = rhoVel; - }; - - void setRho(const MeshFunctionType& rho) + Real operator[]( const IndexType& idx ) const { - this->rho = rho; - }; - - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; + return this->rho[ idx ] / this->rhoVel[ idx ]; + } - __cuda_callable__ - template< typename MeshEntity > - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; + + protected: + + const MeshFunctionType& rho; + + const MeshFunctionType& rhoVel; - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; }; - -#include "EulerVelGetter_impl.h" - #endif /* EulerVelGetter_H */ diff --git a/examples/inviscid-flow/1d/EulerVelGetter_impl.h b/examples/inviscid-flow/1d/EulerVelGetter_impl.h deleted file mode 100644 index 722324c403..0000000000 --- a/examples/inviscid-flow/1d/EulerVelGetter_impl.h +++ /dev/null @@ -1,342 +0,0 @@ -#ifndef EulerVelGetter_IMPL_H -#define EulerVelGetter_IMPL_H - -/**** - * 1D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -tnlString -EulerVelGetter< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return tnlString( "EulerVelGetter< " ) + - MeshType::getType() + ", " + - ::getType< Real >() + ", " + - ::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -EulerVelGetter< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 1, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 1, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - - //vel - const IndexType& center = entity.getIndex(); - return ( rhoVel[ center ] / rho[ center ]); -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -EulerVelGetter< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -EulerVelGetter< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); - matrixRow.setElement( 0, west, - lambdaX ); - matrixRow.setElement( 1, center, 2.0 * lambdaX ); - matrixRow.setElement( 2, east, - lambdaX ); -} - -/**** - * 2D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -tnlString -EulerVelGetter< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return tnlString( "EulerVelGetter< " ) + - MeshType::getType() + ", " + - ::getType< Real >() + ", " + - ::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -EulerVelGetter< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 2, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 2, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); - - const RealType& hxSquareInverse = entity.getMesh().template getSpaceStepsProducts< -2, 0 >(); - const RealType& hySquareInverse = entity.getMesh().template getSpaceStepsProducts< 0, -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse + - ( u[ south ] - 2.0 * u[ center ] + u[ north ] ) * hySquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -EulerVelGetter< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -EulerVelGetter< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2, 0 >(); - const RealType& lambdaY = tau * entity.getMesh().template getSpaceStepsProducts< 0, -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); - matrixRow.setElement( 0, south, -lambdaY ); - matrixRow.setElement( 1, west, -lambdaX ); - matrixRow.setElement( 2, center, 2.0 * ( lambdaX + lambdaY ) ); - matrixRow.setElement( 3, east, -lambdaX ); - matrixRow.setElement( 4, north, -lambdaY ); -} - -/**** - * 3D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -tnlString -EulerVelGetter< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return tnlString( "EulerVelGetter< " ) + - MeshType::getType() + ", " + - ::getType< Real >() + ", " + - ::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -EulerVelGetter< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 3, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 3, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - - const RealType& hxSquareInverse = entity.getMesh().template getSpaceStepsProducts< -2, 0, 0 >(); - const RealType& hySquareInverse = entity.getMesh().template getSpaceStepsProducts< 0, -2, 0 >(); - const RealType& hzSquareInverse = entity.getMesh().template getSpaceStepsProducts< 0, 0, -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse + - ( u[ south ] - 2.0 * u[ center ] + u[ north ] ) * hySquareInverse + - ( u[ up ] - 2.0 * u[ center ] + u[ down ] ) * hzSquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -EulerVelGetter< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -EulerVelGetter< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2, 0, 0 >(); - const RealType& lambdaY = tau * entity.getMesh().template getSpaceStepsProducts< 0, -2, 0 >(); - const RealType& lambdaZ = tau * entity.getMesh().template getSpaceStepsProducts< 0, 0, -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - matrixRow.setElement( 0, down, -lambdaZ ); - matrixRow.setElement( 1, south, -lambdaY ); - matrixRow.setElement( 2, west, -lambdaX ); - matrixRow.setElement( 3, center, 2.0 * ( lambdaX + lambdaY + lambdaZ ) ); - matrixRow.setElement( 4, east, -lambdaX ); - matrixRow.setElement( 5, north, -lambdaY ); - matrixRow.setElement( 6, up, -lambdaZ ); -} - -#endif /* EulerVelGetterIMPL_H */ - diff --git a/examples/inviscid-flow/1d/eulerProblem.h b/examples/inviscid-flow/1d/eulerProblem.h index 849b12d1e3..6bce37a02c 100644 --- a/examples/inviscid-flow/1d/eulerProblem.h +++ b/examples/inviscid-flow/1d/eulerProblem.h @@ -95,7 +95,7 @@ class eulerProblem: MeshFunctionType uRho, uRhoVelocity, uEnergy; MeshFunctionType fuRho, fuRhoVelocity, fuEnergy; - MeshFunctionType rho, rhoVel, energy, pressure, velocity; + MeshFunctionType pressure, velocity, rho, rhoVel, energy; RealType gamma; diff --git a/examples/inviscid-flow/1d/eulerProblem.h~ b/examples/inviscid-flow/1d/eulerProblem.h~ index fc0803a615..cf094a3ba6 100644 --- a/examples/inviscid-flow/1d/eulerProblem.h~ +++ b/examples/inviscid-flow/1d/eulerProblem.h~ @@ -33,23 +33,9 @@ class eulerProblem: typedef typename DifferentialOperator::Velocity Velocity; typedef typename DifferentialOperator::Pressure Pressure; - //definition - tnlVector< RealType, DeviceType, IndexType > _uRho; - tnlVector< RealType, DeviceType, IndexType > _uRhoVelocity; - tnlVector< RealType, DeviceType, IndexType > _uEnergy; - - tnlVector< RealType, DeviceType, IndexType > _fuRho; - tnlVector< RealType, DeviceType, IndexType > _fuRhoVelocity; - tnlVector< RealType, DeviceType, IndexType > _fuEnergy; static tnlString getTypeStatic(); - tnlVector< RealType, DeviceType, IndexType > rho; - tnlVector< RealType, DeviceType, IndexType > rhoVel; - tnlVector< RealType, DeviceType, IndexType > energy; - tnlVector< RealType, DeviceType, IndexType > pressure; - tnlVector< RealType, DeviceType, IndexType > velocity; - double gamma; tnlString getPrologHeader() const; @@ -105,6 +91,15 @@ class eulerProblem: DifferentialOperator differentialOperator; BoundaryCondition boundaryCondition; RightHandSide rightHandSide; + + MeshFunctionType uRho, uRhoVelocity, uEnergy; + MeshFunctionType fuRho, fuRhoVelocity, fuEnergy; + + MeshFunctionType pressure, velocity rho, rhoVel, energy; + + RealType gamma; + + }; #include "eulerProblem_impl.h" diff --git a/examples/inviscid-flow/1d/eulerProblem_impl.h b/examples/inviscid-flow/1d/eulerProblem_impl.h index e97cc2f46e..ea0b55fdb7 100644 --- a/examples/inviscid-flow/1d/eulerProblem_impl.h +++ b/examples/inviscid-flow/1d/eulerProblem_impl.h @@ -100,7 +100,7 @@ setInitialCondition( const tnlParameterContainer& parameters, MeshDependentDataType& meshDependentData ) { typedef typename MeshType::Cell Cell; - double gamma = parameters.getParameter< double >( "gamma" ); + this->gamma = parameters.getParameter< double >( "gamma" ); double rhoL = parameters.getParameter< double >( "left-density" ); double velL = parameters.getParameter< double >( "left-velocity" ); double preL = parameters.getParameter< double >( "left-pressure" ); @@ -112,33 +112,31 @@ setInitialCondition( const tnlParameterContainer& parameters, double x0 = parameters.getParameter< double >( "riemann-border" ); cout << gamma << " " << rhoL << " " << velL << " " << preL << " " << eL << " " << rhoR << " " << velR << " " << preR << " " << eR << " " << x0 << " " << gamma << endl; int count = mesh.template getEntitiesCount< Cell >()/3; - this->rho.bind(dofs,0,count); - this->rhoVel.bind(dofs,count,count); - this->energy.bind(dofs,2 * count,count); - this->data.setSize(2*count); - this->pressure.bind(this->data,0,count); - this->velocity.bind(this->data,count,count); - for(long int i = 0; i < count; i++) + rho.bind(mesh,dofs,0); + rhoVel.bind(mesh,dofs,count); + energy.bind(mesh,dofs,2 * count); +// pressure(mesh); +// velocity(mesh); +/* for(long int i = 0; i < count; i++) if (i < x0 * count ) { - this->rho[i] = rhoL; - this->rhoVel[i] = rhoL * velL; - this->energy[i] = eL; + this->Rho[i] = rhoL; + this->RhoVel[i] = rhoL * velL; + this->Energy[i] = eL; this->velocity[i] = velL; this->pressure[i] = preL; } else { - this->rho[i] = rhoR; - this->rhoVel[i] = rhoR * velR; - this->energy[i] = eR; + this->Rho[i] = rhoR; + this->RhoVel[i] = rhoR * velR; + this->Energy[i] = eR; this->velocity[i] = velR; this->pressure[i] = preR; }; - this->gamma = gamma; cout << "dofs = " << dofs << endl; getchar(); - +*/ /* const tnlString& initialConditionFile = parameters.getParameter< tnlString >( "initial-condition" ); @@ -187,7 +185,7 @@ makeSnapshot( const RealType& time, const MeshType& mesh, DofVectorType& dofs, MeshDependentDataType& meshDependentData ) -{ +{/* cout << endl << "Writing output at time " << time << " step " << step << "." << endl; this->bindDofs( mesh, dofs ); tnlString fileName; @@ -228,7 +226,7 @@ makeSnapshot( const RealType& time, return false; FileNameBaseNumberEnding( "energy-", step, 5, ".tnl", fileName ); if( ! energy.save( fileName ) ) - return false; + return false;*/ return true; } @@ -246,22 +244,17 @@ getExplicitRHS( const RealType& time, MeshDependentDataType& meshDependentData ) { typedef typename MeshType::Cell Cell; - int count = mesh.template getEntitiesCount< Cell >(); ///3; + int count = mesh.template getEntitiesCount< Cell >(); //bind _u - this->uRho.bind(_u,0,count); - this->uRhoVelocity.bind(_u,count,count); - this->uEnergy.bind(_u,2 * count,count); + this->uRho.bind(mesh, _u, 0); + this->uRhoVelocity.bind(mesh, _u ,count); + this->uEnergy.bind(mesh, _u, 2 * count); //bind _fu - this->fuRho.bind(_u,0,count); - this->fuRhoVelocity.bind(_u,count,count); - this->fuEnergy.bind(_u,2 * count,count); - //bind MeshFunctions - this->velocity.setMesh( mesh ); - // TODO: osetrit rychlost podobne jako tlak - this->pressure.setMesh( mesh ); - EulerPressureGetter< Mesh, RealType, IndexType > pressureGetter( velocity, rhoVel, energy, gamma ); - this->pressure = pressureGetter; + this->fuRho.bind(mesh, _u, 0); + this->fuRhoVelocity.bind(mesh, _u, count); + fuEnergy.bind(mesh, _u, 2 * count); + //generating Differential operator object Continuity lF1DContinuity; Momentum lF1DMomentum; @@ -365,37 +358,14 @@ postIterate( const RealType& time, DofVectorType& dofs, MeshDependentDataType& meshDependentData ) { - typedef typename MeshType::Cell Cell; - int count = mesh.template getEntitiesCount< Cell >()/3; - - //bind _u - this->_uRho.bind(dofs, 0, count); - this->_uRhoVelocity.bind(dofs, count, count); - this->_uEnergy.bind(dofs, 2 * count, count); - - //bind meshfunction - MeshFunctionType velocity( mesh, this->velocity ); - MeshFunctionType pressure( mesh, this->pressure ); - MeshFunctionType uRho( mesh, _uRho ); - MeshFunctionType uRhoVel( mesh, _uRhoVelocity ); - MeshFunctionType uEnergy( mesh, _uEnergy ); - //Generating differential operators - Velocity euler1DVelocity; - Pressure euler1DPressure; //velocity - euler1DVelocity.setRhoVel(uRhoVel); - euler1DVelocity.setRho(uRho); -// tnlOperatorFunction< Velocity, MeshFunction, void, true > OFVel; -// velocity = OFVel; - - //pressure - euler1DPressure.setRhoVel(uRhoVel); - euler1DPressure.setVelocity(velocity); - euler1DPressure.setGamma(gamma); - euler1DPressure.setEnergy(uEnergy); -// tnlOperatorFunction< Pressure, MeshFunction, void, true > OFPressure; -// pressure = OFPressure; - + this->velocity.setMesh( mesh ); + Velocity velocityGetter( uRho, uRhoVelocity ); + this->velocity = velocityGetter; +/* //pressure + this->pressure.setMesh( mesh ); + Pressure pressureGetter( velocity, RhoVelocity, Energy, gamma ); + this->pressure = pressureGetter; */ } #endif /* eulerPROBLEM_IMPL_H_ */ diff --git a/examples/inviscid-flow/1d/eulerProblem_impl.h~ b/examples/inviscid-flow/1d/eulerProblem_impl.h~ index 33b5e00a83..4074f7a39a 100644 --- a/examples/inviscid-flow/1d/eulerProblem_impl.h~ +++ b/examples/inviscid-flow/1d/eulerProblem_impl.h~ @@ -100,7 +100,7 @@ setInitialCondition( const tnlParameterContainer& parameters, MeshDependentDataType& meshDependentData ) { typedef typename MeshType::Cell Cell; - double gamma = parameters.getParameter< double >( "gamma" ); + this->gamma = parameters.getParameter< double >( "gamma" ); double rhoL = parameters.getParameter< double >( "left-density" ); double velL = parameters.getParameter< double >( "left-velocity" ); double preL = parameters.getParameter< double >( "left-pressure" ); @@ -112,33 +112,31 @@ setInitialCondition( const tnlParameterContainer& parameters, double x0 = parameters.getParameter< double >( "riemann-border" ); cout << gamma << " " << rhoL << " " << velL << " " << preL << " " << eL << " " << rhoR << " " << velR << " " << preR << " " << eR << " " << x0 << " " << gamma << endl; int count = mesh.template getEntitiesCount< Cell >()/3; - this->rho.bind(dofs,0,count); - this->rhoVel.bind(dofs,count,count); - this->energy.bind(dofs,2 * count,count); - this->data.setSize(2*count); - this->pressure.bind(this->data,0,count); - this->velocity.bind(this->data,count,count); - for(long int i = 0; i < count; i++) + rho.bind(mesh,dofs,0); + rhoVel.bind(mesh,dofs,count); + energy.bind(mesh,dofs,2 * count); + pressure(mesh); +// velocity(mesh); +/* for(long int i = 0; i < count; i++) if (i < x0 * count ) { - this->rho[i] = rhoL; - this->rhoVel[i] = rhoL * velL; - this->energy[i] = eL; + this->Rho[i] = rhoL; + this->RhoVel[i] = rhoL * velL; + this->Energy[i] = eL; this->velocity[i] = velL; this->pressure[i] = preL; } else { - this->rho[i] = rhoR; - this->rhoVel[i] = rhoR * velR; - this->energy[i] = eR; + this->Rho[i] = rhoR; + this->RhoVel[i] = rhoR * velR; + this->Energy[i] = eR; this->velocity[i] = velR; this->pressure[i] = preR; }; - this->gamma = gamma; cout << "dofs = " << dofs << endl; getchar(); - +*/ /* const tnlString& initialConditionFile = parameters.getParameter< tnlString >( "initial-condition" ); @@ -187,7 +185,7 @@ makeSnapshot( const RealType& time, const MeshType& mesh, DofVectorType& dofs, MeshDependentDataType& meshDependentData ) -{ +{/* cout << endl << "Writing output at time " << time << " step " << step << "." << endl; this->bindDofs( mesh, dofs ); tnlString fileName; @@ -228,7 +226,7 @@ makeSnapshot( const RealType& time, return false; FileNameBaseNumberEnding( "energy-", step, 5, ".tnl", fileName ); if( ! energy.save( fileName ) ) - return false; + return false;*/ return true; } @@ -246,51 +244,46 @@ getExplicitRHS( const RealType& time, MeshDependentDataType& meshDependentData ) { typedef typename MeshType::Cell Cell; - int count = mesh.template getEntitiesCount< Cell >()/3; - //bind _u - this->_uRho.bind(_u,0,count); - this->_uRhoVelocity.bind(_u,count,count); - this->_uEnergy.bind(_u,2 * count,count); + int count = mesh.template getEntitiesCount< Cell >(); ///3; +/* //bind _u + this->uRho.bind(_u,0,count); + this->uRhoVelocity.bind(_u,count,count); + this->uEnergy.bind(_u,2 * count,count); //bind _fu - this->_fuRho.bind(_u,0,count); - this->_fuRhoVelocity.bind(_u,count,count); - this->_fuEnergy.bind(_u,2 * count,count); - //bind MeshFunctions - MeshFunctionType velocity( mesh, this->velocity ); - MeshFunctionType pressure( mesh, this->pressure ); - MeshFunctionType uRho( mesh, _uRho ); - MeshFunctionType fuRho( mesh, _fuRho ); - MeshFunctionType uRhoVelocity( mesh, _uRhoVelocity ); - MeshFunctionType fuRhoVelocity( mesh, _fuRhoVelocity ); - MeshFunctionType uEnergy( mesh, _uEnergy ); - MeshFunctionType fuEnergy( mesh, _fuEnergy ); - //generating Differential operator object - Continuity lFContinuity; - Momentum lFMomentum; - Energy lFEnergy; + this->fuRho.bind(_u,0,count); + this->fuRhoVelocity.bind(_u,count,count); + fuEnergy.bind(_u,2 * count,count); +*/ + //generating Differential operator object + Continuity lF1DContinuity; + Momentum lF1DMomentum; + Energy lF1DEnergy; + + + //rho this->bindDofs( mesh, _u ); - lFContinuity.setTau(tau); - lFContinuity.setVelocity(velocity); + lF1DContinuity.setTau(tau); + lF1DContinuity.setVelocity(velocity); tnlExplicitUpdater< Mesh, MeshFunctionType, Continuity, BoundaryCondition, RightHandSide > explicitUpdaterContinuity; explicitUpdaterContinuity.template update< typename Mesh::Cell >( time, mesh, - lFContinuity, + lF1DContinuity, this->boundaryCondition, this->rightHandSide, uRho, fuRho ); //rhoVelocity - lFMomentum.setTau(tau); - lFMomentum.setVelocity(velocity); - lFMomentum.setPressure(pressure); + lF1DMomentum.setTau(tau); + lF1DMomentum.setVelocity(velocity); + lF1DMomentum.setPressure(pressure); tnlExplicitUpdater< Mesh, MeshFunctionType, Momentum, BoundaryCondition, RightHandSide > explicitUpdaterMomentum; explicitUpdaterMomentum.template update< typename Mesh::Cell >( time, mesh, - lFMomentum, + lF1DMomentum, this->boundaryCondition, this->rightHandSide, uRhoVelocity, @@ -298,13 +291,13 @@ getExplicitRHS( const RealType& time, //energy - lFEnergy.setTau(tau); - lFEnergy.setPressure(pressure); - lFEnergy.setVelocity(velocity); + lF1DEnergy.setTau(tau); + lF1DEnergy.setPressure(pressure); + lF1DEnergy.setVelocity(velocity); tnlExplicitUpdater< Mesh, MeshFunctionType, Energy, BoundaryCondition, RightHandSide > explicitUpdaterEnergy; explicitUpdaterEnergy.template update< typename Mesh::Cell >( time, mesh, - lFEnergy, + lF1DEnergy, this->boundaryCondition, this->rightHandSide, uEnergy, @@ -365,37 +358,14 @@ postIterate( const RealType& time, DofVectorType& dofs, MeshDependentDataType& meshDependentData ) { - typedef typename MeshType::Cell Cell; - int count = mesh.template getEntitiesCount< Cell >()/3; - - //bind _u - this->_uRho.bind(dofs, 0, count); - this->_uRhoVelocity.bind(dofs, count, count); - this->_uEnergy.bind(dofs, 2 * count, count); - - //bind meshfunction - MeshFunctionType velocity( mesh, this->velocity ); - MeshFunctionType pressure( mesh, this->pressure ); - MeshFunctionType uRho( mesh, _uRho ); - MeshFunctionType uRhoVel( mesh, _uRhoVelocity ); - MeshFunctionType uEnergy( mesh, _uEnergy ); - //Generating differential operators - Velocity euler1DVelocity; - Pressure euler1DPressure; //velocity - euler1DVelocity.setRhoVel(uRhoVel); - euler1DVelocity.setRho(uRho); -// tnlOperatorFunction< Velocity, MeshFunction, void, true > OFVel; -// velocity = OFVel; - - //pressure - euler1DPressure.setRhoVel(uRhoVel); - euler1DPressure.setVelocity(velocity); - euler1DPressure.setGamma(gamma); - euler1DPressure.setEnergy(uEnergy); -// tnlOperatorFunction< Pressure, MeshFunction, void, true > OFPressure; -// pressure = OFPressure; - + this->velocity.setMesh( mesh ); + Velocity velocityGetter( uRho, uRhoVelocity ); + this->velocity = velocityGetter; +/* //pressure + this->pressure.setMesh( mesh ); + Pressure pressureGetter( velocity, RhoVelocity, Energy, gamma ); + this->pressure = pressureGetter; */ } #endif /* eulerPROBLEM_IMPL_H_ */ -- GitLab