diff --git a/examples/advection/advectionProblem_impl.h b/examples/advection/advectionProblem_impl.h
index bccdeccdfcb44928f4dc961e24981d561638af76..1ae1259782fc5f50b08a705cdf849e9183fd6498 100644
--- a/examples/advection/advectionProblem_impl.h
+++ b/examples/advection/advectionProblem_impl.h
@@ -173,7 +173,7 @@ setInitialCondition( const tnlParameterContainer& parameters,
for (IndexType i = 0; i < inverseSquareCount-1; i++)
for (IndexType j = 0; j < inverseSquareCount-1; j++)
{
- dofs[i * inverseSquareCount + j] = exp(-pow(10/(size*inverseSquareCount)*(size*i-0.2*size*inverseSquareCount),2)-pow(10/(size*inverseSquareCount)*(size*j-0.2*size*inverseSquareCount)*size,2));
+ dofs[i * inverseSquareCount + j] = exp(-pow(10/(size*inverseSquareCount)*(size*i-0.2*size*inverseSquareCount),2)-pow(10/(size*inverseSquareCount)*(size*j-0.2*size*inverseSquareCount),2));
};
};
}
@@ -284,13 +284,14 @@ makeSnapshot( const RealType& time,
cout << endl << "Writing output at time " << time << " step " << step << "." << endl;
this->bindDofs( mesh, dofs );
tnlString fileName;
+ MeshFunctionType dofsh;
+ dofsh.bind(mesh,dofs);
FileNameBaseNumberEnding( "u-", step, 5, ".tnl", fileName );
- if( ! dofs.save( fileName ) )
+ if( ! dofsh.save( fileName ) )
return false;
FileNameBaseNumberEnding( "a-", step, 5, ".tnl", fileName );
if( ! this->analyt.save( fileName ) )
return false;
-cin.ignore();
return true;
}
@@ -317,7 +318,7 @@ getExplicitRHS( const RealType& time,
{
if (dimension == 1)
{
- this->analyt[0];
+ this->analyt[0] = 0;
for (IndexType i = 1; i < count-2; i++)
{
if ((i - step * tau * (count/this->schemeSize) * this -> speedX>0.4*count) && (i - step * tau * (count/this->schemeSize) * this -> speedX<0.5*count)) analyt[i]=1; else analyt[i]=0;
@@ -331,7 +332,7 @@ getExplicitRHS( const RealType& time,
{
if ((i - step * tau * (inverseSquareCount/this->schemeSize) * this -> speedX>0.4*inverseSquareCount) && (i - step * tau * (inverseSquareCount/this->schemeSize) * this -> speedX<0.5*inverseSquareCount)
&& (j - step * tau * (inverseSquareCount/this->schemeSize) * this -> speedY>0.4*inverseSquareCount) && (j - step * tau * (inverseSquareCount/this->schemeSize) * this -> speedY<0.5*inverseSquareCount))
- analyt[i]=1; else analyt[i]=0;
+ analyt[i * inverseSquareCount + j]=1; else analyt[i * inverseSquareCount + j]=0;
};
};
}
@@ -379,8 +380,8 @@ getExplicitRHS( const RealType& time,
else if (dimension == 2)
{
count = sqrt(count);
- for (IndexType i = 1; i < inverseSquareCount-1; i++)
- for (IndexType j = 1; j < inverseSquareCount-1; j++)
+ for (IndexType i = 0; i < inverseSquareCount-1; i++)
+ for (IndexType j = 0; j < inverseSquareCount-1; j++)
{
this->analyt[i * inverseSquareCount + j] = exp(-pow(10/(size*inverseSquareCount)*(size*i-0.2*size*inverseSquareCount)-step * 10 * tau * this->speedX,2)-pow(10/(size*inverseSquareCount)*(size*j-0.2*size*inverseSquareCount)-step * 10 * tau * this->speedY,2));
};
@@ -401,8 +402,8 @@ getExplicitRHS( const RealType& time,
else if (dimension == 2)
{
count = sqrt(count);
- for (IndexType i = 1; i < inverseSquareCount-1; i++)
- for (IndexType j = 1; j < inverseSquareCount-1; j++)
+ for (IndexType i = 0; i < inverseSquareCount-1; i++)
+ for (IndexType j = 0; j < inverseSquareCount-1; j++)
{
if (i - step * tau * (inverseSquareCount/this->schemeSize) * this -> speedX < 0.5*inverseSquareCount &&
j - step * tau * (inverseSquareCount/this->schemeSize) * this -> speedY < 0.5*inverseSquareCount)
diff --git a/examples/inviscid-flow/1d/EulerPressureGetter.h b/examples/inviscid-flow/1d/EulerPressureGetter.h
index 8b4fe917ec8fcebe56db23b736bbd3622555f073..b423baf0d61694b2da8487c10c16884b0e045a3d 100644
--- a/examples/inviscid-flow/1d/EulerPressureGetter.h
+++ b/examples/inviscid-flow/1d/EulerPressureGetter.h
@@ -39,7 +39,9 @@ class EulerPressureGetter
__cuda_callable__
Real operator[]( const IndexType& idx ) const
{
- if (this->rho[ idx ]==0) return 0; else return ( this->gamma - 1.0 ) * ( this->energy[ idx ] - 0.5 * this->rhoVel[ idx ] * this->rhoVel[ idx ] / this->rho[ idx ]);
+ //if (this->rho[ idx ]==0) return 0; else return ( this->gamma - 1.0 ) * ( this->energy[ idx ] - 0.5 * this->rhoVel[ idx ] * this->rhoVel[ idx ] / this->rho[ idx ]);
+ return ( this->gamma - 1.0 ) * this->energy[ idx ] * this->rho[ idx ];
+
}
diff --git a/examples/inviscid-flow/1d/euler.h b/examples/inviscid-flow/1d/euler.h
index ec652677eedecc364b8c7fb707f4ee5b22c6bce8..bd836e7cf21d4c720d797783fb519ac4b5ffe94e 100644
--- a/examples/inviscid-flow/1d/euler.h
+++ b/examples/inviscid-flow/1d/euler.h
@@ -7,6 +7,7 @@
#include "eulerProblem.h"
#include "LaxFridrichs.h"
#include "tnlMyMixedBoundaryConditions.h"
+#include "tnlMyNeumannBoundaryConditions.h"
#include "eulerRhs.h"
#include "eulerBuildConfigTag.h"
@@ -33,6 +34,7 @@ template< typename ConfigTag >class eulerConfig
config.addEntryEnum< tnlString >( "dirichlet" );
config.addEntryEnum< tnlString >( "neumann" );
config.addEntryEnum< tnlString >( "mymixed" );
+ config.addEntryEnum< tnlString >( "myneumann" );
config.addEntry< double >( "boundary-conditions-constant", "This sets a value in case of the constant boundary conditions." );
config.addEntry< double >( "left-density", "This sets a value of left density." );
config.addEntry< double >( "left-velocity", "This sets a value of left velocity." );
@@ -107,6 +109,13 @@ class eulerSetter
SolverStarter solverStarter;
return solverStarter.template run< Problem >( parameters );
}
+ if( boundaryConditionsType == "myneumann" )
+ {
+ typedef tnlMyNeumannBoundaryConditions< MeshType, MeshFunction, MeshType::getMeshDimensions(), Real, Index > BoundaryConditions;
+ typedef eulerProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
+ SolverStarter solverStarter;
+ return solverStarter.template run< Problem >( parameters );
+ }
if( boundaryConditionsType == "neumann" )
{
typedef tnlNeumannBoundaryConditions< MeshType, MeshFunction, Real, Index > BoundaryConditions;
diff --git a/examples/inviscid-flow/1d/eulerProblem_impl.h b/examples/inviscid-flow/1d/eulerProblem_impl.h
index 8a00cd69fcabb981382c2d4245e68c8b79e41a44..8bb06384e3c349857aff7aa8760bb787803a05e1 100644
--- a/examples/inviscid-flow/1d/eulerProblem_impl.h
+++ b/examples/inviscid-flow/1d/eulerProblem_impl.h
@@ -99,29 +99,27 @@ setInitialCondition( const tnlParameterContainer& parameters,
DofVectorType& dofs,
MeshDependentDataType& meshDependentData )
{
- cout << endl << "get conditions from CML";
typedef typename MeshType::Cell Cell;
this->gamma = parameters.getParameter< RealType >( "gamma" );
RealType rhoL = parameters.getParameter< RealType >( "left-density" );
RealType velL = parameters.getParameter< RealType >( "left-velocity" );
RealType preL = parameters.getParameter< RealType >( "left-pressure" );
- RealType eL = ( preL / (gamma - 1) ) + 0.5 * rhoL * velL * velL;
+ RealType eL = ( preL / ( rhoL * (gamma - 1) ) );
+ //RealType eL = ( preL / (gamma - 1) ) + 0.5 * rhoL * velL * velL;
RealType rhoR = parameters.getParameter< RealType >( "right-density" );
RealType velR = parameters.getParameter< RealType >( "right-velocity" );
RealType preR = parameters.getParameter< RealType >( "right-pressure" );
- RealType eR = ( preR / (gamma - 1) ) + 0.5 * rhoR * velR * velR;
+ RealType eR = ( preR / ( rhoR * (gamma - 1) ) );
+ //RealType eR = ( preR / (gamma - 1) ) + 0.5 * rhoR * velR * velR;
RealType x0 = parameters.getParameter< RealType >( "riemann-border" );
- cout <<endl << gamma << " " << rhoL << " " << velL << " " << preL << " " << eL << " " << rhoR << " " << velR << " " << preR << " " << eR << " " << x0 << " " << gamma << endl;
int count = mesh.template getEntitiesCount< Cell >();
-cout << count << endl;
uRho.bind(mesh, dofs, 0);
uRhoVelocity.bind(mesh, dofs, count);
uEnergy.bind(mesh, dofs, 2 * count);
tnlVector < RealType, DeviceType, IndexType > data;
data.setSize(2*count);
velocity.bind( mesh, data, 0);
- pressure.bind( mesh, data, count );
- cout << endl << "set conditions from CML"<< endl;
+ pressure.bind( mesh, data, count );
for(IndexType i = 0; i < count; i++)
if (i < x0 * count )
{
@@ -139,13 +137,6 @@ cout << count << endl;
velocity[i] = velR;
pressure[i] = preR;
};
- cout << "dofs = " << dofs << endl;
- cout << "velocity = " << velocity.getData() << endl;
- cout << "pressure = " << pressure.getData() << endl;
-
- getchar();
-
-
/*
const tnlString& initialConditionFile = parameters.getParameter< tnlString >( "initial-condition" );
if( ! dofs.load( initialConditionFile ) )
diff --git a/examples/inviscid-flow/1d/tnlMyNeumannBoundaryConditions.h b/examples/inviscid-flow/1d/tnlMyNeumannBoundaryConditions.h
new file mode 100644
index 0000000000000000000000000000000000000000..fd41f4a79a574185a7d586294179295f500b392c
--- /dev/null
+++ b/examples/inviscid-flow/1d/tnlMyNeumannBoundaryConditions.h
@@ -0,0 +1,140 @@
+/*** coppied and changed
+/***************************************************************************
+ tnlMyNeumannBoundaryConditions.h - description
+ -------------------
+ begin : Nov 17, 2014
+ copyright : (C) 2014 by oberhuber
+ email : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ ***************************************************************************/
+
+#ifndef TNLMyNeumannBoundaryConditions_H_
+#define TNLMyNeumannBoundaryConditions_H_
+
+#include <operators/tnlOperator.h>
+#include <functions/tnlConstantFunction.h>
+#include <functions/tnlFunctionAdapter.h>
+
+template< typename Mesh,
+ typename Function = tnlConstantFunction< Mesh::getMeshDimensions(), typename Mesh::RealType >,
+ int MeshEntitiesDimensions = Mesh::getMeshDimensions(),
+ typename Real = typename Mesh::RealType,
+ typename Index = typename Mesh::IndexType >
+class tnlMyNeumannBoundaryConditions
+: public tnlOperator< Mesh,
+ MeshBoundaryDomain,
+ MeshEntitiesDimensions,
+ MeshEntitiesDimensions,
+ Real,
+ Index >
+{
+ public:
+
+ typedef Mesh MeshType;
+ typedef Function FunctionType;
+ typedef Real RealType;
+ typedef typename MeshType::DeviceType DeviceType;
+ typedef Index IndexType;
+
+ typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
+ typedef typename MeshType::VertexType VertexType;
+
+ static constexpr int getMeshDimensions() { return MeshType::meshDimensions; }
+
+ static void configSetup( tnlConfigDescription& config,
+ const tnlString& prefix = "" )
+ {
+ Function::configSetup( config, prefix );
+ }
+
+ bool setup( const tnlParameterContainer& parameters,
+ const tnlString& prefix = "" )
+ {
+ return this->function.setup( parameters, prefix );
+ }
+
+ void setFunction( const Function& function )
+ {
+ this->function = function;
+ }
+
+ Function& getFunction()
+ {
+ return this->function;
+ }
+
+ const Function& getFunction() const
+ {
+ return this->function;
+ }
+
+ template< typename EntityType,
+ typename MeshFunction >
+ __cuda_callable__
+ const RealType operator()( const MeshFunction& u,
+ const EntityType& entity,
+ const RealType& time = 0 ) const
+ {
+ const MeshType& mesh = entity.getMesh();
+ const auto& neighbourEntities = entity.getNeighbourEntities();
+ const IndexType& index = entity.getIndex();
+ if( entity.getCoordinates().x() == 0 )
+ return u[ neighbourEntities.template getEntityIndex< 1 >() ];
+ else
+ return u[ neighbourEntities.template getEntityIndex< -1 >() ];
+ //tady se asi delaji okrajove podminky
+ //static_assert( EntityType::getDimensions() == MeshEntitiesDimensions, "Wrong mesh entity dimensions." );
+ }
+
+ template< typename EntityType >
+ __cuda_callable__
+ IndexType getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const EntityType& entity ) const
+ {
+ return 1;
+ }
+
+ template< typename PreimageFunction,
+ typename MeshEntity,
+ typename Matrix,
+ typename Vector >
+ __cuda_callable__
+ void setMatrixElements( const PreimageFunction& u,
+ const MeshEntity& entity,
+ const RealType& time,
+ const RealType& tau,
+ Matrix& matrix,
+ Vector& b ) const
+ {
+ typename Matrix::MatrixRow matrixRow = matrix.getRow( entity.getIndex() );
+ const IndexType& index = entity.getIndex();
+ matrixRow.setElement( 0, index, 1.0 );
+ b[ index ] = tnlFunctionAdapter< MeshType, Function >::getValue( this->function, entity, time );
+ }
+
+
+ protected:
+
+ Function function;
+
+ //static_assert( Device::DeviceType == Function::Device::DeviceType );
+};
+
+template< typename Mesh,
+ typename Function >
+ostream& operator << ( ostream& str, const tnlMyNeumannBoundaryConditions< Mesh, Function >& bc )
+{
+ str << "MyMixed boundary conditions: vector = " << bc.getVector();
+ return str;
+}
+
+#endif /* TNLMyNeumannBoundaryConditions_H_ */
diff --git a/examples/inviscid-flow/2d/EulerPressureGetter.h b/examples/inviscid-flow/2d/EulerPressureGetter.h
index 3af2a4a9c1b03f02ca5c820a61b7bb6f21e33e8c..131db48b9f814f1163a8a484f5bc66f7faa98a3b 100644
--- a/examples/inviscid-flow/2d/EulerPressureGetter.h
+++ b/examples/inviscid-flow/2d/EulerPressureGetter.h
@@ -40,9 +40,10 @@ class EulerPressureGetter
__cuda_callable__
Real operator[]( const IndexType& idx ) const
{
- if (this->rho[ idx ]==0) return 0;
+/* if (this->rho[ idx ]==0) return 0;
else return ( this->gamma - 1.0 ) * ( this->energy[ idx ] - 0.5 * this->rho[ idx ] *
( this->rhoVelX[ idx ] / this->rho[ idx ] + this->rhoVelY[ idx ] / this->rho[ idx ]) );
+*/ return ( this->gamma - 1.0 ) * ( this->energy[ idx ] * this->rho[ idx ] );
}
diff --git a/examples/inviscid-flow/2d/LaxFridrichsMomentumY_impl.h b/examples/inviscid-flow/2d/LaxFridrichsMomentumY_impl.h
index 4574a8d3d2cef8bca3215e093d82d2d9eba560ab..588fab9e048ec32d035b34eb966e149f8e3a83c0 100644
--- a/examples/inviscid-flow/2d/LaxFridrichsMomentumY_impl.h
+++ b/examples/inviscid-flow/2d/LaxFridrichsMomentumY_impl.h
@@ -153,9 +153,9 @@ operator()( const MeshFunction& u,
const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >();
const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >();
return (0.25 / this->tau) * ( u[ west ] + u[ east ] + u[ south ] + u[ north ] - 4.0 * u[ center ] )
- - 0.5 * hyInverse * (( u[ east ] * this->velocityX[ east ] )
+ - 0.5 * hxInverse * (( u[ east ] * this->velocityX[ east ] )
-( u[ west ] * this->velocityX[ west ] ))
- - 0.5 * hxInverse * (( u[ north ] * this->velocityY[ north ] + this->pressure[ north ] )
+ - 0.5 * hyInverse * (( u[ north ] * this->velocityY[ north ] + this->pressure[ north ] )
-( u[ south ] * this->velocityY[ south ] + this->pressure[ south ]));
}
diff --git a/examples/inviscid-flow/2d/euler.h b/examples/inviscid-flow/2d/euler.h
index f1466b9e3254e05c9a9f383eafe4455a1516f9ad..dc29fe2e0cc23c2ee00d0865517755a0ef7b6fdc 100644
--- a/examples/inviscid-flow/2d/euler.h
+++ b/examples/inviscid-flow/2d/euler.h
@@ -9,6 +9,7 @@
#include "eulerRhs.h"
#include "eulerBuildConfigTag.h"
#include "tnlMyMixedBoundaryConditions.h"
+#include "tnlMyNeumannBoundaryConditions.h"
typedef eulerBuildConfigTag BuildConfig;
@@ -32,6 +33,7 @@ template< typename ConfigTag >class eulerConfig
config.addEntryEnum< tnlString >( "dirichlet" );
config.addEntryEnum< tnlString >( "neumann" );
config.addEntryEnum< tnlString >( "mymixed" );
+ config.addEntryEnum< tnlString >( "myneumann" );
config.addEntry< double >( "boundary-conditions-constant", "This sets a value in case of the constant boundary conditions." );
config.addEntry< double >( "left-up-density", "This sets a value of left up density." );
config.addEntry< double >( "left-up-velocityX", "This sets a value of left up x velocity." );
@@ -116,6 +118,13 @@ class eulerSetter
SolverStarter solverStarter;
return solverStarter.template run< Problem >( parameters );
}
+ if( boundaryConditionsType == "myneumann" )
+ {
+ typedef tnlMyNeumannBoundaryConditions< MeshType, MeshFunction, MeshType::getMeshDimensions(), Real, Index > BoundaryConditions;
+ typedef eulerProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
+ SolverStarter solverStarter;
+ return solverStarter.template run< Problem >( parameters );
+ }
if( boundaryConditionsType == "neumann" )
{
typedef tnlNeumannBoundaryConditions< MeshType, MeshFunction, Real, Index > BoundaryConditions;
diff --git a/examples/inviscid-flow/2d/eulerProblem_impl.h b/examples/inviscid-flow/2d/eulerProblem_impl.h
index 60bca6d92ccf8d298dc266bf2a0f0d7d14c80520..ca84918abcb7ee2219d4c74ea7c412c944d946e5 100644
--- a/examples/inviscid-flow/2d/eulerProblem_impl.h
+++ b/examples/inviscid-flow/2d/eulerProblem_impl.h
@@ -108,22 +108,26 @@ setInitialCondition( const tnlParameterContainer& parameters,
RealType velLuX = parameters.getParameter< RealType >( "left-up-velocityX" );
RealType velLuY = parameters.getParameter< RealType >( "left-up-velocityY" );
RealType preLu = parameters.getParameter< RealType >( "left-up-pressure" );
- RealType eLu = ( preLu / (gamma - 1) ) + 0.5 * rhoLu * pow(velLuX,2)+pow(velLuY,2);
+ RealType eLu = ( preLu / ( rhoLu * (gamma - 1) ) );
+ //RealType eLu = ( preLu / (gamma - 1) ) + 0.5 * rhoLu * pow(velLuX,2)+pow(velLuY,2);
RealType rhoLd = parameters.getParameter< RealType >( "left-down-density" );
RealType velLdX = parameters.getParameter< RealType >( "left-down-velocityX" );
RealType velLdY = parameters.getParameter< RealType >( "left-down-velocityY" );
RealType preLd = parameters.getParameter< RealType >( "left-down-pressure" );
- RealType eLd = ( preLd / (gamma - 1) ) + 0.5 * rhoLd * pow(velLdX,2)+pow(velLdY,2);
+ RealType eLd = ( preLd / ( rhoLd * (gamma - 1) ) );
+ //RealType eLd = ( preLd / (gamma - 1) ) + 0.5 * rhoLd * pow(velLdX,2)+pow(velLdY,2);
RealType rhoRu = parameters.getParameter< RealType >( "right-up-density" );
RealType velRuX = parameters.getParameter< RealType >( "right-up-velocityX" );
RealType velRuY = parameters.getParameter< RealType >( "right-up-velocityY" );
RealType preRu = parameters.getParameter< RealType >( "right-up-pressure" );
- RealType eRu = ( preRu / (gamma - 1) ) + 0.5 * rhoRu * pow(velRuX,2)+pow(velRuY,2);
+ RealType eRu = ( preRu / ( rhoRu * (gamma - 1) ) );
+ //RealType eRu = ( preRu / (gamma - 1) ) + 0.5 * rhoRu * pow(velRuX,2)+pow(velRuY,2);
RealType rhoRd = parameters.getParameter< RealType >( "right-down-density" );
RealType velRdX = parameters.getParameter< RealType >( "right-down-velocityX" );
RealType velRdY = parameters.getParameter< RealType >( "right-down-velocityY" );
RealType preRd = parameters.getParameter< RealType >( "right-down-pressure" );
- RealType eRd = ( preRd / (gamma - 1) ) + 0.5 * rhoRd * pow(velRdX,2)+pow(velRdY,2);
+ RealType eRd = ( preRd / ( rhoRd * (gamma - 1) ) );
+ //RealType eRd = ( preRd / (gamma - 1) ) + 0.5 * rhoRd * pow(velRdX,2)+pow(velRdY,2);
RealType x0 = parameters.getParameter< RealType >( "riemann-border" );
int size = mesh.template getEntitiesCount< Cell >();
uRho.bind(mesh, dofs, 0);
@@ -138,7 +142,7 @@ setInitialCondition( const tnlParameterContainer& parameters,
velocityY.bind(mesh, data, 3*size);
for(IndexType j = 0; j < sqrt(size); j++)
for(IndexType i = 0; i < sqrt(size); i++)
- if ((i < x0 * sqrt(size))&&(j < x0 * sqrt(size)) )
+ if ((i <= x0 * sqrt(size))&&(j <= x0 * sqrt(size)) )
{
uRho[j*sqrt(size)+i] = rhoLd;
uRhoVelocityX[j*sqrt(size)+i] = rhoLd * velLdX;
@@ -150,7 +154,7 @@ setInitialCondition( const tnlParameterContainer& parameters,
pressure[j*sqrt(size)+i] = preLd;
}
else
- if ((i >= x0 * sqrt(size))&&(j < x0 * sqrt(size)) )
+ if ((i <= x0 * sqrt(size))&&(j > x0 * sqrt(size)) )
{
uRho[j*sqrt(size)+i] = rhoLu;
uRhoVelocityX[j*sqrt(size)+i] = rhoLu * velLuX;
@@ -162,7 +166,7 @@ setInitialCondition( const tnlParameterContainer& parameters,
pressure[j*sqrt(size)+i] = preLu;
}
else
- if ((i >= x0 * sqrt(size))&&(j >= x0 * sqrt(size)) )
+ if ((i > x0 * sqrt(size))&&(j > x0 * sqrt(size)) )
{
uRho[j*sqrt(size)+i] = rhoRu;
uRhoVelocityX[j*sqrt(size)+i] = rhoRu * velRuX;
diff --git a/examples/inviscid-flow/2d/tnlMyNeumannBoundaryConditions.h b/examples/inviscid-flow/2d/tnlMyNeumannBoundaryConditions.h
new file mode 100644
index 0000000000000000000000000000000000000000..fdc55c33629aad9516a10a3b953d1cf5d87d4784
--- /dev/null
+++ b/examples/inviscid-flow/2d/tnlMyNeumannBoundaryConditions.h
@@ -0,0 +1,152 @@
+/*** coppied and changed
+/***************************************************************************
+ tnlMyNeumannBoundaryConditions.h - description
+ -------------------
+ begin : Nov 17, 2014
+ copyright : (C) 2014 by oberhuber
+ email : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ ***************************************************************************/
+
+#ifndef TNLMyNeumannBOUNDARYCONDITIONS_H_
+#define TNLMyNeumannBOUNDARYCONDITIONS_H_
+
+#include <operators/tnlOperator.h>
+#include <functions/tnlConstantFunction.h>
+#include <functions/tnlFunctionAdapter.h>
+
+template< typename Mesh,
+ typename Function = tnlConstantFunction< Mesh::getMeshDimensions(), typename Mesh::RealType >,
+ int MeshEntitiesDimensions = Mesh::getMeshDimensions(),
+ typename Real = typename Mesh::RealType,
+ typename Index = typename Mesh::IndexType >
+class tnlMyNeumannBoundaryConditions
+: public tnlOperator< Mesh,
+ MeshBoundaryDomain,
+ MeshEntitiesDimensions,
+ MeshEntitiesDimensions,
+ Real,
+ Index >
+{
+ public:
+
+ typedef Mesh MeshType;
+ typedef Function FunctionType;
+ typedef Real RealType;
+ typedef typename MeshType::DeviceType DeviceType;
+ typedef Index IndexType;
+
+ typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
+ typedef typename MeshType::VertexType VertexType;
+
+ static constexpr int getMeshDimensions() { return MeshType::meshDimensions; }
+
+ static void configSetup( tnlConfigDescription& config,
+ const tnlString& prefix = "" )
+ {
+ Function::configSetup( config, prefix );
+ }
+
+ bool setup( const tnlParameterContainer& parameters,
+ const tnlString& prefix = "" )
+ {
+ return this->function.setup( parameters, prefix );
+ }
+
+ void setFunction( const Function& function )
+ {
+ this->function = function;
+ }
+
+ void setX0( const RealType& x0 )
+ {
+ this->x0 = x0;
+ }
+
+ Function& getFunction()
+ {
+ return this->function;
+ }
+
+ const Function& getFunction() const
+ {
+ return this->function;
+ }
+
+ template< typename EntityType,
+ typename MeshFunction >
+ __cuda_callable__
+ const RealType operator()( const MeshFunction& u,
+ const EntityType& entity,
+ const RealType& time = 0 ) const
+ {
+ const MeshType& mesh = entity.getMesh();
+ const auto& neighbourEntities = entity.getNeighbourEntities();
+ typedef typename MeshType::Cell Cell;
+ int count = mesh.template getEntitiesCount< Cell >();
+ count = sqrt(count);
+ if( entity.getCoordinates().x() == 0 )
+ return u[ neighbourEntities.template getEntityIndex< 1, 0 >() ];
+ else if( entity.getCoordinates().x() == count-1 )
+ return u[ neighbourEntities.template getEntityIndex< -1, 0 >() ];
+ else if( entity.getCoordinates().y() == 0 )
+ return u[ neighbourEntities.template getEntityIndex< 0, 1 >() ];
+ else return u[ neighbourEntities.template getEntityIndex< 0, -1 >() ];
+ //tady se asi delaji okrajove podminky
+ //static_assert( EntityType::getDimensions() == MeshEntitiesDimensions, "Wrong mesh entity dimensions." );
+ }
+
+ template< typename EntityType >
+ __cuda_callable__
+ IndexType getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const EntityType& entity ) const
+ {
+ return 1;
+ }
+
+ template< typename PreimageFunction,
+ typename MeshEntity,
+ typename Matrix,
+ typename Vector >
+ __cuda_callable__
+ void setMatrixElements( const PreimageFunction& u,
+ const MeshEntity& entity,
+ const RealType& time,
+ const RealType& tau,
+ Matrix& matrix,
+ Vector& b ) const
+ {
+ typename Matrix::MatrixRow matrixRow = matrix.getRow( entity.getIndex() );
+ const IndexType& index = entity.getIndex();
+ matrixRow.setElement( 0, index, 1.0 );
+ b[ index ] = tnlFunctionAdapter< MeshType, Function >::getValue( this->function, entity, time );
+ }
+
+
+ protected:
+
+ Function function;
+
+ RealType x0 = 0.3;
+
+ //static_assert( Device::DeviceType == Function::Device::DeviceType );
+};
+
+template< typename Mesh,
+ typename Function >
+ostream& operator << ( ostream& str, const tnlMyNeumannBoundaryConditions< Mesh, Function >& bc )
+{
+ str << "MyNeumann boundary conditions: vector = " << bc.getVector();
+ return str;
+}
+
+#endif /* TNLMyNeumannBOUNDARYCONDITIONS_H_ */