This paper introduces a three-dimensional finite element (FE) formulation to accurately model the linear elastic deformation of fractured media under compressive loading. The presented method applies the classic Augmented Lagrangian(AL)-Uzawa method, to evaluate the growth of multiple interacting and intersecting discrete fractures. The volume and surfaces are discretized by unstructured quadratic triangle-tetrahedral meshes; quarter-point triangles and tetrahedra are placed around fracture tips. Frictional contact between crack faces for high contact precisions is modeled using isoparametric integration point-to-integration point contact discretization, and a gap-based augmentation procedure. Contact forces are updated by interpolating tractions over elements that are adjacent to fracture tips, and have boundaries that are excluded from the contact region. Stress intensity factors are computed numerically using the methods of displacement correlation and disk-shaped domain integral. A novel square-root singular variation of the penalty parameter near the crack front is proposed to accurately model the contact tractions near the crack front. Tractions and compressive stress intensity factors are validated against analytical solutions. Numerical examples of cubes containing one, two, twenty four and seventy interacting and intersecting fractures are presented.