CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Special Topics > Mesh Generation & Pre-Processing

Gmsh.jl

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   March 29, 2021, 13:19
Default Gmsh.jl
  #1
Senior Member
 
Domenico Lahaye
Join Date: Dec 2013
Posts: 346
Blog Entries: 1
Rep Power: 10
dlahaye is on a distinguished road
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.
dlahaye is offline   Reply With Quote

Old   March 30, 2021, 09:30
Default
  #2
Senior Member
 
Domenico Lahaye
Join Date: Dec 2013
Posts: 346
Blog Entries: 1
Rep Power: 10
dlahaye is on a distinguished road
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
dlahaye is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On



All times are GMT -4. The time now is 03:22.