Assembler.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 FEM_ASSEMBLER_H
00010 #define FEM_ASSEMBLER_H
00011 
00012 #include <fem/Grid.hpp>
00013 #include <fem/FemKernel.hpp>
00014 
00015 
00016 
00017 
00018 namespace imaging
00019 {
00025   class Assembler
00026   {
00027     static void clear_matrix(ublas::compressed_matrix<float_t> & matrix);
00028     
00029   public:
00030 
00039     template<class fem_types, class equation_t>
00040     void assemble(const equation_t & equation, const Grid<fem_types> & grid,
00041                   ublas::compressed_matrix<float_t> & stiffness_matrix,
00042                   ublas::vector<float_t> & force_vector) const;
00043     
00052     template<class fem_types, class equation_t>
00053     void assemble_stiffness_matrix(const equation_t & equation, const Grid<fem_types> & grid,
00054                   ublas::compressed_matrix<float_t> & stiffness_matrix) const;
00055                   
00064     template<class fem_types, class equation_t>
00065     void assemble_force_vector(const equation_t & equation, const Grid<fem_types> & grid,
00066                   ublas::vector<float_t> & force_vector) const;
00067 
00068   }
00069   ;
00070 
00071   template<class fem_types, class equation_t>
00072   void Assembler::assemble(const equation_t & equation, const Grid<fem_types> & grid,
00073                                        ublas::compressed_matrix<float_t> & stiffness_matrix,
00074                                        ublas::vector<float_t> & force_vector) const
00075   {
00076     typedef typename fem_types::shape_function_t shape_function_t;
00077     typedef typename fem_types::integrator_t integrator_t;
00078     typedef typename fem_types::boundary_integrator_t boundary_integrator_t;
00079     typedef typename fem_types::transform_t transform_t;
00080     
00081     if(stiffness_matrix.size1() != equation.system_size() * grid.n_nodes() &&
00082        stiffness_matrix.size2() != equation.system_size() * grid.n_nodes() )
00083       throw Exception("Exception: Dimension of stiffness matrix does not agree with grid size in Assembler::assemble()");
00084       
00085     clear_matrix(stiffness_matrix);
00086 
00087     force_vector.resize(equation.system_size() * grid.n_nodes(), false);
00088     force_vector.clear();
00089 
00090     integrator_t integrator;
00091     boundary_integrator_t boundary_integrator;
00092 
00093     FemKernel<fem_types> kernel(grid);
00094 
00095     if(grid.is_regular() && grid.n_elements() > 0)
00096       kernel.set_element(0);
00097       
00098     std::string sanity_check_message = "";
00099     if( ! equation.sanity_check_stiffness_matrix(kernel, sanity_check_message) )
00100       throw Exception("Exception: sanity check failed in Assembler::assemble() with message '" + sanity_check_message + "'.");
00101     if( ! equation.sanity_check_force_vector(kernel, sanity_check_message) )
00102       throw Exception("Exception: sanity check failed in Assembler::assemble() with message '" + sanity_check_message + "'.");
00103       
00104     for(size_t element = 0; element < grid.n_elements(); ++element)
00105     {
00106       float_t value;
00107 
00108       if(grid.is_regular())
00109         kernel.lazy_set_element(element);
00110       else
00111         kernel.set_element(element);
00112 
00113       for(size_t i = 0; i < fem_types::shape_function_t::n_element_nodes; ++i)
00114       {
00115         for(size_t l = 0; l < equation.system_size(); ++l)
00116         {
00117           for(size_t j = 0; j < fem_types::shape_function_t::n_element_nodes; ++j)
00118           {
00119             for(size_t m = 0; m < equation.system_size(); ++m)
00120             {
00121               value = 0.0;
00122   
00123               for(size_t k = 0; k < integrator_t::n_nodes; ++k)
00124                 value += equation.stiffness_matrix(l, m, i, j, k, kernel) * 
00125                          kernel.transform_determinant(k) *
00126                          integrator.weight(k);
00127   
00128               stiffness_matrix(equation.system_size() * grid.global_node_index(element, i) + l,
00129                                equation.system_size() * grid.global_node_index(element, j) + m) += value;
00130             }
00131           }
00132 
00133           value = 0.0;
00134 
00135           for(size_t k = 0; k < integrator_t::n_nodes; ++k)
00136             value += equation.force_vector(l, i, k, kernel) *
00137                          kernel.transform_determinant(k) *
00138                          integrator.weight(k);
00139 
00140           force_vector(equation.system_size() * grid.global_node_index(element, i) + l) += value;
00141         }
00142       }
00143     }
00144 
00145 
00146     for(size_t element = 0; element < grid.n_boundary_elements(); ++element)
00147     {
00148       kernel.set_boundary_element(element);
00149       size_t parent_element = grid.parent_element(element);
00150       
00151       for(size_t i = 0; i < fem_types::shape_function_t::n_element_nodes; ++i)
00152       {
00153         for(size_t l = 0; l < equation.system_size(); ++l)
00154         {  
00155           for(size_t j = 0; j < fem_types::shape_function_t::n_element_nodes; ++j)
00156           {
00157             for(size_t m = 0; m < equation.system_size(); ++m)
00158             {
00159               float_t value = 0.0;
00160   
00161               for(size_t k = 0; k < boundary_integrator_t::n_nodes; ++k)
00162                 value += equation.stiffness_matrix_at_boundary(l, m, i, j, k, kernel) *
00163                  kernel.boundary_transform_determinant(k) *
00164                  boundary_integrator.weight(k);
00165                 
00166               stiffness_matrix(equation.system_size() * grid.global_node_index(parent_element, i) + l,
00167                                equation.system_size() * grid.global_node_index(parent_element, j) + m) += value;
00168             }
00169 
00170           }
00171 
00172           float_t value = 0.0;
00173 
00174           for(size_t k = 0; k < boundary_integrator_t::n_nodes; ++k)
00175             value += equation.force_vector_at_boundary(l, i, k, kernel) *
00176                  kernel.boundary_transform_determinant(k) *
00177                  boundary_integrator.weight(k);
00178 
00179           force_vector(equation.system_size() * grid.global_node_index(parent_element, i) + l) += value;
00180         }
00181       }
00182     }
00183 
00184   }
00185 
00186   template<class fem_types, class equation_t>
00187   void Assembler::assemble_stiffness_matrix(const equation_t & equation,
00188                                       const Grid<fem_types> & grid,
00189                                       ublas::compressed_matrix<float_t> & stiffness_matrix) const
00190   {
00191     typedef typename fem_types::shape_function_t shape_function_t;
00192     typedef typename fem_types::integrator_t integrator_t;
00193     typedef typename fem_types::boundary_integrator_t boundary_integrator_t;
00194     typedef typename fem_types::transform_t transform_t;
00195    
00196     if(stiffness_matrix.size1() != equation.system_size() * grid.n_nodes() &&
00197        stiffness_matrix.size2() != equation.system_size() * grid.n_nodes() )
00198       throw Exception("Exception: Dimension of stiffness matrix does not agree with grid size in Assembler::assemble()");
00199 
00200     clear_matrix(stiffness_matrix);
00201     
00202     integrator_t integrator;
00203     boundary_integrator_t boundary_integrator;
00204 
00205     FemKernel<fem_types> kernel(grid);
00206 
00207     if(grid.is_regular() && grid.n_elements() > 0)
00208       kernel.set_element(0);
00209       
00210     std::string sanity_check_message = "";
00211     if( ! equation.sanity_check_stiffness_matrix(kernel, sanity_check_message) )
00212       throw Exception("Exception: sanity check failed in Assembler::assemble() with message '" + sanity_check_message + "'.");  
00213       
00214     for(size_t element = 0; element < grid.n_elements(); ++element)
00215     {
00216       float_t value;
00217 
00218       if(grid.is_regular())
00219         kernel.lazy_set_element(element);
00220       else
00221         kernel.set_element(element);
00222 
00223       for(size_t i = 0; i < fem_types::shape_function_t::n_element_nodes; ++i)
00224       {
00225         for(size_t l = 0; l < equation.system_size(); ++l)
00226         {
00227           for(size_t j = 0; j < fem_types::shape_function_t::n_element_nodes; ++j)
00228           {
00229             for(size_t m = 0; m < equation.system_size(); ++m)
00230             {
00231               value = 0.0;
00232   
00233               for(size_t k = 0; k < integrator_t::n_nodes; ++k)
00234                 value += equation.stiffness_matrix(l, m, i, j, k, kernel) * 
00235                          kernel.transform_determinant(k) *
00236                          integrator.weight(k);
00237   
00238               stiffness_matrix(equation.system_size() * grid.global_node_index(element, i) + l,
00239                                equation.system_size() * grid.global_node_index(element, j) + m) += value;
00240             }
00241           }
00242         }
00243       }
00244     }
00245 
00246 
00247     for(size_t element = 0; element < grid.n_boundary_elements(); ++element)
00248     {
00249       kernel.set_boundary_element(element);
00250       size_t parent_element = grid.parent_element(element);
00251       
00252       for(size_t i = 0; i < fem_types::shape_function_t::n_element_nodes; ++i)
00253       {
00254         for(size_t l = 0; l < equation.system_size(); ++l)
00255         {  
00256           for(size_t j = 0; j < fem_types::shape_function_t::n_element_nodes; ++j)
00257           {
00258             for(size_t m = 0; m < equation.system_size(); ++m)
00259             {
00260               float_t value = 0.0;
00261   
00262               for(size_t k = 0; k < boundary_integrator_t::n_nodes; ++k)
00263                 value += equation.stiffness_matrix_at_boundary(l, m, i, j, k, kernel) *
00264                  kernel.boundary_transform_determinant(k) *
00265                  boundary_integrator.weight(k);
00266                 
00267               stiffness_matrix(equation.system_size() * grid.global_node_index(parent_element, i) + l,
00268                                equation.system_size() * grid.global_node_index(parent_element, j) + m) += value;
00269             }
00270           }
00271         }
00272       }
00273     }
00274 
00275   }
00276 
00277   template<class fem_types, class equation_t>
00278   void Assembler::assemble_force_vector(const equation_t & equation,
00279                                       const Grid<fem_types> & grid,
00280                                       ublas::vector<float_t> & force_vector) const
00281   {
00282     typedef typename fem_types::shape_function_t shape_function_t;
00283     typedef typename fem_types::integrator_t integrator_t;
00284     typedef typename fem_types::boundary_integrator_t boundary_integrator_t;
00285     typedef typename fem_types::transform_t transform_t;
00286    
00287     force_vector.resize(equation.system_size() * grid.n_nodes(), false);
00288     force_vector.clear();
00289 
00290     integrator_t integrator;
00291     boundary_integrator_t boundary_integrator;
00292 
00293     FemKernel<fem_types> kernel(grid);
00294 
00295     if(grid.is_regular() && grid.n_elements() > 0)
00296       kernel.set_element(0);
00297       
00298     std::string sanity_check_message = "";
00299     if( ! equation.sanity_check_force_vector(kernel, sanity_check_message) )
00300       throw Exception("Exception: sanity check failed in Assembler::assemble() with message '" + sanity_check_message + "'.");
00301       
00302     for(size_t element = 0; element < grid.n_elements(); ++element)
00303     {
00304       float_t value;
00305 
00306       if(grid.is_regular())
00307         kernel.lazy_set_element(element);
00308       else
00309         kernel.set_element(element);
00310 
00311       for(size_t i = 0; i < fem_types::shape_function_t::n_element_nodes; ++i)
00312       {
00313         for(size_t l = 0; l < equation.system_size(); ++l)
00314         {
00315           value = 0.0;
00316 
00317           for(size_t k = 0; k < integrator_t::n_nodes; ++k)
00318             value += equation.force_vector(l, i, k, kernel) *
00319                          kernel.transform_determinant(k) *
00320                          integrator.weight(k);
00321 
00322           force_vector(equation.system_size() * grid.global_node_index(element, i) + l) += value;
00323         }
00324       }
00325     }
00326 
00327 
00328     for(size_t element = 0; element < grid.n_boundary_elements(); ++element)
00329     {
00330       kernel.set_boundary_element(element);
00331       size_t parent_element = grid.parent_element(element);
00332       
00333       for(size_t i = 0; i < fem_types::shape_function_t::n_element_nodes; ++i)
00334       {
00335         for(size_t l = 0; l < equation.system_size(); ++l)
00336         {  
00337           float_t value = 0.0;
00338 
00339           for(size_t k = 0; k < boundary_integrator_t::n_nodes; ++k)
00340             value += equation.force_vector_at_boundary(l, i, k, kernel) *
00341                  kernel.boundary_transform_determinant(k) *
00342                  boundary_integrator.weight(k);
00343 
00344           force_vector(equation.system_size() * grid.global_node_index(parent_element, i) + l) += value;
00345         }
00346       }
00347     }
00348 
00349   }
00350 }
00351 
00352 
00353 #endif

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