How to calculate skewness for a mesh?

Question

I am writing a code to calculate mesh quality stats such as: cell volume, face areas and non-orthogonality between faces (basically something like OpenFOAM's checkMesh).
As per F. Moukalled et al, a mesh is skewed when the line connecting adjacent cells centroids does not pass through the centroid of the straddling face connecting the two cells. For example, if face centroid is denoted by $f$ and $f'$ is the intersection between the line connecting the two cells and the face, $f$ and $f'$ coincides for non-skewed meshes.
So, what is the metric to measure skewness?
I found the following code used in OpenFOAM to calculate skewness, but the math behind it is not very clear:
Note: the /* */ comments are mine, however, I am not 100% sure about my interpretation of the variables.
/* fCtrs[facei] is the face centroid of the current straddling face */
/* ownCc is the centroid of the cell that owns facei */
/* neiCc is the centroid of the neighbor cell */

vector Cpf = fCtrs[facei] - ownCc;
vector d = neiCc - ownCc;

// Skewness vector
/* the & operator is an overloaded operator that represents dot product */
/* ROOTVSMALL is a constant, equals "1.0e-18" (defined somewhere else), that prevent errors when dividing by zero */
/* fAreas[facei] returns the area normal vector of the straddling face */
vector sv =
    Cpf
    - ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + ROOTVSMALL))*d;

vector svHat = sv/(mag(sv) + ROOTVSMALL);

Abdullah Ali Sivas · Accepted Answer

From the discussions and the paper, OpenFOAM seems to have implemented a measure of skewness. This answer is not an explanation why the different definitions of skewness might be equivalent, I am just going to justify why this is a measure of skewness. Consider following two elements -for sake of simplicity-

Blue arrow is the outward surface normal fAreas[facei], red dots are, from left to right, ownCc, fCtrs[facei] and neiCc. Now, Cpf is the vector pointing from ownCc to fCtrs[facei] and d is the vector pointing from fCtrs[facei] to neiCc.
This is where I remind that, given two compatible vectors $v,w$: $$vcdot w = |v|  |w|  cos(theta)$$ where $theta$ is the angle between $v$ and $w$.
Let's get back to the formula ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + ROOTVSMALL)). (fAreas[facei] & Cpf) will just give us the norm of fAreas[facei] times Cpf as those two vectors point in the same direction (in this example, if Own was a trapezoid it would not be) thus $theta=0$. (fAreas[facei] & d) may give us a variety of different positive values, but the important points is if fAreas[facei] and d point in the same direction, hence no skewness, it will be the norm of fAreas[facei] times d, e.g. [norm(fAreas[facei])*norm(Cpf)]/[norm(fAreas[facei])*norm(d)] = norm(Cpf)/norm(d). This simplifies
sv = Cpf - ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + ROOTVSMALL))*d;

into
sv = Cpf - norm(Cpf)*d/norm(d); // Note that d/norm(d) is a unit vector pointing
                                // in the same direction as Cpf.

into
sv = Cpf - Cpf; // e.g. zero vector

Hence, if the mesh is not skewed, sv -and as a result svHat- will be zero. If it is skewed, as in the picture, the math is slightly different
sv = Cpf - ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + ROOTVSMALL))*d;

becomes
sv = 
 Cpf - 
  ((norm(fAreas[facei])*norm(Cpf))/(norm(fAreas[facei])*norm(d)*cos(theta) + ROOTVSMALL))*d;

which becomes (disregarding ROOTVSMALL)
sv = Cpf - (norm(Cpf)/(norm(d)*cos(theta) + ROOTVSMALL))*d;

with theta being the angle between d and fAreas[facei]. Let's reorganize (again I am disregarding ROOTVSMALL)
sv = Cpf - norm(Cpf)/norm(d)*d*(1/(cos(theta) + ROOTVSMALL));

This way it is more clear that how this is a measure of skewness. theta can take values in the open interval $(-pi/2,pi/2)$ for meshes without degenerate elements and 1/cos(theta) takes values in the interval $[1,infty)$. At the last step, there is the normalization svHat = sv/(mag(sv) + ROOTVSMALL); which generates a unit vector svHat and gives you the skewness in each direction. 0 means no skewness in the given direction and other values will signify some skewness. I think $-1$ would be the most skewed case and correspond to a degenerate neighbouring element.
Different measures of skewness
As Maxim Umansky mentioned in the comments to the question, there is a wikipedia article which discusses skewness. Those are valid measures of skewness of an element, however, they do not say anything about the skewness of the grid. Except for the one based on equilateral volume. For example, according to those measures meshing of a rhombus domain with rhombus elements would be considered skewed, however, that is not what you want.
Another skewness definition I am familiar with is $$1-frac{||c-d||}{|F|},$$ where $F$ is the face between two neighboring elements, $|F|$ is the area of the face, $c$ is the centroid of the face $F$ and $d$ is the midpoint of the line segment connecting the center of Own element to the center of neighbouring element. In this case, if $c$ and $d$ overlap for each couple of neighboring elements, it means the mesh is not skewed and you get a value of $1$. Hence, this definition of skewness is bounded above by $1$ but it can be an indefinitely large negative number.
Differences between this definition and OpenFOAM measure

The one I am familiar with gives you a scalar, OpenFOAM measure returns a vector and tells you the direction of skewness too
OpenFOAM measure is in $[-1,0]$ (if I am not wrong) and the other one is in $(-infty,1]$.
This one generalizes to polygons and polyhedra (this is a second hand information, i.e. something I heard some time ago), I am not sure about the OpenFOAM one.

Due to these reasons, even though I believe that they are equivalent definitions I cannot prove that they are, e.g. how do I compare a vector to a scalar? However, both would characterize the following two elements as highly skewed so that is my evidence towards their equivalence.

How to calculate skewness for a mesh?

One Answer

Add your own answers!

Ask a Question