Grid.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 GRID_H
#define GRID_H

#include <vector>
#include <map>
#include <set>

#include <fem/FemKernel.hpp>
#include <fem/fem_3d_tetrahedra_types.hpp>
#include <fem/fem_2d_triangle_types.hpp>
#include <fem/fem_1d_types.hpp>
#include <shape/BoundaryDiscretizer.hpp>

// #include <core/cio.hpp>


namespace imaging
{
  template<class fem_types>
  class Image2Grid;
  
  // TODO: clean up accessor functions: void set_vertex() instead of vertex_t & vertex()
  
  template<class fem_types>
  class Grid
  {
  private:
    class BoundaryElement
    {
    public:
      size_t _parent_element;
      size_t _parent_element_face;
      ublas::fixed_vector<float_t, fem_types::data_dimension> _normal;
    };
    
  public:
    static const size_t data_dimension = fem_types::data_dimension;
    
    static const size_t n_element_vertices = fem_types::transform_t::n_element_vertices;
    
    static const size_t n_element_nodes = fem_types::shape_function_t::n_element_nodes;
    
    static const size_t n_element_faces = fem_types::transform_t::n_element_faces;
    
    typedef ublas::fixed_vector<float_t, fem_types::data_dimension> vertex_t;
    
    typedef typename fem_types::transform_t transform_t;
    
    typedef typename fem_types::shape_function_t shape_function_t;
    
    typedef typename fem_types::integrator_t integrator_t;
    
    typedef typename fem_types::boundary_integrator_t boundary_integrator_t;
  
  private:
    typedef ublas::fixed_vector<size_t, n_element_vertices> element_vertices_t;
    typedef ublas::fixed_vector<size_t, n_element_nodes> element_nodes_t;
    
    std::vector<vertex_t> _vertices;
    std::vector<element_vertices_t> _element_vertices;
    std::vector<element_nodes_t> _element_nodes;
    std::vector<BoundaryElement> _boundary_elements;
    std::multimap<size_t, size_t> _boundary_node_boundary_elements;
    std::set<size_t> _boundary_nodes;

    size_t _n_nodes;

    bool _is_regular;

  public:
    
    Grid() : _n_nodes(0), _is_regular(false) {}

    size_t n_elements() const { return _element_vertices.size(); }
    
    size_t n_vertices() const { return _vertices.size(); }
    
    size_t n_nodes() const { return _n_nodes; }
    
    size_t n_boundary_nodes() const { return _boundary_nodes.size(); }
    
    size_t n_boundary_elements() const { return _boundary_elements.size(); }
    
    void set_vertex(size_t global_vertex_index, const vertex_t & vertex) { _vertices[global_vertex_index] = vertex; }
    
    const vertex_t & vertex(size_t global_vertex_index) const { return _vertices[global_vertex_index]; }
    
    void set_vertex(size_t element_index, size_t vertex_index, const vertex_t & vertex) { _vertices[global_vertex_index(element_index, vertex_index)] = vertex; }
    
    const vertex_t & vertex(size_t element_index, size_t vertex_index) const { return _vertices[global_vertex_index(element_index, vertex_index)]; }
    
    size_t global_vertex_index(size_t element_index, size_t vertex_index) const { return _element_vertices[element_index](vertex_index); }
    
    void set_global_vertex_index(size_t element_index, size_t vertex_index, size_t global_vertex_index) { _element_vertices[element_index](vertex_index) = global_vertex_index; }
    
    size_t global_node_index(size_t element_index, size_t node_index) const { return _element_nodes[element_index](node_index); }
    
    void set_global_node_index(size_t element_index, size_t node_index, size_t global_node_index) { _element_nodes[element_index](node_index) = global_node_index; }
    
    void set_boundary_element(size_t boundary_element_index, size_t parent_element_index, size_t parent_element_face)
    { 
      transform_t transform;
      ublas::fixed_vector<float_t, fem_types::data_dimension - 1> in;
      vertex_t out;
      element_transform(parent_element_index, transform);
      transform.boundary_normal(parent_element_face, out);
      
      _boundary_elements[boundary_element_index]._parent_element = parent_element_index; 
      _boundary_elements[boundary_element_index]._parent_element_face = parent_element_face;
      _boundary_elements[boundary_element_index]._normal = out;
      
      for(size_t i = 0; i < shape_function_t::n_shape_face_nodes; ++i)
      {
        _boundary_nodes.insert(global_node_index(parent_element_index, shape_function_t::face_node(parent_element_face, i)));
        _boundary_node_boundary_elements.insert(std::pair<size_t, size_t>(global_node_index(parent_element_index, shape_function_t::face_node(parent_element_face, i)), boundary_element_index));
      }
    }
    
    size_t parent_element(size_t boundary_element_index) const { return _boundary_elements[boundary_element_index]._parent_element; }
    
    size_t parent_element_face(size_t boundary_element_index) const { return _boundary_elements[boundary_element_index]._parent_element_face; }
    
    bool is_boundary_node(size_t global_node_index) const
    {
      return (_boundary_nodes.count(global_node_index) > 0);
    }
    
    bool is_boundary_node(size_t element_index, size_t node_index) const
    {
      return is_boundary_node(global_node_index(element_index, node_index));
    }
    
    const vertex_t boundary_normal(size_t global_node_index) const
    { 
      if( ! _boundary_node_boundary_elements.count(global_node_index))
        throw Exception("Invalid argument 'node' in Grid::boundary_normal().");
        
      std::pair< std::multimap<size_t, size_t>::const_iterator, std::multimap<size_t, size_t>::const_iterator > boundary_elements = _boundary_node_boundary_elements.equal_range(global_node_index);
      
      std::multimap<size_t, size_t>::const_iterator iter = boundary_elements.first;
      
      vertex_t result(0.0);
      for( ; iter != boundary_elements.second; ++iter)
        result += _boundary_elements[iter->second]._normal;
      
      return result / norm_2(result);
    }
  
    const vertex_t boundary_normal(size_t element_index, size_t node_index) const
    { 
      return boundary_normal(global_node_index(element_index, node_index));
    }

    void set_regular(bool is_regular) { _is_regular = is_regular; }
    
    bool is_regular() const { return _is_regular; }

    void set_dimensions(size_t n_vertices, size_t n_elements, size_t n_boundary_elements, size_t n_nodes)
    {
      _vertices.resize(n_vertices);
      _element_vertices.resize(n_elements);
      _element_nodes.resize(n_elements);
      _boundary_elements.resize(n_boundary_elements);

      _n_nodes = n_nodes;
    }

    void element_transform(size_t element_index, transform_t & transform) const
    {
      for(size_t i = 0; i < n_element_vertices; ++i)
        transform.assign(i, _vertices[global_vertex_index(element_index, i)]);
    }
    
    void interpolate_value(size_t integrator_node, const ublas::vector<float_t> & data, const FemKernel<fem_types> & kernel, float_t & value) const
    { 
      value = 0.0;
      
      for(size_t i = 0; i < n_element_nodes; ++i)
        value += data(global_node_index(kernel.current_element(), i)) * kernel.shape_value(integrator_node, i);
    }
    
    template <size_t N>
    void interpolate_value(size_t integrator_node, const ublas::vector<float_t> & data, const FemKernel<fem_types> & kernel, ublas::fixed_vector<float_t, N> & value) const
    { 
      value.assign(0.0);
      
      for(size_t i = 0; i < n_element_nodes; ++i)
        for(size_t j = 0; j < N; ++j)
          value(j) += data(global_node_index(kernel.current_element(), i) * N + j) * kernel.shape_value(integrator_node, i);
    }
    
    void interpolate_gradient(size_t integrator_node, const ublas::vector<float_t> & data, const FemKernel<fem_types> & kernel, ublas::fixed_vector<float_t, data_dimension> & gradient) const
    { 
      gradient.assign(0.0);
      
      for(size_t i = 0; i < n_element_nodes; ++i)
        gradient += data(global_node_index(kernel.current_element(), i)) * kernel.shape_gradient(integrator_node, i);
    }
    
    
    template <size_t N>
    void interpolate_gradient(size_t integrator_node, const ublas::vector<float_t> & data, const FemKernel<fem_types> & kernel, ublas::fixed_matrix<float_t, data_dimension, N> & gradient) const
    { 
      gradient.assign(0.0);
      
      for(size_t i = 0; i < n_element_nodes; ++i)
        for(size_t j = 0; j < N; ++j)
          ublas::matrix_column< ublas::fixed_matrix<float_t, data_dimension, N> >(gradient, j) +=
            data(global_node_index(kernel.current_element(), i) * N + j) * kernel.shape_gradient(integrator_node, i);
    }
    
    void interpolate_boundary_value(size_t boundary_integrator_node, const ublas::mapped_vector<float_t> & data, const FemKernel<fem_types> & kernel, float_t & value) const
    { 
      value = 0.0;
      size_t parent_element = _boundary_elements[kernel.current_boundary_element()]._parent_element;

      for(size_t i = 0; i < n_element_nodes; ++i)
        value += data(global_node_index(parent_element, i)) * kernel.shape_boundary_value(boundary_integrator_node, i);
    }
    
    template <size_t N>
    void interpolate_boundary_value(size_t boundary_integrator_node, const ublas::mapped_vector<float_t> & data, const FemKernel<fem_types> & kernel, ublas::fixed_vector<float_t, N> & value) const
    { 
      value.assign(0.0);
      size_t parent_element = _boundary_elements[kernel.current_boundary_element()]._parent_element;

      for(size_t i = 0; i < n_element_nodes; ++i)
        for(size_t j = 0; j < N; ++j)
          value(j) += data(global_node_index(parent_element, i) * N + j) * kernel.shape_boundary_value(boundary_integrator_node, i);
    }
    
    void interpolate_boundary_derivative(size_t boundary_integrator_node, const ublas::mapped_vector<float_t> & data, const FemKernel<fem_types> & kernel, float_t & value) const
    { 
      value = 0.0;
      size_t parent_element = _boundary_elements[kernel.current_boundary_element()]._parent_element;

      for(size_t i = 0; i < n_element_nodes; ++i)
        value += data(global_node_index(parent_element, i)) * kernel.shape_boundary_derivative(boundary_integrator_node, i);
    }
    
    template <size_t N>
    void interpolate_boundary_derivative(size_t boundary_integrator_node, const ublas::mapped_vector<float_t> & data, const FemKernel<fem_types> & kernel, ublas::fixed_vector<float_t, N> & value) const
    { 
      value.assign(0.0);
      size_t parent_element = _boundary_elements[kernel.current_boundary_element()]._parent_element;
      
      for(size_t i = 0; i < n_element_nodes; ++i)
        for(size_t j = 0; j < N; ++j)
          value(j) += data(global_node_index(parent_element, i) * N + j) * kernel.shape_boundary_derivative(boundary_integrator_node, i);
    }
    
    void global_coordinates(size_t integrator_node, const FemKernel<fem_types> & kernel, ublas::fixed_vector<float_t, data_dimension> & coordinates) const
    { 
      transform_t transform;
      element_transform(kernel.current_element(), transform);
      
      transform.value( integrator_t().node(integrator_node), coordinates);
    }
    
//     void print_grid() const
//     {
//       std::cout << _vertices << std::endl;
//       std::cout << _element_vertices << std::endl;
//       std::cout << _element_nodes << std::endl;
//     }
  }
  ;
}


#endif

---