Mesh from txt of coords and connectivity

I’m trying to generate a mesh from two txt files (nodes.txt and elements.txt ) using the tool FEM_TO_GMSH - Convert Triangulation to GMSH File . nodes.txt is a list of nodal coordinates, elements.txt is the connectivity.

However, the mesh consists of triangles embedded in 3d and this routine only works for solid elements in 3D.

Perhaps there is a simpler way of getting this done using dolfin.mesh directly? Thanks in advance for any suggestions

Have you tried to use mesheditor?

1 Like

Hi @dokken ,
Yes, the mesheditor is really easy… on the surface. Unfortunately, it reorders my nodes in a way that messes up the cell orientations (see attached). If I set close(order=False) I obviously get an UFC error. Is there any way to get around this?

from dolfin import *
import numpy as np

topological_dim = 2
geometrical_dim = 3

COORDS = np.loadtxt("nodes.txt", dtype='float')
CON = np.loadtxt("elements.txt",dtype='uint') - 1  # from MATLAB, start from 0

num_local_vertices = COORDS.shape[0]
num_global_vertices = num_local_vertices  # True if run in serial
num_local_cells = CON.shape[0]
num_global_cells = num_local_cells

# Create mesh object and open editor
mesh = df.Mesh()
editor = df.MeshEditor(), "triangle", topological_dim, geometrical_dim)
editor.init_vertices_global(num_local_vertices, num_global_vertices)
editor.init_cells_global(num_local_cells, num_global_cells)

# Add verticess
for i, coord in enumerate(COORDS):
    editor.add_vertex(i, coord)

# Add cells
for i, cell in enumerate(CON):
    editor.add_cell(i, cell)

# Close editor
# editor.close(order=False)  #  Error:   Unable to create mapping of degrees of freedom.
#*** Reason:  Mesh is not ordered according to the UFC numbering convention. Consider calling mesh.order().

f = df.File('mesh.pvd')
f << mesh
np.testing.assert_equal(mesh.cells()[-1], CON[-1])

Dolfin needs to order the mesh to be able to work with high order function spaces. Why not use order=True? What do you Need the ordering for?
As a slight tangent, you can see the effect of ordering a mesh here: Cell node labeling - #2 by dokken

Ok I see, the nodes were ordered (“upstream”) to produce a consistent outward oriented normal. I do something like this to get my (non-unit) normals:

J = Jacobian(mesh)
G1 = J[:,0]
G2 = J[:,1]
G3 = cross(G1, G2)

Is there another way to achieve this?

As your mesh looks like a sphere, you can easily use geometric information to adjust the normals.
Lets say the center of your sphere is at c[0], c[1], c[2].

J = Jacobian(mesh)
G1 = J[:,0]
G2 = J[:,1]
G3 = cross(G1, G2)
c = as_vector([[0,0,0]])
x = SpatialCoordinate(mesh)
r = x-c
n = G3*dot(G3, r)

gives you another non-unit normal.
If G3 is in outwards direction, dot(G3, r) is positive, otherwise it is negative.
You can then in turn normalize this vector.

Unfortunately I don’t think that works for me. For example if modeling a surface pressure load (psuedo code):

V = VectorFunctionSpace(mesh, 'CG', 1, dim=3)
u = Function(V)
v = TestFunction(V)
gsub1 = Gsub1 + u.dx(0)
gsub2 = Gsub2 + u.dx(1)
gsub3 = cross(gsub1, gsub2) # Deformed normal
f = dot(v, gsub3)*dx(mesh)

Then the deformed normal g_3 will not have the correct (consistent) orientation

(The code snippet assumes a geometric dimension =2 so doesn’t really apply in this case)