Implementing linear diffusion on the deformed grids.

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,

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,

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

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;

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;
 

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();