A variant of the Regularized x-FEM (Rx-FEM) [1,2], where the Heaviside step function is replaced by a continuous or regularized function approximated by the finite element shape functions is further developed in this paper to achieve high crack density with respect to a given mesh. At the same time the method maintains the key feature of Rx-FEM where the elements are not partitioned but duplicated, which allows to maintain the gauss integration schema of the original element. Note that the nature of the regularized step function introduces a gradient zone over which the displacement enrichment is performed and thus generally requires a larger spacing between Rx-FEM cracks than the x-FEM cracks. In this paper the Rx-FEM framework is extended by increasing phantom node multiplicity to eliminate the gap between gradient regions of neighbouring cracks and significantly increase the crack density for a given mesh. In this case multiple phantom nodes are present, whereas the elements are only duplicated once. Dense matrix crack networks in laminated composites are often associated with fatigue loading, which is therefore chosen for the methodology illustration. Three different laminates containing open holes were considered under tension-tension fatigue loading. A single unstructured finite element mesh was used for discrete damage modelling in all laminates. Overall, the analysis results displayed excellent agreement with experimental data in terms of the size and location of damage including fine details of damage distribution such as stitch cracking of small cracks and delamination shapes. A further extension of the Rx-FEM formulation to multiple phantom node and multiple element twinning concept for modelling arbitrary 3D interacting crack networks is also considered. This framework allows for crack crossing and merger. The feature of the proposed method is that the enrichment does not require 3D cutting of an element by multiple cracks. It operates with separate regularized step functions of individual cracks and, as the original formulation, it maintains the Gauss integration schema in all element twins without regard to cracking pattern. Example of crack interaction modelling in 2D and 3D setting will be provided.