diff --git a/Project.toml b/Project.toml
index 815fcfe9..9416e427 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "FinEtools"
 uuid = "91bb5406-6c9a-523d-811d-0644c4229550"
 authors = ["Petr Krysl <pkrysl@ucsd.edu>"]
-version = "8.0.0"
+version = "8.0.1"
 
 [deps]
 ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"
diff --git a/README.md b/README.md
index fe4d0a8b..489b58d8 100644
--- a/README.md
+++ b/README.md
@@ -25,6 +25,7 @@ The package supports application packages, for instance:
 
 ## News
 
+- 02/19/2024: Implement mesh reordering.
 - 02/17/2024: Implement generic parallel matrix assembly using threaded tasks.
 - 12/31/2023: Update for Julia 1.10.
 - 06/19/2023: Introduce DataCache, generic linear and bilinear forms. 
diff --git a/src/FinEtools.jl b/src/FinEtools.jl
index 4682f08f..74df334d 100644
--- a/src/FinEtools.jl
+++ b/src/FinEtools.jl
@@ -297,7 +297,8 @@ using .MeshModificationModule:
     pointpartitioning,
     interior2boundary,
     distortblock,
-    outer_surface_of_solid
+    outer_surface_of_solid,
+    reordermesh
 # Exported: extraction of boundary, fusing of nodes  and merging of meshes, mesh smoothing,  node partitioning
 export meshboundary,
     fusenodes,
@@ -312,7 +313,8 @@ export meshboundary,
     pointpartitioning,
     interior2boundary,
     distortblock,
-    outer_surface_of_solid
+    outer_surface_of_solid,
+    reordermesh
 
 using .MeshImportModule: import_NASTRAN, import_ABAQUS, import_MESH, import_H5MESH
 # Exported: mesh import functions
diff --git a/src/MeshModificationModule.jl b/src/MeshModificationModule.jl
index 878c4fd2..b3946af1 100644
--- a/src/MeshModificationModule.jl
+++ b/src/MeshModificationModule.jl
@@ -565,7 +565,7 @@ numbering for the nodes.
 `new_numbering` = new serial numbers for the nodes.  The connectivity
           should be changed as `conn[j]` --> `new_numbering[conn[j]]`
 
-Let us say there are nodes not connected to any finite element that you would
+Let us say there are nodes not connected to any finite element that we would
 like to remove from the mesh: here is how that would be accomplished.
 ```
 connected = findunconnnodes(fens, fes);
@@ -1238,4 +1238,17 @@ function outer_surface_of_solid(fens::FENodeSet, bdry_fes::ET) where {ET<:Abstra
     return subset(bdry_fes, osfesl)
 end
 
+"""
+    reordermesh(fens, fes, ordering)    
+
+Reorder mesh (reshuffle nodes, renumber connectivities correspondingly).
+
+The ordering may come from Reverse Cuthill-McKey (package SymRCM).
+"""
+function reordermesh(fens, fes, ordering)
+    iordering = collect(1:length(ordering))
+    iordering[ordering] = iordering
+    return FENodeSet(fens.xyz[ordering, :]), renumberconn!(fes, iordering)
+end
+
 end
diff --git a/test/test_meshing_2.jl b/test/test_meshing_2.jl
index a92bd50e..7301ec67 100644
--- a/test/test_meshing_2.jl
+++ b/test/test_meshing_2.jl
@@ -1846,3 +1846,87 @@ function test()
 end
 test()
 end
+
+
+module misc_reorder_1
+using FinEtools
+using SymRCM
+using LinearAlgebra
+using Test
+function test()
+    xs = collect(linearspace(0.0, 5 * pi / 2, 6))
+    ys = collect(linearspace(0.0, 10.0, 7))
+    zs = collect(linearspace(0.0, 10.0, 7))
+
+    for f in (H8blockx, H20blockx, )
+    # for f in (H8blockx, T4blockx, H20blockx, T10blockx)
+        fens, fes = f(xs, ys, zs)
+        geom = NodalField(fens.xyz)
+        femm = FEMMBase(IntegDomain(bfes, GaussRule(3, 3)))
+        V0 = integratefunction(femm, geom, (x) -> 1.0)
+
+        C = connectionmatrix(FEMMBase(IntegDomain(fes)), count(fens))
+            ordering = symrcm(C)
+        fens, fes = reordermesh(fens, fes, ordering)
+        geom = NodalField(fens.xyz)
+        femm = FEMMBase(IntegDomain(bfes, GaussRule(3, 3)))
+        V = integratefunction(femm, geom, (x) -> 1.0)
+        @test abs(V - V0) / V0 < 1.0e-6
+        
+    end
+    nothing
+end
+test()
+end
+
+module misc_reorder_1
+using FinEtools
+using SymRCM
+# using FinEtools.MeshExportModule: VTKWrite
+using LinearAlgebra
+using Test
+function test()
+    xs = collect(linearspace(0.0, 5 * pi / 2, 6))
+    ys = collect(linearspace(0.0, 10.0, 7))
+    zs = collect(linearspace(0.0, 10.0, 7))
+
+    for f in (H8blockx, H20blockx )
+    # for f in (H8blockx, T4blockx, H20blockx, T10blockx)
+        fens, fes = f(xs, ys, zs)
+        geom = NodalField(fens.xyz)
+        femm = FEMMBase(IntegDomain(fes, GaussRule(3, 3)))
+        V0 = integratefunction(femm, geom, (x) -> 1.0)
+
+        C = connectionmatrix(FEMMBase(IntegDomain(fes, GaussRule(3, 3))), count(fens))
+        ordering = symrcm(C)
+        fens, fes = reordermesh(fens, fes, ordering)
+        # File = "misc_reorder_1.vtk"
+        #  VTKWrite.vtkwrite(File, fens, fes)
+        geom = NodalField(fens.xyz)
+        femm = FEMMBase(IntegDomain(fes, GaussRule(3, 3)))
+        V = integratefunction(femm, geom, (x) -> 1.0)
+        @test abs(V - V0) / V0 < 1.0e-6
+    end
+
+    
+    # for f in (H8blockx, H20blockx )
+    for f in (T4blockx, T10blockx)
+        fens, fes = f(xs, ys, zs)
+        geom = NodalField(fens.xyz)
+        femm = FEMMBase(IntegDomain(fes, TetRule(4)))
+        V0 = integratefunction(femm, geom, (x) -> 1.0)
+
+        C = connectionmatrix(FEMMBase(IntegDomain(fes, TetRule(4))), count(fens))
+        ordering = symrcm(C)
+        fens, fes = reordermesh(fens, fes, ordering)
+        # File = "misc_reorder_1.vtk"
+        #  VTKWrite.vtkwrite(File, fens, fes)
+        geom = NodalField(fens.xyz)
+        femm = FEMMBase(IntegDomain(fes, TetRule(4)))
+        V = integratefunction(femm, geom, (x) -> 1.0)
+        @test abs(V - V0) / V0 < 1.0e-6
+    end
+    nothing
+end
+test()
+end
\ No newline at end of file