Fluid–structure interaction formulated in the Lagrangian–Eulerian approach with pressure-based fluid elements, which inherently results in non-symmetric matrices, is commonly treated by either direct un-symmetric solvers or staggered solution methods. Developing an in-house finite element code called SNACS, we implement the pressure-based Eulerian formulation so that a usual symmetric solver can be employed for dam–reservoir interaction analyses. The equilibrium equations of motion for such a coupled system are directly integrated herein by the HHT time marching algorithm in an incremental–iterative solution procedure to tackle a highly nonlinear dynamic analysis of jointed arch dams. The HHT time integration scheme is able to damp out high-frequency noises resulting from sudden changes in the stiffness due to opening/closing of pre-existing joints during an earthquake. First, a Pseudo-Symmetric technique is proposed to store the total matrices and solve the coupled non-symmetric equations in a symmetric manner. Then, the interface element formulation and stress update procedures of two discrete crack constitutive models suitable for contraction and peripheral joints of arch dams are detailed. Afterwards, results of the analysis on a typical arch dam are compared with the case in which the interaction is approximately represented by added masses as a widely-used method in dam engineering practice to illustrate the effects of water compressibility and reservoir bottom absorption. Besides, the straightforward implementation of pressure-based fluid–structure interaction utilizing a symmetric solution strategy is validated herein to be efficient for a 3-D large scale practical usage in other finite element codes even the commercial ones to avoid using un-symmetric solvers. Employing this symmetric implementation for the dam–water interaction drastically reduces the computational effort, particularly when an iterative scheme is required to solve such a highly nonlinear system.