This study introduces a novel control oriented structure-preserving scheme for discretizing a class of multi-dimensional linear port-Hamiltonian systems, preserving their inherent structure while enabling the imposition of diverse combinations of boundary inputs, such as generalized velocities, displacements, and tractions. The proposed approach is grounded on the modified Linked Lagrange Multiplier method and the mixed Finite Element Method (FEM), where Dirichlet and Neumann boundary conditions are weakly enforced. Connections with other standard and mixed FEM approaches are also discussed. The proposed scheme is validated through comparisons with commercial software and simulations using a 2D elasticity model as a demonstrative example.