diff --git a/src/implementation/mesh/tnlGrid2D_impl.h b/src/implementation/mesh/tnlGrid2D_impl.h
index 13a549b57aa91dfae0d6270d0798c318dc7ebfc4..54f32cc510e10ccfbd0744347e95890de6f0eb31 100644
--- a/src/implementation/mesh/tnlGrid2D_impl.h
+++ b/src/implementation/mesh/tnlGrid2D_impl.h
@@ -247,6 +247,16 @@ Real tnlGrid< 2, Real, Device, Index, Geometry > :: getElementMeasure( const Coo
return geometry. getElementMeasure( coordinates );
}
+template< typename Real,
+ typename Device,
+ typename Index,
+ template< int, typename, typename, typename > class Geometry >
+ template< int dx, int dy >
+Real tnlGrid< 2, Real, Device, Index, Geometry > :: getElementCoVolumeMeasure( const CoordinatesType& coordinates ) const
+{
+ return geometry. getElementCoVolumeMeasure< dx, dy >( coordinates );
+}
+
template< typename Real,
typename Device,
typename Index,
@@ -288,7 +298,9 @@ template< typename Real,
void tnlGrid< 2, Real, Device, Index, Geometry > :: getVertex( const CoordinatesType& elementCoordinates,
VertexType& vertex ) const
{
- return geometry. getVertex< dx, dy >( elementCoordinates, vertex );
+ return geometry. getVertex< dx, dy >( elementCoordinates,
+ this -> origin,
+ vertex );
}
template< typename Real,
diff --git a/src/implementation/mesh/tnlIdenticalGridGeometry_impl.h b/src/implementation/mesh/tnlIdenticalGridGeometry_impl.h
index a0c0bc2c78bf380855c6b3eb5925aa387b59f482..20cc9c86456714a180d0339ec7bc7c9440f283e6 100644
--- a/src/implementation/mesh/tnlIdenticalGridGeometry_impl.h
+++ b/src/implementation/mesh/tnlIdenticalGridGeometry_impl.h
@@ -163,6 +163,16 @@ Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementMeasure( co
return elementMeasure;
}
+template< typename Real,
+ typename Device,
+ typename Index >
+ template< int dx, int dy >
+Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementCoVolumeMeasure( const CoordinatesType& coordinates ) const
+{
+ return elementMeasure;
+}
+
+
template< typename Real,
typename Device,
typename Index >
@@ -217,7 +227,21 @@ void tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getEdgeNormal( const
dy == 0 || dy == 1 || dy == -1 ) &&
dx * dy == 0, cerr << " dx = " << dx << " dy = " << dy << endl );
normal. x() = dx * parametricStep. y();
- normal. y() = dy * parametricStep. x();;
+ normal. y() = dy * parametricStep. x();
+}
+
+template< typename Real,
+ typename Device,
+ typename Index >
+ template< Index dx, Index dy >
+void tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getVertex( const CoordinatesType& coordinates,
+ const VertexType& origin,
+ VertexType& vertex ) const
+{
+ tnlAssert( ( dx == 0 || dx == 1 || dx == -1 ||
+ dy == 0 || dy == 1 || dy == -1 ), cerr << " dx = " << dx << " dy = " << dy << endl );
+ vertex. x() = origin. x() + ( coordinates. x() + 0.5 * ( 1 + dx ) ) * parametricStep. x();
+ vertex. y() = origin. y() + ( coordinates. y() + 0.5 * ( 1 + dy ) ) * parametricStep. y();
}
template< typename Real,
diff --git a/src/implementation/schemes/diffusion/tnlLinearDiffusion_impl.h b/src/implementation/schemes/diffusion/tnlLinearDiffusion_impl.h
index 959045105e9284428f4ed19a8885234421c4130f..4c4e5abaece42d5a536f1617543a65fb159e9e56 100644
--- a/src/implementation/schemes/diffusion/tnlLinearDiffusion_impl.h
+++ b/src/implementation/schemes/diffusion/tnlLinearDiffusion_impl.h
@@ -51,8 +51,122 @@ Real tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index > > :: getDiffusion( co
const Index n = mesh -> getElementNeighbour( c, 0, 1 );
const Index s = mesh -> getElementNeighbour( c, 0, -1 );
- return ( f[ e ] - 2.0 * f[ c ] + f[ w ] ) / ( hx * hx ) +
- ( f[ n ] - 2.0 * f[ c ] + f[ s ] ) / ( hy * hy );
+ //return ( f[ e ] - 2.0 * f[ c ] + f[ w ] ) / ( hx * hx ) +
+ // ( f[ n ] - 2.0 * f[ c ] + f[ s ] ) / ( hy * hy );
+
+ /****
+ * Diffusion for the deformed grids
+ */
+ const Index en = mesh -> getElementNeighbour( c, 1, 1 );
+ const Index es = mesh -> getElementNeighbour( c, 1, -1 );
+ const Index wn = mesh -> getElementNeighbour( c, -1, 1 );
+ const Index ws = mesh -> getElementNeighbour( c, -1, -1 );
+
+ CoordinatesType cCoordinates, eCoordinates, nCoordinates, wCoordinates, sCoordinates,
+ enCoordinates, esCoordinates, wnCoordinates, wsCoordinates;
+ this -> mesh -> getElementCoordinates( c, cCoordinates );
+ eCoordinates = nCoordinates = wCoordinates = sCoordinates = cCoordinates;
+ eCoordinates. x() ++;
+ wCoordinates. x() --;
+ nCoordinates. y() ++;
+ sCoordinates. y() --;
+ enCoordinates = esCoordinates = eCoordinates;
+ enCoordinates. y() ++;
+ esCoordinates. y() --;
+ wnCoordinates = wsCoordinates = wCoordinates;
+ wnCoordinates. y() ++;
+ wsCoordinates. y() --;
+
+ const RealType cMeasure = this -> mesh -> getElementMeasure( cCoordinates );
+ const RealType eMeasure = this -> mesh -> getElementMeasure( eCoordinates );
+ const RealType wMeasure = this -> mesh -> getElementMeasure( wCoordinates );
+ const RealType nMeasure = this -> mesh -> getElementMeasure( nCoordinates );
+ const RealType sMeasure = this -> mesh -> getElementMeasure( sCoordinates );
+ const RealType enMeasure = this -> mesh -> getElementMeasure( enCoordinates );
+ const RealType esMeasure = this -> mesh -> getElementMeasure( esCoordinates );
+ const RealType wnMeasure = this -> mesh -> getElementMeasure( wnCoordinates );
+ const RealType wsMeasure = this -> mesh -> getElementMeasure( wsCoordinates );
+
+ const RealType f_en = ( cMeasure * f[ c ] + eMeasure * f[ e ] +
+ nMeasure * f[ n ] + enMeasure * f[ en ] ) /
+ ( cMeasure + eMeasure + nMeasure + enMeasure );
+ const RealType f_es = ( cMeasure * f[ c ] + eMeasure * f[ e ] +
+ sMeasure * f[ s ] + esMeasure * f[ es ] ) /
+ ( cMeasure + eMeasure + sMeasure + esMeasure );
+ const RealType f_wn = ( cMeasure * f[ c ] + wMeasure * f[ w ] +
+ nMeasure * f[ n ] + wnMeasure * f[ wn ] ) /
+ ( cMeasure + wMeasure + nMeasure + wnMeasure );
+ const RealType f_ws = ( cMeasure * f[ c ] + wMeasure * f[ w ] +
+ sMeasure * f[ s ] + wsMeasure * f[ ws ] ) /
+ ( cMeasure + wMeasure + sMeasure + wsMeasure );
+
+ const RealType f_cen = 0.5 * ( f[ c ] + f_en );
+ const RealType f_ces = 0.5 * ( f[ c ] + f_es );
+ const RealType f_cwn = 0.5 * ( f[ c ] + f_wn );
+ const RealType f_cws = 0.5 * ( f[ c ] + f_ws );
+
+ VertexType cCenter, eCenter, nCenter, wCenter, sCenter;
+ this -> mesh -> getElementCenter( cCoordinates, cCenter );
+ this -> mesh -> getElementCenter( eCoordinates, eCenter );
+ this -> mesh -> getElementCenter( wCoordinates, wCenter );
+ this -> mesh -> getElementCenter( nCoordinates, nCenter );
+ this -> mesh -> getElementCenter( sCoordinates, sCenter );
+
+ VertexType enVertex, esVertex, wnVertex, wsVertex;
+ this -> mesh -> template getVertex< 1, 1 >( cCoordinates, enVertex );
+ this -> mesh -> template getVertex< 1, -1 >( cCoordinates, esVertex );
+ this -> mesh -> template getVertex< -1, 1 >( cCoordinates, wnVertex );
+ this -> mesh -> template getVertex< -1, -1 >( cCoordinates, wsVertex );
+
+ const RealType f_x_e = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< 1, 0 >( cCoordinates ) *
+ ( f_cen * ( enVertex. y() - cCenter. y() ) +
+ f_ces * ( esVertex. y() - cCenter. y() ) +
+ 0.5 * ( f_es + f[ e ] ) * ( esVertex. y() - eCenter. y() ) +
+ 0.5 * ( f_en + f[ e ] ) * ( enVertex. y() - eCenter. y() ) );
+ const RealType f_y_e = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< 1, 0 >( cCoordinates ) *
+ ( f_cen * ( enVertex. x() - cCenter. x() ) +
+ f_ces * ( esVertex. x() - cCenter. x() ) +
+ 0.5 * ( f_es + f[ e ] ) * ( esVertex. x() - eCenter. x() ) +
+ 0.5 * ( f_en + f[ e ] ) * ( enVertex. x() - eCenter. x() ) );
+ const RealType f_x_w = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< -1, 0 >( cCoordinates ) *
+ ( f_cwn * ( wnVertex. y() - cCenter. y() ) +
+ f_cws * ( wsVertex. y() - cCenter. y() ) +
+ 0.5 * ( f_ws + f[ w ] ) * ( wsVertex. y() - wCenter. y() ) +
+ 0.5 * ( f_wn + f[ w ] ) * ( wnVertex. y() - wCenter. y() ) );
+ const RealType f_y_w = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< -1, 0 >( cCoordinates ) *
+ ( f_cwn * ( wnVertex. x() - cCenter. x() ) +
+ f_cws * ( wsVertex. x() - cCenter. x() ) +
+ 0.5 * ( f_ws + f[ w ] ) * ( wsVertex. x() - wCenter. x() ) +
+ 0.5 * ( f_wn + f[ w ] ) * ( wnVertex. x() - wCenter. x() ) );
+ const RealType f_x_n = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< 0, 1 >( cCoordinates ) *
+ ( f_cen * ( enVertex. y() - cCenter. y() ) +
+ f_cwn * ( wnVertex. y() - cCenter. y() ) +
+ 0.5 * ( f_en + f[ n ] ) * ( enVertex. y() - nCenter. y() ) +
+ 0.5 * ( f_wn + f[ n ] ) * ( wnVertex. y() - nCenter. y() ) );
+ const RealType f_y_n = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< 0, 1 >( cCoordinates ) *
+ ( f_cen * ( enVertex. x() - cCenter. x() ) +
+ f_cwn * ( wnVertex. x() - cCenter. x() ) +
+ 0.5 * ( f_en + f[ n ] ) * ( enVertex. x() - nCenter. x() ) +
+ 0.5 * ( f_wn + f[ n ] ) * ( wnVertex. x() - nCenter. x() ) );
+ const RealType f_x_s = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< 0, -1 >( cCoordinates ) *
+ ( f_ces * ( esVertex. y() - cCenter. y() ) +
+ f_cws * ( wsVertex. y() - cCenter. y() ) +
+ 0.5 * ( f_es + f[ s ] ) * ( esVertex. y() - sCenter. y() ) +
+ 0.5 * ( f_ws + f[ s ] ) * ( wsVertex. y() - sCenter. y() ) );
+ const RealType f_y_s = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< 0, -1 >( cCoordinates ) *
+ ( f_ces * ( esVertex. x() - cCenter. x() ) +
+ f_cws * ( wsVertex. x() - cCenter. x() ) +
+ 0.5 * ( f_es + f[ s ] ) * ( enVertex. x() - sCenter. x() ) +
+ 0.5 * ( f_ws + f[ s ] ) * ( wnVertex. x() - sCenter. x() ) );
+
+ VertexType eNormal, nNormal, wNormal, sNormal;
+ this -> mesh -> template getEdgeNormal< 1, 0 >( cCoordinates, eNormal );
+ this -> mesh -> template getEdgeNormal< -1, 0 >( cCoordinates, wNormal );
+ this -> mesh -> template getEdgeNormal< 0, 1 >( cCoordinates, nNormal );
+ this -> mesh -> template getEdgeNormal< 0, -1 >( cCoordinates, sNormal );
+ return 1.0 / this -> mesh -> getElementMeasure( cCoordinates ) *
+ ( f_x_e * eNormal. x() + f_x_w * wNormal. x() +
+ f_y_n * nNormal. y() + f_y_s * sNormal. y() );
}
#endif
diff --git a/src/mesh/tnlGrid.h b/src/mesh/tnlGrid.h
index 8611578c72ee7e685c67f843c9b3e23aeec00198..e0b7d2404cf13ba530ea4570fc76f37002b47c59 100644
--- a/src/mesh/tnlGrid.h
+++ b/src/mesh/tnlGrid.h
@@ -161,6 +161,9 @@ class tnlGrid< 2, Real, Device, Index, Geometry > : public tnlObject
Real getElementMeasure( const CoordinatesType& coordinates ) const;
+ template< int dx, int dy >
+ Real getElementCoVolumeMeasure( const CoordinatesType& coordinates ) const;
+
template< int dx, int dy >
void getEdgeNormal( const CoordinatesType& elementCoordinates,
VertexType& normal ) const;
diff --git a/src/mesh/tnlIdenticalGridGeometry.h b/src/mesh/tnlIdenticalGridGeometry.h
index c3b6df7bb7562e577dc6cf753e7d3a3d49731485..13012cfdaec71caa4c2a4f09983730b8546d89ad 100644
--- a/src/mesh/tnlIdenticalGridGeometry.h
+++ b/src/mesh/tnlIdenticalGridGeometry.h
@@ -49,7 +49,6 @@ class tnlIdenticalGridGeometry< 1, Real, Device, Index >
const CoordinatesType& coordinates,
VertexType& center ) const;
-
Real getElementMeasure( const CoordinatesType& i ) const;
template< Index dx >
@@ -105,6 +104,9 @@ class tnlIdenticalGridGeometry< 2, Real, Device, Index >
Real getElementMeasure( const tnlTuple< 2, Index >& coordinates ) const;
+ template< int dx, int dy >
+ Real getElementCoVolumeMeasure( const CoordinatesType& coordinates ) const;
+
Real getElementsDistance( const CoordinatesType& c1,
const CoordinatesType& c2 ) const;
diff --git a/src/schemes/diffusion/tnlLinearDiffusion.h b/src/schemes/diffusion/tnlLinearDiffusion.h
index 7c767e06bbf18fb7638e4331ca96807a7edd99ee..f4589984c67f8c2a89572ebabaaffe5016092bda 100644
--- a/src/schemes/diffusion/tnlLinearDiffusion.h
+++ b/src/schemes/diffusion/tnlLinearDiffusion.h
@@ -34,6 +34,8 @@ class tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index > >
typedef Real RealType;
typedef Device DeviceType;
typedef Index IndexType;
+ typedef typename tnlGrid< 2, Real, Device, Index > :: CoordinatesType CoordinatesType;
+ typedef typename tnlGrid< 2, Real, Device, Index > :: VertexType VertexType;
tnlLinearDiffusion();