# GMSH Mesh Expert to Julia using koehlerson/gmsh.jl
**URL:** https://discourse.julialang.org/t/gmsh-mesh-expert-to-julia-using-koehlerson-gmsh-jl/58169
**Category:** General Usage
**Tags:** question, package, gmsh
**Created:** [March 29, 2021, 4:14pm UTC](https://discourse.julialang.org/t/gmsh-mesh-expert-to-julia-using-koehlerson-gmsh-jl/58169 "2021-03-29T16:14:30Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@ziolai](https://discourse.julialang.org/u/ziolai)
#### Post date: [March 29, 2021, 4:14pm UTC](https://discourse.julialang.org/t/gmsh-mesh-expert-to-julia-using-koehlerson-gmsh-jl/58169/1 "2021-03-29T16:14:30Z")
**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 = x12_y13 - x13_y12; 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.
---
### Author:  [@vavrines](https://discourse.julialang.org/u/vavrines)
#### Post date: [March 29, 2021, 5:18pm UTC](https://discourse.julialang.org/t/gmsh-mesh-expert-to-julia-using-koehlerson-gmsh-jl/58169/2 "2021-03-29T17:18:02Z")
I guess you could use different mesh reader to verify your current code.
A possible choice is [meshio](https://github.com/nschloe/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.
---
### Author:  [@ziolai](https://discourse.julialang.org/u/ziolai)
#### Post date: [March 29, 2021, 5:51pm UTC](https://discourse.julialang.org/t/gmsh-mesh-expert-to-julia-using-koehlerson-gmsh-jl/58169/3 "2021-03-29T17:51:49Z")
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.
---
### Author:  [@ziolai](https://discourse.julialang.org/u/ziolai)
#### Post date: [March 29, 2021, 8:17pm UTC](https://discourse.julialang.org/t/gmsh-mesh-expert-to-julia-using-koehlerson-gmsh-jl/58169/4 "2021-03-29T20:17:13Z")
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
---
### Author:  [@ziolai](https://discourse.julialang.org/u/ziolai)
#### Post date: [March 30, 2021, 11:03am UTC](https://discourse.julialang.org/t/gmsh-mesh-expert-to-julia-using-koehlerson-gmsh-jl/58169/5 "2021-03-30T11:03:03Z")
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]