Finite Element Assembly Using an Embedded Domain Specific Language

Listing 15

Code to build the linear system for in Chorin’s method.

(  // Allow a mix of first and second order shape functions

()  typedef boost::mpl::vector2<mesh::LagrangeP1::Triag2D, mesh::LagrangeP2::Triag2D>

     LagrangeP1P2;

()  // Kinematic viscosity, as a user-configurable constant:

()  PhysicsConstant nu("kinematic_viscosity");

()  // The velocity

()  FieldVariable<0, VectorField> u("u", "navier_stokes_u_velocity", "cf3.mesh.LagrangeP2");

()  // LSSActionUnsteady links with the linear algebra backend and the time tracking

()  Handle<LSSActionUnsteady> auxiliary_lss =

()   create_component<LSSActionUnsteady>("AuxiliaryLSS");

() // The matrix assembly

() Handle<ProtoAction> auxiliary_mat_assembly =

() auxiliary_lss->create_component<ProtoAction>("MatrixAssembly");

() auxiliary_mat_assembly->set_expression(elements_expression(LagrangeP1P2(),

()  group

()  (

() _A(u) = _0, _T(u) = _0,

() element_quadrature

() (

()   _A(u[_i], u[_i]) += transpose(nabla(u))*nabla(u),

()  _T(u[_i], u[_i]) += transpose(N(u))*N(u)

() ),

() auxiliary_lss->system_matrix += auxiliary_lss->invdt()*_T + nu*_A

() )));

()  

() // RHS assembly

() Handle<ProtoAction> auxiliary_rhs_assembly =

() auxiliary_lss->create_component<ProtoAction>("RHSAssembly");

() auxiliary_rhs_assembly->set_expression(elements_expression(LagrangeP1P2(),

() group

() (

() _a[u] = _0,

() element_quadrature

() (

()  _a[u[_i]] += auxiliary_lss->invdt() * transpose(N(u))*u[_i] −

()           transpose(N(u))*(u*nabla(u))*_col(nodal_values(u), _i)

() ),

() auxiliary_lss->system_rhs += _a

() )));