diff --git a/meshes/d3d/d3d-coreMesh-gold.osh/0.osh b/meshes/d3d/d3d-coreMesh-gold.osh/0.osh
index c2c5d5171..d7cf10fad 100644
Binary files a/meshes/d3d/d3d-coreMesh-gold.osh/0.osh and b/meshes/d3d/d3d-coreMesh-gold.osh/0.osh differ
diff --git a/meshes/d3d/d3d-coreMesh-numbered-gold.osh/0.osh b/meshes/d3d/d3d-coreMesh-numbered-gold.osh/0.osh
index 8d23cf490..e202f6122 100644
Binary files a/meshes/d3d/d3d-coreMesh-numbered-gold.osh/0.osh and b/meshes/d3d/d3d-coreMesh-numbered-gold.osh/0.osh differ
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 06dd75828..dc93ab6e2 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -62,6 +62,7 @@ set(Omega_h_SOURCES
   Omega_h_math_lang.cpp
   Omega_h_matrix.cpp
   Omega_h_mesh.cpp
+  Omega_h_mixedMesh.cpp
   Omega_h_rcFields.cpp
   Omega_h_metric.cpp
   Omega_h_metric_input.cpp
@@ -413,11 +414,6 @@ function(test_func_impl TEST_NAME NUM_PROCS)
       ${CMAKE_SOURCE_DIR}/meshes)
 
     osh_add_exe(periodic_test)
-    set(TEST_EXES ${TEST_EXES} periodic_test)
-    test_func(periodic_test 2 ./periodic_test
-      ${CMAKE_SOURCE_DIR}/meshes/wedge_matchZ_12elem.sms
-      ${CMAKE_SOURCE_DIR}/meshes/wedge_match.smd
-      ${CMAKE_SOURCE_DIR}/meshes/wedge_matchZ_12elem_sync_2.osh 2)
   endif()
 
   test_basefunc(run_arrayops 1 ./arrayops_test)
@@ -636,6 +632,7 @@ set(Omega_h_HEADERS
   Omega_h_math_lang.hpp
   Omega_h_matrix.hpp
   Omega_h_mesh.hpp
+  Omega_h_mixedMesh.hpp
   Omega_h_metric.hpp
   Omega_h_mpi.h
   Omega_h_owners.hpp
diff --git a/src/Omega_h_class.cpp b/src/Omega_h_class.cpp
index a8bac30dc..969d7a80a 100644
--- a/src/Omega_h_class.cpp
+++ b/src/Omega_h_class.cpp
@@ -190,32 +190,4 @@ void classify_equal_order(
   mesh->add_tag<ClassId>(ent_dim, "class_id", 1, class_id);
 }
 
-/* TODO this function is included as work in progress and is not tested */
-void classify_equal_order(
-    Mesh* mesh, Topo_type ent_type, LOs eqv2v, Read<ClassId> eq_class_ids) {
-  LOs eq2e;
-  if (ent_type == Topo_type::vertex) {
-    eq2e = eqv2v;
-  } 
-  else {
-    Write<LO> eq2e_w;
-    Write<I8> codes;
-    auto ev2v = mesh->ask_verts_of(ent_type);
-    auto v2e = mesh->ask_up(Topo_type::vertex, ent_type);
-    find_matches(ent_type, eqv2v, ev2v, v2e, &eq2e_w, &codes);
-    eq2e = eq2e_w;
-  }
-
-  auto ent_dim = mesh->ent_dim(ent_type);
-
-  auto neq = eqv2v.size() / (ent_dim + 1);
-  auto eq_class_dim = Read<I8>(neq, I8(ent_dim));
-  auto class_dim =
-      map_onto(eq_class_dim, eq2e, mesh->nents(ent_type), I8(mesh->dim()), 1);
-  auto class_id = map_onto(eq_class_ids, eq2e, mesh->nents(ent_type), -1, 1);
-  mesh->add_tag<I8>(ent_type, "class_dim", 1, class_dim);
-  mesh->add_tag<ClassId>(ent_type, "class_id", 1, class_id);
-}
-/* */
-
 }  // end namespace Omega_h
diff --git a/src/Omega_h_file.hpp b/src/Omega_h_file.hpp
index 7f3755047..e333a382a 100644
--- a/src/Omega_h_file.hpp
+++ b/src/Omega_h_file.hpp
@@ -10,6 +10,7 @@
 #include <Omega_h_defines.hpp>
 #include <Omega_h_filesystem.hpp>
 #include <Omega_h_mesh.hpp>
+#include <Omega_h_mixedMesh.hpp>
 #include <Omega_h_tag.hpp>
 
 namespace Omega_h {
@@ -31,6 +32,13 @@ void write_sol(Mesh* mesh, std::string const& filepath,
 
 #ifdef OMEGA_H_USE_SIMMODSUITE
 namespace meshsim {
+/**
+ * Return if the mesh is mixed or mono topology
+ * @param[in] mesh path to Simmetrix .sms mesh file
+ * @param[in] model path to Simmetrix GeomSim .smd model file
+ */
+bool isMixed(filesystem::path const& mesh, filesystem::path const& model);
+
 /**
  * Convert a serial Simmetrix sms mesh classified on the specified model to an
  * Omega_h mesh instance.
@@ -41,7 +49,8 @@ namespace meshsim {
 Mesh read(filesystem::path const& mesh, filesystem::path const& model,
           CommPtr comm);
 /**
- * Convert a serial Simmetrix sms mesh classified on the specified model to an
+ * Convert a mono topology (i.e., all tets, all triangles, all hex, etc.)
+ * serial Simmetrix sms mesh classified on the specified model to an
  * Omega_h mesh instance and attach a Simmetrix MeshNex vertex numbering.
  * @param[in] mesh path to Simmetrix .sms mesh file
  * @param[in] model path to Simmetrix GeomSim .smd model file
@@ -50,6 +59,15 @@ Mesh read(filesystem::path const& mesh, filesystem::path const& model,
  */
 Mesh read(filesystem::path const& mesh, filesystem::path const& model,
           filesystem::path const& numbering, CommPtr comm);
+/**
+ * Convert a mixed topology (i.e., tet + wedge) serial Simmetrix sms mesh
+ * classified on the specified model to an Omega_h mesh instance.
+ * @param[in] mesh path to Simmetrix .sms mesh file
+ * @param[in] model path to Simmetrix GeomSim .smd model file
+ * @param[in] comm path to Omega_h communicator instance
+ */
+MixedMesh readMixed(filesystem::path const& mesh, filesystem::path const& model,
+          CommPtr comm);
 void matchRead(filesystem::path const& mesh_fname, filesystem::path const& model,
                CommPtr comm, Mesh *mesh, I8 is_in);
 }  // namespace meshsim
@@ -125,7 +143,7 @@ void write_vtu(std::string const& filename, Mesh* mesh, Int cell_dim,
 void write_vtu(std::string const& filename, Mesh* mesh,
     bool compress = OMEGA_H_DEFAULT_COMPRESS);
 
-void write_vtu(filesystem::path const& filename, Mesh* mesh, Topo_type max_type,
+void write_vtu(filesystem::path const& filename, MixedMesh* mesh, Topo_type max_type,
     bool compress = OMEGA_H_DEFAULT_COMPRESS);
 
 void write_parallel(filesystem::path const& path, Mesh* mesh, Int cell_dim,
diff --git a/src/Omega_h_mesh.cpp b/src/Omega_h_mesh.cpp
index 10f7d64e5..785420c2b 100644
--- a/src/Omega_h_mesh.cpp
+++ b/src/Omega_h_mesh.cpp
@@ -17,8 +17,6 @@
 #include "Omega_h_quality.hpp"
 #include "Omega_h_shape.hpp"
 #include "Omega_h_timer.hpp"
-#include "Omega_h_int_scan.hpp"
-#include "Omega_h_atomics.hpp"
 #include "Omega_h_reduce.hpp"
 #include "Omega_h_print.hpp"
 #include "Omega_h_dbg.hpp"
@@ -29,7 +27,6 @@ Mesh::Mesh() {
   family_ = OMEGA_H_SIMPLEX;
   dim_ = -1;
   for (Int i = 0; i <= 3; ++i) nents_[i] = -1;
-  for (Int i = 0; i <= 7; ++i) nents_type_[i] = -1;
   parting_ = -1;
   nghost_layers_ = -1;
   library_ = nullptr;
@@ -88,27 +85,8 @@ void Mesh::set_dim(Int dim_in) {
   dim_ = dim_in;
 }
 
-Int Mesh::ent_dim(Topo_type ent_type) const {
-  Int ent_dim;
-  if (int(ent_type) == 0) {
-    ent_dim = 0;
-  }
-  else if (int(ent_type) == 1) {
-    ent_dim = 1;
-  }
-  else if ((int(ent_type) > 1) && (int(ent_type) < 4)) {      
-    ent_dim = 2;
-  }
-  else {
-    ent_dim = 3;
-  }
-  return ent_dim;
-}
-
 void Mesh::set_verts(LO nverts_in) { nents_[VERT] = nverts_in; }
 
-void Mesh::set_verts_type(LO nverts_in) { nents_type_[int(Topo_type::vertex)] = nverts_in; }
-
 void Mesh::set_ents(Int ent_dim, Adj down) {
   OMEGA_H_TIME_FUNCTION;
   check_dim(ent_dim);
@@ -143,18 +121,6 @@ LOs Mesh::get_model_matches(Int ent_dim) {
   return model_matches_[ent_dim];
 }
 
-void Mesh::set_ents(Topo_type high_type, Topo_type low_type, Adj h2l) {
-  OMEGA_H_TIME_FUNCTION;
-  check_type(high_type);
-  check_type(low_type);
-  if (int(high_type) < 6) {
-    OMEGA_H_CHECK(!has_ents(high_type));
-  }
-  auto deg = element_degree(high_type, low_type);
-  nents_type_[int(high_type)] = divide_no_remainder(h2l.ab2b.size(), deg);
-  add_adj(high_type, low_type, h2l);
-}
-
 void Mesh::set_parents(Int ent_dim, Parents parents) {
   check_dim2(ent_dim);
   parents_[ent_dim] = std::make_shared<Parents>(parents);
@@ -167,11 +133,6 @@ LO Mesh::nents(Int ent_dim) const {
   return nents_[ent_dim];
 }
 
-LO Mesh::nents(Topo_type ent_type) const {
-  check_type2(ent_type);
-  return nents_type_[int(ent_type)];
-}
-
 LO Mesh::nelems() const { return nents(dim()); }
 
 LO Mesh::nregions() const { return nents(REGION); }
@@ -182,34 +143,6 @@ LO Mesh::nedges() const { return nents(EDGE); }
 
 LO Mesh::nverts() const { return nents(VERT); }
 
-LO Mesh::npyrams() const { return nents(Topo_type::pyramid); }
-
-LO Mesh::nwedges() const { return nents(Topo_type::wedge); }
-
-LO Mesh::nhexs() const { return nents(Topo_type::hexahedron); }
-
-LO Mesh::ntets() const { return nents(Topo_type::tetrahedron); }
-
-LO Mesh::nquads() const { return nents(Topo_type::quadrilateral); }
-
-LO Mesh::ntris() const { return nents(Topo_type::triangle); }
-
-LO Mesh::nedges_mix() const { return nents(Topo_type::edge); }
-
-LO Mesh::nverts_mix() const { return nents(Topo_type::vertex); }
-
-LO Mesh::nregions_mix() const { 
-  return (nents(Topo_type::tetrahedron) +
-          nents(Topo_type::hexahedron) +
-          nents(Topo_type::wedge) +
-          nents(Topo_type::pyramid));
-}
-
-LO Mesh::nfaces_mix() const { 
-  return (nents(Topo_type::triangle) +
-          nents(Topo_type::quadrilateral));
-}
-
 GO Mesh::nglobal_ents(Int ent_dim) {
   if (!could_be_shared(ent_dim)) {
     return comm_->allreduce(GO(nents(ent_dim)), OMEGA_H_SUM);
@@ -223,18 +156,6 @@ void Mesh::add_tag(Int ent_dim, std::string const& name, Int ncomps) {
   this->add_tag(ent_dim, name, ncomps, Read<T>(), true);
 }
 
-template <typename T>
-void Mesh::add_tag(Topo_type ent_type, std::string const& name, Int ncomps) {
-  if (has_tag(ent_type, name)) remove_tag(ent_type, name);
-  check_type2(ent_type);
-  check_tag_name(name);
-  OMEGA_H_CHECK(ncomps >= 0);
-  OMEGA_H_CHECK(ncomps <= Int(INT8_MAX));
-  OMEGA_H_CHECK(tags_type_[int(ent_type)].size() < size_t(INT8_MAX));
-  TagPtr ptr(new Tag<T>(name, ncomps));
-  tags_type_[int(ent_type)].push_back(std::move(ptr));
-}
-
 template <typename T>
 void Mesh::add_tag(Int ent_dim, std::string const& name, Int ncomps,
     Read<T> array, bool internal) {
@@ -264,28 +185,6 @@ void Mesh::add_tag(Int ent_dim, std::string const& name, Int ncomps,
   if (!internal) react_to_set_tag(ent_dim, name);
 }
 
-template <typename T>
-void Mesh::add_tag(Topo_type ent_type, std::string const& name, Int ncomps,
-    Read<T> array, bool internal) {
-  check_type2(ent_type);
-  auto it = tag_iter(ent_type, name);
-  auto had_tag = (it != tags_type_[int(ent_type)].end());
-  auto ptr = std::make_shared<Tag<T>>(name, ncomps);
-  ptr->set_array(array);
-  if (had_tag) {
-    OMEGA_H_CHECK(ncomps == ptr->ncomps());
-    *it = std::move(ptr);
-  } else {
-    check_tag_name(name);
-    OMEGA_H_CHECK(ncomps >= 0);
-    OMEGA_H_CHECK(ncomps <= Int(INT8_MAX));
-    OMEGA_H_CHECK(tags_type_[int(ent_type)].size() < size_t(INT8_MAX));
-    tags_type_[int(ent_type)].push_back(std::move(ptr));
-  }
-  OMEGA_H_CHECK(array.size() == nents_type_[int(ent_type)] * ncomps);
-  if (!internal) react_to_set_tag(ent_type, name);
-}
-
 template <typename T>
 void Mesh::set_tag(
     Int ent_dim, std::string const& name, Read<T> array, bool internal) {
@@ -296,17 +195,6 @@ void Mesh::set_tag(
   this->add_tag(ent_dim, name, divide_no_remainder(array.size(), nents(ent_dim)), array, internal);
 }
 
-template <typename T>
-void Mesh::set_tag(
-    Topo_type ent_type, std::string const& name, Read<T> array, bool internal) {
-  if (!has_tag(ent_type, name)) {
-    Omega_h_fail("set_tag(%s, %s): tag doesn't exist (use add_tag first)\n",
-      dimensional_plural_name(ent_type), name.c_str());
-  }
-  auto it = tag_iter(ent_type,name);
-  this->add_tag(ent_type, name, (*it)->ncomps(), array, internal);
-}
-
 void Mesh::react_to_set_tag(Int ent_dim, std::string const& name) {
   /* hardcoded cache invalidations */
   bool is_coordinates = (name == "coordinates");
@@ -319,28 +207,6 @@ void Mesh::react_to_set_tag(Int ent_dim, std::string const& name) {
   }
 }
 
-void Mesh::react_to_set_tag(Topo_type ent_type, std::string const& name) {
-  bool is_coordinates = (name == "coordinates");
-  if ((int(ent_type) == 0) && (is_coordinates || (name == "metric"))) {
-    remove_tag(Topo_type::edge, "length");
-
-    remove_tag(Topo_type::pyramid, "quality");
-    remove_tag(Topo_type::wedge, "quality");
-    remove_tag(Topo_type::hexahedron, "quality");
-    remove_tag(Topo_type::tetrahedron, "quality");
-    remove_tag(Topo_type::quadrilateral, "quality");
-    remove_tag(Topo_type::triangle, "quality");
-  }
-  if ((int(ent_type) == 0) && is_coordinates) {
-    remove_tag(Topo_type::pyramid, "size");
-    remove_tag(Topo_type::wedge, "size");
-    remove_tag(Topo_type::hexahedron, "size");
-    remove_tag(Topo_type::tetrahedron, "size");
-    remove_tag(Topo_type::quadrilateral, "size");
-    remove_tag(Topo_type::triangle, "size");
-  }
-}
-
 TagBase const* Mesh::get_tagbase(Int ent_dim, std::string const& name) const {
   check_dim2(ent_dim);
   auto it = tag_iter(ent_dim, name);
@@ -351,71 +217,33 @@ TagBase const* Mesh::get_tagbase(Int ent_dim, std::string const& name) const {
   return it->get();
 }
 
-TagBase const* Mesh::get_tagbase(Topo_type ent_type, std::string const& name) const {
-  check_type2(ent_type);
-  auto it = tag_iter(ent_type, name);
-  if (it == tags_type_[int(ent_type)].end()) {
-    Omega_h_fail("get_tagbase(%s, %s): doesn't exist\n",
-        dimensional_plural_name(ent_type), name.c_str());
-  }
-  return it->get();
-}
-
 template <typename T>
 Tag<T> const* Mesh::get_tag(Int ent_dim, std::string const& name) const {
   return as<T>(get_tagbase(ent_dim, name));
 }
 
-template <typename T>
-Tag<T> const* Mesh::get_tag(Topo_type ent_type, std::string const& name) const {
-  return as<T>(get_tagbase(ent_type, name));
-}
-
 template <typename T>
 Read<T> Mesh::get_array(Int ent_dim, std::string const& name) const {
   return get_tag<T>(ent_dim, name)->array();
 }
 
-template <typename T>
-Read<T> Mesh::get_array(Topo_type ent_type, std::string const& name) const {
-  return get_tag<T>(ent_type, name)->array();
-}
-
 void Mesh::remove_tag(Int ent_dim, std::string const& name) {
   if (!has_tag(ent_dim, name)) return;
   check_dim2(ent_dim);
   tags_[ent_dim].erase(tag_iter(ent_dim, name));
 }
 
-void Mesh::remove_tag(Topo_type ent_type, std::string const& name) {
-  if (!has_tag(ent_type, name)) return;
-  check_type2(ent_type);
-  OMEGA_H_CHECK(has_tag(ent_type, name));
-  tags_type_[int(ent_type)].erase(tag_iter(ent_type, name));
-}
-
 bool Mesh::has_tag(Int ent_dim, std::string const& name) const {
   check_dim(ent_dim);
   if (!has_ents(ent_dim)) return false;
   return tag_iter(ent_dim, name) != tags_[ent_dim].end();
 }
 
-bool Mesh::has_tag(Topo_type ent_type, std::string const& name) const {
-  check_type(ent_type);
-  if (!has_ents(ent_type)) return false;
-  return tag_iter(ent_type, name) != tags_type_[int(ent_type)].end();
-}
-
 Int Mesh::ntags(Int ent_dim) const {
   check_dim2(ent_dim);
   return static_cast<Int>(tags_[ent_dim].size());
 }
 
-Int Mesh::ntags(Topo_type ent_type) const {
-  check_type2(ent_type);
-  return static_cast<Int>(tags_type_[int(ent_type)].size());
-}
-
 TagBase const* Mesh::get_tag(Int ent_dim, Int i) const {
   check_dim2(ent_dim);
   OMEGA_H_CHECK(0 <= i);
@@ -423,35 +251,17 @@ TagBase const* Mesh::get_tag(Int ent_dim, Int i) const {
   return tags_[ent_dim][static_cast<std::size_t>(i)].get();
 }
 
-TagBase const* Mesh::get_tag(Topo_type ent_type, Int i) const {
-  check_type2(ent_type);
-  OMEGA_H_CHECK(0 <= i);
-  OMEGA_H_CHECK(i <= ntags(ent_type));
-  return tags_type_[int(ent_type)][static_cast<std::size_t>(i)].get();
-}
-
 bool Mesh::has_ents(Int ent_dim) const {
   check_dim(ent_dim);
   return nents_[ent_dim] >= 0;
 }
 
-bool Mesh::has_ents(Topo_type ent_type) const {
-  check_type(ent_type);
-  return nents_type_[int(ent_type)] >= 0;
-}
-
 bool Mesh::has_adj(Int from, Int to) const {
   check_dim(from);
   check_dim(to);
   return bool(adjs_[from][to]);
 }
 
-bool Mesh::has_adj(Topo_type from_type, Topo_type to_type) const {
-  check_type(from_type);
-  check_type(to_type);
-  return bool(adjs_type_[int(from_type)][int(to_type)]);
-}
-
 Adj Mesh::get_adj(Int from, Int to) const {
   check_dim2(from);
   check_dim2(to);
@@ -459,27 +269,13 @@ Adj Mesh::get_adj(Int from, Int to) const {
   return *(adjs_[from][to]);
 }
 
-Adj Mesh::get_adj(Topo_type from_type, Topo_type to_type) const {
-  check_type2(from_type);
-  check_type2(from_type);
-  OMEGA_H_CHECK(has_adj(from_type, to_type));
-  return *(adjs_type_[int(from_type)][int(to_type)]);
-}
-
 Adj Mesh::ask_down(Int from, Int to) {
   OMEGA_H_CHECK(to < from);
   return ask_adj(from, to);
 }
 
-Adj Mesh::ask_down(Topo_type from_type, Topo_type to_type) {
-  OMEGA_H_CHECK(int(to_type) < int(from_type));
-  return ask_adj(from_type, to_type);
-}
-
 LOs Mesh::ask_verts_of(Int ent_dim) { return ask_adj(ent_dim, VERT).ab2b; }
 
-LOs Mesh::ask_verts_of(Topo_type ent_type) { return ask_adj(ent_type, Topo_type::vertex).ab2b; }
-
 LOs Mesh::ask_elem_verts() { return ask_verts_of(dim()); }
 
 Adj Mesh::ask_up(Int from, Int to) {
@@ -487,11 +283,6 @@ Adj Mesh::ask_up(Int from, Int to) {
   return ask_adj(from, to);
 }
 
-Adj Mesh::ask_up(Topo_type from_type, Topo_type to_type) {
-  OMEGA_H_CHECK(int(from_type) < int(to_type));
-  return ask_adj(from_type, to_type);
-}
-
 Graph Mesh::ask_star(Int ent_dim) {
   OMEGA_H_CHECK(ent_dim < dim());
   return ask_adj(ent_dim, ent_dim);
@@ -527,16 +318,6 @@ Mesh::TagCIterResult Mesh::rc_tag_iter(Int ent_dim, std::string const& name) con
   return {found, it};
 }
 
-Mesh::TagIter Mesh::tag_iter(Topo_type ent_type, std::string const& name) {
-  return std::find_if(tags_type_[int(ent_type)].begin(), tags_type_[int(ent_type)].end(),
-      [&](TagPtr const& a) { return a->name() == name; });
-}
-
-Mesh::TagCIter Mesh::tag_iter(Topo_type ent_type, std::string const& name) const {
-  return std::find_if(tags_type_[int(ent_type)].begin(), tags_type_[int(ent_type)].end(),
-      [&](TagPtr const& a) { return a->name() == name; });
-}
-
 void Mesh::check_dim(Int ent_dim) const {
   OMEGA_H_CHECK_OP(0, <=, ent_dim);
   OMEGA_H_CHECK_OP(ent_dim, <=, dim());
@@ -547,16 +328,6 @@ void Mesh::check_dim2(Int ent_dim) const {
   OMEGA_H_CHECK(has_ents(ent_dim));
 }
 
-void Mesh::check_type(Topo_type ent_type) const {
-  OMEGA_H_CHECK(Topo_type::vertex <= ent_type);
-  OMEGA_H_CHECK(ent_type <= Topo_type::pyramid);
-}
-
-void Mesh::check_type2(Topo_type ent_type) const {
-  check_type(ent_type);
-  OMEGA_H_CHECK(has_ents(ent_type));
-}
-
 void Mesh::add_adj(Int from, Int to, Adj adj) {
   check_dim2(from);
   check_dim(to);
@@ -582,34 +353,6 @@ void Mesh::add_adj(Int from, Int to, Adj adj) {
   adjs_[from][to] = std::make_shared<Adj>(adj);
 }
 
-void Mesh::add_adj(Topo_type from_type, Topo_type to_type, Adj adj) {
-  check_type2(from_type);
-  check_type(to_type);
-  OMEGA_H_CHECK(adj.ab2b.exists());
-  const int from = int(from_type);
-  const int to = int(to_type); 
-
-  if (to < from) {
-    OMEGA_H_CHECK(!adj.a2ab.exists());                         
-    if (to_type == Topo_type::vertex) {
-      OMEGA_H_CHECK(!adj.codes.exists());
-    } else {
-      OMEGA_H_CHECK(adj.codes.exists());
-    }
-    OMEGA_H_CHECK(                                            
-        adj.ab2b.size() == nents(from_type) * element_degree(from_type, to_type));
-  } else {
-    if (from < to) {
-      OMEGA_H_CHECK(adj.a2ab.exists());                        
-      OMEGA_H_CHECK(adj.codes.exists());
-      OMEGA_H_CHECK(
-          adj.ab2b.size() == nents(to_type) * element_degree(to_type, from_type)); 
-    }
-    OMEGA_H_CHECK(adj.a2ab.size() == nents(from_type) + 1);         
-  }
-  adjs_type_[from][to] = std::make_shared<Adj>(adj);
-}
-
 Adj Mesh::derive_adj(Int from, Int to) {
   OMEGA_H_TIME_FUNCTION;
   check_dim(from);
@@ -650,44 +393,6 @@ Adj Mesh::derive_adj(Int from, Int to) {
   OMEGA_H_NORETURN(Adj());
 }
 
-Adj Mesh::derive_adj(Topo_type from_type, Topo_type to_type) {
-  OMEGA_H_TIME_FUNCTION;
-  check_type(from_type);
-  check_type2(to_type);
-  const int from = int(from_type);
-  const int to = int(to_type);
-  if (from < to) {
-    Adj down = ask_adj(to_type, from_type);
-    Int nlows_per_high = element_degree(to_type, from_type);
-    LO nlows = nents(from_type);
-    Adj up = invert_adj(down, nlows_per_high, nlows, to_type, from_type);
-    return up;
-  }
-  else if (to < from) {
-    OMEGA_H_CHECK(to + 1 < from);
-    Topo_type mid_type;
-    if (to_type == Topo_type::vertex) {
-      mid_type = Topo_type::edge;
-    }
-    else if ((from_type == Topo_type::tetrahedron) ||
-             (from_type == Topo_type::pyramid)) {
-      mid_type = Topo_type::triangle;
-    }
-    else {
-      mid_type = Topo_type::quadrilateral;
-    }
-    Adj h2m = ask_adj(from_type, mid_type);
-    Adj m2l = ask_adj(mid_type, to_type);
-    Adj h2l = transit(h2m, m2l, from_type, to_type, mid_type); 
-    return h2l;
-  }
-  /* todo: add second order adjacency derivation */
-  Omega_h_fail("can't derive adjacency from %s to %s\n",
-      dimensional_plural_name(from_type),
-      dimensional_plural_name(to_type));
-  OMEGA_H_NORETURN(Adj());
-}
-
 Adj Mesh::ask_adj(Int from, Int to) {
   OMEGA_H_TIME_FUNCTION;
   check_dim2(from);
@@ -700,18 +405,6 @@ Adj Mesh::ask_adj(Int from, Int to) {
   return derived;
 }
 
-Adj Mesh::ask_adj(Topo_type from_type, Topo_type to_type) {
-  OMEGA_H_TIME_FUNCTION;
-  check_type2(from_type);
-  check_type2(to_type);
-  if (has_adj(from_type, to_type)) {
-    return get_adj(from_type, to_type);
-  }
-  Adj derived = derive_adj(from_type, to_type);
-  adjs_type_[int(from_type)][int(to_type)] = std::make_shared<Adj>(derived);
-  return derived;
-}
-
 void Mesh::add_coords(Reals array) {
   add_tag<Real>(0, "coordinates", dim(), array);
 }
@@ -723,12 +416,6 @@ void Mesh::set_coords(Reals const& array) {
   set_tag<Real>(VERT, "coordinates", array);
 }
 
-void Mesh::add_coords_mix(Reals array) {
-  add_tag<Real>(Topo_type::vertex, "coordinates", dim(), array);
-}
-
-Reals Mesh::coords_mix() const { return get_array<Real>(Topo_type::vertex, "coordinates"); }
-
 Read<GO> Mesh::globals(Int ent_dim) const {
   return get_array<GO>(ent_dim, "global");
 }
@@ -1362,13 +1049,6 @@ void get_all_dim_tags(Mesh* mesh, Int dim, TagSet* tags) {
   }
 }
 
-void get_all_type_tags(Mesh* mesh, Int dim, Topo_type ent_type, TagSet* tags) {
-  for (Int j = 0; j < mesh->ntags(ent_type); ++j) {
-    auto tagbase = mesh->get_tag(ent_type, j);
-    (*tags)[size_t(dim)].insert(tagbase->name());
-  }
-}
-
 TagSet get_all_mesh_tags(Mesh* mesh) {
   TagSet out;
   for (Int i = 0; i <= mesh->dim(); ++i) {
@@ -1430,22 +1110,13 @@ __host__
 #define OMEGA_H_INST(T)                                                        \
   template Tag<T> const* Mesh::get_tag<T>(Int dim, std::string const& name)    \
       const;                                                                   \
-  template Tag<T> const* Mesh::get_tag<T>(Topo_type ent_type, std::string const& name)    \
-      const;                                                                   \
   template Read<T> Mesh::get_array<T>(Int dim, std::string const& name) const; \
-  template Read<T> Mesh::get_array<T>(Topo_type ent_type, std::string const& name) const; \
   template void Mesh::add_tag<T>(                                              \
       Int dim, std::string const& name, Int ncomps);                           \
-  template void Mesh::add_tag<T>(                                              \
-      Topo_type ent_type, std::string const& name, Int ncomps);                           \
   template void Mesh::add_tag<T>(Int dim, std::string const& name, Int ncomps, \
       Read<T> array, bool internal);                                           \
-  template void Mesh::add_tag<T>(Topo_type ent_type, std::string const& name, Int ncomps, \
-      Read<T> array, bool internal);                                           \
   template void Mesh::set_tag(                                                 \
       Int dim, std::string const& name, Read<T> array, bool internal);         \
-  template void Mesh::set_tag(                                                 \
-      Topo_type ent_type, std::string const& name, Read<T> array, bool internal);         \
   template Read<T> Mesh::sync_array(Int ent_dim, Read<T> a, Int width);        \
   template Future<T> Mesh::isync_array(Int ent_dim, Read<T> a, Int width);     \
   template Read<T> Mesh::owned_array(Int ent_dim, Read<T> a, Int width);       \
diff --git a/src/Omega_h_mesh.hpp b/src/Omega_h_mesh.hpp
index f39e9dd3e..e5ffae1d9 100644
--- a/src/Omega_h_mesh.hpp
+++ b/src/Omega_h_mesh.hpp
@@ -42,9 +42,7 @@ class Mesh {
   void set_matched(I8 is_matched);
   void set_dim(Int dim_in);
   void set_verts(LO nverts_in);
-  void set_verts_type(LO nverts_in);
   void set_ents(Int ent_dim, Adj down);
-  void set_ents(Topo_type high_type, Topo_type low_type, Adj h2l);
   void set_parents(Int ent_dim, Parents parents);
   Library* library() const;
   CommPtr comm() const;
@@ -56,74 +54,37 @@ class Mesh {
   inline Omega_h_Family family() const { return family_; }
   inline I8 is_matched() const { return matched_; }
   LO nents(Int ent_dim) const;
-  LO nents(Topo_type ent_type) const;
-  Int ent_dim(Topo_type ent_type) const;
   LO nelems() const;
   LO nregions() const;
   LO nfaces() const;
   LO nedges() const;
   LO nverts() const;
-
-  LO npyrams() const;
-  LO nwedges() const;
-  LO nhexs() const;
-  LO ntets() const;
-  LO nquads() const;
-  LO ntris() const;
-  LO nregions_mix() const;
-  LO nfaces_mix() const;
-  LO nedges_mix() const;
-  LO nverts_mix() const;
-
   GO nglobal_ents(Int dim);
   template <typename T>
   void add_tag(Int dim, std::string const& name, Int ncomps);
   template <typename T>
-  void add_tag(Topo_type ent_type, std::string const& name, Int ncomps);
-  template <typename T>
   void add_tag(Int dim, std::string const& name, Int ncomps, Read<T> array,
       bool internal = false);
   template <typename T>
-  void add_tag(Topo_type ent_type, std::string const& name, Int ncomps,
-      Read<T> array, bool internal = false);
-  template <typename T>
   void set_tag(
       Int dim, std::string const& name, Read<T> array, bool internal = false);
-  template <typename T>
-  void set_tag(Topo_type ent_type, std::string const& name, Read<T> array,
-      bool internal = false);
   TagBase const* get_tagbase(Int dim, std::string const& name) const;
-  TagBase const* get_tagbase(Topo_type ent_type, std::string const& name) const;
   template <typename T>
   Tag<T> const* get_tag(Int dim, std::string const& name) const;
   template <typename T>
-  Tag<T> const* get_tag(Topo_type ent_type, std::string const& name) const;
-  template <typename T>
   Read<T> get_array(Int dim, std::string const& name) const;
-  template <typename T>
-  Read<T> get_array(Topo_type ent_type, std::string const& name) const;
   void remove_tag(Int dim, std::string const& name);
-  void remove_tag(Topo_type ent_type, std::string const& name);
   bool has_tag(Int dim, std::string const& name) const;
-  bool has_tag(Topo_type ent_type, std::string const& name) const;
   [[nodiscard]] Int ntags(Int dim) const;
   [[nodiscard]] Int nrctags(Int dim) const;
-  [[nodiscard]] Int ntags(Topo_type ent_type) const;
   TagBase const* get_tag(Int dim, Int i) const;
-  TagBase const* get_tag(Topo_type ent_type, Int i) const;
   bool has_ents(Int dim) const;
-  bool has_ents(Topo_type ent_type) const;
   bool has_adj(Int from, Int to) const;
-  bool has_adj(Topo_type from_type, Topo_type to_type) const;
   Adj get_adj(Int from, Int to) const;
-  Adj get_adj(Topo_type from_type, Topo_type to_type) const;
   Adj ask_down(Int from, Int to);
-  Adj ask_down(Topo_type from_type, Topo_type to_type);
   LOs ask_verts_of(Int dim);
-  LOs ask_verts_of(Topo_type ent_type);
   LOs ask_elem_verts();
   Adj ask_up(Int from, Int to);
-  Adj ask_up(Topo_type from_type, Topo_type to_type);
   Graph ask_star(Int dim);
   Graph ask_dual();
 
@@ -258,22 +219,14 @@ class Mesh {
   };
   TagIter tag_iter(Int dim, std::string const& name);
   TagCIter tag_iter(Int dim, std::string const& name) const;
-  TagIter tag_iter(Topo_type ent_type, std::string const& name);
-  TagCIter tag_iter(Topo_type ent_type, std::string const& name) const;
   TagIterResult rc_tag_iter(Int dim, std::string const& name);
   TagCIterResult rc_tag_iter(Int dim, std::string const& name) const;
   void check_dim(Int dim) const;
   void check_dim2(Int dim) const;
-  void check_type(Topo_type ent_type) const;
-  void check_type2(Topo_type ent_type) const;
   void add_adj(Int from, Int to, Adj adj);
-  void add_adj(Topo_type from_type, Topo_type to_type, Adj adj);
   Adj derive_adj(Int from, Int to);
-  Adj derive_adj(Topo_type from_type, Topo_type to_type);
   Adj ask_adj(Int from, Int to);
-  Adj ask_adj(Topo_type from_type, Topo_type to_type);
   void react_to_set_tag(Int dim, std::string const& name);
-  void react_to_set_tag(Topo_type ent_type, std::string const& name);
   Omega_h_Family family_;
   I8 matched_ = -1;
   Int dim_;
@@ -281,13 +234,10 @@ class Mesh {
   Int parting_;
   Int nghost_layers_;
   LO nents_[DIMS];
-  LO nents_type_[TOPO_TYPES];
   TagVector tags_[DIMS];
-  TagVector tags_type_[TOPO_TYPES];
   // rc field tags stored in "rc" format
   TagVector rc_field_tags_[DIMS];
   AdjPtr adjs_[DIMS][DIMS];
-  AdjPtr adjs_type_[TOPO_TYPES][TOPO_TYPES];
   Remotes owners_[DIMS];
   DistPtr dists_[DIMS];
   RibPtr rib_hints_;
@@ -306,8 +256,6 @@ class Mesh {
   void add_coords(Reals array);
   Reals coords() const;
   void set_coords(Reals const& array);
-  void add_coords_mix(Reals array);
-  Reals coords_mix() const;
   Read<GO> globals(Int dim) const;
   Reals ask_lengths();
   Reals ask_qualities();
@@ -418,7 +366,6 @@ Reals average_field(Mesh* mesh, Int dim, Int ncomps, Reals v2x);
 using TagSet = std::array<std::set<std::string>, DIMS>;
 
 void get_all_dim_tags(Mesh* mesh, Int dim, TagSet* tags);
-void get_all_type_tags(Mesh* mesh, Int dim, Topo_type ent_type, TagSet* tags);
 TagSet get_all_mesh_tags(Mesh* mesh);
 void ask_for_mesh_tags(Mesh* mesh, TagSet const& tags);
 
@@ -445,20 +392,12 @@ __host__
       Int dim, std::string const& name) const;                                 \
   extern template Read<T> Mesh::get_array<T>(Int dim, std::string const& name) \
       const;                                                                   \
-  extern template Read<T> Mesh::get_array<T>(                                  \
-      Topo_type ent_type, std::string const& name) const;                      \
   extern template void Mesh::add_tag<T>(                                       \
       Int dim, std::string const& name, Int ncomps);                           \
-  extern template void Mesh::add_tag<T>(                                       \
-      Topo_type ent_type, std::string const& name, Int ncomps);                \
   extern template void Mesh::add_tag<T>(Int dim, std::string const& name,      \
       Int ncomps, Read<T> array, bool internal);                               \
-  extern template void Mesh::add_tag<T>(Topo_type ent_type,                    \
-      std::string const& name, Int ncomps, Read<T> array, bool internal);      \
   extern template void Mesh::set_tag(                                          \
       Int dim, std::string const& name, Read<T> array, bool internal);         \
-  extern template void Mesh::set_tag(Topo_type ent_type,                       \
-      std::string const& name, Read<T> array, bool internal);                  \
   extern template Read<T> Mesh::sync_array(Int ent_dim, Read<T> a, Int width); \
   extern template Future<T> Mesh::isync_array(                                 \
       Int ent_dim, Read<T> a, Int width);                                      \
diff --git a/src/Omega_h_meshsim.cpp b/src/Omega_h_meshsim.cpp
index e2e445b52..50df67a86 100644
--- a/src/Omega_h_meshsim.cpp
+++ b/src/Omega_h_meshsim.cpp
@@ -4,6 +4,7 @@
 #include "Omega_h_map.hpp"
 #include "Omega_h_vector.hpp"
 #include "Omega_h_mesh.hpp"
+#include "Omega_h_mixedMesh.hpp"
 #include "Omega_h_adj.hpp"
 
 #include "SimModel.h"
@@ -37,28 +38,34 @@ namespace Omega_h {
 
 namespace meshsim {
 
-void read_internal(pMesh m, Mesh* mesh, pMeshNex numbering) {
+struct SimMeshInfo {
+  int count_tet;
+  int count_hex;
+  int count_wedge;
+  int count_pyramid;
+  int count_tri;
+  int count_quad;
+  bool is_simplex;
+  bool is_hypercube;
+};
 
-  (void)mesh;
+SimMeshInfo getSimMeshInfo(pMesh m) {
+  SimMeshInfo info = {0,0,0,0,0,0,false,false};
 
   RIter regions = M_regionIter(m);
-  LO count_tet = 0;
-  LO count_hex = 0;
-  LO count_wedge = 0;
-  LO count_pyramid = 0;
   pRegion rgn;
   while ((rgn = (pRegion) RIter_next(regions))) {
     if (R_topoType(rgn) == Rtet) {
-      count_tet += 1;
+      info.count_tet += 1;
     }
     else if (R_topoType(rgn) == Rhex) {
-      count_hex += 1;
+      info.count_hex += 1;
     }
     else if (R_topoType(rgn) == Rwedge) {
-      count_wedge += 1;
+      info.count_wedge += 1;
     }
     else if (R_topoType(rgn) == Rpyramid) {
-      count_pyramid += 1;
+      info.count_pyramid += 1;
     }
     else {
       Omega_h_fail("Region is not tet, hex, wedge, or pyramid \n");
@@ -68,14 +75,12 @@ void read_internal(pMesh m, Mesh* mesh, pMeshNex numbering) {
 
   FIter faces = M_faceIter(m);
   pFace face;
-  int count_tri = 0;
-  int count_quad = 0;
   while ((face = (pFace) FIter_next(faces))) {
     if (F_numEdges(face) == 3) {
-      count_tri += 1;
+      info.count_tri += 1;
     }
     else if (F_numEdges(face) == 4) {
-      count_quad += 1;
+      info.count_quad += 1;
     }
     else {
       Omega_h_fail ("Face is neither tri nor quad \n");
@@ -83,30 +88,30 @@ void read_internal(pMesh m, Mesh* mesh, pMeshNex numbering) {
   }
   FIter_delete(faces);
 
-  bool is_simplex = 0;
-  bool is_hypercube = 0;
-  if (count_hex == 0 && count_wedge == 0 && count_pyramid == 0) {
-    if (count_tet == 0) {
-      if (count_tri > 0) {
-        is_simplex = 1;
+  info.is_simplex = false;
+  info.is_hypercube = false;
+  if (info.count_hex == 0 && info.count_wedge == 0 && info.count_pyramid == 0) {
+    if (info.count_tet == 0) {
+      if (info.count_tri > 0) {
+        info.is_simplex = true;
       }
     }
-    else if (count_tet > 0) {
-      is_simplex = 1;
+    else if (info.count_tet > 0) {
+      info.is_simplex = true;
     }
     else {
       Omega_h_fail ("Invaild topology type\n");
     }
   }
 
-  if (count_tet == 0 && count_wedge == 0 && count_pyramid == 0) {
-    if (count_hex == 0) {
-      if (count_quad > 0) {
-        is_hypercube = 1;
+  if (info.count_tet == 0 && info.count_wedge == 0 && info.count_pyramid == 0) {
+    if (info.count_hex == 0) {
+      if (info.count_quad > 0) {
+        info.is_hypercube = true;
       }
     }
-    else if (count_hex > 0) {
-      is_hypercube = 1;
+    else if (info.count_hex > 0) {
+      info.is_hypercube = true;
     }
     else {
       Omega_h_fail ("Invaild topology type\n");
@@ -114,526 +119,539 @@ void read_internal(pMesh m, Mesh* mesh, pMeshNex numbering) {
   }
 
   fprintf(stderr, "tet=%d, hex=%d, wedge=%d, pyramid=%d\n",
-         count_tet, count_hex, count_wedge, count_pyramid);
-  fprintf(stderr, "tri=%d, quad=%d\n", count_tri, count_quad);
-
-  const int numVtx = M_numVertices(m);
-  const int numEdges = M_numEdges(m);
-  const int numFaces = M_numFaces(m);
-  const int numRegions = M_numRegions(m);
-
-  std::vector<int> rgn_vertices[4];
-  std::vector<int> face_vertices[2];
-  std::vector<int> edge_vertices[1];
-  std::vector<int> ent_class_ids[4];
-  std::vector<int> ent_class_dim[4];
-  std::vector<int> ent_numbering;
-
-  ent_class_ids[0].reserve(numVtx);
-  ent_class_dim[0].reserve(numVtx);
-  ent_class_ids[1].reserve(numEdges);
-  ent_class_dim[1].reserve(numEdges);
-  ent_class_ids[2].reserve(numFaces);
-  ent_class_dim[2].reserve(numFaces);
-  ent_class_ids[3].reserve(numRegions);
-  ent_class_dim[3].reserve(numRegions);
-  if(numbering) {
-    ent_numbering.reserve(numVtx);
-  }
-
-  Int max_dim;
-  if (numRegions) {
-    max_dim = 3;
-  } else if (numFaces) {
-    max_dim = 2;
-  } else if (numEdges) {
-    max_dim = 1;
-  } else {
-    Omega_h_fail("There were no Elements of dimension higher than zero!\n");
-  }
-
-  HostWrite<Real> host_coords(numVtx*max_dim);
-  VIter vertices = M_vertexIter(m);
-  pVertex vtx;
-  LO v = 0;
-  while ((vtx = (pVertex) VIter_next(vertices))) {
-    double xyz[3];
-    V_coord(vtx,xyz);
-    if( max_dim < 3 && xyz[2] != 0 )
-      Omega_h_fail("The z coordinate must be zero for a 2d mesh!\n");
-    for(int j=0; j<max_dim; j++) {
-      host_coords[v * max_dim + j] = xyz[j];
-    }
-    ent_class_ids[0].push_back(classId(vtx));
-    ent_class_dim[0].push_back(classType(vtx));
-    if(numbering) {
-      ent_numbering.push_back(getNumber(numbering,vtx));
-    }
-    ++v;
-  }
-  VIter_delete(vertices);
-
-
-  HostWrite<LO> host_class_ids_vtx(numVtx);
-  HostWrite<I8> host_class_dim_vtx(numVtx);
-  for (int i = 0; i < numVtx; ++i) {
-    host_class_ids_vtx[i] = ent_class_ids[0][static_cast<std::size_t>(i)];
-    host_class_dim_vtx[i] = ent_class_dim[0][static_cast<std::size_t>(i)];
-  }
-
-  mesh->set_dim(max_dim);
-  if (is_simplex || is_hypercube) {
-    if (is_simplex) {
-      mesh->set_family(OMEGA_H_SIMPLEX);
-    }
-    else if (is_hypercube){
-      mesh->set_family(OMEGA_H_HYPERCUBE);
-    }
-    mesh->set_verts(numVtx);
-    mesh->add_coords(host_coords.write());
-    mesh->add_tag<ClassId>(0, "class_id", 1,
-                           Read<ClassId>(host_class_ids_vtx.write()));
-    mesh->add_tag<I8>(0, "class_dim", 1,
-                      Read<I8>(host_class_dim_vtx.write()));
-    if(numbering) {
-      HostWrite<LO> host_numbering_vtx(numVtx);
-      for (int i = 0; i < numVtx; ++i)
-        host_numbering_vtx[i] = ent_numbering[static_cast<std::size_t>(i)];
-      mesh->add_tag<LO>(0, "simNumbering", 1, Read<LO>(host_numbering_vtx.write()));
-    }
-  }
-  else {
-    mesh->set_family(OMEGA_H_MIXED);
-    mesh->set_verts_type(numVtx);
-    mesh->add_coords_mix(host_coords.write());
-    mesh->add_tag<ClassId>(Topo_type::vertex, "class_id", 1,
-                           Read<ClassId>(host_class_ids_vtx.write()));
-    mesh->add_tag<I8>(Topo_type::vertex, "class_dim", 1,
-                      Read<I8>(host_class_dim_vtx.write()));
-    if(numbering) {
-      fprintf(stderr, "Warning: conversion of Simmetrix MeshNex vertex numbering "
-                      "is not yet supported for mixed meshes\n");
-    }
-  }
+         info.count_tet, info.count_hex, info.count_wedge, info.count_pyramid);
+  fprintf(stderr, "tri=%d, quad=%d\n", info.count_tri, info.count_quad);
+  return info;
+}
 
-  edge_vertices[0].reserve(numEdges*2);
-  EIter edges = M_edgeIter(m);
-  pEdge edge;
-  int count_edge = 0;
-  while ((edge = (pEdge) EIter_next(edges))) {
-    double xyz[3];
-    count_edge += 1;
-    for(int j=0; j<2; ++j) {
-      vtx = E_vertex(edge,j);
-      edge_vertices[0].push_back(EN_id(vtx));
+struct SimMeshEntInfo {
+  int maxDim;
+  bool hasNumbering;
+
+  SimMeshEntInfo(std::array<int,4> numEnts, bool hasNumbering_in) {
+    hasNumbering = hasNumbering_in;
+    maxDim = getMaxDim(numEnts);
+  }
+
+  struct EntClass {
+    HostWrite<LO> id;
+    HostWrite<I8> dim;
+    HostWrite<LO> verts;
+  };
+
+  struct VertexInfo {
+    HostWrite<Real> coords;
+    HostWrite<LO> id;
+    HostWrite<I8> dim;
+    HostWrite<LO> numbering;
+  };
+
+  VertexInfo readVerts(pMesh m,pMeshNex numbering) {
+    const int numVtx = M_numVertices(m);
+    VertexInfo vtxInfo;
+    vtxInfo.coords = HostWrite<Real>(numVtx*maxDim);
+    vtxInfo.id = HostWrite<LO>(numVtx);
+    vtxInfo.dim = HostWrite<I8>(numVtx);
+    if(hasNumbering)
+      vtxInfo.numbering = HostWrite<LO>(numVtx);
+
+    VIter vertices = M_vertexIter(m);
+    pVertex vtx;
+    LO v = 0;
+    while ((vtx = (pVertex) VIter_next(vertices))) {
+      double xyz[3];
       V_coord(vtx,xyz);
-    }
-    ent_class_ids[1].push_back(classId(edge));
-    ent_class_dim[1].push_back(classType(edge));
-  }
-  EIter_delete(edges);
-  
-  HostWrite<LO> host_class_ids_edge(numEdges);
-  HostWrite<I8> host_class_dim_edge(numEdges);
-  for (int i = 0; i < numEdges; ++i) {
-    host_class_ids_edge[i] = ent_class_ids[1][static_cast<std::size_t>(i)];
-    host_class_dim_edge[i] = ent_class_dim[1][static_cast<std::size_t>(i)];
-  }
-
-  HostWrite<LO> host_e2v(numEdges*2);
-  for (Int i = 0; i < numEdges; ++i) {
-    for (Int j = 0; j < 2; ++j) {
-      host_e2v[i*2 + j] =
-          edge_vertices[0][static_cast<std::size_t>(i*2 + j)];
-    }
-  }
-  auto ev2v = Read<LO>(host_e2v.write());
-  if (is_simplex || is_hypercube) {
-    mesh->set_ents(1, Adj(ev2v));
-    mesh->add_tag<ClassId>(1, "class_id", 1,
-                           Read<ClassId>(host_class_ids_edge.write()));
-    mesh->add_tag<I8>(1, "class_dim", 1,
-                      Read<I8>(host_class_dim_edge.write()));
-  }
-  else {
-    mesh->set_ents(Topo_type::edge, Topo_type::vertex, Adj(ev2v));
-    mesh->add_tag<ClassId>(Topo_type::edge, "class_id", 1,
-                           Read<ClassId>(host_class_ids_edge.write()));
-    mesh->add_tag<I8>(Topo_type::edge, "class_dim", 1,
-                      Read<I8>(host_class_dim_edge.write()));
-  }
-
-  face_vertices[0].reserve(count_tri*3);
-  face_vertices[1].reserve(count_quad*4);
-  std::vector<int> face_class_ids[2];
-  std::vector<int> face_class_dim[2];
-  face_class_ids[0].reserve(count_tri);
-  face_class_dim[0].reserve(count_tri);
-  face_class_ids[1].reserve(count_quad);
-  face_class_dim[1].reserve(count_quad);
-
-  faces = M_faceIter(m);
-  while ((face = (pFace) FIter_next(faces))) {
-    if (F_numEdges(face) == 3) {
-      pVertex tri_vertex;
-      pPList tri_vertices = F_vertices(face,1);
-      assert (PList_size(tri_vertices) == 3);
-      void *iter = 0;
-      while ((tri_vertex = (pVertex) PList_next(tri_vertices, &iter))) {
-        face_vertices[0].push_back(EN_id(tri_vertex));
+      if( maxDim < 3 && xyz[2] != 0 )
+        Omega_h_fail("The z coordinate must be zero for a 2d mesh!\n");
+      for(int j=0; j<maxDim; j++) {
+        vtxInfo.coords[v * maxDim + j] = xyz[j];
       }
-      PList_delete(tri_vertices);
-      face_class_ids[0].push_back(classId(face));
-      face_class_dim[0].push_back(classType(face));
-    }
-    else if (F_numEdges(face) == 4) {
-      pVertex quad_vertex;
-      pPList quad_vertices = F_vertices(face,1);
-      assert (PList_size(quad_vertices) == 4);
-      void *iter = 0;
-      while ((quad_vertex = (pVertex) PList_next(quad_vertices, &iter))) {
-        face_vertices[1].push_back(EN_id(quad_vertex));
+      vtxInfo.id[v] = classId(vtx);
+      vtxInfo.dim[v] = classType(vtx);
+      if(hasNumbering) {
+        vtxInfo.numbering[v] = getNumber(numbering,vtx);
+      }
+      ++v;
+    }
+    VIter_delete(vertices);
+    return vtxInfo;
+  }
+
+  void setVtxIds(pPList listVerts, const int vtxPerEnt, const int entIdx, HostWrite<LO>& verts) {
+    assert (PList_size(listVerts) == vtxPerEnt);
+    void *iter = 0;
+    int i = 0;
+    pVertex vtx;
+    while ((vtx = (pVertex) PList_next(listVerts, &iter))) {
+      verts[entIdx*vtxPerEnt+i] = EN_id(vtx);
+      i++;
+    }
+    PList_delete(listVerts);
+  }
+
+  EntClass readEdges(pMesh m) {
+    const int numEdges = M_numEdges(m);
+    EntClass edgeClass;
+    edgeClass.verts = HostWrite<LO>(numEdges*2);
+    edgeClass.id = HostWrite<LO>(numEdges);
+    edgeClass.dim = HostWrite<I8>(numEdges);
+    const int vtxPerEdge = 2;
+    auto edgeIdx = 0;
+    
+    EIter edges = M_edgeIter(m);
+    pEdge edge;
+    while ((edge = (pEdge) EIter_next(edges))) {
+      pPList verts = PList_new();
+      //there is no E_vertices(edge) function that returns a list
+      verts = PList_append(verts, E_vertex(edge,0));
+      verts = PList_append(verts, E_vertex(edge,1));
+      setVtxIds(verts, vtxPerEdge, edgeIdx, edgeClass.verts);
+      edgeClass.id[edgeIdx] = classId(edge);
+      edgeClass.dim[edgeIdx] = classType(edge);
+      edgeIdx++;
+    }
+    EIter_delete(edges);
+    return edgeClass;
+  }
+
+  struct MixedFaceClass {
+    EntClass tri;
+    EntClass quad;
+  };
+
+  MixedFaceClass readMixedFaces(pMesh m, GO count_tri, GO count_quad) {
+    EntClass tri;
+    tri.verts = HostWrite<LO>(count_tri*3);
+    tri.id = HostWrite<LO>(count_tri);
+    tri.dim = HostWrite<I8>(count_tri);
+    const int edgePerTri = 3;
+    const int vtxPerTri = 3;
+    int triIdx = 0;
+
+    EntClass quad;
+    quad.verts = HostWrite<LO>(count_quad*4);
+    quad.id = HostWrite<LO>(count_quad);
+    quad.dim = HostWrite<I8>(count_quad);
+    const int edgePerQuad = 4;
+    const int vtxPerQuad = 4;
+    int quadIdx = 0;
+
+    FIter faces = M_faceIter(m);
+    pFace face;
+    while ((face = (pFace) FIter_next(faces))) {
+      if (F_numEdges(face) == edgePerTri) {
+        pPList verts = F_vertices(face,1);
+        setVtxIds(verts, vtxPerTri, triIdx, tri.verts);
+        tri.id[triIdx] = classId(face);
+        tri.dim[triIdx] = classType(face);
+        triIdx++;
+      }
+      else if (F_numEdges(face) == edgePerQuad) {
+        pPList verts = F_vertices(face,1);
+        setVtxIds(verts, vtxPerQuad, quadIdx, quad.verts);
+        quad.id[quadIdx] = classId(face);
+        quad.dim[quadIdx] = classType(face);
+        quadIdx++;
+      }
+      else {
+        Omega_h_fail ("Face is neither tri nor quad \n");
       }
-      PList_delete(quad_vertices);
-      face_class_ids[1].push_back(classId(face));
-      face_class_dim[1].push_back(classType(face));
-    }
-    else {
-      Omega_h_fail ("Face is neither tri nor quad \n");
-    }
-    ent_class_ids[2].push_back(classId(face));
-    ent_class_dim[2].push_back(classType(face));
-  }
-  FIter_delete(faces);
-
-  HostWrite<LO> host_class_ids_face(numFaces);
-  HostWrite<I8> host_class_dim_face(numFaces);
-  for (int i = 0; i < numFaces; ++i) {
-    host_class_ids_face[i] = ent_class_ids[2][static_cast<std::size_t>(i)];
-    host_class_dim_face[i] = ent_class_dim[2][static_cast<std::size_t>(i)];
-  }
-  HostWrite<LO> host_class_ids_tri(count_tri);
-  HostWrite<I8> host_class_dim_tri(count_tri);
-  for (int i = 0; i < count_tri; ++i) {
-    host_class_ids_tri[i] = face_class_ids[0][static_cast<std::size_t>(i)];
-    host_class_dim_tri[i] = face_class_dim[0][static_cast<std::size_t>(i)];
-  }
-  HostWrite<LO> host_class_ids_quad(count_quad);
-  HostWrite<I8> host_class_dim_quad(count_quad);
-  for (int i = 0; i < count_quad; ++i) {
-    host_class_ids_quad[i] = face_class_ids[1][static_cast<std::size_t>(i)];
-    host_class_dim_quad[i] = face_class_dim[1][static_cast<std::size_t>(i)];
-  }
-
-  Adj edge2vert;
-  Adj vert2edge;
-  if (is_simplex || is_hypercube) {
-    edge2vert = mesh->get_adj(1, 0);
-    vert2edge = mesh->ask_up(0, 1);
-  }
-  else {
-    edge2vert = mesh->get_adj(Topo_type::edge, Topo_type::vertex);
-    vert2edge = mesh->ask_up(Topo_type::vertex, Topo_type::edge);
-  }
-  HostWrite<LO> host_tri2verts(count_tri*3);
-  for (Int i = 0; i < count_tri; ++i) {
-    for (Int j = 0; j < 3; ++j) {
-      host_tri2verts[i*3 + j] =
-          face_vertices[0][static_cast<std::size_t>(i*3 + j)];
-    }
-  }
-  auto tri2verts = Read<LO>(host_tri2verts.write());
-  Adj down;
-  if (is_simplex) {
-    down = reflect_down(tri2verts, edge2vert.ab2b, vert2edge,
-                        OMEGA_H_SIMPLEX, 2, 1);
-    mesh->set_ents(2, down);
-    mesh->add_tag<ClassId>(2, "class_id", 1,
-                           Read<ClassId>(host_class_ids_face.write()));
-    mesh->add_tag<I8>(2, "class_dim", 1,
-                      Read<I8>(host_class_dim_face.write()));
-  }
-  else if (is_hypercube) {
-    //empty since a quad/hex mesh has no triangles and to avoid dropping into
-    //the 'else'
-  }
-  else {
-    down = reflect_down(tri2verts, edge2vert.ab2b, vert2edge,
-                        Topo_type::triangle, Topo_type::edge);
-    mesh->set_ents(Topo_type::triangle, Topo_type::edge, down);
-    mesh->add_tag<ClassId>(Topo_type::triangle, "class_id", 1,
-                           Read<ClassId>(host_class_ids_tri.write()));
-    mesh->add_tag<I8>(Topo_type::triangle, "class_dim", 1,
-                      Read<I8>(host_class_dim_tri.write()));
-  }
-
-  HostWrite<LO> host_quad2verts(count_quad*4);
-  for (Int i = 0; i < count_quad; ++i) {
-    for (Int j = 0; j < 4; ++j) {
-      host_quad2verts[i*4 + j] =
-          face_vertices[1][static_cast<std::size_t>(i*4 + j)];
     }
-  }
-  auto quad2verts = Read<LO>(host_quad2verts.write());
-
-  if (is_hypercube) {
-    down = reflect_down(quad2verts, edge2vert.ab2b, vert2edge,
-                        OMEGA_H_HYPERCUBE, 2, 1);
-    mesh->set_ents(2, down);
-    mesh->add_tag<ClassId>(2, "class_id", 1,
-                           Read<ClassId>(host_class_ids_face.write()));
-    mesh->add_tag<I8>(2, "class_dim", 1,
-                      Read<I8>(host_class_dim_face.write()));
-  }
-  else if (is_simplex) {
-    //empty to avoid dropping into the 'else'
-  }
-  else {
-    down = reflect_down(quad2verts, edge2vert.ab2b, vert2edge,
-                      Topo_type::quadrilateral, Topo_type::edge);
-    mesh->set_ents(Topo_type::quadrilateral, Topo_type::edge, down);
-    mesh->add_tag<ClassId>(Topo_type::quadrilateral, "class_id", 1,
-                           Read<ClassId>(host_class_ids_quad.write()));
-    mesh->add_tag<I8>(Topo_type::quadrilateral, "class_dim", 1,
-                      Read<I8>(host_class_dim_quad.write()));
-  }
-
-  if (!(count_tet == 0 && count_hex == 0 && count_wedge == 0 && 
-        count_pyramid == 0)) {
-    rgn_vertices[0].reserve(count_tet*4);
-    rgn_vertices[1].reserve(count_hex*8);
-    rgn_vertices[2].reserve(count_wedge*6);
-    rgn_vertices[3].reserve(count_pyramid*5);
-    std::vector<int> rgn_class_ids[4];
-    std::vector<int> rgn_class_dim[4];
-    rgn_class_ids[0].reserve(count_tet);
-    rgn_class_dim[0].reserve(count_tet);
-    rgn_class_ids[1].reserve(count_hex);
-    rgn_class_dim[1].reserve(count_hex);
-    rgn_class_ids[2].reserve(count_wedge);
-    rgn_class_dim[2].reserve(count_wedge);
-    rgn_class_ids[3].reserve(count_pyramid);
-    rgn_class_dim[3].reserve(count_pyramid);
-
-    regions = M_regionIter(m);
+    FIter_delete(faces);
+
+    return MixedFaceClass{tri,quad};
+  }
+ 
+  EntClass readMonoTopoFaces(pMesh m, GO numFaces, LO vtxPerFace) {
+    EntClass ents;
+    ents.verts = HostWrite<LO>(numFaces*vtxPerFace);
+    ents.id = HostWrite<LO>(numFaces);
+    ents.dim = HostWrite<I8>(numFaces);
+    int faceIdx = 0;
+
+    FIter faces = M_faceIter(m);
+    pFace face;
+    while ((face = (pFace) FIter_next(faces))) {
+      assert(F_numEdges(face) == vtxPerFace);
+      pPList verts = F_vertices(face,1);
+      setVtxIds(verts, vtxPerFace, faceIdx, ents.verts);
+      ents.id[faceIdx] = classId(face);
+      ents.dim[faceIdx] = classType(face);
+      faceIdx++;
+    }
+    FIter_delete(faces);
+    return ents;
+  }
+
+  struct MixedRgnClass {
+    EntClass tet;
+    EntClass hex;
+    EntClass wedge;
+    EntClass pyramid;
+  };
+
+  MixedRgnClass readMixedTopoRegions(pMesh m,
+      LO numTets, LO numHexs,
+      LO numWedges, LO numPyramids) {
+
+    EntClass tet;
+    const auto vtxPerTet = 4;
+    tet.verts = HostWrite<LO>(numTets*vtxPerTet);
+    tet.id = HostWrite<LO>(numTets);
+    tet.dim = HostWrite<I8>(numTets);
+    auto tetIdx = 0;
+
+    EntClass hex;
+    const auto vtxPerHex = 8;
+    hex.verts = HostWrite<LO>(numHexs*vtxPerHex);
+    hex.id = HostWrite<LO>(numHexs);
+    hex.dim = HostWrite<I8>(numHexs);
+    auto hexIdx = 0;
+
+    EntClass wedge;
+    const auto vtxPerWedge = 6;
+    wedge.verts = HostWrite<LO>(numWedges*vtxPerWedge);
+    wedge.id = HostWrite<LO>(numWedges);
+    wedge.dim = HostWrite<I8>(numWedges);
+    auto wedgeIdx = 0;
+
+    EntClass pyramid;
+    const auto vtxPerPyramid = 5;
+    pyramid.verts = HostWrite<LO>(numPyramids*vtxPerPyramid);
+    pyramid.id = HostWrite<LO>(numPyramids);
+    pyramid.dim = HostWrite<I8>(numPyramids);
+    auto pyramidIdx = 0;
+
+    RIter regions = M_regionIter(m);
+    pRegion rgn;
     while ((rgn = (pRegion) RIter_next(regions))) {
       if (R_topoType(rgn) == Rtet) {
-        pVertex vert;
         pPList verts = R_vertices(rgn,1);
-        assert (PList_size(verts) == 4);
-        void *iter = 0;
-        while ((vert = (pVertex) PList_next(verts, &iter))) {
-          rgn_vertices[0].push_back(EN_id(vert));
-        }
-        PList_delete(verts);
-        rgn_class_ids[0].push_back(classId(rgn));
-        rgn_class_dim[0].push_back(classType(rgn));
+        setVtxIds(verts, vtxPerTet, tetIdx, tet.verts);
+        tet.id[tetIdx] = classId(rgn);
+        tet.dim[tetIdx] = classType(rgn);
+        tetIdx++;
       }
       else if (R_topoType(rgn) == Rhex) {
-        pVertex vert;
         pPList verts = R_vertices(rgn,1);
-        assert (PList_size(verts) == 8);
-        void *iter = 0;
-        while ((vert = (pVertex) PList_next(verts, &iter))) {
-          rgn_vertices[1].push_back(EN_id(vert));
-        }
-        PList_delete(verts);
-        rgn_class_ids[1].push_back(classId(rgn));
-        rgn_class_dim[1].push_back(classType(rgn));
+        setVtxIds(verts, vtxPerHex, hexIdx, hex.verts);
+        hex.id[hexIdx] = classId(rgn);
+        hex.dim[hexIdx] = classType(rgn);
+        hexIdx++;
       }
       else if (R_topoType(rgn) == Rwedge) {
-        pVertex vert;
         pPList verts = R_vertices(rgn,1);
-        assert (PList_size(verts) == 6);
-        void *iter = 0;
-        while ((vert = (pVertex) PList_next(verts, &iter))) {
-          rgn_vertices[2].push_back(EN_id(vert));
-        }
-        PList_delete(verts);
-        rgn_class_ids[2].push_back(classId(rgn));
-        rgn_class_dim[2].push_back(classType(rgn));
+        setVtxIds(verts, vtxPerWedge, wedgeIdx, wedge.verts);
+        wedge.id[wedgeIdx] = classId(rgn);
+        wedge.dim[wedgeIdx] = classType(rgn);
+        wedgeIdx++;
       }
       else if (R_topoType(rgn) == Rpyramid) {
-        pVertex vert;
         pPList verts = R_vertices(rgn,1);
-        assert (PList_size(verts) == 5);
-        void *iter = 0;
-        while ((vert = (pVertex) PList_next(verts, &iter))) {
-          rgn_vertices[3].push_back(EN_id(vert));
-        }
-        PList_delete(verts);
-        rgn_class_ids[3].push_back(classId(rgn));
-        rgn_class_dim[3].push_back(classType(rgn));
+        setVtxIds(verts, vtxPerPyramid, pyramidIdx, pyramid.verts);
+        pyramid.id[pyramidIdx] = classId(rgn);
+        pyramid.dim[pyramidIdx] = classType(rgn);
+        pyramidIdx++;
       }
       else {
         Omega_h_fail ("Region is not tet, hex, wedge, or pyramid \n");
       }
-      ent_class_ids[3].push_back(classId(rgn));
-      ent_class_dim[3].push_back(classType(rgn));
     }
     RIter_delete(regions);
+    return MixedRgnClass({tet,hex,wedge,pyramid});
+  }
 
-    HostWrite<LO> host_class_ids_rgn(numRegions);
-    HostWrite<I8> host_class_dim_rgn(numRegions);
-    for (int i = 0; i < numRegions; ++i) {
-      host_class_ids_rgn[i] = ent_class_ids[3][static_cast<std::size_t>(i)];
-      host_class_dim_rgn[i] = ent_class_dim[3][static_cast<std::size_t>(i)];
-    }
-    HostWrite<LO> host_class_ids_tet(count_tet);
-    HostWrite<I8> host_class_dim_tet(count_tet);
-    for (int i = 0; i < count_tet; ++i) {
-      host_class_ids_tet[i] = rgn_class_ids[0][static_cast<std::size_t>(i)];
-      host_class_dim_tet[i] = rgn_class_dim[0][static_cast<std::size_t>(i)];
-    }
-    HostWrite<LO> host_class_ids_hex(count_hex);
-    HostWrite<I8> host_class_dim_hex(count_hex);
-    for (int i = 0; i < count_hex; ++i) {
-      host_class_ids_hex[i] = rgn_class_ids[1][static_cast<std::size_t>(i)];
-      host_class_dim_hex[i] = rgn_class_dim[1][static_cast<std::size_t>(i)];
-    }
-    HostWrite<LO> host_class_ids_wedge(count_wedge);
-    HostWrite<I8> host_class_dim_wedge(count_wedge);
-    for (int i = 0; i < count_wedge; ++i) {
-      host_class_ids_wedge[i] = rgn_class_ids[2][static_cast<std::size_t>(i)];
-      host_class_dim_wedge[i] = rgn_class_dim[2][static_cast<std::size_t>(i)];
-    }
-    HostWrite<LO> host_class_ids_pyramid(count_pyramid);
-    HostWrite<I8> host_class_dim_pyramid(count_pyramid);
-    for (int i = 0; i < count_pyramid; ++i) {
-      host_class_ids_pyramid[i] = rgn_class_ids[3][static_cast<std::size_t>(i)];
-      host_class_dim_pyramid[i] = rgn_class_dim[3][static_cast<std::size_t>(i)];
+  EntClass readMonoTopoRegions(pMesh m, GO numRgn, LO vtxPerRgn) {
+    assert(vtxPerRgn == 4 || vtxPerRgn == 8);
+    if(vtxPerRgn != 4 && vtxPerRgn != 8) {
+      Omega_h_fail("Mono topology 3d mesh must be all tets or all hex\n");
     }
 
-    Adj tri2vert;
-    Adj vert2tri;
-    Adj quad2vert;
-    Adj vert2quad;
-    if (is_simplex) {
-      tri2vert = mesh->ask_down(2, 0);
-      vert2tri = mesh->ask_up(0, 2);
-    }
-    else if (is_hypercube) {
-      quad2vert = mesh->ask_down(2, 0);
-      vert2quad = mesh->ask_up(0, 2);
-    }
-    else {
-      tri2vert = mesh->ask_down(Topo_type::triangle, Topo_type::vertex);
-      vert2tri = mesh->ask_up(Topo_type::vertex, Topo_type::triangle);
-      quad2vert = mesh->ask_down(Topo_type::quadrilateral, Topo_type::vertex);
-      vert2quad = mesh->ask_up(Topo_type::vertex, Topo_type::quadrilateral);
-    }
+    EntClass rgnClass;
+    rgnClass.verts = HostWrite<LO>(numRgn*vtxPerRgn);
+    rgnClass.id = HostWrite<LO>(numRgn);
+    rgnClass.dim = HostWrite<I8>(numRgn);
 
-    HostWrite<LO> host_tet2verts(count_tet*4);
-    for (Int i = 0; i < count_tet; ++i) {
-      for (Int j = 0; j < 4; ++j) {
-        host_tet2verts[i*4 + j] =
-          rgn_vertices[0][static_cast<std::size_t>(i*4 + j)];
-      }
-    }
-    auto tet2verts = Read<LO>(host_tet2verts.write());
-    if (is_simplex) {
-      down = reflect_down(tet2verts, tri2vert.ab2b, vert2tri,
-          OMEGA_H_SIMPLEX, 3, 2);
-      mesh->set_ents(3, down);
-      mesh->add_tag<ClassId>(3, "class_id", 1,
-          Read<ClassId>(host_class_ids_rgn.write()));
-      mesh->add_tag<I8>(3, "class_dim", 1,
-          Read<I8>(host_class_dim_rgn.write()));
-    }
-    else if (is_hypercube) {
-      //empty to avoid dropping into the 'else'
-    }
-    else {
-      down = reflect_down(tet2verts, tri2vert.ab2b, vert2tri,
-          Topo_type::tetrahedron, Topo_type::triangle);
-      mesh->set_ents(Topo_type::tetrahedron, Topo_type::triangle, down);
-      mesh->add_tag<ClassId>(Topo_type::tetrahedron, "class_id", 1,
-          Read<ClassId>(host_class_ids_tet.write()));
-      mesh->add_tag<I8>(Topo_type::tetrahedron, "class_dim", 1,
-          Read<I8>(host_class_dim_tet.write()));
-    }
+    const auto rgnType = ( vtxPerRgn == 4 ) ? Rtet : Rhex;
 
-    HostWrite<LO> host_hex2verts(count_hex*8);
-    for (Int i = 0; i < count_hex; ++i) {
-      for (Int j = 0; j < 8; ++j) {
-        host_hex2verts[i*8 + j] =
-          rgn_vertices[1][static_cast<std::size_t>(i*8 + j)];
-      }
-    }
-    auto hex2verts = Read<LO>(host_hex2verts.write());
-
-    if (is_hypercube) {
-      down = reflect_down(hex2verts, quad2vert.ab2b, vert2quad,
-          OMEGA_H_HYPERCUBE, 3, 2);
-      mesh->set_ents(3, down);
-      mesh->add_tag<ClassId>(3, "class_id", 1,
-          Read<ClassId>(host_class_ids_rgn.write()));
-      mesh->add_tag<I8>(3, "class_dim", 1,
-          Read<I8>(host_class_dim_rgn.write()));
-    }
-    else if (is_simplex) {
-      //empty to avoid dropping into the 'else'
-    }
-    else {
-      down = reflect_down(hex2verts, quad2vert.ab2b, vert2quad,
-          Topo_type::hexahedron, Topo_type::quadrilateral);
-      mesh->set_ents(Topo_type::hexahedron, Topo_type::quadrilateral, down);
-      mesh->add_tag<ClassId>(Topo_type::hexahedron, "class_id", 1,
-          Read<ClassId>(host_class_ids_hex.write()));
-      mesh->add_tag<I8>(Topo_type::hexahedron, "class_dim", 1,
-          Read<I8>(host_class_dim_hex.write()));
+    RIter regions = M_regionIter(m);
+    pRegion rgn;
+    int rgnIdx = 0;
+    while ((rgn = (pRegion) RIter_next(regions))) {
+      assert(R_topoType(rgn) == rgnType);
+      pPList verts = R_vertices(rgn,1);
+      setVtxIds(verts, vtxPerRgn, rgnIdx, rgnClass.verts);
+      rgnClass.id[rgnIdx] = classId(rgn);
+      rgnClass.dim[rgnIdx] = classType(rgn);
+      rgnIdx++;
     }
+    RIter_delete(regions);
 
-    HostWrite<LO> host_wedge2verts(count_wedge*6);
-    for (Int i = 0; i < count_wedge; ++i) {
-      for (Int j = 0; j < 6; ++j) {
-        host_wedge2verts[i*6 + j] =
-          rgn_vertices[2][static_cast<std::size_t>(i*6 + j)];
-      }
-    }
-    auto wedge2verts = Read<LO>(host_wedge2verts.write());
-    down = reflect_down(wedge2verts, quad2vert.ab2b, vert2quad,
-        Topo_type::wedge, Topo_type::quadrilateral);
-    if ((!is_simplex) && (!is_hypercube)) {
-      mesh->set_ents(Topo_type::wedge, Topo_type::quadrilateral, down);
-    }
+    return rgnClass;
+  }
 
-    down = reflect_down(wedge2verts, tri2vert.ab2b, vert2tri,
-        Topo_type::wedge, Topo_type::triangle);
-    if ((!is_simplex) && (!is_hypercube)) {
-      mesh->set_ents(Topo_type::wedge, Topo_type::triangle, down);
-      mesh->add_tag<ClassId>(Topo_type::wedge, "class_id", 1,
-          Read<ClassId>(host_class_ids_wedge.write()));
-      mesh->add_tag<I8>(Topo_type::wedge, "class_dim", 1,
-          Read<I8>(host_class_dim_wedge.write()));
-    }
 
-    HostWrite<LO> host_pyramid2verts(count_pyramid*5);
-    for (Int i = 0; i < count_pyramid; ++i) {
-      for (Int j = 0; j < 5; ++j) {
-        host_pyramid2verts[i*5 + j] =
-          rgn_vertices[3][static_cast<std::size_t>(i*5 + j)];
-      }
-    }
-    auto pyramid2verts = Read<LO>(host_pyramid2verts.write());
-    down = reflect_down(pyramid2verts, tri2vert.ab2b, vert2tri,
-        Topo_type::pyramid, Topo_type::triangle);
-    if ((!is_simplex) && (!is_hypercube)) {
-      mesh->set_ents(Topo_type::pyramid, Topo_type::triangle, down);
+  private:
+  SimMeshEntInfo();
+  int getMaxDim(std::array<int,4> numEnts) {
+    int max_dim;
+    if (numEnts[3]) {
+      max_dim = 3;
+    } else if (numEnts[2]) {
+      max_dim = 2;
+    } else if (numEnts[1]) {
+      max_dim = 1;
+    } else {
+      Omega_h_fail("There were no Elements of dimension higher than zero!\n");
     }
+    return max_dim;
+  }
+}; //end SimMeshEntInfo
 
-    down = reflect_down(pyramid2verts, quad2vert.ab2b, vert2quad,
-        Topo_type::pyramid, Topo_type::quadrilateral);
-    if ((!is_simplex) && (!is_hypercube)) {
-      mesh->set_ents(Topo_type::pyramid, Topo_type::quadrilateral, down);
-      mesh->add_tag<ClassId>(Topo_type::pyramid, "class_id", 1,
-          Read<ClassId>(host_class_ids_pyramid.write()));
-      mesh->add_tag<I8>(Topo_type::pyramid, "class_dim", 1,
-          Read<I8>(host_class_dim_pyramid.write()));
-    }
+void readMixed_internal(pMesh m, MixedMesh* mesh, SimMeshInfo info) {
+  assert(!info.is_simplex && !info.is_hypercube);
+  if(info.is_simplex || info.is_hypercube) {
+      Omega_h_fail("Attempting to use the mixed mesh reader for a mono topology"
+                   "mesh (family = simplex|hypercube)!\n");
   }
 
-  return;
+  const int numVtx = M_numVertices(m);
+  const int numEdges = M_numEdges(m);
+  const int numFaces = M_numFaces(m);
+  const int numRegions = M_numRegions(m);
+
+  const bool hasNumbering = false;
+
+  SimMeshEntInfo simEnts({{numVtx,numEdges,numFaces,numRegions}}, hasNumbering);
+  mesh->set_dim(simEnts.maxDim);
+
+  // process verts
+  auto vtxInfo = simEnts.readVerts(m,nullptr);
+  mesh->set_verts_type(numVtx);
+  mesh->add_coords_mix(vtxInfo.coords.write());
+  mesh->add_tag<ClassId>(Topo_type::vertex, "class_id", 1,
+                Read<ClassId>(vtxInfo.id.write()));
+  mesh->add_tag<I8>(Topo_type::vertex, "class_dim", 1,
+                    Read<I8>(vtxInfo.dim.write()));
+
+  // process edges
+  auto edges = simEnts.readEdges(m);
+  auto ev2v = Read<LO>(edges.verts.write());
+  mesh->set_ents(Topo_type::edge, Topo_type::vertex, Adj(ev2v));
+  mesh->add_tag<ClassId>(Topo_type::edge, "class_id", 1,
+                         Read<ClassId>(edges.id.write()));
+  mesh->add_tag<I8>(Topo_type::edge, "class_dim", 1,
+                    Read<I8>(edges.dim.write()));
+
+  //process faces
+  auto mixedFaceClass = simEnts.readMixedFaces(m, info.count_tri, info.count_quad);
+  auto edge2vert = mesh->get_adj(Topo_type::edge, Topo_type::vertex);
+  auto vert2edge = mesh->ask_up(Topo_type::vertex, Topo_type::edge);
+
+  //// tris
+  auto tri2verts = Read<LO>(mixedFaceClass.tri.verts.write());
+  auto down = reflect_down(tri2verts, edge2vert.ab2b, vert2edge,
+      Topo_type::triangle, Topo_type::edge);
+  mesh->set_ents(Topo_type::triangle, Topo_type::edge, down);
+  mesh->add_tag<ClassId>(Topo_type::triangle, "class_id", 1,
+      Read<ClassId>(mixedFaceClass.tri.id.write()));
+  mesh->add_tag<I8>(Topo_type::triangle, "class_dim", 1,
+      Read<I8>(mixedFaceClass.tri.dim.write()));
+
+  //// quads
+  auto quad2verts = Read<LO>(mixedFaceClass.quad.verts.write());
+  down = reflect_down(quad2verts, edge2vert.ab2b, vert2edge,
+      Topo_type::quadrilateral, Topo_type::edge);
+  mesh->set_ents(Topo_type::quadrilateral, Topo_type::edge, down);
+  mesh->add_tag<ClassId>(Topo_type::quadrilateral, "class_id", 1,
+      Read<ClassId>(mixedFaceClass.quad.id.write()));
+  mesh->add_tag<I8>(Topo_type::quadrilateral, "class_dim", 1,
+      Read<I8>(mixedFaceClass.quad.dim.write()));
+
+  //process regions
+  if(simEnts.maxDim == 2)
+    return; //there are no regions
+
+  auto mixedRgnClass = simEnts.readMixedTopoRegions(m,
+      info.count_tet, info.count_hex,
+      info.count_wedge, info.count_pyramid);
+  auto tri2vert = mesh->ask_down(Topo_type::triangle, Topo_type::vertex);
+  auto vert2tri = mesh->ask_up(Topo_type::vertex, Topo_type::triangle);
+  auto quad2vert = mesh->ask_down(Topo_type::quadrilateral, Topo_type::vertex);
+  auto vert2quad = mesh->ask_up(Topo_type::vertex, Topo_type::quadrilateral);
+
+  /// tets
+  auto tet2verts = Read<LO>(mixedRgnClass.tet.verts.write());
+  down = reflect_down(tet2verts, tri2vert.ab2b, vert2tri,
+      Topo_type::tetrahedron, Topo_type::triangle);
+  mesh->set_ents(Topo_type::tetrahedron, Topo_type::triangle, down);
+  mesh->add_tag<ClassId>(Topo_type::tetrahedron, "class_id", 1,
+      Read<ClassId>(mixedRgnClass.tet.id.write()));
+  mesh->add_tag<I8>(Topo_type::tetrahedron, "class_dim", 1,
+      Read<I8>(mixedRgnClass.tet.dim.write()));
+
+  /// hexs
+  auto hex2verts = Read<LO>(mixedRgnClass.hex.verts.write());
+  down = reflect_down(hex2verts, quad2vert.ab2b, vert2quad,
+      Topo_type::hexahedron, Topo_type::quadrilateral);
+  mesh->set_ents(Topo_type::hexahedron, Topo_type::quadrilateral, down);
+  mesh->add_tag<ClassId>(Topo_type::hexahedron, "class_id", 1,
+      Read<ClassId>(mixedRgnClass.hex.id.write()));
+  mesh->add_tag<I8>(Topo_type::hexahedron, "class_dim", 1,
+      Read<I8>(mixedRgnClass.hex.dim.write()));
+
+  /// wedges
+  auto wedge2verts = Read<LO>(mixedRgnClass.wedge.verts.write());
+  down = reflect_down(wedge2verts, quad2vert.ab2b, vert2quad,
+      Topo_type::wedge, Topo_type::quadrilateral);
+  mesh->set_ents(Topo_type::wedge, Topo_type::quadrilateral, down);
+  down = reflect_down(wedge2verts, tri2vert.ab2b, vert2tri,
+      Topo_type::wedge, Topo_type::triangle);
+  mesh->set_ents(Topo_type::wedge, Topo_type::triangle, down);
+  mesh->add_tag<ClassId>(Topo_type::wedge, "class_id", 1,
+      Read<ClassId>(mixedRgnClass.wedge.id.write()));
+  mesh->add_tag<I8>(Topo_type::wedge, "class_dim", 1,
+      Read<I8>(mixedRgnClass.wedge.dim.write()));
+
+  /// pyramids
+  auto pyramid2verts = Read<LO>(mixedRgnClass.pyramid.verts.write());
+  down = reflect_down(pyramid2verts, tri2vert.ab2b, vert2tri,
+      Topo_type::pyramid, Topo_type::triangle);
+  mesh->set_ents(Topo_type::pyramid, Topo_type::triangle, down);
+  down = reflect_down(pyramid2verts, quad2vert.ab2b, vert2quad,
+      Topo_type::pyramid, Topo_type::quadrilateral);
+  mesh->set_ents(Topo_type::pyramid, Topo_type::quadrilateral, down);
+  mesh->add_tag<ClassId>(Topo_type::pyramid, "class_id", 1,
+      Read<ClassId>(mixedRgnClass.pyramid.id.write()));
+  mesh->add_tag<I8>(Topo_type::pyramid, "class_dim", 1,
+      Read<I8>(mixedRgnClass.pyramid.dim.write()));
+}
+
+void read_internal(pMesh m, Mesh* mesh, pMeshNex numbering, SimMeshInfo info) {
+  assert(info.is_simplex || info.is_hypercube);
+  if(!info.is_simplex && !info.is_hypercube) {
+      Omega_h_fail("Attempting to use the mono topology reader for a mixed"
+                   "mesh (family = mixed)!\n");
+  }
+
+  const int numVtx = M_numVertices(m);
+  const int numEdges = M_numEdges(m);
+  const int numFaces = M_numFaces(m);
+  const int numRegions = M_numRegions(m);
+
+  const bool hasNumbering = (numbering != NULL);
+
+  SimMeshEntInfo simEnts({{numVtx,numEdges,numFaces,numRegions}}, hasNumbering);
+  mesh->set_dim(simEnts.maxDim);
+  if (info.is_simplex) {
+    mesh->set_family(OMEGA_H_SIMPLEX);
+  } else if (info.is_hypercube){
+    mesh->set_family(OMEGA_H_HYPERCUBE);
+  }
+
+  //process verts
+  auto vtxInfo = simEnts.readVerts(m,numbering);
+  mesh->set_verts(numVtx);
+  mesh->add_coords(vtxInfo.coords.write());
+  mesh->add_tag<ClassId>(0, "class_id", 1,
+      Read<ClassId>(vtxInfo.id.write()));
+  mesh->add_tag<I8>(0, "class_dim", 1,
+      Read<I8>(vtxInfo.dim.write()));
+  if(hasNumbering) {
+    mesh->add_tag<LO>(0, "simNumbering", 1,
+        Read<LO>(vtxInfo.numbering.write()));
+  }
+
+  //process edges
+  auto edges = simEnts.readEdges(m);
+  auto ev2v = Read<LO>(edges.verts.write());
+  mesh->set_ents(1, Adj(ev2v));
+  mesh->add_tag<ClassId>(1, "class_id", 1,
+                Read<ClassId>(edges.id.write()));
+  mesh->add_tag<I8>(1, "class_dim", 1,
+                    Read<I8>(edges.dim.write()));
+
+  //process faces
+  if(info.is_simplex) {
+    const auto vtxPerTri = 3;
+    auto entClass = simEnts.readMonoTopoFaces(m, info.count_tri, vtxPerTri);
+    auto edge2vert = mesh->get_adj(1, 0);
+    auto vert2edge = mesh->ask_up(0, 1);
+    auto tri2verts = Read<LO>(entClass.verts.write());
+    auto down = reflect_down(tri2verts, edge2vert.ab2b, vert2edge,
+                             OMEGA_H_SIMPLEX, 2, 1);
+    mesh->set_ents(2, down);
+    mesh->add_tag<ClassId>(2, "class_id", 1,
+                           Read<ClassId>(entClass.id.write()));
+    mesh->add_tag<I8>(2, "class_dim", 1,
+                      Read<I8>(entClass.dim.write()));
+  } else { // hypercube
+    const auto vtxPerQuad = 4;
+    auto entClass = simEnts.readMonoTopoFaces(m, info.count_quad, vtxPerQuad);
+    auto edge2vert = mesh->get_adj(1, 0);
+    auto vert2edge = mesh->ask_up(0, 1);
+    auto quad2verts = Read<LO>(entClass.verts.write());
+    auto down = reflect_down(quad2verts, edge2vert.ab2b, vert2edge,
+                        OMEGA_H_HYPERCUBE, 2, 1);
+    mesh->set_ents(2, down);
+    mesh->add_tag<ClassId>(2, "class_id", 1,
+                           Read<ClassId>(entClass.id.write()));
+    mesh->add_tag<I8>(2, "class_dim", 1,
+                      Read<I8>(entClass.dim.write()));
+  }
+
+  //process regions
+  if(simEnts.maxDim == 2)
+    return; //there are no regions
+
+  if(info.is_simplex) {
+    const auto vtxPerTet = 4;
+    auto entClass = simEnts.readMonoTopoRegions(m, info.count_tet, vtxPerTet);
+    auto tri2vert = mesh->ask_down(2, 0);
+    auto vert2tri = mesh->ask_up(0, 2);
+    auto tet2verts = Read<LO>(entClass.verts.write());
+    auto down = reflect_down(tet2verts, tri2vert.ab2b, vert2tri,
+        OMEGA_H_SIMPLEX, 3, 2);
+    mesh->set_ents(3, down);
+    mesh->template add_tag<ClassId>(3, "class_id", 1,
+        Read<ClassId>(entClass.id.write()));
+    mesh->template add_tag<I8>(3, "class_dim", 1,
+        Read<I8>(entClass.dim.write()));
+  } else { //hypercube
+    const auto vtxPerHex = 8;
+    auto entClass = simEnts.readMonoTopoRegions(m, info.count_hex, vtxPerHex);
+    auto quad2vert = mesh->ask_down(2, 0);
+    auto vert2quad = mesh->ask_up(0, 2);
+    auto hex2verts = Read<LO>(entClass.verts.write());
+    auto down = reflect_down(hex2verts, quad2vert.ab2b, vert2quad,
+        OMEGA_H_HYPERCUBE, 3, 2);
+    mesh->set_ents(3, down);
+    mesh->template add_tag<ClassId>(3, "class_id", 1,
+        Read<ClassId>(entClass.id.write()));
+    mesh->template add_tag<I8>(3, "class_dim", 1,
+        Read<I8>(entClass.dim.write()));
+  }
+}
+ 
+MixedMesh readMixedImpl(filesystem::path const& mesh_fname,
+    filesystem::path const& mdl_fname,
+    CommPtr comm) {
+  SimModel_start();
+  Sim_readLicenseFile(NULL);
+  SimDiscrete_start(0);
+  pNativeModel nm = NULL;
+  pProgress p = NULL;
+  pGModel g = GM_load(mdl_fname.c_str(), nm, p);
+  pMesh m = M_load(mesh_fname.c_str(), g, p);
+  auto simMeshInfo = getSimMeshInfo(m);
+  auto mesh = MixedMesh(comm->library());
+  mesh.set_comm(comm);
+  meshsim::readMixed_internal(m, &mesh, simMeshInfo);
+  M_release(m);
+  GM_release(g);
+  SimDiscrete_stop(0);
+  SimModel_stop();
+  return mesh;
 }
 
 Mesh readImpl(filesystem::path const& mesh_fname, filesystem::path const& mdl_fname,
@@ -645,12 +663,13 @@ Mesh readImpl(filesystem::path const& mesh_fname, filesystem::path const& mdl_fn
   pProgress p = NULL;
   pGModel g = GM_load(mdl_fname.c_str(), nm, p);
   pMesh m = M_load(mesh_fname.c_str(), g, p);
+  auto simMeshInfo = getSimMeshInfo(m);
   const bool hasNumbering = (numbering_fname.native() != std::string(""));
   pMeshNex numbering = hasNumbering ? MeshNex_load(numbering_fname.c_str(), m) : nullptr;
   auto mesh = Mesh(comm->library());
   mesh.set_comm(comm);
   mesh.set_parting(OMEGA_H_ELEM_BASED);
-  meshsim::read_internal(m, &mesh, numbering);
+  meshsim::read_internal(m, &mesh, numbering, simMeshInfo);
   if(hasNumbering) MeshNex_delete(numbering);
   M_release(m);
   GM_release(g);
@@ -659,6 +678,23 @@ Mesh readImpl(filesystem::path const& mesh_fname, filesystem::path const& mdl_fn
   return mesh;
 }
 
+bool isMixed(filesystem::path const& mesh_fname, filesystem::path const& mdl_fname) {
+  SimModel_start();
+  Sim_readLicenseFile(NULL);
+  SimDiscrete_start(0);
+  pNativeModel nm = NULL;
+  pProgress p = NULL;
+  pGModel g = GM_load(mdl_fname.c_str(), nm, p);
+  pMesh m = M_load(mesh_fname.c_str(), g, p);
+  auto simMeshInfo = getSimMeshInfo(m);
+  M_release(m);
+  GM_release(g);
+  SimDiscrete_stop(0);
+  SimModel_stop();
+  bool isMixed = (!simMeshInfo.is_simplex && !simMeshInfo.is_hypercube);
+  return isMixed;
+}
+
 Mesh read(filesystem::path const& mesh_fname, filesystem::path const& mdl_fname,
     filesystem::path const& numbering_fname, CommPtr comm) {
   return readImpl(mesh_fname, mdl_fname, numbering_fname, comm);
@@ -669,6 +705,11 @@ Mesh read(filesystem::path const& mesh_fname, filesystem::path const& mdl_fname,
   return readImpl(mesh_fname, mdl_fname, std::string(""), comm);
 }
 
+MixedMesh readMixed(filesystem::path const& mesh_fname, filesystem::path const& mdl_fname,
+    CommPtr comm) {
+  return readMixedImpl(mesh_fname, mdl_fname, comm);
+}
+
 }  // namespace meshsim
 
 }  // end namespace Omega_h
diff --git a/src/Omega_h_mixedMesh.cpp b/src/Omega_h_mixedMesh.cpp
new file mode 100644
index 000000000..c72278820
--- /dev/null
+++ b/src/Omega_h_mixedMesh.cpp
@@ -0,0 +1,389 @@
+#include "Omega_h_mixedMesh.hpp"
+
+#include <algorithm>
+#include <cctype>
+#include <iostream>
+
+#include "Omega_h_array_ops.hpp"
+#include "Omega_h_bcast.hpp"
+#include "Omega_h_compare.hpp"
+#include "Omega_h_element.hpp"
+#include "Omega_h_for.hpp"
+#include "Omega_h_ghost.hpp"
+#include "Omega_h_inertia.hpp"
+#include "Omega_h_map.hpp"
+#include "Omega_h_mark.hpp"
+#include "Omega_h_migrate.hpp"
+#include "Omega_h_quality.hpp"
+#include "Omega_h_shape.hpp"
+#include "Omega_h_timer.hpp"
+#include "Omega_h_int_scan.hpp"
+#include "Omega_h_atomics.hpp"
+#include "Omega_h_reduce.hpp"
+#include "Omega_h_print.hpp"
+#include "Omega_h_dbg.hpp"
+
+namespace Omega_h {
+
+MixedMesh::MixedMesh() {
+  family_ = OMEGA_H_MIXED;
+  dim_ = -1;
+  for (Int i = 0; i <= 7; ++i) nents_type_[i] = -1;
+  library_ = nullptr;
+}
+
+MixedMesh::MixedMesh(Library* library_in) : MixedMesh() {
+  set_library(library_in);
+};
+
+void MixedMesh::set_library(Library* library_in) {
+  OMEGA_H_CHECK(library_in != nullptr);
+  library_ = library_in;
+}
+
+void MixedMesh::set_comm(CommPtr const& new_comm) {
+  OMEGA_H_CHECK(new_comm->size() == 1);
+  comm_ = new_comm;
+}
+
+void MixedMesh::set_dim(Int dim_in) {
+  OMEGA_H_CHECK(dim_ == -1);
+  OMEGA_H_CHECK(dim_in >= 1);
+  OMEGA_H_CHECK(dim_in <= 3);
+  dim_ = dim_in;
+}
+
+void MixedMesh::set_verts_type(LO nverts_in) { nents_type_[int(Topo_type::vertex)] = nverts_in; }
+
+void MixedMesh::set_ents(Topo_type high_type, Topo_type low_type, Adj h2l) {
+  OMEGA_H_TIME_FUNCTION;
+  check_type(high_type);
+  check_type(low_type);
+  if (int(high_type) < 6) {
+    OMEGA_H_CHECK(!has_ents(high_type));
+  }
+  auto deg = element_degree(high_type, low_type);
+  nents_type_[int(high_type)] = divide_no_remainder(h2l.ab2b.size(), deg);
+  add_adj(high_type, low_type, h2l);
+}
+
+LO MixedMesh::nents(Topo_type ent_type) const {
+  check_type2(ent_type);
+  return nents_type_[int(ent_type)];
+}
+
+Int MixedMesh::ent_dim(Topo_type ent_type) const {
+  Int ent_dim;
+  if (int(ent_type) == 0) {
+    ent_dim = 0;
+  }
+  else if (int(ent_type) == 1) {
+    ent_dim = 1;
+  }
+  else if ((int(ent_type) > 1) && (int(ent_type) < 4)) {
+    ent_dim = 2;
+  }
+  else {
+    ent_dim = 3;
+  }
+  return ent_dim;
+}
+
+LO MixedMesh::npyrams() const { return nents(Topo_type::pyramid); }
+
+LO MixedMesh::nwedges() const { return nents(Topo_type::wedge); }
+
+LO MixedMesh::nhexs() const { return nents(Topo_type::hexahedron); }
+
+LO MixedMesh::ntets() const { return nents(Topo_type::tetrahedron); }
+
+LO MixedMesh::nquads() const { return nents(Topo_type::quadrilateral); }
+
+LO MixedMesh::ntris() const { return nents(Topo_type::triangle); }
+
+LO MixedMesh::nedges_mix() const { return nents(Topo_type::edge); }
+
+LO MixedMesh::nverts_mix() const { return nents(Topo_type::vertex); }
+
+LO MixedMesh::nregions_mix() const {
+  return (nents(Topo_type::tetrahedron) +
+          nents(Topo_type::hexahedron) +
+          nents(Topo_type::wedge) +
+          nents(Topo_type::pyramid));
+}
+
+LO MixedMesh::nfaces_mix() const {
+  return (nents(Topo_type::triangle) +
+          nents(Topo_type::quadrilateral));
+}
+
+template <typename T>
+void MixedMesh::add_tag(Topo_type ent_type, std::string const& name, Int ncomps) {
+  if (has_tag(ent_type, name)) remove_tag(ent_type, name);
+  check_type2(ent_type);
+  check_tag_name(name);
+  OMEGA_H_CHECK(ncomps >= 0);
+  OMEGA_H_CHECK(ncomps <= Int(INT8_MAX));
+  OMEGA_H_CHECK(tags_type_[int(ent_type)].size() < size_t(INT8_MAX));
+  TagPtr ptr(new Tag<T>(name, ncomps));
+  tags_type_[int(ent_type)].push_back(std::move(ptr));
+}
+
+template <typename T>
+void MixedMesh::add_tag(Topo_type ent_type, std::string const& name, Int ncomps,
+    Read<T> array, bool internal) {
+  check_type2(ent_type);
+  auto it = tag_iter(ent_type, name);
+  auto had_tag = (it != tags_type_[int(ent_type)].end());
+  auto ptr = std::make_shared<Tag<T>>(name, ncomps);
+  ptr->set_array(array);
+  if (had_tag) {
+    OMEGA_H_CHECK(ncomps == ptr->ncomps());
+    *it = std::move(ptr);
+  } else {
+    check_tag_name(name);
+    OMEGA_H_CHECK(ncomps >= 0);
+    OMEGA_H_CHECK(ncomps <= Int(INT8_MAX));
+    OMEGA_H_CHECK(tags_type_[int(ent_type)].size() < size_t(INT8_MAX));
+    tags_type_[int(ent_type)].push_back(std::move(ptr));
+  }
+  OMEGA_H_CHECK(array.size() == nents_type_[int(ent_type)] * ncomps);
+  if (!internal) react_to_set_tag(ent_type, name);
+}
+
+template <typename T>
+void MixedMesh::set_tag(
+    Topo_type ent_type, std::string const& name, Read<T> array, bool internal) {
+  if (!has_tag(ent_type, name)) {
+    Omega_h_fail("set_tag(%s, %s): tag doesn't exist (use add_tag first)\n",
+      dimensional_plural_name(ent_type), name.c_str());
+  }
+  auto it = tag_iter(ent_type,name);
+  this->add_tag(ent_type, name, (*it)->ncomps(), array, internal);
+}
+
+void MixedMesh::react_to_set_tag(Topo_type ent_type, std::string const& name) {
+  bool is_coordinates = (name == "coordinates");
+  if ((int(ent_type) == 0) && (is_coordinates || (name == "metric"))) {
+    remove_tag(Topo_type::edge, "length");
+
+    remove_tag(Topo_type::pyramid, "quality");
+    remove_tag(Topo_type::wedge, "quality");
+    remove_tag(Topo_type::hexahedron, "quality");
+    remove_tag(Topo_type::tetrahedron, "quality");
+    remove_tag(Topo_type::quadrilateral, "quality");
+    remove_tag(Topo_type::triangle, "quality");
+  }
+  if ((int(ent_type) == 0) && is_coordinates) {
+    remove_tag(Topo_type::pyramid, "size");
+    remove_tag(Topo_type::wedge, "size");
+    remove_tag(Topo_type::hexahedron, "size");
+    remove_tag(Topo_type::tetrahedron, "size");
+    remove_tag(Topo_type::quadrilateral, "size");
+    remove_tag(Topo_type::triangle, "size");
+  }
+}
+
+TagBase const* MixedMesh::get_tagbase(Topo_type ent_type, std::string const& name) const {
+  check_type2(ent_type);
+  auto it = tag_iter(ent_type, name);
+  if (it == tags_type_[int(ent_type)].end()) {
+    Omega_h_fail("get_tagbase(%s, %s): doesn't exist\n",
+        dimensional_plural_name(ent_type), name.c_str());
+  }
+  return it->get();
+}
+
+template <typename T>
+Tag<T> const* MixedMesh::get_tag(Topo_type ent_type, std::string const& name) const {
+  return as<T>(get_tagbase(ent_type, name));
+}
+
+template <typename T>
+Read<T> MixedMesh::get_array(Topo_type ent_type, std::string const& name) const {
+  return get_tag<T>(ent_type, name)->array();
+}
+
+void MixedMesh::remove_tag(Topo_type ent_type, std::string const& name) {
+  if (!has_tag(ent_type, name)) return;
+  check_type2(ent_type);
+  OMEGA_H_CHECK(has_tag(ent_type, name));
+  tags_type_[int(ent_type)].erase(tag_iter(ent_type, name));
+}
+
+bool MixedMesh::has_tag(Topo_type ent_type, std::string const& name) const {
+  check_type(ent_type);
+  if (!has_ents(ent_type)) return false;
+  return tag_iter(ent_type, name) != tags_type_[int(ent_type)].end();
+}
+
+Int MixedMesh::ntags(Topo_type ent_type) const {
+  check_type2(ent_type);
+  return static_cast<Int>(tags_type_[int(ent_type)].size());
+}
+
+TagBase const* MixedMesh::get_tag(Topo_type ent_type, Int i) const {
+  check_type2(ent_type);
+  OMEGA_H_CHECK(0 <= i);
+  OMEGA_H_CHECK(i <= ntags(ent_type));
+  return tags_type_[int(ent_type)][static_cast<std::size_t>(i)].get();
+}
+
+bool MixedMesh::has_ents(Topo_type ent_type) const {
+  check_type(ent_type);
+  return nents_type_[int(ent_type)] >= 0;
+}
+
+bool MixedMesh::has_adj(Topo_type from_type, Topo_type to_type) const {
+  check_type(from_type);
+  check_type(to_type);
+  return bool(adjs_type_[int(from_type)][int(to_type)]);
+}
+
+Adj MixedMesh::get_adj(Topo_type from_type, Topo_type to_type) const {
+  check_type2(from_type);
+  check_type2(to_type);
+  OMEGA_H_CHECK(has_adj(from_type, to_type));
+  return *(adjs_type_[int(from_type)][int(to_type)]);
+}
+
+Adj MixedMesh::ask_down(Topo_type from_type, Topo_type to_type) {
+  OMEGA_H_CHECK(int(to_type) < int(from_type));
+  return ask_adj(from_type, to_type);
+}
+
+LOs MixedMesh::ask_verts_of(Topo_type ent_type) { return ask_adj(ent_type, Topo_type::vertex).ab2b; }
+
+Adj MixedMesh::ask_up(Topo_type from_type, Topo_type to_type) {
+  OMEGA_H_CHECK(int(from_type) < int(to_type));
+  return ask_adj(from_type, to_type);
+}
+
+MixedMesh::TagIter MixedMesh::tag_iter(Topo_type ent_type, std::string const& name) {
+  return std::find_if(tags_type_[int(ent_type)].begin(), tags_type_[int(ent_type)].end(),
+      [&](TagPtr const& a) { return a->name() == name; });
+}
+
+MixedMesh::TagCIter MixedMesh::tag_iter(Topo_type ent_type, std::string const& name) const {
+  return std::find_if(tags_type_[int(ent_type)].begin(), tags_type_[int(ent_type)].end(),
+      [&](TagPtr const& a) { return a->name() == name; });
+}
+
+void MixedMesh::check_type(Topo_type ent_type) const {
+  OMEGA_H_CHECK(Topo_type::vertex <= ent_type);
+  OMEGA_H_CHECK(ent_type <= Topo_type::pyramid);
+}
+
+void MixedMesh::check_type2(Topo_type ent_type) const {
+  check_type(ent_type);
+  OMEGA_H_CHECK(has_ents(ent_type));
+}
+
+void MixedMesh::add_adj(Topo_type from_type, Topo_type to_type, Adj adj) {
+  check_type2(from_type);
+  check_type(to_type);
+  OMEGA_H_CHECK(adj.ab2b.exists());
+  const int from = int(from_type);
+  const int to = int(to_type);
+
+  if (to < from) {
+    OMEGA_H_CHECK(!adj.a2ab.exists());
+    if (to_type == Topo_type::vertex) {
+      OMEGA_H_CHECK(!adj.codes.exists());
+    } else {
+      OMEGA_H_CHECK(adj.codes.exists());
+    }
+    OMEGA_H_CHECK(
+        adj.ab2b.size() == nents(from_type) * element_degree(from_type, to_type));
+  } else {
+    if (from < to) {
+      OMEGA_H_CHECK(adj.a2ab.exists());
+      OMEGA_H_CHECK(adj.codes.exists());
+      OMEGA_H_CHECK(
+          adj.ab2b.size() == nents(to_type) * element_degree(to_type, from_type));
+    }
+    OMEGA_H_CHECK(adj.a2ab.size() == nents(from_type) + 1);
+  }
+  adjs_type_[from][to] = std::make_shared<Adj>(adj);
+}
+
+Adj MixedMesh::derive_adj(Topo_type from_type, Topo_type to_type) {
+  OMEGA_H_TIME_FUNCTION;
+  check_type(from_type);
+  check_type2(to_type);
+  const int from = int(from_type);
+  const int to = int(to_type);
+  if (from < to) {
+    Adj down = ask_adj(to_type, from_type);
+    Int nlows_per_high = element_degree(to_type, from_type);
+    LO nlows = nents(from_type);
+    Adj up = invert_adj(down, nlows_per_high, nlows, to_type, from_type);
+    return up;
+  }
+  else if (to < from) {
+    OMEGA_H_CHECK(to + 1 < from);
+    Topo_type mid_type;
+    if (to_type == Topo_type::vertex) {
+      mid_type = Topo_type::edge;
+    }
+    else if ((from_type == Topo_type::tetrahedron) ||
+             (from_type == Topo_type::pyramid)) {
+      mid_type = Topo_type::triangle;
+    }
+    else {
+      mid_type = Topo_type::quadrilateral;
+    }
+    Adj h2m = ask_adj(from_type, mid_type);
+    Adj m2l = ask_adj(mid_type, to_type);
+    Adj h2l = transit(h2m, m2l, from_type, to_type, mid_type);
+    return h2l;
+  }
+  /* todo: add second order adjacency derivation */
+  Omega_h_fail("can't derive adjacency from %s to %s\n",
+      dimensional_plural_name(from_type),
+      dimensional_plural_name(to_type));
+  OMEGA_H_NORETURN(Adj());
+}
+
+Adj MixedMesh::ask_adj(Topo_type from_type, Topo_type to_type) {
+  OMEGA_H_TIME_FUNCTION;
+  check_type2(from_type);
+  check_type2(to_type);
+  if (has_adj(from_type, to_type)) {
+    return get_adj(from_type, to_type);
+  }
+  Adj derived = derive_adj(from_type, to_type);
+  adjs_type_[int(from_type)][int(to_type)] = std::make_shared<Adj>(derived);
+  return derived;
+}
+
+void MixedMesh::add_coords_mix(Reals array) {
+  add_tag<Real>(Topo_type::vertex, "coordinates", dim(), array);
+}
+
+Reals MixedMesh::coords_mix() const { return get_array<Real>(Topo_type::vertex, "coordinates"); }
+
+void get_all_type_tags(MixedMesh* mesh, Int dim, Topo_type ent_type, TagSet* tags) {
+  for (Int j = 0; j < mesh->ntags(ent_type); ++j) {
+    auto tagbase = mesh->get_tag(ent_type, j);
+    (*tags)[size_t(dim)].insert(tagbase->name());
+  }
+}
+
+#define OMEGA_H_INST(T)                                                        \
+  template Tag<T> const* MixedMesh::get_tag<T>(Topo_type ent_type, std::string const& name)    \
+      const;                                                                   \
+  template Read<T> MixedMesh::get_array<T>(Topo_type ent_type, std::string const& name) const; \
+  template void MixedMesh::add_tag<T>(                                              \
+      Topo_type ent_type, std::string const& name, Int ncomps);                           \
+  template void MixedMesh::add_tag<T>(Topo_type ent_type, std::string const& name, Int ncomps, \
+      Read<T> array, bool internal);                                           \
+  template void MixedMesh::set_tag(                                                 \
+      Topo_type ent_type, std::string const& name, Read<T> array, bool internal);
+OMEGA_H_INST(I8)
+OMEGA_H_INST(I32)
+OMEGA_H_INST(I64)
+OMEGA_H_INST(Real)
+#undef OMEGA_H_INST
+
+}  // end namespace Omega_h
diff --git a/src/Omega_h_mixedMesh.hpp b/src/Omega_h_mixedMesh.hpp
new file mode 100644
index 000000000..3f5edd08a
--- /dev/null
+++ b/src/Omega_h_mixedMesh.hpp
@@ -0,0 +1,125 @@
+#ifndef OMEGA_H_MIXEDMESH_HPP
+#define OMEGA_H_MIXEDMESH_HPP
+
+#include <Omega_h_adj.hpp>
+#include <Omega_h_comm.hpp>
+#include <Omega_h_dist.hpp>
+#include <Omega_h_library.hpp>
+#include <Omega_h_tag.hpp>
+#include <array>
+#include <map>
+#include <set>
+#include <string>
+#include <vector>
+
+namespace Omega_h {
+
+class MixedMesh {
+ public:
+  MixedMesh();
+  MixedMesh(Library* library);
+  void set_library(Library* library);
+  void set_comm(CommPtr const& new_comm);
+  void set_dim(Int dim_in);
+  void set_verts(LO nverts_in);
+  void set_ents(Int ent_dim, Adj down);
+  Library* library() const;
+  CommPtr comm() const;
+  inline Int dim() const {
+    OMEGA_H_CHECK(0 <= dim_ && dim_ <= 3);
+    return dim_;
+  }
+  inline Omega_h_Family family() const { return family_; }
+
+  void set_verts_type(LO nverts_in);
+  void set_ents(Topo_type high_type, Topo_type low_type, Adj h2l);
+  LO nents(Topo_type ent_type) const;
+  Int ent_dim(Topo_type ent_type) const;
+
+  LO npyrams() const;
+  LO nwedges() const;
+  LO nhexs() const;
+  LO ntets() const;
+  LO nquads() const;
+  LO ntris() const;
+  LO nregions_mix() const;
+  LO nfaces_mix() const;
+  LO nedges_mix() const;
+  LO nverts_mix() const;
+
+  template <typename T>
+  void add_tag(Topo_type ent_type, std::string const& name, Int ncomps);
+  template <typename T>
+  void add_tag(Topo_type ent_type, std::string const& name, Int ncomps,
+      Read<T> array, bool internal = false);
+  template <typename T>
+  void set_tag(Topo_type ent_type, std::string const& name, Read<T> array,
+      bool internal = false);
+  TagBase const* get_tagbase(Topo_type ent_type, std::string const& name) const;
+  template <typename T>
+  Tag<T> const* get_tag(Topo_type ent_type, std::string const& name) const;
+  template <typename T>
+  Read<T> get_array(Topo_type ent_type, std::string const& name) const;
+  void remove_tag(Topo_type ent_type, std::string const& name);
+  bool has_tag(Topo_type ent_type, std::string const& name) const;
+  [[nodiscard]] Int ntags(Topo_type ent_type) const;
+  TagBase const* get_tag(Topo_type ent_type, Int i) const;
+  bool has_ents(Topo_type ent_type) const;
+  bool has_adj(Topo_type from_type, Topo_type to_type) const;
+  Adj get_adj(Topo_type from_type, Topo_type to_type) const;
+  Adj ask_down(Topo_type from_type, Topo_type to_type);
+  LOs ask_verts_of(Topo_type ent_type);
+  Adj ask_up(Topo_type from_type, Topo_type to_type);
+
+  typedef std::shared_ptr<const TagBase> TagPtr;
+  typedef std::shared_ptr<const Adj> AdjPtr;
+
+ private:
+  typedef std::vector<TagPtr> TagVector;
+  typedef TagVector::iterator TagIter;
+  typedef TagVector::const_iterator TagCIter;
+  TagIter tag_iter(Topo_type ent_type, std::string const& name);
+  TagCIter tag_iter(Topo_type ent_type, std::string const& name) const;
+  void check_type(Topo_type ent_type) const;
+  void check_type2(Topo_type ent_type) const;
+  void add_adj(Topo_type from_type, Topo_type to_type, Adj adj);
+  Adj derive_adj(Topo_type from_type, Topo_type to_type);
+  Adj ask_adj(Topo_type from_type, Topo_type to_type);
+  void react_to_set_tag(Topo_type ent_type, std::string const& name);
+  Library* library_;
+  Omega_h_Family family_;
+  Int dim_;
+  CommPtr comm_;
+  LO nents_[DIMS];
+  TagVector tags_[DIMS];
+  LO nents_type_[TOPO_TYPES];
+  TagVector tags_type_[TOPO_TYPES];
+  AdjPtr adjs_type_[TOPO_TYPES][TOPO_TYPES];
+
+ public:
+  void add_coords_mix(Reals array);
+  Reals coords_mix() const;
+};
+
+using TagSet = std::array<std::set<std::string>, DIMS>;
+  
+void get_all_type_tags(MixedMesh* mesh, Int dim, Topo_type ent_type, TagSet* tags);
+
+#define OMEGA_H_EXPL_INST_DECL(T)                                              \
+  extern template Read<T> MixedMesh::get_array<T>(                                  \
+      Topo_type ent_type, std::string const& name) const;                      \
+  extern template void MixedMesh::add_tag<T>(                                       \
+      Topo_type ent_type, std::string const& name, Int ncomps);                \
+  extern template void MixedMesh::add_tag<T>(Topo_type ent_type,                    \
+      std::string const& name, Int ncomps, Read<T> array, bool internal);      \
+  extern template void MixedMesh::set_tag(Topo_type ent_type,                       \
+      std::string const& name, Read<T> array, bool internal);
+OMEGA_H_EXPL_INST_DECL(I8)
+OMEGA_H_EXPL_INST_DECL(I32)
+OMEGA_H_EXPL_INST_DECL(I64)
+OMEGA_H_EXPL_INST_DECL(Real)
+#undef OMEGA_H_EXPL_INST_DECL
+
+}  // namespace Omega_h
+
+#endif
diff --git a/src/Omega_h_vtk.cpp b/src/Omega_h_vtk.cpp
index 9dfcc2d09..2bc8103a6 100644
--- a/src/Omega_h_vtk.cpp
+++ b/src/Omega_h_vtk.cpp
@@ -19,6 +19,7 @@
 #include "Omega_h_file.hpp"
 #include "Omega_h_for.hpp"
 #include "Omega_h_mesh.hpp"
+#include "Omega_h_mixedMesh.hpp"
 #include "Omega_h_tag.hpp"
 #include "Omega_h_xml_lite.hpp"
 
@@ -42,7 +43,7 @@ TagSet get_all_vtk_tags(Mesh* mesh, Int cell_dim) {
   return tags;
 }
 
-TagSet get_all_vtk_tags_mix(Mesh* mesh, Int cell_dim) {
+TagSet get_all_vtk_tags_mix(MixedMesh* mesh, Int cell_dim) {
   TagSet tags;
   get_all_type_tags(mesh, VERT, Topo_type::vertex, &tags);
   if (cell_dim == 3) {
@@ -299,7 +300,7 @@ void write_tag_impl<Real>(
 }  // namespace detail
 
 void write_tag(std::ostream& stream, TagBase const* tag, Int space_dim,
-    Int ent_dim, Mesh* mesh, bool compress) {
+    Int ent_dim, bool compress) {
   OMEGA_H_TIME_FUNCTION;
   const auto name = tag->name();
   const auto class_ids = tag->class_ids();
@@ -458,7 +459,7 @@ static void write_piece_start_tag(
 }
 
 static void write_piece_start_tag_mix(
-    std::ostream& stream, Mesh const* mesh, Int cell_dim) {
+    std::ostream& stream, MixedMesh const* mesh, Int cell_dim) {
   stream << "<Piece NumberOfPoints=\"" << mesh->nverts_mix() << "\"";
   if (cell_dim == 3) {
     stream << " NumberOfCells=\"" << mesh->nregions_mix() << "\">\n";
@@ -490,7 +491,7 @@ static void write_connectivity(
   write_array(stream, "offsets", 1, ends, compress);
 }
 
-static void write_connectivity(std::ostream& stream, Mesh* mesh, Int cell_dim,
+static void write_connectivity(std::ostream& stream, MixedMesh* mesh, Int cell_dim,
     Topo_type max_type, bool compress) {
   if (cell_dim == 3) {
     Read<I8> types_t(mesh->nents(Topo_type::tetrahedron),
@@ -721,8 +722,8 @@ filesystem::path get_pvd_path(filesystem::path const& root_path) {
   return result;
 }
 
-static void default_dim(Mesh* mesh, Int* cell_dim) {
-  if (*cell_dim == -1) *cell_dim = mesh->dim();
+static void default_dim(Int meshDim, Int* cell_dim) {
+  if (*cell_dim == -1) *cell_dim = meshDim;
 }
 
 static void verify_vtk_tagset(Mesh* mesh, Int cell_dim, TagSet const& tags) {
@@ -768,7 +769,7 @@ void write_vtkfile_vtu_start_tag(std::ostream& stream, bool compress) {
 void write_vtu(std::ostream& stream, Mesh* mesh, Int cell_dim,
     TagSet const& tags, bool compress) {
   OMEGA_H_TIME_FUNCTION;
-  default_dim(mesh, &cell_dim);
+  default_dim(mesh->dim(), &cell_dim);
   verify_vtk_tagset(mesh, cell_dim, tags);
   write_vtkfile_vtu_start_tag(stream, compress);
   stream << "<UnstructuredGrid>\n";
@@ -785,14 +786,14 @@ void write_vtu(std::ostream& stream, Mesh* mesh, Int cell_dim,
   /* globals go first so read_vtu() knows where to find them */
   if (mesh->has_tag(VERT, "global") && tags[VERT].count("global")) {
     write_tag(stream, mesh->get_tag<GO>(VERT, "global"), mesh->dim(), VERT,
-        mesh, compress);
+        compress);
   }
   write_locals_and_owners(stream, mesh, VERT, tags, compress);
   for (Int i = 0; i < mesh->ntags(VERT); ++i) {
     auto tag = mesh->get_tag(VERT, i);
     if (tag->name() != "coordinates" && tag->name() != "global" &&
         tags[VERT].count(tag->name())) {
-      write_tag(stream, tag, mesh->dim(), VERT, mesh, compress);
+      write_tag(stream, tag, mesh->dim(), VERT, compress);
     }
   }
   stream << "</PointData>\n";
@@ -801,7 +802,7 @@ void write_vtu(std::ostream& stream, Mesh* mesh, Int cell_dim,
   if (mesh->has_tag(cell_dim, "global") &&
       tags[size_t(cell_dim)].count("global")) {
     write_tag(stream, mesh->get_tag<GO>(cell_dim, "global"), mesh->dim(),
-        cell_dim, mesh, compress);
+        cell_dim, compress);
   }
   write_locals_and_owners(stream, mesh, cell_dim, tags, compress);
   if (tags[size_t(cell_dim)].count("vtkGhostType")) {
@@ -810,7 +811,7 @@ void write_vtu(std::ostream& stream, Mesh* mesh, Int cell_dim,
   for (Int i = 0; i < mesh->ntags(cell_dim); ++i) {
     auto tag = mesh->get_tag(cell_dim, i);
     if (tag->name() != "global" && tags[size_t(cell_dim)].count(tag->name())) {
-      write_tag(stream, tag, mesh->dim(), cell_dim, mesh, compress);
+      write_tag(stream, tag, mesh->dim(), cell_dim, compress);
     }
   }
   stream << "</CellData>\n";
@@ -819,14 +820,14 @@ void write_vtu(std::ostream& stream, Mesh* mesh, Int cell_dim,
   stream << "</VTKFile>\n";
 }
 
-void write_vtu(filesystem::path const& filename, Mesh* mesh, Topo_type max_type,
+void write_vtu(filesystem::path const& filename, MixedMesh* mesh, Topo_type max_type,
     bool compress) {
   auto tags = get_all_vtk_tags_mix(mesh, mesh->dim());
   OMEGA_H_TIME_FUNCTION;
   std::ofstream stream(filename.c_str());
   OMEGA_H_CHECK(stream.is_open());
   auto cell_dim = mesh->ent_dim(max_type);
-  default_dim(mesh, &cell_dim);
+  default_dim(mesh->dim(), &cell_dim);
   write_vtkfile_vtu_start_tag(stream, compress);
   stream << "<UnstructuredGrid>\n";
   write_piece_start_tag_mix(stream, mesh, cell_dim);
@@ -840,33 +841,26 @@ void write_vtu(filesystem::path const& filename, Mesh* mesh, Topo_type max_type,
       compress);
   stream << "</Points>\n";
   stream << "<PointData>\n";
-  if (mesh->has_tag(VERT, "global") && tags[VERT].count("global")) {
-    write_tag(stream, mesh->get_tag<GO>(VERT, "global"), mesh->dim(), 0, mesh,
+  if (mesh->has_tag(Topo_type::vertex, "global") && tags[VERT].count("global")) {
+    write_tag(stream, mesh->get_tag<GO>(Topo_type::vertex, "global"), mesh->dim(), 0,
         compress);
   }
   for (Int i = 0; i < mesh->ntags(Topo_type::vertex); ++i) {
     auto tag = mesh->get_tag(Topo_type::vertex, i);
     if (tag->name() != "coordinates" && tag->name() != "global" &&
         tags[VERT].count(tag->name())) {
-      write_tag(stream, tag, mesh->dim(), 0, mesh, compress);
+      write_tag(stream, tag, mesh->dim(), 0, compress);
     }
   }
   stream << "</PointData>\n";
   stream << "<CellData>\n";
-  if (mesh->has_tag(cell_dim, "global") &&
-      tags[size_t(cell_dim)].count("global")) {
-    write_tag(stream, mesh->get_tag<GO>(cell_dim, "global"), mesh->dim(),
-        cell_dim, mesh, compress);
-  }
-  if (tags[size_t(cell_dim)].count("vtkGhostType")) {
-    write_vtk_ghost_types(stream, mesh, cell_dim, compress);
-  }
+  OMEGA_H_CHECK(cell_dim <= 3);
   if (cell_dim == 3) {
     for (Int i = 0; i < mesh->ntags(Topo_type::tetrahedron); ++i) {
       auto tag = mesh->get_tag(Topo_type::tetrahedron, i);
       if (tag->name() != "global" &&
           tags[size_t(cell_dim)].count(tag->name())) {
-        write_tag(stream, tag, mesh->dim(), cell_dim, mesh, compress);
+        write_tag(stream, tag, mesh->dim(), cell_dim, compress);
       }
     }
 
@@ -874,7 +868,7 @@ void write_vtu(filesystem::path const& filename, Mesh* mesh, Topo_type max_type,
       auto tag = mesh->get_tag(Topo_type::hexahedron, i);
       if (tag->name() != "global" &&
           tags[size_t(cell_dim)].count(tag->name())) {
-        write_tag(stream, tag, mesh->dim(), cell_dim, mesh, compress);
+        write_tag(stream, tag, mesh->dim(), cell_dim, compress);
       }
     }
 
@@ -882,7 +876,7 @@ void write_vtu(filesystem::path const& filename, Mesh* mesh, Topo_type max_type,
       auto tag = mesh->get_tag(Topo_type::wedge, i);
       if (tag->name() != "global" &&
           tags[size_t(cell_dim)].count(tag->name())) {
-        write_tag(stream, tag, mesh->dim(), cell_dim, mesh, compress);
+        write_tag(stream, tag, mesh->dim(), cell_dim, compress);
       }
     }
 
@@ -890,7 +884,7 @@ void write_vtu(filesystem::path const& filename, Mesh* mesh, Topo_type max_type,
       auto tag = mesh->get_tag(Topo_type::pyramid, i);
       if (tag->name() != "global" &&
           tags[size_t(cell_dim)].count(tag->name())) {
-        write_tag(stream, tag, mesh->dim(), cell_dim, mesh, compress);
+        write_tag(stream, tag, mesh->dim(), cell_dim, compress);
       }
     }
   } else if (cell_dim == 2) {
@@ -898,7 +892,7 @@ void write_vtu(filesystem::path const& filename, Mesh* mesh, Topo_type max_type,
       auto tag = mesh->get_tag(Topo_type::triangle, i);
       if (tag->name() != "global" &&
           tags[size_t(cell_dim)].count(tag->name())) {
-        write_tag(stream, tag, mesh->dim(), cell_dim, mesh, compress);
+        write_tag(stream, tag, mesh->dim(), cell_dim, compress);
       }
     }
 
@@ -906,15 +900,16 @@ void write_vtu(filesystem::path const& filename, Mesh* mesh, Topo_type max_type,
       auto tag = mesh->get_tag(Topo_type::quadrilateral, i);
       if (tag->name() != "global" &&
           tags[size_t(cell_dim)].count(tag->name())) {
-        write_tag(stream, tag, mesh->dim(), cell_dim, mesh, compress);
+        write_tag(stream, tag, mesh->dim(), cell_dim, compress);
       }
     }
-  } else {
-    for (Int i = 0; i < mesh->ntags(cell_dim); ++i) {
-      auto tag = mesh->get_tag(cell_dim, i);
+  } else { //cell_dim == 1 or 0
+    auto cellTopoType = static_cast<Topo_type>(cell_dim);
+    for (Int i = 0; i < mesh->ntags(cellTopoType); ++i) {
+      auto tag = mesh->get_tag(cellTopoType, i);
       if (tag->name() != "global" &&
           tags[size_t(cell_dim)].count(tag->name())) {
-        write_tag(stream, tag, mesh->dim(), cell_dim, mesh, compress);
+        write_tag(stream, tag, mesh->dim(), cell_dim, compress);
       }
     }
   }
@@ -1003,7 +998,7 @@ void write_vtu(filesystem::path const& filename, Mesh* mesh, Int cell_dim,
 
 void write_vtu(
     std::string const& filename, Mesh* mesh, Int cell_dim, bool compress) {
-  default_dim(mesh, &cell_dim);
+  default_dim(mesh->dim(), &cell_dim);
   write_vtu(
       filename, mesh, cell_dim, get_all_vtk_tags(mesh, cell_dim), compress);
 }
@@ -1129,7 +1124,7 @@ void read_pvtu(filesystem::path const& pvtupath, CommPtr comm, I32* npieces_out,
 void write_parallel(filesystem::path const& path, Mesh* mesh, Int cell_dim,
     TagSet const& tags, bool compress) {
   ScopedTimer timer("vtk::write_parallel");
-  default_dim(mesh, &cell_dim);
+  default_dim(mesh->dim(), &cell_dim);
   ask_for_mesh_tags(mesh, tags);
   auto const rank = mesh->comm()->rank();
   if (rank == 0) {
@@ -1152,7 +1147,7 @@ void write_parallel(filesystem::path const& path, Mesh* mesh, Int cell_dim,
 
 void write_parallel(
     std::string const& path, Mesh* mesh, Int cell_dim, bool compress) {
-  default_dim(mesh, &cell_dim);
+  default_dim(mesh->dim(), &cell_dim);
   ScopedChangeRCFieldsToMesh rc_to_mesh(*mesh);
   write_parallel(
       path, mesh, cell_dim, get_all_vtk_tags(mesh, cell_dim), compress);
@@ -1288,7 +1283,7 @@ Writer::Writer(filesystem::path const& root_path, Mesh* mesh, Int cell_dim,
       compress_(compress),
       step_(0),
       pvd_pos_(0) {
-  default_dim(mesh_, &cell_dim_);
+  default_dim(mesh_->dim(), &cell_dim_);
   auto const comm = mesh->comm();
   auto const rank = comm->rank();
   if (rank == 0) {
diff --git a/src/meshsim2osh.cpp b/src/meshsim2osh.cpp
index 38929dfb0..7077fc265 100644
--- a/src/meshsim2osh.cpp
+++ b/src/meshsim2osh.cpp
@@ -25,10 +25,14 @@ int main(int argc, char** argv) {
     std::cout << "attaching numbering...\n";
     numbering_in = cmdline.get<std::string>("-numbering", "numbering-in");
   }
-  auto mesh = Omega_h::meshsim::read(mesh_in, model_in, numbering_in, comm);
-  auto family = mesh.family();
-  if ((family == OMEGA_H_SIMPLEX) || family == OMEGA_H_HYPERCUBE) {
+  auto isMixed = Omega_h::meshsim::isMixed(mesh_in, model_in);
+  std::cerr << "isMixed " << isMixed << "\n";
+  //TODO - call the correct reader (mixed vs mono)
+  if( !isMixed ) {
+    auto mesh = Omega_h::meshsim::read(mesh_in, model_in, numbering_in, comm);
     Omega_h::binary::write(mesh_out, &mesh);
+  } else {
+    auto mesh = Omega_h::meshsim::readMixed(mesh_in, model_in, comm);
   }
 
   return 0;
diff --git a/src/mixed_test.cpp b/src/mixed_test.cpp
index 4fae8f437..8f6269e08 100644
--- a/src/mixed_test.cpp
+++ b/src/mixed_test.cpp
@@ -2,7 +2,7 @@
 #include <string>
 
 #include "Omega_h_file.hpp"
-#include "Omega_h_mesh.hpp"
+#include "Omega_h_mixedMesh.hpp"
 
 #include "Omega_h_element.hpp"
 #include "Omega_h_array_ops.hpp"
@@ -17,7 +17,7 @@ void test_degree() {
   OMEGA_H_CHECK(element_degree(Topo_type::pyramid, Topo_type::quadrilateral) == 1);
 }
 
-void test_tags(Mesh* mesh) {
+void test_tags(MixedMesh* mesh) {
   auto num_wedge = mesh->nwedges();
   mesh->add_tag<LO>(Topo_type::wedge, "gravity", 1);
   mesh->set_tag<LO>(Topo_type::wedge, "gravity", LOs(num_wedge,10));
@@ -32,7 +32,7 @@ void test_tags(Mesh* mesh) {
   mesh->remove_tag(Topo_type::pyramid, "density");
 }
 
-void test_adjs(Mesh* mesh) {
+void test_adjs(MixedMesh* mesh) {
 
   //get number of entities
   auto num_vertex = mesh->nverts_mix();
@@ -157,29 +157,26 @@ void test_adjs(Mesh* mesh) {
 void test_finer_meshes(CommPtr comm, std::string mesh_dir) {
 
   //tet=4609, hex=249, wedge=262, pyramid=313
-  std::string mesh_in = std::string(mesh_dir) +
-                        "/localconcave_tutorial_mixedvol_geomsim3-case1.sms";
-  std::string model_in = std::string(mesh_dir) +
-                         "/localconcave_tutorial_mixedvol_geomsim3.smd";
-  auto mesh = meshsim::read(mesh_in, model_in, comm);
+  auto name = std::string("/localconcave_tutorial_mixedvol_geomsim3");
+  std::string mesh_in = std::string(mesh_dir) + name + "-case1.sms";
+  std::string model_in = std::string(mesh_dir) + name + ".smd";
+  auto mesh = meshsim::readMixed(mesh_in, model_in, comm);
   test_adjs(&mesh);
   test_tags(&mesh);
 
   //tet=1246, hex=8, wedge=0, pyramid=229
-  mesh_in = std::string(mesh_dir) +
-            "/localconcave_tutorial_mixedvol_geomsim2-case1.sms";
-  model_in = std::string(mesh_dir) +
-             "/localconcave_tutorial_mixedvol_geomsim2.smd";
-  mesh = meshsim::read(mesh_in, model_in, comm);
+  name = std::string("/localconcave_tutorial_mixedvol_geomsim2");
+  mesh_in = std::string(mesh_dir) + name + "-case1.sms";
+  model_in = std::string(mesh_dir) + name + ".smd";
+  mesh = meshsim::readMixed(mesh_in, model_in, comm);
   test_adjs(&mesh);
   test_tags(&mesh);
 
   //tet=3274, hex=68, wedge=0, pyramid=171
-  mesh_in = std::string(mesh_dir) +
-            "/localconcave_tutorial_mixedvol_geomsim-case1.sms";
-  model_in = std::string(mesh_dir) +
-             "/localconcave_tutorial_mixedvol_geomsim.smd";
-  mesh = meshsim::read(mesh_in, model_in, comm);
+  name = std::string("/localconcave_tutorial_mixedvol_geomsim");
+  mesh_in = std::string(mesh_dir) + name + "-case1.sms";
+  model_in = std::string(mesh_dir) + name + ".smd";
+  mesh = meshsim::readMixed(mesh_in, model_in, comm);
   test_adjs(&mesh);
   test_tags(&mesh);
 
@@ -188,7 +185,7 @@ void test_finer_meshes(CommPtr comm, std::string mesh_dir) {
             "/localconcave_tutorial-case1_v7.sms";
   model_in = std::string(mesh_dir) +
              "/localconcave_tutorial_geomsim.smd";
-  mesh = meshsim::read(mesh_in, model_in, comm);
+  auto simplexMesh = meshsim::read(mesh_in, model_in, comm);
 
 }
 
@@ -204,37 +201,39 @@ int main(int argc, char** argv) {
 
   /* Generate these meshes using the mixed_writeMesh ctest */
   //TODO: add support for writing mixed meshes containing tags(class info) to vtk
+  {
   std::string mesh_in = std::string(mesh_dir) + "/Example_hex.sms";
   std::string model_in = std::string(mesh_dir) + "/Example_hex.smd";
   auto mesh = meshsim::read(mesh_in, model_in, comm);
+  }
 
-  mesh_in = std::string(mesh_dir) + "/Example_wedge.sms";
-  model_in = std::string(mesh_dir) + "/Example_wedge.smd";
-  mesh = meshsim::read(mesh_in, model_in, comm);
+  auto mesh_in = std::string(mesh_dir) + "/Example_wedge.sms";
+  auto model_in = std::string(mesh_dir) + "/Example_wedge.smd";
+  auto mesh = meshsim::readMixed(mesh_in, model_in, comm);
   test_adjs(&mesh);
   test_tags(&mesh);
 
   mesh_in = std::string(mesh_dir) + "/Example_pym.sms";
   model_in = std::string(mesh_dir) + "/Example_pym.smd";
-  mesh = meshsim::read(mesh_in, model_in, comm);
+  mesh = meshsim::readMixed(mesh_in, model_in, comm);
   test_adjs(&mesh);
   test_tags(&mesh);
 
   mesh_in = std::string(mesh_dir) + "/Example_tet_wedge.sms";
   model_in = std::string(mesh_dir) + "/Example_tet_wedge.smd";
-  mesh = meshsim::read(mesh_in, model_in, comm);
+  mesh = meshsim::readMixed(mesh_in, model_in, comm);
   test_adjs(&mesh);
   test_tags(&mesh);
 
   mesh_in = std::string(mesh_dir) + "/Example_pym_hex.sms";
   model_in = std::string(mesh_dir) + "/Example_pym_hex.smd";
-  mesh = meshsim::read(mesh_in, model_in, comm);
+  mesh = meshsim::readMixed(mesh_in, model_in, comm);
   test_adjs(&mesh);
   test_tags(&mesh);
 
   mesh_in = std::string(mesh_dir) + "/Example_allType.sms";
   model_in = std::string(mesh_dir) + "/Example_allType.smd";
-  mesh = meshsim::read(mesh_in, model_in, comm);
+  mesh = meshsim::readMixed(mesh_in, model_in, comm);
   test_adjs(&mesh);
   test_tags(&mesh);