Research Article

Finite Element Assembly Using an Embedded Domain Specific Language

Listing 19

Code for the specialized user-defined terminal for the Poisson problem, valid for linear shape functions over a triangle.
(1) // Specialized code for triangles
(2) struct PoissonTriagAssembly
(3) {
(4)  typedef void result_type;
(5)  // Functor that takes: source term f, Linear system lss, element matrix and vector acc
(6)  template  <typename FT, typename LSST>
(7)  void operator()(const FT& f, LSST& lss, math::LSS::BlockAccumulator& acc) const
(8)  {
(9)   typedef mesh::LagrangeP1::Triag2D ElementT;
(10)   // Get the coordinates of the element nodes
(11)   const ElementT::NodesT& nodes = f.support().nodes();
(12)
(13)   // Compute normals
(14)   ElementT::NodesT normals;
(15)   normals(0, XX) = nodes(1, YY) nodes(2, YY);
(16)   // repetitive code omitted
(17)   // Jacobian determinant
(18)   const Real det_jac = normals(2, YY)*normals(1, XX) normals(1, YY)*normals(2, XX);
(19)   const Real c = 1. / (2.*det_jac);
(20)
(21)   // Indices of the nodes of the current element
(22)   acc.neighbour_indices(f.support().element_connectivity());
(23)
(24)   for(Uint i = 0; i != 3; ++i)
(25)    for(Uint j = 0; j != 3; ++j)
(26)     acc.mat(i, j) = c * (normals(i, XX)*normals(j, XX) + normals(i, YY)*normals(j, YY));
(27)
(28)   // Get the values of the source term
(29)   const Real f0 = f.value()[];
(30)   // f1 and f2
(31)   acc.rhs[] = (2*f0 + f1 + f2);
(32)   // acc.rhsand acc.rhs
(33)   acc.rhs *= det_jac/24.;
(34)   lss.matrix().add_values(acc);
(35)   lss.rhs().add_rhs_values(acc);
(36)  }
(37) };
(38) // Create an terminal that can be used as a function in a proto expression
(39) static MakeSFOp  <PoissonTriagAssembly>::type const assemble_triags = ;
(40) // Usage example:
(41) assembly->set_expression(elements_expression
(42) (
(43)  boost::mpl::vector1< mesh::LagrangeP1::Triag2D>(),
(44)  assemble_triags(f, system_matrix, m_block_accumulator)
(45) ));