# How to calculate skewness for a mesh?

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);

• You can have mesh skewed (non-orthogonal) but the line connecting centers of adjacent cells would pass through the center of the face between those cells. Just imagine each cell is a rhombus of the same dimensions. – Maxim Umansky Aug 14 at 17:50
• @MaximUmansky But aren't skewness and non-orthogonality two different things? If the line connecting centers of adjacent cells passes through the center of the face, the mesh is non-skewed but may be non-orthogonal (depending on the angle between face normal vector and the line joining the cells centroids). – Algo Aug 14 at 18:33
• I thought skewness is the same as non-orthogonality. If you google "skewed grid" it would show images that to me look like non-orthogonal grids. But what is your definition of grid skewness? – Maxim Umansky Aug 14 at 19:05
• @MaximUmansky Well, this is actually my question. The book I referenced describes skewness as the line connecting the cells centers not passing through the face center, but no metric or whatsoever. – Algo Aug 14 at 19:13
• There are some metrics of skewness on en.wikipedia.org/wiki/Types_of_mesh#Skewness. I don't know what those are based on, if there are real references; but at least a starting point. – Maxim Umansky Aug 14 at 19:27

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$$: $$v\cdot 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. • Thanks, Abdullah, great code explaination. You mentioned different definitions of skewness, can you expand the answer and list some? – Algo Aug 16 at 9:49
• @Algo, thanks a lot for the compliments. I expanded my answer to include a small discussion on skewness definitions. – Abdullah Ali Sivas Aug 16 at 21:05
• Can't upvote enough, many thanks. – Algo Aug 16 at 21:55
• @Algo, I don't know any references. I have it in my lecture notes from years ago but I do not have a cited paper there. When googled, I found that Ansys implements it as ||c-d||/|F|\$, so maybe you can refer to that. Probably, it is just considered common knowledge and noone knows the original paper. – Abdullah Ali Sivas Aug 17 at 19:03
• Looking for a reference, I found this "Dose, B. (2013). CFD Simulations of a 2.5 MW wind turbine using ANSYS CFX and OpenFOAM (Doctoral dissertation, MSc Thesis, UAS Kiel and FhG IWES, Germany)." which has a different metric for skewness. They attribute their definition to OpenFOAM but it is definitely not what OpenFOAM implements according to the code above. – Abdullah Ali Sivas Aug 17 at 19:11