diff --git a/apps/bingham/src/bingham_vector.cpp b/apps/bingham/src/bingham_vector.cpp
index 3b16ac23..1472e713 100644
--- a/apps/bingham/src/bingham_vector.cpp
+++ b/apps/bingham/src/bingham_vector.cpp
@@ -182,7 +182,8 @@ int main(int argc, char **argv)
     if (std::regex_match(input_mesh, std::regex(".*\\.medit2d$")))
     {
         std::cout << "Guessed mesh format: Medit format" << std::endl;
-        auto msh = disk::load_medit_2d_mesh<T>(input_mesh.c_str());
+        disk::generic_mesh<T, 2> msh;
+        disk::load_mesh_medit<T>(input_mesh.c_str(), msh);
         if(!run_bingham(msh, hdi, vp))
             std::cout << "No convergence" << std::endl;
 
@@ -193,7 +194,8 @@ int main(int argc, char **argv)
     if (std::regex_match(input_mesh, std::regex(".*\\.mesh2d$") ))
     {
         std::cout << "Guessed mesh format: Netgen 2D" << std::endl;
-        auto msh = disk::load_netgen_2d_mesh<T>(input_mesh.c_str());
+        disk::simplicial_mesh<T, 2> msh;
+        disk::load_mesh_netgen<T>(input_mesh.c_str(), msh);
         if(!run_bingham(msh, hdi, vp))
             std::cout << "No convergence" << std::endl;
 
@@ -204,7 +206,8 @@ int main(int argc, char **argv)
     if (std::regex_match(input_mesh, std::regex(".*\\.quad$") ))
     {
         std::cout << "Guessed mesh format: DiSk++ Cartesian 2D" << std::endl;
-        auto msh = disk::load_cartesian_2d_mesh<T>(input_mesh.c_str());
+        disk::cartesian_mesh<T, 2> msh;
+        disk::load_mesh_diskpp_cartesian<T>(input_mesh.c_str(), msh);
         if(!run_bingham(msh, hdi, vp))
             std::cout << "No convergence" << std::endl;
 
diff --git a/apps/bingham/src/bingham_vector_higher.cpp b/apps/bingham/src/bingham_vector_higher.cpp
index ce2bd5a7..fa470785 100644
--- a/apps/bingham/src/bingham_vector_higher.cpp
+++ b/apps/bingham/src/bingham_vector_higher.cpp
@@ -181,7 +181,9 @@ int main(int argc, char **argv)
     if (std::regex_match(input_mesh, std::regex(".*\\.medit2d$")))
     {
         std::cout << "Guessed mesh format: Medit format" << std::endl;
-        auto msh = disk::load_medit_2d_mesh<T>(input_mesh.c_str());
+        disk::generic_mesh<T, 2> msh;
+        disk::load_mesh_medit<T>(input_mesh.c_str(), msh);
+
         if(!run_bingham(msh, hdi, vp))
             std::cout << "No convergence" << std::endl;
 
@@ -192,7 +194,9 @@ int main(int argc, char **argv)
     if (std::regex_match(input_mesh, std::regex(".*\\.mesh2d$") ))
     {
         std::cout << "Guessed mesh format: Netgen 2D" << std::endl;
-        auto msh = disk::load_netgen_2d_mesh<T>(input_mesh.c_str());
+        disk::simplicial_mesh<T, 2> msh;
+        disk::load_mesh_netgen<T>(input_mesh.c_str(), msh);
+
         if(!run_bingham(msh, hdi, vp))
             std::cout << "No convergence" << std::endl;
 
@@ -203,7 +207,9 @@ int main(int argc, char **argv)
     if (std::regex_match(input_mesh, std::regex(".*\\.quad$") ))
     {
         std::cout << "Guessed mesh format: DiSk++ Cartesian 2D" << std::endl;
-        auto msh = disk::load_cartesian_2d_mesh<T>(input_mesh.c_str());
+        disk::cartesian_mesh<T, 2> msh;
+        disk::load_mesh_diskpp_cartesian<T>(input_mesh.c_str(), msh);
+
         if(!run_bingham(msh, hdi, vp))
             std::cout << "No convergence" << std::endl;
 
diff --git a/apps/contact/src/diffusion_nitsche.cpp b/apps/contact/src/diffusion_nitsche.cpp
index def48cc7..f0e009ab 100644
--- a/apps/contact/src/diffusion_nitsche.cpp
+++ b/apps/contact/src/diffusion_nitsche.cpp
@@ -100,7 +100,8 @@ int main(int argc, char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.typ1$") ))
     {
         std::cout << "Guessed mesh format: FVCA5 2D" << std::endl;
-        auto msh = disk::load_fvca5_2d_mesh<T>(mesh_filename);
+        disk::generic_mesh<T,2> msh;
+        disk::load_mesh_fvca5_2d<T>(mesh_filename, msh);
         auto error = run_diffusion_solver(msh, ap);
         return 0;
     }
@@ -109,7 +110,8 @@ int main(int argc, char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.mesh2d$") ))
     {
         std::cout << "Guessed mesh format: Netgen 2D" << std::endl;
-        auto msh = disk::load_netgen_2d_mesh<T>(mesh_filename);
+        disk::simplicial_mesh<T, 2> msh;
+        disk::load_mesh_netgen<T>(mesh_filename, msh);
         auto error = run_diffusion_solver(msh, ap);
         return 0;
     }
@@ -118,7 +120,8 @@ int main(int argc, char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.quad$") ))
     {
         std::cout << "Guessed mesh format: DiSk++ Cartesian 2D" << std::endl;
-        auto msh = disk::load_cartesian_2d_mesh<T>(mesh_filename);
+        disk::cartesian_mesh<T, 2> msh;
+        disk::load_mesh_diskpp_cartesian<T>(mesh_filename, msh);
         auto error = run_diffusion_solver(msh, ap);
         return 0;
     }
diff --git a/apps/contact/src/signorini_newton.cpp b/apps/contact/src/signorini_newton.cpp
index 77b8b8a0..4ee035b2 100644
--- a/apps/contact/src/signorini_newton.cpp
+++ b/apps/contact/src/signorini_newton.cpp
@@ -107,17 +107,9 @@ int main(int argc, char **argv)
     /* Netgen 2d*/
     if (std::regex_match(filename, std::regex(".*\\.mesh2d$") ))
     {
-        typedef disk::simplicial_mesh<T, 2>  mesh_type;
         std::cout << "Guessed mesh format: Netgen 2D" << std::endl;
-        mesh_type msh;
-        disk::netgen_mesh_loader<T, 2> loader;
-        if (!loader.read_mesh(filename))
-        {
-            std::cout << "Problem loading mesh." << std::endl;
-            return 1;
-        }
-        loader.populate_mesh(msh);
-
+        disk::simplicial_mesh<T, 2> msh;
+        disk::load_mesh_netgen<T>(filename, msh);
         run_signorini(msh, ap, run_exact);
         return 0;
     }
@@ -126,7 +118,8 @@ int main(int argc, char **argv)
     if (std::regex_match(filename, std::regex(".*\\.typ1$") ))
     {
         std::cout << "Guessed mesh format: FVCA5 2D" << std::endl;
-        auto msh = disk::load_fvca5_2d_mesh<T>(filename);
+        disk::generic_mesh<T,2> msh;
+        disk::load_mesh_fvca5_2d<T>(filename, msh);
         run_signorini(msh, ap, run_exact);
         return 0;
     }
@@ -135,22 +128,11 @@ int main(int argc, char **argv)
     if (std::regex_match(filename, std::regex(".*\\.quad$") ))
     {
         std::cout << "Guessed mesh format: DiSk++ Cartesian 2D" << std::endl;
-        auto msh = disk::load_cartesian_2d_mesh<T>(filename);
+        disk::cartesian_mesh<T, 2> msh;
+        disk::load_mesh_diskpp_cartesian<T>(filename, msh);
         run_signorini(msh, ap, run_exact);
         return 0;
     }
 
-
-    #if 0
-    /* Medit 2d*/
-    if (std::regex_match(filename, std::regex(".*\\.medit2d$")))
-    {
-        std::cout << "Guessed mesh format: Medit format" << std::endl;
-        typedef disk::generic_mesh<T, 2>  mesh_type;
-        mesh_type msh = disk::load_medit_2d_mesh<T>(filename);
-        run_signorini(msh, ap);
-    }
-    #endif
-
     return 0;
 }
diff --git a/apps/diffusion/src/diffusion_hho.cpp b/apps/diffusion/src/diffusion_hho.cpp
index 37ced564..dc531cc8 100644
--- a/apps/diffusion/src/diffusion_hho.cpp
+++ b/apps/diffusion/src/diffusion_hho.cpp
@@ -66,9 +66,9 @@ int main(int argc, char **argv)
                 tmpdeg = std::stoi(optarg);
                 if (not lflag)
                     hdi.cell_degree(tmpdeg);
-                
+
                 hdi.face_degree(tmpdeg);
-                
+
                 if (not rflag)
                     hdi.reconstruction_degree(tmpdeg+1);
                 break;
@@ -105,9 +105,7 @@ int main(int argc, char **argv)
     {
         std::cout << "Mesh format: 1D uniform" << std::endl;
 
-        typedef disk::generic_mesh<T, 1>  mesh_type;
-
-        mesh_type msh;
+        disk::generic_mesh<T, 1>  msh;
         disk::uniform_mesh_loader<T, 1> loader(0, 1, num_elems);
         loader.populate_mesh(msh);
 
@@ -121,7 +119,8 @@ int main(int argc, char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.typ1$") ))
     {
         std::cout << "Guessed mesh format: FVCA5 2D" << std::endl;
-        auto msh = disk::load_fvca5_2d_mesh<T>(mesh_filename);
+        disk::generic_mesh<T, 2> msh;
+        disk::load_mesh_fvca5_2d<T>(mesh_filename, msh);
         run_hho_diffusion_solver(msh, hdi, stat_cond, stab_diam_F);
         return 0;
     }
@@ -130,10 +129,9 @@ int main(int argc, char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.mesh2d$") ))
     {
         std::cout << "Guessed mesh format: Netgen 2D" << std::endl;
-        auto msh = disk::load_netgen_2d_mesh<T>(mesh_filename);
-
+        disk::simplicial_mesh<T, 2> msh;
+        disk::load_mesh_netgen<T>(mesh_filename, msh);
         std::cout << msh.faces_size() << std::endl;
-
         run_hho_diffusion_solver(msh, hdi, stat_cond, stab_diam_F);
         return 0;
     }
@@ -142,17 +140,18 @@ int main(int argc, char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.quad$") ))
     {
         std::cout << "Guessed mesh format: DiSk++ Cartesian 2D" << std::endl;
-        auto msh = disk::load_cartesian_2d_mesh<T>(mesh_filename);
+        disk::cartesian_mesh<T, 2> msh;
+        disk::load_mesh_diskpp_cartesian<T>(mesh_filename, msh);
         run_hho_diffusion_solver(msh, hdi, stat_cond, stab_diam_F);
         return 0;
     }
 
-
     /* Netgen 3D */
     if (std::regex_match(mesh_filename, std::regex(".*\\.mesh$") ))
     {
         std::cout << "Guessed mesh format: Netgen 3D" << std::endl;
-        auto msh = disk::load_netgen_3d_mesh<T>(mesh_filename);
+        disk::simplicial_mesh<T, 3> msh;
+        disk::load_mesh_netgen<T>(mesh_filename, msh);
         run_hho_diffusion_solver(msh, hdi, stat_cond, stab_diam_F);
         return 0;
     }
@@ -161,7 +160,8 @@ int main(int argc, char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.hex$") ))
     {
         std::cout << "Guessed mesh format: DiSk++ Cartesian 3D" << std::endl;
-        auto msh = disk::load_cartesian_3d_mesh<T>(mesh_filename);
+        disk::cartesian_mesh<T, 3> msh;
+        disk::load_mesh_diskpp_cartesian<T>(mesh_filename, msh);
         run_hho_diffusion_solver(msh, hdi, stat_cond, stab_diam_F);
         return 0;
     }
@@ -171,11 +171,8 @@ int main(int argc, char **argv)
     {
         std::cout << "Guessed mesh format: FVCA6 3D" << std::endl;
         disk::generic_mesh<T,3> msh;
-
         disk::load_mesh_fvca6_3d<T>(mesh_filename, msh);
-
         run_hho_diffusion_solver(msh, hdi, stat_cond, stab_diam_F);
-
         return 0;
     }
 }
diff --git a/apps/diskpp_tutorials/diskpp_tutorial_3.cpp b/apps/diskpp_tutorials/diskpp_tutorial_3.cpp
index d8f203b3..7e106002 100644
--- a/apps/diskpp_tutorials/diskpp_tutorial_3.cpp
+++ b/apps/diskpp_tutorials/diskpp_tutorial_3.cpp
@@ -75,11 +75,15 @@ int main(int argc, const char **argv)
     const char *mesh_filename = argv[1];
 
     using T = double;
+
     /* FVCA5 2D */
     if (std::regex_match(mesh_filename, std::regex(".*\\.typ1$") ))
     {
         std::cout << "Guessed mesh format: FVCA5 2D" << std::endl;
-        auto msh = disk::load_fvca5_2d_mesh<T>(mesh_filename);
+        disk::generic_mesh<T, 2> msh;
+        auto success = disk::load_mesh_fvca5_2d<T>(mesh_filename, msh);
+        if (!success)
+            return 1;
         run(msh);
         return 0;
     }
@@ -88,7 +92,10 @@ int main(int argc, const char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.mesh2d$") ))
     {
         std::cout << "Guessed mesh format: Netgen 2D" << std::endl;
-        auto msh = disk::load_netgen_2d_mesh<T>(mesh_filename);
+        disk::simplicial_mesh<T, 2> msh;
+        auto success = disk::load_mesh_netgen<T>(mesh_filename, msh);
+        if (!success)
+            return 1;
         run(msh);
         return 0;
     }
@@ -97,17 +104,22 @@ int main(int argc, const char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.quad$") ))
     {
         std::cout << "Guessed mesh format: DiSk++ Cartesian 2D" << std::endl;
-        auto msh = disk::load_cartesian_2d_mesh<T>(mesh_filename);
+        disk::cartesian_mesh<T, 2> msh;
+        auto success = disk::load_mesh_diskpp_cartesian<T>(mesh_filename, msh);
+        if (!success)
+            return 1;
         run(msh);
         return 0;
     }
 
-
     /* Netgen 3D */
     if (std::regex_match(mesh_filename, std::regex(".*\\.mesh$") ))
     {
         std::cout << "Guessed mesh format: Netgen 3D" << std::endl;
-        auto msh = disk::load_netgen_3d_mesh<T>(mesh_filename);
+        disk::simplicial_mesh<T, 3> msh;
+        auto success = disk::load_mesh_netgen<T>(mesh_filename, msh);
+        if (!success)
+            return 1;
         run(msh);
         return 0;
     }
@@ -116,7 +128,10 @@ int main(int argc, const char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.hex$") ))
     {
         std::cout << "Guessed mesh format: DiSk++ Cartesian 3D" << std::endl;
-        auto msh = disk::load_cartesian_3d_mesh<T>(mesh_filename);
+        disk::cartesian_mesh<T, 3> msh;
+        auto success = disk::load_mesh_diskpp_cartesian<T>(mesh_filename, msh);
+        if (!success)
+            return 1;
         run(msh);
         return 0;
     }
@@ -126,13 +141,14 @@ int main(int argc, const char **argv)
     {
         std::cout << "Guessed mesh format: FVCA6 3D" << std::endl;
         disk::generic_mesh<T,3> msh;
-        disk::load_mesh_fvca6_3d<T>(mesh_filename, msh);
+        auto success = disk::load_mesh_fvca6_3d<T>(mesh_filename, msh);
+        if (!success)
+            return 1;
         run(msh);
         return 0;
     }
 
 
     return 0;
-}
-
 
+}
diff --git a/apps/helmholtz/helmholtz.cpp b/apps/helmholtz/helmholtz.cpp
index d19cb4d6..dd881f11 100644
--- a/apps/helmholtz/helmholtz.cpp
+++ b/apps/helmholtz/helmholtz.cpp
@@ -198,8 +198,8 @@ run_hho_helmholtz_solver(Mesh& msh, size_t degree, const typename Mesh::coordina
     }
 
     std::cout << "h = " << disk::average_diameter(msh) << std::endl;
-    std::cout << "L2 error:     " << std::sqrt(l2_error) << std::endl; 
-    std::cout << "Energy error: " << std::sqrt(nrg_error) << std::endl; 
+    std::cout << "L2 error:     " << std::sqrt(l2_error) << std::endl;
+    std::cout << "Energy error: " << std::sqrt(nrg_error) << std::endl;
 
     disk::silo_database silo_db;
     silo_db.create("helmholtz.silo");
@@ -248,8 +248,8 @@ int main(int argc, char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.mesh2d$") ))
     {
         std::cout << "Guessed mesh format: Netgen 2D" << std::endl;
-        auto msh = disk::load_netgen_2d_mesh<T>(mesh_filename);
-
+        disk::simplicial_mesh<T, 2> msh;
+        disk::load_mesh_netgen<T>(mesh_filename, msh);
         run_hho_helmholtz_solver(msh, degree, stab_param);
 
         return 0;
@@ -259,7 +259,8 @@ int main(int argc, char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.typ1$") ))
     {
         std::cout << "Guessed mesh format: FVCA5 2D" << std::endl;
-        auto msh = disk::load_fvca5_2d_mesh<T>(mesh_filename);
+        disk::generic_mesh<T, 2> msh;
+        disk::load_mesh_fvca5_2d<T>(mesh_filename, msh);
 
         run_hho_helmholtz_solver(msh, degree, stab_param);
 
@@ -270,7 +271,8 @@ int main(int argc, char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.mesh$") ))
     {
         std::cout << "Guessed mesh format: Netgen 3D" << std::endl;
-        auto msh = disk::load_netgen_3d_mesh<T>(mesh_filename);
+        disk::simplicial_mesh<T, 3> msh;
+        disk::load_mesh_netgen<T>(mesh_filename, msh);
 
         run_hho_helmholtz_solver(msh, degree, stab_param);
 
@@ -281,10 +283,11 @@ int main(int argc, char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.hex$") ))
     {
         std::cout << "Guessed mesh format: DiSk++ Cartesian 3D" << std::endl;
-        auto msh = disk::load_cartesian_3d_mesh<T>(mesh_filename);
-        
+        disk::cartesian_mesh<T, 3> msh;
+        disk::load_mesh_diskpp_cartesian<T>(mesh_filename, msh);
+
         run_hho_helmholtz_solver(msh, degree, stab_param);
-        
+
         return 0;
     }
 
@@ -293,11 +296,10 @@ int main(int argc, char **argv)
     {
         std::cout << "Guessed mesh format: FVCA6 3D" << std::endl;
         disk::generic_mesh<T,3> msh;
-        
         disk::load_mesh_fvca6_3d<T>(mesh_filename, msh);
-        
+
         run_hho_helmholtz_solver(msh, degree, stab_param);
-        
+
         return 0;
     }
 }
\ No newline at end of file
diff --git a/apps/linear_elasticity/src/linear_elasticity.cpp b/apps/linear_elasticity/src/linear_elasticity.cpp
index 43bb12d5..37a26cf9 100644
--- a/apps/linear_elasticity/src/linear_elasticity.cpp
+++ b/apps/linear_elasticity/src/linear_elasticity.cpp
@@ -267,7 +267,8 @@ main(int argc, char** argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.typ1$")))
     {
         std::cout << "Guessed mesh format: FVCA5 2D" << std::endl;
-        auto msh = disk::load_fvca5_2d_mesh<RealType>(mesh_filename);
+        disk::generic_mesh<RealType,2> msh;
+        disk::load_mesh_fvca5_2d<RealType>(mesh_filename, msh);
         run_linear_elasticity_solver(msh, rp, material_data);
         return 0;
     }
@@ -276,7 +277,8 @@ main(int argc, char** argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.medit2d$")))
     {
         std::cout << "Guessed mesh format: Medit format" << std::endl;
-        auto msh = disk::load_medit_2d_mesh<RealType>(mesh_filename);
+        disk::generic_mesh<RealType, 2> msh;
+        disk::load_mesh_medit<RealType>(mesh_filename, msh);
         run_linear_elasticity_solver(msh, rp, material_data);
         return 0;
     }
@@ -285,7 +287,8 @@ main(int argc, char** argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.mesh2d$")))
     {
         std::cout << "Guessed mesh format: Netgen 2D" << std::endl;
-        auto msh = disk::load_netgen_2d_mesh<RealType>(mesh_filename);
+        disk::simplicial_mesh<RealType, 2> msh;
+        disk::load_mesh_netgen<RealType>(mesh_filename, msh);
         run_linear_elasticity_solver(msh, rp, material_data);
         return 0;
     }
@@ -294,7 +297,8 @@ main(int argc, char** argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.quad$")))
     {
         std::cout << "Guessed mesh format: DiSk++ Cartesian 2D" << std::endl;
-        auto msh = disk::load_cartesian_2d_mesh<RealType>(mesh_filename);
+        disk::cartesian_mesh<RealType, 2> msh;
+        disk::load_mesh_diskpp_cartesian<RealType>(mesh_filename, msh);
         run_linear_elasticity_solver(msh, rp, material_data);
         return 0;
     }
@@ -303,16 +307,18 @@ main(int argc, char** argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.mesh$")))
     {
         std::cout << "Guessed mesh format: Netgen 3D" << std::endl;
-        auto msh = disk::load_netgen_3d_mesh<RealType>(mesh_filename);
+        disk::simplicial_mesh<RealType, 3> msh;
+        disk::load_mesh_netgen<RealType>(mesh_filename, msh);
         run_linear_elasticity_solver(msh, rp, material_data);
         return 0;
     }
 
-    /* FVCA5 2D */
+    /* FVCA6 3D */
     if (std::regex_match(mesh_filename, std::regex(".*\\.msh$")))
     {
         std::cout << "Guessed mesh format: FVCA6 3D" << std::endl;
-        auto msh = disk::load_fvca6_3d_mesh<RealType>(mesh_filename);
+        disk::generic_mesh<RealType,3> msh;
+        disk::load_mesh_fvca6_3d<RealType>(mesh_filename, msh);
         run_linear_elasticity_solver(msh, rp, material_data);
         return 0;
     }
@@ -321,7 +327,8 @@ main(int argc, char** argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.hex$")))
     {
         std::cout << "Guessed mesh format: DiSk++ Cartesian 3D" << std::endl;
-        auto msh = disk::load_cartesian_3d_mesh<RealType>(mesh_filename);
+        disk::cartesian_mesh<RealType, 3> msh;
+        disk::load_mesh_diskpp_cartesian<RealType>(mesh_filename, msh);
         run_linear_elasticity_solver(msh, rp, material_data);
         return 0;
     }
diff --git a/apps/linear_elasticity/src/linear_elasticity_test.cpp b/apps/linear_elasticity/src/linear_elasticity_test.cpp
index 2c0da490..4ed48e4f 100644
--- a/apps/linear_elasticity/src/linear_elasticity_test.cpp
+++ b/apps/linear_elasticity/src/linear_elasticity_test.cpp
@@ -318,7 +318,8 @@ test_triangles_fvca5(const run_params& rp, const ElasticityParameters material_d
 
     for (int i = 0; i < runs; i++)
     {
-        auto msh = disk::load_fvca5_2d_mesh<T>(paths[i].c_str());
+        disk::generic_mesh<T, 2> msh;
+        disk::load_mesh_fvca5_2d<T>(paths[i].c_str(), msh);
         error_sumup.push_back(run_linear_elasticity_solver(msh, rp, material_data));
     }
     printResults(error_sumup);
@@ -341,7 +342,8 @@ test_triangles_netgen(const run_params& rp, const ElasticityParameters material_
 
     for (int i = 0; i < runs; i++)
     {
-        auto msh = disk::load_netgen_2d_mesh<T>(paths[i].c_str());
+        disk::simplicial_mesh<T, 2> msh;
+        disk::load_mesh_netgen<T>(paths[i].c_str(), msh);
         error_sumup.push_back(run_linear_elasticity_solver(msh, rp, material_data));
     }
     printResults(error_sumup);
@@ -364,7 +366,8 @@ test_hexagons(const run_params& rp, const ElasticityParameters material_data)
 
     for (int i = 0; i < runs; i++)
     {
-        auto msh = disk::load_fvca5_2d_mesh<T>(paths[i].c_str());
+        disk::generic_mesh<T, 2> msh;
+        disk::load_mesh_fvca5_2d<T>(paths[i].c_str(), msh);
         error_sumup.push_back(run_linear_elasticity_solver(msh, rp, material_data));
     }
     printResults(error_sumup);
@@ -387,7 +390,8 @@ test_kershaws(const run_params& rp, const ElasticityParameters material_data)
 
     for (int i = 0; i < runs; i++)
     {
-        auto msh = disk::load_fvca5_2d_mesh<T>(paths[i].c_str());
+        disk::generic_mesh<T, 2> msh;
+        disk::load_mesh_fvca5_2d<T>(paths[i].c_str(), msh);
         error_sumup.push_back(run_linear_elasticity_solver(msh, rp, material_data));
     }
     printResults(error_sumup);
@@ -410,7 +414,8 @@ test_quads_fvca5(const run_params& rp, const ElasticityParameters material_data)
 
     for (int i = 0; i < runs; i++)
     {
-        auto msh = disk::load_fvca5_2d_mesh<T>(paths[i].c_str());
+        disk::generic_mesh<T, 2> msh;
+        disk::load_mesh_fvca5_2d<T>(paths[i].c_str(), msh);
         error_sumup.push_back(run_linear_elasticity_solver(msh, rp, material_data));
     }
     printResults(error_sumup);
@@ -433,7 +438,10 @@ test_quads_diskpp(const run_params& rp, const ElasticityParameters material_data
 
     for (int i = 0; i < runs; i++)
     {
-        auto msh = disk::load_cartesian_2d_mesh<T>(paths[i].c_str());
+        disk::cartesian_mesh<T, 2> msh;
+        disk::load_mesh_diskpp_cartesian<T>(paths[i].c_str(), msh);
+
+
         error_sumup.push_back(run_linear_elasticity_solver(msh, rp, material_data));
     }
     printResults(error_sumup);
@@ -456,7 +464,8 @@ test_hexahedra_diskpp(const run_params& rp, const ElasticityParameters material_
 
     for (int i = 0; i < runs; i++)
     {
-        auto msh = disk::load_cartesian_3d_mesh<T>(paths[i].c_str());
+        disk::cartesian_mesh<T, 3> msh;
+        disk::load_mesh_diskpp_cartesian<T>(paths[i].c_str(), msh);
         error_sumup.push_back(run_linear_elasticity_solver(msh, rp, material_data));
     }
     printResults(error_sumup);
@@ -479,7 +488,9 @@ test_hexahedra_fvca6(const run_params& rp, const ElasticityParameters material_d
 
     for (int i = 0; i < runs; i++)
     {
-        auto msh = disk::load_fvca6_3d_mesh<T>(paths[i].c_str());
+        disk::generic_mesh<T,3> msh;
+        disk::load_mesh_fvca6_3d<T>(paths[i].c_str(), msh);
+
         error_sumup.push_back(run_linear_elasticity_solver(msh, rp, material_data));
     }
     printResults(error_sumup);
@@ -502,7 +513,8 @@ test_tetrahedra_netgen(const run_params& rp, const ElasticityParameters material
 
     for (int i = 0; i < runs; i++)
     {
-        auto msh = disk::load_netgen_3d_mesh<T>(paths[i].c_str());
+        disk::simplicial_mesh<T, 3> msh;
+        disk::load_mesh_netgen<T>(paths[i].c_str(), msh);
         error_sumup.push_back(run_linear_elasticity_solver(msh, rp, material_data));
     }
     printResults(error_sumup);
@@ -524,7 +536,9 @@ test_polyhedra_fvca6(const run_params& rp, const ElasticityParameters material_d
 
     for (int i = 0; i < runs; i++)
     {
-        auto msh = disk::load_fvca6_3d_mesh<T>(paths[i].c_str());
+        disk::generic_mesh<T,3> msh;
+        disk::load_mesh_fvca6_3d<T>(paths[i].c_str(), msh);
+
         error_sumup.push_back(run_linear_elasticity_solver(msh, rp, material_data));
     }
     printResults(error_sumup);
@@ -547,7 +561,8 @@ test_tetrahedra_fvca6(const run_params& rp, const ElasticityParameters material_
 
     for (int i = 0; i < runs; i++)
     {
-        auto msh = disk::load_fvca6_3d_mesh<T>(paths[i].c_str());
+        disk::generic_mesh<T,3> msh;
+        disk::load_mesh_fvca6_3d<T>(paths[i].c_str(), msh);
         error_sumup.push_back(run_linear_elasticity_solver(msh, rp, material_data));
     }
     printResults(error_sumup);
diff --git a/apps/maxwell/src/maxwell.cpp b/apps/maxwell/src/maxwell.cpp
index 57f32725..6e358264 100644
--- a/apps/maxwell/src/maxwell.cpp
+++ b/apps/maxwell/src/maxwell.cpp
@@ -1,11 +1,11 @@
 /*
  * DISK++, a template library for DIscontinuous SKeletal methods.
- *  
+ *
  * Matteo Cicuttin (C) 2020,2021
  * matteo.cicuttin@uliege.be
  *
  * University of Liège - Montefiore Institute
- * Applied and Computational Electromagnetics group  
+ * Applied and Computational Electromagnetics group
  */
 
 #include <iostream>
@@ -45,7 +45,7 @@ class maxwell_assembler
     std::vector<size_t>     expand_table;
 
 public:
-    
+
     SparseMatrix<T>         LHS;
     Matrix<T, Dynamic, 1>   RHS;
 
@@ -212,7 +212,7 @@ class maxwell_assembler
         LHS.setFromTriplets(triplets.begin(), triplets.end());
         triplets.clear();
         std::cout << "Solving for " << syssz << " DOFs. ";
-        std::cout << "Matrix has " << LHS.nonZeros() << " nonzeros." << std::endl; 
+        std::cout << "Matrix has " << LHS.nonZeros() << " nonzeros." << std::endl;
     }
 };
 
@@ -264,7 +264,7 @@ class maxwell_assembler_condensed
     }
 
 public:
-    
+
     SparseMatrix<ScalT>         LHS;
     Matrix<ScalT, Dynamic, 1>   RHS;
 
@@ -333,7 +333,7 @@ class maxwell_assembler_condensed
                 auto lofsi = fi*fbs;
                 auto lofsj = fj*fbs;
 
-                if ( not is_in_system(msh, fcs[fj]) ) 
+                if ( not is_in_system(msh, fcs[fj]) )
                 {
                     RHS.segment(cofsi, fbs) += -lhsc.block(lofsi, lofsj, fbs, fbs)*dirichlet_data.segment(lofsj, fbs);
                     continue;
@@ -354,7 +354,7 @@ class maxwell_assembler_condensed
         LHS.setFromTriplets(triplets.begin(), triplets.end());
         triplets.clear();
         std::cout << "Solving for " << syssz << " DOFs. ";
-        std::cout << "Matrix has " << LHS.nonZeros() << " nonzeros." << std::endl; 
+        std::cout << "Matrix has " << LHS.nonZeros() << " nonzeros." << std::endl;
     }
 
     disk::dynamic_vector<ScalT>
@@ -517,7 +517,7 @@ class maxwell_eigenvalue_assembler
             {
                 if ( msh.is_boundary(fcs[fj]) )
                     continue;
-            
+
                 auto dofsi = cbs*msh.cells_size() + fbs * compress_table.at( offset(msh, fcs[fi]) );
                 auto dofsj = cbs*msh.cells_size() + fbs * compress_table.at( offset(msh, fcs[fj]) );
 
@@ -526,7 +526,7 @@ class maxwell_eigenvalue_assembler
                         trpK.push_back( Triplet<T>(dofsi+i, dofsj+j, K(cbs+fi*fbs+i, cbs+fj*fbs+j)) );
 
             }
-        }  
+        }
     }
 
     void finalize(void)
@@ -538,7 +538,7 @@ class maxwell_eigenvalue_assembler
 
         std::cout << "Solving for " << syssz << " DOFs. ";
         std::cout << "K has " << gK.nonZeros() << " nonzeros, M has ";
-        std::cout << gM.nonZeros() << " nonzeros." << std::endl; 
+        std::cout << gM.nonZeros() << " nonzeros." << std::endl;
     }
 };
 
@@ -556,7 +556,7 @@ vector_wave_solver(Mesh<T,2,Storage>& msh, size_t order,
     typedef Mesh<T,2,Storage>                   mesh_type;
     typedef typename mesh_type::coordinate_type scalar_type;
     typedef typename mesh_type::point_type      point_type;
-    
+
     disk::hho_degree_info chdi( { .rd = (size_t) order,
                                   .cd = (size_t) order,
                                   .fd = (size_t) order } );
@@ -632,7 +632,7 @@ vector_wave_solver(Mesh<T,2,Storage>& msh, size_t order,
 #ifdef USE_STATIC_CONDENSATION
         auto CR = disk::curl_reconstruction(msh, cl, chdi);
         auto ST = disk::curl_hdg_stabilization(msh, cl, chdi);
-        
+
         Matrix<T, Dynamic, Dynamic> MM = Matrix<T, Dynamic, Dynamic>::Zero(ST.rows(), ST.cols());
         MM.block(0,0,cb.size(),cb.size()) = make_mass_matrix(msh, cl, cb);
 
@@ -648,7 +648,7 @@ vector_wave_solver(Mesh<T,2,Storage>& msh, size_t order,
         Matrix<T, Dynamic, 1> prj = project_tangent(msh, cl, chdi, sol_fun);
 
         Matrix<T, Dynamic, 1> hsol = Matrix<T, Dynamic, 1>::Zero(rb.size());
-        
+
         hsol.segment(1, rb.size()-1) = CR.first * esol;
 
         energy_err += (prj-esol).dot(lhs*(prj-esol));
@@ -664,7 +664,7 @@ vector_wave_solver(Mesh<T,2,Storage>& msh, size_t order,
         for (auto& qp : qps)
         {
             Matrix<T, Dynamic, 2> hphi  = rb.eval_curls2(qp.point());
-            Matrix<T, Dynamic, 2> hphi2 = hphi.block(0,0,cb.size(),2); 
+            Matrix<T, Dynamic, 2> hphi2 = hphi.block(0,0,cb.size(),2);
             Matrix<T, Dynamic, 1> ephi  = cb.eval_functions(qp.point());
             //Matrix<T, Dynamic, 3> rphi = rb.eval_functions(qp.point());
 
@@ -715,7 +715,7 @@ vector_wave_solver(Mesh<T,3,Storage>& msh, size_t order)
     typedef Mesh<T,3,Storage>                   mesh_type;
     typedef typename mesh_type::coordinate_type scalar_type;
     typedef typename mesh_type::point_type      point_type;
-    
+
     disk::hho_degree_info chdi( { .rd = (size_t) order,
                                   .cd = (size_t) order+1,
                                   .fd = (size_t) order } );
@@ -797,7 +797,7 @@ vector_wave_solver(Mesh<T,3,Storage>& msh, size_t order)
         auto CR = disk::curl_reconstruction_pk(msh, cl, chdi);
 
         Matrix<T, Dynamic, Dynamic> ST = disk::curl_hdg_stabilization(msh, cl, chdi);
-        
+
         auto MM = make_vector_mass_oper(msh, cl, chdi);
         Matrix<T, Dynamic, Dynamic> lhs = CR.second + omega*(ST - omega*MM);
         Matrix<T, Dynamic, 1> rhs = Matrix<T, Dynamic, 1>::Zero(lhs.rows());
@@ -899,7 +899,7 @@ vector_wave_solver_complex(Mesh<CoordT,3,Storage>& msh, parameter_loader<Mesh<Co
     auto jw = scalar_type(0,omega);
     auto jwmu0 = scalar_type(0,omega*mu0);
     auto ksq = (eps0*omega)*(mu0*omega);
-    
+
     auto& lua = pl.state();
 
     double wgz = lua["wgz"];
@@ -994,7 +994,7 @@ vector_wave_solver_complex(Mesh<CoordT,3,Storage>& msh, parameter_loader<Mesh<Co
         {
             disk::dynamic_vector<scalar_type> tfsf =
                         disk::dynamic_vector<scalar_type>::Zero(lhs.rows());
-            
+
             auto qps = disk::integrate(msh, cl, 2*chdi.cell_degree());
             for (const auto& qp : qps)
             {
@@ -1021,7 +1021,7 @@ vector_wave_solver_complex(Mesh<CoordT,3,Storage>& msh, parameter_loader<Mesh<Co
                 if (di.tag() == 14 and bi.tag() == 1188)
                 //if (di.tag() == 5)
                 {
-                    const auto fb = disk::make_vector_monomial_tangential_basis<mesh_type, face_type, scalar_type>(msh, fc, chdi.face_degree()); 
+                    const auto fb = disk::make_vector_monomial_tangential_basis<mesh_type, face_type, scalar_type>(msh, fc, chdi.face_degree());
 
                     disk::dynamic_matrix<scalar_type> tfsf_mass =
                         disk::dynamic_matrix<scalar_type>::Zero(fbs, fbs);
@@ -1042,11 +1042,11 @@ vector_wave_solver_complex(Mesh<CoordT,3,Storage>& msh, parameter_loader<Mesh<Co
                     }
                     //if (bi.tag() == 7)
                     {
-                        disk::dynamic_vector<scalar_type> s = disk::dynamic_vector<scalar_type>::Zero(lhs.rows()); 
+                        disk::dynamic_vector<scalar_type> s = disk::dynamic_vector<scalar_type>::Zero(lhs.rows());
                         s.segment(cbs+i*fbs, fbs) = tfsf_mass.ldlt().solve(tfsf_rhs);
                         rhs += lhst*s;
                         //auto Y = disk::make_impedance_term(msh, fc, chdi.face_degree());
-                        rhs.segment(cbs+i*fbs, fbs) -= (jwmu0/wgz)*tfsf_rhs; //Y*tfsf_mass.ldlt().solve(tfsf_rhs); 
+                        rhs.segment(cbs+i*fbs, fbs) -= (jwmu0/wgz)*tfsf_rhs; //Y*tfsf_mass.ldlt().solve(tfsf_rhs);
                         //rhs.segment(cbs+i*fbs, fbs) += (omega*mu0/wgz)*tfsf_rhs;
                     }
 
@@ -1272,7 +1272,7 @@ vector_wave_solver_complex(Mesh<CoordT,3,Storage>& msh, parameter_loader<Mesh<Co
 
     std::vector<scalar_type> nodal_ex( data_ex.size() );
     std::transform(data_ex.begin(), data_ex.end(), nodal_ex.begin(), tran);
-    
+
     std::vector<scalar_type> nodal_ey( data_ey.size() );
     std::transform(data_ey.begin(), data_ey.end(), nodal_ey.begin(), tran);
 
@@ -1420,7 +1420,7 @@ void maxwell_eigenvalue_solver(Mesh& msh, size_t order, const typename Mesh::coo
 
         disk::silo_zonal_variable<T> uz("uz", data_uz);
         silo_db.add_variable("mesh", uz);
-    
+
         silo_db.add_expression("u", "{ux, uy, uz}", DB_VARTYPE_VECTOR);
         silo_db.add_expression("mag_u", "magnitude(u)", DB_VARTYPE_SCALAR);
     }
@@ -1430,7 +1430,7 @@ void maxwell_eigenvalue_solver(Mesh& msh, size_t order, const typename Mesh::coo
         std::cout << std::setprecision(8) << hho_eigvals(i) << " -> ";
         std::cout << std::sqrt( hho_eigvals(i) ) << std::endl;
     }
-    
+
 }
 #endif
 
@@ -1438,7 +1438,7 @@ void maxwell_eigenvalue_solver(Mesh& msh, size_t order, const typename Mesh::coo
 void autotest_alpha(size_t order)
 {
     using T = double;
-    
+
     using Mesh = disk::simplicial_mesh<T,3>;
 
     std::stringstream ssl2;
@@ -1454,7 +1454,7 @@ void autotest_alpha(size_t order)
     {
         ofs_l2  << sp << " ";
         ofs_nrg << sp << " ";
-        
+
         for (size_t i = 1; i < 5; i++)
         {
             Mesh msh;
@@ -1465,7 +1465,7 @@ void autotest_alpha(size_t order)
 
             auto diam = disk::average_diameter(msh);
 
-            auto ci = vector_wave_solver(msh, order, sp, M_PI, false);  
+            auto ci = vector_wave_solver(msh, order, sp, M_PI, false);
             ofs_l2  << ci.l2_error_e << " " << ci.l2_error_h << " ";
             ofs_nrg << ci.nrg_error << " ";
         }
@@ -1480,7 +1480,7 @@ void autotest_alpha(size_t order)
 void autotest_convergence(size_t order_min, size_t order_max)
 {
     using T = double;
-    
+
     using Mesh = disk::simplicial_mesh<T,3>;
 
     std::ofstream ofs( "hho_convergence_convt.txt" );
@@ -1626,7 +1626,8 @@ int main(int argc, char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.mesh$") ))
     {
         std::cout << "Guessed mesh format: Netgen 3D" << std::endl;
-        auto msh = disk::load_netgen_3d_mesh<coord_T>(mesh_filename);
+        disk::simplicial_mesh<coord_T, 3> msh;
+        disk::load_mesh_netgen<coord_T>(mesh_filename, msh);
 
         //if (solve_eigvals)
         //    maxwell_eigenvalue_solver(msh, degree, stab_param);
@@ -1645,7 +1646,7 @@ int main(int argc, char **argv)
         //disk::gmsh_geometry_loader< disk::simplicial_mesh<T,3> > loader;
         disk::generic_mesh<T,3> msh;
         disk::gmsh_geometry_loader< disk::generic_mesh<T,3> > loader;
-        
+
         loader.read_mesh(mesh_filename);
         loader.populate_mesh(msh);
 
@@ -1662,7 +1663,7 @@ int main(int argc, char **argv)
         disk::gmsh_geometry_loader< disk::simplicial_mesh<T,3> > loader;
         //disk::generic_mesh<T,3> msh;
         //disk::gmsh_geometry_loader< disk::generic_mesh<T,3> > loader;
-        
+
         loader.read_mesh(mesh_filename);
         loader.populate_mesh(msh);
 
@@ -1676,14 +1677,15 @@ int main(int argc, char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.hex$") ))
     {
         std::cout << "Guessed mesh format: DiSk++ Cartesian 3D" << std::endl;
-        auto msh = disk::load_cartesian_3d_mesh<coord_T>(mesh_filename);
-        
+        disk::cartesian_mesh<coord_T, 3> msh;
+        disk::load_mesh_diskpp_cartesian<coord_T>(mesh_filename, msh);
+
         if (solve_eigvals)
             maxwell_eigenvalue_solver(msh, degree, stab_param);
         else
             //vector_wave_solver(msh, degree, stab_param, omega);
             vector_wave_solver_complex(msh, degree);
-        
+
         return 0;
     }
 
@@ -1692,19 +1694,18 @@ int main(int argc, char **argv)
     {
         std::cout << "Guessed mesh format: FVCA6 3D" << std::endl;
         disk::generic_mesh<coord_T,3> msh;
-        
+
         disk::load_mesh_fvca6_3d<coord_T>(mesh_filename, msh);
-        
+
         if (solve_eigvals)
             maxwell_eigenvalue_solver(msh, degree, stab_param);
         else
             //vector_wave_solver(msh, degree, stab_param, omega);
             vector_wave_solver_complex(msh, degree);
-        
+
         return 0;
     }
 #endif
 
     return 0;
 }
-
diff --git a/apps/maxwell/src/maxwell_sip_dg.cpp b/apps/maxwell/src/maxwell_sip_dg.cpp
index 2d26bc8c..abd4866b 100644
--- a/apps/maxwell/src/maxwell_sip_dg.cpp
+++ b/apps/maxwell/src/maxwell_sip_dg.cpp
@@ -1,11 +1,11 @@
 /*
  * DISK++, a template library for DIscontinuous SKeletal methods.
- *  
+ *
  * Matteo Cicuttin (C) 2020
  * matteo.cicuttin@uliege.be
  *
  * University of Liège - Montefiore Institute
- * Applied and Computational Electromagnetics group  
+ * Applied and Computational Electromagnetics group
  */
 
 #include <iostream>
@@ -215,7 +215,7 @@ run_maxwell_solver(Mesh& msh, size_t degree, const typename Mesh::coordinate_typ
 
         auto tbasis = disk::make_vector_monomial_basis(msh, tcl, degree);
         auto qps = disk::integrate(msh, tcl, 2*degree);
-        
+
         matrix_type M = matrix_type::Zero(tbasis.size(), tbasis.size());
         matrix_type K = matrix_type::Zero(tbasis.size(), tbasis.size());
         vector_type loc_rhs = vector_type::Zero(tbasis.size());
@@ -223,7 +223,7 @@ run_maxwell_solver(Mesh& msh, size_t degree, const typename Mesh::coordinate_typ
         {
             auto phi = tbasis.eval_functions(qp.point());
             auto cphi = tbasis.eval_curls2(qp.point());
-            
+
             M += qp.weight() * phi * phi.transpose();
             K += qp.weight() * cphi * cphi.transpose();
             //loc_rhs += qp.weight() * phi * f(qp.point());
@@ -240,12 +240,12 @@ run_maxwell_solver(Mesh& msh, size_t degree, const typename Mesh::coordinate_typ
 
             matrix_type Att = matrix_type::Zero(tbasis.size(), tbasis.size());
             matrix_type Atn = matrix_type::Zero(tbasis.size(), tbasis.size());
-            
+
             auto nv = cvf.neighbour_via(msh, tcl, fc);
             auto ncl = nv.first;
             auto nbasis = disk::make_vector_monomial_basis(msh, ncl, degree);
             assert(tbasis.size() == nbasis.size());
-            
+
             auto n     = normal(msh, tcl, fc);
             auto eta_l = eta / diameter(msh, fc);
             auto f_qps = disk::integrate(msh, fc, 2*degree);
@@ -256,7 +256,7 @@ run_maxwell_solver(Mesh& msh, size_t degree, const typename Mesh::coordinate_typ
                 auto tcphi      = tbasis.eval_curls2(fqp.point());
                 auto n_x_tphi   = disk::vcross(n, tphi);
                 auto n_x_tphi_x_n = disk::vcross(n_x_tphi, n);
-                
+
                 if (nv.second)
                 {   /* NOT on a boundary */
                     Att += + fqp.weight() * eta_l * n_x_tphi * n_x_tphi.transpose();
@@ -273,7 +273,7 @@ run_maxwell_solver(Mesh& msh, size_t degree, const typename Mesh::coordinate_typ
                         Att += - fqp.weight() * n_x_tphi * tcphi.transpose();
                     }
 
-                    
+
                     //loc_rhs -= fqp.weight() * tcphi_x_n;
                     //loc_rhs += fqp.weight() * eta_l * tphi;
 
@@ -285,16 +285,16 @@ run_maxwell_solver(Mesh& msh, size_t degree, const typename Mesh::coordinate_typ
 
                     continue;
                 }
-                
+
                 auto nphi       = nbasis.eval_functions(fqp.point());
                 auto ncphi      = nbasis.eval_curls2(fqp.point());
                 auto n_x_nphi   = disk::vcross(n, nphi);
-                
+
                 Atn += - fqp.weight() * eta_l * n_x_tphi * n_x_nphi.transpose();
                 Atn += - fqp.weight() * 0.5 * n_x_tphi * ncphi.transpose();
                 Atn += + fqp.weight() * 0.5 * tcphi * n_x_nphi.transpose();
             }
-            
+
             assm.assemble(msh, tcl, tcl, Att);
             if (nv.second)
                 assm.assemble(msh, tcl, ncl, Atn);
@@ -418,14 +418,14 @@ run_maxwell_solver(Mesh& msh, size_t degree, const typename Mesh::coordinate_typ
 
                 auto cb = make_vector_monomial_basis(msh, cl, degree);
                 Matrix<scalar_type, Dynamic, 1> lsol = sol.segment(cell_i*cb.size(), cb.size());
-            
+
                 auto phi = cb.eval_functions(pt);
                 auto ls = phi.transpose()*lsol;
 
                 voltage += ls(0)*eval_increment.x() +
                            ls(1)*eval_increment.y() +
                            ls(2)*eval_increment.z();
-            
+
                 continue;
             }
             cell_i++;
@@ -445,7 +445,7 @@ run_maxwell_solver(Mesh& msh, size_t degree, const typename Mesh::coordinate_typ
 
     std::vector<scalar_type> nodal_ex( data_ex.size() );
     std::transform(data_ex.begin(), data_ex.end(), nodal_ex.begin(), tran);
-    
+
     std::vector<scalar_type> nodal_ey( data_ey.size() );
     std::transform(data_ey.begin(), data_ey.end(), nodal_ey.begin(), tran);
 
@@ -491,12 +491,12 @@ run_maxwell_solver(Mesh& msh, size_t degree, const typename Mesh::coordinate_typ
 template<typename Mesh>
 void
 run_maxwell_eigenvalue_solver(Mesh& msh, size_t degree, const typename Mesh::coordinate_type eta)
-{   
+{
     auto cvf = connectivity_via_faces(msh);
     using T = typename Mesh::coordinate_type;
     typedef Matrix<T, Dynamic, Dynamic> matrix_type;
     typedef Matrix<T, Dynamic, 1>       vector_type;
-    
+
     auto f = make_rhs_function(msh);
 
     auto cbs = disk::vector_basis_size(degree, Mesh::dimension, Mesh::dimension);
@@ -506,7 +506,7 @@ run_maxwell_eigenvalue_solver(Mesh& msh, size_t degree, const typename Mesh::coo
     {
         auto tbasis = disk::make_vector_monomial_basis(msh, tcl, degree);
         auto qps = disk::integrate(msh, tcl, 2*degree);
-        
+
         matrix_type M = matrix_type::Zero(tbasis.size(), tbasis.size());
         matrix_type K = matrix_type::Zero(tbasis.size(), tbasis.size());
         vector_type loc_rhs = vector_type::Zero(tbasis.size());
@@ -514,7 +514,7 @@ run_maxwell_eigenvalue_solver(Mesh& msh, size_t degree, const typename Mesh::coo
         {
             auto phi = tbasis.eval_functions(qp.point());
             auto cphi = tbasis.eval_curls2(qp.point());
-            
+
             M += qp.weight() * phi * phi.transpose();
             K += qp.weight() * cphi * cphi.transpose();
         }
@@ -527,22 +527,22 @@ run_maxwell_eigenvalue_solver(Mesh& msh, size_t degree, const typename Mesh::coo
         {
             matrix_type Att = matrix_type::Zero(tbasis.size(), tbasis.size());
             matrix_type Atn = matrix_type::Zero(tbasis.size(), tbasis.size());
-            
+
             auto nv = cvf.neighbour_via(msh, tcl, fc);
             auto ncl = nv.first;
             auto nbasis = disk::make_vector_monomial_basis(msh, ncl, degree);
             assert(tbasis.size() == nbasis.size());
-            
+
             auto n     = normal(msh, tcl, fc);
             auto eta_l = eta / diameter(msh, fc);
             auto f_qps = disk::integrate(msh, fc, 2*degree);
-            
+
             for (auto& fqp : f_qps)
             {
                 auto tphi       = tbasis.eval_functions(fqp.point());
                 auto tcphi      = tbasis.eval_curls2(fqp.point());
                 auto n_x_tphi   = disk::vcross(n, tphi);
-                
+
                 if (nv.second)
                 {   /* NOT on a boundary */
                     Att += + fqp.weight() * eta_l * n_x_tphi * n_x_tphi.transpose();
@@ -556,21 +556,21 @@ run_maxwell_eigenvalue_solver(Mesh& msh, size_t degree, const typename Mesh::coo
                     Att += - fqp.weight() * tcphi * n_x_tphi.transpose();
                     continue;
                 }
-                
+
                 auto nphi       = nbasis.eval_functions(fqp.point());
                 auto ncphi      = nbasis.eval_curls2(fqp.point());
                 auto n_x_nphi   = disk::vcross(n, nphi);
-                
+
                 Atn += - fqp.weight() * eta_l * n_x_tphi * n_x_nphi.transpose();
                 Atn += - fqp.weight() * 0.5 * n_x_tphi * ncphi.transpose();
                 Atn += + fqp.weight() * 0.5 * tcphi * n_x_nphi.transpose();
             }
-            
+
             assm.assemble(msh, tcl, tcl, Att);
             if (nv.second)
                 assm.assemble(msh, tcl, ncl, Atn);
         }
-        
+
     }
 
     assm.finalize();
@@ -632,7 +632,7 @@ run_maxwell_eigenvalue_solver(Mesh& msh, size_t degree, const typename Mesh::coo
 
         disk::silo_zonal_variable<T> uz("uz", data_uz);
         silo_db.add_variable("mesh", uz);
-    
+
         silo_db.add_expression("u", "{ux, uy, uz}", DB_VARTYPE_VECTOR);
         silo_db.add_expression("mag_u", "magnitude(u)", DB_VARTYPE_SCALAR);
     }
@@ -650,7 +650,7 @@ run_maxwell_eigenvalue_solver(Mesh& msh, size_t degree, const typename Mesh::coo
 void autotest_convergence(size_t order_min, size_t order_max)
 {
     using T = double;
-    
+
     using Mesh = disk::simplicial_mesh<T,3>;
 
     std::ofstream ofs( "dg_convergence_convt.txt" );
@@ -728,7 +728,9 @@ int main(int argc, char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.mesh$") ))
     {
         std::cout << "Guessed mesh format: Netgen 3D" << std::endl;
-        auto msh = disk::load_netgen_3d_mesh<T>(mesh_filename);
+        disk::simplicial_mesh<T, 3> msh;
+        disk::load_mesh_netgen<T>(mesh_filename, msh);
+
         run_maxwell_solver(msh, degree, stab_param, param_filename);
         return 0;
     }
@@ -740,7 +742,7 @@ int main(int argc, char **argv)
         using Mesh = disk::simplicial_mesh<T,3>;
         Mesh msh;
         disk::gmsh_geometry_loader< Mesh > loader;
-        
+
         loader.read_mesh(mesh_filename);
         loader.populate_mesh(msh);
 
@@ -756,4 +758,3 @@ int main(int argc, char **argv)
         return 0;
     }
 }
-
diff --git a/apps/maxwell/src/maxwell_td.cpp b/apps/maxwell/src/maxwell_td.cpp
index 579b349d..0ce81169 100644
--- a/apps/maxwell/src/maxwell_td.cpp
+++ b/apps/maxwell/src/maxwell_td.cpp
@@ -1,11 +1,11 @@
 /*
  * DISK++, a template library for DIscontinuous SKeletal methods.
- *  
+ *
  * Matteo Cicuttin (C) 2020,2021
  * matteo.cicuttin@uliege.be
  *
  * University of Liège - Montefiore Institute
- * Applied and Computational Electromagnetics group  
+ * Applied and Computational Electromagnetics group
  */
 
 #include <iostream>
@@ -84,7 +84,7 @@ class maxwell_assembler_condensed
     }
 
 public:
-    
+
     SparseMatrix<ScalT>         LHS;
     SparseMatrix<ScalT>         Cond1;
     SparseMatrix<ScalT>         Cond2;
@@ -115,7 +115,7 @@ class maxwell_assembler_condensed
 
         sysfcs = std::count_if(msh.faces_begin(), msh.faces_end(), in_system);
         syssz = fbs*sysfcs;
-        fullsz = cbs*msh.cells_size() + syssz; 
+        fullsz = cbs*msh.cells_size() + syssz;
 
         LHS = SparseMatrix<ScalT>(syssz, syssz);
         Cond1 = SparseMatrix<ScalT>(syssz, fullsz);
@@ -187,7 +187,7 @@ class maxwell_assembler_condensed
                 auto lofsi = fi*fbs;
                 auto lofsj = fj*fbs;
 
-                if ( not is_in_system(msh, fcs[fj]) ) 
+                if ( not is_in_system(msh, fcs[fj]) )
                 {
                     //RHS.segment(cofsi, fbs) += -lhsc.block(lofsi, lofsj, fbs, fbs)*dirichlet_data.segment(lofsj, fbs);
                     continue;
@@ -234,11 +234,11 @@ class maxwell_assembler_condensed
 
         for (size_t fj = 0; fj < fcs.size(); fj++)
         {
-            if ( not is_in_system(msh, fcs[fj]) ) 
+            if ( not is_in_system(msh, fcs[fj]) )
                 continue;
-            
+
             auto cofsj = get_system_offset(msh, fcs[fj]);
-            
+
             for (size_t i = 0; i < cbs; i++)
             {
                 for (size_t j = 0; j < fbs; j++)
@@ -267,7 +267,7 @@ class maxwell_assembler_condensed
         triplets_decond2.clear();
         triplets_decond3.clear();
         std::cout << "Solving for " << syssz << " DOFs. ";
-        std::cout << "Matrix has " << LHS.nonZeros() << " nonzeros." << std::endl; 
+        std::cout << "Matrix has " << LHS.nonZeros() << " nonzeros." << std::endl;
     }
 
     disk::dynamic_vector<ScalT>
@@ -329,7 +329,7 @@ vector_wave_solver(Mesh<T,3,Storage>& msh, size_t order, T alpha, size_t num_ts)
     typedef typename mesh_type::cell_type       cell_type;
     typedef typename mesh_type::coordinate_type scalar_type;
     typedef typename mesh_type::point_type      point_type;
-    
+
     disk::hho_degree_info chdi( disk::priv::hdi_named_args {
                                   .rd = (size_t) order,
                                   .cd = (size_t) order,
@@ -365,7 +365,7 @@ vector_wave_solver(Mesh<T,3,Storage>& msh, size_t order, T alpha, size_t num_ts)
         auto fcs = faces(msh, cl);
         auto num_faces = fcs.size();
         auto fbs = disk::vector_basis_size(chdi.face_degree(), 2, 2);
-        auto nfd = num_faces * fbs; 
+        auto nfd = num_faces * fbs;
         auto CR = disk::curl_reconstruction_pk(msh, cl, chdi);
         Matrix<T, Dynamic, Dynamic> ST = disk::curl_hdg_stabilization(msh, cl, chdi);
 
@@ -377,7 +377,7 @@ vector_wave_solver(Mesh<T,3,Storage>& msh, size_t order, T alpha, size_t num_ts)
 
         Matrix<T, Dynamic, Dynamic> mA = beta*dt2*mK + gamma*dt*mC + mM;
         Matrix<T, Dynamic, Dynamic> mP = 2.0*mM + (2.0*gamma - 1.0)*dt*mC - (0.5 - 2.0*beta + gamma)*dt2*mK;
-        Matrix<T, Dynamic, Dynamic> mQ = -mM - (gamma - 1.0)*dt*mC - (0.5 + beta - gamma)*dt2*mK; 
+        Matrix<T, Dynamic, Dynamic> mQ = -mM - (gamma - 1.0)*dt*mC - (0.5 + beta - gamma)*dt2*mK;
 
         Matrix<T, Dynamic, Dynamic> mA_TT = mA.block(  0,   0, cbs, cbs);
         Matrix<T, Dynamic, Dynamic> mA_TF = mA.block(  0, cbs, cbs, nfd);
@@ -448,7 +448,7 @@ vector_wave_solver(Mesh<T,3,Storage>& msh, size_t order, T alpha, size_t num_ts)
     std::chrono::duration<double> elapsed_fact = fact_end - fact_start;
     std::cout << sgr::Yellowfg << "Factorization time: ";
     std::cout << elapsed_fact.count() << " s\n";
-    
+
     disk::dynamic_vector<T> RHS = disk::dynamic_vector<T>::Zero(assm.RHS.rows());
     disk::dynamic_vector<T> X_prev = assm.X_prev;
     disk::dynamic_vector<T> X_curr = assm.X_curr;
@@ -471,11 +471,11 @@ vector_wave_solver(Mesh<T,3,Storage>& msh, size_t order, T alpha, size_t num_ts)
             //+ dt2*(0.5 + gamma - 2.0*beta)*B_curr
             //+ dt2*(0.5 - gamma + beta)*B_prev;
 
-        
+
         X_next = solver.solve(RHS);
         auto bs_end = std::chrono::steady_clock::now();
 
-        X_next_decond = assm.Decond1*X_curr 
+        X_next_decond = assm.Decond1*X_curr
                       + assm.Decond2*X_prev
                       + assm.Decond3*X_next;
         auto decond_end = std::chrono::steady_clock::now();
@@ -627,7 +627,8 @@ int main(int argc, char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.mesh$") ))
     {
         std::cout << "Guessed mesh format: Netgen 3D" << std::endl;
-        auto msh = disk::load_netgen_3d_mesh<coord_T>(mesh_filename);
+        disk::simplicial_mesh<T, 3> msh;
+        disk::load_mesh_netgen<coord_T>(mesh_filename, msh);
 
         vector_wave_solver(msh, degree, stab_param, num_ts);
 
@@ -642,7 +643,7 @@ int main(int argc, char **argv)
         disk::gmsh_geometry_loader< disk::simplicial_mesh<T,3> > loader;
         //disk::generic_mesh<T,3> msh;
         //disk::gmsh_geometry_loader< disk::generic_mesh<T,3> > loader;
-        
+
         loader.read_mesh(mesh_filename);
         loader.populate_mesh(msh);
 
@@ -656,11 +657,10 @@ int main(int argc, char **argv)
     {
         std::cout << "Guessed mesh format: FVCA6 3D" << std::endl;
         disk::generic_mesh<coord_T,3> msh;
-        
         disk::load_mesh_fvca6_3d<coord_T>(mesh_filename, msh);
-        
+
         vector_wave_solver(msh, degree, stab_param, num_ts);
-        
+
         return 0;
     }
 
diff --git a/apps/maxwell/src/rectest.cpp b/apps/maxwell/src/rectest.cpp
index 7f5e25ee..c461bed0 100644
--- a/apps/maxwell/src/rectest.cpp
+++ b/apps/maxwell/src/rectest.cpp
@@ -104,7 +104,7 @@ void rectest(Mesh& msh, size_t order)
         //std::cout << " v'Sv = " << vSv << " " << std::sqrt(vSv) << std::endl;
 
         JacobiSVD<MatrixXd> svd(GR.second + ST);
-        double cond = svd.singularValues()(0) 
+        double cond = svd.singularValues()(0)
             / svd.singularValues()(svd.singularValues().size()-1);
 
         //std::cout << cond << std::endl;
@@ -116,7 +116,7 @@ void rectest(Mesh& msh, size_t order)
         auto CR = disk::curl_reconstruction(msh, cl, chdi);
 
         Matrix<T, Dynamic, 1> pfun = disk::project_tangent(msh, cl, chdi, fun, 1);
-        
+
         Matrix<T, Dynamic, 1> rfun = Matrix<T, Dynamic, 1>::Zero(rb.size());
         rfun.segment(3, rb.size()-3) = CR.first * pfun;
 
@@ -176,7 +176,7 @@ int main(int argc, char **argv)
             case 'k':
                 degree = std::stoi(optarg);
                 break;
-            
+
             case '?':
             default:
                 std::cout << "Invalid option" << std::endl;
@@ -188,7 +188,8 @@ int main(int argc, char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.mesh2d$") ))
     {
         std::cout << "Guessed mesh format: Netgen 2D" << std::endl;
-        auto msh = disk::load_netgen_2d_mesh<T>(mesh_filename);
+        disk::simplicial_mesh<T, 2> msh;
+        disk::load_mesh_netgen<T>(mesh_filename, msh);
 
         rectest(msh, degree);
 
@@ -199,10 +200,11 @@ int main(int argc, char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.quad$") ))
     {
         std::cout << "Guessed mesh format: DiSk++ Cartesian 2D" << std::endl;
-        auto msh = disk::load_cartesian_2d_mesh<T>(mesh_filename);
-        
+        disk::cartesian_mesh<T, 2> msh;
+        disk::load_mesh_diskpp_cartesian<T>(mesh_filename, msh);
+
         rectest(msh, degree);
-        
+
         return 0;
     }
 
@@ -210,7 +212,8 @@ int main(int argc, char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.typ1$") ))
     {
         std::cout << "Guessed mesh format: FVCA5 2D" << std::endl;
-        auto msh = disk::load_fvca5_2d_mesh<T>(mesh_filename);
+        disk::generic_mesh<T, 2> msh;
+        disk::load_mesh_fvca5_2d<T>(mesh_filename, msh);
 
         rectest(msh, degree);
 
@@ -221,7 +224,8 @@ int main(int argc, char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.mesh$") ))
     {
         std::cout << "Guessed mesh format: Netgen 3D" << std::endl;
-        auto msh = disk::load_netgen_3d_mesh<T>(mesh_filename);
+        disk::simplicial_mesh<T, 3> msh;
+        disk::load_mesh_netgen<T>(mesh_filename, msh);
 
         rectest(msh, degree);
 
@@ -232,10 +236,11 @@ int main(int argc, char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.hex$") ))
     {
         std::cout << "Guessed mesh format: DiSk++ Cartesian 3D" << std::endl;
-        auto msh = disk::load_cartesian_3d_mesh<T>(mesh_filename);
-        
+        disk::cartesian_mesh<T, 3> msh;
+        disk::load_mesh_diskpp_cartesian<T>(mesh_filename, msh);
+
         rectest(msh, degree);
-        
+
         return 0;
     }
 
@@ -243,9 +248,11 @@ int main(int argc, char **argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.msh$") ))
     {
         std::cout << "Guessed mesh format: FVCA6 3D" << std::endl;
-        auto msh = disk::load_fvca6_3d_mesh<T>(mesh_filename);
+        disk::generic_mesh<T,3> msh;
+        disk::load_mesh_fvca6_3d<T>(mesh_filename, msh);
+
         rectest(msh, degree);
-        
+
         return 0;
     }
 }
\ No newline at end of file
diff --git a/apps/stokes/stokes_driven_cavity.cpp b/apps/stokes/stokes_driven_cavity.cpp
index 52a86feb..53620f02 100644
--- a/apps/stokes/stokes_driven_cavity.cpp
+++ b/apps/stokes/stokes_driven_cavity.cpp
@@ -376,7 +376,7 @@ auto run_stokes(const Mesh &msh, size_t degree, bool use_sym_grad = true) {
 
 // #if 0
 int main(int argc, char **argv) {
-  using RealType = double;
+  using T = double;
   bool use_sym_grad = false;
 
   char *filename = nullptr;
@@ -409,42 +409,23 @@ int main(int argc, char **argv) {
 
   if (std::regex_match(filename, std::regex(".*\\.typ1$"))) {
     std::cout << "Guessed mesh format: FVCA5 2D" << std::endl;
-
-    typedef disk::generic_mesh<RealType, 2> mesh_type;
-
-    mesh_type msh;
-    disk::fvca5_mesh_loader<RealType, 2> loader;
-    if (!loader.read_mesh(filename)) {
-      std::cout << "Problem loading mesh." << std::endl;
-      return 1;
-    }
-    loader.populate_mesh(msh);
-
+    disk::generic_mesh<T, 2> msh;
+    disk::load_mesh_fvca5_2d<T>(filename, msh);
     run_stokes(msh, degree, use_sym_grad);
-    std::cout << "fini" << std::endl;
   }
 
   if (std::regex_match(filename, std::regex(".*\\.mesh2d$"))) {
     std::cout << "Guessed mesh format: Netgen 2D" << std::endl;
-
-    typedef disk::simplicial_mesh<RealType, 2> mesh_type;
-
-    mesh_type msh;
-    disk::netgen_mesh_loader<RealType, 2> loader;
-    if (!loader.read_mesh(filename)) {
-      std::cout << "Problem loading mesh." << std::endl;
-      return 1;
-    }
-    loader.populate_mesh(msh);
-
+    disk::simplicial_mesh<T, 2> msh;
+    disk::load_mesh_netgen<T>(filename, msh);
     run_stokes(msh, degree, use_sym_grad);
   }
 
   /* Medit 2d*/
   if (std::regex_match(filename, std::regex(".*\\.medit2d$"))) {
     std::cout << "Guessed mesh format: Medit format" << std::endl;
-    typedef disk::generic_mesh<RealType, 2> mesh_type;
-    mesh_type msh = disk::load_medit_2d_mesh<RealType>(filename);
+    disk::generic_mesh<T, 2> msh;
+    disk::load_mesh_medit<T>(filename, msh);
     run_stokes(msh, degree, use_sym_grad);
   }
 
diff --git a/apps/unsteady_laplacian/src/heat_dG.cpp b/apps/unsteady_laplacian/src/heat_dG.cpp
index 9e41c6a4..ae503060 100644
--- a/apps/unsteady_laplacian/src/heat_dG.cpp
+++ b/apps/unsteady_laplacian/src/heat_dG.cpp
@@ -593,7 +593,6 @@ unsteady_laplacian_solver(const Mesh& msh, size_t degree, size_t time_steps, siz
 int main(int argc, char **argv)
 {
     using T = double;
-    disk::cartesian_mesh<T, 2> msh;
 
     size_t      degree = 1;
     size_t      time_degree = 0;
@@ -626,7 +625,8 @@ int main(int argc, char **argv)
 	    }
     }
 
-    msh = load_cartesian_2d_mesh<T>(mesh_filename);
+    disk::cartesian_mesh<T, 2> msh;
+    disk::load_mesh_diskpp_cartesian<T>(mesh_filename, msh);
 
     std::cout << "Mesh loaded ..." << std::endl;
     unsteady_laplacian_solver(msh, degree, N, time_degree);
diff --git a/apps/unsteady_laplacian/src/unsteady_laplacian.cpp b/apps/unsteady_laplacian/src/unsteady_laplacian.cpp
index 4134ce97..d2111c11 100644
--- a/apps/unsteady_laplacian/src/unsteady_laplacian.cpp
+++ b/apps/unsteady_laplacian/src/unsteady_laplacian.cpp
@@ -243,7 +243,6 @@ unsteady_laplacian_solver(const Mesh& msh, size_t degree, typename Mesh::coordin
 int main(int argc, char **argv)
 {
     using T = double;
-    disk::cartesian_mesh<T, 2> msh;
 
     T           dt = 0.01;
     size_t      degree = 1;
@@ -251,7 +250,7 @@ int main(int argc, char **argv)
     int ch;
     while ( (ch = getopt(argc, argv, "k:m:t:")) != -1 )
     {
-	switch(ch)
+	    switch(ch)
         {
             case 'k':
                 degree = std::stoi(optarg);
@@ -268,10 +267,11 @@ int main(int argc, char **argv)
             default:
                 std::cout << "Invalid option" << std::endl;
                 return 1;
-	}
+	    }
     }
 
-    msh = load_cartesian_2d_mesh<T>(mesh_filename);
+    disk::cartesian_mesh<T, 2> msh;
+    disk::load_mesh_diskpp_cartesian<T>(mesh_filename, msh);
 
     std::cout << "Mesh loaded ..." << std::endl;
     unsteady_laplacian_solver(msh, degree, dt, 1.0);
diff --git a/apps/utils/mesh2silo.cpp b/apps/utils/mesh2silo.cpp
index fd18f6c3..bd764cc8 100644
--- a/apps/utils/mesh2silo.cpp
+++ b/apps/utils/mesh2silo.cpp
@@ -38,7 +38,8 @@ int main(int argc, const char *argv[])
     if (std::regex_match(mesh_filename, std::regex(".*\\.typ1$") ))
     {
         std::cout << "Guessed mesh format: FVCA5 2D" << std::endl;
-        auto msh = disk::load_fvca5_2d_mesh<T>(mesh_filename);
+        disk::generic_mesh<T,2> msh;
+        disk::load_mesh_fvca5_2d<T>(mesh_filename, msh);
         export_mesh_to_silo(msh, silo_filename);
         return 0;
     }
@@ -47,7 +48,8 @@ int main(int argc, const char *argv[])
     if (std::regex_match(mesh_filename, std::regex(".*\\.mesh2d$") ))
     {
         std::cout << "Guessed mesh format: Netgen 2D" << std::endl;
-        auto msh = disk::load_netgen_2d_mesh<T>(mesh_filename);
+        disk::simplicial_mesh<T, 2> msh;
+        disk::load_mesh_netgen<T>(mesh_filename, msh);
         export_mesh_to_silo(msh, silo_filename);
         return 0;
     }
@@ -56,7 +58,8 @@ int main(int argc, const char *argv[])
     if (std::regex_match(mesh_filename, std::regex(".*\\.quad$") ))
     {
         std::cout << "Guessed mesh format: DiSk++ Cartesian 2D" << std::endl;
-        auto msh = disk::load_cartesian_2d_mesh<T>(mesh_filename);
+        disk::cartesian_mesh<T, 2> msh;
+        disk::load_mesh_diskpp_cartesian<T>(mesh_filename, msh);
         export_mesh_to_silo(msh, silo_filename);
         return 0;
     }
@@ -66,7 +69,8 @@ int main(int argc, const char *argv[])
     if (std::regex_match(mesh_filename, std::regex(".*\\.mesh$") ))
     {
         std::cout << "Guessed mesh format: Netgen 3D" << std::endl;
-        auto msh = disk::load_netgen_3d_mesh<T>(mesh_filename);
+        disk::simplicial_mesh<T, 3> msh;
+        disk::load_mesh_netgen<T>(mesh_filename, msh);
         export_mesh_to_silo(msh, silo_filename);
         return 0;
     }
@@ -75,7 +79,8 @@ int main(int argc, const char *argv[])
     if (std::regex_match(mesh_filename, std::regex(".*\\.hex$") ))
     {
         std::cout << "Guessed mesh format: DiSk++ Cartesian 3D" << std::endl;
-        auto msh = disk::load_cartesian_3d_mesh<T>(mesh_filename);
+        disk::cartesian_mesh<T, 3> msh;
+        disk::load_mesh_diskpp_cartesian<T>(mesh_filename, msh);
         export_mesh_to_silo(msh, silo_filename);
         return 0;
     }
@@ -89,7 +94,6 @@ int main(int argc, const char *argv[])
         export_mesh_to_silo(msh, silo_filename);
         return 0;
     }
-    
+
     return 0;
 }
-
diff --git a/apps/vector_laplacian/src/vector_laplacian.cpp b/apps/vector_laplacian/src/vector_laplacian.cpp
index 957bd44f..6db07ab6 100755
--- a/apps/vector_laplacian/src/vector_laplacian.cpp
+++ b/apps/vector_laplacian/src/vector_laplacian.cpp
@@ -259,7 +259,8 @@ main(int argc, char** argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.typ1$")))
     {
         std::cout << "Guessed mesh format: FVCA5 2D" << std::endl;
-        auto msh = disk::load_fvca5_2d_mesh<RealType>(mesh_filename);
+        disk::generic_mesh<RealType,2> msh;
+        disk::load_mesh_fvca5_2d<RealType>(mesh_filename, msh);
         run_vector_laplacian_solver(msh, rp, material_data);
         return 0;
     }
@@ -268,7 +269,8 @@ main(int argc, char** argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.medit2d$")))
     {
         std::cout << "Guessed mesh format: Medit format" << std::endl;
-        auto msh = disk::load_medit_2d_mesh<RealType>(mesh_filename);
+        disk::generic_mesh<RealType, 2> msh;
+        disk::load_mesh_medit<RealType>(mesh_filename, msh);
         run_vector_laplacian_solver(msh, rp, material_data);
         return 0;
     }
@@ -277,7 +279,8 @@ main(int argc, char** argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.mesh2d$")))
     {
         std::cout << "Guessed mesh format: Netgen 2D" << std::endl;
-        auto msh = disk::load_netgen_2d_mesh<RealType>(mesh_filename);
+        disk::simplicial_mesh<RealType, 2> msh;
+        disk::load_mesh_netgen<RealType>(mesh_filename, msh);
         run_vector_laplacian_solver(msh, rp, material_data);
         return 0;
     }
@@ -286,7 +289,8 @@ main(int argc, char** argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.quad$")))
     {
         std::cout << "Guessed mesh format: DiSk++ Cartesian 2D" << std::endl;
-        auto msh = disk::load_cartesian_2d_mesh<RealType>(mesh_filename);
+        disk::cartesian_mesh<RealType, 2> msh;
+        disk::load_mesh_diskpp_cartesian<RealType>(mesh_filename, msh);
         run_vector_laplacian_solver(msh, rp, material_data);
         return 0;
     }
@@ -295,7 +299,8 @@ main(int argc, char** argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.mesh$")))
     {
         std::cout << "Guessed mesh format: Netgen 3D" << std::endl;
-        auto msh = disk::load_netgen_3d_mesh<RealType>(mesh_filename);
+        disk::simplicial_mesh<RealType, 3> msh;
+        disk::load_mesh_netgen<RealType>(mesh_filename, msh);
         run_vector_laplacian_solver(msh, rp, material_data);
         return 0;
     }
@@ -304,7 +309,8 @@ main(int argc, char** argv)
     if (std::regex_match(mesh_filename, std::regex(".*\\.hex$")))
     {
         std::cout << "Guessed mesh format: DiSk++ Cartesian 3D" << std::endl;
-        auto msh = disk::load_cartesian_3d_mesh<RealType>(mesh_filename);
+        disk::cartesian_mesh<RealType, 3> msh;
+        disk::load_mesh_diskpp_cartesian<RealType>(mesh_filename, msh);
         run_vector_laplacian_solver(msh, rp, material_data);
         return 0;
     }
diff --git a/apps/vector_laplacian/src/vector_laplacian_test.cpp b/apps/vector_laplacian/src/vector_laplacian_test.cpp
index 9c31202a..ecc9af77 100644
--- a/apps/vector_laplacian/src/vector_laplacian_test.cpp
+++ b/apps/vector_laplacian/src/vector_laplacian_test.cpp
@@ -247,7 +247,8 @@ test_triangles_fvca5(const run_params& rp, const Parameters material_data)
 
     for (size_t i = 0; i < runs; i++)
     {
-        auto msh = disk::load_fvca5_2d_mesh<T>(paths[i].c_str());
+        disk::generic_mesh<T, 2> msh;
+        disk::load_mesh_fvca5_2d<T>(paths[i].c_str(), msh);
         error_sumup.push_back(run_vector_laplacian_solver(msh, rp, material_data));
     }
     printResults(error_sumup);
@@ -270,7 +271,8 @@ test_triangles_netgen(const run_params& rp, const Parameters material_data)
 
     for (size_t i = 0; i < runs; i++)
     {
-        auto msh = disk::load_netgen_2d_mesh<T>(paths[i].c_str());
+        disk::simplicial_mesh<T, 2> msh;
+        disk::load_mesh_netgen<T>(paths[i].c_str(), msh);
         error_sumup.push_back(run_vector_laplacian_solver(msh, rp, material_data));
     }
     printResults(error_sumup);
@@ -293,7 +295,8 @@ test_hexagons(const run_params& rp, const Parameters material_data)
 
     for (size_t i = 0; i < runs; i++)
     {
-        auto msh = disk::load_fvca5_2d_mesh<T>(paths[i].c_str());
+        disk::generic_mesh<T, 2> msh;
+        disk::load_mesh_fvca5_2d<T>(paths[i].c_str(), msh);
         error_sumup.push_back(run_vector_laplacian_solver(msh, rp, material_data));
     }
     printResults(error_sumup);
@@ -316,7 +319,8 @@ test_kershaws(const run_params& rp, const Parameters material_data)
 
     for (size_t i = 0; i < runs; i++)
     {
-        auto msh = disk::load_fvca5_2d_mesh<T>(paths[i].c_str());
+        disk::generic_mesh<T, 2> msh;
+        disk::load_mesh_fvca5_2d<T>(paths[i].c_str(), msh);
         error_sumup.push_back(run_vector_laplacian_solver(msh, rp, material_data));
     }
     printResults(error_sumup);
@@ -339,7 +343,8 @@ test_quads_fvca5(const run_params& rp, const Parameters material_data)
 
     for (size_t i = 0; i < runs; i++)
     {
-        auto msh = disk::load_fvca5_2d_mesh<T>(paths[i].c_str());
+        disk::generic_mesh<T, 2> msh;
+        disk::load_mesh_fvca5_2d<T>(paths[i].c_str(), msh);
         error_sumup.push_back(run_vector_laplacian_solver(msh, rp, material_data));
     }
     printResults(error_sumup);
@@ -362,7 +367,9 @@ test_quads_diskpp(const run_params& rp, const Parameters material_data)
 
     for (size_t i = 0; i < runs; i++)
     {
-        auto msh = disk::load_cartesian_2d_mesh<T>(paths[i].c_str());
+        disk::cartesian_mesh<T, 2> msh;
+        disk::load_mesh_diskpp_cartesian<T>(paths[i].c_str(), msh);
+
         error_sumup.push_back(run_vector_laplacian_solver(msh, rp, material_data));
     }
     printResults(error_sumup);
@@ -385,7 +392,8 @@ test_hexahedra_diskpp(const run_params& rp, const Parameters material_data)
 
     for (int i = 0; i < runs; i++)
     {
-        auto msh = disk::load_cartesian_3d_mesh<T>(paths[i].c_str());
+        disk::cartesian_mesh<T, 3> msh;
+        disk::load_mesh_diskpp_cartesian<T>(paths[i].c_str(), msh);
         error_sumup.push_back(run_vector_laplacian_solver(msh, rp, material_data));
     }
     printResults(error_sumup);
@@ -408,7 +416,9 @@ test_hexahedra_fvca6(const run_params& rp, const Parameters material_data)
 
     for (int i = 0; i < runs; i++)
     {
-        auto msh = disk::load_fvca6_3d_mesh<T>(paths[i].c_str());
+        disk::generic_mesh<T,3> msh;
+        disk::load_mesh_fvca6_3d<T>(paths[i].c_str(), msh);
+
         error_sumup.push_back(run_vector_laplacian_solver(msh, rp, material_data));
     }
     printResults(error_sumup);
@@ -431,7 +441,8 @@ test_tetrahedra_netgen(const run_params& rp, const Parameters material_data)
 
     for (int i = 0; i < runs; i++)
     {
-        auto msh = disk::load_netgen_3d_mesh<T>(paths[i].c_str());
+        disk::simplicial_mesh<T, 3> msh;
+        disk::load_mesh_netgen<T>(paths[i].c_str(), msh);
         error_sumup.push_back(run_vector_laplacian_solver(msh, rp, material_data));
     }
     printResults(error_sumup);
@@ -453,7 +464,8 @@ test_polyhedra_fvca6(const run_params& rp, const Parameters material_data)
 
     for (int i = 0; i < runs; i++)
     {
-        auto msh = disk::load_fvca6_3d_mesh<T>(paths[i].c_str());
+        disk::generic_mesh<T,3> msh;
+        disk::load_mesh_fvca6_3d<T>(paths[i].c_str(), msh);
         error_sumup.push_back(run_vector_laplacian_solver(msh, rp, material_data));
     }
     printResults(error_sumup);
@@ -476,7 +488,8 @@ test_tetrahedra_fvca6(const run_params& rp, const Parameters material_data)
 
     for (int i = 0; i < runs; i++)
     {
-        auto msh = disk::load_fvca6_3d_mesh<T>(paths[i].c_str());
+        disk::generic_mesh<T,3> msh;
+        disk::load_mesh_fvca6_3d<T>(paths[i].c_str(), msh);
         error_sumup.push_back(run_vector_laplacian_solver(msh, rp, material_data));
     }
     printResults(error_sumup);
diff --git a/libdiskpp/include/diskpp/loaders/loader.hpp b/libdiskpp/include/diskpp/loaders/loader.hpp
index d5306ddb..22dcfc0f 100644
--- a/libdiskpp/include/diskpp/loaders/loader.hpp
+++ b/libdiskpp/include/diskpp/loaders/loader.hpp
@@ -1776,6 +1776,7 @@ class medit_mesh_loader<T, 2> : public mesh_loader<generic_mesh<T, 2>>
    }
 
  public:
+   static const char constexpr *expected_extension = "medit2d";
    medit_mesh_loader() = default;
 
    bool
@@ -2077,6 +2078,7 @@ class medit_mesh_loader<T, 3> : public mesh_loader<generic_mesh<T, 3>>
    }
 
  public:
+   static const char constexpr *expected_extension = "medit3d";
    medit_mesh_loader() = default;
 
    bool
@@ -2446,7 +2448,8 @@ class poly_mesh_loader<T, 2> : public mesh_loader<generic_mesh<T, 2>>
         std::vector<size_t> conv_table;
         /* Sort the edges in lexicographical order, remember their original
          * position to convert the pointers in the faces */
-        auto comp_edges = [](const std::pair<size_t, edge_type>& e1, const std::pair<size_t, edge_type>& e2) {
+        auto comp_edges = [](const std::pair<size_t, edge_type>& e1, const std::pair<size_t, edge_type>& e2)
+        {
             return e1.second < e2.second;
         };
         std::sort(m_edges.begin(), m_edges.end(), comp_edges);
@@ -2480,7 +2483,8 @@ class poly_mesh_loader<T, 2> : public mesh_loader<generic_mesh<T, 2>>
 
         /* Sort in lexicographical order and remember original position */
         auto comp_vecs = [](const std::pair<size_t, std::vector<size_t>>& e1,
-                            const std::pair<size_t, std::vector<size_t>>& e2) { return e1.second < e2.second; };
+                            const std::pair<size_t, std::vector<size_t>>& e2)
+                            { return e1.second < e2.second; };
         std::sort(faces_to_edges.begin(), faces_to_edges.end(), comp_vecs);
 
         std::vector<surface_type> faces;
@@ -2722,7 +2726,8 @@ class poly_mesh_loader<T, 3> : public mesh_loader<generic_mesh<T, 3>>
         std::vector<size_t> conv_table;
         /* Sort the edges in lexicographical order, remember their original
          * position to convert the pointers in the faces */
-        auto comp_edges = [](const std::pair<size_t, edge_type>& e1, const std::pair<size_t, edge_type>& e2) {
+        auto comp_edges = [](const std::pair<size_t, edge_type>& e1, const std::pair<size_t, edge_type>& e2)
+        {
             return e1.second < e2.second;
         };
         std::sort(m_edges.begin(), m_edges.end(), comp_edges);
@@ -2745,7 +2750,8 @@ class poly_mesh_loader<T, 3> : public mesh_loader<generic_mesh<T, 3>>
 
         /* Sort in lexicographical order and remember original position */
         auto comp_vecs = [](const std::pair<size_t, std::vector<size_t>>& e1,
-                            const std::pair<size_t, std::vector<size_t>>& e2) { return e1.second < e2.second; };
+                            const std::pair<size_t, std::vector<size_t>>& e2)
+                            { return e1.second < e2.second; };
         std::sort(faces_to_edges.begin(), faces_to_edges.end(), comp_vecs);
 
         std::vector<surface_type> faces;
@@ -2804,6 +2810,7 @@ class poly_mesh_loader<T, 3> : public mesh_loader<generic_mesh<T, 3>>
     }
 };
 
+#if 0
 /* Helper to load uniform 1D meshes. */
 template<typename T>
 [[deprecated("DiSk++ deprecation: The load_mesh_*() functions should be preferred")]]
@@ -2819,6 +2826,7 @@ load_uniform_1d_mesh(T min, T max, size_t cells)
     return msh;
 }
 
+
 /* Helper to load 2D meshes in FVCA5 format */
 template<typename T>
 [[deprecated("DiSk++ deprecation: The load_mesh_*() functions should be preferred")]]
@@ -2916,6 +2924,7 @@ load_cartesian_3d_mesh(const char *filename)
     return msh;
 }
 
+
 /* Helper to load 2D meshes in Medit format */
 template<typename T>
 [[deprecated("DiSk++ deprecation: The load_mesh_*() functions should be preferred")]]
@@ -2947,7 +2956,7 @@ load_medit_3d_mesh(const char* filename)
 
    return msh;
 }
-
+#endif
 /**************************************************************************/
 // New mesh loader helpers
 namespace priv {
@@ -3047,6 +3056,23 @@ load_mesh_poly3d(const char* filename, disk::generic_mesh<T, 3>& msh, bool verbo
     return priv::load_mesh(filename, loader, msh);
 }
 
+
+template<typename T>
+bool
+load_mesh_medit(const char *filename, disk::generic_mesh<T, 2>& msh)
+{
+    disk::medit_mesh_loader<T, 2> loader;
+    return priv::load_mesh(filename, loader, msh);
+}
+
+template<typename T>
+bool
+load_mesh_medit(const char *filename, disk::generic_mesh<T, 3>& msh)
+{
+    disk::medit_mesh_loader<T, 3> loader;
+    return priv::load_mesh(filename, loader, msh);
+}
+
 } // namespace disk
 
 
@@ -3061,21 +3087,21 @@ namespace disk {
 /**
  * @brief Deduce mesh type from filename extension, create appropriate mesh object
  *        and dispatch a function on it
- * 
+ *
  * If you have the function `process_mesh(msh, p1, p2, ..., pn)` you can
- * 
+ *
  *  ```
  *  disk::dispatch_all_meshes(mesh_filename,
  *        [](auto ...args) { process_mesh(args...); },
  *        p1, p2, ..., pn);
  *  ```
- * 
- * @tparam Function 
- * @tparam Args 
+ *
+ * @tparam Function
+ * @tparam Args
  * @param mesh_filename Name of the file containing the mesh.
  * @param func Function to dispatch. The mesh object is passed as first parameter.
  * @param args Arguments to the function to dispatch.
- * @return int 
+ * @return int
  */
 
 template<typename Function, typename... Args>
@@ -3191,4 +3217,3 @@ dispatch_all_meshes(const char *mesh_filename, Function func, Args && ...args)
 }
 
 } // namespace disk
-
diff --git a/libdiskpp/include/diskpp/mesh/meshgen_fvca5hex.hpp b/libdiskpp/include/diskpp/mesh/meshgen_fvca5hex.hpp
index cf165980..49b37408 100644
--- a/libdiskpp/include/diskpp/mesh/meshgen_fvca5hex.hpp
+++ b/libdiskpp/include/diskpp/mesh/meshgen_fvca5hex.hpp
@@ -16,6 +16,10 @@
 #include <vector>
 #include <map>
 
+#include <cassert>
+#include <cmath>
+#include <algorithm>
+
 namespace disk::fvca5 {
 
 enum class boundary {
@@ -28,22 +32,22 @@ enum class boundary {
 struct point {
     size_t x;
     size_t y;
-    
+
     point(size_t px, size_t py)
         : x(px), y(py)
     {}
-    
+
     void make_absolute(size_t base_x, size_t base_y)
     {
         x += base_x;
         y += base_y;
     }
-    
+
     bool operator<(const point& other) const
     {
         return (x < other.x) or (x == other.x and y < other.y);
     }
-    
+
     bool operator==(const point& other) const
     {
         return x == other.x and y == other.y;
@@ -53,7 +57,7 @@ struct point {
 struct edge {
     point   a;
     point   b;
-    
+
     edge(const point& pa, const point& pb)
         : a(pa), b(pb)
     {
@@ -61,18 +65,18 @@ struct edge {
         if (b < a)
             std::swap(a, b);
     }
-    
+
     void make_absolute(size_t base_x, size_t base_y)
     {
         a.make_absolute(base_x, base_y);
         b.make_absolute(base_x, base_y);
     }
-    
+
     bool operator<(const edge& other) const
     {
         return (a < other.a) or (a == other.a and b < other.b);
     }
-    
+
     bool operator==(const edge& other) const
     {
         return a == other.a and b == other.b;
@@ -82,24 +86,24 @@ struct edge {
 struct polygon {
     std::vector<point>  points;
     size_t              domain_id;
-    
+
     void add_point(const point& pt) {
         points.push_back(pt);
     }
-    
+
     void make_absolute(size_t base_x, size_t base_y)
     {
         for (auto& point : points)
             point.make_absolute(base_x, base_y);
     }
-    
+
     bool operator<(const polygon& other) const {
         if (points.size() < other.points.size())
             return true;
-        
+
         if (points.size() > other.points.size())
             return false;
-        
+
         return std::lexicographical_compare(points.begin(), points.end(),
             other.points.begin(), other.points.end());
     }
@@ -122,7 +126,7 @@ void add_bottom(std::vector<polygon>& local_polygons,
 {
     bndedges[boundary::BOTTOM].push_back( edge(point(0,0), point(3,0)) );
     bndedges[boundary::BOTTOM].push_back( edge(point(3,0), point(4,0)) );
-    
+
     polygon bottom_quad;
     bottom_quad.add_point( point(0,0) );
     bottom_quad.add_point( point(3,0) );
@@ -136,7 +140,7 @@ void add_top(std::vector<polygon>& local_polygons,
 {
     bndedges[boundary::TOP].push_back( edge(point(0,4), point(3,4)) );
     bndedges[boundary::TOP].push_back( edge(point(3,4), point(4,4)) );
-    
+
     polygon top_quad;
     top_quad.add_point( point(0,4) );
     top_quad.add_point( point(1,3) );
@@ -162,13 +166,13 @@ void add_left(std::vector<polygon>& local_polygons,
 {
     bndedges[boundary::LEFT].push_back( edge(point(0,0), point(0,2)) );
     bndedges[boundary::LEFT].push_back( edge(point(0,2), point(0,4)) );
-    
+
     polygon tri_lower;
     tri_lower.add_point( point(0,0) );
     tri_lower.add_point( point(1,1) );
     tri_lower.add_point( point(0,2) );
     local_polygons.push_back(tri_lower);
-    
+
     polygon tri_upper;
     tri_upper.add_point( point(0,2) );
     tri_upper.add_point( point(1,3) );
@@ -181,7 +185,7 @@ void add_right(std::vector<polygon>& local_polygons,
 {
     bndedges[boundary::RIGHT].push_back( edge(point(4,0), point(4,2)) );
     bndedges[boundary::RIGHT].push_back( edge(point(4,2), point(4,4)) );
-    
+
     polygon penta_lower;
     penta_lower.add_point( point(4,0) );
     penta_lower.add_point( point(4,2) );
@@ -189,7 +193,7 @@ void add_right(std::vector<polygon>& local_polygons,
     penta_lower.add_point( point(2,1) );
     penta_lower.add_point( point(3,0) );
     local_polygons.push_back(penta_lower);
-    
+
     polygon penta_upper;
     penta_upper.add_point( point(4,2) );
     penta_upper.add_point( point(4,4) );
@@ -209,7 +213,7 @@ void add_right(std::vector<polygon>& local_polygons)
     hex_lower.add_point( point(2,1) );
     hex_lower.add_point( point(3,0) );
     local_polygons.push_back(hex_lower);
-    
+
     polygon hex_upper;
     hex_upper.add_point( point(4,2) );
     hex_upper.add_point( point(5,3) ); // on the neighboring tile
@@ -226,7 +230,7 @@ struct hexmesh {
     std::vector<edge>                       edges;
     std::vector<polygon>                    polygons;
     std::map<boundary, std::vector<edge>>   boundary_edges;
-    
+
     using optn = std::optional<size_t>;
     using npair = std::pair<optn, optn>;
     std::map<edge, npair>                   edge_owners;
@@ -241,7 +245,7 @@ update_geometry(hexmesh& hm, const std::vector<polygon>& polys,
         auto end = hm.boundary_edges[boundary].end();
         hm.boundary_edges[boundary].insert(end, edges.begin(), edges.end());
     }
-    
+
     for (const auto& poly : polys)
     {
         const auto& points = poly.points;
@@ -267,37 +271,37 @@ void hexmesh_work(hexmesh& hm, size_t xmin, size_t xmax, size_t ymin, size_t yma
 
         std::vector<polygon>                    local_polygons;
         std::map<boundary, std::vector<edge>>   local_boundary_edges;
-        
+
         make_main_hex(local_polygons);
 
         if (bottom)
             add_bottom(local_polygons, local_boundary_edges);
-        
+
         if (top)
             add_top(local_polygons, local_boundary_edges);
         else
             add_top(local_polygons);
-        
+
         if (left)
             add_left(local_polygons, local_boundary_edges);
-        
+
         if (right)
             add_right(local_polygons, local_boundary_edges);
         else
             add_right(local_polygons);
-        
+
         for(auto& polygon : local_polygons)
             polygon.make_absolute(xmin, ymin);
-        
+
         for(auto& [boundary,edges] : local_boundary_edges)
             for (auto& edge : edges)
                 edge.make_absolute(xmin, ymin);
-        
+
         update_geometry(hm, local_polygons, local_boundary_edges);
-        
+
         return;
     }
-    
+
     hexmesh_work(hm, xmin, (xmax+xmin)/2, ymin, (ymax+ymin)/2);
     hexmesh_work(hm, (xmax+xmin)/2, xmax, ymin, (ymax+ymin)/2);
     hexmesh_work(hm, xmin, (xmax+xmin)/2, (ymax+ymin)/2, ymax);
@@ -318,18 +322,18 @@ void make_hexmesh(hexmesh& hm, size_t l)
     hm.coord_max = ipow(2,l+2);
     std::cout << "coord_max: " << hm.coord_max << std::endl;
     hexmesh_work(hm, 0, hm.coord_max, 0, hm.coord_max);
-    
+
     std::cout << "Sorting" << std::endl;
     std::sort(hm.points.begin(), hm.points.end());
     hm.points.erase(
         std::unique(hm.points.begin(), hm.points.end()), hm.points.end() );
-    
+
     std::sort(hm.edges.begin(), hm.edges.end());
     hm.edges.erase(
         std::unique(hm.edges.begin(), hm.edges.end()), hm.edges.end() );
-    
+
     std::sort(hm.polygons.begin(), hm.polygons.end());
-    
+
     auto add_edge_owner = [&](const edge& edge, size_t polyofs) {
         if (hm.edge_owners[edge].first)
         {
@@ -341,7 +345,7 @@ void make_hexmesh(hexmesh& hm, size_t l)
         else
             hm.edge_owners[edge].first = polyofs;
     };
-    
+
     for (size_t polyofs = 0; polyofs < hm.polygons.size(); polyofs++)
     {
         const auto& polygon = hm.polygons[polyofs];
@@ -353,7 +357,7 @@ void make_hexmesh(hexmesh& hm, size_t l)
             add_edge_owner( edge(pa, pb), polyofs);
         }
     }
-    
+
     std::cout << "Points: " << hm.points.size() << ", Edges: ";
     std::cout << hm.edges.size() << ", Polygons: " << hm.polygons.size();
     std::cout << std::endl;
@@ -367,10 +371,10 @@ offset(const hexmesh& hm, const point& pt, bool one_based)
     auto ipt = std::lower_bound(hm.points.begin(), hm.points.end(), pt);
     if (ipt == hm.points.end() or not(*ipt == pt))
         throw std::logic_error("point not found, this is a bug.");
-    
+
     if (one_based)
         return std::distance(hm.points.begin(), ipt) + 1;
-    
+
     return std::distance(hm.points.begin(), ipt);
 }