00001
00002
00003
00004
00005
00006
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
