Mesh generation


2D Quadrilateral Meshing 

Automatic generation of 2D quadrilateral meshes is based on the subdivision of an arbrarily complex domain into simpler meshable subregions using the Medial Axis Transform or skeleton of the region (Tam & Armstrong, 1991). These simple subregions, which have at most 6 sides, are meshed with quad elements using a technique called midpoint subdivision. Compatibility between meshes in adjacent subrergions is enforced by Integer Programming (Tam & Armstrong, 1993). 

3D Hexahedral Mesh Generation


In 3D, hexahedral mesh generation is accomplished using the Medial Surface of a solid to subdivide a complex body into a collection of simpler domains. Each face, edge and vertex of the Medial Surface corresponds to a meshable domain (Price, Armstrong & Sabin, 1994). A 3D midpoint subdivision algorithm is used in conjunction with Integer Programming to create a compatible mesh (Li, Armstrong & McKeag, 1994) Adaptive analysis based on this technique has been demonstrated (Armstrong et al, 1993).


Designers draft and analyse simple conceptual layouts to prove an idea before expending resources on a full-blown 3D finite element analysis model. Once an acceptable solution has been found at the conceptual stage, the original CAD representation may then be re-constructed as a 3D model. A considerable amount of unnecessary duplication of endeavours occurs during this process. Dimensional addition, the reverse process of the more established idea of dimensional reduction , is therefore advocated as an effective saving for the designer, whereby the transition from 1D skeleton to 3D solid at is least semi-automating and so minimising what would otherwise be a tedious and costly drafting stage. 

Geometric operations are necessary to generate a coherent solid model from the conceptual layout. Modern CAD systems are highly productive tools for concept embodiment, detailing, and preproduction engineering. Systems with fully integrated behaviour modelling and simulation functionality provide the engineer with tools capable of supporting a number of engineering steps in all life-cycle phases of a product. The simplified operations developed at QUB FEM facilitate rapid and full re-use of idealised conceptual model details to create full solid CAD and analysis models. They provide functionality for adding complexity to simple models are dimensional addition, detail insertion and symmetry addition, all of which have been demonstrated in the examples shown. 
The integration of analysis modelling into CAD systems has enormous potential benefits. A huge range of geometry construction and modification tools are available, whilst access to the feature tree used to create the model means that optimisation studies can modify design geometry directly rather than the finite element mesh. Thus analysis-driven design is facilitated.


In 2-D the medial axis is the locus of the centre of an inscribed disc of maximal diameter as it rolls around the domain interior expanding and contracting to maintain contact with the domain boundary. The combination of the medial axis and the radius function, which describes the radius of the inscribed disc at any point on the medial axis, is known as the Medial Axis Transform (MAT). In 3-D the equivalent construction is the locus of the centres of all inscribed spheres of maximal diameter. This is also known as the medial axis, though perhaps the medial surface would be a more appropriate description.

The medial axis captures the geometric proximity of the boundary elements in a simple form and therefore provides a complimentary representation of physical objects in computer aided design systems. It is obvious that the effectiveness of the medial axis to capture an object's geometric characteristics influences its ability to serve these purposes, e.g. meshing, features recognition, object decomposition, path planning etc.

Currently there are various ways to construct the medial axis topology. Typically these algorithms generate excessive points such that the medial axis patch is geometrically over approximated, i.e. the number of generated points, which define a particular medial axis patch, is more than what is practically required. Hence it appears that an algorithm for the construction of medial axis geometry, which will be carried out after the topological structure has been constructed, is necessary for medial axis patch approximation. With an efficient algorithm to successfully approximate the medial axis, the subsequent processes will run more effectively.



Develop an adaptive curvature-sensitive mesh generator for ordinary surfaces. For an ordinary surface in 3D space, it has one, or any combination, of the following natures: elliptic, parabolic, hyperbolic and planar. The differences of these surface natures are based on the magnitude and directions of the curvature vectors. The working principle of this mesh generator lies on the first and second fundamental forms of the surface.


Develop a set of theories and formulae so that the geometric properties of the medial axis of 2D planar object can be obtained. Derived from the equal distance criterion of medial axis, this set of theories and formulae is able to generate all the necessary geometric information, such as the tangent and curvature vectors, of a particular point at the medial axis.



Extent the theories and formulae found in stage 2 to the mid surface of 3D solid object. All the geometric information required to do curvature-sensitive meshing on ordinary surfaces is equally essential for mesh generation on the mid surfaces. Some principles found in stage 2 can be extended to the mid surface of 3D solid object, and the coefficients of the first and second fundamental forms of the mid surface have been derived.



Combine all the results to produce an adaptive curvature-sensitive mesh generator for mid surface. By dropping all the equations and functions of stage 3 into the mesh generator developed at stage 1, the mesh generator is able to do adaptive curvature-sensitive meshing on mid surface.




CAD is traditionally seen as the carrier of information in design. However, creating a finite element analysis model of a design invariably requires simplification. Therefore, the analysis model is often rebuilt from scratch, relying upon the judgment of skilled analysts, much of the work in creating the geometry being duplicated. As well as removing details from the design geometry, plate, shell and beam element types are also often used where they can give more accurate and cost effective results.

Procedures have been developed for the automatic dimensional reduction of a 2D geometric model to an equivalent 1D beam analysis model. This was achieved using the medial axis transform, an alternative skeleton-like representation of the geometric model having properties relevant to model simplification. Mixed dimensional models are created implicitly with the required coupling automatically applied.

Operations have also been defined and implemented for the expansion of these operations to the dimensional reduction of 3D solid models to equivalent solid/shell/beam analysis models. These operations are interactive, with appropriate physical properties such as shell thickness, beam section moment of areas and torsion constants calculated automatically. 

The automatic dimensional reduction of thin walled solids to equivalent shell models has been implemented using a geometric measure similar to an aspect ratio, but extended for solid to shell dimensional reduction. This process uses the operations described earlier directly, with shell thicknesses being calculated automatically.


The first stage is to subdivide the 2D model into subregions which are suitable for representation by 1D beam elements, and those which are not to either by 2D elements, or by an equivalent 0D point mass.

A threshold value of aspect ratio can then be introduced, below which a region is not suitable for representation as 1D beam elements and must be represented either by the 2D region or an equivalent 0D point mass. All other regions can then be represented by 1D beams.

Progressive dimensional reduction based on aspect ratio

For the case of mixed 1D/0D abstracted models, the 1D beams can either be connected together with rigid links, or extrapolated to meet at their closest point inside the neighbouring 0D region.


The diagram below shows a 2D axisymmetric section of part of a turbine casing. They were used here to validate the dimensionally reduced beam models created. A natural frequency analysis of both the 2D sections and the abstracted beam models was carried out. The results for the first three modes are presented below. Three analyses were carried out, namely for a full 2D mesh, a 1D model with rigid link coupling and an extrapolated 1D beam model. Also shown is the extrapolated beam model (top left), the beams together with 2D regions represented as point masses in the analysis model (centre left) and the beam model coupled with rigid links (bottom left).

From the analysis results it can be seen that the frequencies obtained from the 2D models are bounded by the reduced model results. The frequencies for the rigid link model are higher, as might be expected since the full mass of the structure is represented, but some of the flexibility has been replaced by rigid links of infinite stiffness. The extrapolated beam models give lower natural frequency results than the full 2D models. In both cases the mode shapes obtained are representative of the 2D model, and the analysis results are good given the small number of degrees of freedom in the reduced models.


A lot of 2D planar models used are cross sections of a full 3D model which has shown some form of of symmetry. Therefore it was logical to extend the creation of analysis models to swept 2D shell / 1D beam stiffener models. These models are created automatically, with all required physical properties and coupling between element types applied. Some simple models are shown below. 

Section swept in Z-direction, coupled with rigid links and extrapolation 

Section rotationally swept about x-axis, rigid links and extrapolation





This model was dimensionally reduced to the equivalent shell model below automatically, no user interaction was required.  

The thin walled solid model below was also dimensionally reduced automatically to a shell model, and meshed using triangular elements. 

Mixed Dimensional Coupling


The complex structure on the left is part of a facility used for offshore oil exploration. As can be seen, this is a framed structure. Frames are major structural components in many engineering structures. Applications vary from simple roof trusses to very complex structures such as the one illustrated.


The design and assessment of such components require consideration in strength, fatigue and fracture. The finite element method is used to simulate the behaviour of joints of complex geometry and welded joints with material heterogeneity. Analysis of complete frame structures is desirable so that load paths and magnitudes are computed correctly. Local stress-strain history, and failure modes and their locations can then be evaluated precisely.

Accurate predictions of static collapse behaviour for aging structures are essential so that reserve strength and likelihood of damage can be assessed. Full 3D elastic-plastic analyses are currently used to predict the responses of framed structures to the point of collapse. Due to the structural complexity, these impractical analyses result in long computing times and hence solution can only be obtained by using expensive supercomputers.


Substantial reductions in analysis times are achieved if idealisation techniques are employed. Due to discontinuities in geometry, loading and material behaviour, the highest stresses are produced at the joints of these structures, and hence failure and collapse usually occurs there. This non-linear behaviour requires full 3D analysis of each joint. Between the frame joints, the structures are long and slender, and can be replaced by beam elements. Beam elements are appropiate as they accurately represent the slender regions and are computationally efficient.


In order to achieve a state of compatibility of displacements and stress equilibrium at the junction or interface between the differing element types, it is important to integrate into the analysis some scheme for coupling the element types in a way that conforms with the theory of elasticity.


The work done at either side of the interface is equal, and so we can write:

  • (force.distance)1d = (force.distance)3d
  • (force.distance)1d = (stress.area.distance)3d



So the essence of the problem is to determine the correct stress distributions on the beam cross section. Extensive research at QUB has resulted in the correct evaluation of these stress distributions for arbitrary sections with any number of holes. The interaction of either side of the interface is achieved by using multi-point constraint equations. These equations are a way of connecting the two sides of the interface, so that when one side moves, the other side moves in a manner that is defined by the constraint equation.

By adopting the above idea, a marked reduction in analysis times can be achieved without any loss of accuracy. At the QUB Finite Element Research Group, focus in multi-dimensional coupling has been on 3D-1D coupling, the mathematical coupling of 1D beam elements to 3D continuum elements at a common interface. Coupling of elements of other dimensions have been considered, and have yet to be developed to full automation.

The beam to solid coupling procedures have been designed to facilitate multi-cellular sections at any orientation.  They also include the ability to considir problems such as asymmetric bending and shear centers. Software has been developed to couple between beams and standard solid mesh types (3 noded & 6 noded triangular based 3D meshes, and 4 & 8 noded quadrilateral based 3D meshes). Click here for comparisons between the coupling technique developed at Queen's and the others that are available. Coupling of shells and solids has also been considered and developed to maturity. The procedure in generating the constraint equations again follows the aforementioned idea, i.e. equating the work done on either side of the dimensional interface. This method of coupling also proves itself to be applicable to problems involving composite laminates.