To co jsme tady udelali

examples/incompressible-navier-stokes/tnlLinearDiffusion.h

+#ifndef TNLLINEARDIFFUSION_H
+#define	TNLLINEARDIFFUSION_H
+
+#include <core/vectors/tnlVector.h>
+#include <mesh/tnlGrid.h>
+
+template< typename Mesh,
+          typename Real = typename Mesh::RealType,
+          typename Index = typename Mesh::IndexType >
+class tnlLinearDiffusion
+{
+ 
+};
+
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+class tnlLinearDiffusion< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index >
+{
+   public: 
+   
+   typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType;
+   typedef typename MeshType::CoordinatesType CoordinatesType;
+   typedef Real RealType;
+   typedef Device DeviceType;
+   typedef Index IndexType;
+
+   static tnlString getType();
+   
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   Real getValue( const MeshType& mesh,
+                  const IndexType cellIndex,
+                  const CoordinatesType& coordinates,
+                  const Vector& u,
+                  const RealType& time ) const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   Index getLinearSystemRowLength( const MeshType& mesh,
+                                   const IndexType& index,
+                                   const CoordinatesType& coordinates ) const;
+
+   template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+      void updateLinearSystem( const RealType& time,
+                               const RealType& tau,
+                               const MeshType& mesh,
+                               const IndexType& index,
+                               const CoordinatesType& coordinates,
+                               Vector& u,
+                               Vector& b,
+                               MatrixRow& matrixRow ) const;
+
+};
+
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+class tnlLinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >
+{
+   public: 
+   
+   typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType;
+   typedef typename MeshType::CoordinatesType CoordinatesType;
+   typedef Real RealType;
+   typedef Device DeviceType;
+   typedef Index IndexType;
+
+   static tnlString getType();
+   
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   Real getValue( const MeshType& mesh,
+                  const IndexType cellIndex,
+                  const CoordinatesType& coordinates,
+                  const Vector& u,
+                  const Real& time ) const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   Index getLinearSystemRowLength( const MeshType& mesh,
+                                   const IndexType& index,
+                                   const CoordinatesType& coordinates ) const;
+
+   template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+      void updateLinearSystem( const RealType& time,
+                               const RealType& tau,
+                               const MeshType& mesh,
+                               const IndexType& index,
+                               const CoordinatesType& coordinates,
+                               Vector& u,
+                               Vector& b,
+                               MatrixRow& matrixRow ) const;
+};
+
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+class tnlLinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >
+{
+   public: 
+   
+   typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType;
+   typedef typename MeshType::CoordinatesType CoordinatesType;
+   typedef Real RealType;
+   typedef Device DeviceType;
+   typedef Index IndexType;
+
+   static tnlString getType();
+   
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   Real getValue( const MeshType& mesh,
+                  const IndexType cellIndex,
+                  const CoordinatesType& coordinates,
+                  const Vector& u,
+                  const Real& time ) const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   Index getLinearSystemRowLength( const MeshType& mesh,
+                                   const IndexType& index,
+                                   const CoordinatesType& coordinates ) const;
+
+   template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+      void updateLinearSystem( const RealType& time,
+                               const RealType& tau,
+                               const MeshType& mesh,
+                               const IndexType& index,
+                               const CoordinatesType& coordinates,
+                               Vector& u,
+                               Vector& b,
+                               MatrixRow& matrixRow ) const;
+
+};
+
+
+#include <operators/diffusion/tnlLinearDiffusion_impl.h>
+
+
+#endif	/* TNLLINEARDIFFUSION_H */

examples/incompressible-navier-stokes/tnlLinearDiffusion_impl.h

+
+#ifndef TNLLINEARDIFFUSION_IMP_H
+#define	TNLLINEARDIFFUSION_IMP_H
+
+#include <operators/diffusion/tnlLinearDiffusion.h>
+#include <mesh/tnlGrid.h>
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+tnlString
+tnlLinearDiffusion< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+   return tnlString( "tnlLinearDiffusion< " ) +
+          MeshType::getType() + ", " +
+          ::getType< Real >() + ", " +
+          ::getType< Index >() + " >";
+}
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlLinearDiffusion< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
+getValue( const MeshType& mesh,
+          const IndexType cellIndex,
+          const CoordinatesType& coordinates,
+          const Vector& u,
+          const Real& time ) const
+{
+   return ( u[ mesh.getCellXPredecessor( cellIndex ) ]
+            - 2.0 * u[ cellIndex ]
+            + u[ mesh.getCellXSuccessor( cellIndex ) ] ) * mesh.getHxSquareInverse();
+}
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+Index
+tnlLinearDiffusion< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
+getLinearSystemRowLength( const MeshType& mesh,
+                          const IndexType& index,
+                          const CoordinatesType& coordinates ) const
+{
+   return 3;
+}
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+   template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+void
+tnlLinearDiffusion< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
+updateLinearSystem( const RealType& time,
+                    const RealType& tau,
+                    const MeshType& mesh,
+                    const IndexType& index,
+                    const CoordinatesType& coordinates,
+                    Vector& u,
+                    Vector& b,
+                    MatrixRow& matrixRow ) const
+{
+   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.getCellXPredecessor( index ),     - lambdaX );
+   matrixRow.setElement( 1, index,                             2.0 * lambdaX );
+   matrixRow.setElement( 2, mesh.getCellXSuccessor( index ),       - lambdaX );
+   //printf( "Linear diffusion index %d columns %d %d %d \n", index, columns[ 0 ], columns[ 1 ], columns[ 2 ] );
+}
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+tnlString
+tnlLinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+   return tnlString( "tnlLinearDiffusion< " ) +
+          MeshType::getType() + ", " +
+          ::getType< Real >() + ", " +
+          ::getType< Index >() + " >";
+}
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+Index
+tnlLinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+getLinearSystemRowLength( const MeshType& mesh,
+                          const IndexType& index,
+                          const CoordinatesType& coordinates ) const
+{
+   return 5;
+}
+
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlLinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+getValue( const MeshType& mesh,
+          const IndexType cellIndex,
+          const CoordinatesType& coordinates,
+          const Vector& u,
+          const Real& time ) const
+{
+   return ( u[ mesh.getCellXPredecessor( cellIndex ) ]
+            - 2.0 * u[ cellIndex ]
+            + u[ mesh.getCellXSuccessor( cellIndex ) ] ) * mesh.getHxSquareInverse() +
+           ( u[ mesh.getCellYPredecessor( cellIndex ) ]
+             - 2.0 * u[ cellIndex ]
+             + u[ mesh.getCellYSuccessor( cellIndex ) ] ) * mesh.getHySquareInverse();
+}
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+   template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+void
+tnlLinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+updateLinearSystem( const RealType& time,
+                    const RealType& tau,
+                    const MeshType& mesh,
+                    const IndexType& index,
+                    const CoordinatesType& coordinates,
+                    Vector& u,
+                    Vector& b,
+                    MatrixRow& matrixRow ) const
+{
+   const RealType lambdaX = tau * mesh.getHxSquareInverse();
+   const RealType lambdaY = tau * mesh.getHySquareInverse();
+   matrixRow.setElement( 0, mesh.getCellYPredecessor( index ), -lambdaY );
+   matrixRow.setElement( 1, mesh.getCellXPredecessor( index ), -lambdaX );
+   matrixRow.setElement( 2, index,                             2.0 * ( lambdaX + lambdaY ) );
+   matrixRow.setElement( 3, mesh.getCellXSuccessor( index ),   -lambdaX );
+   matrixRow.setElement( 4, mesh.getCellYSuccessor( index ),   -lambdaY );
+}
+
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+tnlString
+tnlLinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+   return tnlString( "tnlLinearDiffusion< " ) +
+          MeshType::getType() + ", " +
+          ::getType< Real >() + ", " +
+          ::getType< Index >() + " >";
+}
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlLinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
+getValue( const MeshType& mesh,
+          const IndexType cellIndex,
+          const CoordinatesType& coordinates,
+          const Vector& u,
+          const Real& time ) const
+{
+   return ( u[ mesh.getCellXPredecessor( cellIndex ) ]
+            - 2.0 * u[ cellIndex ]
+            + u[ mesh.getCellXSuccessor( cellIndex ) ] ) * mesh.getHxSquareInverse() +
+          ( u[ mesh.getCellYPredecessor( cellIndex ) ]
+            - 2.0 * u[ cellIndex ]
+            + u[ mesh.getCellYSuccessor( cellIndex ) ] ) * mesh.getHySquareInverse() +
+          ( u[ mesh.getCellZPredecessor( cellIndex ) ]
+            - 2.0 * u[ cellIndex ]
+            + u[ mesh.getCellZSuccessor( cellIndex ) ] ) * mesh.getHzSquareInverse();
+}
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+Index
+tnlLinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
+getLinearSystemRowLength( const MeshType& mesh,
+                          const IndexType& index,
+                          const CoordinatesType& coordinates ) const
+{
+   return 7;
+}
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+   template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+void
+tnlLinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
+updateLinearSystem( const RealType& time,
+                    const RealType& tau,
+                    const MeshType& mesh,
+                    const IndexType& index,
+                    const CoordinatesType& coordinates,
+                    Vector& u,
+                    Vector& b,
+                    MatrixRow& matrixRow ) const
+{
+   const RealType lambdaX = tau * mesh.getHxSquareInverse();
+   const RealType lambdaY = tau * mesh.getHySquareInverse();
+   const RealType lambdaZ = tau * mesh.getHzSquareInverse();
+   matrixRow.setElement( 0, mesh.getCellZPredecessor( index ), -lambdaZ );
+   matrixRow.setElement( 1, mesh.getCellYPredecessor( index ), -lambdaY );
+   matrixRow.setElement( 2, mesh.getCellXPredecessor( index ), -lambdaX );
+   matrixRow.setElement( 3, index,                             2.0 * ( lambdaX + lambdaY + lambdaZ ) );
+   matrixRow.setElement( 4, mesh.getCellXSuccessor( index ),   -lambdaX );
+   matrixRow.setElement( 5, mesh.getCellYSuccessor( index ),   -lambdaY );
+   matrixRow.setElement( 6, mesh.getCellZSuccessor( index ),   -lambdaZ );
+}
+
+
+
+#endif	/* TNLLINEARDIFFUSION_IMP_H */

examples/incompressible-navier-stokes/tnlNSProblem.h

+/***************************************************************************
+                          tnlHeatEquationProblem.h  -  description
+                             -------------------
+    begin                : Feb 23, 2013
+    copyright            : (C) 2013 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef TNLHEATEQUATIONPROBLEM_H_
+#define TNLHEATEQUATIONPROBLEM_H_
+
+#include <problems/tnlPDEProblem.h>
+#include <operators/diffusion/tnlLinearDiffusion.h>
+#include <array/tnlStaticArray.h>
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator = tnlLinearDiffusion< Mesh,
+                                                              typename BoundaryCondition::RealType > >
+class tnlNSProblem : public tnlPDEProblem< Mesh,
+                                                     typename DifferentialOperator::RealType,
+                                                     typename Mesh::DeviceType,
+                                                     typename DifferentialOperator::IndexType  >
+{
+   public:
+
+      typedef typename DifferentialOperator::RealType RealType;
+      typedef typename Mesh::DeviceType DeviceType;
+      typedef typename DifferentialOperator::IndexType IndexType;
+      typedef tnlPDEProblem< Mesh, RealType, DeviceType, IndexType > BaseType;
+
+      using typename BaseType::MeshType;
+      using typename BaseType::DofVectorType;
+	  
+	  enum { Dimensions = Mesh::Dimensions };
+
+      static tnlString getTypeStatic() {return tnlString( "tnlNSProblem< " ) + Mesh :: getTypeStatic() + " >";}
+
+      tnlString getPrologHeader() const{return tnlString( "NS equation" );}
+
+      void writeProlog( tnlLogger& logger,
+                        const tnlParameterContainer& parameters ) const {}
+
+      bool setup( const tnlParameterContainer& parameters );
+
+      bool setInitialCondition( const tnlParameterContainer& parameters,
+                                const MeshType& mesh,
+                                DofVectorType& dofs,
+                                DofVectorType& auxDofs );
+
+      template< typename MatrixType >
+      bool setupLinearSystem( const MeshType& mesh,
+                              MatrixType& matrix );
+
+      bool makeSnapshot( const RealType& time,
+                         const IndexType& step,
+                         const MeshType& mesh,
+                         DofVectorType& dofs,
+                         DofVectorType& auxDofs );
+
+      IndexType getDofs( const MeshType& mesh ) const;
+
+      void bindDofs( const MeshType& mesh,
+                     DofVectorType& dofs );
+
+      void getExplicitRHS( const RealType& time,
+                           const RealType& tau,
+                           const MeshType& mesh,
+                           DofVectorType& _u,
+                           DofVectorType& _fu );
+
+      template< typename MatrixType >
+      void assemblyLinearSystem( const RealType& time,
+                                 const RealType& tau,
+                                 const MeshType& mesh,
+                                 DofVectorType& dofs,
+                                 DofVectorType& auxDofs,
+                                 MatrixType& matrix,
+                                 DofVectorType& rightHandSide );
+
+
+      protected:
+
+      tnlSharedVector< RealType, DeviceType, IndexType > p;
+	  
+	  tnlStaticArray< Dimensions, tnlSharedVector< RealType, DeviceType, IndexType > > u;
+
+      DifferentialOperator differentialOperator;
+
+      BoundaryCondition boundaryCondition;
+   
+      RightHandSide rightHandSide;
+};
+
+#include <problems/tnlHeatEquationProblem_impl.h>
+
+#endif /* TNLHEATEQUATIONPROBLEM_H_ */

examples/incompressible-navier-stokes/tnlNSProblem_impl.h

+/***************************************************************************
+                          tnlHeatEquationProblem_impl.h  -  description
+                             -------------------
+    begin                : Mar 10, 2013
+    copyright            : (C) 2013 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef TNLHEATEQUATIONPROBLEM_IMPL_H_
+#define TNLHEATEQUATIONPROBLEM_IMPL_H_
+
+#include <core/mfilename.h>
+#include <matrices/tnlMatrixSetter.h>
+#include <matrices/tnlMultidiagonalMatrixSetter.h>
+#include <core/tnlLogger.h>
+#include <solvers/pde/tnlExplicitUpdater.h>
+#include <solvers/pde/tnlLinearSystemAssembler.h>
+
+
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator >
+bool
+tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+setup( const tnlParameterContainer& parameters )
+{
+   if( ! this->boundaryCondition.setup( parameters, "boundary-conditions-" ) ||
+       ! this->rightHandSide.setup( parameters, "right-hand-side-" ) )
+      return false;
+   return true;
+}
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator >
+typename tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::IndexType
+tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+getDofs( const MeshType& mesh ) const
+{
+   /****
+    * Set-up DOFs and supporting grid functions
+    */
+   return mesh.getNumberOfFaces();
+}
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator >
+typename tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::IndexType
+tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+getAuxiliaryDofs( const MeshType& mesh ) const
+{
+   /****
+    * Set-up DOFs and supporting grid functions
+    */
+   return mesh.getNumberOfCells();
+}
+
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator >
+void
+tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+bindDofs( const MeshType& mesh,
+          DofVectorType& dofVector )
+{
+   const IndexType uDofs = mesh.getNumberOfFaces() / 2;
+   this->u[ 0 ].bind( dofVector.getData(), uDofs );
+   this->u[ 1 ].bind( &dofVector.getData()[ uDofs ], uDofs );
+}
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator >
+void
+tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+bindAuxiliaryDofs( const MeshType& mesh,
+          DofVectorType& auxDofs )
+{
+   const IndexType pDofs = mesh.getNumberOfCells();
+   this->p.bind( auxDofs.getData(), pDofs );
+}
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator >
+bool
+tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+setInitialCondition( const tnlParameterContainer& parameters,
+                     const MeshType& mesh,
+                     DofVectorType& dofs,
+                     DofVectorType& auxiliaryDofs )
+{
+   this->bindDofs( mesh, dofs );
+   this->bindAuxiliaryDofs( mesh, dofs );
+   this->p.setValue( 0.0 );
+   this->u[ 0 ].setValue( 0.0 );
+   this->u[ 1 ].setValue( 0.0 );
+   
+   /*const tnlString& initialConditionFile = parameters.getParameter< tnlString >( "initial-condition" );
+   if( ! this->solution.load( initialConditionFile ) )
+   {
+      cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." << endl;
+      return false;
+   }*/
+   return true;
+}
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator >
+   template< typename MatrixType >          
+bool
+tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+setupLinearSystem( const MeshType& mesh,
+                   MatrixType& matrix )
+{
+   const IndexType dofs = this->getDofs( mesh );
+   typedef typename MatrixType::RowLengthsVector RowLengthsVectorType;
+   RowLengthsVectorType rowLengths;
+   if( ! rowLengths.setSize( dofs ) )
+      return false;
+   tnlMatrixSetter< MeshType, DifferentialOperator, BoundaryCondition, RowLengthsVectorType > matrixSetter;
+   matrixSetter.template getRowLengths< Mesh::Dimensions >( mesh,
+                                                            differentialOperator,
+                                                            boundaryCondition,
+                                                            rowLengths );
+   matrix.setDimensions( dofs, dofs );
+   if( ! matrix.setRowLengths( rowLengths ) )
+      return false;
+   return true;
+   //return tnlMultidiagonalMatrixSetter< Mesh >::setupMatrix( mesh, matrix );
+}
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator >
+bool
+tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+makeSnapshot( const RealType& time,
+              const IndexType& step,
+              const MeshType& mesh,
+              DofVectorType& dofs,
+              DofVectorType& auxiliaryDofs )
+{
+   cout << endl << "Writing output at time " << time << " step " << step << "." << endl;
+
+   this->bindDofs( mesh, dofs );
+   this->bindAuxiliaryDofs( mesh, auxiliaryDofs );
+   //cout << "dofs = " << dofs << endl;
+   tnlString fileName;
+   FileNameBaseNumberEnding( "u-", step, 5, ".tnl", fileName );
+   if( ! this->solution.save( fileName ) )
+      return false;
+   return true;
+}
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator >
+void
+tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+getExplicitRHS( const RealType& time,
+                const RealType& tau,
+                const MeshType& mesh,
+                DofVectorType& u,
+                DofVectorType& fu )
+{
+   /****
+    * If you use an explicit solver like tnlEulerSolver or tnlMersonSolver, you
+    * need to implement this method. Compute the right-hand side of
+    *
+    *   d/dt u(x) = fu( x, u )
+    *
+    * You may use supporting vectors again if you need.
+    */
+
+   //cout << "u = " << u << endl;
+   this->bindDofs( mesh, u );
+   
+   tnlExplicitUpdater< Mesh, DofVectorType, DifferentialOperator, BoundaryCondition, RightHandSide > explicitUpdater;
+   explicitUpdater.template update< Mesh::Dimensions >( time,
+                                                        mesh,
+                                                        this->differentialOperator,
+                                                        this->boundaryCondition,
+                                                        this->rightHandSide,
+                                                        u,
+                                                        fu );
+   //cout << "u = " << u << endl;
+   //cout << "fu = " << fu << endl;
+   //_u.save( "u.tnl" );
+   //_fu.save( "fu.tnl" );
+   //getchar();
+}
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator >
+    template< typename MatrixType >          
+void
+tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+assemblyLinearSystem( const RealType& time,
+                      const RealType& tau,
+                      const MeshType& mesh,
+                      DofVectorType& u,
+                      DofVectorType& auxDofs,
+                      MatrixType& matrix,
+                      DofVectorType& b )
+{
+   tnlLinearSystemAssembler< Mesh, DofVectorType, DifferentialOperator, BoundaryCondition, RightHandSide, MatrixType > systemAssembler;
+   systemAssembler.template assembly< Mesh::Dimensions >( time,
+                                                          tau,
+                                                          mesh,
+                                                          this->differentialOperator,
+                                                          this->boundaryCondition,
+                                                          this->rightHandSide,
+                                                          u,
+                                                          matrix,
+                                                          b );
+   /*matrix.print( cout );
+   cout << endl << b << endl;
+   cout << endl << u << endl;
+   abort();*/
+}
+
+#endif /* TNLHEATEQUATIONPROBLEM_IMPL_H_ */