1. Introduction
===============

The recent progress in computers has enabled a number of complicated natural phenomena to be accurately simulated. Among various computer simulation techniques, the finite element method (FEM) has been most widely used since the invention of the method in late 1950s in a variety of practical fields, such as mechanical, aerospace, nuclear, chemical and civil engineering. It can be applied flexibly to complex boundary shapes using unstructured computational grids, and results, accurate enough for engineering purposes, are obtainable at reasonable cost.^[@r01],[@r02])^

When using FEM, one needs to divide a continuum into a number of elements. The individual elements are connected together by a topological map, which is usually called a mesh. The finite element interpolation functions are then built on the mesh. Not only in FEM,^[@r01],[@r02])^ but in the Finite Difference Method (FDM),^[@r03])^ and the Finite Volume Method (FVM),^[@r04])^ the problem spatial domain is discretized into meshes. The mesh must be predefined to provide a certain relationship between the nodes, which becomes the building blocks of the formulation procedure of these conventional numerical methods.

It is well known that the mesh generation process often becomes a serious bottleneck in large-scale computing,^[@r05]--[@r07])^ especially for problems that require frequent mesh refinement, such as the moving boundary problems. The limitations of the mesh-based methods, particularly those of FEM are summarized as follows,

1.  The analyst spends much time in creating the mesh, which becomes a major part of the cost of a simulation process. As the cost of central processing unit (CPU) time is drastically decreasing, the concern is more the manpower time, and less the computer time.

2.  It is rather difficult to simulate both crack growth with arbitrary and complex paths and phase transformations due to discontinuities that do not coincide with the original nodal lines.

3.  Remeshing methods have been proposed for handling these types of problems, where the problem domain is remeshed at steps during the simulation process to prevent the severe distortion of meshes and to allow the nodal lines to remain coincident with the discontinuity boundaries. For this purpose, complex, robust, and adaptive mesh generation processors have to be developed. However, these processors are only workable for 2D problems.

It is natural that the concept of meshless method has evolved,^[@r08]--[@r12])^ which is the idea of eliminating or reducing the reliance on the elements and more flexible ways to make use of mesh. In this concept, the domain of the problem is represented, ideally, only by a set of arbitrarily distributed nodes, where the physical field under consideration is approximated by a set of points. Smoothed particle hydrodynamics (SPH),^[@r13]--[@r16])^ diffuse element method (DEM),^[@r17])^ gridless Euler/Navier--Stokes solution,^[@r18],[@r19])^ element-free Galerkin method (EFGM),^[@r20]--[@r22])^ reproducing kernel particle method (RKPM),^[@r23])^ hp-meshless cloud method,^[@r24],[@r25])^ finite point method (FPM),^[@r26])^ moving-particle semi-implicit method (MPS),^[@r27])^ meshless local boundary integral equation (MLBIE)^[@r28])^ and meshless local Petrov--Galerkin method (MLPG)^[@r29])^ are considered as typical meshless methods.

Among these methods, SPH was originally developed for modeling astronomical phenomena using particles, with RKPM and MPS being improvements to this development. These particle type meshless methods require a predefinition of particles for their volumes or masses. The algorithm will then carry out the analyses even if the problem domain undergoes extremely large deformation or separation. This type of method suffers from problems in the imposition of boundary conditions. In addition, predefining the particles still technically requires some kind of mesh. SPH simulates well the overall behaviors of certain class of problems such as highly nonlinear and momentum-driven problems. These particle-based methods mainly aim at numerical simulation of phenomena such as splitting, scattering and uniting of bodies.

In general, methods that do not require a mesh at all are less stable and less accurate. Local point collocation methods and FDMs using irregular grids are among this category. Selection of nodes based on the type of a physical problem can be important for obtaining stable and accurate results. Automation of nodal selection and improving the stability of the solution are still some of the challenges in these kinds of methods. This type of method has a very significant advantage: It is very easy to implement, because no integration is required. There are, however, vital instability issues that require special treatments.^[@r30],[@r31])^

In contrast, DEM, EFGM, MLBIE and MLPG are classified as meshless methods in a narrow sense, and the major aim of these methods is the numerical simulation of continuum mechanics without finite element mesh, while the nature of the approximation is very similar to FEM in that they are based on the Galerkin method. While these meshless methods exhibit excellent performance in some special fields, they have not always succeeded in replacing the conventional FEM. The primary reason is that the meshless methods are not effective in representing the three-dimensional complex shapes required in practical engineering problems.

In this context, a number of FEM-based meshless approaches have been proposed to overcome the difficulty of mesh generation. Manifold method,^[@r32])^ voxel finite element method,^[@r33])^ generalized finite element method (GFEM),^[@r34])^ partition of unity finite element method (PUFEM),^[@r35])^ extended finite element method (X-FEM)^[@r36],[@r37])^ and finite cover method (FCM)^[@r38])^ can be classified as FEM-based meshless method. Even though these methods employ elements or a mesh to approximate the physical field, they are still regarded as meshless methods, because an explicit mesh is not required in the input data, or else a flexible mesh, that can be prepared much more easily than a conventional FEM mesh, is available as input data.

The free mesh method (FMM)^[@r39]--[@r45])^ is one of the earliest FEM-based meshless methods intended for particle-like finite element analysis of problems that are difficult to handle using global mesh generation, such as moving boundary problems, large deformation problems, crack propagation analysis, the separation or unification of bodies, as well as adaptive mesh refinement analysis, especially on parallel processors. FMM is a node-based finite element method that employs a local mesh generation technique and a node-by-node algorithm for the finite element calculations. The method can overcome difficulties arising from the distortion of elements by employing a node-based approach, with maintaining the superior theoretical foundations and the ability to handle complicated boundary shapes in accordance with finite elements. In FMM, both pre-processing and main-processing of analysis can easily be parallelized in terms of nodes, where the pre-processing involves local mesh generation, and the main-processing the construction and solution of a system of equations. The method is quite suitable for parallel environments.

Reviewed in this paper are FMM as well as the Enriched Free Mesh Method (EFMM):^[@r46]--[@r53])^ a new version of FMM, focusing on their fundamental conception, algorithms and accuracy. In the following section, the fundamental concept of FMM is reviewed, and the third section deals with EFMM. Finally, the present paper is summarized by some concluding remarks.

2. Basic conception of Free Mesh Method (FMM)^[@r39],[@r48])^
=============================================================

2.1. Outline.
-------------

FMM is a node-based finite element method, characterized by both the node-based local mesh generation and the node-by-node finite element calculation, which starts with distributing the nodes across the analysis domain. Data to input to the method are surface patches and nodes as shown in Fig. [1](#fig01){ref-type="fig"}, where surface patches are sets of segments in two-dimensional space or planes in three-dimensional space, which surround the analysis domain and define the shape of the domain.

Figure [2](#fig02){ref-type="fig"} illustrates the basic concept of the method. First, a single node *P*~*i*~ among the distributed nodes is designated as a central node. A local mesh is then generated around *P*~*i*~ using any local mesh generation technique. A group of elements around *P*~*i*~, which consist of several local meshes, are called the satellite elements associated with *P*~*i*~. Nodes of satellite elements, other than the central node, are referred to as satellite nodes. Local meshes must be generated under the constraints that the edges of the satellite elements do not cross each other, and that the analysis domain is covered fully by local meshes. If this condition fails, the satellite elements are judged to be inconsistent (see Fig. [3](#fig03){ref-type="fig"}). After the above local mesh generation, the components of the global coefficient matrix associated with *P*~*i*~ are calculated by integrating all the satellite elements around *P*~*i*~. The above process is repeated for all the nodes in the domain, and a global system of equations is obtained after all the nodes have been integrated. The flowchart of the method is shown in Fig. [4](#fig04){ref-type="fig"}. Several techniques have been proposed for the above local mesh generation.

2.2. Algorithms for local mesh generation.
------------------------------------------

We summarizes three types of local mesh generation algorithms in what follows.

### Algorithm I: segment comparison algorithm.^[@r39],[@r45])^

The segment comparison algorithm is proposed in the early stage of FMM. Satellite nodes for each central node are found by comparing the lengths of segments formed by nominated nodes existing in a local area called the searching circle, which is defined around the central node with a radius determined by the density of nodes at the central node. The process of searching for satellite nodes is performed using local information about nodal co-ordinates, and can be parallelized in a node-based manner.

### Algorithm II: local Delaunay triangulation algorithm with searching circle.^[@r52])^

A disadvantage of the segment comparison algorithm described above is related to the determination of an appropriate radius of the searching circle. If the radius is too large or too small, the algorithm may fail to determine appropriate satellite nodes, so that inconsistent overlapping of satellite elements occurs as shown in Fig. [3](#fig03){ref-type="fig"}. This problem is referred to as inconsistency in satellite elements.^[@r45])^ The presence of only a few inconsistencies in satellite elements is known to have little influence on the solution accuracy in the case of linear and stationary problems.^[@r46])^ However, too many inconsistencies may cause unsatisfactory results in non-linear or time-dependent problems, such as fluid dynamics problems. Therefore, a local Delaunay triangulation algorithm is developed as a substitute for the segment comparison algorithm, where the Delaunay triangulation^[@r54]--[@r57])^ is performed locally around each node. However, if a normal Delaunay triangulation procedure is used to search for satellite nodes, inconsistency problems can still occur when degeneracy exists in the Delaunay triangulation (see Fig. [5](#fig05){ref-type="fig"}), because satellite nodes are selected independently at each node in a node-by-node algorithm. To overcome this problem, the Delaunay triangulation algorithm has been modified using nodal identification number information. This algorithm has the advantage that inconsistencies in satellite elements can be avoided in almost any distribution of nodes, as long as the radius of the searching circle is sufficiently large.

### Algorithm III: local Delaunay triangulation algorithm with gift-wrapping method.^[@r47],[@r58])^

The local Delaunay triangulation algorithm mentioned above employs a searching circle. This sometimes causes a critical failure in the algorithm, because no theoretical basis exists for determining the radius of the searching circle. If the radius is too small, the algorithm cannot provide a proper number of satellite elements. On the other hand, if the radius is too large, too many unnecessary tentative elements are generated, and the algorithm needs extra computer power. In particular, this is a serious problem if the density of nodes varies rapidly over a domain to be solved. To overcome this problem, alternative local Delaunay triangulation algorithm is developed based on the gift-wrapping method,^[@r59],[@r60])^ which is an algorithm that calculates a triangle of Delaunay triangulation (DT) in a one-by-one manner. The algorithm is well suited to parallel computing, because the algorithm is highly localized spatially. Local meshes can always be generated by minimum searching of nodes without searching circles. Furthermore, the algorithm can generate local meshes robustly, even when the object contains knife-edge geometries.

Consider a triangle *abc* in a two-dimensional space (see Fig. [6](#fig06){ref-type="fig"}). To obtain an adjacent Delaunay triangle on the left side of segment *bc* (see Fig. [6](#fig06){ref-type="fig"}(a)), the circumcircle of the triangle *abc* is first drawn as shown in Fig. [6](#fig06){ref-type="fig"}(b). The circle is then moved as shown in Fig. [6](#fig06){ref-type="fig"}(c), until the circle meets a node (see Fig. [6](#fig06){ref-type="fig"}(d)). When the circle meets node *d*, a Delaunay triangle is formed between nodes *b*, *c* and *d* because the circle is nothing but the largest empty circle at this moment. Briefly, a new Delaunay triangle is obtained by using segment *bc* of the settled Delaunay triangle and the nearest node in the outside of the segment *bc*. It is noted that this Delaunay triangle is determined uniquely.

Furthermore, the gift-wrapping method for obtaining Delaunay triangle can easily be extended to the three-dimensional space as shown in Fig. [7](#fig07){ref-type="fig"}. It is noted here that Figs. [7](#fig07){ref-type="fig"}(a) through [7](#fig07){ref-type="fig"}(d) correspond to Figs. [6](#fig06){ref-type="fig"}(a) through [6](#fig06){ref-type="fig"}(d), respectively. However, the Delaunay triangulation described above examines only node co-ordinates data, so that certain boundary shapes, such as reentrant angles, cannot be treated. In order to deal with these shapes, the constrained Delaunay triangulation (CDT),^[@r61]--[@r63])^ which divides the given domain into Delaunay-like triangulations using the conforming constraint of surface patches, is proposed: a new algorithm for obtaining constrained Delaunay triangulations based on the gift-wrapping in a node-based manner^[@r47],[@r58])^ (see Fig. [8](#fig08){ref-type="fig"}). Another difficulty in the three-dimensional mesh generation is related to the so-called slivers, which can be avoided by changing the connectivity of nodes (see Fig. [9](#fig09){ref-type="fig"}). However, this operation reduces the efficiency of parallel mesh generation considerably. The remedy to avoid slivers will be to regularly arrange nodes as much as possible, which can be done by generating input nodal data using the Centroidal Voronoi Tessellation (CVT).^[@r64],[@r65])^

2.3. Node-by-node algorithm and parallelization based on nodes.^[@r43])^
------------------------------------------------------------------------

The present method consists of a node-by-node procedure for local mesh generation and construction of the global matrix. In other words, the global matrix is constructed in a row-by-row manner (see Figs. [2](#fig02){ref-type="fig"} and [4](#fig04){ref-type="fig"}), where the non-zero components of the global matrix in the row of a node correspond to the components derived from the satellite nodes associated with the node.

The final system of equations obtained by the present node-by-node procedure is equivalent to that obtained by the element-by-element procedure of the conventional FEM, while difference between the two procedures appears in the parallel computing. Figure [10](#fig10){ref-type="fig"}(a) shows the element-based domain decomposition for parallel computing, where the global matrix is not constructed explicitly, but parallel computing is performed instead in an element-by-element manner.^[@r66],[@r67])^ On the other hand, the node-based domain decomposition is employed in the proposed method (see Fig. [10](#fig10){ref-type="fig"}(b)), where nodes in the analysis domain are classified into two types: communication-independent nodes and communication-dependent nodes. The former ones are defined as the nodes, which do not require communication between processors when parallel computing is performed, whereas the latter ones are defined as those which require communication between processors. Figure [11](#fig11){ref-type="fig"} shows both types of nodes regarding PE1 in Fig. [10](#fig10){ref-type="fig"}. Parallel efficiency can be improved by dividing the nodes into these two categories.

The node-based partition of computational loads across multiprocessors shown in Fig. [10](#fig10){ref-type="fig"}(b) relates to row-based division for parallel calculation of the product of a matrix and a vector (see Fig. [12](#fig12){ref-type="fig"}). Communication between processors is required when satellite nodes associated with a central node are located on other processors than that in which the central node is located (see Fig. [13](#fig13){ref-type="fig"}), where the satellite nodes *S*~1~ and *S*~6~ are handled by a processor different from that handles the central node *P*~*i*~. In other words, inter-processor communication is required for nodes *S*~1~ and *S*~6~ when the product of a matrix by a vector is calculated. These satellite nodes are referred to as external satellite nodes, whereas other satellite nodes are called as internal satellite nodes. In Fig. [12](#fig12){ref-type="fig"}, the internal satellite nodes correspond to those marked by open circles, whereas the external satellite nodes correspond to those marked by gray triangles and the central nodes to the diagonal components of the matrix. Inter-processor communication is carried out regarding the components marked by gray triangles in Fig. [12](#fig12){ref-type="fig"}. The amount of communication between processors depends on the numbering of nodal identifiers. In this regards, a parallel graph-partitioning library ParMETIS^[@r68])^ is effectively used to renumber nodal identification numbers to minimize the amount of communication.^[@r47])^

2.4. Results of local mesh generation.
--------------------------------------

In Fig. [14](#fig14){ref-type="fig"}, the meshes generated by the above algorithms are compared, where Figs. [14](#fig14){ref-type="fig"}(a)--(c) show the results of Algorithms I, II and III, respectively. Bold lines in Fig. [14](#fig14){ref-type="fig"}(a) depict the inconsistencies in satellite elements. As can be seen from these figures, satisfactory results are achieved except Algorithm I. As another comparison of these three algorithms, meshes with combination of two different nodal densities are shown in Fig. [15](#fig15){ref-type="fig"}, where Algorithm I is seen to fail again to generate appropriate local meshes if the radius of the searching circle is too large or too small. In Algorithm II, the inconsistencies can be avoided if the radius of the searching circle is sufficiently large. However, using too large a searching circle in Algorithm II results in many useless elements, being computationally wasteful. On the other hand, Algorithm III can generate local meshes robustly by a minimal search of neighboring nodes as demonstrated in Figs. [14](#fig14){ref-type="fig"}(c) and [15](#fig15){ref-type="fig"}(c), respectively. Furthermore, Algorithm III is sufficiently robust to treat singular boundary shapes such as cracks: Fig. [16](#fig16){ref-type="fig"}(a) shows input surface patch and nodal data as an example of a knife-edge boundary shape and Fig. [16](#fig16){ref-type="fig"}(b) the mesh generated by Algorithm III.

Figures [17](#fig17){ref-type="fig"} (a)--(f) show the process of mesh generation applied to a three-dimensional crack model using Algorithm III, in which central nodes are selected at random and all the local meshes are generated independently around each node in a robust manner.

Figures [18](#fig18){ref-type="fig"} (a)--(f) and Figs. [19](#fig19){ref-type="fig"} (a)--(f) show the processes of Algorithm III applied to the three dimensional mesh generation for a model of the "Pantheon" building, which are viewed from the outside (see Fig. [18](#fig18){ref-type="fig"}) and from the inside (see Fig. [19](#fig19){ref-type="fig"}), respectively. Figures [20](#fig20){ref-type="fig"} (a)--(f) show the process of parallel mesh generation with eight processors. As can be seen in these figures, the three-dimensional meshes are robustly generated. It is noted that the present mesh generation technique is favored with parallel environments: Fig. [21](#fig21){ref-type="fig"} shows a large-scale mesh generated by Algorithm III, where the mesh shown in Fig. [21](#fig21){ref-type="fig"}(a) has approximately 1.2 × 10^5^ nodes and 4.7 × 10^5^ elements and that shown in Fig. [21](#fig21){ref-type="fig"}(b) approximately 2.1 × 10^7^ nodes and 1.2 × 10^8^ elements, in which virtually parallel technique is employed with single processor.^[@r58])^

3. Enriched Free Mesh Method (EFMM)^[@r50],[@r51])^
===================================================

3.1. Outline.
-------------

"Assumed strain on the clustered local elements" characterizes EFMM (see Fig. [22](#fig22){ref-type="fig"}), where the strain field on the clustered local elements and the displacement field of each local element are assumed independently. Coupling these two independent fields, the following approaches have been proposed: the localized least square method and the method based on the Hellinger--Reissner variational principle.

3.2. EFMM based on localized least square method.
-------------------------------------------------

EFMM based on the localized least square method (EFMM-LS) assumes the strain field on the clustered local elements aswhere {ε(***x***)}={ε~*xx*~,ε~*yy*~,γ~*xy*~} is the strain field defined on the clustered local elements and each component of {ε~*xx*~,ε~*yy*~,γ~*xy*~} is assumed independently, and \[***N***^ε^\] is a matrix, which consists of arbitrary polynomials as follows,where *p*^*t*^(**x**) is assumed on the clustered local elements asThe coefficients vector {**a**} in Eq. \[[1](#e01){ref-type="disp-formula"}\] is obtained by minimizing the *L*~2~ norm *J* as follows,where *n*~*e*~ is the number of local elements with *c* (=1,2,...*n*~*e*~) being current local element, *p* the number of points, which are called as the "strain monitoring points" on the clustered local elements with *i* (=1,2,...*p*) being the current strain monitoring point and {ε~*i*~^*c*^} the strain vector of *i*-th strain monitoring point on the *c*-th local element, which is called as the "mother element". The stationary condition of Eq. \[[4](#e04){ref-type="disp-formula"}\] iswhich yields the coefficients vector {**a**} as follows, Let us consider a simple constant strain triangle as the mother element, in which the displacement field of each local element is defined bywhere {**u**} is the displacement field of the local element, {**u**~*i*~} the nodal displacement, and ζ~*i*~ the area-coordinate.^[@r69])^ Thus, the strain on the strain monitoring points is given bywhereBy substituting Eq. \[[8](#e08){ref-type="disp-formula"}\] into Eq. \[[6](#e06){ref-type="disp-formula"}\], the unknown coefficient {**a**} is determined asSubstituting Eq. \[[10](#e10){ref-type="disp-formula"}\] into Eq. \[[1](#e01){ref-type="disp-formula"}\], we obtainwhere In the elasticity problem, the stress vector {σ} and the strain vector {ε} have the relation as follows,where \[**D**\] is the symmetric matrix of material stiffness. With using Eq. \[[12](#e12){ref-type="disp-formula"}\], the stiffness matrix based on the localized least square method is calculated on the clustered local elements aswhere Ω is the clustered local elements. It is important to note that the above stiffness matrix is computed in a node-wise manner.

It is also noted that EFMM-LS is closely related to the superconvergent patch recovery proposed by Zienkiewicz and Zhu (Z--Z).^[@r70],[@r71])^ In the adaptive finite element method,^[@r72],[@r73])^ the Z--Z error estimator has been most widely used to estimate the error. The error estimator requires an exact solution, but generally it is impossible to have it. The Z--Z technique then obtains the recovered solution in a post processing stage. The clustered local elements in the present method are equivalent to the superconvergent patch used in the Z--Z technique. The difference between these appears in the fact that the recovering procedure in EFMM-LS is in a main process stage when computing element stiffness matrices. The use of the assumed strain is, in some sense, equivalent to the "post-process" of the Z--Z superconvergent patch recovery.

3.3. EFMM based on Hellinger--Reissner principle.
-------------------------------------------------

The Hellinger--Reissner (HR) variational principle is employed in EFMM based on the Hellinger--Reissner principle (EFMM-HR)^[@r74])^ to improve the accuracy of solutions. Let the HR principle of a linear elastic body be defined on the clustered local elements bywherewith {**b**} being the applied body force per unit mass, and ![](pjab-87-115-e026.jpg) the applied traction on the boundary *S*~σ~, {**ū**} the unknown nodal displacement and {![](pjab-87-115-e027.jpg)} the unknown nodal strain. The unknown values (**ū**, ![](pjab-87-115-e027.jpg)) of the HR principle satisfy the following equations in a weak manner, and

It is noted here that the strain field is defined on the clustered local elements by a node-wise manner, whereas the displacement field is defined on each element by an element-wise manner. Equations \[[17](#e17){ref-type="disp-formula"}\] and \[[18](#e18){ref-type="disp-formula"}\] yield the linear matrix equation as follows,whereBy condensing the coefficient matrix of Eq. \[[19](#e19){ref-type="disp-formula"}\], we obtain the following equation:where the condensation should be performed on the clustered local elements. Thus, the stiffness matrix based on the HR principle is computed on the clustered local elements as follows,It is noted that we can obtain the enriched stiffness matrix without increasing the degrees of freedom of nodes.

3.4. Convergence study.
-----------------------

Two kinds of error norms for a beam bending problem as shown in Fig. [23](#fig23){ref-type="fig"} ^[@r75])^ are taken, which are, respectively, given as andwhere \|*E*\|~2~ is the displacement error norm and \|*E*\|~*e*~ the energy error norm, respectively. **u**, ε and σ are, respectively, the numerical results of displacement, strain and stress, whereas ***u***^*exact*^, ε^*exact*^ and σ^*exact*^ the exact solutions. Assume that the beam of length L = 10, height D = 2 and thickness t = 1 is subjected to a shear load in plane stress. The material parameters are given as Young's modulus *E* = 1000 and Poisson's ratio *v* = 0.25. The above displacement and energy convergence norms are plotted against the DOFs in Figs. [24](#fig24){ref-type="fig"} and [25](#fig25){ref-type="fig"}, respectively, where the meshes employed are 1 × 1 (see Fig. [23](#fig23){ref-type="fig"}), 2 × 2, 4 × 4, 8 × 8, 16 × 16, 32 × 32 and 64 × 64, respectively. It can be seen from these figures that

1.  The error norms of displacement of EFMMs are between those of the linear and the quadratic FEMs (see Fig. [24](#fig24){ref-type="fig"}). However, the convergence rates of EFMMs are almost equal to that of the quadratic FEM.

2.  The error norms of energy of the quadratic EFMMs are almost the same as that of the quadratic FEM, whereas those of the linear EFMMs are between the linear and the quadratic FEMs.

The patch test is finally performed using the three models of patch as shown in Fig. [26](#fig26){ref-type="fig"}, where the displacement field:is applied at the boundary. Table [1](#tbl01){ref-type="table"} shows the results for FEMs and EFMMs. As shown in the table, all of the methods pass the patch test for the Model A, which is a regular mesh division model. However, for the Models B and C, which are irregular ones, EFMM-LSs fail to pass the test. Here, "Pass" means that the displacement of the internal node (M1 or M2) satisfies Eq. \[[25](#e25){ref-type="disp-formula"}\]. This implies that EFMM-LSs are of non-conforming type for irregular mesh.

Summarized here are that the accuracy of EFMM are equivalent to that of FEM (quadratic), which is considered to be the most accurate, although the CPU time and the required memory size are found much less compared with those of the latter.

4. Concluding remarks
=====================

In the present paper, the fundamental conception, the algorithms and the accuracy of the Free Mesh Method (FMM), are reviewed, which aims at a particle-like finite element analysis of problems that are difficult to handle with using global mesh generation, such as moving boundaries problems, large deformation problems, crack propagation analysis, the separation or unification of bodies and adaptive mesh refinement analysis, especially on parallel processors. In the method, local finite elements are generated around each node, with local mesh data structures and a system of equations being created based on these nodes. As such, FMM can overcome difficulties associated with the distortion of elements by simply adding or deleting nodes, similar to particle methods. This characteristic is particularly advantageous when a parallel computer is employed. While parallelization of mesh generation is generally difficult, only the distribution of nodes needs to be considered in this method to perform parallel remeshing that favors excellent load balancing between processors. This node-based finite element computation is realized by a robust local mesh generation technique based on the gift-wrapping method. However, the study that still remains is the fast and robust generation of nodes, especially on massively parallel processors. A probabilistic node generation approach with centroidal Voronoi tessellation under a parallel environment will be the most promising.

Next, a new Free Mesh Method called Enriched Free Mesh Method (EFMM) is taken, by which a high accuracy can be obtained without explicitly increasing the degrees of freedom. The central idea of EFMM is that the strain field is assumed on clustered local elements in addition to the usual FEM displacement field on each element. To relate the above two fields, the localized least square method or the Hellinger--Reissner principle are independently employed. The convergence characteristics of the displacement error norms are between that of FEM with the linear displacement field and that with the quadratic one, whereas that of the energy error norms with the quadratic strain field for the clustered elements is equivalent to that of FEM with the quadratic displacement field. EFMM based on the Hellinger--Reissner principle passes the patch test, whereas that based on the localized least square method does not for irregular nodal arrangements. This would be an open question and there is a room for future research.

I would like to thank all the researchers and students, who have been involved in the study of FMM with the present author, which has been performed under the sponsorship of University of Tokyo, Toyo University and Sungkyunkwan University in Korea. I am grateful to Professor Tamio Yamakawa, Editor-in-Chief for Proceedings of the Japan Academy, Series B, for his encouragement to prepare this review paper.

Genki Yagawa was born in 1942 and started his research career in 1965 with studies on the finite element method superposed on the analytical solution, applying to linear and nonlinear fracture mechanics, after graduating from the Faculty of Engineering at the University of Tokyo. Among the most important works he has accomplished is a massively parallel computing method, where the hierarchical domain decomposition together with the variational method to satisfy the boundary conditions between the sub-domains is employed. It has been proved that the method can solve huge scale problems of billions-degrees of freedom with high accuracy and high performance when applied to computational solid and fluid mechanics problems. He has also created a new method called the Free Mesh Method, where both pre-processing and main-processing of the finite element analysis can be parallelized easily in terms of nodes. These methods, being quite suitable for massively parallel environments and embedded within the internationally well-known codes, have been employed by many researchers and engineers worldwide.

He became Lecturer in 1970, Associate Professor in 1971, Professor in 1984 and Council Member in 1998 at the University of Tokyo. He is currently Emeritus Professor at the University of Tokyo, Professor and Director at the Computational Mechanics Research Center of Toyo University, Visiting Professor at Sungkyunkwan University in Korea, Council Member of Science Council of Japan, and President of International Association for Computational Mechanics.

Among his awards and honors are the Japan Academy Prize (2009), the Toray Science and Technology Prize (2009), the International Association for Computational Mechanics (IACM) Award (2008), the Asia-Pacific Association for Computational Mechanics (APACM) Zienkiewicz Medal (2007), the Japan Association for Computational Mechanics (JACM) Computational Mechanics Award (2004), and Honorary Doctor endowed from Iasi Technical University (Romania) (2010). He is also an Honorary Member of JSME, and a Fellow of IACM, ASME, JSIAM and AESJ.

![Surface patches and nodes.](pjab-87-115-g001){#fig01}

![Conception of FMM and node-based data structure.](pjab-87-115-g002){#fig02}

![Inconsistency in satellite elements.](pjab-87-115-g003){#fig03}

![Node-by-node assembling procedure.](pjab-87-115-g004){#fig04}

![Inconsistency caused by degeneracy.](pjab-87-115-g005){#fig05}

![Scheme of Delaunay triangulation based on gift-wrapping method in two-dimensional space.](pjab-87-115-g006){#fig06}

![Scheme of Delaunay triangulation based on gift-wrapping method in three-dimensional space.](pjab-87-115-g007){#fig07}

![Scheme of constrained Delaunay triangulation based on gift-wrapping method in two-dimensional space.](pjab-87-115-g008){#fig08}

![Treatment of slivers: (a) before treatment, (b) after treatment.](pjab-87-115-g009){#fig09}

![Domain decompositions for parallel computing: (a) element-based domain decomposition, (b) node-based domain decomposition.](pjab-87-115-g010){#fig10}

![Communication-independent and communication-dependent nodes.](pjab-87-115-g011){#fig11}

![Parallelization of calculation of product of matrix by vector.](pjab-87-115-g012){#fig12}

![Internal and external satellite nodes.](pjab-87-115-g013){#fig13}

![Comparison of meshes generated for a plate with a central hole: (a) Algorithm I, (b) Algorithm II, and (c) Algorithm III.](pjab-87-115-g014){#fig14}

![Comparison of meshes for combination of two different nodal densities: (a) Algorithm I, (b) Algorithm II, and (c) Algorithm III.](pjab-87-115-g015){#fig15}

![Mesh for a crack model using Algorithm: (a) surface patches and nodes given as input data, and (b) resulted mesh.](pjab-87-115-g016){#fig16}

![Process of mesh generation on a three-dimensional crack model using Algorithm III.](pjab-87-115-g017){#fig17}

![Process of three-dimensional mesh generation for a Pantheon model by Algorithm III (viewed from the outside).](pjab-87-115-g018){#fig18}

![Process of three-dimensional mesh generation for a Pantheon model by Algorithm III (viewed from the inside).](pjab-87-115-g019){#fig19}

![Process of three-dimensional mesh generation for a Pantheon model by Algorithm III when central node are selected in the order determined by multi-level buckets (viewed from the outside).](pjab-87-115-g020){#fig20}

![Large-scale finite element meshes for a Pantheon model generated by Algorithm III: (a) approximately 1.2 × 10^5^ nodes and 4.7 × 10^5^ elements, and (b) approximately 2.1 × 10^7^ nodes and 1.2 × 10^8^ elements.](pjab-87-115-g021){#fig21}

![Conception of enriched free mesh method.](pjab-87-115-g022){#fig22}

![Beam bending model for error norm study: only the simplest case of 2 elements shown here as an example.](pjab-87-115-g023){#fig23}

![Error norm of displacement vs. DOF.](pjab-87-115-g024){#fig24}

![Error norm of energy vs. DOF.](pjab-87-115-g025){#fig25}

![Models for patch test.](pjab-87-115-g026){#fig26}

###### 

Results of patch test: displacements at internal nodes M1 and M2

                        Model A   Model B   Model C                                         
  --------------------- --------- --------- --------- --------- -------- --------- -------- ---------
  FEM (linear)          0.9000    −2.4000   1.2000    −3.0000   0.6000   −2.4000   1.3000   −1.8000
  FEM (quadratic)       0.9000    −2.4000   1.2000    −3.0000   0.6000   −2.4000   1.3000   −1.8000
  EFMM-LS (linear)      0.9000    −2.4000   1.1984    −3.1698   0.6113   −2.3718   1.2849   −1.8990
  EFMM-LS (quadratic)   0.9000    −2.4000   1.1975    −3.1409   0.6108   −2.3681   1.2895   −1.8588
  EFMM-HR (linear)      0.9000    −2.4000   1.2000    −3.0000   0.6000   −2.4000   1.3000   −1.8000
  EFMM-HR (quadratic)   0.9000    −2.4000   1.2000    −3.0000   0.6000   −2.4000   1.3000   −1.8000
                                                                                            
  Exact                 0.9000    −2.4000   1.2000    −3.0000   0.6000   −2.4000   1.3000   −1.8000

[^1]: (Communicated by Takao INUI, M.J.A.)
