Finite Element Assembly Using an Embedded Domain Specific Language

Listing 8

The element looping algorithm.

(  /// MPL functor to loop over elements

()  template<typename ExprT>

()  struct element_looper

()  {

()   element_looper(const fem::mesh_data& m, const ExprT& e): mesh(m), expr(e)

()   {

()   }

()  

()   template < typename ElemT >

() void operator()(ElemT) const

() {

()  // Bail out if the shape function doesn't match

()  if(ElemT::dimension != mesh.coordinates.shape()[] || ElemT::nb_nodes != mesh.

      connectivity.shape()[])

()     return;

()  

()  // Construct helper data

()  dsl_data<ElemT> data(mesh);

()  

()  // Evaluate for all elements

()  const int nb_elems = mesh.connectivity.size();

()  for(int i = 0; i != nb_elems; ++i)

()  {

()     data.set_element(i);

()     fem_grammar()(expr, 0, data);

()  }

() }

()  

() const fem::mesh_data& mesh;

() const ExprT& expr;

() };

()  

() /// Execute the given expression for every element in the mesh

() template<typename ExprT>

() void for_each_element(const fem::mesh_data& mesh, const ExprT& expr)

() {

() // Allowed element types

() typedef mpl::vector2<fem::line1d, fem::triag2d> element_types;

()  mpl::for_each<element_types>(element_looper<ExprT>(mesh, expr));

() }