Commit 6200b55b authored by Tomáš Oberhuber's avatar Tomáš Oberhuber
Browse files

Implementing the tnlGrid and the Lax-Fridrichs scheme.

parent e42fb244
Loading
Loading
Loading
Loading
+13 −9
Original line number Diff line number Diff line
@@ -66,7 +66,7 @@ void tnlGrid< 1, Real, Device, Index, Geometry > :: setDimensions( const Index x
   this -> dimensions. x() = xSize;
   dofs = xSize;

   tnlTuple< 1, Real > parametricStep;
   VertexType parametricStep;
   parametricStep. x() = proportions. x() / xSize;
   geometry. setParametricStep( parametricStep );
}
@@ -75,7 +75,7 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry  >
void tnlGrid< 1, Real, Device, Index, Geometry > :: setDimensions( const tnlTuple< 1, Index >& dimensions )
void tnlGrid< 1, Real, Device, Index, Geometry > :: setDimensions( const CoordinatesType& dimensions )
{
   this -> setDimensions( dimensions. x() );
}
@@ -84,7 +84,8 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry  >
const tnlTuple< 1, Index >& tnlGrid< 1, Real, Device, Index, Geometry > :: getDimensions() const
const typename tnlGrid< 1, Real, Device, Index, Geometry > :: CoordinatesType& 
   tnlGrid< 1, Real, Device, Index, Geometry > :: getDimensions() const
{
   return this -> dimensions;
}
@@ -93,7 +94,7 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry  >
void tnlGrid< 1, Real, Device, Index, Geometry > :: setOrigin( const tnlTuple< 1, Real >& origin )
void tnlGrid< 1, Real, Device, Index, Geometry > :: setOrigin( const VertexType& origin )
{
   this -> origin = origin;
}
@@ -102,7 +103,8 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry  >
const tnlTuple< 1, Real >& tnlGrid< 1, Real, Device, Index, Geometry > :: getOrigin() const
const typename tnlGrid< 1, Real, Device, Index, Geometry > :: VertexType& 
  tnlGrid< 1, Real, Device, Index, Geometry > :: getOrigin() const
{
   return this -> origin;
}
@@ -111,7 +113,7 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry  >
void tnlGrid< 1, Real, Device, Index, Geometry > :: setProportions( const tnlTuple< 1, Real >& proportions )
void tnlGrid< 1, Real, Device, Index, Geometry > :: setProportions( const VertexType& proportions )
{
   this -> proportions = proportions;
   this -> setDimensions( this -> dimensions );
@@ -121,7 +123,8 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry  >
const tnlTuple< 1, Real >& tnlGrid< 1, Real, Device, Index, Geometry > :: getProportions() const
const typename tnlGrid< 1, Real, Device, Index, Geometry > :: VertexType& 
   tnlGrid< 1, Real, Device, Index, Geometry > :: getProportions() const
{
   return this -> proportions;
}
@@ -130,7 +133,7 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry  >
void tnlGrid< 1, Real, Device, Index, Geometry > :: setParametricStep( const tnlTuple< 1, Real >& parametricStep )
void tnlGrid< 1, Real, Device, Index, Geometry > :: setParametricStep( const VertexType& parametricStep )
{
   this -> proportions. x() = this -> dimensions. x() * parametricStep. x();
   geometry. setParametricStep( parametricStep );
@@ -140,7 +143,8 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry  >
const tnlTuple< 1, Real >& tnlGrid< 1, Real, Device, Index, Geometry > :: getParametricStep() const
const typename tnlGrid< 1, Real, Device, Index, Geometry > :: VertexType& 
   tnlGrid< 1, Real, Device, Index, Geometry > :: getParametricStep() const
{
   return geometry. getParametricStep();
}
+37 −22
Original line number Diff line number Diff line
@@ -68,7 +68,7 @@ void tnlGrid< 2, Real, Device, Index, Geometry > :: setDimensions( const Index x
   this -> dimensions. x() = xSize;
   this -> dimensions. y() = ySize;
   dofs = ySize * xSize;
   tnlTuple< 2, Real > parametricStep;
   VertexType parametricStep;
   parametricStep. x() = proportions. x() / xSize;
   parametricStep. y() = proportions. y() / ySize;
   geometry. setParametricStep( parametricStep );
@@ -78,7 +78,7 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry >
void tnlGrid< 2, Real, Device, Index, Geometry > :: setDimensions( const tnlTuple< 2, Index >& dimensions )
void tnlGrid< 2, Real, Device, Index, Geometry > :: setDimensions( const CoordinatesType& dimensions )
{
   this -> setDimensions( dimensions. x(), dimensions. y() );
}
@@ -87,7 +87,8 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry >
const tnlTuple< 2, Index >& tnlGrid< 2, Real, Device, Index, Geometry > :: getDimensions() const
const typename tnlGrid< 2, Real, Device, Index, Geometry > :: CoordinatesType& 
   tnlGrid< 2, Real, Device, Index, Geometry > :: getDimensions() const
{
   return this -> dimensions;
}
@@ -96,7 +97,7 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry >
void tnlGrid< 2, Real, Device, Index, Geometry > :: setOrigin( const tnlTuple< 2, Real >& origin )
void tnlGrid< 2, Real, Device, Index, Geometry > :: setOrigin( const VertexType& origin )
{
   this -> origin = origin;
}
@@ -105,7 +106,8 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry >
const tnlTuple< 2, Real >& tnlGrid< 2, Real, Device, Index, Geometry > :: getOrigin() const
const typename tnlGrid< 2, Real, Device, Index, Geometry > :: VertexType& 
   tnlGrid< 2, Real, Device, Index, Geometry > :: getOrigin() const
{
   return this -> origin;
}
@@ -114,10 +116,10 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry >
void tnlGrid< 2, Real, Device, Index, Geometry > :: setProportions( const tnlTuple< 2, Real >& proportions )
void tnlGrid< 2, Real, Device, Index, Geometry > :: setProportions( const VertexType& proportions )
{
   this -> proportions = proportions;
   tnlTuple< 2, Real > parametricStep;
   VertexType parametricStep;
   parametricStep. x() = proportions. x() / ( this -> dimensions. x() - 1 );
   parametricStep. y() = proportions. y() / ( this -> dimensions. y() - 1 );
   geometry. setParametricStep( parametricStep );
@@ -127,7 +129,8 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry >
const tnlTuple< 2, Real >& tnlGrid< 2, Real, Device, Index, Geometry > :: getProportions() const
const typename tnlGrid< 2, Real, Device, Index, Geometry > :: VertexType& 
   tnlGrid< 2, Real, Device, Index, Geometry > :: getProportions() const
{
   return this -> proportions;
}
@@ -136,7 +139,7 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry >
void tnlGrid< 2, Real, Device, Index, Geometry > :: setParametricStep( const tnlTuple< 2, Real >& spaceStep )
void tnlGrid< 2, Real, Device, Index, Geometry > :: setParametricStep( const VertexType& spaceStep )
{
   this -> proportions. x() = this -> dimensions. x() *
                              geometry. getParametricStep(). x();
@@ -148,7 +151,8 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry >
const tnlTuple< 2, Real >& tnlGrid< 2, Real, Device, Index, Geometry > :: getParametricStep() const
const typename tnlGrid< 2, Real, Device, Index, Geometry > :: VertexType& 
   tnlGrid< 2, Real, Device, Index, Geometry > :: getParametricStep() const
{
   /*tnlAssert( dimensions. x() > 0,
              cerr << "Cannot get the space step hx since number of Elements along the x axis is not known in tnlGrid "
@@ -188,7 +192,7 @@ template< typename Real,
          typename Index,
          template< int, typename, typename, typename > class Geometry >
void tnlGrid< 2, Real, Device, Index, Geometry > :: getElementCoordinates( const Index element,
                                                                           tnlTuple< 2, Index >& coordinates ) const
                                                                           CoordinatesType& coordinates ) const
{
   tnlAssert( i >= 0 && i < dofs,
              cerr << " i = " << i << " dofs = " << dofs
@@ -203,8 +207,8 @@ template< typename Real,
          typename Index,
          template< int, typename, typename, typename > class Geometry >
Index tnlGrid< 2, Real, Device, Index, Geometry > :: getElementNeighbour( const Index Element,
                                                            const Index dy,
                                                            const Index dx ) const
                                                                          const Index dx,
                                                                          const Index dy ) const
{
   tnlAssert( Element + dy * this -> dimensions. x() + dx < getDofs(),
              cerr << "Index of neighbour with dx = " << dx
@@ -228,8 +232,8 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry >
void tnlGrid< 2, Real, Device, Index, Geometry > :: getElementCenter( const tnlTuple< 2, Index >& coordinates,
                                                                      tnlTuple< 2, Real >& center ) const
void tnlGrid< 2, Real, Device, Index, Geometry > :: getElementCenter( const CoordinatesType& coordinates,
                                                                      VertexType& center ) const
{
      geometry. getElementCenter( origin, coordinates, center );
}
@@ -238,7 +242,7 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry >
Real tnlGrid< 2, Real, Device, Index, Geometry > :: getElementMeasure( const tnlTuple< 2, Index >& coordinates ) const
Real tnlGrid< 2, Real, Device, Index, Geometry > :: getElementMeasure( const CoordinatesType& coordinates ) const
{
   return geometry. getElementMeasure( coordinates );
}
@@ -263,18 +267,29 @@ Real tnlGrid< 2, Real, Device, Index, Geometry > :: getEdgeLength( const Index j
                                                                   const Index i ) const
{
   return geometry. getEdgeLength< dy, dx >( j, i );
}*/

template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry >
template< int dx, int dy >
void tnlGrid< 2, Real, Device, Index, Geometry > :: getEdgeNormal( const CoordinatesType& coordinates,
                                                                   VertexType& normal ) const
{
   return geometry. getEdgeNormal< dx, dy >( coordinates, normal );
}

template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry >
template< int dy, int dx >
tnlTuple< 2, Real > tnlGrid< 2, Real, Device, Index, Geometry > :: getEdgeNormal( const Index j,
                                                                                  const Index i ) const
   template< int dx, int dy >
void tnlGrid< 2, Real, Device, Index, Geometry > :: getVertex( const CoordinatesType& elementCoordinates,
                                                               VertexType& vertex ) const
{
   return geometry. getEdgeNormal< dy, dx >( j, i );
}*/
   return geometry. getVertex< dx, dy >( elementCoordinates, vertex );
}

template< typename Real,
          typename Device,
@@ -317,7 +332,7 @@ bool tnlGrid< 2, Real, Device, Index, Geometry > :: load( tnlFile& file )
      return false;
   this -> dofs = this -> getDimensions(). x() *
                  this -> getDimensions(). y();
   tnlTuple< 2, Real > parametricStep;
   VertexType parametricStep;
   parametricStep. x() = proportions. x() / ( this -> dimensions. x() - 1 );
   parametricStep. y() = proportions. y() / ( this -> dimensions. y() - 1 );
   geometry. setParametricStep( parametricStep );
+13 −27
Original line number Diff line number Diff line
@@ -74,7 +74,7 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry >
void tnlGrid< 3, Real, Device, Index, Geometry > :: setDimensions( const tnlTuple< 3, Index >& dimensions )
void tnlGrid< 3, Real, Device, Index, Geometry > :: setDimensions( const CoordinatesType& dimensions )
{
   this -> setDimensions( this -> dimensions. x(), this -> dimensions. y(), this -> dimensions. z() );
}
@@ -83,7 +83,8 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry >
const tnlTuple< 3, Index >& tnlGrid< 3, Real, Device, Index, Geometry > :: getDimensions() const
const typename tnlGrid< 3, Real, Device, Index, Geometry > :: CoordinatesType& 
   tnlGrid< 3, Real, Device, Index, Geometry > :: getDimensions() const
{
   return this -> dimensions;
}
@@ -92,7 +93,7 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry >
void tnlGrid< 3, Real, Device, Index, Geometry > :: setOrigin( const tnlTuple< 3, Real >& origin )
void tnlGrid< 3, Real, Device, Index, Geometry > :: setOrigin( const VertexType& origin )
{
   this -> origin = origin;
}
@@ -101,7 +102,8 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry >
const tnlTuple< 3, Real >& tnlGrid< 3, Real, Device, Index, Geometry > :: getOrigin() const
const typename tnlGrid< 3, Real, Device, Index, Geometry > :: VertexType&
   tnlGrid< 3, Real, Device, Index, Geometry > :: getOrigin() const
{
   return this -> origin;
}
@@ -110,7 +112,7 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry >
void tnlGrid< 3, Real, Device, Index, Geometry > :: setProportions( const tnlTuple< 3, Real >& proportions )
void tnlGrid< 3, Real, Device, Index, Geometry > :: setProportions( const VertexType& proportions )
{
   this -> proportions = proportions;
}
@@ -119,7 +121,8 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry >
const tnlTuple< 3, Real >& tnlGrid< 3, Real, Device, Index, Geometry > :: getProportions() const
const typename tnlGrid< 3, Real, Device, Index, Geometry > :: VertexType&
   tnlGrid< 3, Real, Device, Index, Geometry > :: getProportions() const
{
   return this -> proportions;
}
@@ -128,7 +131,7 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry >
void tnlGrid< 3, Real, Device, Index, Geometry > :: setParametricStep( const tnlTuple< 3, Real >& spaceStep )
void tnlGrid< 3, Real, Device, Index, Geometry > :: setParametricStep( const VertexType& spaceStep )
{
   this -> proportions. x() = this -> dimensions. x() *
                              spaceStep. x();
@@ -142,27 +145,10 @@ template< typename Real,
          typename Device,
          typename Index,
          template< int, typename, typename, typename > class Geometry >
tnlTuple< 3, Real > tnlGrid< 3, Real, Device, Index, Geometry > :: getParametricStep() const
const typename tnlGrid< 3, Real, Device, Index, Geometry > :: VertexType&
   tnlGrid< 3, Real, Device, Index, Geometry > :: getParametricStep() const
{
   tnlAssert( dimensions. x() > 0,
              cerr << "Cannot get the space step hx since number of Elements along the x axis is not known in tnlGrid "
                   << this -> getName() );
   tnlAssert( dimensions. y() > 0,
              cerr << "Cannot get the space step hy since number of Elements along the y axis is not known in tnlGrid "
                   << this -> getName() );
   tnlAssert( dimensions. z() > 0,
              cerr << "Cannot get the space step hz since number of Elements along the z axis is not known in tnlGrid "
                   << this -> getName() );

   tnlTuple< 3, RealType > spaceStep;
   spaceStep. x() =
            this -> proportions. x() / ( Real ) ( this -> dimensions. x() - 1 );
   spaceStep. y() =
            this -> proportions. y() / ( Real ) ( this -> dimensions. y() - 1 );
   spaceStep. z() =
            this -> proportions. z() / ( Real ) ( this -> dimensions. z() - 1 );

   return spaceStep;
   //return geometry. getParametricStep();
}

template< typename Real,
+25 −24
Original line number Diff line number Diff line
@@ -28,7 +28,7 @@
template< typename Real,
          typename Device,
          typename Index >
void tnlIdenticalGridGeometry< 1, Real, Device, Index > :: setParametricStep( const tnlTuple< 1, Real >& parametricStep )
void tnlIdenticalGridGeometry< 1, Real, Device, Index > :: setParametricStep( const VertexType& parametricStep )
{
   this -> parametricStep = parametricStep;
   this -> elementMeasure = this -> parametricStep. x();
@@ -37,7 +37,8 @@ void tnlIdenticalGridGeometry< 1, Real, Device, Index > :: setParametricStep( co
template< typename Real,
          typename Device,
          typename Index >
const tnlTuple< 1, Real >& tnlIdenticalGridGeometry< 1, Real, Device, Index > :: getParametricStep() const
const typename tnlIdenticalGridGeometry< 1, Real, Device, Index > :: VertexType& 
   tnlIdenticalGridGeometry< 1, Real, Device, Index > :: getParametricStep() const
{
   return this -> parametricStep;
}
@@ -45,9 +46,9 @@ const tnlTuple< 1, Real >& tnlIdenticalGridGeometry< 1, Real, Device, Index > ::
template< typename Real,
          typename Device,
          typename Index >
void tnlIdenticalGridGeometry< 1, Real, Device, Index > :: getElementCenter( const tnlTuple< 1, Real >& origin,
                                                                             const tnlTuple< 1, Index >& coordinates,
                                                                             tnlTuple< 1, Real >& center ) const
void tnlIdenticalGridGeometry< 1, Real, Device, Index > :: getElementCenter( const VertexType& origin,
                                                                             const CoordinatesType& coordinates,
                                                                             VertexType& center ) const
{
   center. x() = ( coordinates. x() + 0.5 ) * parametricStep. x();
}
@@ -55,7 +56,7 @@ void tnlIdenticalGridGeometry< 1, Real, Device, Index > :: getElementCenter( con
template< typename Real,
          typename Device,
          typename Index >
Real tnlIdenticalGridGeometry< 1, Real, Device, Index > :: getElementMeasure( const Index i ) const
Real tnlIdenticalGridGeometry< 1, Real, Device, Index > :: getElementMeasure( const CoordinatesType& i ) const
{
   return elementMeasure;
}
@@ -74,8 +75,8 @@ template< typename Real,
          typename Index >
template< Index dx >
void tnlIdenticalGridGeometry< 1, Real, Device, Index > :: getEdgeCoordinates( const Index i,
                                                                               const tnlTuple< 1, Real >& origin,
                                                                               tnlTuple< 1, Real >& coordinates ) const
                                                                               const VertexType& origin,
                                                                               VertexType& coordinates ) const
{
   coordinates. x() = origin. x() + ( i + 0.5 * ( 1.0 + dx ) ) * parametricStep. x();
}
@@ -94,7 +95,7 @@ template< typename Real,
          typename Index >
template< Index dx >
void tnlIdenticalGridGeometry< 1, Real, Device, Index > :: getEdgeNormal( const Index i,
                                                                          tnlTuple< 1, Real >& normal ) const
                                                                          VertexType& normal ) const
{
   tnlAssert( dx == 1 || dx == -1, cerr << " dx = " << dx << " dy = " << dy << endl );
   normal. x() = dx;
@@ -128,7 +129,7 @@ bool tnlIdenticalGridGeometry< 1, Real, Device, Index > :: load( tnlFile& file )
template< typename Real,
          typename Device,
          typename Index >
void tnlIdenticalGridGeometry< 2, Real, Device, Index > :: setParametricStep( const tnlTuple< 2, Real >& parametricStep )
void tnlIdenticalGridGeometry< 2, Real, Device, Index > :: setParametricStep( const VertexType& parametricStep )
{
   this -> parametricStep = parametricStep;
   this -> elementMeasure = this -> parametricStep. x() * this -> parametricStep. y();
@@ -137,7 +138,8 @@ void tnlIdenticalGridGeometry< 2, Real, Device, Index > :: setParametricStep( co
template< typename Real,
          typename Device,
          typename Index >
const tnlTuple< 2, Real >& tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getParametricStep() const
const typename tnlIdenticalGridGeometry< 2, Real, Device, Index > :: VertexType&
   tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getParametricStep() const
{
   return this -> parametricStep;
}
@@ -145,9 +147,9 @@ const tnlTuple< 2, Real >& tnlIdenticalGridGeometry< 2, Real, Device, Index > ::
template< typename Real,
          typename Device,
          typename Index >
void tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementCenter( const tnlTuple< 2, Real >& origin,
                                                                             const tnlTuple< 2, Index >& coordinates,
                                                                             tnlTuple< 2, Real >& center ) const
void tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementCenter( const VertexType& origin,
                                                                             const CoordinatesType& coordinates,
                                                                             VertexType& center ) const
{
   center. x() = ( coordinates. x() + 0.5 ) * parametricStep. x();
   center. y() = ( coordinates. y() + 0.5 ) * parametricStep. y();
@@ -156,12 +158,12 @@ void tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementCenter( con
template< typename Real,
          typename Device,
          typename Index >
Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementMeasure( const tnlTuple< 2, Index >& coordinates ) const
Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementMeasure( const CoordinatesType& coordinates ) const
{
   return elementMeasure;
}

template< typename Real,
/*template< typename Real,
          typename Device,
          typename Index >
   template< Index dx, Index dy >
@@ -183,8 +185,8 @@ template< typename Real,
template< Index dx, Index dy >
void tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getEdgeCoordinates( const Index i,
                                                                               const Index j,
                                                                               const tnlTuple< 2, Real >& origin,
                                                                               tnlTuple< 2, Real >& coordinates ) const
                                                                               const VertexType& origin,
                                                                               VertexType& coordinates ) const
{
   coordinates. x() = origin. x() + ( i + 0.5 * ( 1.0 + dx ) ) * parametricStep. x();
   coordinates. y() = origin. y() + ( j + 0.5 * ( 1.0 + dy ) ) * parametricStep. y();
@@ -202,21 +204,20 @@ Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getEdgeLength( const
   if( dy == 1 && dx == 0 )
      return parametricStep. x();
   tnlAssert( false, cerr << "Bad values of dx and dy - dx = " << dx << " dy = " << dy );
}
}*/

template< typename Real,
          typename Device,
          typename Index >
template< Index dx, Index dy >
void tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getEdgeNormal( const Index i,
                                                                          const Index j,
                                                                          tnlTuple< 2, Real >& normal ) const
void tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getEdgeNormal( const CoordinatesType& coordinates,                                                                     
                                                                          VertexType& normal ) const
{
   tnlAssert( ( dx == 0 || dx == 1 || dx == -1 ||
                dy == 0 || dy == 1 || dy == -1 ) &&
               dx * dy == 0, cerr << " dx = " << dx << " dy = " << dy << endl );
   normal. x() = dx;
   normal. y() = dy;
   normal. x() = dx * parametricStep. y();
   normal. y() = dy * parametricStep. x();;
}

template< typename Real,
+4 −4
Original line number Diff line number Diff line
@@ -46,10 +46,10 @@ Real tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index > > :: getDiffusion( co
   const Real hx = mesh -> getParametricStep(). x();
   const Real hy = mesh -> getParametricStep(). y();

   const Index e = mesh -> getElementNeighbour( c,  0,  1 );
   const Index w = mesh -> getElementNeighbour( c,  0, -1 );
   const Index n = mesh -> getElementNeighbour( c,  1,  0 );
   const Index s = mesh -> getElementNeighbour( c, -1,  0 );
   const Index e = mesh -> getElementNeighbour( c,  1,  0 );
   const Index w = mesh -> getElementNeighbour( c, -1,  0 );
   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 );
Loading