Skip to content

Commit

Permalink
Add more gmsh stuff
Browse files Browse the repository at this point in the history
* handle geo files
* allow for `simplexgrid("my.msh")` and `simplexgrid("my.geo")`
* allow for `write("my.msh",grid)
* file handling with gmsh initializes and finalizes
  • Loading branch information
j-fu committed Nov 28, 2023
1 parent 502b627 commit eea840c
Show file tree
Hide file tree
Showing 6 changed files with 205 additions and 74 deletions.
92 changes: 55 additions & 37 deletions ext/ExtendableGridsGmshExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -117,18 +117,29 @@ end
gmshfile_to_simplexgrid(filename::String, Tc, Ti)
````
This function just reads an msh file, and creates a gmsh.model and then calls the 'mod_to_simplexgrid' function
(This function has to be called with an initialized gmsh environment)
This function reads a .msh or a .geo file, and creates a `gmsh.model`
If it is a .geo file, `gmsh.model.mesh.generate()` is called.
Finally, it calls the 'mod_to_simplexgrid' function.
This function is called in 'simplexgrid_from_gmsh'
`Tc` is the type of coordinates, `Ti` is the index type.
The function initializes and finalized the `gmsh` module.
"""
function gmshfile_to_simplexgrid(filename::String; incomplete = false, Tc, Ti)
gmsh.initialize()
base, ext = splitext(filename)
gmsh.open(filename)
if ext == ".geo"
gmsh.model.mesh.generate()
end
if !incomplete
return mod_to_simplexgrid(gmsh.model, Tc, Ti)
grid = mod_to_simplexgrid(gmsh.model, Tc, Ti)
else
return incomplete_mod_to_simplexgrid(gmsh.model, Tc, Ti)
grid = incomplete_mod_to_simplexgrid(gmsh.model, Tc, Ti)
end
gmsh.finalize()
return grid
end

"""
Expand All @@ -137,13 +148,18 @@ gmshfile_to_mixedgrid(filename::String, Tc, Ti)
````
This function just reads an msh file, and creates a gmsh.model and then calls the 'mod_to_mixedgrid' function
(This function has to be called with an initialized gmsh environment)
This function is called in 'mixedgrid_from_gmsh'
`Tc` is the type of coordinates, `Ti` is the index type.
This function initalizes and finalized gmsh.
"""
function gmshfile_to_mixedgrid(filename::String, Tc, Ti)
gmsh.initialize()
gmsh.open(filename)
return mod_to_mixedgrid(gmsh.model, Tc, Ti)
grid = mod_to_mixedgrid(gmsh.model, Tc, Ti)
gmsh.finalize()
return grid
end

"""
Expand All @@ -152,11 +168,14 @@ function simplexgrid_to_gmshfile(grid::ExtendableGrid, filename::String)
````
This function takes a simplexgrid, uses 'grid_to_mod' to create a corresponding gmsh module
Then it writes the module to a file
(this function has to be called with an initialized gmsh environment)
This function is called in 'write_gmsh'
Then it writes the module to a file.
This function initalizes and finalized gmsh.
"""
function simplexgrid_to_gmshfile(grid::ExtendableGrid; filename::String = "")
gmsh.initialize()

if VERSION >= v"1.9"
# This possibility is new in 1.9, see
# https://github.com/JuliaLang/julia/blob/release-1.9/NEWS.md#new-language-features
Expand All @@ -169,7 +188,7 @@ function simplexgrid_to_gmshfile(grid::ExtendableGrid; filename::String = "")
if filename != ""
gmsh.write(filename)
end
return mod
gmsh.finalize()
end

"""
Expand All @@ -179,12 +198,12 @@ mixedgrid_to_gmshfile(grid::ExtendableGrid, filename::String)
This function takes a mixed grid, uses 'grid_to_mod' to create a corresponding gmsh module
Then it writes the module to a file
(this function has to be called with an initialized gmsh environment)
This function is called in 'write_gmsh'
grid[CellNodes] must be a VariableTargetAdjacency structure
This function initializes and finalized gmsh.
"""
function mixedgrid_to_gmshfile(grid::ExtendableGrid; filename::String = "")
gmsh.initialize()
if VERSION >= v"1.9"
# This possibility is new in 1.9, see
# https://github.com/JuliaLang/julia/blob/release-1.9/NEWS.md#new-language-features
Expand All @@ -197,7 +216,7 @@ function mixedgrid_to_gmshfile(grid::ExtendableGrid; filename::String = "")
if filename != ""
gmsh.write(filename)
end
return mod
gmsh.finalize()
end

#---------------------------------------------------------------------------------------------
Expand Down Expand Up @@ -326,13 +345,13 @@ function mod_to_simplexgrid(model::Module, Tc, Ti)
cell_types, element_tags_cells, cell_node_tags = model.mesh.getElements(dim, -1)
face_types, element_tags_faces, face_node_tags = model.mesh.getElements(dim - 1, -1)

unique_cell_types=unique(cell_types)
tag_nodes=invperm(node_tags)
unique_cell_types = unique(cell_types)
tag_nodes = invperm(node_tags)

if length(unique_cell_types)>1
if length(unique_cell_types) > 1
@warn "mesh contains different cell types"
end

#check whether cells are tetrahedrons in 3d or triangles in 2d:
if dim == 3
if unique_cell_types[1] != 4
Expand All @@ -351,36 +370,36 @@ function mod_to_simplexgrid(model::Module, Tc, Ti)

#if dim=3, the coordinates (of the nodes) just have to be reshaped
#for dim=2, the z-coordinate has to be deleted
ncoord=Int(length(coords) / 3)
ncoord = Int(length(coords) / 3)
if dim == 3
coords_new = Tc.(reshape(coords, (3,ncoord)))
coords_new = Tc.(reshape(coords, (3, ncoord)))
else
coords_new = Tc.(reshape(coords, (3,ncoord))[1:2,:])
coords_new = Tc.(reshape(coords, (3, ncoord))[1:2, :])
end

#number of cells
ncells = Int(length(cell_node_tags[1])/(dim+1))
ncells = Int(length(cell_node_tags[1]) / (dim + 1))

#the nodes making up the cells is stored in "cell_node_tags",
#just in the wrong format and permuted
simplices=zeros(Ti,dim+1, ncells)
simplices = zeros(Ti, dim + 1, ncells)

for i in eachindex(simplices)
simplices[i]=tag_nodes[cell_node_tags[1][i]]
simplices[i] = tag_nodes[cell_node_tags[1][i]]
end

#the physicalnames are currently unused
cellregion_to_physicalname = Bijection{Ti, String}()
pgnum_to_physcialname = Dict()
cr_count = 1

#the cellregions correspond to the physical groups in which the cells are
cellregions = ones(Ti, ncells)
pgs_data=model.getPhysicalGroups(dim)
if length(pgs_data)>0

pgs_data = model.getPhysicalGroups(dim)
if length(pgs_data) > 0
pgs = take_second(pgs_data)

for pg in pgs
name = model.getPhysicalName(dim, pg)
if length(name) == 0
Expand All @@ -397,31 +416,30 @@ function mod_to_simplexgrid(model::Module, Tc, Ti)
if entitytag in model.getEntitiesForPhysicalGroup(dim, pg)
cellregions[i] = cellregion_to_physicalname(pgnum_to_physcialname[pg]) #pg
break
end
end
end
end
end
# assemble the boundary faces, just reads the faces stored in the msh file
# for incomplete boundaries, there will be a function to seal them

nfaces = length(element_tags_faces[1])
bfaces=zeros(Ti,dim, nfaces)
bfaces = zeros(Ti, dim, nfaces)

for i in eachindex(bfaces)
bfaces[i]=tag_nodes[face_node_tags[1][i]]
bfaces[i] = tag_nodes[face_node_tags[1][i]]
end


# physical groups for bfaces
bfaceregions = ones(Ti, nfaces)
pgs_data=model.getPhysicalGroups(dim - 1)
if length(pgs_data)>0
pgs_data = model.getPhysicalGroups(dim - 1)
if length(pgs_data) > 0
pgs = take_second(pgs_data)

bfaceregion_to_physicalname = Bijection{Ti, String}()
pgnum_to_physcialname = Dict()
fr_count = 1

for pg in pgs
name = model.getPhysicalName(dim - 1, pg)
if length(name) == 0
Expand All @@ -431,7 +449,7 @@ function mod_to_simplexgrid(model::Module, Tc, Ti)
bfaceregion_to_physicalname[fr_count] = name
fr_count += 1
end

for i = 1:nfaces
_, _, _, entitytag = model.mesh.getElement(element_tags_faces[1][i])
for pg in pgs
Expand Down
22 changes: 20 additions & 2 deletions src/io.jl
Original file line number Diff line number Diff line change
Expand Up @@ -42,13 +42,21 @@ end
"""
$(TYPEDSIGNATURES)
Write grid to file. Currently for pdelib sg format only
Write grid to file. Currently for pdelib sg and Gmsh formats.
"""
function Base.write(fname::String, g::ExtendableGrid; format = "")
(fbase, fext) = splitext(fname)
if format == ""
format = fext[2:end]
end
if format == "msh"
try
simplexgrid_to_gmsh(g; filename = fname)
catch e
throw(ErrorException("Missing Gmsh extension. Add Gmsh.jl to your environment and import it to write $(fext) files."))
end
return
end
@assert format == "sg"

dim_g = dim_grid(g)
Expand Down Expand Up @@ -107,14 +115,24 @@ end
"""
$(TYPEDSIGNATURES)
Read grid from file. Currently for pdelib sg format only.
Read grid from file. Currently for pdelib sg and Gmsh formats.
"""
function simplexgrid(file::String; format = "")
Ti = Cint
(fbase, fext) = splitext(file)
if format == ""
format = fext[2:end]
end
if format == "msh" || format == "geo"
grid = nothing
try
grid = simplexgrid_from_gmsh(file)
catch e
throw(ErrorException("Missing Gmsh extension. Add Gmsh.jl to your environment and import it to read $(fext) files."))
end
return grid
end

@assert format == "sg"

tks = TokenStream(file)
Expand Down
74 changes: 74 additions & 0 deletions test/cube6.geo
Original file line number Diff line number Diff line change
@@ -0,0 +1,74 @@
//
// Copied and distributed with permission from
// https://www.math.univ-paris13.fr/~cuvelier/software/gmshgeo.html
//

DefineConstant[
N = {10, Name "Input/1Points "}
];
L=1;
h=L/N;

Point (1) = {0, 0, 0, h};
Point (2) = {L, 0, 0, h};
Point (3) = {L, L, 0, h};
Point (4) = {0, L, 0, h};
Point (5) = {0, 0, L, h};
Point (6) = {L, 0, L, h};
Point (7) = {L, L, L, h};
Point (8) = {0, L, L, h};
Line (1) = {1, 2};
Line (2) = {2, 3};
Line (3) = {3, 4};
Line (4) = {4, 1};
Line (5) = {5, 6};
Line (6) = {6, 7};
Line (7) = {7, 8};
Line (8) = {8, 5};
Line (9) = {1, 5};
Line (10) = {2, 6};
Line (11) = {3, 7};
Line (12) = {4, 8};

Line Loop(13) = {3, 12, -7, -11};
Plane Surface(4) = {13};
Line Loop(15) = {4, 9, -8, -12};
Plane Surface(1) = {15};
Line Loop(17) = {5, -10, -1, 9};
Plane Surface(3) = {17};
Line Loop(19) = {3, 4, 1, 2};
Plane Surface(5) = {19};
Line Loop(21) = {11, -6, -10, 2};
Plane Surface(2) = {21};
Line Loop(23) = {7, 8, 5, 6};
Plane Surface(6) = {23};
Surface Loop(25) = {4, 5, 1, 3, 6, 2};
Volume(1) = {25};

Physical Line(1) = {1};
Physical Line(2) = {2};
Physical Line(3) = {3};
Physical Line(4) = {4};
Physical Line(5) = {5};
Physical Line(6) = {6};
Physical Line(7) = {7};
Physical Line(8) = {8};
Physical Line(9) = {9};
Physical Line(10) = {10};
Physical Line(11) = {11};
Physical Line(12) = {12};
Physical Surface(1) = {1};
Physical Surface(4) = {4};
Physical Surface(5) = {5};
Physical Surface(2) = {2};
Physical Surface(6) = {6};
Physical Surface(3) = {3};
Physical Volume(1) = {1};
Physical Point(101) = {1};
Physical Point(102) = {2};
Physical Point(103) = {3};
Physical Point(104) = {4};
Physical Point(201) = {5};
Physical Point(202) = {6};
Physical Point(203) = {7};
Physical Point(204) = {8};
34 changes: 34 additions & 0 deletions test/disk1hole.geo
Original file line number Diff line number Diff line change
@@ -0,0 +1,34 @@
//
// Copied and distributed with permission from
// https://www.math.univ-paris13.fr/~cuvelier/software/gmshgeo.html
//

N=20;
h = 1/N;
R=1.;
r=0.45;
// Center Points of big circle
Point(1) = {0, 0, 0, h};
Point(2) = {R, 0, 0, h};
Point(3) = {-R, 0, 0, h/2};
Point(4) = {0,R, 0, h};
Point(5) = {0,-R, 0, h};

Point(6) = {-r, 0, 0, h}; // center
Point(7) = {-r, r, 0, h};
Point(8) = {-r, -r, 0, h};
Point(9) = {-2*r, 0, 0, h/2};
Circle(1) = {2, 1, 4};
Circle(2) = {4, 1, 3};
Circle(3) = {3, 1, 5};
Circle(4) = {5, 1, 2};
Circle(5) = {1, 6, 7};
Circle(6) = {7, 6, 9};
Circle(7) = {9, 6, 8};
Circle(8) = {8, 6, 1};
Line Loop(9) = {1, 2, 3, 4};
Line Loop(10) = {5, 6, 7, 8};
Plane Surface(11) = {9, 10};
Physical Line(1) = {1, 2, 3, 4};
Physical Line(2) = {5, 6, 7, 8};
Physical Surface(1) = {11};
Loading

0 comments on commit eea840c

Please sign in to comment.