Grid.hpp

00001 // This file is part of the imaging2 class library.
00002 //
00003 // University of Innsbruck, Infmath Imaging, 2009.
00004 // http://infmath.uibk.ac.at
00005 //
00006 // All rights reserved.
00007 
00008 
00009 #ifndef GRID_H
00010 #define GRID_H
00011 
00012 #include <vector>
00013 #include <map>
00014 #include <set>
00015 
00016 #include <fem/FemKernel.hpp>
00017 #include <fem/fem_3d_tetrahedra_types.hpp>
00018 #include <fem/fem_2d_triangle_types.hpp>
00019 #include <fem/fem_1d_types.hpp>
00020 #include <shape/BoundaryDiscretizer.hpp>
00021 
00022 // #include <core/cio.hpp>
00023 
00024 
00025 namespace imaging
00026 {
00027   template<class fem_types>
00028   class Image2Grid;
00029   
00030   // TODO: clean up accessor functions: void set_vertex() instead of vertex_t & vertex()
00031   
00045   template<class fem_types>
00046   class Grid
00047   {
00048   private:
00049     class BoundaryElement
00050     {
00051     public:
00052       size_t _parent_element;
00053       size_t _parent_element_face;
00054       ublas::fixed_vector<float_t, fem_types::data_dimension> _normal;
00055     };
00056     
00057   public:
00059     static const size_t data_dimension = fem_types::data_dimension;
00060     
00062     static const size_t n_element_vertices = fem_types::transform_t::n_element_vertices;
00063     
00065     static const size_t n_element_nodes = fem_types::shape_function_t::n_element_nodes;
00066     
00068     static const size_t n_element_faces = fem_types::transform_t::n_element_faces;
00069     
00071     typedef ublas::fixed_vector<float_t, fem_types::data_dimension> vertex_t;
00072     
00074     typedef typename fem_types::transform_t transform_t;
00075     
00077     typedef typename fem_types::shape_function_t shape_function_t;
00078     
00080     typedef typename fem_types::integrator_t integrator_t;
00081     
00083     typedef typename fem_types::boundary_integrator_t boundary_integrator_t;
00084   
00085   private:
00086     typedef ublas::fixed_vector<size_t, n_element_vertices> element_vertices_t;
00087     typedef ublas::fixed_vector<size_t, n_element_nodes> element_nodes_t;
00088     
00089     std::vector<vertex_t> _vertices;
00090     std::vector<element_vertices_t> _element_vertices;
00091     std::vector<element_nodes_t> _element_nodes;
00092     std::vector<BoundaryElement> _boundary_elements;
00093     std::multimap<size_t, size_t> _boundary_node_boundary_elements;
00094     std::set<size_t> _boundary_nodes;
00095 
00096     size_t _n_nodes;
00097 
00098     bool _is_regular;
00099 
00100   public:
00101     
00103     Grid() : _n_nodes(0), _is_regular(false) {}
00104 
00106     size_t n_elements() const { return _element_vertices.size(); }
00107     
00109     size_t n_vertices() const { return _vertices.size(); }
00110     
00112     size_t n_nodes() const { return _n_nodes; }
00113     
00115     size_t n_boundary_nodes() const { return _boundary_nodes.size(); }
00116     
00118     size_t n_boundary_elements() const { return _boundary_elements.size(); }
00119     
00121     void set_vertex(size_t global_vertex_index, const vertex_t & vertex) { _vertices[global_vertex_index] = vertex; }
00122     
00124     const vertex_t & vertex(size_t global_vertex_index) const { return _vertices[global_vertex_index]; }
00125     
00127     void set_vertex(size_t element_index, size_t vertex_index, const vertex_t & vertex) { _vertices[global_vertex_index(element_index, vertex_index)] = vertex; }
00128     
00130     const vertex_t & vertex(size_t element_index, size_t vertex_index) const { return _vertices[global_vertex_index(element_index, vertex_index)]; }
00131     
00133     size_t global_vertex_index(size_t element_index, size_t vertex_index) const { return _element_vertices[element_index](vertex_index); }
00134     
00136     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; }
00137     
00139     size_t global_node_index(size_t element_index, size_t node_index) const { return _element_nodes[element_index](node_index); }
00140     
00142     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; }
00143     
00145     void set_boundary_element(size_t boundary_element_index, size_t parent_element_index, size_t parent_element_face)
00146     { 
00147       transform_t transform;
00148       ublas::fixed_vector<float_t, fem_types::data_dimension - 1> in;
00149       vertex_t out;
00150       element_transform(parent_element_index, transform);
00151       transform.boundary_normal(parent_element_face, out);
00152       
00153       _boundary_elements[boundary_element_index]._parent_element = parent_element_index; 
00154       _boundary_elements[boundary_element_index]._parent_element_face = parent_element_face;
00155       _boundary_elements[boundary_element_index]._normal = out;
00156       
00157       for(size_t i = 0; i < shape_function_t::n_shape_face_nodes; ++i)
00158       {
00159         _boundary_nodes.insert(global_node_index(parent_element_index, shape_function_t::face_node(parent_element_face, i)));
00160         _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));
00161       }
00162     }
00163     
00165     size_t parent_element(size_t boundary_element_index) const { return _boundary_elements[boundary_element_index]._parent_element; }
00166     
00168     size_t parent_element_face(size_t boundary_element_index) const { return _boundary_elements[boundary_element_index]._parent_element_face; }
00169     
00171     bool is_boundary_node(size_t global_node_index) const
00172     {
00173       return (_boundary_nodes.count(global_node_index) > 0);
00174     }
00175     
00177     bool is_boundary_node(size_t element_index, size_t node_index) const
00178     {
00179       return is_boundary_node(global_node_index(element_index, node_index));
00180     }
00181     
00183     const vertex_t boundary_normal(size_t global_node_index) const
00184     { 
00185       if( ! _boundary_node_boundary_elements.count(global_node_index))
00186         throw Exception("Invalid argument 'node' in Grid::boundary_normal().");
00187         
00188       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);
00189       
00190       std::multimap<size_t, size_t>::const_iterator iter = boundary_elements.first;
00191       
00192       vertex_t result(0.0);
00193       for( ; iter != boundary_elements.second; ++iter)
00194         result += _boundary_elements[iter->second]._normal;
00195       
00196       return result / norm_2(result);
00197     }
00198   
00200     const vertex_t boundary_normal(size_t element_index, size_t node_index) const
00201     { 
00202       return boundary_normal(global_node_index(element_index, node_index));
00203     }
00204 
00206     void set_regular(bool is_regular) { _is_regular = is_regular; }
00207     
00209     bool is_regular() const { return _is_regular; }
00210 
00212     void set_dimensions(size_t n_vertices, size_t n_elements, size_t n_boundary_elements, size_t n_nodes)
00213     {
00214       _vertices.resize(n_vertices);
00215       _element_vertices.resize(n_elements);
00216       _element_nodes.resize(n_elements);
00217       _boundary_elements.resize(n_boundary_elements);
00218 
00219       _n_nodes = n_nodes;
00220     }
00221 
00223     void element_transform(size_t element_index, transform_t & transform) const
00224     {
00225       for(size_t i = 0; i < n_element_vertices; ++i)
00226         transform.assign(i, _vertices[global_vertex_index(element_index, i)]);
00227     }
00228     
00232     void interpolate_value(size_t integrator_node, const ublas::vector<float_t> & data, const FemKernel<fem_types> & kernel, float_t & value) const
00233     { 
00234       value = 0.0;
00235       
00236       for(size_t i = 0; i < n_element_nodes; ++i)
00237         value += data(global_node_index(kernel.current_element(), i)) * kernel.shape_value(integrator_node, i);
00238     }
00239     
00243     template <size_t N>
00244     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
00245     { 
00246       value.assign(0.0);
00247       
00248       for(size_t i = 0; i < n_element_nodes; ++i)
00249         for(size_t j = 0; j < N; ++j)
00250           value(j) += data(global_node_index(kernel.current_element(), i) * N + j) * kernel.shape_value(integrator_node, i);
00251     }
00252     
00256     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
00257     { 
00258       gradient.assign(0.0);
00259       
00260       for(size_t i = 0; i < n_element_nodes; ++i)
00261         gradient += data(global_node_index(kernel.current_element(), i)) * kernel.shape_gradient(integrator_node, i);
00262     }
00263     
00264     
00268     template <size_t N>
00269     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
00270     { 
00271       gradient.assign(0.0);
00272       
00273       for(size_t i = 0; i < n_element_nodes; ++i)
00274         for(size_t j = 0; j < N; ++j)
00275           ublas::matrix_column< ublas::fixed_matrix<float_t, data_dimension, N> >(gradient, j) +=
00276             data(global_node_index(kernel.current_element(), i) * N + j) * kernel.shape_gradient(integrator_node, i);
00277     }
00278     
00282     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
00283     { 
00284       value = 0.0;
00285       size_t parent_element = _boundary_elements[kernel.current_boundary_element()]._parent_element;
00286 
00287       for(size_t i = 0; i < n_element_nodes; ++i)
00288         value += data(global_node_index(parent_element, i)) * kernel.shape_boundary_value(boundary_integrator_node, i);
00289     }
00290     
00294     template <size_t N>
00295     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
00296     { 
00297       value.assign(0.0);
00298       size_t parent_element = _boundary_elements[kernel.current_boundary_element()]._parent_element;
00299 
00300       for(size_t i = 0; i < n_element_nodes; ++i)
00301         for(size_t j = 0; j < N; ++j)
00302           value(j) += data(global_node_index(parent_element, i) * N + j) * kernel.shape_boundary_value(boundary_integrator_node, i);
00303     }
00304     
00308     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
00309     { 
00310       value = 0.0;
00311       size_t parent_element = _boundary_elements[kernel.current_boundary_element()]._parent_element;
00312 
00313       for(size_t i = 0; i < n_element_nodes; ++i)
00314         value += data(global_node_index(parent_element, i)) * kernel.shape_boundary_derivative(boundary_integrator_node, i);
00315     }
00316     
00320     template <size_t N>
00321     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
00322     { 
00323       value.assign(0.0);
00324       size_t parent_element = _boundary_elements[kernel.current_boundary_element()]._parent_element;
00325       
00326       for(size_t i = 0; i < n_element_nodes; ++i)
00327         for(size_t j = 0; j < N; ++j)
00328           value(j) += data(global_node_index(parent_element, i) * N + j) * kernel.shape_boundary_derivative(boundary_integrator_node, i);
00329     }
00330     
00332     void global_coordinates(size_t integrator_node, const FemKernel<fem_types> & kernel, ublas::fixed_vector<float_t, data_dimension> & coordinates) const
00333     { 
00334       transform_t transform;
00335       element_transform(kernel.current_element(), transform);
00336       
00337       transform.value( integrator_t().node(integrator_node), coordinates);
00338     }
00339     
00340 //     void print_grid() const
00341 //     {
00342 //       std::cout << _vertices << std::endl;
00343 //       std::cout << _element_vertices << std::endl;
00344 //       std::cout << _element_nodes << std::endl;
00345 //     }
00346   }
00347   ;
00348 }
00349 
00350 
00351 #endif

Generated on Tue Feb 10 10:01:30 2009 for imaging2 by  doxygen 1.5.5