refactor(io): remove bundled Gmsh reader and document I/O policy

The core package keeps the mesh and assembly stack type-stable without
shipping a Gmsh dependency or an always-loaded ASCII parser, and the IO
README now states where format readers should live going forward.

- Delete `src/io/gmsh_reader.jl` (`GmshReader`, `read_gmsh_mesh`, `GmshMesh`) including its Dict-heavy parsing path
- Rewrite `src/io/README.md` to describe policy-only `src/io/` (weakdeps/extensions or separate packages), point legacy `.inp`/`.med` readers at `src/legacy/io/` behind `JULIAFEM_ENABLE_LEGACY=1`, and warn against unconditional `include` growth for heavy I/O

src/io/README.md

@@ -1,24 +1,19 @@
 # src/io/
 
-Mesh readers that ship with the current 0.x package surface.
+Policy for **external mesh formats** (Gmsh, Netgen, Abaqus, Code Aster, …).
 
-The Abaqus `.inp` and Code Aster `.med` readers were built around the
+The core package owns the type-stable mesh model (`Mesh{N, T}` under
-older `Element(Poi1, ...)` constructors and the Dict-based field
+`src/mesh/`) and assembly; it does **not** ship format-specific parsers in
-system; they have been moved into the optional `JuliaFEM.Legacy`
+`src/io/`. Converters that read third-party files and build `Mesh{…}`,
-submodule (see `src/legacy/io/`). Set the environment variable
+`element_sets`, and `node_sets` should live in **separate packages** or in
-`JULIAFEM_ENABLE_LEGACY=1` before `using JuliaFEM` to load them.
+**JuliaFEM package extensions** wired through `[weakdeps]` / `[extensions]` in
+`Project.toml`, mirroring the MPI extension pattern (`JuliaFEMMPIExt`).
 
-## Files
+Legacy `.inp` / `.med` readers tied to the pre-reset API remain under
+`src/legacy/io/` and load only with `JULIAFEM_ENABLE_LEGACY=1`.
 
-- `gmsh_reader.jl`
+## This directory
-  Self-contained Gmsh `.msh` (ASCII format 4.1) reader.
-  Defines its own `JuliaFEM.GmshReader` submodule with `GmshMesh`,
-  `read_gmsh_mesh`. Has no dependency on legacy types and is loaded
-  unconditionally.
 
-## Future work
+This folder holds **documentation only** until a new extension or companion
-
+package is added. Do not grow unconditional `include("io/*.jl")` paths in
-A current-API replacement for the Abaqus reader (returning
+`JuliaFEM.jl` for heavy or format-specific I/O.
-`Mesh{T<:AbstractTopology, N}` and elements built via `create_elements!`)
-is on the roadmap. When it lands it will live here next to
-`gmsh_reader.jl`.

src/io/gmsh_reader.jl

@@ -1,173 +0,0 @@
-"""
-Simple Gmsh mesh reader for tetrahedral meshes
-
-Reads .msh format (ASCII version 4.1) and extracts:
-- Nodes (coordinates)
-- Elements (Tet4 connectivity)
-- Physical groups (for boundary conditions)
-"""
-
-module GmshReader
-
-export read_gmsh_mesh, GmshMesh
-
-using LinearAlgebra
-
-"""
-Mesh data structure
-"""
-struct GmshMesh
-    nodes::Matrix{Float64}              # 3 × n_nodes
-    elements::Matrix{Int}               # 4 × n_elems (Tet4)
-    physical_groups::Dict{String,Vector{Int}}  # Name → element/node indices
-end
-
-"""
-Read Gmsh mesh file (ASCII format 4.1)
-"""
-function read_gmsh_mesh(filename::String)
-    @info "Reading Gmsh mesh" filename
-
-    # Storage
-    nodes = Dict{Int,Vector{Float64}}()
-    elements = Vector{Tuple{Int,Int,Int,Int}}()
-    physical_names = Dict{Int,String}()
-    element_groups = Dict{Int,Vector{Int}}()
-
-    open(filename, "r") do io
-        while !eof(io)
-            line = strip(readline(io))
-
-            # Parse physical names
-            if line == "\$PhysicalNames"
-                n_names = parse(Int, readline(io))
-                for _ in 1:n_names
-                    parts = split(readline(io))
-                    dim = parse(Int, parts[1])
-                    tag = parse(Int, parts[2])
-                    name = strip(parts[3], '"')
-                    physical_names[tag] = name
-                end
-                readline(io)  # \$EndPhysicalNames
-            end
-
-            # Parse nodes (format 4.1)
-            if line == "\$Nodes"
-                header = split(readline(io))
-                numEntityBlocks = parse(Int, header[1])
-                numNodes = parse(Int, header[2])
-
-                for _ in 1:numEntityBlocks
-                    # Entity block header: entityDim entityTag parametric numNodesInBlock
-                    block_header = split(readline(io))
-                    numNodesInBlock = parse(Int, block_header[4])
-
-                    # Read node IDs
-                    node_ids = Int[]
-                    for _ in 1:numNodesInBlock
-                        push!(node_ids, parse(Int, readline(io)))
-                    end
-
-                    # Read coordinates
-                    for node_id in node_ids
-                        coords = split(readline(io))
-                        x = parse(Float64, coords[1])
-                        y = parse(Float64, coords[2])
-                        z = parse(Float64, coords[3])
-                        nodes[node_id] = [x, y, z]
-                    end
-                end
-                readline(io)  # \$EndNodes
-            end
-
-            # Parse elements (format 4.1)
-            if line == "\$Elements"
-                header = split(readline(io))
-                numEntityBlocks = parse(Int, header[1])
-                numElements = parse(Int, header[2])
-
-                for _ in 1:numEntityBlocks
-                    # Entity block header: entityDim entityTag elementType numElementsInBlock
-                    block_header = split(readline(io))
-                    entityTag = parse(Int, block_header[2])
-                    elementType = parse(Int, block_header[3])
-                    numElementsInBlock = parse(Int, block_header[4])
-
-                    for _ in 1:numElementsInBlock
-                        parts = split(readline(io))
-                        elem_id = parse(Int, parts[1])
-
-                        # Tet4 element (type 4)
-                        if elementType == 4
-                            n1 = parse(Int, parts[2])
-                            n2 = parse(Int, parts[3])
-                            n3 = parse(Int, parts[4])
-                            n4 = parse(Int, parts[5])
-                            push!(elements, (n1, n2, n3, n4))
-
-                            # Track physical group (entity tag is the physical group)
-                            if entityTag > 0
-                                if !haskey(element_groups, entityTag)
-                                    element_groups[entityTag] = Int[]
-                                end
-                                push!(element_groups[entityTag], length(elements))
-                            end
-                        end
-                    end
-                end
-                readline(io)  # \$EndElements
-            end
-        end
-    end
-
-    # Convert nodes to matrix (3 × n_nodes)
-    n_nodes = length(nodes)
-    node_matrix = zeros(3, n_nodes)
-    for (node_id, coords) in nodes
-        node_matrix[:, node_id] = coords
-    end
-
-    # Convert elements to matrix (4 × n_elems)
-    n_elems = length(elements)
-    elem_matrix = zeros(Int, 4, n_elems)
-    for (i, elem) in enumerate(elements)
-        elem_matrix[:, i] = [elem[1], elem[2], elem[3], elem[4]]
-    end
-
-    # Build physical groups with names
-    groups = Dict{String,Vector{Int}}()
-    for (tag, elem_ids) in element_groups
-        name = get(physical_names, tag, "Group_$tag")
-        groups[name] = elem_ids
-    end
-
-    mesh = GmshMesh(node_matrix, elem_matrix, groups)
-
-    @info "Gmsh mesh loaded" nodes = size(node_matrix, 2) elements = size(elem_matrix, 2) physical_groups = collect(keys(groups))
-
-    return mesh
-end
-
-"""
-Extract nodes on a boundary surface
-"""
-function get_surface_nodes(mesh::GmshMesh, surface_name::String)
-    # For now, return nodes where X ≈ 0 (fixed end) or Z ≈ max (top)
-    # This is a placeholder - proper implementation would track surface elements
-
-    if surface_name == "FixedEnd"
-        # Find nodes with X ≈ 0
-        x_min = minimum(mesh.nodes[1, :])
-        tol = 1e-6
-        return findall(abs.(mesh.nodes[1, :] .- x_min) .< tol)
-    elseif surface_name == "PressureSurface"
-        # Find nodes with Z ≈ max
-        z_max = maximum(mesh.nodes[3, :])
-        tol = 1e-6
-        return findall(abs.(mesh.nodes[3, :] .- z_max) .< tol)
-    else
-        return Int[]
-    end
-end
-
-end  # module