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

Adding total variation minimization operator.

parent cc518b8f
Loading
Loading
Loading
Loading
+6 −0
Original line number Diff line number Diff line
TODO:
 - implementovat tnlMixedGridBoundaryConditions, kde by se pro kazdou stranu gridu definoval jiny zvlastni typ
   okrajovych podminek
 - dalo by se tim resit i skladani zpetnych a doprednych diferenci u nelinearni difuze, kdy je potreba napr. dopredne diference
   vycislit i na leve a dolni hranici 2D gridu

TODO:
 - implementovat tuple pro snazsi a snad efektoivnejsi prenos dat na GPU
 - nebylo by nutne definovat pomocne datove structury pro traverser
+21 −29
Original line number Diff line number Diff line
@@ -49,8 +49,7 @@ class tnlOneSidedNonlinearDiffusion< tnlGrid< 1,MeshReal, Device, MeshIndex >, N
      typedef Index IndexType;
      typedef Nonlinearity NonlinearityType;
      typedef typename MeshType::template MeshEntity< MeshType::getMeshDimensions() > CellType;
      template< typename NonlinearityFunction >
      using ExactOperatorType = tnlExactNonlinearDiffusion< NonlinearityFunction, MeshType::getMeshDimensions() > ExactOperatorType;
      typedef tnlExactNonlinearDiffusion< typename Nonlinearity::ExactOperatorType, MeshType::getMeshDimensions() > ExactOperatorType;

      tnlOneSidedNonlinearDiffusion( const Nonlinearity& nonlinearity )
      : nonlinearity( nonlinearity ){}
@@ -59,17 +58,17 @@ class tnlOneSidedNonlinearDiffusion< tnlGrid< 1,MeshReal, Device, MeshIndex >, N
      {
         return tnlString( "tnlOneSidedNonlinearDiffusion< " ) +
            MeshType::getType() + ", " +
            NonlinearDiffusionOperator::getType() + "," +
            Nonlinearity::getType() + "," +
            ::getType< Real >() + ", " +
            ::getType< Index >() + " >";         
      }     

      template< typename MeshEntity,
                typename Vector >
      template< typename MeshFunction,
                typename MeshEntity >
      __cuda_callable__
      Real operator()( const MeshEntity& entity,
                       const Vector& u,
                       const RealType& time) const
      Real operator()( const MeshFunction& u,
                       const MeshEntity& entity,
                       const RealType& time = 0.0 ) const
      {
         const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities();
         const typename MeshEntity::MeshType& mesh = entity.getMesh();
@@ -146,9 +145,7 @@ class tnlOneSidedNonlinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >,
      typedef Device DeviceType;
      typedef Index IndexType;
      typedef Nonlinearity NonlinearityType;
      template< typename NonlinearityFunction >
      using ExactOperatorType = tnlExactNonlinearDiffusion< NonlinearityFunction, MeshType::getMeshDimensions() > ExactOperatorType;
      
      typedef tnlExactNonlinearDiffusion< typename Nonlinearity::ExactOperatorType, MeshType::getMeshDimensions() > ExactOperatorType;      

      tnlOneSidedNonlinearDiffusion( const Nonlinearity& nonlinearity )
      : nonlinearity( nonlinearity ){}      
@@ -157,17 +154,17 @@ class tnlOneSidedNonlinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >,
      {
         return tnlString( "tnlOneSidedNonlinearDiffusion< " ) +
            MeshType::getType() + ", " +
            NonlinearDiffusionOperator::getType() + "," +
            Nonlinearity::getType() + "," +
            ::getType< Real >() + ", " +
            ::getType< Index >() + " >";         
      }      

      template< typename MeshEntity,
                typename Vector >
      template< typename MeshFunction,
                typename MeshEntity >
      __cuda_callable__
      Real operator()( const MeshEntity& entity,
                       const Vector& u,
                       const RealType& time) const
      Real operator()( const MeshFunction& u,
                       const MeshEntity& entity,
                       const RealType& time = 0.0 ) const
      {
         const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities();
         const typename MeshEntity::MeshType& mesh = entity.getMesh();
@@ -259,8 +256,7 @@ class tnlOneSidedNonlinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >,
      typedef Device DeviceType;
      typedef Index IndexType;
      typedef Nonlinearity NonlinearityType;
      template< typename NonlinearityFunction >
      using ExactOperatorType = tnlExactNonlinearDiffusion< NonlinearityFunction, MeshType::getMeshDimensions() > ExactOperatorType;      
      typedef tnlExactNonlinearDiffusion< typename Nonlinearity::ExactOperatorType, MeshType::getMeshDimensions() > ExactOperatorType;

      tnlOneSidedNonlinearDiffusion( const Nonlinearity& nonlinearity )
      : nonlinearity( nonlinearity ){}
@@ -269,17 +265,17 @@ class tnlOneSidedNonlinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >,
      {
         return tnlString( "tnlOneSidedNonlinearDiffusion< " ) +
            MeshType::getType() + ", " +
            NonlinearDiffusionOperator::getType() + "," +
            Nonlinearity::getType() + "," +
            ::getType< Real >() + ", " +
            ::getType< Index >() + " >";         
      }

      template< typename MeshEntity,
                typename MeshFunction >
      template< typename MeshFunction,
                typename MeshEntity >
      __cuda_callable__
      Real operator()( const MeshEntity& entity,
                       const MeshFunction& u,
                       const RealType& time) const
      Real operator()( const MeshFunction& u,
                       const MeshEntity& entity,
                       const RealType& time = 0.0 ) const
      {
         const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities();
         const typename MeshEntity::MeshType& mesh = entity.getMesh();
@@ -373,8 +369,4 @@ class tnlOneSidedNonlinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >,
      Nonlinearity& nonlinearity;
};


#include "tnlOneSidedNonlinearDiffusion_impl.h"


#endif	/* TNLONESIDEDNONLINEARDIFFUSION_H */
+119 −0
Original line number Diff line number Diff line
/***************************************************************************
                          tnlOneSidedTotalVariationMinimization.h  -  description
                             -------------------
    begin                : Feb 17, 2016
    copyright            : (C) 2016 by 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 TNLONESIDEDTOTALVARIATIONMINIMIZATION_H
#define	TNLONESIDEDTOTALVARIATIONMINIMIZATION_H

#include <operators/tnlOperator.h>
#include <operators/tnlFunctionInverseOperator.h>
#include <operators/geometric/tnlFDMGradientNorm.h>
#include <operators/tnlNeumannBoundaryConditions.h>
#include <operators/diffusion/tnlOneSidedNonlinearDiffusion.h>
#include <functions/tnlOperatorFunction.h>
#include <functions/tnlConstantFunction.h>

template< typename Mesh,
          typename Real = typename Mesh::RealType,
          typename Index = typename Mesh::IndexType,
          bool EvaluateNonlinearityOnFly = false >
class tnlOneSidedTotalVariationMinimization
   : public tnlOperator< Mesh, MeshInteriorDomain, Mesh::getMeshDimensions(), Mesh::getMeshDimensions(), Real, Index >
{
   public:
            
      typedef Mesh MeshType;
      typedef Real RealType;
      typedef Index IndexType;
      typedef tnlFDMGradientNorm< MeshType, tnlForwardFiniteDifference, RealType, IndexType > GradientNorm;
      typedef tnlFunctionInverseOperator< GradientNorm > NonlinearityOperator;
      typedef tnlMeshFunction< MeshType, MeshType::getMeshDimensions(), RealType > NonlinearityMeshFunction;
      typedef tnlConstantFunction< MeshType::getMeshDimensions(), RealType > NonlinearityBoundaryConditionsFunction;
      typedef tnlNeumannBoundaryConditions< MeshType, NonlinearityBoundaryConditionsFunction > NonlinearityBoundaryConditions;
      typedef tnlOperatorFunction< NonlinearityOperator, NonlinearityMeshFunction, NonlinearityBoundaryConditions, EvaluateNonlinearityOnFly > Nonlinearity;
      typedef tnlOneSidedNonlinearDiffusion< Mesh, Nonlinearity, RealType, IndexType > NonlinearDiffusion;
      
      tnlOneSidedTotalVariationMinimization()
      : nonlinearity( gradientNorm, nonlinearityBoundaryConditions ),
        nonlinearDiffusion( nonlinearity ){}
      
      static tnlString getType()
      {
         return tnlString( "tnlOneSidedTotalVariationMinimization< " ) +
            MeshType::getType() + ", " +
            ::getType< Real >() + ", " +
            ::getType< Index >() + " >";         
      }
      
      bool refresh( const RealType& time = 0.0 )
      {
         return this->nonlinearity.refresh( time );
      }
      
      bool deepRefresh( const RealType& time = 0.0 )
      {
         return this->nonlinearity.deepRefresh( time );
      }      
      
      template< typename MeshFunction,
                typename MeshEntity >
      __cuda_callable__
      Real operator()( const MeshFunction& u,
                       const MeshEntity& entity,
                       const RealType& time = 0.0 ) const
      {
         return this->nonlinearDiffusion( u, entity, time );
      }

      template< typename MeshEntity >
      __cuda_callable__
      Index getLinearSystemRowLength( const MeshType& mesh,
                                      const IndexType& index,
                                      const MeshEntity& entity ) const
      {
         return this->nonlinearDiffusion.getLinearSystemRowLength( mesh, index, entity );
      }

      template< typename MeshEntity,
                typename MeshFunction,
                typename Vector,
                typename Matrix >
      __cuda_callable__
      void updateLinearSystem( const RealType& time,
                               const RealType& tau,
                               const MeshType& mesh,
                               const IndexType& index,
                               const MeshEntity& entity,
                               const MeshFunction& u,
                               Vector& b,
                               Matrix& matrix ) const
      {
         this->nonlinearDiffusion.updateLinearSystem( time, tau, mesh, index, entity, u, b, matrix );
      }            
      
   protected:      
      
      NonlinearityBoundaryConditions nonlinearityBoundaryConditions;
      
      GradientNorm gradientNorm;
      
      Nonlinearity nonlinearity;
      
      NonlinearDiffusion nonlinearDiffusion;            
};

#endif	/* TNLONESIDEDTOTALVARIATIONMINIMIZATION_H */
+61 −0
Original line number Diff line number Diff line
/***************************************************************************
                          tnlFunctionInverseOperator.h  -  description
                             -------------------
    begin                : Feb 17, 2016
    copyright            : (C) 2016 by 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 TNLFUNCTIONINVERSEOPERATOR_H
#define	TNLFUNCTIONINVERSEOPERATOR_H

template< typename Operator >
class tnlFunctionInverseOperator
{
   public:
      
      typedef Operator OperatorType;
      typedef typename OperatorType::RealType RealType;
      typedef typename OperatorType::IndexType IndexType;
      typedef tnlFunctionInverseOperator< Operator > ThisType;
      typedef ThisType ExactOperatorType;
      
      
      tnlFunctionInverseOperator( const OperatorType& operator_ )
      : operator_( operator_ ) {};
      
      static tnlString getType()
      {
         return tnlString( "tnlFunctionInverseOperator< " ) + Operator::getType() + " >";         
      }
      
      const OperatorType& getOperator() const { return this->operator_; }
      
      template< typename MeshFunction,
                typename MeshEntity >
      __cuda_callable__
      typename MeshFunction::RealType
      operator()( const MeshFunction& u,
                  const MeshEntity& entity,
                  const typename MeshFunction::RealType& time = 0.0 ) const
      {
         return 1.0 / operator_( u, entity, time );
      }
      
   protected:
      
      const OperatorType& operator_;
};


#endif	/* TNLINVERSEFUNCTIONOPERATOR_H */
+1 −0
Original line number Diff line number Diff line
@@ -3,6 +3,7 @@ TARGET_LINK_LIBRARIES( tnlLinearDiffusionTest${mpiExt}${debugExt} ${CPPUNIT_LIBR
                                                                 tnl${mpiExt}${debugExt}-0.1 )
ADD_EXECUTABLE( tnlOneSidedNonlinearDiffusionTest${mpiExt}${debugExt} tnlOneSidedNonlinearDiffusionTest.cpp )
TARGET_LINK_LIBRARIES( tnlOneSidedNonlinearDiffusionTest${mpiExt}${debugExt} ${CPPUNIT_LIBRARIES}
                                                                             tnl${mpiExt}${debugExt}-0.1 )


if( BUILD_CUDA )                                                           
Loading