FemKernel.hpp

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


#ifndef FEM_FEMKERNEL_H
#define FEM_FEMKERNEL_H

#include <core/imaging2.hpp>
#include <core/matrix_utilities.hpp>

namespace imaging
{

  template<class fem_types> class Grid;

  template<std::size_t M, std::size_t N>
  float_t transform_det(const ublas::fixed_matrix<float_t, M, N> & derivative)
  {
    return fabs(determinant(derivative));
  }

  template<class fem_types>
  class FemKernel
  {
  private:
    static const std::size_t data_dimension = fem_types::data_dimension;

    static const std::size_t n_element_nodes = fem_types::shape_function_t::n_element_nodes;
    static const std::size_t n_element_vertices = fem_types::transform_t::n_element_vertices;

    typedef typename fem_types::shape_function_t shape_function_t;

    typedef typename fem_types::transform_t transform_t;
    typedef typename fem_types::integrator_t integrator_t;
    typedef typename fem_types::boundary_integrator_t boundary_integrator_t;
    
    ublas::fixed_vector< ublas::fixed_vector<float_t, n_element_nodes>, integrator_t::n_nodes> _shape_values;
    ublas::fixed_vector< ublas::fixed_vector< ublas::fixed_vector<float_t, data_dimension>, n_element_nodes >, integrator_t::n_nodes > _shape_gradients;
    ublas::fixed_vector<float_t, integrator_t::n_nodes> _transform_determinants;

    ublas::fixed_vector< ublas::fixed_vector<float_t, n_element_nodes>, boundary_integrator_t::n_nodes > _shape_boundary_values;
    ublas::fixed_vector< ublas::fixed_vector<float_t, n_element_nodes>, boundary_integrator_t::n_nodes > _shape_boundary_derivatives;
    ublas::fixed_vector<float_t, boundary_integrator_t::n_nodes> _boundary_transform_determinants;
    
    ublas::fixed_vector<float_t, data_dimension> _boundary_normal;
    ublas::fixed_vector< ublas::fixed_vector<float_t, data_dimension>, shape_function_t::n_shape_face_nodes > _boundary_normals_at_nodes;
    

    const Grid<fem_types> & _grid;
    std::size_t _current_element;
    std::size_t _current_boundary_element;

  public:

    FemKernel(const Grid<fem_types> & grid) : _grid(grid) {}

    const Grid<fem_types> & grid() const { return _grid; };
    
    void set_element(std::size_t element);
    
    void set_boundary_element(std::size_t boundary_element);

    void lazy_set_element(std::size_t element) { _current_element = element; }

    float_t shape_value(std::size_t integrator_node, std::size_t element_node) const
      { return _shape_values(integrator_node)(element_node); }
    
    const ublas::fixed_vector<float_t, fem_types::data_dimension> & shape_gradient(std::size_t integrator_node, std::size_t element_node) const
      { return _shape_gradients(integrator_node)(element_node); }
    
    float_t transform_determinant(std::size_t integrator_node) const
      { return _transform_determinants(integrator_node); }
  
    float_t shape_boundary_value(std::size_t integrator_node, std::size_t element_node) const
      { return _shape_boundary_values(integrator_node)(element_node); }

    float_t shape_boundary_derivative(std::size_t integrator_node, std::size_t element_node) const
      { return _shape_boundary_derivatives(integrator_node)(element_node); }

    float_t boundary_transform_determinant(std::size_t integrator_node) const
      { return _boundary_transform_determinants(integrator_node); }
      
    const ublas::fixed_vector<float_t, fem_types::data_dimension> & boundary_normal() const
      { return _boundary_normal; }

    std::size_t current_element() const { return _current_element; }
    
    std::size_t current_boundary_element() const { return _current_boundary_element; }
    
    bool is_boundary_node(std::size_t element_node) const { return _grid.is_boundary_node(_current_element,element_node); }
    
    const ublas::fixed_vector<float_t, fem_types::data_dimension> boundary_normal_at_node(std::size_t element_node) const { return _grid.boundary_normal(_current_element, element_node); }
  }
  ;

  template <class fem_types>
  void FemKernel<fem_types>::set_element(std::size_t element)
  {
    transform_t transform;
    shape_function_t shape_function;
    integrator_t integrator;

    _current_element = element;

    _grid.element_transform(element, transform);

    for(std::size_t integrator_node = 0; integrator_node < integrator_t::n_nodes; ++integrator_node)
    {
      ublas::fixed_matrix<float_t, data_dimension, data_dimension> derivative, inverse_derivative;
      ublas::fixed_vector<float_t, data_dimension> temp;

      transform.derivative(integrator.node(integrator_node), derivative);
      inverse_derivative = inverse(derivative);

      _transform_determinants(integrator_node) = transform_det(derivative);

      for(std::size_t element_node = 0; element_node < n_element_nodes; ++element_node)
      {
        _shape_values(integrator_node)(element_node) = shape_function.value(element_node, integrator.node(integrator_node));
        shape_function.gradient(element_node, integrator.node(integrator_node), temp);
        _shape_gradients(integrator_node)(element_node) = prod(temp, inverse_derivative);
      }
    }

  }



  template <class fem_types>
  void FemKernel<fem_types>::set_boundary_element(std::size_t element)
  {
    transform_t transform;
    shape_function_t shape_function;
    boundary_integrator_t integrator;
    std::size_t parent_element = _grid.parent_element(element);
    std::size_t parent_element_face = _grid.parent_element_face(element);

    _current_boundary_element = element;

    _grid.element_transform(parent_element, transform);

    transform.boundary_normal(parent_element_face, _boundary_normal);
    
    for(std::size_t integrator_node = 0; integrator_node < boundary_integrator_t::n_nodes; ++integrator_node)
    {
      ublas::fixed_matrix<float_t, data_dimension, data_dimension> derivative, inverse_derivative;
      ublas::fixed_matrix<float_t, data_dimension, data_dimension - 1> boundary_transform_derivative;
      ublas::fixed_vector<float_t, data_dimension> element_coordinates, gradient, boundary_normal, transformed_boundary_normal, temp;
      
      transform.boundary2element(parent_element_face, integrator.node(integrator_node), element_coordinates);

      transform.derivative(element_coordinates, derivative);
      inverse_derivative = inverse(derivative);
      transform.boundary_derivative(parent_element_face, integrator.node(integrator_node), boundary_transform_derivative);
      
      _boundary_transform_determinants(integrator_node) = transform_det(boundary_transform_derivative);

      for(std::size_t element_node = 0; element_node < n_element_nodes; ++element_node)
      {
        _shape_boundary_values(integrator_node)(element_node) =
          shape_function.value(element_node, element_coordinates);

        shape_function.gradient(element_node, element_coordinates, temp);
        gradient = prod(temp, inverse_derivative);

        _shape_boundary_derivatives(integrator_node)(element_node) = inner_prod(gradient, transformed_boundary_normal);
      }
    }
  }
  
  template <>
  float_t transform_det(const ublas::fixed_matrix<float_t, 1, 0> & derivative);

  template <>
  float_t transform_det(const ublas::fixed_matrix<float_t, 1, 1> & derivative);
   
  template <>
  float_t transform_det(const ublas::fixed_matrix<float_t, 2, 1> & derivative);
   
  template <>
  float_t transform_det(const ublas::fixed_matrix<float_t, 3, 2> & derivative);

}


#endif

---