Accurate simulations of cracks remain a challenging task in the computational mechanics community. The main purpose of this work is to develop a nodal-based Lagrange multiplier/cohesive zone (LM/CZ) method for arbitrary three-dimensional (3D) crack simulations in the finite element (FE) models of quasi-brittle materials. The main idea of the developed method is to fulfill the displacement compatibility condition before crack onset via LMs enforced at nodes and describe crack behaviors with the aid of the cohesive zone model (CZM). The main appeal of our method lies in effectively addressing the so-called artificial compliance numerical issue arising in the conventional intrinsic CZMs, thereby yielding more accurate numerical results than intrinsic CZMs. Node-pair groups are generated via a robust data structure. In each group, the LM forces for all nodal pairs are calculated in a predictor–corrector form, which can achieve convergence within a limited number of iterations using the Gauss–Seidel method. Two loop types for nodal pairs are compared to find a more efficient scheme with fewer iterations. The nodal constraints are smoothly switched from LMs to cohesive forces after a defined fracture criterion is satisfied. The developed method is implemented in conjunction with tetrahedral elements. Several quasi-static numerical examples, including modes I and I/II cracks, are performed to demonstrate the numerical accuracy of our developed method. Finally, this work simulates a Kolthoff plate impact problem, where the computational accuracy and efficiency are further investigated and discussed by comparing our results with those calculated via the intrinsic CZMs.