###### Abstract
Typically, the analysis of heterogeneous materials involves the determination of effective properties obtained through homogenising the response at microscopic scales, which are subsequently extrapolated to macroscopic analyses. Numerous analytical models have been developed within the framework of small deformation linear elasticity to obtain the homogenised constitutive response of heterogeneous materials at the macroscopic level, including the attributes of their microstructure [10]. These models are founded on the Hill-Mandel condition for homogeneity [11; 12], which posits that the volume-averaged strain energy within a representative volume element (RVE) can be expressed as the product of the volume-averaged stress and strain fields within the same RVE, thereby demonstrating energy equivalence between homogeneous and heterogeneous materials [11; 12]. Micromechanical analysis of composite materials frequently employs discretisation-based approaches to extrapolate the overall response from their microstructure [13]. Addressing this, the Generalised Method of Cells (GMC) offers a solution by defining a repeating unit cell (RUC) within the periodic composite structure and further subdividing it into orthogonal subcells [14]. As a mature numerical technique, the GMC framework has been effectively applied across various composite types, including fibre-reinforced composites [14], metal matrix composites [15], and woven polymer matrix composites [16]. Within computational mechanics, the finite element method (FEM) is the most commonly used discretisation technique, possessing a long history of successful applications across various engineering domains [17]. Nonetheless, the computational mechanics research community continuously investigates and explores advanced new discretisation methodologies capable of offering enhanced efficiency and accuracy. Meshless methods, also known as meshfree methods, emerge as promising alternatives capable of supplanting FEM in numerous applications [18; 19]. Unlike FEM, which relies on a standard element mesh for domain discretisation, meshless methods employ an unstructured nodal set to discretise the solid domain [20; 21]. While FEM establishes nodal connectivity through the element concept, meshless methods achieve this connectivity via the influence domain concept [22]. The Diffuse Element Method (DEM), introduced by [PERSON] et al. [23], emerged as a pioneering meshless technique. It applies moving least square (MLS) approximants, offering a generalisation of FEM but without the constraints of a pre-defined mesh [23]. [PERSON] et al. further developed DEM, refining it and extending its applicability to elasticity problems [24]. Their key innovation involved the incorporation of a background integration mesh based on the Gauss-Legendre quadrature scheme. [PERSON] et al. also coupled DEM with Lagrange multipliers for enforcing boundary conditions, forming the well-known Element Free Galerkin Method (EFGM) [24], the most popular meshless method ever created. The late 1990s witnessed a surge in meshless method development. The meshless method's portfolio was enriched with robust formulations, such as the Reproducing Kernel Particle Method (RKPM) [25] or the Meshless Local Petrov-Galerkin (MLPG) method [25]. These methods, by exploiting the advantages of approximation shape functions and higher nodal connectivity, were able to deliver more accurate solutions and smoother variable fields [22]. Despite their advantages, meshless methods based on approximation shape functions face a fundamental impairment: the absence of the delta Kronecker property. Unlike the FEM, this characteristic impedes the straightforward imposition of essential (displacement) and natural (traction) boundary conditions [22]. Thus, in approximation meshless methods, Lagrange multipliers are the primary numerical tool for handling boundary conditions.