From b8f31f4b2efa025519862974ccb2bbc8c2816642 Mon Sep 17 00:00:00 2001
From: Tomas Oberhuber <tomas.oberhuber@fjfi.cvut.cz>
Date: Tue, 2 Jun 2015 20:19:43 +0200
Subject: [PATCH] Fixing use of MatrixRow in matrxi assembler.

---
 .../diffusion/tnlLinearDiffusion_impl.h       | 15 ++++++++------
 ...AnalyticDirichletBoundaryConditions_impl.h |  5 +++--
 ...nlAnalyticNeumannBoundaryConditions_impl.h | 15 ++++++++------
 .../tnlDirichletBoundaryConditions_impl.h     |  5 +++--
 .../tnlNeumannBoundaryConditions_impl.h       | 15 ++++++++------
 src/solvers/pde/tnlLinearSystemAssembler.h    | 20 +++++++------------
 6 files changed, 40 insertions(+), 35 deletions(-)

diff --git a/src/operators/diffusion/tnlLinearDiffusion_impl.h b/src/operators/diffusion/tnlLinearDiffusion_impl.h
index 07abc05d24..c7fcc8a0fd 100644
--- a/src/operators/diffusion/tnlLinearDiffusion_impl.h
+++ b/src/operators/diffusion/tnlLinearDiffusion_impl.h
@@ -76,7 +76,7 @@ template< typename MeshReal,
           typename MeshIndex,
           typename Real,
           typename Index >
-   template< typename Vector, typename MatrixRow >
+   template< typename Vector, typename Matrix >
 __cuda_callable__
 void
 tnlLinearDiffusion< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
@@ -87,8 +87,9 @@ updateLinearSystem( const RealType& time,
                     const CoordinatesType& coordinates,
                     Vector& u,
                     Vector& b,
-                    MatrixRow& matrixRow ) const
+                    Matrix& matrix ) const
 {
+   typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
    const RealType lambdaX = tau * mesh.getHxSquareInverse();
    //printf( "tau = %f lambda = %f dx_sqr = %f dx = %f, \n", tau, lambdaX, mesh.getHxSquareInverse(), mesh.getHx() );
    matrixRow.setElement( 0, mesh.template getCellNextToCell< -1 >( index ),     - lambdaX );
@@ -156,7 +157,7 @@ template< typename MeshReal,
           typename MeshIndex,
           typename Real,
           typename Index >
-   template< typename Vector, typename MatrixRow >
+   template< typename Vector, typename Matrix >
 __cuda_callable__
 void
 tnlLinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
@@ -167,8 +168,9 @@ updateLinearSystem( const RealType& time,
                     const CoordinatesType& coordinates,
                     Vector& u,
                     Vector& b,
-                    MatrixRow& matrixRow ) const
+                    Matrix& matrix ) const
 {
+   typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
    const RealType lambdaX = tau * mesh.getHxSquareInverse();
    const RealType lambdaY = tau * mesh.getHySquareInverse();
    matrixRow.setElement( 0, mesh.template getCellNextToCell< 0, -1 >( index ), -lambdaY );
@@ -240,7 +242,7 @@ template< typename MeshReal,
           typename MeshIndex,
           typename Real,
           typename Index >
-   template< typename Vector, typename MatrixRow >
+   template< typename Vector, typename Matrix >
 __cuda_callable__
 void
 tnlLinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
@@ -251,8 +253,9 @@ updateLinearSystem( const RealType& time,
                     const CoordinatesType& coordinates,
                     Vector& u,
                     Vector& b,
-                    MatrixRow& matrixRow ) const
+                    Matrix& matrix ) const
 {
+   typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
    const RealType lambdaX = tau * mesh.getHxSquareInverse();
    const RealType lambdaY = tau * mesh.getHySquareInverse();
    const RealType lambdaZ = tau * mesh.getHzSquareInverse();
diff --git a/src/operators/tnlAnalyticDirichletBoundaryConditions_impl.h b/src/operators/tnlAnalyticDirichletBoundaryConditions_impl.h
index f744de07d6..da90d69aa2 100644
--- a/src/operators/tnlAnalyticDirichletBoundaryConditions_impl.h
+++ b/src/operators/tnlAnalyticDirichletBoundaryConditions_impl.h
@@ -135,7 +135,7 @@ template< int Dimensions,
           typename Function,
           typename Real,
           typename Index >
-   template< typename MatrixRow >          
+   template< typename Matrix >          
 __cuda_callable__
 void
 tnlAnalyticDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Function, Real, Index >::
@@ -145,8 +145,9 @@ updateLinearSystem( const RealType& time,
                     const CoordinatesType& coordinates,
                     DofVectorType& u,
                     DofVectorType& b,
-                    MatrixRow& matrixRow ) const
+                    Matrix& matrix ) const
 {
+   typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
    matrixRow.setElement( 0, index, 1.0 );
    b[ index ] = function.getValue( mesh.template getCellCenter< VertexType >( coordinates ), time );
 }
diff --git a/src/operators/tnlAnalyticNeumannBoundaryConditions_impl.h b/src/operators/tnlAnalyticNeumannBoundaryConditions_impl.h
index 42f15cf38f..5e65dde3f5 100644
--- a/src/operators/tnlAnalyticNeumannBoundaryConditions_impl.h
+++ b/src/operators/tnlAnalyticNeumannBoundaryConditions_impl.h
@@ -112,7 +112,7 @@ template< typename MeshReal,
           typename Function,
           typename Real,
           typename Index >
-   template< typename MatrixRow >
+   template< typename Matrix >
 __cuda_callable__
 void
 tnlAnalyticNeumannBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >::
@@ -122,8 +122,9 @@ updateLinearSystem( const RealType& time,
                     const CoordinatesType& coordinates,
                     DofVectorType& u,
                     DofVectorType& b,
-                    MatrixRow& matrixRow ) const
+                    Matrix& matrix ) const
 {
+   typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
    const Real functionValue = this->function.getValue( mesh.template getCellCenter< VertexType >( coordinates ), time );
    if( coordinates.x() == 0 )
    {
@@ -205,7 +206,7 @@ template< typename MeshReal,
           typename Function,
           typename Real,
           typename Index >
-   template< typename MatrixRow >
+   template< typename Matrix >
 __cuda_callable__
 void
 tnlAnalyticNeumannBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >::
@@ -215,8 +216,9 @@ updateLinearSystem( const RealType& time,
                     const CoordinatesType& coordinates,
                     DofVectorType& u,
                     DofVectorType& b,
-                    MatrixRow& matrixRow ) const
+                    Matrix& matrix ) const
 {
+   typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
    const Real functionValue = this->function.getValue( mesh.template getCellCenter< VertexType >( coordinates ), time );
    if( coordinates.x() == 0 )
    {
@@ -320,7 +322,7 @@ template< typename MeshReal,
           typename Function,
           typename Real,
           typename Index >
-   template< typename MatrixRow >
+   template< typename Matrix >
 __cuda_callable__
 void
 tnlAnalyticNeumannBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >::
@@ -330,8 +332,9 @@ updateLinearSystem( const RealType& time,
                     const CoordinatesType& coordinates,
                     DofVectorType& u,
                     DofVectorType& b,
-                    MatrixRow& matrixRow ) const
+                    Matrix& matrix ) const
 {
+   typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
    const Real functionValue = this->function.getValue( mesh.template getCellCenter< VertexType >( coordinates ), time );
    if( coordinates.x() == 0 )
    {
diff --git a/src/operators/tnlDirichletBoundaryConditions_impl.h b/src/operators/tnlDirichletBoundaryConditions_impl.h
index 1fc9b245a8..c392491807 100644
--- a/src/operators/tnlDirichletBoundaryConditions_impl.h
+++ b/src/operators/tnlDirichletBoundaryConditions_impl.h
@@ -128,7 +128,7 @@ template< int Dimensions,
           typename Vector,
           typename Real,
           typename Index >
-   template< typename MatrixRow >
+   template< typename Matrix >
 __cuda_callable__
 void
 tnlDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Vector, Real, Index >::
@@ -138,8 +138,9 @@ updateLinearSystem( const RealType& time,
                     const CoordinatesType& coordinates,
                     DofVectorType& u,
                     DofVectorType& b,
-                    MatrixRow& matrixRow ) const
+                    Matrix& matrix ) const
 {
+   typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
    matrixRow.setElement( 0, index, 1.0 );
    b[ index ] = this->vector[ index ];
 }
diff --git a/src/operators/tnlNeumannBoundaryConditions_impl.h b/src/operators/tnlNeumannBoundaryConditions_impl.h
index 333bf09da0..a31150150d 100644
--- a/src/operators/tnlNeumannBoundaryConditions_impl.h
+++ b/src/operators/tnlNeumannBoundaryConditions_impl.h
@@ -89,7 +89,7 @@ template< typename MeshReal,
           typename Vector,
           typename Real,
           typename Index >
-   template< typename MatrixRow >
+   template< typename Matrix >
 __cuda_callable__
 void
 tnlNeumannBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >, Vector, Real, Index >::
@@ -99,8 +99,9 @@ updateLinearSystem( const RealType& time,
                     const CoordinatesType& coordinates,
                     DofVectorType& u,
                     DofVectorType& b,
-                    MatrixRow& matrixRow ) const
+                    Matrix& matrix ) const
 {
+   typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
    if( coordinates.x() == 0 )
    {
       matrixRow.setElement( 0, index, 1.0 );
@@ -179,7 +180,7 @@ template< typename MeshReal,
           typename Vector,
           typename Real,
           typename Index >
-   template< typename MatrixRow >
+   template< typename Matrix >
 __cuda_callable__
 void
 tnlNeumannBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >, Vector, Real, Index >::
@@ -189,8 +190,9 @@ updateLinearSystem( const RealType& time,
                     const CoordinatesType& coordinates,
                     DofVectorType& u,
                     DofVectorType& b,
-                    MatrixRow& matrixRow ) const
+                    Matrix& matrix ) const
 {
+   typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
    if( coordinates.x() == 0 )
    {
       matrixRow.setElement( 0, index,                            1.0 );
@@ -291,7 +293,7 @@ template< typename MeshReal,
           typename Vector,
           typename Real,
           typename Index >
-   template< typename MatrixRow >
+   template< typename Matrix >
 __cuda_callable__
 void
 tnlNeumannBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >, Vector, Real, Index >::
@@ -301,8 +303,9 @@ updateLinearSystem( const RealType& time,
                     const CoordinatesType& coordinates,
                     DofVectorType& u,
                     DofVectorType& b,
-                    MatrixRow& matrixRow ) const
+                    Matrix& matrix ) const
 {
+   typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
    if( coordinates.x() == 0 )
    {
       matrixRow.setElement( 0, index,                            1.0 );
diff --git a/src/solvers/pde/tnlLinearSystemAssembler.h b/src/solvers/pde/tnlLinearSystemAssembler.h
index 4ac6c99d3e..065f2cf3d9 100644
--- a/src/solvers/pde/tnlLinearSystemAssembler.h
+++ b/src/solvers/pde/tnlLinearSystemAssembler.h
@@ -118,13 +118,12 @@ class tnlLinearSystemAssembler
                                     TraverserUserData& userData,
                                     const IndexType index )
          {
-            typename MatrixType::MatrixRow matrixRow = userData.matrix->getRow( index );
             userData.boundaryConditions->updateLinearSystem( *userData.time + *userData.tau,
                                                              mesh,
                                                              index,
                                                              *userData.u,
                                                              *userData.b,
-                                                             matrixRow );
+                                                             *userData.matrix );
          }
 
    };
@@ -148,14 +147,13 @@ class tnlLinearSystemAssembler
                                                                     index,
                                                                     *userData.time );*/
 
-            typename MatrixType::MatrixRow matrixRow = userData.matrix->getRow( index );
             userData.differentialOperator->updateLinearSystem( *userData.time,
                                                                *userData.tau,
                                                                mesh,
                                                                index,
                                                                *userData.u,
                                                                *userData.b,
-                                                               matrixRow );
+                                                               *userData.matrix );
             //userData.matrix->addElement( index, index, 1.0, 1.0 );
             const RealType& rhs = FunctionAdapter::getValue( mesh,
                                                              *userData.rightHandSide,
@@ -239,15 +237,14 @@ class tnlLinearSystemAssembler< tnlGrid< Dimensions, Real, Device, Index >,
                                   const CoordinatesType& coordinates )
          {
             //printf( "index = %d \n", index );
-             ( *userData.b )[ index ] = 0.0;
-            typename MatrixType::MatrixRow matrixRow = userData.matrix->getRow( index );
+             ( *userData.b )[ index ] = 0.0;           
             userData.boundaryConditions->updateLinearSystem( *userData.time + *userData.tau,
                                                              mesh,
                                                              index,
                                                              coordinates,
                                                              *userData.u,
                                                              *userData.b,
-                                                             matrixRow );
+                                                             *userData.matrix );
          }
 
 #ifdef HAVE_CUDA
@@ -261,14 +258,13 @@ class tnlLinearSystemAssembler< tnlGrid< Dimensions, Real, Device, Index >,
             //printf( "index = %d \n", index );
             // printf("Matrix assembler: Index = %d \n", index );
             ( *userData.b )[ index ] = 0.0;
-            typename MatrixType::MatrixRow matrixRow = userData.matrix->getRow( index );
             userData.boundaryConditions->updateLinearSystem( *userData.time,
                                                              mesh,
                                                              index,
                                                              coordinates,
                                                              *userData.u,
                                                              *userData.b,
-                                                             matrixRow );
+                                                             *userData.matrix );
             //printf( "BC: index = %d, b = %f \n", index, ( *userData.b )[ index ] );
          }
 
@@ -301,7 +297,6 @@ class tnlLinearSystemAssembler< tnlGrid< Dimensions, Real, Device, Index >,
                                                              coordinates,
                                                              *userData.time );*/
             
-            typename MatrixType::MatrixRow matrixRow = userData.matrix->getRow( index );
             userData.differentialOperator->updateLinearSystem( *userData.time,
                                                                *userData.tau,
                                                                mesh,
@@ -309,7 +304,7 @@ class tnlLinearSystemAssembler< tnlGrid< Dimensions, Real, Device, Index >,
                                                                coordinates,
                                                                *userData.u,
                                                                *userData.b,
-                                                               matrixRow );
+                                                               *userData.matrix );
             /*if( *userData.timeDiscretisationCoefficient != 0.0 )
                userData.matrix->addElementFast( index,
                                                 index,
@@ -348,7 +343,6 @@ class tnlLinearSystemAssembler< tnlGrid< Dimensions, Real, Device, Index >,
                                                              coordinates,
                                                              *userData.time );*/
 
-            typename MatrixType::MatrixRow matrixRow = userData.matrix->getRow( index );
             userData.differentialOperator->updateLinearSystem( *userData.time,
                                                                *userData.tau,
                                                                mesh,
@@ -356,7 +350,7 @@ class tnlLinearSystemAssembler< tnlGrid< Dimensions, Real, Device, Index >,
                                                                coordinates,
                                                                *userData.u,
                                                                *userData.b,
-                                                               matrixRow );
+                                                               *userData.matrix );
             /*if( *userData.timeDiscretisationCoefficient != 0.0 )
                userData.matrix->addElementFast( index,
                                                 index,
-- 
GitLab