From 3aa7debdeba8a71b7933364059273cde314ee9c5 Mon Sep 17 00:00:00 2001
From: Tomas Oberhuber <tomas.oberhuber@fjfi.cvut.cz>
Date: Wed, 8 May 2013 22:32:37 +0200
Subject: [PATCH] Implementing linear diffusion on the deformed grids.
---
src/implementation/mesh/tnlGrid2D_impl.h | 14 ++-
.../mesh/tnlIdenticalGridGeometry_impl.h | 26 +++-
.../diffusion/tnlLinearDiffusion_impl.h | 118 +++++++++++++++++-
src/mesh/tnlGrid.h | 3 +
src/mesh/tnlIdenticalGridGeometry.h | 4 +-
src/schemes/diffusion/tnlLinearDiffusion.h | 2 +
6 files changed, 162 insertions(+), 5 deletions(-)
diff --git a/src/implementation/mesh/tnlGrid2D_impl.h b/src/implementation/mesh/tnlGrid2D_impl.h
index 13a549b57a..54f32cc510 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 a0c0bc2c78..20cc9c8645 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 959045105e..4c4e5abaec 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 8611578c72..e0b7d2404c 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 c3b6df7bb7..13012cfdae 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 7c767e06bb..f4589984c6 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();
--
GitLab