From a variety of yield stress fluid models, an elastoviscoplastic Herschel-Bulkley (EVPHB) model written in 3D is selected and coupled with a Finite Element Method with Lagrangian Integration Points (FEMLIP) to solve boundary value problems with large deformation process. By tracing the historical variables of a material point, it is verified that in a time-independent flow the elastic strain and viscous strain rate could be accurately reproduced by EVPHB model. For a time-dependent flow, because of the addition of elasticity, the EVPHB model makes the material experience a deformation process which is significantly distinctive from that produced by a pure regularized Herschel-Bulkley model. Benchmarks also show that in FEMLIP the yielded and unyielded zones could be easily defined by EVPHB model according to the stress of a material point. Lastly, it is shown that EVPHB model also induces a stress relaxation process for materials under constant strain. The suitability of FEMLIP to model elastoviscoplastic fluid is verified.