Simulations of biomembranes have gained an increasing interest in the past years. specificities of these membranes propose new challenges for the numerics. In particular, vesicle dynamics are governed by bending forces as well as a surface incompressibility constraint. A method to compute the bending force density resultant onto piecewise linearly triangulated surface meshes is described. This method is coupled with a boundary element method solver for inner and outer fluids, to compute vesicle dynamics under external flows. The surface incompressibility constraint is satisfied by the construction of a projection operator.