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

Implementing linear grid geometry.

parent ebcf1213

src/core/mfuncs.h

+1 −0
Original line number Diff line number Diff line
@@ -77,4 +77,5 @@ void swap( T& a, T& b) 
   b = aux;

}*/





#endif

src/core/tnlTuple.h

+44 −0
Original line number Diff line number Diff line
@@ -637,4 +637,48 @@ template< int Size, typename Real > tnlString getParameterType() 
          tnlString( " >" );

};



template< typename Real >

tnlTuple< 3, Real > tnlVectorProduct( const tnlTuple< 3, Real >& u,

                                      const tnlTuple< 3, Real >& v )

{

   tnlTuple< 3, Real > p;

   p[ 0 ] = u[ 1 ] * v[ 2 ] - u[ 2 ] * v[ 1 ];

   p[ 1 ] = u[ 2 ] * v[ 0 ] - u[ 0 ] * v[ 2 ];

   p[ 2 ] = u[ 0 ] * v[ 1 ] - u[ 1 ] * v[ 0 ];

   return p;

};



template< typename Real >

Real tnlScalarProduct( const tnlTuple< 2, Real >& u,

                       const tnlTuple< 2, Real >& v )

{

   return u[ 0 ] * v[ 0 ] + u[ 1 ] * v[ 1 ];

};



template< typename Real >

Real tnlScalarProduct( const tnlTuple< 3, Real >& u,

                       const tnlTuple< 3, Real >& v )

{

   return u[ 0 ] * v[ 0 ] + u[ 1 ] * v[ 1 ] + u[ 2 ] * v[ 2 ];

};



template< typename Real >

Real tnlTriangleArea( const tnlTuple< 2, Real >& a,

                      const tnlTuple< 2, Real >& b,

                      const tnlTuple< 2, Real >& c )

{

   tnlTuple< 3, Real > u1, u2;

   u1. x() = b. x() - a. x();

   u1. y() = b. y() - a. y();

   u1. z() = 0.0;

   u2. x() = c. x() - a. x();

   u2. y() = c. y() - a. y();

   u2. z() = 0;



   const tnlTuple< 3, Real > v = tnlVectorProduct( u1, u2 );

   return 0.5 * sqrt( tnlScalarProduct( v, v ) );

};







#endif

src/implementation/mesh/tnlLinearGridGeometry_impl.h

+20 −11
Original line number Diff line number Diff line
@@ -96,7 +96,7 @@ template< Index dx > 
void tnlLinearGridGeometry< 1, Real, Device, Index > :: getEdgeNormal( const Index i,

                                                                          VertexType& normal ) const

{

   tnlAssert( dx == 1 || dx == -1, cerr << " dx = " << dx << " dy = " << dy << endl );

   tnlAssert( dx == 1 || dx == -1, cerr << " dx = " << dx << endl );

   normal. x() = dx;

}
@@ -159,7 +159,16 @@ template< typename Real, 
          typename Index >

Real tnlLinearGridGeometry< 2, Real, Device, Index > :: getElementMeasure( const CoordinatesType& coordinates ) const

{

   return elementMeasure;

   VertexType v0, v1, v2, v3;

   const VertexType origin( 0.0 );

   this -> template getVertex< -1, -1 >( coordinates, origin, v0 );

   this -> template getVertex<  1, -1  >( coordinates, origin, v1 );

   this -> template getVertex<  1,  1 >( coordinates, origin, v2 );

   this -> template getVertex< -1,  1 >( coordinates, origin, v3 );



   return tnlTriangleArea( v0, v1, v3 ) + tnlTriangleArea( v2, v3, v1 );



   //return elementMeasure;

}



template< typename Real,
@@ -168,6 +177,10 @@ template< typename Real, 
   template< int dx, int dy >

Real tnlLinearGridGeometry< 2, Real, Device, Index > :: getElementCoVolumeMeasure( const CoordinatesType& coordinates ) const

{

   tnlAssert( ( dx == 0 && ( dy == 1 || dy == -1 ) ) ||

              ( dy == 0 && ( dx == 1 || dx == -1 ) ),

              cerr << " dx = " << dx << " dy = " << dy << endl );



   return 0.5 * elementMeasure;

}
@@ -222,10 +235,10 @@ template< Index dx, Index dy > 
void tnlLinearGridGeometry< 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 && dx + dy != 0 , cerr << " dx = " << dx << " dy = " << dy << endl );

   /*VertexType v1, v2, origin( 0 );

   tnlAssert( ( dx == 0 && ( dy == 1 || dy == -1 ) ) ||

              ( dy == 0 && ( dx == 1 || dx == -1 ) ),

              cerr << " dx = " << dx << " dy = " << dy << endl );

   VertexType v1, v2, origin( 0.0 );

   if( dy == 0 )

   {

      if( dx == 1 )
@@ -253,11 +266,7 @@ void tnlLinearGridGeometry< 2, Real, Device, Index > :: getEdgeNormal( const Coo 
      }

   }

   normal. x() = v1. y() - v2. y();

   normal. y() = v2. x() - v1. x();*/



   normal. x() = dx * parametricStep. y();

   normal. y() = dy * parametricStep. x();



   normal. y() = v2. x() - v1. x();

}



template< typename Real,

src/implementation/schemes/gradient/tnlCentralFDMGradient_impl.h

+27 −8
Original line number Diff line number Diff line
@@ -51,17 +51,17 @@ template< typename Real, 
          typename Device,

          typename Index,

          template< int, typename, typename, typename > class GridGeometry >

void tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index, GridGeometry > > :: getGradient( const Index& i,

void tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index, GridGeometry > > :: getGradient( const Index& c,

                                                                                              VertexType& grad_f  ) const

{

   tnlAssert( this -> mesh, cerr << "No mesh was set in tnlCentralFDMGradient. Use the bindMesh method." );



   const Index e = mesh -> getElementNeighbour( i,  1,  0 );

   const Index w = mesh -> getElementNeighbour( i, -1,  0 );

   const Index n = mesh -> getElementNeighbour( i,  0,  1 );

   const Index s = mesh -> getElementNeighbour( i,  0, -1 );

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

   CoordinatesType cCoordinates;

   mesh -> getElementCoordinates( i, cCoordinates );

   mesh -> getElementCoordinates( c, cCoordinates );

   CoordinatesType eCoordinates( cCoordinates ),

                   wCoordinates( cCoordinates ),

                   nCoordinates( cCoordinates ),
@@ -71,8 +71,27 @@ void tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index, GridGeometry > > :: 
   nCoordinates. y() ++;

   sCoordinates. y() --;



   grad_f. x() = ( f[ e ] - f[ w ] ) / ( mesh -> getElementsDistance( eCoordinates, wCoordinates ) );

   grad_f. y() = ( f[ n ] - f[ s ] ) / ( mesh -> getElementsDistance( nCoordinates, sCoordinates ) );



   VertexType eNormal, nNormal, wNormal, sNormal;

   mesh -> template getEdgeNormal<  1,  0 >( cCoordinates, eNormal );

   mesh -> template getEdgeNormal<  0,  1 >( cCoordinates, nNormal );

   mesh -> template getEdgeNormal< -1,  0 >( cCoordinates, wNormal );

   mesh -> template getEdgeNormal<  0, -1 >( cCoordinates, sNormal );

   const RealType cMu = mesh -> getElementMeasure( cCoordinates );

   const RealType eMu = mesh -> getElementMeasure( eCoordinates );

   const RealType nMu = mesh -> getElementMeasure( nCoordinates );

   const RealType wMu = mesh -> getElementMeasure( wCoordinates );

   const RealType sMu = mesh -> getElementMeasure( sCoordinates );

   const RealType f_e = ( cMu * f[ c ] + eMu * f[ e ] ) / ( cMu + eMu );

   const RealType f_n = ( cMu * f[ c ] + nMu * f[ n ] ) / ( cMu + nMu );

   const RealType f_w = ( cMu * f[ c ] + wMu * f[ w ] ) / ( cMu + wMu );

   const RealType f_s = ( cMu * f[ c ] + sMu * f[ s ] ) / ( cMu + sMu );

   grad_f. x() = 1.0 / cMu * ( f_e * eNormal. x() + f_n * nNormal. x() + f_w * wNormal. x() + f_s * sNormal. x() );

   grad_f. y() = 1.0 / cMu * ( f_e * eNormal. y() + f_n * nNormal. y() + f_w * wNormal. y() + f_s * sNormal. y() );





   //grad_f. x() = ( f[ e ] - f[ w ] ) / ( mesh -> getElementsDistance( eCoordinates, wCoordinates ) );

   //grad_f. y() = ( f[ n ] - f[ s ] ) / ( mesh -> getElementsDistance( nCoordinates, sCoordinates ) );

}



/****

src/matrix/tnlCSRMatrix.h

+1 −1
Original line number Diff line number Diff line
@@ -1064,7 +1064,7 @@ bool tnlCSRMatrix< Real, Device, Index > :: read( istream& file, 
      }

      Index I = atoi( parsed_line[ 0 ]. getString() ) - 1;

      Index J = atoi( parsed_line[ 1 ]. getString() ) - 1;

      double A = atof( parsed_line[ 2 ]. getString() );

      Real A = ( Real ) atof( parsed_line[ 2 ]. getString() );



      if( I < 0 || I >= size || J < 0 || J >= size )

      {