SimpleEquationAdaptor.hpp

00001 #ifndef EQUATION_SIMPLEEQUATIONADAPTOR_H
00002 #define EQUATION_SIMPLEEQUATIONADAPTOR_H
00003 
00004 
00005 #include <fem/equation/EquationInterface.hpp>
00006 #include <fem/equation/SimpleEquationInterface.hpp>
00007 #include <fem/FemKernel.hpp>
00008 
00009 
00010 namespace imaging
00011 {
00034   template <class equation_t>
00035   class SimpleEquationAdaptor : public EquationInterface<typename equation_t::fem_types>
00036   {
00037     const equation_t & _equation;
00038     typedef typename equation_t::fem_types fem_types;
00039     
00040   public:
00041     const std::size_t system_size() const { return 1; }
00042     
00044     SimpleEquationAdaptor(const equation_t & equation) : _equation(equation) {}
00045     
00046     float_t stiffness_matrix(std::size_t equation, std::size_t component,
00047                              std::size_t i, std::size_t j,
00048                              std::size_t integrator_node,
00049                              const FemKernel<fem_types> & kernel) const
00050     {
00051       float_t value = 0.0;
00052       
00053       typename equation_t::matrix_coefficient_t A;
00054       ublas::fixed_vector<float_t, fem_types::data_dimension> a, b;
00055       float_t c;
00056       
00057       if( ! kernel.is_boundary_node(i) ||
00058             _equation.boundary_data_type == SimpleEquationInterface<fem_types>::NO_BOUNDARY_DATA ||
00059             _equation.boundary_data_type == SimpleEquationInterface<fem_types>::IMPLICIT_NEUMANN_DATA )
00060       {
00061         _equation.stiffness_matrix(integrator_node, kernel, A, a, b, c);
00062         
00063         value += inner_prod( prod(A, kernel.shape_gradient(integrator_node, i)),
00064                                 kernel.shape_gradient(integrator_node, j) );
00065 
00066         if(_equation.a_active)
00067           value += inner_prod(a, kernel.shape_gradient(integrator_node, j) ) *
00068                   kernel.shape_value(integrator_node, i);
00069 
00070         if(_equation.b_active)
00071           value += inner_prod(b, kernel.shape_gradient(integrator_node, i) ) *
00072                   kernel.shape_value(integrator_node, j);
00073 
00074         if(_equation.c_active)
00075           value += c *
00076                    kernel.shape_value(integrator_node, i) *
00077                    kernel.shape_value(integrator_node, j);
00078       }
00079       
00080       return value;
00081     }
00082     
00083     
00084     float_t force_vector(std::size_t equation,
00085                          std::size_t i,
00086                          std::size_t integrator_node,
00087                          const FemKernel<fem_types> & kernel) const
00088     {
00089       float_t value = 0.0;
00090       
00091       float_t f;
00092       ublas::fixed_vector<float_t, fem_types::data_dimension> g;
00093       
00094       if( ! kernel.is_boundary_node(i) ||
00095             _equation.boundary_data_type == equation_t::NO_BOUNDARY_DATA ||
00096             _equation.boundary_data_type == equation_t::IMPLICIT_NEUMANN_DATA )
00097       {
00098         _equation.force_vector(integrator_node, kernel, f, g);
00099         
00100         if(_equation.f_active)
00101           value += f *
00102                    kernel.shape_value(integrator_node, i);
00103 
00104         if(_equation.g_active)
00105           value += inner_prod(g, kernel.shape_gradient(integrator_node, i) );
00106       }
00107                 
00108       return value;
00109     }
00110     
00111     
00112     float_t stiffness_matrix_at_boundary(std::size_t equation, std::size_t component,
00113                                          std::size_t i, std::size_t j,
00114                                          std::size_t integrator_node,
00115                                          const FemKernel<fem_types> & kernel) const
00116     {
00117       float_t value = 0.0;
00118       
00119       float_t h;
00120       
00121       if(_equation.boundary_data_type != equation_t::NO_BOUNDARY_DATA)
00122       {
00123         if(_equation.boundary_data_type == equation_t::MIXED_DATA)
00124         {
00125           _equation.stiffness_matrix_at_boundary(integrator_node, kernel, h);
00126           
00127           value += kernel.shape_boundary_derivative(integrator_node, i) * kernel.shape_boundary_value(integrator_node, j) *
00128                    kernel.boundary_transform_determinant(integrator_node);
00129 
00130           value += h *
00131                    kernel.shape_boundary_value(integrator_node, i) * kernel.shape_boundary_value(integrator_node, j);
00132         }
00133 
00134         if(_equation.boundary_data_type == equation_t::NEUMANN_DATA)
00135           value += kernel.shape_boundary_derivative(integrator_node, i) * kernel.shape_boundary_value(integrator_node, j);
00136     
00137         if(_equation.boundary_data_type == equation_t::DIRICHLET_DATA)
00138           value += kernel.shape_boundary_value(integrator_node, i) * kernel.shape_boundary_value(integrator_node, j);
00139       }
00140                 
00141       return value;
00142     }
00143     
00144     
00145     float_t force_vector_at_boundary(std::size_t equation,
00146                                      std::size_t i,
00147                                      std::size_t integrator_node,
00148                                      const FemKernel<fem_types> & kernel) const
00149     {
00150       float_t value = 0.0;
00151       
00152       float_t v;
00153       
00154       if(_equation.boundary_data_type != equation_t::NO_BOUNDARY_DATA)
00155       {
00156         _equation.force_vector_at_boundary(integrator_node, kernel, v);
00157         
00158         if(_equation.boundary_data_type == equation_t::IMPLICIT_NEUMANN_DATA)
00159           value += -1.0 *
00160                    v *
00161                    kernel.shape_boundary_value(integrator_node, i);
00162         else
00163           value += v *
00164                    kernel.shape_boundary_value(integrator_node, i);
00165       }      
00166                 
00167       return value;
00168     }
00169     
00170     bool sanity_check_stiffness_matrix(const FemKernel<fem_types> & kernel, std::string & error_message) const
00171     { 
00172       return _equation.sanity_check_stiffness_matrix(kernel, error_message);
00173     }
00174     
00175     bool sanity_check_force_vector(const FemKernel<fem_types> & kernel, std::string & error_message) const
00176     { 
00177       return _equation.sanity_check_force_vector(kernel, error_message);
00178     }
00179   };
00180 }
00181 
00182 #endif
00183 

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