This is a powerful 2D orthogonal elliptic mesh (grid) generator which works by solving the Winslow partial differential equations. 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 (original discovery). 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. Then one must specify the coordinates of the mesh domain (warning: 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 course 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 will also be provided.
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.
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 described 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 gij to be:
where each of the partial derivatives comes from the definition of the covariant base vectors. Each of these base bectors 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.
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. Pij 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 gij 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.
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 (x1, y1), (x2, y2) 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 analytical 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.
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. Upon implementation, the Winslow system becomes the Poisson system
where f1 and f2 are the stretching functions of and respectively. 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. In the context of this program, 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.
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 height to the length of the cell.
- JMathPlot by Yann Richet
© Chaitanya Varier 2017. This software is protected under the MIT License.
I wrote this software independently for Prof. Jonathon Sumner at Dawson College, who provided me with invaluable mentorship and advice throughout the course of my research.
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.