diff --git a/src/io/abaqus_reader.jl b/src/io/abaqus_reader.jl
new file mode 100644
index 0000000..9c066b0
--- /dev/null
+++ b/src/io/abaqus_reader.jl
@@ -0,0 +1,67 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+# AbaqusReader consolidated into src/readers/ - no longer a separate package
+
+"""
+    abaqus_read_mesh(fn::String)
+
+Read and parse ABAQUS `.inp` file.
+
+`fn` (filename) is the name of the file to parse.
+"""
+function abaqus_read_mesh(fn::String)
+    m = AbaqusReader.abaqus_read_mesh(fn)
+    return Mesh(m)
+end
+
+""" 
+    create_surface_elements(mesh::Mesh, surface_name::Symbol)
+
+Create a set of surface elements from solid elements.
+
+Notation follow what is defined in ABAQUS. For example, if solid elements
+are Tet10, surface elements will be Tri6 and they can be used to define
+boundary conditions.
+"""
+function create_surface_elements(mesh::Mesh, surface_name::Symbol)
+    elements = Element[]
+    for (elid, elsi) in mesh.surface_sets[surface_name]
+        elty = mesh.element_types[elid]
+        elco = mesh.elements[elid]
+        chel, chcon = AbaqusReader.create_surface_element(elty, elsi, elco)
+        ch = Element(getfield(JuliaFEM, chel), chcon)
+        push!(elements, ch)
+    end
+    update!(elements, "geometry", mesh.nodes)
+    return elements
+end
+
+""" 
+    create_surface_elements(mesh::Mesh, surface_name::String)
+
+Create a set of surface elements from solid elements.
+
+Notation follow what is defined in ABAQUS. For example, if solid elements
+are Tet10, surface elements will be Tri6 and they can be used to define
+boundary conditions.
+"""
+function create_surface_elements(mesh::Mesh, surface_name::String)
+    return create_surface_elements(mesh, Symbol(surface_name))
+end
+
+"""
+    create_nodal_elements(mesh::Mesh, node_set_name::String)
+
+Create a set of "nodal" elements from node set.
+
+Nodal elements are of type Poi1. In principle they don't have volume and it's
+not possible to integrate over them. It is however possible to add boundary
+conditions to them.
+"""
+function create_nodal_elements(mesh::Mesh, node_set_name::String)
+    node_ids = mesh.node_sets[Symbol(node_set_name)]
+    elements = [Element(Poi1, [j]) for j in node_ids]
+    update!(elements, "geometry", mesh.nodes)
+    return elements
+end
diff --git a/src/io/aster_reader.jl b/src/io/aster_reader.jl
new file mode 100644
index 0000000..04aef0b
--- /dev/null
+++ b/src/io/aster_reader.jl
@@ -0,0 +1,82 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+# using AsterReader  # Consolidated into src/readers/, functions are now directly available
+
+"""
+Nodes are different order in Code Aster compared to ABAQUS. This is yet
+incomplete mapping between the permutations. Most of this is still a
+great mystery.
+
+# References
+- http://onelab.info/pipermail/gmsh/2008/003850.html
+- http://caelinux.org/wiki/index.php/Proj:UNVConvert
+
+"""
+const med_connectivity = Dict{Symbol,Vector{Int}}(
+    :Tet4 => [4, 3, 1, 2],
+    :Tet10 => [4, 3, 1, 2, 10, 7, 8, 9, 6, 5],
+    :Pyr5 => [1, 4, 3, 2, 5],
+    :Wedge6 => [4, 5, 6, 1, 2, 3],
+    :Hex8 => [4, 8, 7, 3, 1, 5, 6, 2],
+    :Hex20 => [4, 8, 7, 3, 1, 5, 6, 2, 20, 15, 19, 11, 12, 16, 14, 10, 17, 13, 18, 9],
+    :Hex27 => [4, 8, 7, 3, 1, 5, 6, 2, 20, 15, 19, 11, 12, 16, 14, 10, 17, 13, 18, 9, 24, 25, 26, 23, 21, 22, 27])
+
+"""
+Map element names used in in Code Aster to element names used in in JuliaFEM
+"""
+const med_element_names = Dict{Symbol,Symbol}(
+    :PO1 => :Poi1,
+    :SE2 => :Seg2,
+    :SE3 => :Seg3,
+    :SE4 => :Seg4,
+    :TR3 => :Tri3,
+    :TR6 => :Tri6,
+    :TR7 => :Tri7,
+    :QU4 => :Quad4,
+    :QU8 => :Quad8,
+    :QU9 => :Quad9,
+    :TE4 => :Tet4,
+    :T10 => :Tet10,
+    :PE6 => :Wedge6,
+    :P15 => :Wedge15,
+    :P18 => :Wedge18,
+    :HE8 => :Hex8,
+    :H20 => :Hex20,
+    :H27 => :Hex27,
+    :PY5 => :Pyr5,
+    :P13 => :Pyr13)
+
+"""
+    aster_read_mesh(filename, mesh_name=nothing; reorder_element_connectivity=true)
+
+Read code aster mesh from file and return `Mesh` structure.
+
+If mesh file contains several meshes, a name of mesh must be given. By default,
+elements are reordered so that they match to the conventions used in JuliaFEM.
+"""
+function aster_read_mesh(filename::String, mesh_name=nothing; reorder_element_connectivity=true)
+    m = aster_read_mesh_raw(filename, mesh_name)
+    mesh = Mesh(m)
+    for (elid, eltype) in mesh.element_types
+        mesh.element_types[elid] = med_element_names[eltype]
+    end
+    if reorder_element_connectivity
+        reorder_element_connectivity!(mesh, med_connectivity)
+    end
+    nnodes = length(mesh.nodes)
+    nelements = length(mesh.elements)
+    @info("Mesh parsed from Code Aster file $filename.")
+    @info("Mesh contains $nnodes nodes and $nelements elements.")
+    for (elset_name, elset_elids) in mesh.element_sets
+        content = Dict{Symbol,Int}()
+        for elid in elset_elids
+            eltype = mesh.element_types[elid]
+            content[eltype] = get(content, eltype, 0) + 1
+        end
+        s = join(("$v x $k" for (k, v) in content), ", ")
+        nels = length(elset_elids)
+        @info("Element set $elset_name contains $nels elements ($s).")
+    end
+    return mesh
+end
diff --git a/src/io/io.jl b/src/io/io.jl
new file mode 100644
index 0000000..5e57684
--- /dev/null
+++ b/src/io/io.jl
@@ -0,0 +1,43 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+"""
+    JuliaFEM.IO
+
+Input/Output submodule for reading mesh files and writing results.
+
+# Available (zero external dependencies)
+- `abaqus_read_mesh(filename)`: Read ABAQUS .inp files
+- `create_surface_elements(mesh, name)`: Create surface elements
+- `create_nodal_elements(mesh, name)`: Create nodal elements
+
+# Future (requires optional dependencies)
+- Code Aster .med files: `aster_read_mesh` (requires HDF5.jl)
+- Xdmf output format: `Xdmf` writer (requires HDF5.jl + LightXML.jl)
+- VTK output: Native Julia implementation (no dependencies)
+
+# Usage
+```julia
+using JuliaFEM
+mesh = abaqus_read_mesh("model.inp")
+body = create_elements(mesh, "BODY")
+```
+"""
+module IO
+
+using ..JuliaFEM: Element, Mesh, add_node!, add_element_to_element_set!,
+    add_node_to_node_set!
+
+# Re-export for convenience
+import ..JuliaFEM.Preprocess
+using SparseArrays, LinearAlgebra
+using Logging
+
+include("abaqus_reader.jl")
+export abaqus_read_mesh, create_surface_elements, create_nodal_elements
+
+# AsterReader would go here when HDF5 is added as optional dependency
+# include("aster_reader.jl")
+# export aster_read_mesh
+
+end # module IO