GMSH Mesh Expert to Julia using koehlerson/gmsh.jl

General Usage
, ,
1

Short Version of Query:

Consider a unit square (2D) domain and a triangular mesh on this geometry. I would like to perform a loop over all the elements in the mesh, compute the area of each triangle and verify that the sum of all areas equals one. I I am not sure that I am using connectivity between the elements and the nodes correctly. That is, I have doubt on the functions getNodes() and getElements() in Julia interface to gmsh correctly. I am obtaining elements with zero area.

Longer Version of Query:

using Plots

using Gmsh

gmsh.initialize()

#…define the model

model = gmsh.model

model.add(“Square”)

#…define mesh density near a point

cl = 0.5;

#…define four points in the geometry

gmsh.model.geo.addPoint(0.5, 0.5, 0., cl, 1)

gmsh.model.geo.addPoint(-0.5, 0.5, 0., cl, 2)

gmsh.model.geo.addPoint(-0.5, -0.5,0., cl, 3)

gmsh.model.geo.addPoint(0.5, -0.5, 0., cl, 4)

#…define four edges in the geometry

gmsh.model.geo.addLine(1, 2, 1)

gmsh.model.geo.addLine(2, 3, 2)

gmsh.model.geo.addLine(3, 4, 3)

gmsh.model.geo.addLine(4, 1, 4)

#…define outer boundary

gmsh.model.geo.addCurveLoop([1, 2, 3, 4], 1)

#…define planar surface

gmsh.model.geo.addPlaneSurface([1], 1)

#…define physics

gmsh.model.addPhysicalGroup(1, [1], 1)

#…synchronize the model

gmsh.model.geo.synchronize()

#…generate the mesh in 2D

model.mesh.generate(2)

#…save the mesh to file for future reference

gmsh.write(“square.msh”)

Finalize GMSH

#gmsh.finalize()

node_ids, node_coord, _ = gmsh.model.mesh.getNodes()

nnodes = length(node_ids)

xnode = node_coord[1:3:end]

ynode = node_coord[2:3:end]

#…Plot the mesh nodes

#…plots nodes only

scatter(xnode, ynode)

#…or alternatively

using GR

z = ones(length(xnode))

trisurf(xnode,ynode,z)

#…Retrieve the set of mesh nodes and the number of nodes (nnodes)

Observe that although the mesh is a two-dimensional mesh, the z-coordinate that is equal to zero is stored as well. Observe that the coordinates are stored contiguously

node_ids, node_coord, _ = gmsh.model.mesh.getNodes()

nnodes = length(node_ids)

xnode = node_coord[1:3:end]

ynode = node_coord[2:3:end]

#…Retrieve the set of mesh elements and the number of elements in the mesh…

element_types, element_ids, element_connectivity = gmsh.model.mesh.getElements(2,1)

nelements = length(element_ids[1])

for element_id in 1:nelements

node1_id = element_connectivity[1][3*(element_id-1)+1]

node2_id = element_connectivity[1][3*(element_id-1)+2]

node3_id = element_connectivity[1][3*(element_id-1)+3]

xnode1 = xnode[node1_id]; xnode2 = xnode[node2_id]; xnode3 = xnode[node3_id];

ynode1 = ynode[node1_id]; ynode2 = ynode[node2_id]; ynode3 = ynode[node3_id];

x12 = xnode2 - xnode1; x13 = xnode3-xnode1;

y12 = ynode2 - ynode1; y13 = ynode3-ynode1;

area_id = x12y13 - x13y12; area_id = abs(area_id)/2

println("on element ", element_id, " node-1 has global number ", node1_id)

println("on element ", element_id, " node-2 has global number ", node2_id)

println("on element ", element_id, " node-3 has global number ", node3_id)

println("on element ", element_id, " area = ", area_id)

println(" ")

end

Components Used:

I am using Julia Version 1.3.1 and koehlerson/gmsh.jl .

Thank you in advance for your input, Domenico Lahaye.

2

I guess you could use different mesh reader to verify your current code.
A possible choice is meshio that can be called via PyCall in Julia.
If the only task is to calculate areas, no complex connectivity information is needed. You just need to know the locations of nodes and their affiliations to cell ids.

3

Thanks! My understanding is that affiliation of nodes to cell ids is encoded in connectivity. I will see whether the mesh.cell field after mesh import gives me this information.

4

Any ideas why the following happens? Sincere thanks for the additional input.

Generating mesh using gmsh.jl

julia> #…generate the mesh in 2D
model.mesh.generate(2)
Info : Meshing 1D…
Info : [ 0%] Meshing curve 1 (Line)
Info : [ 30%] Meshing curve 2 (Line)
Info : [ 50%] Meshing curve 3 (Line)
Info : [ 80%] Meshing curve 4 (Line)
Info : Done meshing 1D (Wall 0.00127763s, CPU 0.000458s)
Info : Meshing 2D…
Info : Meshing surface 1 (Plane, Frontal-Delaunay)
Info : Done meshing 2D (Wall 0.00805564s, CPU 0.000929s)
Info : 12 nodes 26 elements

julia> #…save the mesh to file for future reference
gmsh.write(“square.msh”)
Info : Writing ‘square.msh’…
Info : Done writing ‘square.msh’

Checking mesh using meshio-info

ziolai$ meshio-info square.msh

Number of points: 3
Number of cells:
line: 2
Cell sets: gmsh:bounding_entities
Point data: gmsh:dim_tags
Cell data: gmsh:physical, gmsh:geometrical

5

Solved!

Below code that sort the node list prior to loop over elements. This sorting does the trick.

#…sort the node coordinates by ID, such that Node one sits at row 1

tosort = [node_ids node_coord[1:3:end] node_coord[2:3:end]];

sorted = sortslices(tosort , dims = 1);

node_ids = sorted[:,1]

xnode = sorted[:,2]

ynode = sorted[:,3]