Finite Element Assembly Using an Embedded Domain Specific Language

Listing 6

Functor to perform numerical integration over an element and assemble the result into the global matrix.

(struct eval_element_quadrature: proto::callable

()  {

()   template<typename ExprT, typename DataT>

()   void operator()(Eigen::MatrixXd& mat, const ExprT& expr, DataT& data) const

()   {

()    // Temporary storage for the result from the RHS

()    Eigen::Matrix<double, DataT::element_t::nb_nodes, DataT::element_t::nb_nodes> rhs_result

      ;

()    rhs_result.setZero();

()  

()  // Loop over the gauss points

()  const typename DataT::element_t::gauss_points_t gauss_points = DataT::element_t::

       gauss_points();

()  const typename DataT::element_t::gauss_weights_t gauss_weights = DataT::element_t::

      gauss_weights();

()  for(int i = 0; i != DataT::element_t::nb_gauss_points; ++i)

()  {

()   data.shape_func = DataT::element_t::shape_function(gauss_points.row(i));

()   data.det_jacobian = DataT::element_t::jacobian_determinant(gauss_points.row(i), data.

       coord_mat);

()   rhs_result += gauss_weights[i] * data.det_jacobian * fem_grammar()(expr, 0, data);

()  }

()  

()  // Place the result in the global matrix

()  for(int i = 0; i != DataT::element_t::nb_nodes; ++i)

()   for(int j = 0; j != DataT::element_t::nb_nodes; ++j)

()    mat(data.node_indices[i], data.node_indices[j]) += rhs_result(i,j);

()   }

() };