00001 #ifndef EQUATION_TVFLOWSTEP_H
00002 #define EQUATION_TVFLOWSTEP_H
00003
00004 #include <fem/equation/SimpleEquationInterface.hpp>
00005 #include <boost/shared_ptr.hpp>
00006
00007 namespace imaging
00008 {
00018 template<class fem_types>
00019 class TvFlowStep: public SimpleEquationInterface<fem_types>
00020 {
00021 boost::shared_ptr< ublas::vector<float_t> > _input_ptr;
00022 float_t _step_size;
00023 float_t _epsilon;
00024
00025 public:
00026 typedef ublas::fixed_matrix<float_t, fem_types::data_dimension, fem_types::data_dimension> matrix_coefficient_t;
00027 typedef typename ublas::fixed_vector<float_t, fem_types::data_dimension> vector_coefficient_t;
00028
00029 TvFlowStep() : _step_size(0.0), _epsilon(1.0) {}
00030
00032 const ublas::vector<float_t> & input() const { return *_input_ptr; }
00033
00035 void set_input(boost::shared_ptr< ublas::vector<float_t> > input_ptr) { _input_ptr = input_ptr; }
00036
00038 float_t step_size() const { return _step_size; }
00039
00041 void set_step_size(float_t step_size) { _step_size = step_size; }
00042
00044 float_t epsilon() const { return _epsilon; }
00045
00047 void set_epsilon(float_t epsilon) { _epsilon = epsilon; }
00048
00049 static const bool a_active = false;
00050 static const bool b_active = false;
00051 static const bool c_active = true;
00052 static const bool f_active = true;
00053 static const bool g_active = false;
00054
00055 static const size_t boundary_data_type = SimpleEquationInterface<fem_types>::NO_BOUNDARY_DATA;
00056
00057 void stiffness_matrix(std::size_t integrator_node,
00058 const FemKernel<fem_types> & kernel,
00059 matrix_coefficient_t & A,
00060 ublas::fixed_vector<float_t, fem_types::data_dimension> & a,
00061 ublas::fixed_vector<float_t, fem_types::data_dimension> & b,
00062 float_t & c) const
00063 {
00064 ublas::fixed_vector<float_t, fem_types::data_dimension> local_gradient;
00065 kernel.grid().interpolate_gradient(integrator_node, *_input_ptr, kernel, local_gradient);
00066 A = _step_size / sqrt(inner_prod(local_gradient, local_gradient) + square(_epsilon)) *
00067 ublas::identity_matrix<float_t>(fem_types::data_dimension);
00068 c = 1.0;
00069 }
00070
00071 void force_vector(std::size_t integrator_node,
00072 const FemKernel<fem_types> & kernel,
00073 float_t & f,
00074 ublas::fixed_vector<float_t, fem_types::data_dimension> & g) const
00075 {
00076 kernel.grid().interpolate_value(integrator_node, *_input_ptr, kernel, f);
00077 }
00078
00079 bool sanity_check_stiffness_matrix(const FemKernel<fem_types> & kernel, std::string & error_message) const
00080 {
00081 if(! _input_ptr)
00082 return false;
00083
00084 if(_input_ptr->size() != kernel.grid().n_nodes())
00085 return false;
00086
00087 return true;
00088 }
00089
00090 bool sanity_check_force_vector(const FemKernel<fem_types> & kernel, std::string & error_message) const
00091 {
00092 if(! _input_ptr)
00093 return false;
00094
00095 if(_input_ptr->size() != kernel.grid().n_nodes())
00096 return false;
00097
00098 return true;
00099 }
00100
00101 };
00102
00103 }
00104
00105
00106 #endif
