From e39433287397845fa07038c80a45a34e5dd9faa5 Mon Sep 17 00:00:00 2001
From: Ondrej Szekely <ondra.szekely@gmail.com>
Date: Mon, 8 Jun 2015 08:48:05 +0200
Subject: [PATCH] Adding Finite Volume method for 2D mesh
---
.../tnl-mean-curvature-flow-eoc.h | 8 +-
.../CMakeLists.txt | 4 +-
.../tnlFiniteVolumeNonlinearOperator.h | 187 +++++
.../tnlFiniteVolumeNonlinearOperator_impl.h | 307 +++++++++
src/operators/operator-Q/CMakeLists.txt | 4 +-
.../operator-Q/tnlFiniteVolumeOperatorQ.h | 387 +++++++++++
.../tnlFiniteVolumeOperatorQ_impl.h | 636 ++++++++++++++++++
7 files changed, 1528 insertions(+), 5 deletions(-)
create mode 100644 src/operators/diffusion/nonlinear-diffusion-operators/tnlFiniteVolumeNonlinearOperator.h
create mode 100644 src/operators/diffusion/nonlinear-diffusion-operators/tnlFiniteVolumeNonlinearOperator_impl.h
create mode 100644 src/operators/operator-Q/tnlFiniteVolumeOperatorQ.h
create mode 100644 src/operators/operator-Q/tnlFiniteVolumeOperatorQ_impl.h
diff --git a/examples/mean-curvature-flow/tnl-mean-curvature-flow-eoc.h b/examples/mean-curvature-flow/tnl-mean-curvature-flow-eoc.h
index 483bff18cf..eb79ecc722 100644
--- a/examples/mean-curvature-flow/tnl-mean-curvature-flow-eoc.h
+++ b/examples/mean-curvature-flow/tnl-mean-curvature-flow-eoc.h
@@ -29,8 +29,10 @@
#include <operators/diffusion/tnlExactNonlinearDiffusion.h>
#include <operators/diffusion/tnlNonlinearDiffusion.h>
#include <operators/operator-Q/tnlOneSideDiffOperatorQ.h>
+#include <operators/operator-Q/tnlFiniteVolumeOperatorQ.h>
#include <operators/diffusion/tnlExactNonlinearDiffusion.h>
#include <operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator.h>
+#include <operators/diffusion/nonlinear-diffusion-operators/tnlFiniteVolumeNonlinearOperator.h>
#include <operators/operator-Q/tnlExactOperatorQ.h>
//typedef tnlDefaultConfigTag BuildConfig;
@@ -45,7 +47,7 @@ class meanCurvatureFlowEocConfig
config.addDelimiter( "Mean Curvature Flow EOC settings:" );
config.addEntry< double >( "eps", "This sets a eps in operator Q.", 1.0 );
config.addDelimiter( "Tests setting::" );
- tnlTestFunction< 2, double >::configSetup( config );
+ tnlTestFunction< 3, double >::configSetup( config );
}
};
@@ -68,8 +70,8 @@ class meanCurvatureFlowEocSetter
static bool run( const tnlParameterContainer& parameters )
{
enum { Dimensions = MeshType::Dimensions };
- typedef tnlOneSideDiffOperatorQ<MeshType, Real, Index, 0> OperatorQ;
- typedef tnlOneSideDiffNonlinearOperator<MeshType, OperatorQ, Real, Index > NonlinearOperator;
+ typedef tnlFiniteVolumeOperatorQ<MeshType, Real, Index, 0> OperatorQ;
+ typedef tnlFiniteVolumeNonlinearOperator<MeshType, OperatorQ, Real, Index > NonlinearOperator;
typedef tnlNonlinearDiffusion< MeshType, NonlinearOperator, Real, Index > ApproximateOperator;
typedef tnlExactNonlinearDiffusion< tnlExactOperatorQ<Dimensions>, Dimensions > ExactOperator;
typedef tnlTestFunction< MeshType::Dimensions, Real, Device > TestFunction;
diff --git a/src/operators/diffusion/nonlinear-diffusion-operators/CMakeLists.txt b/src/operators/diffusion/nonlinear-diffusion-operators/CMakeLists.txt
index 0962ac61ee..41886e78d4 100755
--- a/src/operators/diffusion/nonlinear-diffusion-operators/CMakeLists.txt
+++ b/src/operators/diffusion/nonlinear-diffusion-operators/CMakeLists.txt
@@ -1,4 +1,6 @@
SET( headers tnlOneSideDiffNonlinearOperator.h
- tnlOneSideDiffNonlinearOperator_impl.h )
+ tnlOneSideDiffNonlinearOperator_impl.h
+ tnlFiniteVolumeNonlinearOperator.h
+ tnlFiniteVolumeNonlinearOperator_impl.h )
INSTALL( FILES ${headers} DESTINATION include/tnl-${tnlVersion}/operators/diffusion/nonlinear-diffusion-operators )
diff --git a/src/operators/diffusion/nonlinear-diffusion-operators/tnlFiniteVolumeNonlinearOperator.h b/src/operators/diffusion/nonlinear-diffusion-operators/tnlFiniteVolumeNonlinearOperator.h
new file mode 100644
index 0000000000..82557f4f56
--- /dev/null
+++ b/src/operators/diffusion/nonlinear-diffusion-operators/tnlFiniteVolumeNonlinearOperator.h
@@ -0,0 +1,187 @@
+#ifndef TNLFINITEVOLUMENONLINEAROPERATOR_H
+#define TNLFINITEVOLUMENONLINEAROPERATOR_H
+
+#include <core/vectors/tnlVector.h>
+#include <mesh/tnlGrid.h>
+
+template< typename Mesh,
+ typename NonlinearDiffusionOperator,
+ typename OperatorQ,
+ typename Real = typename Mesh::RealType,
+ typename Index = typename Mesh::IndexType >
+class tnlFiniteVolumeNonlinearOperator
+{
+
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+class tnlFiniteVolumeNonlinearOperator< tnlGrid< 1,MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >
+{
+ public:
+
+ typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef OperatorQ OperatorQType;
+
+ static tnlString getType();
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const RealType& time) const;
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const;
+
+ template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const;
+
+ public:
+
+ OperatorQ operatorQ;
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+class tnlFiniteVolumeNonlinearOperator< tnlGrid< 2, MeshReal, Device, MeshIndex >, OperatorQ, 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 OperatorQ OperatorQType;
+
+
+ static tnlString getType();
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const RealType& time) const;
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const;
+
+ template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const;
+
+ public:
+
+ OperatorQ operatorQ;
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+class tnlFiniteVolumeNonlinearOperator< tnlGrid< 3, MeshReal, Device, MeshIndex >, OperatorQ, 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 OperatorQ OperatorQType;
+
+ static tnlString getType();
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const RealType& time) const;
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const;
+
+ template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const;
+
+ public:
+
+ OperatorQ operatorQ;
+};
+
+
+#include "tnlFiniteVolumeNonlinearOperator_impl.h"
+
+
+#endif /* TNLFINITEVOLUMENONLINEAROPERATOR_H */
diff --git a/src/operators/diffusion/nonlinear-diffusion-operators/tnlFiniteVolumeNonlinearOperator_impl.h b/src/operators/diffusion/nonlinear-diffusion-operators/tnlFiniteVolumeNonlinearOperator_impl.h
new file mode 100644
index 0000000000..730adff00a
--- /dev/null
+++ b/src/operators/diffusion/nonlinear-diffusion-operators/tnlFiniteVolumeNonlinearOperator_impl.h
@@ -0,0 +1,307 @@
+
+#ifndef TNLFINITEVOLUMENONLINEAROPERATOR__IMPL_H
+#define TNLFINITEVOLUMENONLINEAROPERATOR__IMPL_H
+
+#include "tnlFiniteVolumeNonlinearOperator.h"
+
+#include <mesh/tnlGrid.h>
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+tnlString
+tnlFiniteVolumeNonlinearOperator< tnlGrid< 1, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+getType()
+{
+ return tnlString( "tnlFiniteVolumeNonlinearOperator< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + ", " +
+ OperatorQ::getType() + " >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlFiniteVolumeNonlinearOperator< tnlGrid< 1, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const
+{
+ return 0.0;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+Index
+tnlFiniteVolumeNonlinearOperator< tnlGrid< 1, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const
+{
+ return 1;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+void
+tnlFiniteVolumeNonlinearOperator< tnlGrid< 1, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const
+{
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+tnlString
+tnlFiniteVolumeNonlinearOperator< tnlGrid< 2, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+getType()
+{
+ return tnlString( "tnlFiniteVolumeNonlinearOperator< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + ", " +
+ OperatorQ::getType() + " >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlFiniteVolumeNonlinearOperator< tnlGrid< 2, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const
+{
+ return operatorQ.getValue(mesh,cellIndex, coordinates, u, time)*( (u[ mesh.template getCellNextToCell< 1,0 >( cellIndex ) ] - u[ cellIndex ])
+ * mesh.getHxSquareInverse() / operatorQ.getValue(mesh,cellIndex, coordinates, u, time, 1 )
+ + ( u[ mesh.template getCellNextToCell< 0,1 >( cellIndex ) ] - u[ cellIndex ]) * mesh.getHySquareInverse()/
+ operatorQ.getValue(mesh,cellIndex, coordinates, u, time, 0, 1 )
+ - ( - u[ mesh.template getCellNextToCell< -1,0 >( cellIndex ) ] + u[ cellIndex ])
+ * mesh.getHxSquareInverse()/operatorQ.getValue(mesh,cellIndex, coordinates, u, time, -1)
+ -( - u[ mesh.template getCellNextToCell< 0,-1 >( cellIndex ) ] + u[ cellIndex ]) * mesh.getHySquareInverse()
+ /operatorQ.getValue(mesh,cellIndex,coordinates, u, time, 0, -1) );
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+Index
+tnlFiniteVolumeNonlinearOperator< tnlGrid< 2, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const
+{
+ return 5;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+void
+tnlFiniteVolumeNonlinearOperator< tnlGrid< 2, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const
+{
+ const RealType aCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time ) * mesh.getHySquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time, 0, -1 );
+ const RealType bCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time ) * mesh.getHxSquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time, -1 );
+ const RealType cCoef = tau * operatorQ.getValue(mesh, index, coordinates, u, time ) * ( mesh.getHxSquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time, 1 ) + mesh.getHySquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time, 0, 1 )
+ + mesh.getHxSquareInverse() / operatorQ.getValue(mesh, index, coordinates, u, time, -1 ) +
+ mesh.getHySquareInverse() / operatorQ.getValue(mesh, index, coordinates, u, time, 0, -1 ) );
+ const RealType dCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time ) * mesh.getHxSquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time, 1 );
+ const RealType eCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time ) * mesh.getHySquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time, 0, 1 );
+ matrixRow.setElement( 0, mesh.template getCellNextToCell< 0,-1 >( index ), aCoef );
+ matrixRow.setElement( 1, mesh.template getCellNextToCell< -1,0 >( index ), bCoef );
+ matrixRow.setElement( 2, index , cCoef );
+ matrixRow.setElement( 3, mesh.template getCellNextToCell< 1,0 >( index ), dCoef );
+ matrixRow.setElement( 4, mesh.template getCellNextToCell< 0,1 >( index ), eCoef );
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+tnlString
+tnlFiniteVolumeNonlinearOperator< tnlGrid< 3, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+getType()
+{
+ return tnlString( "tnlFiniteVolumeNonlinearOperator< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + ", " +
+ OperatorQ::getType() + " >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlFiniteVolumeNonlinearOperator< tnlGrid< 3, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const
+{
+ return operatorQ.getValue(mesh,cellIndex, coordinates, u, time)*( (u[ mesh.template getCellNextToCell< 1,0,0 >( cellIndex ) ] - u[ cellIndex ])
+ * mesh.getHxSquareInverse() / operatorQ.getValue(mesh,cellIndex, coordinates, u, time, 1 )
+ + ( u[ mesh.template getCellNextToCell< 0,1,0 >( cellIndex ) ] - u[ cellIndex ]) * mesh.getHySquareInverse()/
+ operatorQ.getValue(mesh,cellIndex, coordinates, u, time, 0, 1 )
+ + ( u[ mesh.template getCellNextToCell< 0,0,1 >( cellIndex ) ] - u[ cellIndex ]) * mesh.getHzSquareInverse()/
+ operatorQ.getValue(mesh,cellIndex, coordinates, u, time, 0, 0, 1 )
+ - ( - u[ mesh.template getCellNextToCell< -1,0,0 >( cellIndex ) ] + u[ cellIndex ])
+ * mesh.getHxSquareInverse()/operatorQ.getValue(mesh,cellIndex, coordinates, u, time, -1)
+ -( - u[ mesh.template getCellNextToCell< 0,-1,0 >( cellIndex ) ] + u[ cellIndex ]) * mesh.getHySquareInverse()
+ /operatorQ.getValue(mesh,cellIndex,coordinates, u, time, 0, -1)
+ -( - u[ mesh.template getCellNextToCell< 0,0,-1 >( cellIndex ) ] + u[ cellIndex ]) * mesh.getHzSquareInverse()
+ /operatorQ.getValue(mesh,cellIndex,coordinates, u, time, 0, 0, -1) );
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+Index
+tnlFiniteVolumeNonlinearOperator< tnlGrid< 3, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const
+{
+ return 7;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index,
+ typename OperatorQ >
+template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+void
+tnlFiniteVolumeNonlinearOperator< tnlGrid< 3, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::
+updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const
+{
+ const RealType aCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time ) * mesh.getHzSquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time, 0, 0, -1 );
+ const RealType bCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time ) * mesh.getHySquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time, 0, -1, 0 );
+ const RealType cCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time ) * mesh.getHxSquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time, -1, 0, 0 );
+ const RealType dCoef = tau * operatorQ.getValue(mesh, index, coordinates, u, time ) * ( mesh.getHxSquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time, 1, 0, 0 ) + mesh.getHySquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time, 0, 1, 0 )
+ + mesh.getHzSquareInverse() / operatorQ.getValue(mesh, index, coordinates, u, time, 0, 0, 1 ) +
+ mesh.getHxSquareInverse() / operatorQ.getValue(mesh, index, coordinates, u, time, -1, 0, 0 )
+ + mesh.getHySquareInverse() / operatorQ.getValue(mesh, index, coordinates, u, time, 0, -1, 0 ) +
+ mesh.getHzSquareInverse() / operatorQ.getValue(mesh, index, coordinates, u, time, 0, 0, -1 ) );
+ const RealType eCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time, 1, 0, 0 ) * mesh.getHxSquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time );
+ const RealType fCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time, 0, 1, 0 ) * mesh.getHySquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time );
+ const RealType gCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time, 0, 0, 1 ) * mesh.getHzSquareInverse() /
+ operatorQ.getValue(mesh, index, coordinates, u, time );
+ matrixRow.setElement( 0, mesh.template getCellNextToCell< 0,0,-1 >( index ), aCoef );
+ matrixRow.setElement( 1, mesh.template getCellNextToCell< 0,-1,0 >( index ), bCoef );
+ matrixRow.setElement( 2, mesh.template getCellNextToCell< -1,0,0 >( index ), cCoef );
+ matrixRow.setElement( 3, index, dCoef );
+ matrixRow.setElement( 4, mesh.template getCellNextToCell< 1,0,0 >( index ), eCoef );
+ matrixRow.setElement( 5, mesh.template getCellNextToCell< 0,1,0 >( index ), fCoef );
+ matrixRow.setElement( 6, mesh.template getCellNextToCell< 0,0,1 >( index ), gCoef );
+}
+#endif /* TNLFINITEVOLUMENONLINEAROPERATOR__IMPL_H */
diff --git a/src/operators/operator-Q/CMakeLists.txt b/src/operators/operator-Q/CMakeLists.txt
index 747d01a753..192c55d961 100755
--- a/src/operators/operator-Q/CMakeLists.txt
+++ b/src/operators/operator-Q/CMakeLists.txt
@@ -1,6 +1,8 @@
SET( headers tnlExactOperatorQ.h
tnlExactOperatorQ_impl.h
tnlOneSideDiffOperatorQ.h
- tnlOneSideDiffOperatorQ_impl.h )
+ tnlOneSideDiffOperatorQ_impl.h
+ tnlFiniteVolumeOperatorQ.h
+ tnlFiniteVolumeOperatorQ_impl.h )
INSTALL( FILES ${headers} DESTINATION include/tnl-${tnlVersion}/operators/operator-Q )
diff --git a/src/operators/operator-Q/tnlFiniteVolumeOperatorQ.h b/src/operators/operator-Q/tnlFiniteVolumeOperatorQ.h
new file mode 100644
index 0000000000..a3188559e2
--- /dev/null
+++ b/src/operators/operator-Q/tnlFiniteVolumeOperatorQ.h
@@ -0,0 +1,387 @@
+#ifndef TNLFINITEVOLUMEOPERATORQ_H
+#define TNLFINITEVOLUMEOPERATORQ_H
+
+#include <core/vectors/tnlVector.h>
+#include <core/vectors/tnlSharedVector.h>
+#include <mesh/tnlGrid.h>
+
+template< typename Mesh,
+ typename Real = typename Mesh::RealType,
+ typename Index = typename Mesh::IndexType,
+ int Precomputation = 0 >
+class tnlFiniteVolumeOperatorQ
+{
+
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+class tnlFiniteVolumeOperatorQ< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index, 0 >
+{
+ public:
+
+ typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ static tnlString getType();
+
+ template< typename Vector >
+ IndexType bind( Vector& u)
+ { return 0; }
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void update( const MeshType& mesh, const RealType& time )
+ {}
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType& dx = 0,
+ const IndexType& dy = 0,
+ const IndexType& dz = 0 ) const;
+
+ bool setEps(const Real& eps);
+
+ private:
+
+ template< typename Vector, int AxeX = 0, int AxeY = 0, int AxeZ = 0 >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real boundaryDerivative( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType& dx = 0,
+ const IndexType& dy = 0,
+ const IndexType& dz = 0 ) const;
+
+ RealType eps;
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+class tnlFiniteVolumeOperatorQ< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index, 0 >
+{
+ public:
+
+ typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ static tnlString getType();
+
+ template< typename Vector >
+ IndexType bind( Vector& u)
+ { return 0; }
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void update( const MeshType& mesh, const RealType& time )
+ {}
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType& dx = 0,
+ const IndexType& dy = 0,
+ const IndexType& dz = 0 ) const;
+
+ bool setEps(const Real& eps);
+
+ private:
+
+ template< typename Vector, int AxeX = 0, int AxeY = 0, int AxeZ = 0 >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real boundaryDerivative( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType& dx = 0,
+ const IndexType& dy = 0,
+ const IndexType& dz = 0 ) const;
+
+ RealType eps;
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+class tnlFiniteVolumeOperatorQ< tnlGrid< 3,MeshReal, Device, MeshIndex >, Real, Index, 0 >
+{
+ public:
+
+ typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ static tnlString getType();
+
+ template< typename Vector >
+ IndexType bind( Vector& u)
+ { return 0; }
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void update( const MeshType& mesh, const RealType& time )
+ {}
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType& dx = 0,
+ const IndexType& dy = 0,
+ const IndexType& dz = 0 ) const;
+
+ bool setEps(const Real& eps);
+
+ private:
+
+ template< typename Vector, int AxeX = 0, int AxeY = 0, int AxeZ = 0 >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real boundaryDerivative( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType& dx = 0,
+ const IndexType& dy = 0,
+ const IndexType& dz = 0 ) const;
+
+ RealType eps;
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+class tnlFiniteVolumeOperatorQ< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index, 1 >
+{
+ public:
+
+ typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ static tnlString getType();
+
+ template< typename Vector >
+ Index bind( Vector& u);
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void update( const MeshType& mesh, const RealType& time );
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType& dx = 0,
+ const IndexType& dy = 0,
+ const IndexType& dz = 0 ) const;
+
+ bool setEps(const Real& eps);
+
+ private:
+
+ template< typename Vector, int AxeX = 0, int AxeY = 0, int AxeZ = 0 >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real boundaryDerivative( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType& dx = 0,
+ const IndexType& dy = 0,
+ const IndexType& dz = 0 ) const;
+
+ tnlSharedVector< RealType, DeviceType, IndexType > u;
+ tnlVector< RealType, DeviceType, IndexType> q;
+ RealType eps;
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+class tnlFiniteVolumeOperatorQ< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index, 1 >
+{
+ public:
+
+ typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef tnlSharedVector< RealType, DeviceType, IndexType > DofVectorType;
+
+ static tnlString getType();
+
+ template< typename Vector >
+ Index bind( Vector& u);
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void update( const MeshType& mesh, const RealType& time );
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType& dx = 0,
+ const IndexType& dy = 0,
+ const IndexType& dz = 0 ) const;
+
+ bool setEps(const Real& eps);
+
+ private:
+
+ template< typename Vector, int AxeX = 0, int AxeY = 0, int AxeZ = 0 >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real boundaryDerivative( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType& dx = 0,
+ const IndexType& dy = 0,
+ const IndexType& dz = 0 ) const;
+
+ tnlSharedVector< RealType, DeviceType, IndexType > u;
+ tnlVector< RealType, DeviceType, IndexType> q;
+ RealType eps;
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+class tnlFiniteVolumeOperatorQ< tnlGrid< 3,MeshReal, Device, MeshIndex >, Real, Index, 1 >
+{
+ public:
+
+ typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ static tnlString getType();
+
+ template< typename Vector >
+ Index bind( Vector& u);
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void update( const MeshType& mesh, const RealType& time );
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType& dx = 0,
+ const IndexType& dy = 0,
+ const IndexType& dz = 0 ) const;
+
+ bool setEps(const Real& eps);
+
+ private:
+
+ template< typename Vector, int AxeX = 0, int AxeY = 0, int AxeZ = 0 >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real boundaryDerivative( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType& dx = 0,
+ const IndexType& dy = 0,
+ const IndexType& dz = 0 ) const;
+
+ tnlSharedVector< RealType, DeviceType, IndexType > u;
+ tnlVector< RealType, DeviceType, IndexType> q;
+ RealType eps;
+};
+
+#include <operators/operator-Q/tnlFiniteVolumeOperatorQ_impl.h>
+
+
+#endif /* TNLFINITEVOLUMEOPERATORQ_H */
diff --git a/src/operators/operator-Q/tnlFiniteVolumeOperatorQ_impl.h b/src/operators/operator-Q/tnlFiniteVolumeOperatorQ_impl.h
new file mode 100644
index 0000000000..dbc4133329
--- /dev/null
+++ b/src/operators/operator-Q/tnlFiniteVolumeOperatorQ_impl.h
@@ -0,0 +1,636 @@
+
+#ifndef TNLFINITEVOLUMEOPERATORQ_IMPL_H
+#define TNLFINITEVOLUMEOPERATORQ_IMPL_H
+
+#include <operators/operator-Q/tnlFiniteVolumeOperatorQ.h>
+#include <mesh/tnlGrid.h>
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+tnlString
+tnlFiniteVolumeOperatorQ< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index, 0 >::
+getType()
+{
+ return tnlString( "tnlFiniteVolumeOperatorQ< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + ", 0 >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+tnlString
+tnlFiniteVolumeOperatorQ< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+getType()
+{
+ return tnlString( "tnlFiniteVolumeOperatorQ< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + ", 1 >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+bool tnlFiniteVolumeOperatorQ< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index, 0 >::setEps( const Real& eps )
+{
+ this->eps = eps;
+
+ return true;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+bool tnlFiniteVolumeOperatorQ< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index, 1 >::setEps( const Real& eps )
+{
+ this->eps = eps;
+
+ return true;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+void
+tnlFiniteVolumeOperatorQ< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+update( const MeshType& mesh, const RealType& time )
+{
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector, int AxeX, int AxeY, int AxeZ >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlFiniteVolumeOperatorQ< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index, 0 >::
+boundaryDerivative( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType& dx,
+ const IndexType& dy,
+ const IndexType& dz ) const
+{
+ return 0.0;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlFiniteVolumeOperatorQ< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index, 0 >::
+getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType& dx,
+ const IndexType& dy,
+ const IndexType& dz ) const
+{
+ return 0.0;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector, int AxeX, int AxeY, int AxeZ >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlFiniteVolumeOperatorQ< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+boundaryDerivative( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType& dx,
+ const IndexType& dy,
+ const IndexType& dz ) const
+{
+ return 0.0;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlFiniteVolumeOperatorQ< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType& dx,
+ const IndexType& dy,
+ const IndexType& dz ) const
+{
+ return 0.0;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+tnlString
+tnlFiniteVolumeOperatorQ< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index, 0 >::
+getType()
+{
+ return tnlString( "tnlFiniteVolumeOperatorQ< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + ", 0 >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+tnlString
+tnlFiniteVolumeOperatorQ< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+getType()
+{
+ return tnlString( "tnlFiniteVolumeOperatorQ< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + ", 1 >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+bool tnlFiniteVolumeOperatorQ< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index, 1 >::setEps( const Real& eps )
+{
+ this->eps = eps;
+
+ return true;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+bool tnlFiniteVolumeOperatorQ< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index, 0 >::setEps( const Real& eps )
+{
+ this->eps = eps;
+
+ return true;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector >
+Index
+tnlFiniteVolumeOperatorQ< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+bind( Vector& u)
+{
+ return 0;
+}
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+void
+tnlFiniteVolumeOperatorQ< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+update( const MeshType& mesh, const RealType& time )
+{
+ CoordinatesType dimensions = mesh.getDimensions();
+ CoordinatesType coordinates;
+
+ for( coordinates.x()=1; coordinates.x() < dimensions.x()-1; coordinates.x()++ )
+ for( coordinates.y()=1; coordinates.y() < dimensions.y()-1; coordinates.y()++ )
+ {
+ q.setElement( mesh.getCellIndex(coordinates), getValue( mesh, mesh.getCellIndex(coordinates), coordinates, u, time ) );
+ }
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector, int AxeX, int AxeY, int AxeZ >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlFiniteVolumeOperatorQ< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index, 0 >::
+boundaryDerivative( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType& dx,
+ const IndexType& dy,
+ const IndexType& dz ) const
+{
+ if ( ( AxeX == 1 ) && ( AxeY == 0 ) && ( AxeZ == 0 ) )
+ {
+ if ( ( dx == 1 ) && ( dy == 0 ) && ( dz == 0 ) )
+ return mesh.getHxInverse() * ( u[ mesh.template getCellNextToCell< 1,0 >( cellIndex ) ] - u[ cellIndex ] );
+ if ( ( dx == -1 ) && ( dy == 0 ) && ( dz == 0 ) )
+ return mesh.getHxInverse() * ( u[ cellIndex ] - u[ mesh.template getCellNextToCell< -1,0 >( cellIndex ) ] );
+ if ( ( dx == 0 ) && ( dy == 1 ) && ( dz == 0 ) )
+ return mesh.getHxInverse() * 0.25 * ( u[ mesh.template getCellNextToCell< 1,0 >( cellIndex ) ] +
+ u[ mesh.template getCellNextToCell< 1,1 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< -1,0 >( cellIndex ) ] -
+ u[ mesh.template getCellNextToCell< -1,1 >( cellIndex ) ] );
+ if ( ( dx == 0 ) && ( dy == -1 ) && ( dz == 0 ) )
+ return mesh.getHxInverse() * 0.25 * ( u[ mesh.template getCellNextToCell< 1,0 >( cellIndex ) ] +
+ u[ mesh.template getCellNextToCell< 1,-1 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< -1,0 >( cellIndex ) ] -
+ u[ mesh.template getCellNextToCell< -1,-1 >( cellIndex ) ] );
+ }
+ if ( ( AxeX == 0 ) && ( AxeY == 1 ) && ( AxeZ == 0 ) )
+ {
+ if ( ( dx == 0 ) && ( dy == 1 ) && ( dz == 0 ) )
+ return mesh.getHyInverse() * ( u[ mesh.template getCellNextToCell< 0,1 >( cellIndex ) ] - u[ cellIndex ] );
+ if ( ( dx == 0 ) && ( dy == -1 ) && ( dz == 0 ) )
+ return mesh.getHyInverse() * ( u[ cellIndex ] - u[ mesh.template getCellNextToCell< 0,-1 >( cellIndex ) ] );
+ if ( ( dx == 1 ) && ( dy == 0 ) && ( dz == 0 ) )
+ return mesh.getHyInverse() * 0.25 * ( u[ mesh.template getCellNextToCell< 0,1 >( cellIndex ) ] +
+ u[ mesh.template getCellNextToCell< 1,1 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< 0,-1 >( cellIndex ) ] -
+ u[ mesh.template getCellNextToCell< 1,-1 >( cellIndex ) ] );
+ if ( ( dx == -1 ) && ( dy == 0 ) && ( dz == 0 ) )
+ return mesh.getHyInverse() * 0.25 * ( u[ mesh.template getCellNextToCell< 0,1 >( cellIndex ) ] +
+ u[ mesh.template getCellNextToCell< -1,1 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< 0,-1 >( cellIndex ) ] -
+ u[ mesh.template getCellNextToCell< -1,-1 >( cellIndex ) ] );
+ }
+ return 0.0;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlFiniteVolumeOperatorQ< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index, 0 >::
+getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType& dx,
+ const IndexType& dy,
+ const IndexType& dz ) const
+{
+ if ( ( dx == 0 ) && ( dy == 0 ) && ( dz == 0 ) )
+ return sqrt( this->eps + ( u[ mesh.template getCellNextToCell< 0,1 >( cellIndex ) ] - u[ cellIndex ] ) *
+ ( u[ mesh.template getCellNextToCell< 0,1 >( cellIndex ) ] - u[ cellIndex ] )
+ * mesh.getHyInverse() * mesh.getHyInverse() + ( u[ mesh.template getCellNextToCell< 1,0 >( cellIndex ) ] - u[ cellIndex ] )
+ * ( u[ mesh.template getCellNextToCell< 1,0 >( cellIndex ) ] - u[ cellIndex ] ) * mesh.getHxInverse() * mesh.getHxInverse() );
+ if ( ( dx == 1 ) && ( dy == 0 ) && ( dz == 0 ) )
+ return sqrt( this->eps + this->template boundaryDerivative< Vector,1,0 >( mesh, cellIndex, coordinates, u, time, 1, 0 ) *
+ this->template boundaryDerivative< Vector,1,0 >( mesh, cellIndex, coordinates, u, time, 1, 0 ) +
+ this->template boundaryDerivative< Vector,0,1 >( mesh, cellIndex, coordinates, u, time, 1, 0 ) *
+ this->template boundaryDerivative< Vector,0,1 >( mesh, cellIndex, coordinates, u, time, 1, 0 ) );
+ if ( ( dx == -1 ) && ( dy == 0 ) && ( dz == 0 ) )
+ return sqrt( this->eps + this->template boundaryDerivative< Vector,1,0 >( mesh, cellIndex, coordinates, u, time, -1, 0 ) *
+ this->template boundaryDerivative< Vector,1,0 >( mesh, cellIndex, coordinates, u, time, -1, 0 ) +
+ this->template boundaryDerivative< Vector,0,1 >( mesh, cellIndex, coordinates, u, time, -1, 0 ) *
+ this->template boundaryDerivative< Vector,0,1 >( mesh, cellIndex, coordinates, u, time, -1, 0 ) );
+ if ( ( dx == 0 ) && ( dy == 1 ) && ( dz == 0 ) )
+ return sqrt( this->eps + this->template boundaryDerivative< Vector,1,0 >( mesh, cellIndex, coordinates, u, time, 0, 1 ) *
+ this->template boundaryDerivative< Vector,1,0 >( mesh, cellIndex, coordinates, u, time, 0, 1 ) +
+ this->template boundaryDerivative< Vector,0,1 >( mesh, cellIndex, coordinates, u, time, 0, 1 ) *
+ this->template boundaryDerivative< Vector,0,1 >( mesh, cellIndex, coordinates, u, time, 0, 1 ) );
+ if ( ( dx == 0 ) && ( dy == -1 ) && ( dz == 0 ) )
+ return sqrt( this->eps + this->template boundaryDerivative< Vector,1,0 >( mesh, cellIndex, coordinates, u, time, 0, -1 ) *
+ this->template boundaryDerivative< Vector,1,0 >( mesh, cellIndex, coordinates, u, time, 0, -1 ) +
+ this->template boundaryDerivative< Vector,0,1 >( mesh, cellIndex, coordinates, u, time, 0, -1 ) *
+ this->template boundaryDerivative< Vector,0,1 >( mesh, cellIndex, coordinates, u, time, 0, -1 ) );
+ return 0.0;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlFiniteVolumeOperatorQ< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType& dx,
+ const IndexType& dy,
+ const IndexType& dz ) const
+{
+ return q.getElement(cellIndex);
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+tnlString
+tnlFiniteVolumeOperatorQ< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index, 0 >::
+getType()
+{
+ return tnlString( "tnlFiniteVolumeOperatorQ< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + ", 0 >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+tnlString
+tnlFiniteVolumeOperatorQ< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+getType()
+{
+ return tnlString( "tnlFiniteVolumeOperatorQ< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + ", 1 >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+bool tnlFiniteVolumeOperatorQ< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index, 1 >::setEps( const Real& eps )
+{
+ this->eps = eps;
+
+ return true;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+bool tnlFiniteVolumeOperatorQ< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index, 0 >::setEps( const Real& eps )
+{
+ this->eps = eps;
+
+ return true;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector >
+Index
+tnlFiniteVolumeOperatorQ< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+bind( Vector& u)
+{
+ this->u.bind(u);
+ if(q.setSize(u.getSize()))
+ return 1;
+ q.setValue(0);
+ return 0;
+}
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+void
+tnlFiniteVolumeOperatorQ< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index, 1 >::
+update( const MeshType& mesh, const RealType& time )
+{
+ CoordinatesType dimensions = mesh.getDimensions();
+ CoordinatesType coordinates;
+
+ for( coordinates.x()=1; coordinates.x() < dimensions.x()-1; coordinates.x()++ )
+ for( coordinates.y()=1; coordinates.y() < dimensions.y()-1; coordinates.y()++ )
+ for( coordinates.z()=1; coordinates.z() < dimensions.z()-1; coordinates.z()++ )
+ q.setElement( mesh.getCellIndex(coordinates), getValue( mesh, mesh.getCellIndex(coordinates), coordinates, u, time ) );
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector, int AxeX, int AxeY, int AxeZ >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlFiniteVolumeOperatorQ< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index, 0 >::
+boundaryDerivative( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType& dx,
+ const IndexType& dy,
+ const IndexType& dz ) const
+{
+ if ( ( AxeX == 1 ) && ( AxeY == 0 ) && ( AxeZ == 0 ) )
+ {
+ if ( ( dx == 1 ) && ( dy == 0 ) && ( dz == 0 ) )
+ return mesh.getHxInverse() * ( u[ mesh.template getCellNextToCell< 1,0,0 >( cellIndex ) ] - u[ cellIndex ] );
+ if ( ( dx == -1 ) && ( dy == 0 ) && ( dz == 0 ) )
+ return mesh.getHxInverse() * ( u[ cellIndex ] - u[ mesh.template getCellNextToCell< -1,0,0 >( cellIndex ) ] );
+ if ( ( dx == 0 ) && ( dy == 1 ) && ( dz == 0 ) )
+ return mesh.getHxInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 1,0,0 >( cellIndex ) ] +
+ u[ mesh.template getCellNextToCell< 1,1,0 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< -1,0,0 >( cellIndex ) ] -
+ u[ mesh.template getCellNextToCell< -1,1,0 >( cellIndex ) ] );
+ if ( ( dx == 0 ) && ( dy == -1 ) && ( dz == 0 ) )
+ return mesh.getHxInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 1,0,0 >( cellIndex ) ] +
+ u[ mesh.template getCellNextToCell< 1,-1,0 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< -1,0,0 >( cellIndex ) ] -
+ u[ mesh.template getCellNextToCell< -1,-1,0 >( cellIndex ) ] );
+ if ( ( dx == 0 ) && ( dy == 0 ) && ( dz == 1 ) )
+ return mesh.getHxInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 1,0,0 >( cellIndex ) ] +
+ u[ mesh.template getCellNextToCell< 1,0,1 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< -1,0,0 >( cellIndex ) ] -
+ u[ mesh.template getCellNextToCell< -1,0,1 >( cellIndex ) ] );
+ if ( ( dx == 0 ) && ( dy == 0 ) && ( dz == -1 ) )
+ return mesh.getHxInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 1,0,0 >( cellIndex ) ] +
+ u[ mesh.template getCellNextToCell< 1,0,-1 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< -1,0,0 >( cellIndex ) ] -
+ u[ mesh.template getCellNextToCell< -1,0,-1 >( cellIndex ) ] );
+ }
+ if ( ( AxeX == 0 ) && ( AxeY == 1 ) && ( AxeZ == 0 ) )
+ {
+ if ( ( dx == 0 ) && ( dy == 1 ) && ( dz == 0 ) )
+ return mesh.getHyInverse() * ( u[ mesh.template getCellNextToCell< 0,1,0 >( cellIndex ) ] - u[ cellIndex ] );
+ if ( ( dx == 0 ) && ( dy == -1 ) && ( dz == 0 ) )
+ return mesh.getHyInverse() * ( u[ cellIndex ] - u[ mesh.template getCellNextToCell< 0,-1,0 >( cellIndex ) ] );
+ if ( ( dx == 1 ) && ( dy == 0 ) && ( dz == 0 ) )
+ return mesh.getHyInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 0,1,0 >( cellIndex ) ] +
+ u[ mesh.template getCellNextToCell< 1,1,0 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< 0,-1,0 >( cellIndex ) ] -
+ u[ mesh.template getCellNextToCell< 1,-1,0 >( cellIndex ) ] );
+ if ( ( dx == -1 ) && ( dy == 0 ) && ( dz == 0 ) )
+ return mesh.getHyInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 0,1,0 >( cellIndex ) ] +
+ u[ mesh.template getCellNextToCell< -1,1,0 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< 0,-1,0 >( cellIndex ) ] -
+ u[ mesh.template getCellNextToCell< -1,-1,0 >( cellIndex ) ] );
+ if ( ( dx == 0 ) && ( dy == 0 ) && ( dz == 1 ) )
+ return mesh.getHyInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 0,1,0 >( cellIndex ) ] +
+ u[ mesh.template getCellNextToCell< 0,1,1 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< 0,-1,0 >( cellIndex ) ] -
+ u[ mesh.template getCellNextToCell< 0,-1,1 >( cellIndex ) ] );
+ if ( ( dx == 0 ) && ( dy == 0 ) && ( dz == -1 ) )
+ return mesh.getHyInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 0,1,0 >( cellIndex ) ] +
+ u[ mesh.template getCellNextToCell< 0,1,-1 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< 0,-1,0 >( cellIndex ) ] -
+ u[ mesh.template getCellNextToCell< 0,-1,-1 >( cellIndex ) ] );
+ }
+ if ( ( AxeX == 0 ) && ( AxeY == 0 ) && ( AxeZ == 1 ) )
+ {
+ if ( ( dx == 0 ) && ( dy == 0 ) && ( dz == 1 ) )
+ return mesh.getHzInverse() * ( u[ mesh.template getCellNextToCell< 0,0,1 >( cellIndex ) ] - u[ cellIndex ] );
+ if ( ( dx == 0 ) && ( dy == 0 ) && ( dz == -1 ) )
+ return mesh.getHzInverse() * ( u[ cellIndex ] - u[ mesh.template getCellNextToCell< 0,0,-1 >( cellIndex ) ] );
+ if ( ( dx == 1 ) && ( dy == 0 ) && ( dz == 0 ) )
+ return mesh.getHzInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 0,0,1 >( cellIndex ) ] +
+ u[ mesh.template getCellNextToCell< 1,0,1 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< 0,0,-1 >( cellIndex ) ] -
+ u[ mesh.template getCellNextToCell< 1,0,-1 >( cellIndex ) ] );
+ if ( ( dx == -1 ) && ( dy == 0 ) && ( dz == 0 ) )
+ return mesh.getHzInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 0,0,1 >( cellIndex ) ] +
+ u[ mesh.template getCellNextToCell< -1,0,1 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< 0,0,-1 >( cellIndex ) ] -
+ u[ mesh.template getCellNextToCell< -1,0,-1 >( cellIndex ) ] );
+ if ( ( dx == 0 ) && ( dy == 1 ) && ( dz == 0 ) )
+ return mesh.getHzInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 0,0,1 >( cellIndex ) ] +
+ u[ mesh.template getCellNextToCell< 0,1,1 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< 0,0,-1 >( cellIndex ) ] -
+ u[ mesh.template getCellNextToCell< 0,1,-1 >( cellIndex ) ] );
+ if ( ( dx == 0 ) && ( dy == -1 ) && ( dz == 0 ) )
+ return mesh.getHzInverse() * 0.125 * ( u[ mesh.template getCellNextToCell< 0,0,1 >( cellIndex ) ] +
+ u[ mesh.template getCellNextToCell< 0,-1,1 >( cellIndex ) ] - u[ mesh.template getCellNextToCell< 0,0,-1 >( cellIndex ) ] -
+ u[ mesh.template getCellNextToCell< 0,-1,-1 >( cellIndex ) ] );
+ }
+ return 0.0;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlFiniteVolumeOperatorQ< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index, 0 >::
+getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType& dx,
+ const IndexType& dy,
+ const IndexType& dz ) const
+{
+ if ( ( dx == 0 ) && ( dy == 0 ) && ( dz == 0 ) )
+ return sqrt( this->eps + ( u[ mesh.template getCellNextToCell< 0,1,0 >( cellIndex ) ] - u[ cellIndex ] ) *
+ ( u[ mesh.template getCellNextToCell< 0,1,0 >( cellIndex ) ] - u[ cellIndex ] )
+ * mesh.getHyInverse() * mesh.getHyInverse() + ( u[ mesh.template getCellNextToCell< 1,0,0 >( cellIndex ) ] - u[ cellIndex ] )
+ * ( u[ mesh.template getCellNextToCell< 1,0,0 >( cellIndex ) ] - u[ cellIndex ] ) * mesh.getHxInverse() * mesh.getHxInverse()
+ + ( u[ mesh.template getCellNextToCell< 0,0,1 >( cellIndex ) ] - u[ cellIndex ] )
+ * ( u[ mesh.template getCellNextToCell< 0,0,1 >( cellIndex ) ] - u[ cellIndex ] ) * mesh.getHzInverse() * mesh.getHzInverse() );
+ if ( ( dx == 1 ) && ( dy == 0 ) && ( dz == 0 ) )
+ return sqrt( this->eps + this->template boundaryDerivative< Vector,1,0,0 >( mesh, cellIndex, coordinates, u, time, 1, 0, 0 ) *
+ this->template boundaryDerivative< Vector,1,0,0 >( mesh, cellIndex, coordinates, u, time, 1, 0, 0 ) +
+ this->template boundaryDerivative< Vector,0,1,0 >( mesh, cellIndex, coordinates, u, time, 1, 0, 0 ) *
+ this->template boundaryDerivative< Vector,0,1,0 >( mesh, cellIndex, coordinates, u, time, 1, 0, 0 ) +
+ this->template boundaryDerivative< Vector,0,0,1 >( mesh, cellIndex, coordinates, u, time, 1, 0, 0 ) *
+ this->template boundaryDerivative< Vector,0,0,1 >( mesh, cellIndex, coordinates, u, time, 1, 0, 0 ) );
+ if ( ( dx == -1 ) && ( dy == 0 ) && ( dz == 0 ) )
+ return sqrt( this->eps + this->template boundaryDerivative< Vector,1,0,0 >( mesh, cellIndex, coordinates, u, time, -1, 0, 0 ) *
+ this->template boundaryDerivative< Vector,1,0,0 >( mesh, cellIndex, coordinates, u, time, -1, 0, 0 ) +
+ this->template boundaryDerivative< Vector,0,1,0 >( mesh, cellIndex, coordinates, u, time, -1, 0, 0 ) *
+ this->template boundaryDerivative< Vector,0,1,0 >( mesh, cellIndex, coordinates, u, time, -1, 0, 0 ) +
+ this->template boundaryDerivative< Vector,0,0,1 >( mesh, cellIndex, coordinates, u, time, -1, 0, 0 ) *
+ this->template boundaryDerivative< Vector,0,0,1 >( mesh, cellIndex, coordinates, u, time, -1, 0, 0 ) );
+ if ( ( dx == 0 ) && ( dy == 1 ) && ( dz == 0 ) )
+ return sqrt( this->eps + this->template boundaryDerivative< Vector,1,0,0 >( mesh, cellIndex, coordinates, u, time, 0, 1, 0 ) *
+ this->template boundaryDerivative< Vector,1,0,0 >( mesh, cellIndex, coordinates, u, time, 0, 1, 0 ) +
+ this->template boundaryDerivative< Vector,0,1,0 >( mesh, cellIndex, coordinates, u, time, 0, 1, 0 ) *
+ this->template boundaryDerivative< Vector,0,1,0 >( mesh, cellIndex, coordinates, u, time, 0, 1, 0 ) +
+ this->template boundaryDerivative< Vector,0,0,1 >( mesh, cellIndex, coordinates, u, time, 0, 1, 0 ) *
+ this->template boundaryDerivative< Vector,0,0,1 >( mesh, cellIndex, coordinates, u, time, 0, 1, 0 ));
+ if ( ( dx == 0 ) && ( dy == -1 ) && ( dz == 0 ) )
+ return sqrt( this->eps + this->template boundaryDerivative< Vector,1,0,0 >( mesh, cellIndex, coordinates, u, time, 0, -1, 0 ) *
+ this->template boundaryDerivative< Vector,1,0,0 >( mesh, cellIndex, coordinates, u, time, 0, -1, 0 ) +
+ this->template boundaryDerivative< Vector,0,1,0 >( mesh, cellIndex, coordinates, u, time, 0, -1, 0 ) *
+ this->template boundaryDerivative< Vector,0,1,0 >( mesh, cellIndex, coordinates, u, time, 0, -1, 0 ) +
+ this->template boundaryDerivative< Vector,0,0,1 >( mesh, cellIndex, coordinates, u, time, 0, -1, 0 ) *
+ this->template boundaryDerivative< Vector,0,0,1 >( mesh, cellIndex, coordinates, u, time, 0, -1, 0 ) );
+ if ( ( dx == 0 ) && ( dy == 0 ) && ( dz == 1 ) )
+ return sqrt( this->eps + this->template boundaryDerivative< Vector,1,0,0 >( mesh, cellIndex, coordinates, u, time, 0, 0, 1 ) *
+ this->template boundaryDerivative< Vector,1,0,0 >( mesh, cellIndex, coordinates, u, time, 0, 0, 1 ) +
+ this->template boundaryDerivative< Vector,0,1,0 >( mesh, cellIndex, coordinates, u, time, 0, 0, 1 ) *
+ this->template boundaryDerivative< Vector,0,1,0 >( mesh, cellIndex, coordinates, u, time, 0, 0, 1 ) +
+ this->template boundaryDerivative< Vector,0,0,1 >( mesh, cellIndex, coordinates, u, time, 0, 0, 1 ) *
+ this->template boundaryDerivative< Vector,0,0,1 >( mesh, cellIndex, coordinates, u, time, 0, 0, 1 ));
+ if ( ( dx == 0 ) && ( dy == 0 ) && ( dz == -1 ) )
+ return sqrt( this->eps + this->template boundaryDerivative< Vector,1,0,0 >( mesh, cellIndex, coordinates, u, time, 0, 0, -1 ) *
+ this->template boundaryDerivative< Vector,1,0,0 >( mesh, cellIndex, coordinates, u, time, 0, 0, -1 ) +
+ this->template boundaryDerivative< Vector,0,1,0 >( mesh, cellIndex, coordinates, u, time, 0, 0, -1 ) *
+ this->template boundaryDerivative< Vector,0,1,0 >( mesh, cellIndex, coordinates, u, time, 0, 0, -1 ) +
+ this->template boundaryDerivative< Vector,0,0,1 >( mesh, cellIndex, coordinates, u, time, 0, 0, -1 ) *
+ this->template boundaryDerivative< Vector,0,0,1 >( mesh, cellIndex, coordinates, u, time, 0, 0, -1 ) );
+ return 0.0;
+}
+#endif /* TNLFINITEVOLUMEOPERATORQ_IMPL_H */
--
GitLab