DiffusionStep.hpp

// This file is part of the imaging2 class library.
//
// Fuchs Matthias
// University of Innsbruck, Infmath Imaging, 2008.
// http://infmath.uibk.ac.at
//
// All rights reserved.


#ifndef EQUATION_DIFFUSIONSTEP_H
#define EQUATION_DIFFUSIONSTEP_H

#include <fem/equation/SimpleEquationInterface.hpp>
#include <boost/shared_ptr.hpp>


namespace imaging
{
  template<class fem_types>
  class DiffusionStep : public SimpleEquationInterface<fem_types>
  {
    boost::shared_ptr< ublas::vector<float_t> > _input_ptr;
    boost::shared_ptr< ublas::vector<float_t> > _diffusion_ptr;
    float_t _step_size;

  public:
    typedef ublas::fixed_matrix<float_t, fem_types::data_dimension, fem_types::data_dimension> matrix_coefficient_t;
    typedef typename ublas::fixed_vector<float_t, fem_types::data_dimension> vector_coefficient_t;

    DiffusionStep() : _step_size(0.0) {}

    const ublas::vector<float_t> & input() const { return *_input_ptr; }
    
    void set_input(boost::shared_ptr< ublas::vector<float_t> > input_ptr) { _input_ptr = input_ptr; }
    
    const ublas::vector<float_t> & diffusion() const { return *_diffusion_ptr; }
    
    void set_diffusion(boost::shared_ptr< ublas::vector<float_t> > diffusion_ptr) { _diffusion_ptr = diffusion_ptr; }
    
    float_t step_size() const { return _step_size; }
    
    void set_step_size(float_t step_size) { _step_size = step_size; }
    
    static const size_t boundary_data_type = SimpleEquationInterface<fem_types>::NO_BOUNDARY_DATA;

    static const bool a_active = false;
    static const bool b_active = false;
    static const bool c_active = true;
    static const bool f_active = true;
    static const bool g_active = false;
    
    void stiffness_matrix(std::size_t integrator_node,
                          const FemKernel<fem_types> & kernel,
                          matrix_coefficient_t & A,
                          ublas::fixed_vector<float_t, fem_types::data_dimension> & a,
                          ublas::fixed_vector<float_t, fem_types::data_dimension> & b,
                          float_t & c) const
    { 
      float_t local_diffusion;
      kernel.grid().interpolate_value(integrator_node, *_diffusion_ptr, kernel, local_diffusion);
      A = _step_size * local_diffusion * ublas::identity_matrix<float_t>(fem_types::data_dimension);
      c = 1.0;
    }
    
    void force_vector(std::size_t integrator_node,
                      const FemKernel<fem_types> & kernel,
                      float_t & f,
                      ublas::fixed_vector<float_t, fem_types::data_dimension> & g) const
    {
      kernel.grid().interpolate_value(integrator_node, *_input_ptr, kernel, f);
    }
    
    bool sanity_check_stiffness_matrix(const FemKernel<fem_types> & kernel, std::string & error_message) const
    { 
      if(! _diffusion_ptr)
        return false;
        
      if(_diffusion_ptr->size() != kernel.grid().n_nodes())
        return false;
        
      return true;
    }
    
    bool sanity_check_force_vector(const FemKernel<fem_types> & kernel, std::string & error_message) const
    { 
      if(! _input_ptr)
        return false;
        
      if(_input_ptr->size() != kernel.grid().n_nodes())
        return false;
        
      return true;
    }
  };

}


#endif

---