Gmsh.jl
 

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.

edit: issue resolved by reordering nodes prior to loop over elements.

#..1/ Define geometry and mesh

gmsh.initialize()
#..define the model
model = gmsh.model
model.add("Square")
#..define mesh density near a point
cl = 0.1;
#..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()
# gmsh.option.setNumber("Mesh.MeshSizeMin", 5)
# gmsh.option.setNumber("Mesh.MeshSizeMax", 10)
#..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()

#..2/ Get and sort the mesh nodes
#..Observe that although the mesh is two-dimensional,
#..the z-coordinate that is equal to zero is stored as well.
#..Observe that the coordinates are stored contiguously for computational
#..efficiency

node_ids, node_coord, _ = gmsh.model.mesh.getNodes()
nnodes = length(node_ids)
#..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]

#..3/ Plotting the mesh

#..first alternative of two: plotting the nodes only
scatter(xnode, ynode)

#..second alternative of two: plotting the elements
z = ones(length(xnode))
trisurf(xnode,ynode,z)

#..4/ Get the mesh elements

element_types, element_ids, element_connectivity = gmsh.model.mesh.getElements(2)
nelements = length(element_ids[1])

#..5/ Perform a loop over the elements

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