In the following, we take solid mechanics as a demonstration example for setting up the weighted residual formulation. To do this, multiplying the PDE (3) on both sides by a weight function w and integrating it through the computational domain Ω , it follows that ∫ Ω w x ∂ ∂ x l D i j k l x ∂ u k x ∂ x j d Ω + ∫ Ω w x b i x d Ω = 0 (5)

Taking integration by parts and applying Gauss’ divergence theorem to the first domain integral of Eq. 5, the above equation becomes: ∫ Ω ∂ w x ∂ x l D i j k l x ∂ u k x ∂ x j d Ω = ∫ Γ w x t i x d Γ + ∫ Ω w x b i x d Ω (6) where t i is the traction component on the boundary Γ of the domain Ω , which has the relationship with the displacement gradient shown in Eq. 4.

In Eq. 6, the basic physical variable u k is mainly included in the volume integral of the left-hand side; therefore, it is called the volume-based weighted residual formulation. Taking integration by parts to the first domain integral of Eq. 6 and applying Gauss’ divergence theorem again, the following equation can be obtained: ∫ Γ ∂ w ∂ x l D i j k l n j u k d Γ − ∫ Ω ∂ ∂ x j D i j k l ∂ w ∂ x l u k d Ω = ∫ Γ w t i d Γ + ∫ Ω w b i d Ω (7)

In Eq. 7, the basic physical variable u k is included in both the surface and volume integrals of the left-hand side; therefore, it is called the surface-volume–based weighted residual formulation.

It is noted that Eqs 6, 7 are valid for any sized closed domain Ω , and from this feature various weak-form solution algorithms can be generated, such as FEMs and BEMs, by taking different kinds of the weigh function w in an element or in the whole domain.

In FEM, the computational domain is discretized into a series of elements [6, 7] with a certain number of nodes. Usually, the nodes on the element interfaces should be linked point-to-point, as shown in Figure 1, for a 2D computation FEM model. Over each element, the displacement u k is approximated using its nodal values u k α of the element by the shape function N α as follows: u k = N α u k α (8) where the repeated index α represents summation through all element nodes.

The above-described Galerkin FEM and SVFEM are weak-form algorithms, which require integration over elements to form the system of equations. In recent years, new types of strong-form FEM-like volume-operation–based methods have been proposed, which belong to a type of element collocation method and do not need integration computations. However, the stability of the strong-form algorithms is usually not as good as the weak-form algorithms, although for general problems these algorithms can still give satisfactory results.

Wen and Li et al. [8, 24] proposed the finite block method (FBM) in 2014, in which isoparametric element-like blocks are used to compute the first-order partial derivative of physical variables with respect to the global coordinates. FBM has the advantage of simple coding, and since element-like blocks are used, the stability of the solution is usually good. On the other hand, since all nodal values of physical variables over each block are independently inserted into the system of equations by introducing a consistent condition of physical variables and an equilibrium condition of the physical variable gradient in the system, there are more unknowns in the formed final system of equations than other frequently used numerical methods in the case of the same number of total nodes. In view of this issue, as few blocks as possible should be used when solving a problem using FBM to ensure that the final system of equations is not so large.

In the same period, Fantuzzi et al. [25, 26] proposed another type of strong-form FEM (SFEM), in which a set of formulations computing the first- and second-order spatial partial derivatives are derived for 2D problems and are used to collocate the governing PDEs in solid mechanics. In SFEM, the continuity condition among elements is determined by the compatibility, and a mapping technique is used to transform both the governing differential equations and the compatibility conditions between two adjacent sub-domains into the regular master element in the computational space. As in FBM, the treatment of the compatibility and equilibrium conditions between elements is still complicated; this makes SFEM not as flexible as GFEM when solving complicated engineering problems.

In 2017, Gao et al. [27] proposed a new type of strong-form FEM, called the Element Differential Method (EDM), for solving heat conduction problems, and later it was successfully used to solve solid mechanics [4, 28], electromagnetic [29], and thermo-mechanical-seepage coupled [30] problems. As in FBM, Lagrange polynomials are used to construct high-order elements in EDM. The essential difference between EDM and FBM is that both the first- and second-order partial derivatives were derived for 2D and 3D problems in EDM. The following is a brief review of EDM.

Looking back at Eq. 8 for physical variable interpolation over an element, the global coordinates can also be expressed by their nodal values and shape functions, as follows: x i = N α x i α (16)

Based on Eqs 8, 16, the following expressions can be derived for the first- and second-order partial derivatives [4, 27]: ∂ u k ∂ x i = d i c α ′ u k α ′ (17) ∂ 2 u k ∂ x i ∂ x j = d i j c α ″ u k α ″ (18) where d i c α ′ = ∂ N α ∂ x i = J i k − 1 ∂ N α ∂ ξ k (19a) d i j c α ″ = ∂ 2 N α ∂ x i ∂ x j = J i k − 1 ∂ 2 N α ∂ ξ k ∂ ξ l + ∂ J i k − 1 ∂ ξ l ∂ N α ∂ ξ k ∂ ξ l ∂ x j (19b) where J = ∂ x / ∂ ξ is the Jacobian matrix between the global coordinate x and the local coordinate ξ of the element. The detailed expressions for each term in Eqs 19a, 19b can be found in [27, 28].

The main advantage of the strong-form FEMs over weak-form FEMs is that the derived spatial partial derivatives can be directly substituted into the problem’s PEDs and B.C. to set up the system of equations. For example, by using Eqs 17, 18, the PDE and B.C. for the solid mechanics shown in Eqs 3, 4 can be directly used to generate the following discretized equations: d l c β ′ D i j k l x β ′ d j β ′ α ′ u k α ′ + b i x c = 0 , x c ∈ Ω e (20) u i x c = u ¯ i , x c ∈ Γ u D i j k l x c n j x c d l c α ′ u k α ′ = t ¯ i x c , x c ∈ Γ t (21)

In the above equations, Γ u and Γ t are the out boundaries of the problem, over which the displacement and traction boundary conditions are specified. The big issue in the strong-form FEM is how to set up the discretized equation when the collocation point x c is located on the interface Γ I between elements. To solve this issue, Gao et al. [4, 27] proposed the summation-equilibrium technique for all related tractions, that is, ∑ e = 1 , s = 1 t i e s = 0 , to form a single set of equations at an element interface node, which can be expressed as ∑ e = 1 , s = 1 E c D i j k l e x c n j e s x c d l e c α ′ u k e α ′ = 0 , x c ∈ Γ I (22) where E c is the number of all elements connected with the interface collocation point c. Various numerical examples [27–31] have proved that the above equation can give correct results. The important point is that Eq. 22 allows the final system of equations to have the same size as the conventional FEM, which is much smaller than those in FBM [8] and SFEM [25].

The FVM looks like a volume-based method [32–35]. However, in this article, it is classified into the category of surface-operation–based methods. This is because its main operation is over the surfaces of the control volume, rather than on the volume itself. To see this, let us take the weight function w to be 1 in Eq. 6. This results in the following: ∫ Γ t i d Γ + ∫ Ω b i d Ω = 0 (23)

In FVM, the computational domain is discretized into a series of control volumes [32]. Applying Eq. 23 to each control volume, say volume Ω c , and dividing its boundary ∂ Ω c into two parts, the inner boundary Γ I c and outer boundary Γ o c , Eq. 23 can be written as: ∫ Γ I c D i j k l n j ∂ u k ∂ x l d Γ + ∫ Γ o c t ¯ i d Γ + ∫ Ω c b i d Ω = 0 (24) where Γ I c ∪ Γ o c = ∂ Ω c .

Equation 24 is a typical formulation of FVM, from which we can see that the main computation is over the control surfaces of a control volume. The key work in FVM is the evaluation of the physical variable gradient ∂ u k / ∂ x l included in the first control surface integral of Eq. 24 [36, 37]. In the conventional FVM, the interface Γ I c is taken as the mid-surface connected by the collocation point c and around-neighbour points; thus, ∂ u k / ∂ x l at the mid-surface can be easily computed using the values of u k between c and the neighbour points [38–40]. However, only the linear variation of u k over the operation surface can be easily achieved. It is difficult to construct a high-order scheme to compute the value of ∂ u k / ∂ x l on the operation surface. To overcome this problem, the free element [15, 31] can be used in FVM analysis.

In Eq. 7, if the weight function w is taken as the displacement fundamental solution u m i * [9, 10], i.e., w = u m i * , the following BEM integral equation can be established: ∫ Γ ∂ u m i * ∂ x l D i j k l n j u k d Γ − ∫ Ω ∂ ∂ x j D i j k l ∂ u m i * ∂ x l u k d Ω = ∫ Γ u m i * t i d Γ + ∫ Ω u m i * b i d Ω (31)

Recalling that the displacement fundamental solution u m i * satisfies the following equation: ∂ ∂ x j D i j k l ∂ u m i * ∂ x l + δ m k δ p , q = 0 (32) Equation 31 becomes u m + ∫ Γ t m k * u k d Γ = ∫ Γ u m i * t i d Γ + ∫ Ω u m i * b i d Ω (33) where t m k * = D i j k l n j ∂ u m i * ∂ x l (34)

To evaluate the domain integral appearing in Eq. 33, the conventional technique is to discretize the domain into internal cells [9]; however, this eliminates the advantage of BEM where only the boundary of the problem needs to be discretized into elements. To overcome this drawback, a transformation technique is usually employed to transform the domain integral into an equivalent integral. The most extensively used transformation technique is the Dual Reciprocity Method [DRM] [43]. Another technique used is the Radial Integration Method (RIM) [42], which can give more accurate results than DRM.

Using RIM, the domain integral in Eq. 33 can be expressed as ∫ Ω u m i * b i d Ω = ∫ Γ F m i r n p , q ∂ r ∂ n   d Γ q (35) where the radial integral is F m i = ∫ 0 r p , q u m i * b i r n d r (36) where n = 1 for 2D problems and n = 2 for 3D problems. When b i is a known function, Eqs 35, 36 can give rise to a very accurate result. On the other hand, for the case where b i is not very complicated, Eq. 36 can be analytically integrated [42]. However, when b i is a complicated function, numerical integration should be performed [44].

In FrEM, a free element is independently formed for each collocation point c [31], with the domain denoted by Ω c . The shape function shown in Eq. 8 is still employed for the formed free element. Let us apply the weighted residual formulation (6) to Ω c and take the weight function as the shape function of the collocation point c, i.e., w = N c . Thus, the following equation can be obtained: ∫ Ω c ∂ N c ∂ x l D i j k l ∂ u k ∂ x j d Ω = ∫ ∂ Ω c N c t i x d Γ + ∫ Ω c N c b i d Ω (44) where the derivatives of the shape functions are computed using Eq. 19a.

Dividing the boundary ∂ Ω c of Ω c into two parts; the inner boundary Γ I c , which is located within the problem and the outer boundary Γ o c , which is located on the outer surface of the problem. Remembering that N c is zero on the surfaces excluding point c, making the integral over Γ I c zero, as a result, Eq. 44 can be written as: ∫ Ω c ∂ N c ∂ x l D i j k l ∂ N α ∂ x j d Ω   u k α = ∫ Γ o c N c t ¯ i d Γ + ∫ Ω c N c b i d Ω (45) where Γ o c = ∂ Ω c ∩ ∂ Ω , which is the outer boundary containing c.

From Eq. 45, it can be seen that the form of the basic equation in WFrEM is similar to that in the conventional FEM. The essential difference between them is that the element in WFrEM is freely formed at each collocation point, the nodes of which are not restricted to any particular nodes of adjacent elements. It is also noteworthy that the free elements formed by around-collocation points are overlapped in FrEM since they are formed locally and independently at each point.

The SFrEM is a type of collocation method [15]. To achieve a highly accurate result, the collocation point c should be placed inside the formed free element. For this reason, the free elements used should have at least one internal node. In principle, any type of isoparametric elements with internal nodes can be utilized in SFrEM analysis [54–56]. For example, Figure 8 shows new types of quadratic triangular and tetrahedral elements [55] and Figure 9 shows a 21-noded block element [15].

The shape functions for the above triangular and tetrahedral elements can be found in [55] and that for the 21-noded quadratic block element in [15].

For higher-order elements, the best method is to use Lagrange elements. For example, Figure 10 shows a 16-noded 2D third-order Lagrange element.

The shape functions of Lagrange elements for 2D and 3D problems can be constructed as follows [4]: N α ξ , η = L I ξ L J η , f o r   2 D N α ξ , η , ζ = L I ξ L J η L K ζ , f o r   3 D (46) where L I , L J , a n d   L K are determined by Eq. 39 and the superscript α is determined by the permutation of subscripts I, J, and K sequentially.

From the shape functions shown in Eq. 46, the analytical expressions for computing the first- and second-order partial derivatives can be derived, which are the same as those shown in Eqs 17–19a, 19b. The collocation scheme to form the system of equations for the governing PDEs is the same as that in EDM, shown in Eqs 20, 21, the difference being that in SFrEM, only the interior and out boundary nodes are used.