diff --git a/src/TNL/Problems/HeatEquationProblem_impl.h b/src/TNL/Problems/HeatEquationProblem_impl.h
index 8215e915b843af3fe55ebc10701e567af69e418d..84492612b30af9354403633a2b0e50636fd1ce7f 100644
--- a/src/TNL/Problems/HeatEquationProblem_impl.h
+++ b/src/TNL/Problems/HeatEquationProblem_impl.h
@@ -227,11 +227,11 @@ getExplicitRHS( const RealType& time,
       this->rightHandSidePointer,
       this->uPointer,
       fuPointer );
-   /*BoundaryConditionsSetter< MeshFunctionType, BoundaryCondition > boundaryConditionsSetter;
+   Solvers::PDE::BoundaryConditionsSetter< MeshFunctionType, BoundaryCondition > boundaryConditionsSetter;
    boundaryConditionsSetter.template apply< typename Mesh::Cell >(
-      this->boundaryCondition,
+      this->boundaryConditionPointer,
       time + tau,
-      this->u );*/
+      this->uPointer );
    
    //fu.write( "fu.txt", "gnuplot" );
    //this->u.write( "u.txt", "gnuplot");
diff --git a/src/TNL/Solvers/PDE/BoundaryConditionsSetter.h b/src/TNL/Solvers/PDE/BoundaryConditionsSetter.h
index 037d5a05b39b22da86a71eca1e019f2663d9c93e..277d39c5b10c1d5ec6fbb465ce749fc62fbc6ac6 100644
--- a/src/TNL/Solvers/PDE/BoundaryConditionsSetter.h
+++ b/src/TNL/Solvers/PDE/BoundaryConditionsSetter.h
@@ -13,6 +13,7 @@
 
 #include <TNL/Devices/Cuda.h>
 #include <TNL/Functions/FunctionAdapter.h>
+#include <TNL/SharedPointer.h>
 
 namespace TNL {
 namespace Solvers {
@@ -55,11 +56,14 @@ class BoundaryConditionsSetter
          RealType,
          MeshFunction,
          BoundaryConditions > TraverserUserData;
+      typedef SharedPointer< MeshType, DeviceType > MeshPointer;
+      typedef SharedPointer< BoundaryConditions, DeviceType > BoundaryConditionsPointer;
+      typedef SharedPointer< MeshFunction, DeviceType > MeshFunctionPointer;
 
       template< typename EntityType = typename MeshType::Cell >
-      static void apply( const BoundaryConditions& boundaryConditions,
+      static void apply( const BoundaryConditionsPointer& boundaryConditions,
                          const RealType& time,
-                         MeshFunction& u );
+                         MeshFunctionPointer& u );
  
       class TraverserBoundaryEntitiesProcessor
       {
diff --git a/src/TNL/Solvers/PDE/BoundaryConditionsSetter_impl.h b/src/TNL/Solvers/PDE/BoundaryConditionsSetter_impl.h
index 26ea5a4a7c531182f6cfe797be16dd67d68db9db..779a09e543677fa480078685867c8d8711e94896 100644
--- a/src/TNL/Solvers/PDE/BoundaryConditionsSetter_impl.h
+++ b/src/TNL/Solvers/PDE/BoundaryConditionsSetter_impl.h
@@ -22,20 +22,22 @@ template< typename MeshFunction,
    template< typename EntityType >
 void
 BoundaryConditionsSetter< MeshFunction, BoundaryConditions >::
-apply( const BoundaryConditions& boundaryConditions,
+apply( const BoundaryConditionsPointer& boundaryConditions,
        const RealType& time,
-       MeshFunction& u )
+       MeshFunctionPointer& u )
 {
-   if( std::is_same< DeviceType, Devices::Host >::value )
+   //if( std::is_same< DeviceType, Devices::Host >::value )
    {
-      TraverserUserData userData( time, boundaryConditions, u );
+      TraverserUserData userData( time,
+                                  boundaryConditions.template getData< DeviceType >(),
+                                  u.template modifyData< DeviceType >() );
       Meshes::Traverser< MeshType, EntityType > meshTraverser;
       meshTraverser.template processBoundaryEntities< TraverserUserData,
                                                       TraverserBoundaryEntitiesProcessor >
-                                                    ( u.getMeshPointer(),
+                                                    ( u->getMeshPointer(),
                                                       userData );
    }
-   if( std::is_same< DeviceType, Devices::Cuda >::value )
+   /*if( std::is_same< DeviceType, Devices::Cuda >::value )
    {
       RealType* kernelTime = Devices::Cuda::passToDevice( time );
       BoundaryConditions* kernelBoundaryConditions = Devices::Cuda::passToDevice( boundaryConditions );
@@ -51,7 +53,7 @@ apply( const BoundaryConditions& boundaryConditions,
       Devices::Cuda::freeFromDevice( kernelBoundaryConditions );
       Devices::Cuda::freeFromDevice( kernelU );
       checkCudaDevice;
-   }
+   }*/
 }
 
 } // namespace PDE
diff --git a/src/Tools/tnl-quickstart/operator-grid-specialization.h.in b/src/Tools/tnl-quickstart/operator-grid-specialization.h.in
index a65f3f9548d3908f1fc940f6755a6c37f0397a8b..cac8dbb2aa04435388ccccf4bac205d9bae55dd6 100644
--- a/src/Tools/tnl-quickstart/operator-grid-specialization.h.in
+++ b/src/Tools/tnl-quickstart/operator-grid-specialization.h.in
@@ -30,14 +30,14 @@ class {operatorName}< TNL::Meshes::Grid< {meshDimensions}, MeshReal, Device, Mes
 
       template< typename MeshEntity, typename Vector, typename MatrixRow >
       __cuda_callable__
-      void updateLinearSystem( const RealType& time,
-                               const RealType& tau,
-                               const MeshType& mesh,
-                               const IndexType& index,
-                               const MeshEntity& entity,
-                               const MeshFunctionType& u,
-                               Vector& b,
-                               MatrixRow& matrixRow ) const;
+      void setMatrixElements( const RealType& time,
+                              const RealType& tau,
+                              const MeshType& mesh,
+                              const IndexType& index,
+                              const MeshEntity& entity,
+                              const MeshFunctionType& u,
+                              Vector& b,
+                              MatrixRow& matrixRow ) const;
 }};
     
 
diff --git a/src/Tools/tnl-quickstart/operator-grid-specialization_impl.h.in b/src/Tools/tnl-quickstart/operator-grid-specialization_impl.h.in
index 2e975772fce8d3966c1ba6b102c0bcccee1b062c..d9e0a7f856644140e088b13ea0074e57ac1787d1 100644
--- a/src/Tools/tnl-quickstart/operator-grid-specialization_impl.h.in
+++ b/src/Tools/tnl-quickstart/operator-grid-specialization_impl.h.in
@@ -73,14 +73,14 @@ template< typename MeshReal,
 __cuda_callable__        
 void
 {operatorName}< TNL::Meshes::Grid< {meshDimensions}, MeshReal, Device, MeshIndex >, Real, Index >::
-updateLinearSystem( const RealType& time,
-                    const RealType& tau,
-                    const MeshType& mesh,
-                    const IndexType& index,
-                    const MeshEntity& entity,
-                    const MeshFunctionType& u,
-                    Vector& b,
-                    MatrixRow& matrixRow ) const
+setMatrixElements( const RealType& time,
+                   const RealType& tau,
+                   const MeshType& mesh,
+                   const IndexType& index,
+                   const MeshEntity& entity,
+                   const MeshFunctionType& u,
+                   Vector& b,
+                   MatrixRow& matrixRow ) const
 {{
    /****
     * Setup the non-zero elements of the linear system here.
diff --git a/src/Tools/tnl-quickstart/problem_impl.h.in b/src/Tools/tnl-quickstart/problem_impl.h.in
index 34aeb6dd6ae8fda9c4ab08aac89d6f9752d92e21..0a495676127f6e5a5f85402ebc7c994d19c79c36 100644
--- a/src/Tools/tnl-quickstart/problem_impl.h.in
+++ b/src/Tools/tnl-quickstart/problem_impl.h.in
@@ -53,8 +53,8 @@ setup( const MeshPointer& meshPointer,
        const TNL::Config::ParameterContainer& parameters,
        const TNL::String& prefix )
 {{
-   if( ! this->boundaryCondition.setup( meshPointer, parameters, "boundary-conditions-" ) ||
-       ! this->rightHandSide.setup( parameters, "right-hand-side-" ) )
+   if( ! this->boundaryCondition->setup( meshPointer, parameters, "boundary-conditions-" ) ||
+       ! this->rightHandSide->setup( parameters, "right-hand-side-" ) )
       return false;
    return true;
 }}
@@ -179,8 +179,8 @@ getExplicitRHS( const RealType& time,
 
    this->bindDofs( mesh, _u );
    TNL::Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, DifferentialOperator, BoundaryCondition, RightHandSide > explicitUpdater;
-   MeshFunctionType u( mesh, _u ); 
-   MeshFunctionType fu( mesh, _fu ); 
+   TNL::SharedPointer< MeshFunctionType > u( mesh, _u ); 
+   TNL::SharedPointer< MeshFunctionType > fu( mesh, _fu ); 
    explicitUpdater.template update< typename Mesh::Cell >( time,
                                                            mesh,
                                                            this->differentialOperator,
@@ -219,7 +219,7 @@ assemblyLinearSystem( const RealType& time,
                              typename MatrixPointer::ObjectType,
                              DofVectorType > systemAssembler;
 
-   TNL::Functions::MeshFunction< Mesh > u( mesh, _u );
+   TNL::SharedPointer< TNL::Functions::MeshFunction< Mesh > > u( mesh, _u );
    systemAssembler.template assembly< typename Mesh::Cell >( time,
                                                              tau,
                                                              mesh,