(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.rhs and 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) )); |