Mesh and Normals

[JarosT55]

Hi

'Mesh and Normals' … I am referring here to normal vectors of mesh's polygonal faces (in the case of skfem: quads and triangles)
I might have missed some skfem functionality of controlling mesh faces normal vectors direction … or might not?
The meshes (Quad and Tri's) normals are not unified. ***.
It is not so important in the context of the main skfem functionality, which is concentrated on solving linear equations in the form $Ab=x$
Changing the order of how meshe's triangles are enumerated ($mesh.t[:,i]$) does not affect the equation to be solved nor the result ($x$).
However, it does affect how a surface object (e.g. STL), the topological representation of a mesh, is constructed.
I have prepared a basic example file ( ex02_unify_normals - included ) showing the issue … once again, which is out of the main scope of the skfem package.
Nevertheless, although it has minor significance now, …. the feature of 'unified normals' could be had in the future.
Such functionality can be implemented 'externally' to skfem classes … but it would be more structural to have it 'inside'… and check the orientation of successive triangles (quads) in the run time, as the process should not be that complex. Kwarging of the functionality would ensure the backward functionality is intact.
What about other meshes, in particular 3D ones like MeshTet's? The answer is in the file!

Regards
MichaelT

` # ex02_unify_normals

from skfem import *
from skfem.helpers import dot, grad
import meshio
from skfem.io.meshio import to_meshio, from_meshio

# create the mesh
m = MeshTri().refined(6)
# or, with your own points and cells:
# m = MeshTri(points, cells)

e = ElementTriP1()
basis = Basis(m, e)

# this method could also be imported from skfem.models.laplace
@BilinearForm
def laplace(u, v, _):
return dot(grad(u), grad(v))

# this method could also be imported from skfem.models.unit_load
@LinearForm
def rhs(v, _):
return 1.0 * v

A = asm(laplace, basis)
b = asm(rhs, basis)
# enforce Dirichlet boundary conditions
A, b = enforce(A, b, D=m.boundary_nodes())
# solve -- can be anything that takes a sparse matrix and a right-hand side
x = solve(A, b)

def visualize():
from skfem.visuals.matplotlib import plot
return plot(m, x, shading='gouraud', colorbar=True)

visualize().show()

# #####################  Mesh to STL test section  ##########################

def unify_normals(m):
""" Unify normal directions of mesh cells. WARNING: input argument m is modified """
def crossproduct(u, v):
    return u[0]*v[1] - u[1]*v[0]

for e in range(m.nelements):
    if  crossproduct(m.p[:,m.t[1,e]] - m.p[:,m.t[0,e]], m.p[:,m.t[2,e]] - m.p[:,m.t[0,e]]) < 0.0:
        m.t[:,e] = m.t[::-1,e]
return m                                        # not required, but....

# Define your path_results, possibly with a dedicated folder
path_results = '?????'
subfolder = '????\\'
fpath     = os.path.join(path_results, subfolder )  
f_prefix  = 'unify_normals_Tri' 

# TEST A - MeshTri as they are created
testName  = 'no'
fpath_STL = os.path.join(fpath, f'{f_prefix}_{testName}.stl' ) 
m.save(fpath_STL)  

# TEST B - MeshTri after unification of 'normals' of triangles
testName  = 'yes'
fpath_STL = os.path.join(fpath, f'{f_prefix}_{testName}.stl' ) 
m_n = MeshTri(m.p, m.t)
m_n = unify_normals(m_n)
m_n.save(fpath_STL)  

# Compare created STLs e.g. in Meshmixer or other STL viewers showing 'normals'
# And assemble/resolve for m_n mesh
basis_n = Basis(m_n, e)
A_n = asm(laplace, basis_n)
b_n = asm(rhs, basis_n)
A_n, b_n = enforce(A_n, b, D=m_n.boundary_nodes())
x_n = solve(A_n, b_n)
visualize().show()

print( 'A=A_n =', np.allclose(A.toarray(), A_n.toarray(), rtol=1e-05, atol=1e-08))
print( 'b=b_n =', np.allclose(b,           b_n, rtol=1e-05, atol=1e-08))
print( 'x=x_n =', np.allclose(x,           x_n, rtol=1e-05, atol=1e-08))

# TEST C - Checking up what is going with MeshTet_s
testName  = 'no_corrupted_STL'
f_prefix  = 'unify_normals_Tet' 
fpath_STL = os.path.join(fpath, f'{f_prefix}_{testName}.stl' ) 
m_t = MeshTet.init_tensor(np.linspace(0, 1, 5), np.linspace(0, 1, 5), np.linspace(0, 1, 5))
m_t.save(fpath_STL)  

testName  = 'no_also_corrupted'
fpath_STL = os.path.join(fpath, f'{f_prefix}_{testName}.stl' ) 
meshio.write(fpath_STL, to_meshio(m_t))    

# EOF`

[kinnala]

There are several different reasons to reorder indices when implementing the finite element method. Usually it is related to edge orientation which can be communicated implicitly through the ordering of indices. The edge orientation can be used when implementing some elements - in fact, we have many elements that assume a specific order. The edge orientation is also used in some refinement algorithms.

This implementation technique has the unfortunate downside that the ordering is not suitable for the implicit communication of normal vectors. We have a completely different scheme for normals within scikit-fem which is not based on the ordering of the nodes.

I do not have plans to change this, it would be rather nontrivial task. If we want to save stl-files or some other files with normal vector information, it should be done differently, reordering during export.

[JarosT55]

Hi Mr Kinnala,

“I do not have plans to change this, it would be rather nontrivial task.”
Fair enough, … certainly, even a simple change in as complex software as skfem … is not trivial.

“We have a completely different scheme for normals within scikit-fem, which is not based on the ordering of the nodes.”
I understand this. Immersing homotopicaly $R^n$ space in vector/functional ($Sobolev$) space with their own norms and metrics $(L^p-norms$) does not guarantee that ‘normal’ in $R^n$ is the same as in what we call in skfem ElementSpace(). Element()s define their own metrics/norms, which can be pretty diverse. Thus, a ‘normal’ in $R^n$ is generally not the same as a ‘normal’ in space, called $H(R^n)$.

“This implementation technique has the unfortunate downside that the ordering is not suitable for the implicit communication of normal vectors.”
I do not know what you mean, particularly in the context of your following sentence. If the mapping of $R^n$ to $H(R^n)$ is homomorphic, then by definition … order(ing) of nodes in Rn does not affect the ‘shape’ of the final space (in the case of skfem, a assembly object).
However, as you clarified out …

“There are several different reasons to reorder indices when implementing the finite element method. Usually it is related to edge orientation which can be communicated implicitly through the ordering of indices. The edge orientation can be used when implementing some elements - in fact, we have many elements that assume a specific order. The edge orientation is also used in some refinement algorithms.”
Doing 'puristic mathematics’ exclusively in $H(R^n)$ is not always ‘ergonomic’, and some ‘shortcuts’ might be necessary (like, "... communicated implicitly through the ordering of indices"). I recognize this, … you certainly know better what is under skfem hood!

Regards
MichaelT
PS. … still, it is a pity that you have decided to leave $R^n$ space (the agnostic, innocent and having so many friends, standing patiently in line before Assembly() immersion) … unorganized!