diff --git a/src/operators/diffusion/tnlExactNonlinearDiffusion.h b/src/operators/diffusion/tnlExactNonlinearDiffusion.h
index 1522343b231f2735e99d54d463191c218f9e65ab..b8e093da7c87c2a1d1545756087a81896946a815 100644
--- a/src/operators/diffusion/tnlExactNonlinearDiffusion.h
+++ b/src/operators/diffusion/tnlExactNonlinearDiffusion.h
@@ -20,71 +20,97 @@
 
 #include <functions/tnlDomain.h>
 
-template< typename OperatorQ, int Dimensions >
+template< typename Nonlinearity, int Dimensions >
 class tnlExactNonlinearDiffusion
 {};
 
-template< typename OperatorQ >
-class tnlExactNonlinearDiffusion< OperatorQ, 1 > : public tnlDomain< 1, SpaceDomain >
+template< typename Nonlinearity >
+class tnlExactNonlinearDiffusion< Nonlinearity, 1 > 
+   : public tnlDomain< 1, SpaceDomain >
 {
    public:
 
-      enum { Dimensions = 1 };
-
-      static tnlString getType();
-   
-#ifdef HAVE_NOT_CXX11      
-      template< typename Function, typename Vertex, typename Real >
-#else   
-      template< typename Function, typename Vertex, typename Real = typename Vertex::RealType >
-#endif      
+      static tnlString getType()
+      {
+         return "tnlExactNonlinearDiffusion< " + Nonlinearity::getType() + ", 1 >";
+      };
       
+      void setNonlinearity( const Nonlinearity& nonlinearity )
+      {
+         this->nonlinearity = nonlinearity;
+      }
+   
+      template< typename Function >
       __cuda_callable__
       Real operator()( const Function& function,
                        const Vertex& v,
-                       const Real& time = 0.0 ) const;
+                       const Real& time = 0.0 ) const
+      {         
+         const Real u_x = function.template getPartialDerivative< 1, 0, 0 >( v, time );
+         const Real u_xx = function.template getPartialDerivative< 2, 0, 0 >( v, time );
+         const Real g = nonlinearity( function, v, time ) 
+         const Real g_x = nonlinearity::template getPartialDerivative< Function, 1, 0, 0 >( function, v, time );
+         return u_xx - u_x * g_x / g;          
+      }
+   
+      protected:
+         
+         Nonlinearity nonlinearity;
 };
 
-template< typename OperatorQ >
-class tnlExactNonlinearDiffusion< OperatorQ, 2 > : public tnlDomain< 2, SpaceDomain >
+template< typename Nonlinearity >
+class tnlExactNonlinearDiffusion< Nonlinearity, 2 >
+   : public tnlDomain< 2, SpaceDomain >
 {
    public:
 
-      enum { Dimensions = 2 };
-
-      static tnlString getType();
+      static tnlString getType()
+      {
+         return "tnlExactNonlinearDiffusion< " + Nonlinearity::getType() + ", 2 >";
+      };
 
-#ifdef HAVE_NOT_CXX11      
-      template< typename Function, typename Vertex, typename Real >
-#else   
-      template< typename Function, typename Vertex, typename Real = typename Vertex::RealType >
-#endif 
+      void setNonlinearity( const Nonlinearity& nonlinearity )
+      {
+         this->nonlinearity = nonlinearity;
+      }
 
+      template< typename Function >
       __cuda_callable__
       Real operator()( const Function& function,
-                       const Vertex& v,
+                       const VertexType& v,
                        const Real& time = 0.0 ) const;
+
+      protected:
+         
+         Nonlinearity nonlinearity;
+      
 };
 
-template< typename OperatorQ >
-class tnlExactNonlinearDiffusion< OperatorQ, 3 > : public tnlDomain< 3, SpaceDomain >
+template< typename Nonlinearity >
+class tnlExactNonlinearDiffusion< Nonlinearity, 3 >
+   : public tnlDomain< 3, SpaceDomain >
 {
    public:
 
-      enum { Dimensions = 3 };
-
-      static tnlString getType();
-
-#ifdef HAVE_NOT_CXX11      
-      template< typename Function, typename Vertex, typename Real >
-#else   
-      template< typename Function, typename Vertex, typename Real = typename Vertex::RealType >
-#endif 
+      static tnlString getType()
+      {
+         return "tnlExactNonlinearDiffusion< " + Nonlinearity::getType() + ", 3 >";
+      }
 
+      void setNonlinearity( const Nonlinearity& nonlinearity )
+      {
+         this->nonlinearity = nonlinearity;
+      }      
+      
+      template< typename Function >
       __cuda_callable__
       Real operator()( const Function& function,
-                       const Vertex& v,
+                       const VertexType& v,
                        const Real& time = 0.0 ) const;
+      
+      protected:
+         
+         Nonlinearity nonlinearity;
 };
 
 #include "tnlExactNonlinearDiffusion_impl.h"
diff --git a/src/operators/diffusion/tnlExactNonlinearDiffusion_impl.h b/src/operators/diffusion/tnlExactNonlinearDiffusion_impl.h
index 1444f124505e48bee71dfd82132a7819582eb71f..293e8ab52c37fc6b6ecaa9e425ebb46cb16c1ea0 100644
--- a/src/operators/diffusion/tnlExactNonlinearDiffusion_impl.h
+++ b/src/operators/diffusion/tnlExactNonlinearDiffusion_impl.h
@@ -18,73 +18,73 @@
 #ifndef TNLEXACTNONLINEARDIFFUSION_IMPL_H_
 #define TNLEXACTNONLINEARDIFFUSION_IMPL_H_
 
-template< typename OperatorQ >
+template< typename Nonlinearity >
 tnlString
-tnlExactNonlinearDiffusion< OperatorQ, 1 >::
+tnlExactNonlinearDiffusion< Nonlinearity, 1 >::
 getType()
 {
-   return "tnlExactNonlinearDiffusion< " + OperatorQ::getType() + ", 1 >";
+   return "tnlExactNonlinearDiffusion< " + Nonlinearity::getType() + ", 1 >";
 }
 
-template< typename OperatorQ >
+template< typename Nonlinearity >
 template< typename Function, typename Vertex, typename Real >
 __cuda_callable__
 Real
-tnlExactNonlinearDiffusion< OperatorQ, 1 >::
+tnlExactNonlinearDiffusion< Nonlinearity, 1 >::
 operator()( const Function& function,
           const Vertex& v,
           const Real& time ) const
 {
    return function.template getPartialDerivative< 2, 0, 0 >( v, time ) - function.template getPartialDerivative< 1, 0, 0 >( v, time ) 
-          * OperatorQ::template getPartialDerivative<1, 0, 0>(function, v, time ) / OperatorQ::template getPartialDerivative<0, 0, 0>(function, v, time );
+          * Nonlinearity::template getPartialDerivative< Function, 1, 0, 0>(function, v, time ) / Nonlinearity( function, v, time );
 }
 
-template< typename OperatorQ >
+template< typename Nonlinearity >
 tnlString
-tnlExactNonlinearDiffusion< OperatorQ, 2 >::
+tnlExactNonlinearDiffusion< Nonlinearity, 2 >::
 getType()
 {
-   return "tnlExactNonlinearDiffusion< " + OperatorQ::getType() + ", 2 >";
+   return "tnlExactNonlinearDiffusion< " + Nonlinearity::getType() + ", 2 >";
 }
 
-template< typename OperatorQ >
+template< typename Nonlinearity >
 template< typename Function, typename Vertex, typename Real >
 __cuda_callable__
 Real
-tnlExactNonlinearDiffusion< OperatorQ, 2 >::
+tnlExactNonlinearDiffusion< Nonlinearity, 2 >::
 operator()( const Function& function,
           const Vertex& v,
           const Real& time ) const
 {
    return  function.template getPartialDerivative< 2, 0, 0 >( v, time ) +  function.template getPartialDerivative< 0, 2, 0 >( v, time )
-           -( OperatorQ::template getPartialDerivative<1, 0, 0> (function, v, time) * function.template getPartialDerivative< 1, 0, 0 >( v, time ) 
-           + OperatorQ::template getPartialDerivative<0, 1, 0> (function, v, time) * function.template getPartialDerivative< 0, 1, 0 >( v, time ) ) 
-           / OperatorQ::template getPartialDerivative<0, 0, 0> (function, v, time);
+           -( Nonlinearity::template getPartialDerivative<1, 0, 0> (function, v, time) * function.template getPartialDerivative< 1, 0, 0 >( v, time ) 
+           + Nonlinearity::template getPartialDerivative<0, 1, 0> (function, v, time) * function.template getPartialDerivative< 0, 1, 0 >( v, time ) ) 
+           / Nonlinearity::template getPartialDerivative<0, 0, 0> (function, v, time);
 }
 
-template< typename OperatorQ >
+template< typename Nonlinearity >
 tnlString
-tnlExactNonlinearDiffusion< OperatorQ, 3 >::
+tnlExactNonlinearDiffusion< Nonlinearity, 3 >::
 getType()
 {
-   return "tnlExactNonlinearDiffusion< " + OperatorQ::getType() + ", 3 >";
+   return "tnlExactNonlinearDiffusion< " + Nonlinearity::getType() + ", 3 >";
 }
 
-template< typename OperatorQ >
+template< typename Nonlinearity >
 template< typename Function, typename Vertex, typename Real >
 __cuda_callable__
 Real
-tnlExactNonlinearDiffusion< OperatorQ, 3 >::
+tnlExactNonlinearDiffusion< Nonlinearity, 3 >::
 operator()( const Function& function,
           const Vertex& v,
           const Real& time ) const
 {
    return  function.template getPartialDerivative< 2, 0, 0 >( v, time ) +  function.template getPartialDerivative< 0, 2, 0 >( v, time )
            +  function.template getPartialDerivative< 0, 0, 2 >( v, time )
-           -( OperatorQ::template getPartialDerivative<1, 0, 0> (function, v, time) * function.template getPartialDerivative< 1, 0, 0 >( v, time ) 
-           + OperatorQ::template getPartialDerivative<0, 1, 0> (function, v, time) * function.template getPartialDerivative< 0, 1, 0 >( v, time )
-           + OperatorQ::template getPartialDerivative<0, 0, 1> (function, v, time) * function.template getPartialDerivative< 0, 0, 1 >( v, time ) )
-           / OperatorQ::template getPartialDerivative<0, 0, 0> (function, v, time);
+           -( Nonlinearity::template getPartialDerivative<1, 0, 0> (function, v, time) * function.template getPartialDerivative< 1, 0, 0 >( v, time ) 
+           + Nonlinearity::template getPartialDerivative<0, 1, 0> (function, v, time) * function.template getPartialDerivative< 0, 1, 0 >( v, time )
+           + Nonlinearity::template getPartialDerivative<0, 0, 1> (function, v, time) * function.template getPartialDerivative< 0, 0, 1 >( v, time ) )
+           / Nonlinearity::template getPartialDerivative<0, 0, 0> (function, v, time);
 }
 
 #endif /* TNLEXACTNONLINEARDIFFUSION_IMPL_H_ */