2D Elliptic Mesh Generator

This is a 2D orthogonal elliptic mesh (grid) generator which works by solving the Winslow partial differential equations (Elliptic PDEs). It is capable of modifying the meshes with stretching functions and an orthogonality adjustment algorithm. This algorithm works by calculating curve slopes using a tilted parabola tangent line fitter. The mesh generator is packaged as a Java program which can be compiled and executed via the command line.

The program allows one to choose from six different boundary types: rectangular, Gaussian, absolute value, greatest-integer, forwards step and semi-ellipse. One can then specify the coordinates of the mesh domain (note: the domain must be perfectly square). Finally, one can choose to add refinements to the mesh, such as orthogonality adjustment and stretching functions. The program will then generate an initial coarse mesh and iteratively refine it to produce a smooth mesh with the given parameters and refinement options.

A distinct feature of the elliptic mesh solver is that it corrects overlapping and misplaced grid lines very well. A detailed analysis of the quality of the resulting mesh will also be provided as part of the program output.

Running the Program

To run the pre-built JAR file of this project, ensure you have Java version >= 1.8 installed and execute the following on the command line: java -jar 2DEllipticMeshGenerator.jar. From here, simply follow the prompts to start the generation of a mesh. As the solver runs, you will see a new grid frame popup for every few iterations completed. This allows you to see the progress of the mesh solution as it converges in real time.

Once the mesh has been generated, you will find 4 new files in the root directory: animatedMesh.gif, Initial_Grid.txt, Final_Grid.txt and Grid_Info.txt.

• animatedMesh.gif will show you the time-evolution of the mesh solution from the initial grid to the final grid in its converged state.

• Initial_Grid.txt contains the physical XY-coordinates for each point in the initial mesh.

• Final_Grid.txt contains the physical XY-coordinates for each point in the final solution.

• Grid_Info.txt contains the statistical quality analysis report for the initial and final meshes.

Contributing

If you'd like to contribute to the development of this project, you may want to build the source files from scratch and test your modifications. I've structured this project using Maven to help simplify the build process and manage dependencies more easily. To build the project, follow these steps:

1. Ensure you have the latest version of Maven installed (https://maven.apache.org/download.cgi).
2. Once you're done editing the source files, within the root directory of the repo, run mvn package.
3. If the build completes successfuly, then you will find an updated executable JAR with the following path: target/2DEllipticMeshGenerator-1.0-SNAPSHOT-jar-with-dependencies. To run your changes, execute java -jar target/2DEllipticMeshGenerator-1.0-SNAPSHOT-jar-with-dependencies.

Screenshots

Here are some examples of meshes generated with the program (initial in blue and final in green):

                            

A more complete collection can be found within the Screenshots folder.

Animations

The application automatically generates an animation of the mesh as it converges in real-time. This is contained within the produced file animatedMesh.gif. By default, the animation shows the solution at every 5 iterations of the solver per frame. This can be adjusted in the main class, EllipticMeshGenerator2D through the constant animationFPS.

Here are some examples of animations that are generated by the program:

A more complete collection can be found within the Animations folder.

Mathematical Framework

The algorithms implemented in this project are mainly founded upon the principles of differential geometry and tensor calculus. The most fundamental concept behind the mathematics governing this project is the transition between coordinate systems in order to obtain the solution to partial differential equations in the most efficient manner possible. More specifically, in the context of this project, this implies transforming a set of equations written in Cartesian coordinates to curvilinear coordinates. This concept can be extended to any number of spatial dimensions, which will later be shown. However, a two-dimensional solution was developed in order to demonstrate the feasibility of generating smooth elliptic grids with the prescribed specifications.

Consider a system of n dimensions which can be represented by the set of Cartesian coordinates and the set of curvilinear coordinates . We define the covariant metric tensor g_ij to be:

where each of the partial derivatives comes from the definition of the covariant base vectors. Each of these base vectors describe how one coordinate system changes with respect to another, when any particular coordinate is held fixed. For our two-dimensional problem, we can expand these sums to yield the following expressions for the metric tensors:

   

   

where (x, y) represents our Cartesian coordinate system and represents our curvilinear coordinate system. We could also define the contravariant base vectors and metric tensors if we wished to solve the inverse problem, which would be transitioning from curvilinear to Cartesian coordinates.

Elliptic Mesh Generation Algorithm

Firstly to construct an initial mesh, the Transfinite Interpolation algorithm is applied to the given domain constrained by the specified boundary conditions. This algorithm is implemented by mapping each point within the domain (regardless of the boundaries) to a new domain existing within the boundaries. This algorithm works by iteratively solving the parametric vector equation

where and represent parameters in the original domain and , , and represent the curves defining the left, top, bottom and right boundaries. P_ij represents the point of intersection between curve and .

At the heart of the solver is the mesh smoothing algorithm, which at a high level, works by solving the pair of Laplace equations

and

where and represent the x and y coordinates of every point in the target domain, mapped to a transformed, computational space using the change of variables method. This renders the calculations simpler and faster to compute. However, we wish to solve the inverse problem, where we transition from the computational space to the curvilinear solution space. Using tensor mathematics, it can be shown that this problem entails solving the equations

and

where g_ij is the covariant metric tensor at entry (i,j) within the matrix of covariant tensor components defining the mapping of the computational space coordinates onto the physical solution space coordinates (x,y). In this model, x and y are computed as functions of and .

This set of equations are the elliptic PDEs known as the Winslow equations. These are applied to the mesh using the method of mixed-order finite differences on the partial derivatives (and tensor coefficients, as they are a function of these derivatives), thereby resulting in the equations (for a single node):

and

,

where i and j are the coordinates of a node in the mesh in computational space. Here and are equal increments in and respectively.

The coefficients for these equations can be generated for each point to form a system of linear equations, which is then modeled in matrix representation, resulting in a tri-diagonal matrix. This matrix is then solved iteratively using the Thomas Tri-Diagonal Matrix Algorithm line-by-line by traversing from the bottom up.

The solution to the matrix generated from a single iteration can then be further processed by the orthogonality adjustment algorithm and stretching function methods as necessary. The solver then calculates the solution for all other node lines and repeats the process until the difference between adjacent nodes meets a threshold convergence criteria.

Orthogonality Adjustment Algorithm

In several computational fluid dynamics applications, an orthogonal mesh is necessary in certain regions to ensure a high enough accuracy when performing calculations. However, it is not always possible to achieve a fully orthogonal solution, and thus the problem becomes finding a nearly-orthogonal solution to an arbitrarily defined domain.

The implemented solution uses an iterative approach to find the angles of intersection and adjust the position of the nodes until their respective angles of intersection converge to a reasonable threshold value from 90 degrees. The exact method makes use of the linear approximation of the grid lines intersecting at each node within the mesh.

A remarkable result from the research was the development of an accurate method for obtaining these linear approximations. This method consists of fitting a tilted parabola to three adjacent nodes, which are defined as (x₁, y₁), (x₂, y₂) and the node in between these two. By applying coordinate transformations, we can obtain the trigonometric function

whose roots can be solved for using the bisection method, which represent the angular position of the parabola (can be improved with Newton's method).

The same process is applied to the three oppositely adjacent nodes. From this, a suitable linear approximation is obtained, and the adjustment is determined by plugging the slopes of the two linear functions into the linear equation relating the two derived using basic analytic geometry. This describes the system of equations

for the vertical grid line V and the horizontal grid line H at a given iteration k. Thus, this system is solved iteratively for each mesh node in the same line-by-line fashion as the Winslow system solver.

Stretching Functions

In order to further improve the quality of the mesh, one can introduce univariate stretching functions to either compress or expand grid lines in order to correct non-uniformity where grid lines are more or less dense. These functions are arbitrarily chosen and only reflect the distribution of grid lines.

We can derive a new set of equations by combining our previously established differential model for grid generation and a set of univariate stretching functions of our choice. In order to do so in a straightforward manner, we can transform our Cartesian coordinates (x, y) to a new set of coordinates (χ, σ) which exists in a different space, χσ, called the parameter space. We define our stretching functions, f₁ and f₂ as bijectional univariate functions of ξ and η respectively. They are described as follows:

and

If we define the mapping from physical space to parameter space as the same elliptic system as before, we get the equations:

and

To obtain an expression for this system, we can get the first order partial derivative with respect to x as so:

Using the chain rule, we can get the second order partial derivative:

Similarly for y, we get:

Adding the previous two equations, and rearranging, we get:

Similarly for η, we get:

However, we wish to solve the inverse problem as before, and thus in finding the inverse of these equations, we obtain the Poisson system:

However, if we wish to maintain orthogonality everywhere in the grid, then we must eliminate the g₁₂ terms to get the equations:

This slightly modifies the solution process by changing the values of the matrix coefficients in the TDMA setup.

The discretized version of the new system is

and

where the metric tensor coefficients are the same as before.

The default stretching functions used in the program are

and

where |∝| << 1 and |β| << 1. The parameters ∝ and β, when positive, indicate how much stretching occurs in the x and y directions respectively. When these values are negative, compression of grid lines will occur instead.

Extension to Higher Dimensions

If we wished to extend the elliptic solver to 3D, we would need to develop equations for transitioning from a curvilinear coordinate system (ξ, η, ς) to a 3D Cartesian coordinate system (x, y, z). Using the same elliptic model as before, we get the 3D version of the Winslow equations

where each of the metric tensor coefficients is determined by taking the cofactors of the contravariant tensor matrix. The contravariant tensor matrix is used to obtain the coefficients for the Winslow equations, which are the inverse of the Laplace equations as stated before.

In general, if we wish to extend our elliptic mesh solver to n dimensions, then we will have n sets of equations each with n!/(2(n-2)!) + n terms. This renders the problem gradually more and more difficult to solve for higher dimensions with the existing elliptic scheme, implying that a different type of PDE might be needed in these cases.

Another complication that arises in higher dimensions is adjusting grid lines to enforce orthogonality. Using the aforementioned algorithm for adjusting grid lines to achieve either complete or partial orthogonality on the boundary, we would need to iteratively solve three sets of two linear equations for each node in the mesh, as well as solve three trigonometric equations per iteration to compute the tangents.

For two dimensions, if our grid domain was of size n x n, and we assumed an average of k iterations per mesh node, then our runtime complexity would be approximately . However for three dimensions, our runtime complexity would be . If our domain had a side length of 100 and for each node, the solver took 10 iterations in the 2D version and 25 iterations in the 3D version, then the 3D version would take over 1000 times longer to return a solution. This means we would definitely need to seek a more efficient and powerful algorithm. However, another issue that may arise is that for three dimensions, orthogonal solutions to smooth meshes are very rare, and even getting a solution with the most efficient algorithm might be impossible for many cases.

Mesh Quality Analysis Report

In order to determine the quality of the resulting mesh, it was necessary to construct an objective means of quality measurement. Therefore, several statistical procedures were implemented in the program to produce a meaningful mesh quality analysis report. The metrics which are presented are divided into the following categories:

• Orthogonality Metrics

  • Standard deviation of angles
  • Mean angle
  • Maximum deviation from 90 degrees
  • Percentage of angles within x degrees from 90 degrees (x can be set as a constant in the code)

• Cell Quality Metrics

  • Average aspect ratio of all cells
  • Standard deviation of all aspect ratios

In the above list, "angle" refers to the angle of intersection of grid lines at a node and "aspect ratio" refers to the skewness of a grid cell measured as a ratio of the cell's longest side to its shortest side.

Libraries Used

• JMathPlot by Yann Richet

Recognition

• Cited in the CFD Open Series textbook, Structure Meshing for CFD

License

© Chaitanya Varier 2017. This software is protected under the MIT License.

Context and Acknowledgement

I wrote this software and paper independently for Prof. Jonathon Sumner at Dawson College, who provided me with invaluable mentorship and advice throughout the course of my research.

References

I gleaned all the necessary mathematical knowledge for completing this project from the following sources:

• Farrashkhalvat, M., and J. P. Miles. Basic structured grid generation with an introduction to unstructured grid generation. Oxford: Butterworth Heinemann, 2003. Print.

• Versteeg, H.K, and W. Malalasekera. An introduction to computational fluid dynamics: the finite volume method. Harlow: Pearson Education, 2011. Print.

• Knupp, Patrick M., and Stanly Steinberg. Fundamentals of grid generation. Boca Raton: CRC Press, 1994. Print.