The present invention is concerned with a computer implemented method of remeshing patches in a triangular meshed surface, in particular, a method that results in creating a remeshed surface.

Computer-Aided Engineering (CAE) is used increasingly in manufacturing operations to aid in the creation, modification or optimisation of a design. This is a critical step in product design, and CAE offers the opportunity to analyse the design at all stages of the Product Lifecycle Management (PLM). CAE software almost always requires mesh generation software, where this is used as part of finite element analysis or computational fluid dynamics as required by the design. Mesh generation (also known as grid generation, meshing or gridding) takes input from various types of models, including solid modelling (parametric or direct), geometric modelling, NURBS (non-uniform rational basis spline), B-rep (boundary representation), STL (a stereolithography file format) or a point cloud. The mesh is, in its basic form, a subdivision of a continuous geometric space into discrete geometric and topological cells, creating a discrete local approximation of a larger domain.

A mesh is typically formed from simple cells such as triangles or quadrilaterals, with many algorithms proposed over the last few decades for surface meshing and remeshing. However, few of these algorithms are designed to align triangles to both surface curvature and boundaries, thereby minimizing the deviation from the original surface being modelled. There may also be issues with the uniformity of the triangles within the mesh, and lack of conformity with the original imposed variable target size of the surface. One solution proposed to alleviate some of these issues is the use of an advancing front approach for generating triangular surface meshes, such as that discussed in "High Quality Mesh Generation Using cross and Asterisk Fields: Application on Coastal Domains", Georgiadis et al, 26 ^th International Meshing Roundtable, https://imr.sandia.gov/papers/imr26/2007 imr26RN Georgiadis.pdf. This paper proposes the use of direction fields and a frontal point insertion strategy. Asterisk fields are used to generate triangulations and cross fields are used to generate right-angled triangulations optimal for transformation into quadrilateral meshes. However, this paper only considers spherical surfaces, where curvature is equivalent in all directions, and whilst vertices are defined, no information is given on the generation of the triangles in the mesh subsequent to the positioning of an apex. Consequently, whilst this method offers a solution for the spherical surfaces illustrated there are many situations where non-uniform surfaces need to be modelled and meshes generated accurately and with high quality to enable this, such as in finite element analysis (FEA) and finite volume analysis (FVA). There therefore exists a need to be able to create quality meshed surfaces that closely resemble an original surface being simulated as part of an FEA/FVA process.

The present invention aims to address these issues by providing, in a first aspect, a computer implemented method of remeshing patches in a triangular meshed surface, wherein an initial surface IS comprises a triangular surface mesh M and a target size field TS that specifies an optimal triangle edge length at each position on the triangular surface mesh M defined over every point in M, the triangular surface mesh M being partitioned into patches, wherein each patch comprises a contiguous set of adjacent faces delimited by closed loops of boundary or feature edges and has principal surface curvatures, the method comprising the steps of: a) meshing each patch using an advancing front process, wherein the advancing front process comprises: b) computing an asterisk field AF over the patch; c) initialising a front F comprising one or more closed chains of connected edges by the boundary or the features of the patch; d) generating a new triangle by selecting a front edge AB from the edges of the front and determining a position of the apex C of the new triangle by: i) calculating the radius of a first circle based on the target size at point A, CA = TS(A), calculating the radius of a second circle based on the target size at point B, re = TS(B); ii) calculating the growth direction gd for the front edge AB from the asterisk field AF by calculating the asterisk field AF(X) at the midpoint X of the front edge AB, calculating the conjugate asterisk field CAF(X) ofAF(X) as the field obtained from a 30 degree rotation of AF(x), calculating normal n as a vector orthogonal to front edge AB and aligned with the direction of the advancing font, obtaining the growth direction gd by determining the direction of the conjugate asterisk field CAF(X) having the greatest dot product with the normal n; iii) determining the intersection Co of two circles CA centred on A with radius CA and CB centred on B with radius re that lies in the positive half plane of the edge AB with respect to growth direction gd; iv) generating the circumcircle cc of triangle ABCo; v) determining the point O of intersection between the circumcircle cc and the axis of the front edge AB; vi) constructing a line r passing through the point O with the orientation of the growth direction gd; vii) obtaining the position of an apex C of the new triangle at the point where r intersects the circumcircle cc other than at point O; e) if the surface mesh M is not planar, projecting the apex C onto the initial surface IS; and f) inserting point C and enforcing the edges AC and BC into the mesh and replacing any internal triangle in the loop AB, BC, CA, with a single triangle ABC; wherein the asterisk field AF is generated from a cross field CF aligned to the surface principal curvature directions and boundary orientation of the patch.

By generating an asterisk field AFfrom the cross field CFan advancing front approach can be used to remesh a surface, where the asterisk field AF is aligned to the surface principal curvature directions and the boundary orientation of the patch in question. This leads to an increase in mesh quality when compared with prior art approaches.

In one case neither of the contiguous front edges of front edge AB lies at an angle & of less than 4/9n radians (80°) to front edge AB. In another case, if one of the two front edges contiguous with front edge AB, edge s, forms an angle of less than 4/9n radians (80°) with front edge AB, the first approximation of the apex C of the new triangle is the vertex of the end of edge s that is neither A nor B; and the method further comprises: optimising the position of the final apex C by taking the average of all of the apex positions C generated in the iteration of step d). Alternatively, when both of the contiguous front edges of front edge AB lie at an angle of less than 4/9n radians (80°) to the front edge AB, the front edge laying at the smallest angle to front edge AB is selected as the edge s, and the method further comprises: optimising the position of the final apex C by taking the average of all of the apex positions C generated in the iteration of step d).

The method preferably further comprises the step of: carrying out a preliminary edge check on the front edge AB comprising determining the length LAB of the front edge AB and comparing with the optimal length of the front edge AB with respect to the target size field TS, LTS, and either if the length LAB is too short, the front edge AB is collapsed, or if the length LAB of the front edge AB is too long, the front edge AB is split.

Preferably, the method further comprises the step of: placing all edges into a priority queue to determine the order in which new triangles are generated, and sorting the order of the edges in the priority queue by: determining if the angle between the edge and at least one of is contiguous edges is less than 4/9n radians (80°); determining the number of rows of triangles created whilst moving towards the centre of the patch; and determining the iteration number of the front edge, wherein the iteration number of the front edge is zero.

Preferably, the method further comprises the steps of: if the internal angles of the new triangle are measured and one of these internal angles is below a pre-determined threshold, collapsing the triangle.

Preferably, the generation of the asterisk field Ffrom the cross field CF comprises the steps of: generating the cross field using the boundary and feature edge orientations of the patch as constraints and the directions of the principal surface curvature of the patch as the internal guiding field; selecting a primary and a secondary axis from each pair of axes of the cross field CF calculated across the patch; and designating the primary cross field CF axis as a first asterisk field AF axis direction; determining a second asterisk field Faxis direction by rotating 30° counter-clockwise from the secondary cross field CF axis; and determining a third asterisk field AF axis direction by rotating 30° clockwise from the primary cross field CFaxis.

Preferably, if the target size TS is uniform over the triangular surface mesh M, then rA = r _B .

Preferably, if the value of the asterisk field is zero, then apex C= Co, where Co is the intersection of the two circles CA centred on A with radius rA and c _B centred on B with radius re

In another aspect, the present invention also provides a computer-implemented method of creating a remeshed surface, wherein an initial surface IS comprises a triangular surface mesh M and a target size field TS that specifies an optimal triangle edge length at each position on the triangular surface mesh M defined over every point in M, comprising; partitioning the mesh M into patches, wherein each patch comprises a contiguous set of adjacent faces delimited by closed loops of boundary or feature edges; remeshing each of the patches using the method above; and once all patches have been meshed, smoothing the final mesh.

In yet another aspect, the present invention also provides a computer program comprising instructions, which, when the program is executed on a computer, causes the computer to carry out the steps of the method. In yet a further aspect, the present invention also provides a non-transient computer program storage medium, comprising code that when executed on a computer causes the computer to carry out the steps of the method.

The invention will now be illustrated by way of example only, and with reference to the accompanying drawings, in which:

Figure 1 is a diagrammatic illustration of the derivation of an asterisk field from a cross field;

Figure 2a is a diagrammatic illustration of possible front edge positions in the generation of a triangle within a mesh when neither angle ABD nor CAB are less than 4n/9 radians (80°);

Figure 2b is a diagrammatic illustration of possible front edge positions in the generation of a triangle within a mesh when one of angle ABD or CAB is less than 4n/9 radians (80°);

Figure 2c is a diagrammatic illustration of possible front edge positions in the generation of a triangle within a mesh when both angle ABD and CAB are less than 4n/9 radians (80°);

Figure 3 is a diagrammatic illustration of apex position during the process of generating a triangle in accordance with embodiments of the present invention having Status 0:

Figure 4a illustrates an example leading edge AB and triangle apex C; and

Figure 4b illustrates the position of a new triangle apex, C‘.

Embodiments of the present invention offer advantages over the prior art as they are able to produce a high quality, uniformly sized mesh with triangles aligned to surface boundaries and curvature thus minimising any deviation from the original surface. In addition, the generated mesh size conforms closely to the imposed variable target size - the desired size of the edges in the mesh. This is due to the adoption of an advancing front approach for mesh generation driven by an asterisk field aligned to the surface principal curvature directions and the boundary orientation of a patch formed by partitioning the mesh. To achieve such an alignment the present invention takes the approach of converting a cross field into an asterisk field. This is achieved by using a computer implemented method of remeshing patches in a triangular meshed surface, wherein an initial surface IS comprises a triangular surface mesh M and a target size field TS that specifies an optimal triangle edge length at each position on the triangular surface mesh M defined over every point in M. The triangular surface mesh M is partitioned into patches, wherein each patch comprises a contiguous set of adjacent faces delimited by closed loops of boundary or feature edges and has principal surface curvatures. The method comprises a number of steps, beginning initially with meshing each patch using an advancing front process. The advancing front process involves the computation of an asterisk field AF over the patch and then initialising a front F comprising one or more closed chains of connected edges by the boundary or the features of the patch. Once this has been done, new triangles are iteratively generated by selecting a front edge AB from the edges of the initial front as a basis, the position of the apex of the new triangle C is determined using a circumcircle method. Initially, the radius of a first circle based on the target size at point A, CA = TS(A) is calculated, and the radius of a second circle based on the target size at point B, re = TS(B) is calculated. A growth direction gd is then calculated for the front edge AB from the asterisk field AF. This is done by calculating the asterisk field AF(X) at the midpoint X of the front edge AB, calculating the conjugate asterisk field CAF(X) ofAF(X) as the field obtained from a 30 degree rotation of AF(X), and then calculating normal n as a vector orthogonal to front edge AB and aligned with the direction of the advancing font. The growth direction gd s then obtained by determining the direction of the conjugate asterisk field CAF(X) having the greatest dot product with the normal n. The intersection Co of two circles CA centred on A with radius CA and CB centred on B with radius re that lies in the positive half plane of the edge AB with respect to growth direction gd is determined. The circumcircle cc of triangle ABCo is created and the point O of the intersection between the circumcircle cc and the axis of the front edge AB is determined. At this point a line r passing through the point O with the orientation of the growth direction gd is determined such that apex C of the new triangle at the point where r intersects the circumcircle cc other than at point O is found. If the surface mesh M is not planar, point C is projected onto the initial surface IS. The point C is inserted into the mesh M and the edges AC and BC are enmeshed into the mesh M so that any internal triangle in the loop AB, BC, CA are replaced with a single triangle ABC. The key feature of this method is that the asterisk field AF is generated from a cross field CF aligned to the surface principal curvature directions and boundary orientation of the patch. It is this, in conjunction with the methodology by which the individual triangles are calculated in the remeshing process that results in the quality of the meshes generated. The detailed embodiments of the invention described below all show such improved remeshing behaviour.

Asterisk Field Generation

As described above, the embodiments of the present invention employ an asterisk field AF generated from a cross field CF. An initial surface IS comprising a triangular surface mesh M and a target size TS that specifies an optimal triangle edge length at each position on the triangular surface mesh M defined over every point in M is created. The triangular surface mesh M is partitioned into patches, wherein each patch comprises a contiguous set of adjacent faces delimited by closed loops of boundary or feature edges. Each patch also has principal surface curvatures. To generate a cross field CF, a complex linear system is assembled, where the boundary or feature edges and principal surface curvatures of the patch are used to define the alignment energy for the interior of the patch domain. The cross field CF is then computed for the vertices of the patch using the smooth n-direction fields theory, where n = 4, as described in "Globally Optimal Direction Fields", F. Knbppel et al, http://www.cs.cmu.edu/~kmcrane/Proiects/GloballyOptimalDirec tionFields/paper.pdf, and transferred to the centroid of the patch. Whilst a cross field CF can be used to drive the advancing front of a quad mesh, it is not suitable for triangle meshing as it is only able to drive the advancing front in four directions. Triangle meshing instead requires the use of an asterisk field AF, since this is able to drive the advancing front in six directions. In a perfect equilateral triangle mesh every internal vertex is attached to six edges and six faces, hence the need for an asterisk field AF with six directions to drive the advancing front. Whilst it is also possible to apply the "Globally Optimal Direction Fields" theory to create an asterisk field AF this does not lead to expected results. This is because an asterisk field AF cannot be aligned to the surface principal curvature directions as there are only two axes (four directions) and not three axes (six directions). Furthermore, an asterisk field AF does not align well to patches having rectangular boundaries (and, in general, boundaries with 90° angles), which are by far the most commonly used in CAE applications. It is for these reasons that the present invention adopts the approach of firstly, computing a cross field CF and secondly, converting to an asterisk field AF.

Figure 1 is a diagrammatic illustration of the derivation of an asterisk field AF from a cross field CF. Initially, for a chosen patch, a cross field CF is generated using the boundary and feature edge orientations of the patch as constraints and the directions of the principal surface curvatures of the patch are used as the internal guiding field. This creates a pair of orthogonal, cross field CF axes as illustrated in Figure 1. To create the asterisk field AFthe one of the axes from the pair of cross field (Taxes calculated across the patch is selected as a primary axis P and the other as a secondary axis S. The primary axis P is designated as a first asterisk field AF axis direction AFi. The second asterisk field AFaxis direction AF2 is determined by rotating the 30° counter-clockwise from the secondary cross field CF axis S. The third and final asterisk field AF axis direction AF3 is then determined by rotating 30° from the secondary cross field ( axis S. For curved surfaces, the minimum absolute curvature is preferably chosen as the principal direction 4, however in some cases it may alternatively be desirable to choose the maximum absolute curvature as the primary cross field CF axis P. For planar surfaces, where there is no curvature information, primary cross field CF axis P is imposed by the boundary orientation, and this information is propagated into the interior of the domain using a flood fill algorithm. Once the asterisk field AF has been created, the remeshing of the initial surface IS can begin by generating new triangles to replace those in the triangular meshed surface M.

Triangle Generation

In triangular mesh generation, each time the front edge is processed the mesh is modified by building a new triangle that has the front edge as one of its edges. This is illustrated further in Figure 2. Figure 2 is a diagrammatic illustration of possible front edge positions in the generation of a triangle within a triangular surface mesh M. The front edge is denoted as AB, with the previous edge in the loop denoted as CA and the next edge in the loop denoted BD. This illustrates that a front F com prises one of more closed chains of connected edges by the boundary or the features of the patch. Depending on angles between the leading edge AB and its two contiguous edges BD and

CA there are three possible cases: i) Status 0: Neither angle ABD and CAB is less than 4n/9 radians (80°) (shown in Figure 2a); ii) Status 1: The angle ABD is less than 4n/9 radians (80°), the angle CAB is greater than 4n/9 radians (80°) or vice versa (shown in Figure 2b); iii) Status 2: the angles ABD and CAB are both less than 4n/9 radians (80°) (shown in Figure 2c).

The angle & that the contiguous edges make with the leading edge AB is key in determining the complexity of the calculation of the position C of the apex of the new triangle to be generated by the advancing front. These are outlined below.

Status 0

Figure 3 is a diagrammatic illustration of apex position during the process of generating a triangle in accordance with embodiments of the present invention. When neither of the contiguous front edges of front edge AB lies at an angle tf of less than 4/9n radians (80°) to front edge AB, the position of apex of the new triangle C is determined by the following steps.

1. calculating the radius of a first circle based on the target size at point A, rA = TS(A), calculating the radius of a second circle based on the target size at point B, rB = TS(B);

2) calculating the growth direction gd for the front edge AB from the asterisk field AF by calculating the asterisk field AF(X) at the midpoint Xof the front edge AB, calculating the conjugate asterisk field CAF(X) ofAF(X) as the field obtained from a 30 degree rotation of AF(x), calculating normal n as a vector orthogonal to front edge AB and aligned with the direction of the advancing font, obtaining the growth direction gd by determining the direction of the conjugate asterisk field CAF(X) having the greatest dot product with the normal n;

3) determining the intersection Co of two circles CA centred on A with radius rA and CB centred on B with radius re that lies in the positive half plane of the edge AB with respect to growth direction gd;

4) generating the circumcircle cc of triangle ABCo; 5) determining the point O of intersection between the circumcircle cc and the axis of the front edge AB;

6) constructing a line r passing through the point O with the orientation of the growth direction gd;

7) obtaining the apex C of the new triangle at the point where r intersects the circumcircle cc other than at point O;

8) if the surface mesh M is not planar, projecting point C onto the initial surface IS; and

9) inserting point C and enforcing the edges AC and BC into the mesh and replacing any internal triangle in the loop AB, BC, CA, with a single triangle ABC.

Status 1

In the case of Status 1 only one of the contiguous front edges of front edge AB lies at an angle & of less than 4/9n radians (80°) to front edge AB. This edge is nominated as the edge s that will be used to calculated the apex of the new triangle. Since there are already two edges (the leading edge AB and the edge s) that are suitable for forming a new triangle, the vertex of this new triangle C is given by the vertex of the end of edge s that is neither/5 nor B. Once chosen, the apex position is optimised by calculating the average position of all apex locations constructed for the triangles attached to C that have been generated during step d) above.

Status 2

When both of the contiguous front edges of front edge AB lie at an angle of less than 4/9n radians (80°) to the front edge AB, the front edge laying at the smallest angle to front edge AB is selected as the edge s, and the position of the apex C is generated and optimised as in Status 1.

In the above examples the threshold chosen to determine whether a front edge is 4/9n radians, or 80°. However, other angles may be suitable or desirable, and preferably fall into the range of being less than n/2 radians, or 90°. The reasoning for this is as follows. The simplest example to consider is that shown in Figures 4a and 4b. Figure 4a illustrates an example leading edge AB and triangle apex C, and Figure 4b illustrates the position of a new triangle apex, C‘. In Figure 4a, a represents the angle ABC, with the values of a determining whether the side edge BC is the most appropriate edge to use for a new triangle or whether it is better to discard it. If the side edge BC is suitable, A will be connected with C, creating a triangle ABC. However, if not, a new edge BC' will be created with an angle of a/2, resulting in a triangle ABC' on the first iteration, and on the next iteration a triangle BCC', also with an angle of a/2, will be created.

An ideal triangle, e.g. an equilateral triangle, has an internal angle of n/3 radians, or 60°. Therefore, to choose between the two possibilities in Figures 4a and 4b the difference between the internal angle on the vertex B of the generated triangle and n/3 (60°) should be minimised. The difference between these two values in the case of Figure 4a is given by: d(a) = | a-n/31 (1) whilst in the case of Figure 4b, the difference is given by: d'(a) = |a/2-n/3 | (2)

In order to choose the simplest case to generate triangle ABC (Figure 4a) d(a) < d'(a) (3) In order to solve these equations for a (1) and (2) are substituted into (3) to give:

| a-n/3 | < | a/2-n/3 | (4)

For n/3 < a < 2n/3 a - n/3 < n/3 - a/2 (5) leading to: a < 4n/9 = 80°

Whilst the above considers the effect of the angle & on the generation of the triangle, it is also necessary to consider the effect of the target size TS of the patch on the triangle generated. The target size TS may be uniform over the triangular surface mesh M, resulting in the values of rA and re being equal and the intersection Co of circles CA, CB forms an equilateral triangle ABCo. If the target size TS is variable, the method continues as above with individual values of rA and re. Similarly, if there is no imposed asterisk filed AF such that the value of the asterisk field AF is zero, the apex C= Co where Co is the intersection of the two circles: CA centred on A with radius rA, and CB centred on B with radius re.

In use, it may be desirable to perform further operations as part of the method of the present invention. For example, before beginning any calculations it may be desirable to carry out a preliminary edge check on the front edge AB. This comprises determining the length LAB of the front edge AB and comparing with the optimal length of the front edge AB with respect to the target size field TS, LTS. If the length LAB is too short, the front edge AB is collapsed, and if the length LAB of the front edge AB is too long, the front edge AB is split. It may also be necessary or desirable to consider the internal angles of the new triangle in the triangular surface mesh M that has been generated. If the internal angles of the new triangle are measured and one of the internal angles is below a predetermined threshold, the triangle is collapsed. This prevents the generation of triangles where the quality of the resulting remeshed surface is compromised. Both of these steps may be combined or alternated with a step that determines the priority in which leading edges will be processed by the advancing front. For example, all the edges of the front F are placed into a priority queue to determine the order in which the new triangles are generated. The order of the edges in the priority queue is given by looking at various criteria. Firstly, the status: a determination of whether the angle between the edge and at least one of is contiguous edges is less than 4/9n radians (80°). A higher status has a higher priority, such that Status 2 takes priority over Status 1, which takes priority over Status 0. Secondly, the level: the number of rows of triangles created whilst moving towards the centre of the patch is calculated. Thirdly, the age of the front edge: an integer value that specifies the iteration in which the front edge was created, wherein initial edges have an age equal to zero. Of course, other methods of determining the priority of the edges may be used in place of these, if desired.

The method outlined above finds particular use in computer implemented methods of creating a meshed surface. This is where the initial surface IS comprises a triangular surface mesh M and a target size field TS that specifies an optimal triangle edge length at each position on the triangular surface mesh M defined over every point in M. To begin with, the method comprises partitioning the mesh M into patches, wherein each patch comprises a contiguous set of adjacent faces delimited by closed loops of boundary or feature edges before remeshing each of the patches using the method outlined in the embodiments described above iteratively. Then, once all patches have been remeshed, a step smoothing the final mesh using known techniques can be applied. To enable this, both a computer program comprising instructions, which, when the program is executed on a computer, causes the computer to carry out the steps of the methods, or a nontransient computer program storage medium, comprising code that when executed on a computer causes the computer to carry out the steps of the methods can be provided.