From 6657ec8419248d4606bc09727861899bde6bdf89 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Mon, 11 May 2026 02:57:38 +0300 Subject: [PATCH] feat(ext): add Gmsh and WriteVTK package extensions Gmsh and WriteVTK remain optional weakdeps; these extensions implement the typed mesh import and VTU export entry points used by examples and tests. - Import linear Tris, Quads, Tets, Hexes from GMSH MSH with physical groups - Export Mesh to VTU with optional point and cell named tuples --- ext/JuliaFEMGmshExt.jl | 162 +++++++++++++++++++++++++++++++++++++ ext/JuliaFEMWriteVTKExt.jl | 119 +++++++++++++++++++++++++++ 2 files changed, 281 insertions(+) create mode 100644 ext/JuliaFEMGmshExt.jl create mode 100644 ext/JuliaFEMWriteVTKExt.jl diff --git a/ext/JuliaFEMGmshExt.jl b/ext/JuliaFEMGmshExt.jl new file mode 100644 index 0000000..a53093b --- /dev/null +++ b/ext/JuliaFEMGmshExt.jl @@ -0,0 +1,162 @@ +# SPDX-FileCopyrightText: 2015-2026 Jukka Aho +# SPDX-License-Identifier: MIT + +module JuliaFEMGmshExt + +using JuliaFEM +using Gmsh +using Tensors: Vec + +const _GMSH_LINEAR_TYPE = Dict{Int32,DataType}( + Int32(2) => Tri3, + Int32(3) => Quad4, + Int32(4) => Tet4, + Int32(5) => Hex8, +) + +function _physical_symbol(name::AbstractString, dim::Integer, tag::Integer)::Symbol + s = strip(String(name)) + if isempty(s) + return Symbol("physical_", dim, "_", tag) + end + symstr = replace(s, r"[^0-9a-zA-Z_]+" => "_") + if isempty(symstr) + return Symbol("physical_", dim, "_", tag) + end + return Symbol(symstr) +end + +function _mesh_dim!(gmsh, dim::Union{Nothing,Int})::Int32 + if dim === nothing + et3, _, _ = gmsh.model.mesh.getElements(3, -1) + if !isempty(et3) + return Int32(3) + end + et2, _, _ = gmsh.model.mesh.getElements(2, -1) + if !isempty(et2) + return Int32(2) + end + throw(ArgumentError("Gmsh model has no 2D or 3D mesh elements")) + end + dim == 2 || dim == 3 || throw(ArgumentError("dim must be 2 or 3, got $dim")) + return Int32(dim) +end + +function _topology_from_elem_types(elem_types::Vector{Int32})::DataType + isempty(elem_types) && throw(ArgumentError("Gmsh returned no element types for the requested dimension")) + Ts = DataType[] + for gtyp in elem_types + T = get(_GMSH_LINEAR_TYPE, gtyp, nothing) + T === nothing && throw( + ArgumentError( + "Unsupported Gmsh element type $gtyp (only linear types " * + "$(join(sort!(collect(keys(_GMSH_LINEAR_TYPE))), ", ")) are supported)", + ), + ) + push!(Ts, T) + end + all(==(Ts[1]), Ts) || + throw(ArgumentError("Mixed Gmsh element types in one dimension are not supported: $Ts")) + return Ts[1] +end + +function _import_mesh(gmsh; dim::Union{Nothing,Int}, quiet::Bool) + mesh_dim = _mesh_dim!(gmsh, dim) + elem_types, elem_tags, elem_node_tags = gmsh.model.mesh.getElements(mesh_dim, -1) + Ttop = _topology_from_elem_types(elem_types) + + node_tags, coord, _ = gmsh.model.mesh.getNodes() + length(node_tags) * 3 == length(coord) || + throw(ArgumentError("Unexpected Gmsh node coordinate layout")) + + n_nodes = length(node_tags) + tag_to_idx = Dict{UInt64,UInt32}() + nodes = Vector{Vec{3,Float64}}(undef, n_nodes) + @inbounds for i in 1:n_nodes + tag_to_idx[node_tags[i]] = UInt32(i) + base = 3 * (i - 1) + nodes[i] = Vec(coord[base + 1], coord[base + 2], coord[base + 3]) + end + + N = nnodes(Ttop()) + connectivity = NTuple{N,UInt32}[] + elem_tag_to_julia = Dict{UInt64,UInt32}() + @inbounds for it in eachindex(elem_types) + tags_i = elem_tags[it] + nt_i = elem_node_tags[it] + n_elem = length(tags_i) + expected_len = n_elem * N + length(nt_i) == expected_len || + throw(ArgumentError("Gmsh node tag list length mismatch for element type $(elem_types[it])")) + for e in 1:n_elem + julia_e = UInt32(length(connectivity) + 1) + elem_tag_to_julia[tags_i[e]] = julia_e + offs = N * (e - 1) + conn = ntuple(N) do k + tag = nt_i[offs + k] + idx = get(tag_to_idx, tag, nothing) + idx === nothing && throw(ArgumentError("Unknown node tag $tag in element connectivity")) + idx + end + push!(connectivity, conn) + end + end + + element_sets = Dict{Symbol,Set{UInt32}}(:all => Set(UInt32(1):UInt32(length(connectivity)))) + node_sets = Dict{Symbol,Set{UInt32}}(:all => Set(UInt32(1):UInt32(n_nodes))) + + for (pdim, ptag) in gmsh.model.getPhysicalGroups() + name = gmsh.model.getPhysicalName(Int(pdim), Int(ptag)) + sym = _physical_symbol(name, pdim, ptag) + if pdim == mesh_dim + acc = get!(Set{UInt32}, element_sets, sym) + for ent in gmsh.model.getEntitiesForPhysicalGroup(Int(pdim), Int(ptag)) + _, e_tags_b, _ = gmsh.model.mesh.getElements(Int(pdim), Int(ent)) + for it in eachindex(e_tags_b) + for etag in e_tags_b[it] + ji = get(elem_tag_to_julia, etag, nothing) + ji === nothing || push!(acc, ji) + end + end + end + elseif pdim == mesh_dim - 1 + acc = get!(Set{UInt32}, node_sets, sym) + for ent in gmsh.model.getEntitiesForPhysicalGroup(Int(pdim), Int(ptag)) + ntags_b, _, _ = gmsh.model.mesh.getNodes(Int(pdim), Int(ent)) + for t in ntags_b + ni = get(tag_to_idx, t, nothing) + ni === nothing || push!(acc, ni) + end + end + end + end + + for d in (element_sets, node_sets) + for k in collect(keys(d)) + k === :all && continue + isempty(d[k]) && delete!(d, k) + end + end + + return Mesh{N,Ttop}(nodes, connectivity, element_sets, node_sets) +end + +function JuliaFEM.read_gmsh_msh(path::AbstractString; dim=nothing, quiet::Bool=true) + argv = quiet ? String["-v", "0"] : String[] + started = Gmsh.initialize(argv; finalize_atexit=false) + try + Gmsh.gmsh.open(String(path)) + return _import_mesh(Gmsh.gmsh; dim, quiet) + finally + started && Gmsh.finalize() + end +end + +function JuliaFEM.mesh_from_current_gmsh_model(; dim=nothing, quiet::Bool=true) + Bool(Gmsh.gmsh.isInitialized()) || + throw(ArgumentError("Gmsh is not initialized; call Gmsh.initialize first")) + _ = quiet + return _import_mesh(Gmsh.gmsh; dim, quiet) +end + +end # module diff --git a/ext/JuliaFEMWriteVTKExt.jl b/ext/JuliaFEMWriteVTKExt.jl new file mode 100644 index 0000000..437558b --- /dev/null +++ b/ext/JuliaFEMWriteVTKExt.jl @@ -0,0 +1,119 @@ +# SPDX-FileCopyrightText: 2015-2026 Jukka Aho +# SPDX-License-Identifier: MIT + +module JuliaFEMWriteVTKExt + +using JuliaFEM +using WriteVTK + +const _VTK = WriteVTK.VTKCellTypes + +function _vtk_cell_type(::Type{T}) where {T} + if T === Hex8 + return _VTK.VTK_HEXAHEDRON + elseif T === Tet4 + return _VTK.VTK_TETRA + elseif T === Quad4 + return _VTK.VTK_QUAD + elseif T === Tri3 + return _VTK.VTK_TRIANGLE + elseif T === Seg2 + return _VTK.VTK_LINE + else + throw( + ArgumentError( + "VTK export does not support topology $T. " * + "Supported: Hex8, Tet4, Quad4, Tri3, Seg2.", + ), + ) + end +end + +function _points_matrix(mesh::Mesh{N,T}) where {N,T} + n = nnodes_total(mesh) + pts = Matrix{Float64}(undef, 3, n) + @inbounds for i in 1:n + v = mesh.nodes[i] + pts[1, i] = v[1] + pts[2, i] = v[2] + pts[3, i] = v[3] + end + return pts +end + +function _mesh_cells(mesh::Mesh{N,T}) where {N,T} + vtk_ct = _vtk_cell_type(T) + ne = nelements(mesh) + cells = Vector{MeshCell}(undef, ne) + @inbounds for e in 1:ne + conn = mesh.connectivity[e] + cells[e] = MeshCell(vtk_ct, collect(Int, conn)) + end + return cells +end + +function _strip_vtu_extension(basepath::AbstractString) + path = String(basepath) + if endswith(lowercase(path), ".vtu") + return path[1:(end - 4)] + end + return path +end + +function _validate_point_data(mesh::Mesh, point_data::NamedTuple) + n = nnodes_total(mesh) + for (k, v) in pairs(point_data) + if v isa AbstractVector + length(v) == n || + throw(ArgumentError("point_data.$k: expected length $n (nnodes), got $(length(v))")) + elseif v isa AbstractMatrix + size(v, 2) == n || + throw(ArgumentError("point_data.$k: expected matrix with $n columns, got $(size(v))")) + size(v, 1) in (1, 2, 3) || + throw(ArgumentError("point_data.$k: expected 1×n, 2×n, or 3×n matrix, got $(size(v))")) + else + throw(ArgumentError("point_data.$k: unsupported type $(typeof(v))")) + end + end + return nothing +end + +function _validate_cell_data(mesh::Mesh, cell_data::NamedTuple) + ne = nelements(mesh) + for (k, v) in pairs(cell_data) + v isa AbstractVector || + throw(ArgumentError("cell_data.$k: expected AbstractVector, got $(typeof(v))")) + length(v) == ne || + throw(ArgumentError("cell_data.$k: expected length $ne (nelements), got $(length(v))")) + end + return nothing +end + +function JuliaFEM.write_vtu_mesh( + basepath::AbstractString, + mesh::Mesh{N,T}; + point_data::NamedTuple = (;), + cell_data::NamedTuple = (;), +) where {N,T} + _vtk_cell_type(T) + root = _strip_vtu_extension(basepath) + pts = _points_matrix(mesh) + cells = _mesh_cells(mesh) + isempty(point_data) || _validate_point_data(mesh, point_data) + isempty(cell_data) || _validate_cell_data(mesh, cell_data) + + vtk = vtk_grid(root, pts, cells) + for (name, vals) in pairs(point_data) + if vals isa AbstractVector + vtk[String(name), VTKPointData()] = collect(Float64, vals) + else + vtk[String(name), VTKPointData()] = Float64.(vals) + end + end + for (name, vals) in pairs(cell_data) + vtk[String(name), VTKCellData()] = collect(Float64, vals) + end + return first(vtk_save(vtk)) +end + +end # module