diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9b48f1d71..621afc77b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -105,14 +105,18 @@ set(Trilinos_PREFIX "" CACHE STRING "Trilinos installation directory")
 option(SKIP_SIMMETRIX_VERSION_CHECK "enable at your own risk; it may result in undefined behavior" OFF)
 option(ENABLE_SIMMETRIX "Build with Simmetrix support" OFF)
 message(STATUS "ENABLE_SIMMETRIX: ${ENABLE_SIMMETRIX}")
-option(ENABLE_CAPSTONE "Build with Capstone support" OFF)
-message(STATUS "ENABLE_CAPSTONE: ${ENABLE_CAPSTONE}")
+option(PUMI_ENABLE_CAPSTONE "Build PUMI with Capstone support" OFF)
+# For compatibility and to inherit from owning projects:
+if(ENABLE_CAPSTONE)
+  set(PUMI_ENABLE_CAPSTONE ON)
+endif()
+message(STATUS "PUMI_ENABLE_CAPSTONE: ${PUMI_ENABLE_CAPSTONE}")
 
 if(ENABLE_SIMMETRIX)
   add_definitions(-DHAVE_SIMMETRIX)
 endif()
-if(ENABLE_CAPSTONE)
-  add_definitions(-DHAVE_CAPSTONE)
+if(PUMI_ENABLE_CAPSTONE)
+  add_definitions(-DPUMI_HAS_CAPSTONE)
 endif()
 
 option(ENABLE_FPP "Build with snapping to first problem plane" OFF)
@@ -144,8 +148,8 @@ if(ENABLE_SIMMETRIX)
   find_package(SimModSuite MODULE REQUIRED)
 endif()
 
-if(ENABLE_CAPSTONE)
-  set(SCOREC_USE_CreateMG_DEFAULT ${ENABLE_CAPSTONE})
+if(PUMI_ENABLE_CAPSTONE)
+  set(SCOREC_USE_CreateMG_DEFAULT ${PUMI_ENABLE_CAPSTONE})
   bob_public_dep(CreateMG)
 endif()
 
diff --git a/apf/apfMesh.cc b/apf/apfMesh.cc
index 7f57fc983..80131a6e4 100644
--- a/apf/apfMesh.cc
+++ b/apf/apfMesh.cc
@@ -51,7 +51,7 @@ int const quad_edge_verts[4][2] =
 */
 
 int const tet_edge_verts[6][2] =
-#ifndef HAVE_CAPSTONE
+#ifndef PUMI_HAS_CAPSTONE
 {{0,1}
 ,{1,2}
 ,{2,0}
@@ -78,7 +78,7 @@ int const pyramid_edge_verts[8][2] =
 ,{0,4},{1,4},{2,4},{3,4}};
 
 int const tet_tri_verts[4][3] =
-#ifndef HAVE_CAPSTONE
+#ifndef PUMI_HAS_CAPSTONE
 {{0,1,2}
 ,{0,1,3}
 ,{1,2,3}
diff --git a/apf_cap/CMakeLists.txt b/apf_cap/CMakeLists.txt
index 785a38a6b..7ae827496 100644
--- a/apf_cap/CMakeLists.txt
+++ b/apf_cap/CMakeLists.txt
@@ -3,7 +3,7 @@ if(DEFINED TRIBITS_PACKAGE)
   return()
 endif()
 
-if(NOT ENABLE_CAPSTONE)
+if(NOT PUMI_ENABLE_CAPSTONE)
   return()
 endif()
 
@@ -12,25 +12,35 @@ include(CheckIncludeFileCXX)
 cmake_policy(SET CMP0075 NEW) # Observe CMAKE_REQUIRED_LIBRARIES.
 
 #Sources & Headers
-set(SOURCES apfCAP.cc)
+set(SOURCES apfCAP.cc apfCAPsizing.cc)
 set(HEADERS apfCAP.h)
 
 #Library
 add_library(apf_cap ${SOURCES})
 target_link_libraries(apf_cap PUBLIC apf gmi_cap)
-target_link_libraries(apf_cap PUBLIC capstone_module
-				     framework_testing)
+target_link_libraries(apf_cap PRIVATE framework_mesh framework_application)
 
 set(CMAKE_CXX_OLD_STANDARD "${CMAKE_CXX_STANDARD}")
 cmake_push_check_state(RESET)
 set(CMAKE_CXX_STANDARD 14)
+
+# CreateFT_Vector.h via CreateMG_SizingMetricTool.h is missing
+# `#include <limits>` in Capstone 12.5.2 so add it explicitly in this check.
+if(MSVC)
+  set(CMAKE_REQUIRED_FLAGS "/FI limits")
+else()
+  set(CMAKE_REQUIRED_FLAGS "-include limits")
+endif()
+
 set(CMAKE_REQUIRED_LIBRARIES "framework_meshing")
-check_include_file_cxx("CreateMG_SizingMetricTool.h" HAVE_CAPSTONE_SIZINGMETRICTOOL)
+check_include_file_cxx(
+  "CreateMG_SizingMetricTool.h" PUMI_HAS_CAPSTONE_SIZINGMETRICTOOL
+)
 cmake_pop_check_state()
 set(CMAKE_CXX_STANDARD "${CMAKE_CXX_OLD_STANDARD}")
 
-if(HAVE_CAPSTONE_SIZINGMETRICTOOL)
-target_compile_definitions(apf_cap PRIVATE HAVE_CAPSTONE_SIZINGMETRICTOOL)
+if(PUMI_HAS_CAPSTONE_SIZINGMETRICTOOL)
+target_compile_definitions(apf_cap PRIVATE PUMI_HAS_CAPSTONE_SIZINGMETRICTOOL)
 target_compile_features(framework_meshing INTERFACE cxx_std_14)
 target_link_libraries(apf_cap PRIVATE framework_meshing)
 endif()
diff --git a/apf_cap/apfCAP.cc b/apf_cap/apfCAP.cc
index 13919376c..88ce841e2 100644
--- a/apf_cap/apfCAP.cc
+++ b/apf_cap/apfCAP.cc
@@ -1,17 +1,25 @@
 #include "apfCAP.h"
+
+#include <algorithm>
+#include <cstdlib>
+
 #include <apf.h>
+#include <apfMesh2.h>
 #include <apfShape.h>
 #include <gmi.h>
 #include <gmi_cap.h>
-#include <cstdlib>
+#include <lionPrint.h>
 #include <mth.h>
 #include <pcu_util.h>
-#include <algorithm>
-#include <lionPrint.h>
+#include <PCU.h>
 
-#ifdef HAVE_CAPSTONE_SIZINGMETRICTOOL
-#include <CreateMG_SizingMetricTool.h>
-#endif
+#include <CreateMG_AppProcessor.h>
+#include <CreateMG_Framework_Geometry.h>
+#include <CreateMG_Framework_Mesh.h>
+#include <CreateMG_Function.h>
+
+using namespace CreateMG;
+namespace MeshMG = CreateMG::Mesh;
 
 namespace apf {
 
@@ -63,8 +71,9 @@ int const degree[Mesh::TYPES][4] =
 ,{5, 8,5,1} /* MDS_PYRAMID */
 };
 
-
-
+bool hasCAP() noexcept {
+  return true;
+}
 
 MeshEntity* toEntity(M_MTopo topo)
 {
@@ -81,21 +90,182 @@ M_MTopo fromEntity(MeshEntity* e)
   return topo;
 }
 
-MeshCAP::MeshCAP(MeshDatabaseInterface* mdb, GeometryDatabaseInterface* gdb):
-  meshInterface(mdb), geomInterface(gdb)
+class TagCAP;
+
+class MeshCAP : public Mesh2
 {
+  public:
+    MeshCAP(MDBI* mdb, GDBI* gdb);
+    MeshCAP(gmi_model* mdl, MDBIP mdb);
+    virtual ~MeshCAP();
+    /* --------------------------------------------------------------------- */
+    /* Category 00: General Mesh APIs */
+    /* --------------------------------------------------------------------- */
+    // REQUIRED Member Functions //
+    int getDimension();
+    std::size_t count(int dimension);
+    Type getType(MeshEntity* e);
+    void verify();
+    // OPTIONAL Member Functions //
+    void writeNative(const char* fileName);
+    void destroyNative();
+
+    /* --------------------------------------------------------------------- */
+    /* Category 01: General Getters and Setters for vertex coordinates */
+    /* --------------------------------------------------------------------- */
+    // REQUIRED Member Functions //
+    void getPoint_(MeshEntity* e, int, Vector3& point);
+    void setPoint_(MeshEntity* e, int, Vector3 const& p);
+    void getParam(MeshEntity* e, Vector3& p);
+    void setParam(MeshEntity* e, Vector3 const& point);
+
+    /* --------------------------------------------------------------------- */
+    /* Category 02: Iterators */
+    /* --------------------------------------------------------------------- */
+    // REQUIRED Member Functions //
+    MeshIterator* begin(int dimension);
+    MeshEntity* iterate(MeshIterator* it);
+    void end(MeshIterator* it);
+    void increment(MeshIterator* it);
+    bool isDone(MeshIterator* it);
+    MeshEntity* deref(MeshIterator* it);
+
+    /* --------------------------------------------------------------------- */
+    /* Category 03: Adjacencies */
+    /* --------------------------------------------------------------------- */
+    // REQUIRED Member Functions //
+    void getAdjacent(MeshEntity* e, int dimension, Adjacent& adjacent);
+    int getDownward(MeshEntity* e, int dimension, MeshEntity** adjacent);
+    MeshEntity* getUpward(MeshEntity* e, int i);
+    bool hasUp(MeshEntity* e);
+    // OPTIONAL Member Functions //
+    bool hasAdjacency(int from_dim, int to_dim);
+    void createAdjacency(int from_dim, int to_dim);
+    void deleteAdjacency(int from_dim, int to_dim);
+    void getUp(MeshEntity* e, Up& up);
+    int countUpward(MeshEntity* e);
+
+    /* --------------------------------------------------------------------- */
+    /* Category 04: CAD model inquires */
+    /* --------------------------------------------------------------------- */
+    // REQUIRED Member Functions //
+    ModelEntity* toModel(MeshEntity* e);
+    // OPTIONAL Member Functions //
+    gmi_model* getModel();
+    void setModel(gmi_model* newModel);
+    void setModelEntity(MeshEntity* e, ModelEntity* me);
+
+    /* --------------------------------------------------------------------- */
+    /* Category 05: Entity Creation/Deletion */
+    /* --------------------------------------------------------------------- */
+    // REQUIRED Member Functions //
+    MeshEntity* createVert_(ModelEntity* me);
+    MeshEntity* createEntity_(int type, ModelEntity* me, MeshEntity** down);
+    void destroy_(MeshEntity* e);
+
+    /* --------------------------------------------------------------------- */
+    /* Category 06: Attachable Data Functionality */
+    /* --------------------------------------------------------------------- */
+    // REQUIRED Member Functions //
+    MeshTag* createDoubleTag(const char* name, int size);
+    MeshTag* createIntTag(const char* name, int size);
+    MeshTag* createLongTag(const char* name, int size);
+    MeshTag* findTag(const char* name);
+    void destroyTag(MeshTag* t);
+    void renameTag(MeshTag* t, const char* name);
+    void getTags(DynamicArray<MeshTag*>& tags);
+    /* void getTag(MeshEntity* e, MeshTag* t, void* data); */
+    /* void setTag(MeshEntity* e, MeshTag* t, void const* data); */
+    void getDoubleTag(MeshEntity* e, MeshTag* tag, double* data);
+    void setDoubleTag(MeshEntity* e, MeshTag* tag, double const* data);
+    void getIntTag(MeshEntity* e, MeshTag* tag, int* data);
+    void setIntTag(MeshEntity* e, MeshTag* tag, int const* data);
+    void getLongTag(MeshEntity* e, MeshTag* tag, long* data);
+    void setLongTag(MeshEntity* e, MeshTag* tag, long const* data);
+    void removeTag(MeshEntity* e, MeshTag* t);
+    bool hasTag(MeshEntity* e, MeshTag* t);
+    int getTagType(MeshTag* t);
+    int getTagSize(MeshTag* t);
+    const char* getTagName(MeshTag* t);
+    unsigned getTagChecksum(MeshTag*,int);
+
+
+    /* --------------------------------------------------------------------- */
+    /* Category 07: Distributed Meshes */
+    /* --------------------------------------------------------------------- */
+    // REQUIRED Member Functions //
+    bool isShared(MeshEntity* e);
+    bool isGhost(MeshEntity*) { return false; }
+    bool isGhosted(MeshEntity*) { return false; }
+    bool isOwned(MeshEntity* e);
+    int getOwner(MeshEntity* e);
+    void getRemotes(MeshEntity* e, Copies& remotes);
+    void getResidence(MeshEntity* e, Parts& residence);
+    int getId();
+    void setResidence(MeshEntity*, Parts&) {}
+    void acceptChanges() {}
+    // OPTIONAL Member Functions //
+    void deleteGhost(MeshEntity*) {}
+    void addGhost(MeshEntity*, int, MeshEntity*) {}
+    int getGhosts(MeshEntity*, Copies&) { return 0; }
+    void migrate(Migration* plan);
+    void setRemotes(MeshEntity*, Copies&) {}
+    void addRemote(MeshEntity*, int, MeshEntity*) {}
+    void clearRemotes(MeshEntity*) {}
+
+
+    /* --------------------------------------------------------------------- */
+    /* Category 08: Periodic Meshes */
+    /* --------------------------------------------------------------------- */
+    // REQUIRED Member Functions //
+    bool hasMatching() { return false; }
+    void getMatches(MeshEntity* e, Matches& m);
+    // OPTIONAL Member Functions //
+    void addMatch(MeshEntity*, int, MeshEntity*) {}
+    void clearMatches(MeshEntity*) {}
+    void clear_() {}
+    void getDgCopies(MeshEntity* e, DgCopies& dgCopies, ModelEntity* me);
+
+    MDBI* getMesh() noexcept { return meshInterface; }
+    void disownModel() noexcept { ownsModel = false; }
+  private:
+    void setupAdjacencies();
+    gmi_model* model{nullptr};
+    MDBIP meshOwner;
+    MDBI* meshInterface;
+    int d;
+    bool ownsModel{true};
+    std::vector<TagCAP*> tags;
+};
+
+void MeshCAP::setupAdjacencies() {
   PCU_ALWAYS_ASSERT(meshInterface);
+  MG_API_CALL(meshInterface, get_dimension(d));
+  MG_API_CALL(meshInterface, set_adjacency_state(
+    REGION2FACE|REGION2EDGE|REGION2VERTEX| FACE2EDGE|FACE2VERTEX
+  ));
+  MG_API_CALL(meshInterface, set_reverse_states());
+  MG_API_CALL(meshInterface, compute_adjacency());
+}
 
-  std::size_t numRegions = 0;
-  meshInterface->get_num_topos(TOPO_REGION, numRegions);
-  d = numRegions ? 3 : 2;
-  iterDim = -1;
-  model = gmi_import_cap(geomInterface);
+MeshCAP::MeshCAP(MDBI* mdb, GDBI* gdb): meshInterface(mdb) {
+  model = gmi_import_cap(gdb);
+  setupAdjacencies();
+}
+
+MeshCAP::MeshCAP(gmi_model* mdl, MDBIP mdb):
+  model(mdl), meshOwner(mdb), meshInterface(mdb.get())
+{
+  PCU_ALWAYS_ASSERT(mdl);
+  PCU_ALWAYS_ASSERT(mdb.is_valid());
+  setupAdjacencies();
 }
 
 MeshCAP::~MeshCAP()
 {
-  if (model) destroyNative();
+  if (ownsModel && model) gmi_destroy(model);
+  model = nullptr;
+  meshInterface = nullptr;
 }
 
 int MeshCAP::getDimension()
@@ -107,37 +277,29 @@ std::size_t MeshCAP::count(int dimension)
 {
   std::size_t count = 0;
   if (dimension == 0)
-    meshInterface->get_num_topos(TOPO_VERTEX, count);
+    meshInterface->get_num_topos(MeshMG::TOPO_VERTEX, count);
   if (dimension == 1)
-    meshInterface->get_num_topos(TOPO_EDGE, count);
+    meshInterface->get_num_topos(MeshMG::TOPO_EDGE, count);
   if (dimension == 2)
-    meshInterface->get_num_topos(TOPO_FACE, count);
+    meshInterface->get_num_topos(MeshMG::TOPO_FACE, count);
   if (dimension == 3)
-    meshInterface->get_num_topos(TOPO_REGION, count);
+    meshInterface->get_num_topos(MeshMG::TOPO_REGION, count);
   return count;
 }
 
 Mesh::Type MeshCAP::getType(MeshEntity* e)
 {
   M_MTopo topo = fromEntity(e);
-  MeshShape topoShape;
+  MeshMG::MeshShape topoShape;
   meshInterface->get_topo_shape(topo, topoShape);
-  if (topoShape == SHAPE_NODE)
-    return Mesh::VERTEX;
-  else if (topoShape == SHAPE_SEGMENT)
-    return Mesh::EDGE;
-  else if (topoShape == SHAPE_TRIANGLE)
-    return Mesh::TRIANGLE;
-  else if (topoShape == SHAPE_QUAD)
-    return Mesh::QUAD;
-  else if (topoShape == SHAPE_TETRA)
-    return Mesh::TET;
-  else if (topoShape == SHAPE_HEX)
-    return Mesh::HEX;
-  else if (topoShape == SHAPE_PRISM)
-    return Mesh::PRISM;
-  else if (topoShape == SHAPE_PYRAMID)
-    return Mesh::PYRAMID;
+  if (topoShape == MeshMG::SHAPE_NODE) return VERTEX;
+  else if (topoShape == MeshMG::SHAPE_SEGMENT) return EDGE;
+  else if (topoShape == MeshMG::SHAPE_TRIANGLE) return TRIANGLE;
+  else if (topoShape == MeshMG::SHAPE_QUAD) return QUAD;
+  else if (topoShape == MeshMG::SHAPE_TETRA) return TET;
+  else if (topoShape == MeshMG::SHAPE_HEX) return HEX;
+  else if (topoShape == MeshMG::SHAPE_PRISM) return PRISM;
+  else if (topoShape == MeshMG::SHAPE_PYRAMID) return PYRAMID;
   else
     apf::fail("MeshCAP::getType encountered an unknown entity type!\n");
 }
@@ -184,7 +346,7 @@ void MeshCAP::getParam(MeshEntity* e, Vector3& point)
 {
   M_MTopo topo = fromEntity(e);
   double u, v;
-  GeometryTopoType gtype;
+  MeshMG::GeometryTopoType gtype;
   meshInterface->get_vertex_uv_parameters(topo, u, v, gtype);
   point = Vector3(u, v, 0.);
 }
@@ -198,15 +360,15 @@ void MeshCAP::setParam(MeshEntity* e, Vector3 const& point)
 
 MeshIterator* MeshCAP::begin(int dimension)
 {
-  MeshSmartIterator* miter = new MeshSmartIterator(meshInterface);
+  auto miter = new MeshMG::MeshSmartIterator(meshInterface);
   if (dimension == 0)
-    meshInterface->get_topo_iterator(TOPO_VERTEX, *miter);
+    meshInterface->get_topo_iterator(MeshMG::TOPO_VERTEX, *miter);
   if (dimension == 1)
-    meshInterface->get_topo_iterator(TOPO_EDGE, *miter);
+    meshInterface->get_topo_iterator(MeshMG::TOPO_EDGE, *miter);
   if (dimension == 2)
-    meshInterface->get_topo_iterator(TOPO_FACE, *miter);
+    meshInterface->get_topo_iterator(MeshMG::TOPO_FACE, *miter);
   if (dimension == 3)
-    meshInterface->get_topo_iterator(TOPO_REGION, *miter);
+    meshInterface->get_topo_iterator(MeshMG::TOPO_REGION, *miter);
   meshInterface->iterator_begin(*miter);
   return reinterpret_cast<MeshIterator*>(miter);
 }
@@ -216,7 +378,7 @@ MeshIterator* MeshCAP::begin(int dimension)
  */
 MeshEntity* MeshCAP::iterate(MeshIterator* it)
 {
-  MeshSmartIterator* miter = reinterpret_cast<MeshSmartIterator*>(it);
+  auto miter = reinterpret_cast<MeshMG::MeshSmartIterator*>(it);
 
   M_MTopo topo = meshInterface->iterator_value(*miter);
 
@@ -230,26 +392,26 @@ MeshEntity* MeshCAP::iterate(MeshIterator* it)
 
 void MeshCAP::end(MeshIterator* it)
 {
-  MeshSmartIterator* miter = reinterpret_cast<MeshSmartIterator*>(it);
+  auto miter = reinterpret_cast<MeshMG::MeshSmartIterator*>(it);
   delete miter;
 }
 
 void MeshCAP::increment(MeshIterator* it)
 {
-  MeshSmartIterator* miter = reinterpret_cast<MeshSmartIterator*>(it);
+  auto miter = reinterpret_cast<MeshMG::MeshSmartIterator*>(it);
   meshInterface->iterator_next(*miter);
   it = reinterpret_cast<MeshIterator*>(miter);
 }
 
 bool MeshCAP::isDone(MeshIterator* it)
 {
-  MeshSmartIterator* miter = reinterpret_cast<MeshSmartIterator*>(it);
+  auto miter = reinterpret_cast<MeshMG::MeshSmartIterator*>(it);
   return meshInterface->iterator_end(*miter);
 }
 
 MeshEntity* MeshCAP::deref(MeshIterator* it)
 {
-  MeshSmartIterator* miter = reinterpret_cast<MeshSmartIterator*>(it);
+  auto miter = reinterpret_cast<MeshMG::MeshSmartIterator*>(it);
   M_MTopo topo = meshInterface->iterator_value(*miter);
   return toEntity(topo);
 }
@@ -259,7 +421,7 @@ void MeshCAP::getAdjacent(MeshEntity* e,
     DynamicArray<MeshEntity*>& adjacent)
 {
   M_MTopo topo = fromEntity(e);
-  MeshTopo type;
+  MeshMG::MeshTopo type;
   meshInterface->get_topo_type(topo, type);
   std::vector<M_MTopo> adjTopos;
   if (apf::getDimension(this, e) == dimension)
@@ -268,41 +430,41 @@ void MeshCAP::getAdjacent(MeshEntity* e,
     adjacent[0] = e;
     return;
   }
-  if (type == TOPO_VERTEX)
+  if (type == MeshMG::TOPO_VERTEX)
   {
     if (dimension == 1)
-      meshInterface->get_adjacency_vector(topo, TOPO_EDGE, adjTopos);
+      meshInterface->get_adjacency_vector(topo, MeshMG::TOPO_EDGE, adjTopos);
     if (dimension == 2)
-      meshInterface->get_adjacency_vector(topo, TOPO_FACE, adjTopos);
+      meshInterface->get_adjacency_vector(topo, MeshMG::TOPO_FACE, adjTopos);
     if (dimension == 3)
-      meshInterface->get_adjacency_vector(topo, TOPO_REGION, adjTopos);
+      meshInterface->get_adjacency_vector(topo, MeshMG::TOPO_REGION, adjTopos);
   }
-  if (type == TOPO_EDGE)
+  if (type == MeshMG::TOPO_EDGE)
   {
     if (dimension == 0)
-      meshInterface->get_adjacency_vector(topo, TOPO_VERTEX, adjTopos);
+      meshInterface->get_adjacency_vector(topo, MeshMG::TOPO_VERTEX, adjTopos);
     if (dimension == 2)
-      meshInterface->get_adjacency_vector(topo, TOPO_FACE, adjTopos);
+      meshInterface->get_adjacency_vector(topo, MeshMG::TOPO_FACE, adjTopos);
     if (dimension == 3)
-      meshInterface->get_adjacency_vector(topo, TOPO_REGION, adjTopos);
+      meshInterface->get_adjacency_vector(topo, MeshMG::TOPO_REGION, adjTopos);
   }
-  if (type == TOPO_FACE)
+  if (type == MeshMG::TOPO_FACE)
   {
     if (dimension == 0)
-      meshInterface->get_adjacency_vector(topo, TOPO_VERTEX, adjTopos);
+      meshInterface->get_adjacency_vector(topo, MeshMG::TOPO_VERTEX, adjTopos);
     if (dimension == 1)
-      meshInterface->get_adjacency_vector(topo, TOPO_EDGE, adjTopos);
+      meshInterface->get_adjacency_vector(topo, MeshMG::TOPO_EDGE, adjTopos);
     if (dimension == 3)
-      meshInterface->get_adjacency_vector(topo, TOPO_REGION, adjTopos);
+      meshInterface->get_adjacency_vector(topo, MeshMG::TOPO_REGION, adjTopos);
   }
-  if (type == TOPO_REGION)
+  if (type == MeshMG::TOPO_REGION)
   {
     if (dimension == 0)
-      meshInterface->get_adjacency_vector(topo, TOPO_VERTEX, adjTopos);
+      meshInterface->get_adjacency_vector(topo, MeshMG::TOPO_VERTEX, adjTopos);
     if (dimension == 1)
-      meshInterface->get_adjacency_vector(topo, TOPO_EDGE, adjTopos);
+      meshInterface->get_adjacency_vector(topo, MeshMG::TOPO_EDGE, adjTopos);
     if (dimension == 2)
-      meshInterface->get_adjacency_vector(topo, TOPO_FACE, adjTopos);
+      meshInterface->get_adjacency_vector(topo, MeshMG::TOPO_FACE, adjTopos);
   }
   adjacent.setSize(adjTopos.size());
   for (size_t i = 0; i < adjTopos.size(); i++) {
@@ -315,15 +477,15 @@ int MeshCAP::getDownward(MeshEntity* e,
     MeshEntity** down)
 {
   M_MTopo topo = fromEntity(e);
-  MeshTopo type;
+  MeshMG::MeshTopo type;
   meshInterface->get_topo_type(topo, type);
 
   int from;
-  if (type == TOPO_VERTEX)
+  if (type == MeshMG::TOPO_VERTEX)
     from = 0;
-  else if(type == TOPO_EDGE)
+  else if(type == MeshMG::TOPO_EDGE)
     from = 1;
-  else if(type == TOPO_FACE)
+  else if(type == MeshMG::TOPO_FACE)
     from = 2;
   else
     from = 3;
@@ -337,11 +499,11 @@ int MeshCAP::getDownward(MeshEntity* e,
 
   PCU_ALWAYS_ASSERT(from > dimension);
   if (dimension == 0)
-    meshInterface->get_adjacency_vector(topo, TOPO_VERTEX, adjTopos);
+    meshInterface->get_adjacency_vector(topo, MeshMG::TOPO_VERTEX, adjTopos);
   if (dimension == 1)
-    meshInterface->get_adjacency_vector(topo, TOPO_EDGE, adjTopos);
+    meshInterface->get_adjacency_vector(topo, MeshMG::TOPO_EDGE, adjTopos);
   if (dimension == 2)
-    meshInterface->get_adjacency_vector(topo, TOPO_FACE, adjTopos);
+    meshInterface->get_adjacency_vector(topo, MeshMG::TOPO_FACE, adjTopos);
 
   for (std::size_t i = 0; i < adjTopos.size(); i++)
     down[i] = toEntity(adjTopos[i]);
@@ -353,25 +515,25 @@ int MeshCAP::getDownward(MeshEntity* e,
 MeshEntity* MeshCAP::getUpward(MeshEntity* e, int i)
 {
   M_MTopo topo = fromEntity(e);
-  MeshTopo type;
+  MeshMG::MeshTopo type;
   meshInterface->get_topo_type(topo, type);
   std::vector<std::size_t> adjId;
-  if (type == TOPO_VERTEX)
+  if (type == MeshMG::TOPO_VERTEX)
   {
-    meshInterface->get_adjacency_id_vector(topo, TOPO_EDGE, adjId);
-    M_MTopo topo = meshInterface->get_topo_by_id(TOPO_EDGE, adjId[i]);
+    meshInterface->get_adjacency_id_vector(topo, MeshMG::TOPO_EDGE, adjId);
+    M_MTopo topo = meshInterface->get_topo_by_id(MeshMG::TOPO_EDGE, adjId[i]);
     return toEntity(topo);
   }
-  if (type == TOPO_EDGE)
+  if (type == MeshMG::TOPO_EDGE)
   {
-    meshInterface->get_adjacency_id_vector(topo, TOPO_FACE, adjId);
-    M_MTopo topo = meshInterface->get_topo_by_id(TOPO_FACE, adjId[i]);
+    meshInterface->get_adjacency_id_vector(topo, MeshMG::TOPO_FACE, adjId);
+    M_MTopo topo = meshInterface->get_topo_by_id(MeshMG::TOPO_FACE, adjId[i]);
     return toEntity(topo);
   }
-  if (type == TOPO_FACE)
+  if (type == MeshMG::TOPO_FACE)
   {
-    meshInterface->get_adjacency_id_vector(topo, TOPO_REGION, adjId);
-    M_MTopo topo = meshInterface->get_topo_by_id(TOPO_REGION, adjId[i]);
+    meshInterface->get_adjacency_id_vector(topo, MeshMG::TOPO_REGION, adjId);
+    M_MTopo topo = meshInterface->get_topo_by_id(MeshMG::TOPO_REGION, adjId[i]);
     return toEntity(topo);
   }
   return 0;
@@ -404,15 +566,15 @@ void MeshCAP::deleteAdjacency(int from_dim, int to_dim)
 void MeshCAP::getUp(MeshEntity* e, Up& up)
 {
   M_MTopo topo = fromEntity(e);
-  MeshTopo type;
+  MeshMG::MeshTopo type;
   meshInterface->get_topo_type(topo, type);
   std::vector<M_MTopo> adjTopos;
-  if (type == TOPO_VERTEX)
-    meshInterface->get_adjacency_vector(topo, TOPO_EDGE, adjTopos);
-  if (type == TOPO_EDGE)
-    meshInterface->get_adjacency_vector(topo, TOPO_FACE, adjTopos);
-  if (type == TOPO_FACE)
-    meshInterface->get_adjacency_vector(topo, TOPO_REGION, adjTopos);
+  if (type == MeshMG::TOPO_VERTEX)
+    meshInterface->get_adjacency_vector(topo, MeshMG::TOPO_EDGE, adjTopos);
+  if (type == MeshMG::TOPO_EDGE)
+    meshInterface->get_adjacency_vector(topo, MeshMG::TOPO_FACE, adjTopos);
+  if (type == MeshMG::TOPO_FACE)
+    meshInterface->get_adjacency_vector(topo, MeshMG::TOPO_REGION, adjTopos);
   up.n = adjTopos.size();
   for (int i = 0; i < up.n; i++) {
     up.e[i] = toEntity(adjTopos[i]);
@@ -422,15 +584,15 @@ void MeshCAP::getUp(MeshEntity* e, Up& up)
 int MeshCAP::countUpward(MeshEntity* e)
 {
   M_MTopo topo = fromEntity(e);
-  MeshTopo type;
+  MeshMG::MeshTopo type;
   meshInterface->get_topo_type(topo, type);
   std::vector<std::size_t> adjTopos;
-  if (type == TOPO_VERTEX)
-    meshInterface->get_adjacency_id_vector(topo, TOPO_EDGE, adjTopos);
-  if (type == TOPO_EDGE)
-    meshInterface->get_adjacency_id_vector(topo, TOPO_FACE, adjTopos);
-  if (type == TOPO_FACE)
-    meshInterface->get_adjacency_id_vector(topo, TOPO_REGION, adjTopos);
+  if (type == MeshMG::TOPO_VERTEX)
+    meshInterface->get_adjacency_id_vector(topo, MeshMG::TOPO_EDGE, adjTopos);
+  if (type == MeshMG::TOPO_EDGE)
+    meshInterface->get_adjacency_id_vector(topo, MeshMG::TOPO_FACE, adjTopos);
+  if (type == MeshMG::TOPO_FACE)
+    meshInterface->get_adjacency_id_vector(topo, MeshMG::TOPO_REGION, adjTopos);
   return (int)adjTopos.size();
 }
 
@@ -438,7 +600,7 @@ ModelEntity* MeshCAP::toModel(MeshEntity* e)
 {
   M_MTopo topo = fromEntity(e);
   M_GTopo gtopo;
-  GeometryTopoType gtype;
+  MeshMG::GeometryTopoType gtype;
   meshInterface->get_geom_entity(topo, gtype, gtopo);
   gmi_ent* g = toGmiEntity(gtopo);
   return reinterpret_cast<ModelEntity*>(g);
@@ -461,21 +623,21 @@ void MeshCAP::setModelEntity(MeshEntity* e, ModelEntity* me)
   apf::fail("MeshCAP::setModelEntity called!\n");
 }
 
-static GeometryTopoType getCapGeomType(int d)
+static MeshMG::GeometryTopoType getCapGeomType(int d)
 {
-  GeometryTopoType gtype = GVERTEX;
+  MeshMG::GeometryTopoType gtype = MeshMG::GVERTEX;
   switch (d) {
     case 0:
-      gtype = GVERTEX;
+      gtype = MeshMG::GVERTEX;
       break;
     case 1:
-      gtype = GEDGE;
+      gtype = MeshMG::GEDGE;
       break;
     case 2:
-      gtype = GFACE;
+      gtype = MeshMG::GFACE;
       break;
     case 3:
-      gtype = GREGION;
+      gtype = MeshMG::GREGION;
       break;
     default:
       break;
@@ -494,38 +656,38 @@ MeshEntity* MeshCAP::createVert_(ModelEntity* me)
   gmi_ent* g = reinterpret_cast<gmi_ent*>(me);
   M_GTopo gtopo = fromGmiEntity(g);
   int d = getModelType(me);
-  GeometryTopoType gtype = getCapGeomType(d);
+  MeshMG::GeometryTopoType gtype = getCapGeomType(d);
   meshInterface->create_vertex(xyz, vertex, gtype, gtopo);
   return toEntity(vertex);
 }
 
-static MeshShape getCapShape(int type)
+static MeshMG::MeshShape getCapShape(int type)
 {
-  MeshShape shape = SHAPE_UNKNOWN;
+  MeshMG::MeshShape shape = MeshMG::SHAPE_UNKNOWN;
   switch (type) {
     case Mesh::VERTEX:
-      shape = SHAPE_NODE;
+      shape = MeshMG::SHAPE_NODE;
       break;
     case Mesh::EDGE:
-      shape = SHAPE_SEGMENT;
+      shape = MeshMG::SHAPE_SEGMENT;
       break;
     case Mesh::TRIANGLE:
-      shape = SHAPE_TRIANGLE;
+      shape = MeshMG::SHAPE_TRIANGLE;
       break;
     case Mesh::QUAD:
-      shape = SHAPE_QUAD;
+      shape = MeshMG::SHAPE_QUAD;
       break;
     case Mesh::TET:
-      shape = SHAPE_TETRA;
+      shape = MeshMG::SHAPE_TETRA;
       break;
     case Mesh::HEX:
-      shape = SHAPE_HEX;
+      shape = MeshMG::SHAPE_HEX;
       break;
     case Mesh::PRISM:
-      shape = SHAPE_PRISM;
+      shape = MeshMG::SHAPE_PRISM;
       break;
     case Mesh::PYRAMID:
-      shape = SHAPE_PYRAMID;
+      shape = MeshMG::SHAPE_PYRAMID;
       break;
     default:
       break;
@@ -592,14 +754,14 @@ MeshEntity* MeshCAP::createEntity_(int type, ModelEntity* me, MeshEntity** down)
   }
 
   M_MTopo topo;
-  MeshShape shape = getCapShape(type);
+  MeshMG::MeshShape shape = getCapShape(type);
   if ( !me ) {
     meshInterface->create_topo(shape, mtopos, topo);
   }
   // if model is not 0 figure out its type
   else {
     int d = getModelType(me);
-    GeometryTopoType gtype = getCapGeomType(d);
+    MeshMG::GeometryTopoType gtype = getCapGeomType(d);
     gmi_ent* g = reinterpret_cast<gmi_ent*>(me);
     M_GTopo gtopo = fromGmiEntity(g);
     meshInterface->create_topo(shape, mtopos, topo, gtype, gtopo);
@@ -621,15 +783,12 @@ void MeshCAP::destroy_(MeshEntity* e)
 class TagCAP
 {
   public:
-    TagCAP(MeshDatabaseInterface* m,
-     const char* n,
-           int c):
+    TagCAP(MDBI* m, const char* n, int c):
       mesh(m),
       count(c),
       name(n)
     {}
-    virtual ~TagCAP()
-    {}
+    virtual ~TagCAP() {}
     virtual void* allocate() = 0;
     virtual void deallocate(void* p) = 0;
     virtual int getType() = 0;
@@ -658,7 +817,7 @@ class TagCAP
     {
       this->name = n;
     }
-    MeshDatabaseInterface* mesh;
+    MDBI* mesh;
     int count;
     std::string name;
     std::map<MeshEntity*, void*> tagContainer;
@@ -667,7 +826,7 @@ class TagCAP
 class DoubleTagCAP : public TagCAP
 {
   public:
-    DoubleTagCAP(MeshDatabaseInterface* m, const char* name, int c):
+    DoubleTagCAP(MDBI* m, const char* name, int c):
       TagCAP(m, name,c)
     {}
     virtual void* allocate()
@@ -700,7 +859,7 @@ class DoubleTagCAP : public TagCAP
 class IntTagCAP : public TagCAP
 {
   public:
-    IntTagCAP(MeshDatabaseInterface* m, const char* name, int c):
+    IntTagCAP(MDBI* m, const char* name, int c):
       TagCAP(m, name,c)
     {}
     virtual void* allocate()
@@ -925,109 +1084,129 @@ void MeshCAP::getDgCopies(MeshEntity* e, DgCopies& dgCopies, ModelEntity* me)
   apf::fail("MeshCAP::getDgCopies called!\n");
 }
 
-Mesh2* createMesh(capMesh* mesh)
-{
-  (void)mesh;
-  apf::fail("MeshCAP::createMesh called!\n");
-  return 0;
+Mesh2* createCapMesh(MDBI* mdb, GDBI* gdb, pcu::PCU *PCUObj) {
+  MeshCAP* m = new MeshCAP(mdb, gdb);
+  m->init(getLagrange(1), PCUObj);
+  return m;
 }
 
-MeshEntity* castEntity(capEntity* entity)
-{
-  (void)entity;
-  apf::fail("MeshCAP::castEntity called!\n");
-  return 0;
+Mesh2* createCapMesh(
+  gmi_model* model, const char* meshname, pcu::PCU* PCUObj
+) {
+  GDBI* gdb = gmi_export_cap(model);
+  AppContext* ctx = gdb->get_context();
+  MDBIP mdp = get_context_mesh_database_interface(ctx);
+  M_MModel mmodel;
+  MG_API_CALL(mdp.get(), get_model_by_name(meshname, mmodel));
+  MG_API_CALL(mdp.get(), set_current_model(mmodel));
+  MeshCAP* m = new MeshCAP(model, mdp);
+  m->init(getLagrange(1), PCUObj);
+  return m;
 }
 
-Mesh2* createMesh(MeshDatabaseInterface* mdb, GeometryDatabaseInterface* gdb, pcu::PCU *PCUObj)
-{
-  MeshCAP* m = new MeshCAP(mdb, gdb);
+Mesh2* createCapMesh(gmi_model* model, pcu::PCU* PCUObj) {
+  GDBI* gdb = gmi_export_cap(model);
+  AppContext* ctx = gdb->get_context();
+  MDBIP mdp = get_context_mesh_database_interface(ctx);
+  M_MModel mmodel;
+  MG_API_CALL(mdp.get(), get_model_by_index(0, mmodel));
+  MG_API_CALL(mdp.get(), set_current_model(mmodel));
+  MeshCAP* m = new MeshCAP(model, mdp);
   m->init(getLagrange(1), PCUObj);
   return m;
 }
 
-bool has_smoothCAPAnisoSizes(void) noexcept {
-#ifdef HAVE_CAPSTONE_SIZINGMETRICTOOL
-  return true;
-#else
-  return false;
-#endif
+Mesh2* makeEmptyCapMesh(
+  gmi_model* model, const char* meshname, pcu::PCU* PCUObj
+) {
+  PCU_ALWAYS_ASSERT(model);
+  GDBI* gdb = gmi_export_cap(model);
+  M_GModel gmodel;
+  MG_API_CALL(gdb, get_current_model(gmodel));
+  AppContext* ctx = gdb->get_context();
+  MDBIP mdp = get_context_mesh_database_interface(ctx);
+  M_MModel mmodel;
+  MG_API_CALL(mdp.get(), create_associated_model(mmodel, gmodel, meshname));
+  MeshCAP* m = new MeshCAP(model, mdp);
+  m->init(getLagrange(1), PCUObj);
+  return m;
 }
 
-bool smoothCAPAnisoSizes(apf::Mesh2* mesh, std::string analysis,
-  apf::Field* scales, apf::Field* frames) {
-#ifdef HAVE_CAPSTONE_SIZINGMETRICTOOL
-  // Ensure input is a MeshCAP.
-  apf::MeshCAP* m = dynamic_cast<apf::MeshCAP*>(mesh);
-  if (!m) {
-    lion_eprint(1, "ERROR: smoothCAPAnisoSizes: mesh is not an apf::MeshCAP*\n");
-    return false;
-  }
-
-  // Extract metric tensors from MeshAdapt frames and scales.
-  std::vector<Metric6> sizing6(m->count(0));
-  apf::Matrix3x3 Q;
-  apf::Vector3 H;
-  apf::MeshIterator* it = m->begin(0);
-  for (apf::MeshEntity* e = m->iterate(it); e; e = m->iterate(it)) {
-    apf::getVector(scales, e, 0, H); // Desired element lengths.
-    apf::getMatrix(frames, e, 0, Q); // MeshAdapt uses column vectors.
-    apf::Matrix3x3 L(1.0/(H[0]*H[0]), 0, 0,
-      0, 1.0/(H[1]*H[1]), 0,
-      0, 0, 1.0/(H[2]*H[2]));
-    apf::Matrix3x3 t = Q * L * apf::transpose(Q); // Invert orthogonal frames.
-    size_t id;
-    MG_API_CALL(m->getMesh(), get_id(fromEntity(e), id));
-    PCU_DEBUG_ASSERT(id != 0);
-    --id;
-    sizing6[id][0] = t[0][0];
-    sizing6[id][1] = t[0][1];
-    sizing6[id][2] = t[0][2];
-    sizing6[id][3] = t[1][1];
-    sizing6[id][4] = t[1][2];
-    sizing6[id][5] = t[2][2];
+Mesh2* generateCapMesh(
+  gmi_model* model, int dimension, pcu::PCU* PCUObj
+) {
+  PCU_ALWAYS_ASSERT(model);
+  PCU_ALWAYS_ASSERT(0 <= dimension && dimension <= 3);
+  auto gdb = gmi_export_cap(model);
+  auto ctx = gdb->get_context();
+  FunctionPtr fn;
+  switch (dimension) {
+  case 0:
+    fn = FunctionPtr(get_function(ctx, "GenerateVertexMesh"));
+    break;
+  case 1:
+    fn = FunctionPtr(get_function(ctx, "GenerateEdgeMesh"));
+    break;
+  case 2:
+    fn = FunctionPtr(get_function(ctx, "GenerateFaceMesh"));
+    break;
+  case 3:
+    fn = FunctionPtr(get_function(ctx, "GenerateRegionMesh"));
+    break;
   }
-  m->end(it);
-  auto smooth_tool = get_sizing_metric_tool(m->getMesh()->get_context(),
-    "CreateSmoothingBase");
-  if (smooth_tool == nullptr) {
-    lion_eprint(1, "ERROR: Unable to find \"CreateSmoothingBase\"\n");
-    return false;
+  M_GModel gmodel;
+  MG_API_CALL(gdb, get_current_model(gmodel));
+  set_input(fn, "Model", gmodel);
+  auto proc = get_context_processor(ctx);
+  if (proc->execute(fn) != STATUS_OK) {
+    fail("gmi_cap_gen_mesh: failed to mesh the model");
   }
-  smooth_tool->set_context(m->getMesh()->get_context());
   M_MModel mmodel;
-  MG_API_CALL(m->getMesh(), get_current_model(mmodel));
-  smooth_tool->set_metric(mmodel, "sizing6", sizing6);
-  std::vector<Metric6> ometric;
-  smooth_tool->smooth_metric(mmodel, analysis, "sizing6", ometric);
-  it = m->begin(0);
-  for (apf::MeshEntity* e = m->iterate(it); e; e = m->iterate(it)) {
-    size_t id;
-    MG_API_CALL(m->getMesh(), get_id(fromEntity(e), id));
-    PCU_DEBUG_ASSERT(id != 0);
-    --id;
-    const Metric6& m = ometric[id];
-    apf::Matrix3x3 t(m[0], m[1], m[2],
-      m[1], m[3], m[4],
-      m[2], m[4], m[5]);
-    int n = apf::eigen(t, &Q[0], &H[0]); // Eigenvectors in rows of Q.
-    PCU_DEBUG_ASSERT(n == 3);
-    Q = apf::transpose(Q); // Put eigenvectors back into columns for MeshAdapt.
-    for (int i = 0; i < 3; ++i) {
-      H[i] = 1.0/sqrt(H[i]);
-    }
-    apf::setMatrix(frames, e, 0, Q);
-    apf::setVector(scales, e, 0, H);
+  get_output(fn, "MeshModel", mmodel);
+  auto mdp = get_context_mesh_database_interface(ctx, mmodel);
+  MeshCAP* m = new MeshCAP(model, mdp);
+  m->init(getLagrange(1), PCUObj);
+  return m;
+}
+
+void disownCapModel(Mesh2* capMesh) {
+  auto m = dynamic_cast<MeshCAP*>(capMesh);
+  if (!m) fail("disownCapModel: not a Capstone mesh");
+  m->disownModel();
+}
+
+
+size_t getCapId(Mesh2* capMesh, MeshEntity* e) {
+  auto m = dynamic_cast<MeshCAP*>(capMesh);
+  if (!m) fail("getCapId: not a Capstone mesh");
+  MDBI* mdbi = m->getMesh();
+  auto mtopo = fromEntity(e);
+  size_t id;
+  MG_API_CALL(mdbi, get_id(mtopo, id));
+  return id;
+}
+
+MeshEntity* getCapEntity(Mesh2* capMesh, int dimension, size_t id) {
+  auto m = dynamic_cast<MeshCAP*>(capMesh);
+  if (!m) fail("getCapEntity: not a Capstone mesh");
+  MDBI* mdbi = m->getMesh();
+  MeshMG::MeshTopo topo = MeshMG::TOPO_UNKNOWN;
+  switch (dimension) {
+    case 0: topo = MeshMG::TOPO_VERTEX; break;
+    case 1: topo = MeshMG::TOPO_EDGE; break;
+    case 2: topo = MeshMG::TOPO_FACE; break;
+    case 3: topo = MeshMG::TOPO_REGION; break;
+    default: return nullptr;
   }
-  m->end(it);
-  return true;
-#else
-  (void) mesh;
-  (void) analysis;
-  (void) scales;
-  (void) frames;
-  apf::fail("smoothCAPAnisoSizes: Capstone does not have SizingMetricTool.");
-#endif
+  auto mtopo = mdbi->get_topo_by_id(topo, id);
+  return mtopo.is_valid() ? toEntity(mtopo) : nullptr;
+}
+
+
+CreateMG::MDBI* exportCapNative(Mesh2* capMesh) {
+  auto m = dynamic_cast<MeshCAP*>(capMesh);
+  if (!m) fail("exportCapNative: not a Capstone mesh");
+  return m->getMesh();
 }
 
-}//namespace apf
+} // namespace apf
diff --git a/apf_cap/apfCAP.h b/apf_cap/apfCAP.h
index 849c90d22..59feba225 100644
--- a/apf_cap/apfCAP.h
+++ b/apf_cap/apfCAP.h
@@ -1,187 +1,296 @@
-#ifndef APFCAP
-#define APFCAP
+/****************************************************************************** 
 
-#include <apfMesh2.h>
-#include <PCU.h>
+  Copyright 2025 Scientific Computation Research Center,
+      Rensselaer Polytechnic Institute. All rights reserved.
 
-#include "CapstoneModule.h"
-#include "CreateMG_Framework_Geometry.h"
-#include "CreateMG_Framework_Mesh.h"
-
-using namespace CreateMG;
-using namespace CreateMG::Attribution;
-using namespace CreateMG::Mesh;
-using namespace CreateMG::Geometry;
+  This work is open source software, licensed under the terms of the
+  BSD license as described in the LICENSE file in the top-level directory.
 
+******************************************************************************/
+#ifndef APF_CAP_H
+#define APF_CAP_H
+/**
+ * \file apfCAP.h
+ * \brief Capstone apf::Mesh2 implementation and interface.
+ *
+ * Like gmi_cap, the interface is used in two ways:
+ *
+ * 1. to import an existing Capstone mesh; and
+ * 2. to load or generate a mesh associated with a model which was previously
+ *    loaded by gmi_cap.
+ *
+ * \note Files which `#include` apfCAP.h should also `#include`
+ * CreateMG_Framework_Mesh.h (or another Capstone header with the full
+ * MeshDatabaseInterface definition) to import meshes or use exportCapNative.
+ */
 
+#include <string>
+#include <vector>
 
-class capEntity;
-class capMesh;
+/**
+ * \cond
+ * Forward declarations
+ */
+struct gmi_model;
+namespace pcu {
+  class PCU;
+}
+namespace CreateMG {
+  namespace Geometry { class GeometryDatabaseInterface; }
+  namespace Mesh { class MeshDatabaseInterface; }
+  typedef Geometry::GeometryDatabaseInterface GDBI;
+  typedef Mesh::MeshDatabaseInterface MDBI;
+  class Metric6;
+}
+/** \endcond */
 
 namespace apf {
 
+class Field;
+class Mesh2;
+class MeshEntity;
+
+/**
+ * \defgroup apf_cap Capstone APF mesh interface
+ *
+ * apf_cap provides access to an implementation of apf::Mesh2. The model must
+ * be loaded by gmi_cap. The interface in apfCAP.h provides additional
+ * functions to simplify loading and interacting with the underlying mesh.
+ *
+ * \{
+ */
+
+/**
+ * \brief Test for compiled Capstone library support.
+ * \return true if apf_cap was compiled with Capstone. otherwise all routines
+ *         will call apf::fail.
+ */
+bool hasCAP() noexcept;
+
+/**
+ * \brief Create an apf::Mesh2 object from a Capstone mesh database.
+ *
+ * This object should be destroyed by apf::destroyMesh. Since Capstone meshes
+ * are serial only right now, PCUObj will not be directly used, but may have an
+ * impact on collective calls involving the mesh (e.g. a PCU::Add which uses
+ * Mesh::getPCU()).
+ *
+ * \param mdb A Capstone MeshDatabaseInterface with a current M_MModel.
+ * \param gdb A Capstone GeometryDatabaseInterface with a current M_GModel.
+ * \param PCUObj The PCU communicator to define the mesh over.
+ */
+Mesh2* createCapMesh(
+  CreateMG::MDBI* mdb, CreateMG::GDBI* gdb, pcu::PCU *PCUObj
+);
+
 /**
- * \brief Creates an apf::Mesh from a CapStone mesh.
- *
- * \details This object should be destroyed by apf::destroyMesh.
- */
-Mesh2* createMesh(MeshDatabaseInterface* mdb, GeometryDatabaseInterface* gdb, pcu::PCU *PCUObj);
-
-/**
-  * \brief Casts a CapStone entity to an apf::MeshEntity.
-  *
-  * \details This does not create any objects, use freely.
-  */
-MeshEntity* castEntity(capEntity* entity);
-
-class TagCAP;
-
-class MeshCAP : public Mesh2
-{
-  public:
-    MeshCAP(MeshDatabaseInterface* mdb, GeometryDatabaseInterface* gdb);
-    virtual ~MeshCAP();
-    /* --------------------------------------------------------------------- */
-    /* Category 00: General Mesh APIs */
-    /* --------------------------------------------------------------------- */
-    // REQUIRED Member Functions //
-    int getDimension();
-    std::size_t count(int dimension);
-    Type getType(MeshEntity* e);
-    void verify();
-    // OPTIONAL Member Functions //
-    void writeNative(const char* fileName);
-    void destroyNative();
-
-    /* --------------------------------------------------------------------- */
-    /* Category 01: General Getters and Setters for vertex coordinates */
-    /* --------------------------------------------------------------------- */
-    // REQUIRED Member Functions //
-    void getPoint_(MeshEntity* e, int, Vector3& point);
-    void setPoint_(MeshEntity* e, int, Vector3 const& p);
-    void getParam(MeshEntity* e, Vector3& p);
-    void setParam(MeshEntity* e, Vector3 const& point);
-
-    /* --------------------------------------------------------------------- */
-    /* Category 02: Iterators */
-    /* --------------------------------------------------------------------- */
-    // REQUIRED Member Functions //
-    MeshIterator* begin(int dimension);
-    MeshEntity* iterate(MeshIterator* it);
-    void end(MeshIterator* it);
-    void increment(MeshIterator* it);
-    bool isDone(MeshIterator* it);
-    MeshEntity* deref(MeshIterator* it);
-
-    /* --------------------------------------------------------------------- */
-    /* Category 03: Adjacencies */
-    /* --------------------------------------------------------------------- */
-    // REQUIRED Member Functions //
-    void getAdjacent(MeshEntity* e, int dimension, Adjacent& adjacent);
-    int getDownward(MeshEntity* e, int dimension, MeshEntity** adjacent);
-    MeshEntity* getUpward(MeshEntity* e, int i);
-    bool hasUp(MeshEntity* e);
-    // OPTIONAL Member Functions //
-    bool hasAdjacency(int from_dim, int to_dim);
-    void createAdjacency(int from_dim, int to_dim);
-    void deleteAdjacency(int from_dim, int to_dim);
-    void getUp(MeshEntity* e, Up& up);
-    int countUpward(MeshEntity* e);
-
-    /* --------------------------------------------------------------------- */
-    /* Category 04: CAD model inquires */
-    /* --------------------------------------------------------------------- */
-    // REQUIRED Member Functions //
-    ModelEntity* toModel(MeshEntity* e);
-    // OPTIONAL Member Functions //
-    gmi_model* getModel();
-    void setModel(gmi_model* newModel);
-    void setModelEntity(MeshEntity* e, ModelEntity* me);
-
-    /* --------------------------------------------------------------------- */
-    /* Category 05: Entity Creation/Deletion */
-    /* --------------------------------------------------------------------- */
-    // REQUIRED Member Functions //
-    MeshEntity* createVert_(ModelEntity* me);
-    MeshEntity* createEntity_(int type, ModelEntity* me, MeshEntity** down);
-    void destroy_(MeshEntity* e);
-
-    /* --------------------------------------------------------------------- */
-    /* Category 06: Attachable Data Functionality */
-    /* --------------------------------------------------------------------- */
-    // REQUIRED Member Functions //
-    MeshTag* createDoubleTag(const char* name, int size);
-    MeshTag* createIntTag(const char* name, int size);
-    MeshTag* createLongTag(const char* name, int size);
-    MeshTag* findTag(const char* name);
-    void destroyTag(MeshTag* t);
-    void renameTag(MeshTag* t, const char* name);
-    void getTags(DynamicArray<MeshTag*>& tags);
-    /* void getTag(MeshEntity* e, MeshTag* t, void* data); */
-    /* void setTag(MeshEntity* e, MeshTag* t, void const* data); */
-    void getDoubleTag(MeshEntity* e, MeshTag* tag, double* data);
-    void setDoubleTag(MeshEntity* e, MeshTag* tag, double const* data);
-    void getIntTag(MeshEntity* e, MeshTag* tag, int* data);
-    void setIntTag(MeshEntity* e, MeshTag* tag, int const* data);
-    void getLongTag(MeshEntity* e, MeshTag* tag, long* data);
-    void setLongTag(MeshEntity* e, MeshTag* tag, long const* data);
-    void removeTag(MeshEntity* e, MeshTag* t);
-    bool hasTag(MeshEntity* e, MeshTag* t);
-    int getTagType(MeshTag* t);
-    int getTagSize(MeshTag* t);
-    const char* getTagName(MeshTag* t);
-    unsigned getTagChecksum(MeshTag*,int);
-
-
-    /* --------------------------------------------------------------------- */
-    /* Category 07: Distributed Meshes */
-    /* --------------------------------------------------------------------- */
-    // REQUIRED Member Functions //
-    bool isShared(MeshEntity* e);
-    bool isGhost(MeshEntity*) { return false; }
-    bool isGhosted(MeshEntity*) { return false; }
-    bool isOwned(MeshEntity* e);
-    int getOwner(MeshEntity* e);
-    void getRemotes(MeshEntity* e, Copies& remotes);
-    void getResidence(MeshEntity* e, Parts& residence);
-    int getId();
-    void setResidence(MeshEntity*, Parts&) {}
-    void acceptChanges() {}
-    // OPTIONAL Member Functions //
-    void deleteGhost(MeshEntity*) {}
-    void addGhost(MeshEntity*, int, MeshEntity*) {}
-    int getGhosts(MeshEntity*, Copies&) { return 0; }
-    void migrate(Migration* plan);
-    void setRemotes(MeshEntity*, Copies&) {}
-    void addRemote(MeshEntity*, int, MeshEntity*) {}
-    void clearRemotes(MeshEntity*) {}
-
-
-    /* --------------------------------------------------------------------- */
-    /* Category 08: Periodic Meshes */
-    /* --------------------------------------------------------------------- */
-    // REQUIRED Member Functions //
-    bool hasMatching() { return false; }
-    void getMatches(MeshEntity* e, Matches& m);
-    // OPTIONAL Member Functions //
-    void addMatch(MeshEntity*, int, MeshEntity*) {}
-    void clearMatches(MeshEntity*) {}
-    void clear_() {}
-    void getDgCopies(MeshEntity* e, DgCopies& dgCopies, ModelEntity* me);
-
-
-
-    MeshDatabaseInterface* getMesh() { return meshInterface; }
-  protected:
-    /* CapstoneModule capModule; */
-    MeshDatabaseInterface* meshInterface;
-    GeometryDatabaseInterface  *geomInterface;
-    /* AppContext                 *c; */
-    int iterDim;
-    int d;
-    gmi_model* model;
-    std::vector<TagCAP*> tags;
-};
+ * \brief Create an apf::Mesh2 object from the first mesh linked to a loaded
+ * geometry model.
+ *
+ * The gmi_model should be loaded previously by gmi_load (on a .cre file),
+ * gmi_cap_load, or gmi_cap_load_selective. Try to load the first mesh (by
+ * index).
+ *
+ * \param model A gmi_model associated with a Capstone geometry.
+ * \param PCUObj The PCU communicator to define the mesh over
+ * \return an apf::Mesh2 interface to the Capstone mesh.
+ */
+Mesh2* createCapMesh(gmi_model* model, pcu::PCU* PCUObj);
+
+/**
+ * \brief Create an apf::Mesh2 object from a mesh linked to a loaded
+ * geometry model.
+ *
+ * The gmi_model should be loaded previously by gmi_load (on a .cre file),
+ * gmi_cap_load, or gmi_cap_load_selective. The list of acceptable values for
+ * meshname can be found by using gmi_cap_probe or direct Capstone interfaces.
+ *
+ * \param model A gmi_model associated with a Capstone geometry.
+ * \param meshname The name of a mesh associated with
+ * \param PCUObj The PCU communicator to define the mesh over
+ * \return an apf::Mesh2 interface to the Capstone mesh.
+ */
+Mesh2* createCapMesh(gmi_model* model, const char* meshname, pcu::PCU* PCUObj);
+
+/**
+ * \brief Generate Capstone mesh object on a model linked to a Capstone model.
+ *
+ * \param model A Capstone GMI model to generate a mesh on.
+ * \param dimension The dimension of mesh to generate
+ * \param PCUObj The PCU object to link to the mesh.
+ * \return An apf::Mesh2 object with the mesh.
+ */
+Mesh2* generateCapMesh(
+  gmi_model* model, int dimension, pcu::PCU* PCUObj
+);
+
+/**
+ * \brief Make an empty Capstone M_MModel associated with the model.
+ *
+ * \param model Previously loaded Capstone gmi_model.
+ * \param meshname Name for new mesh in the Capstone database.
+ * \param PCUObj The PCU object to associate the new apf::Mesh2 with.
+ * \return a new apf::Mesh2 object.
+ */
+Mesh2* makeEmptyCapMesh(
+  gmi_model* model, const char* meshname, pcu::PCU* PCUObj
+);
+
+/**
+ * \brief Disown capMesh's gmi_model.
+ *
+ * Mark the gmi_model as non-owned so that the destructor does not call
+ * gmi_destroy.
+ *
+ * \param capMesh A Capstone mesh wrapper.
+ */
+void disownCapModel(Mesh2* capMesh);
+
+/**
+ * \brief Get the native Capstone id of an APF entity.
+ *
+ * \param m A Capstone mesh
+ * \param e A MeshEntity on m
+ * \return Unique id associated with the underlying Capstone Topo.
+ */
+size_t getCapId(Mesh2* m, MeshEntity* e);
+
+/**
+ * \brief Get an entity from a Capstone mesh by dimension and native id.
+ *
+ * \param m A Capstone mesh
+ * \param dimension The dimension of the entity to retrieve
+ * \param id The native Capstone id
+ * \return The corresponding MeshEntity or nullptr if no such entity exists
+ */
+MeshEntity* getCapEntity(Mesh2* m, int dimension, size_t id);
+
+/**
+ * \brief Get native Capstone mesh database interface.
+ * \param capMesh Previously loaded capstone mesh.
+ * \return Underlying Capstone mesh database interface from capMesh.
+ */
+CreateMG::MDBI* exportCapNative(Mesh2* capMesh);
+
+/**
+ * \defgroup apf_cap_sizing Capstone APF mesh sizing utilities
+ * \{
+ */
+
+/**
+ * \brief Load Capstone bulk sizing into apf::Fields.
+ *
+ * \param[in] m An apf Capstone mesh
+ * \param[in] sizing A bulk sizing vector of Metric6 from Capstone
+ * \param[out] scales An apf::Field of apf::Vector3's on vertices for
+ *                    anisotropic sizing scales
+ * \param[out] frames An apf::Field of apf::Matrix3x3's on vertices for
+ *                    anisotropic sizing frames
+ * \return true on success, false otherwise
+ */
+bool loadCapSizing(
+  apf::Mesh2* m, const std::vector<CreateMG::Metric6>& sizing,
+  apf::Field* scales, apf::Field* frames
+);
+
+/**
+ * \brief Load Capstone bulk sizing into apf::Fields.
+ *
+ * This overload is a convenience wrapper to simplify Field creation.
+ *
+ * \param[in] m An apf Capstone mesh
+ * \param[in] sizing A bulk sizing vector of Metric6 from Capstone
+ * \param[out] scales The name of a new apf::Field to create with anisotropic
+ *                    sizing scales
+ * \param[out] frames The name of a new apf::Field to create with anisotropic
+ *                    sizing frames
+ * \return true on success, false otherwise
+ */
+bool loadCapSizing(
+  apf::Mesh2* m, const std::vector<CreateMG::Metric6>& sizing,
+  const char* scales, const char* frames
+);
+
+/**
+ * \brief Load metrics from a sizing file and vmap file.
+ *
+ * \param m An apf Capstone mesh
+ * \param sizingFile The name of the bulk sizing file
+ * \param vmapFile The name of the vmap file
+ * \param sizing6 A vector to fill with bulk sizing metrics
+ */
+void loadCapSizingFileMetrics(
+  apf::Mesh2* m, const std::string& sizingFile, const std::string& vmapFile,
+  std::vector<CreateMG::Metric6>& sizing6
+);
+
+/**
+ * \brief Load Capstone bulk sizing into apf::Fields from a sizing file.
+ *
+ * The bulk sizing file contains 6 doubles for each vertex. The vmap file
+ * contains 1 size_t for each vertex, plus an extra one with vertex count
+ * information. The vmap contains the info to convert solver ids to mesh ids.
+ *
+ * \param[in] m An apf Capstone mesh
+ * \param[in] sizingFile The name of the bulk sizing file.
+ * \param[in] vmapFile The name of the vmap file.
+ * \param[out] scales An apf::Field of apf::Vector3's on vertices with
+ *                    anisotropic sizing scales
+ * \param[out] frames An apf::Field of apf::Matrix3x3's on vertices with
+ *                    anisotropic sizing frames
+ * \param[in] smooth A boolean to request smoothing before loading to
+ *                   apf::Fields
+ * \param[in] analysis The analysis which may contain additional sizing
+ *                     information (passed to Capstone smoothing routine)
+ * \return true on success, false on failure or if smoothing is requested but
+ *         not supported
+ */
+bool loadCapSizingFile(
+  apf::Mesh2* m, const std::string& sizingFile, const std::string& vmapFile,
+  apf::Field* scales, apf::Field* frames,
+  bool smooth = false, const std::string& analysis = ""
+);
+
+/**
+ * \brief Load Capstone bulk sizing into apf::Fields from a sizing file.
+ *
+ * The bulk sizing file contains 6 doubles for each vertex. The vmap file
+ * contains 1 size_t for each vertex, plus an extra one with vertex count
+ * information. The vmap contains the info to convert solver ids to mesh ids.
+ *
+ * This overload is a convenience wrapper to simplify Field creation.
+ *
+ * \param[in] m An apf Capstone mesh
+ * \param[in] sizingFile The name of the bulk sizing file.
+ * \param[in] vmapFile The name of the vmap file.
+ * \param[out] scales The name of a new apf::Field to create with anisotropic
+ *                    sizing scales
+ * \param[out] frames The name of a new apf::Field to create with anisotropic
+ *                    sizing frames
+ * \param[in] smooth Whether to request smoothing before loading to apf::Fields
+ * \param[in] analysis The analysis which may contain additional sizing
+ *                     information (passed to Capstone smoothing routine)
+ * \return true on success, false on failure or if smoothing is requested but
+ *         not supported
+ */
+bool loadCapSizingFile(
+  apf::Mesh2* m, const std::string& sizingFile, const std::string& vmapFile,
+  const char* scales, const char* frames,
+  bool smooth = false, const std::string& analysis = ""
+);
+
+/**
+ * \brief Extract Capstone metric tensors from MeshAdapt frames and scales.
+ *
+ * \param[in] m An apf Capstone mesh
+ * \param[in] scales An apf::Field with anisotropic sizing scales
+ * \param[in] frames An apf::Field with anisotropic sizing frames
+ * \param[out] sizing A vector to output metric tensors to
+ */
+void extractCapSizing(
+  apf::Mesh2* m, apf::Field* scales, apf::Field* frames,
+  std::vector<CreateMG::Metric6>& sizing
+);
 
 /**
  * \brief Test for smoothCAPAnisoSizes support.
@@ -194,7 +303,7 @@ class MeshCAP : public Mesh2
  *          against. Otherwise the call will always apf::fail. Use this
  *          function to programmatically test for the capability.
  */
-bool has_smoothCAPAnisoSizes(void) noexcept;
+bool has_smoothCapAnisoSizes(void) noexcept;
 
 /**
  * \brief Use the SizingMetricTool to smooth a size field on a Capstone mesh.
@@ -204,11 +313,16 @@ bool has_smoothCAPAnisoSizes(void) noexcept;
  * \param frames An apf::Field of apf::Matrix3x3 with orthogonal basis frames.
  * \param scales An apf::Field of apf::Vector3 with frame scales (eigenvalues).
  * \return A boolean indicating success.
- * \pre m must be an apf::MeshCAP.
+ * \pre m must be a Capstone mesh.
  */
-bool smoothCAPAnisoSizes(apf::Mesh2* m, std::string analysis,
-  apf::Field* scales, apf::Field* frames);
+bool smoothCapAnisoSizes(
+  apf::Mesh2* m, std::string analysis, apf::Field* scales, apf::Field* frames
+);
+
+/** \} apf_cap_sizing */
+
+/** \} apf_cap */
 
-}//namespace apf
+} // namespace apf
 
-#endif
+#endif // APF_CAP_H
diff --git a/apf_cap/apfCAPsizing.cc b/apf_cap/apfCAPsizing.cc
new file mode 100644
index 000000000..38b08e5e7
--- /dev/null
+++ b/apf_cap/apfCAPsizing.cc
@@ -0,0 +1,223 @@
+#include "apfCAP.h"
+
+#include <exception>
+#include <fstream>
+#include <stdexcept>
+
+#include <apf.h>
+#include <apfMesh2.h>
+#include <lionPrint.h>
+#include <pcu_util.h>
+
+#include <CreateMG_Framework_Mesh.h>
+#ifdef PUMI_HAS_CAPSTONE_SIZINGMETRICTOOL
+#include <CreateMG_SizingMetricTool.h>
+#endif
+
+namespace {
+
+/** \brief Load bulk file of one type into std::vector. */
+template <class T>
+typename std::enable_if<std::is_trivially_copyable<T>::value, void>::type
+load_file(const std::string& fname, std::vector<T>& data) {
+  try {
+    std::ifstream file(fname, std::ios::in | std::ios::binary);
+    if (!file.is_open()) {
+      throw std::runtime_error("failed to open file `" + fname + "`");
+    }
+    T val;
+    while (file) {
+      file.read(reinterpret_cast<char*>(std::addressof(val)), sizeof(T));
+      if (file.gcount() == sizeof(T)) data.push_back(val);
+    }
+  } catch (...) {
+    std::throw_with_nested(
+      std::runtime_error("failed to load file `" + fname + "`")
+    );
+  }
+}
+
+bool doSmoothing(
+  apf::Mesh2* m, const std::string& analysis,
+  const std::vector<CreateMG::Metric6>& in, std::vector<CreateMG::Metric6>& out
+) {
+  auto smooth_tool = CreateMG::get_sizing_metric_tool(
+    apf::exportCapNative(m)->get_context(), "CreateSmoothingBase"
+  );
+  if (smooth_tool == nullptr) {
+    lion_eprint(1, "ERROR: Unable to find \"CreateSmoothingBase\"\n");
+    return false;
+  }
+  smooth_tool->set_context(apf::exportCapNative(m)->get_context());
+  CreateMG::M_MModel mmodel;
+  MG_API_CALL(apf::exportCapNative(m), get_current_model(mmodel));
+  out.resize(in.size());
+  smooth_tool->set_metric(mmodel, "apf_cap_sizing6", in);
+  smooth_tool->smooth_metric(mmodel, analysis, "apf_cap_sizing6", out);
+  smooth_tool->remove(mmodel, "apf_cap_sizing6");
+  return true;
+}
+
+} // namespace
+
+namespace apf {
+
+bool loadCapSizing(
+  apf::Mesh2* m, const std::vector<CreateMG::Metric6>& sizing,
+  apf::Field* scales, apf::Field* frames
+) {
+  PCU_DEBUG_ASSERT(m);
+  PCU_DEBUG_ASSERT(scales);
+  PCU_DEBUG_ASSERT(frames);
+  apf::Vector3   scale;
+  apf::Matrix3x3 frame;
+  for (std::size_t i = 0; i < sizing.size(); ++i) {
+    // field eigen values
+    apf::Matrix3x3 mat(
+      sizing[i][0], sizing[i][1], sizing[i][2],
+      sizing[i][1], sizing[i][3], sizing[i][4],
+      sizing[i][2], sizing[i][4], sizing[i][5]
+    );
+    int n = apf::eigen(mat, &frame[0], &scale[0]);
+    PCU_DEBUG_ASSERT(n == 3);
+    for (int i = 0; i < 3; i++) {
+      scale[i] = 1.0/sqrt(scale[i]);
+    }
+    // MeshAdapt wants frames on the columns.
+    frame = apf::transpose(frame);
+    apf::MeshEntity* e = getCapEntity(m, 0, i + 1);
+    apf::setVector(scales, e, 0, scale);
+    apf::setMatrix(frames, e, 0, frame);
+  }
+  return true;
+}
+
+bool loadCapSizing(
+  apf::Mesh2* m, const std::vector<CreateMG::Metric6>& sizing,
+  const char* scales, const char* frames
+) {
+  apf::Field* scalesField = apf::createFieldOn(m, scales, VECTOR);
+  apf::Field* framesField = apf::createFieldOn(m, frames, MATRIX);
+  return loadCapSizing(m, sizing, scalesField, framesField);
+}
+
+void loadCapSizingFileMetrics(
+  apf::Mesh2* m, const std::string& sizingFile, const std::string& vmapFile,
+  std::vector<CreateMG::Metric6>& sizing6
+) {
+  std::vector<size_t> vmap;
+  std::vector<double> sizing;
+  load_file(vmapFile, vmap);
+  load_file(sizingFile, sizing);
+  // Validate sizing
+  if (sizing.size() / 6 != m->count(0)) {
+    std::string msg = "loadCapSizingFile: sizing file size does not match mesh"
+      " vertex count: " + std::to_string(sizing.size() / 6) + " / " +
+      std::to_string(m->count(0));
+    fail(msg.c_str());
+  }
+  // Validate vmap
+  if (
+    vmap.size() < 1 || vmap[0] != vmap.size() - 1 || vmap[0] != m->count(0)
+  ) {
+    fail("loadCapSizingFile: loaded invalid vmap");
+  }
+  // Copy sizefield to Metric6
+  sizing6.resize(vmap[0]);
+  for (size_t si = 0, vi = 1; vi < vmap.size(); ++vi) {
+    size_t tid = vmap[vi];
+    PCU_DEBUG_ASSERT(tid != 0);
+    for (size_t j = 0; j < 6; ++j, ++si) {
+      sizing6[tid - 1][j] = sizing[si];
+    }
+  }
+}
+
+bool loadCapSizingFile(
+  apf::Mesh2* m, const std::string& sizingFile, const std::string& vmapFile,
+  apf::Field* scales, apf::Field* frames,
+  bool smooth, const std::string& analysis
+) {
+  if (smooth && !has_smoothCapAnisoSizes()) {
+    lion_eprint(1,
+      "WARNING: loadCapSizingFile: smoothing requested but apf_cap was"
+      " compiled without support for smoothing\n");
+    return false;
+  }
+  std::vector<CreateMG::Metric6> sizing6;
+  loadCapSizingFileMetrics(m, sizingFile, vmapFile, sizing6);
+  if (smooth) {
+    std::vector<CreateMG::Metric6> ometric;
+    if (!doSmoothing(m, analysis, sizing6, ometric)) return false;
+    std::swap(sizing6, ometric);
+  }
+  return loadCapSizing(m, sizing6, scales, frames);
+}
+
+bool loadCapSizingFile(
+  apf::Mesh2* m, const std::string& sizingFile, const std::string& vmapFile,
+  const char* scales, const char* frames,
+  bool smooth, const std::string& analysis
+) {
+  apf::Field* scalesField = apf::createFieldOn(m, scales, VECTOR);
+  apf::Field* framesField = apf::createFieldOn(m, frames, MATRIX);
+  return loadCapSizingFile(
+    m, sizingFile, vmapFile, scalesField, framesField, smooth, analysis
+  );
+}
+
+void extractCapSizing(
+  apf::Mesh2* mesh, apf::Field* scales, apf::Field* frames,
+  std::vector<CreateMG::Metric6>& sizing
+) {
+  apf::Matrix3x3 Q;
+  apf::Vector3 H;
+  sizing.resize(mesh->count(0));
+  apf::MeshIterator* it = mesh->begin(0);
+  for (apf::MeshEntity* e = mesh->iterate(it); e; e = mesh->iterate(it)) {
+    apf::getVector(scales, e, 0, H); // Desired element lengths.
+    apf::getMatrix(frames, e, 0, Q); // MeshAdapt uses column vectors.
+    apf::Matrix3x3 L(
+      1.0/(H[0]*H[0]), 0, 0,
+      0, 1.0/(H[1]*H[1]), 0,
+      0, 0, 1.0/(H[2]*H[2])
+    );
+    apf::Matrix3x3 t = Q * L * apf::transpose(Q); // Invert orthogonal frames.
+    size_t id = getCapId(mesh, e);
+    PCU_DEBUG_ASSERT(id != 0);
+    --id;
+    sizing[id][0] = t[0][0];
+    sizing[id][1] = t[0][1];
+    sizing[id][2] = t[0][2];
+    sizing[id][3] = t[1][1];
+    sizing[id][4] = t[1][2];
+    sizing[id][5] = t[2][2];
+  }
+  mesh->end(it);
+}
+
+bool has_smoothCapAnisoSizes(void) noexcept {
+#ifdef PUMI_HAS_CAPSTONE_SIZINGMETRICTOOL
+  return true;
+#else
+  return false;
+#endif
+}
+
+bool smoothCapAnisoSizes(apf::Mesh2* mesh, std::string analysis,
+  apf::Field* scales, apf::Field* frames) {
+#ifdef PUMI_HAS_CAPSTONE_SIZINGMETRICTOOL
+  std::vector<CreateMG::Metric6> sizing6, ometric;
+  extractCapSizing(mesh, scales, frames, sizing6);
+  if (!doSmoothing(mesh, analysis, sizing6, ometric)) return false;
+  return loadCapSizing(mesh, ometric, frames, scales);
+#else
+  (void) mesh;
+  (void) analysis;
+  (void) scales;
+  (void) frames;
+  apf::fail("smoothCAPAnisoSizes: Capstone does not have SizingMetricTool.");
+#endif
+}
+
+} // namespace apf
diff --git a/gmi_cap/CMakeLists.txt b/gmi_cap/CMakeLists.txt
index 3fa431f34..8df250997 100644
--- a/gmi_cap/CMakeLists.txt
+++ b/gmi_cap/CMakeLists.txt
@@ -3,7 +3,7 @@ if(DEFINED TRIBITS_PACKAGE)
   return()
 endif()
 
-if(NOT ENABLE_CAPSTONE)
+if(NOT PUMI_ENABLE_CAPSTONE)
   return()
 endif()
 
@@ -13,7 +13,8 @@ set(HEADERS gmi_cap.h)
 
 add_library(gmi_cap ${SOURCES})
 
-target_link_libraries(gmi_cap PUBLIC gmi pcu lion capstone_module framework_testing)
+target_link_libraries(gmi_cap PUBLIC gmi framework_core)
+target_link_libraries(gmi_cap PRIVATE pcu lion capstone_module)
 
 # Include directories
 target_include_directories(gmi_cap PUBLIC
diff --git a/gmi_cap/gmi_cap.cc b/gmi_cap/gmi_cap.cc
index 64b23e26a..9996caed6 100644
--- a/gmi_cap/gmi_cap.cc
+++ b/gmi_cap/gmi_cap.cc
@@ -8,12 +8,44 @@
 
 *******************************************************************************/
 #include "gmi_cap.h"
+
+#include <memory>
 #include <stdlib.h>
-#include <gmi.h>
 #include <vector>
-#include <pcu_util.h>
+
+#include <gmi.h>
 #include <lionPrint.h>
+#include <pcu_util.h>
+
+#include <CapstoneModule.h>
+#include <CreateMG_AppProcessor.h>
+#include <CreateMG_Framework_Geometry.h>
+#include <CreateMG_Function.h>
+#include <CreateMG_Reader.h>
+
+using namespace CreateMG;
+
+static std::unique_ptr<CapstoneModule> cs_module;
+
+void gmi_cap_start(void) {
+  if (!cs_module) {
+    cs_module.reset(new CapstoneModule(
+      "SCOREC/gmi_cap",
+      "Geometry Database : SMLIB",
+      "Mesh Database : Create",
+      "Attribution Database : Create"
+    ));
+    PCU_ALWAYS_ASSERT(cs_module->get_context());
+    PCU_ALWAYS_ASSERT(cs_module->get_geometry());
+    PCU_ALWAYS_ASSERT(cs_module->get_mesh());
+  }
+}
 
+void gmi_cap_stop(void) {
+  if (!cs_module)
+    gmi_fail("gmi_cap_stop called before gmi_cap_start");
+  cs_module.reset();
+}
 
 gmi_ent* toGmiEntity(M_GTopo topo)
 {
@@ -30,10 +62,10 @@ M_GTopo fromGmiEntity(gmi_ent* g)
   return topo;
 }
 
-
 struct cap_model {
   struct gmi_model model;
-  GeometryDatabaseInterface* geomInterface;
+  GDBI* geomInterface;
+  M_GModel gmodel;
   bool owned;
 };
 
@@ -42,17 +74,17 @@ static gmi_iter* begin(gmi_model* m, int dim)
   cap_model* cm = (cap_model*)m;
   M_GBRep brep;
   cm->geomInterface->get_brep_by_index(0, brep);
-  GeometrySmartIterator* giter = new GeometrySmartIterator(cm->geomInterface);
+  auto giter = new Geometry::GeometrySmartIterator(cm->geomInterface);
   if (dim == 0)
-    cm->geomInterface->get_topo_iterator(brep, VERTEX, *giter);
+    cm->geomInterface->get_topo_iterator(brep, Geometry::VERTEX, *giter);
   if (dim == 1)
-    cm->geomInterface->get_topo_iterator(brep, EDGE, *giter);
+    cm->geomInterface->get_topo_iterator(brep, Geometry::EDGE, *giter);
   if (dim == 2)
-    cm->geomInterface->get_topo_iterator(brep, FACE, *giter);
+    cm->geomInterface->get_topo_iterator(brep, Geometry::FACE, *giter);
   if (dim == 3)
-    cm->geomInterface->get_topo_iterator(brep, REGION, *giter);
+    cm->geomInterface->get_topo_iterator(brep, Geometry::REGION, *giter);
   cm->geomInterface->iterator_begin(*giter);
-  return (gmi_iter*)(giter);
+  return reinterpret_cast<gmi_iter*>(giter);
 }
 
 /* NOTE: giter is located at the first item in the list, therefore
@@ -61,7 +93,7 @@ static gmi_iter* begin(gmi_model* m, int dim)
 static gmi_ent* next(gmi_model*m, gmi_iter* i)
 {
   cap_model* cm = (cap_model*)m;
-  GeometrySmartIterator* giter = (GeometrySmartIterator*)i;
+  auto giter = reinterpret_cast<Geometry::GeometrySmartIterator*>(i);
 
   M_GTopo topo = cm->geomInterface->iterator_value(*giter);
 
@@ -73,10 +105,8 @@ static gmi_ent* next(gmi_model*m, gmi_iter* i)
   return toGmiEntity(topo);
 }
 
-static void end(gmi_model*, gmi_iter* i)
-{
-  GeometrySmartIterator* giter = (GeometrySmartIterator*)i;
-  delete giter;
+static void end(gmi_model*, gmi_iter* i) {
+  delete reinterpret_cast<Geometry::GeometrySmartIterator*>(i);
 }
 
 static int get_dim(gmi_model* m, gmi_ent* e)
@@ -108,13 +138,13 @@ static gmi_ent* find(gmi_model* m, int dim, int tag)
   cap_model* cm = (cap_model*)m;
   M_GTopo topo;
   if (dim == 0)
-    topo = cm->geomInterface->get_topo_by_id(VERTEX, tag);
+    topo = cm->geomInterface->get_topo_by_id(Geometry::VERTEX, tag);
   else if (dim == 1)
-    topo = cm->geomInterface->get_topo_by_id(EDGE, tag);
+    topo = cm->geomInterface->get_topo_by_id(Geometry::EDGE, tag);
   else if (dim == 2)
-    topo = cm->geomInterface->get_topo_by_id(FACE, tag);
+    topo = cm->geomInterface->get_topo_by_id(Geometry::FACE, tag);
   else if (dim == 3)
-    topo = cm->geomInterface->get_topo_by_id(REGION, tag);
+    topo = cm->geomInterface->get_topo_by_id(Geometry::REGION, tag);
   else
     gmi_fail("input dim is out of range.");
   return toGmiEntity(topo);
@@ -136,21 +166,21 @@ static gmi_set* adjacent(gmi_model* m, gmi_ent* e, int dim)
   int edim = gmi_dim(m, e);
   std::vector<M_GTopo> gtopos;
   if (edim == 0 && dim == 1)
-    cm->geomInterface->get_adjacency(gtopo, EDGE, gtopos);
+    cm->geomInterface->get_adjacency(gtopo, Geometry::EDGE, gtopos);
   else if (edim == 1 && dim == 0)
-    cm->geomInterface->get_adjacency(gtopo, VERTEX, gtopos);
+    cm->geomInterface->get_adjacency(gtopo, Geometry::VERTEX, gtopos);
   else if (edim == 1 && dim == 2)
-    cm->geomInterface->get_adjacency(gtopo, FACE, gtopos);
+    cm->geomInterface->get_adjacency(gtopo, Geometry::FACE, gtopos);
   else if (edim == 2 && dim == 0)
-    cm->geomInterface->get_adjacency(gtopo, VERTEX, gtopos);
+    cm->geomInterface->get_adjacency(gtopo, Geometry::VERTEX, gtopos);
   else if (edim == 2 && dim == 1)
-    cm->geomInterface->get_adjacency(gtopo, EDGE, gtopos);
+    cm->geomInterface->get_adjacency(gtopo, Geometry::EDGE, gtopos);
   else if (edim == 2 && dim == 3)
-    cm->geomInterface->get_adjacency(gtopo, REGION, gtopos);
+    cm->geomInterface->get_adjacency(gtopo, Geometry::REGION, gtopos);
   else if (edim == 1 && dim == 3)
-    cm->geomInterface->get_adjacency(gtopo, REGION, gtopos);
+    cm->geomInterface->get_adjacency(gtopo, Geometry::REGION, gtopos);
   else if (edim == 3 && dim == 2)
-    cm->geomInterface->get_adjacency(gtopo, FACE, gtopos);
+    cm->geomInterface->get_adjacency(gtopo, Geometry::FACE, gtopos);
   else
     gmi_fail("requested adjacency not available\n");
   return vector_to_set(gtopos);
@@ -285,7 +315,7 @@ static int is_in_closure_of(struct gmi_model* m, struct gmi_ent* e,
     cap_model* cm = (cap_model*)m;
     M_GTopo ce = fromGmiEntity(e);
     M_GTopo cet = fromGmiEntity(et);
-    GeometryTopo ce_type;
+    Geometry::GeometryTopo ce_type;
     MG_API_CALL(cm->geomInterface, get_topo_type(ce, ce_type));
     std::vector<M_GTopo> adj;
     MG_API_CALL(cm->geomInterface, get_adjacency(cet, ce_type, adj));
@@ -325,28 +355,78 @@ static int is_discrete_ent(struct gmi_model*, struct gmi_ent* e)
   return 0;
 }
 
-static void destroy(gmi_model* m)
-{
+static void destroy(gmi_model* m) {
   cap_model* cm = (cap_model*)m;
+  if (cm->owned) {
+    if (!cs_module) gmi_fail(
+      "gmi_cap destroy: called before gmi_cap_start or after gmi_cap_stop"
+    );
+    auto ctx = cm->geomInterface->get_context();
+    FunctionPtr fn = get_function(ctx, "DeleteGeometryModel");
+    set_input(fn, "Model", cm->gmodel);
+    auto proc = get_context_processor(ctx);
+    if (proc->execute(fn) != STATUS_OK)
+      gmi_fail("gmi_cap destroy: failed to delete Capstone geometry model");
+  }
   free(cm);
 }
 
-/* static gmi_model* create_cre(const char* filename) */
-/* { */
-/*   return gmi_cap_load(filename); */
-/* } */
-
-static struct gmi_model_ops ops;
-
+static gmi_model* create_cre(const char* filename) {
+  return gmi_cap_load(filename);
+}
 
-void gmi_cap_start(void)
-{
+void gmi_cap_probe(
+  const char* creFileName, std::string& model_content,
+  std::vector<std::string>& mesh_names, std::vector<std::string>& mesh_contents
+) {
+  model_content.clear();
+  mesh_names.clear();
+  mesh_contents.clear();
+	Reader *reader = get_reader(
+    cs_module->get_context(), "Create Native Reader"
+  );
+  DataFileInfo info;
+  if (!reader->probe(creFileName, info))
+    gmi_fail("gmi_cap_probe: failed to read file");
+  for (std::size_t i = 0; i < info.get_num_sections(); i++) {
+    std::string secType, secName;
+    info.get_section(i, secType, secName);
+    if (secType == "mmodel") mesh_names.push_back(secName);
+    v_string infoNames, infoValues;
+    info.get_section_info(i, infoNames, infoValues);
+    PCU_DEBUG_ASSERT(infoNames.size() == infoValues.size());
+    for (std::size_t j = 0; j < infoNames.size(); ++j) {
+      if (secType == "gmodel" && infoNames[j] == "content") {
+        PCU_DEBUG_ASSERT(model_content.empty()); // should only be one gmodel
+        model_content = infoValues[j];
+      } else if (secType == "mmodel" && infoNames[j] == "content") {
+        mesh_contents.push_back(infoValues[j]);
+      }
+    }
+  }
+  PCU_DEBUG_ASSERT(!model_content.empty());
+  PCU_DEBUG_ASSERT(mesh_names.size() == mesh_contents.size());
 }
 
-void gmi_cap_stop(void)
-{
+void gmi_cap_probe(
+  const char* creFileName, std::vector<std::string>& mesh_names
+) {
+  mesh_names.clear();
+	Reader *reader = get_reader(
+    cs_module->get_context(), "Create Native Reader"
+  );
+  DataFileInfo info;
+  if (!reader->probe(creFileName, info))
+    gmi_fail("gmi_cap_probe: failed to read file");
+  for (std::size_t i = 0; i < info.get_num_sections(); i++) {
+    std::string secType, secName;
+    info.get_section(i, secType, secName);
+    if (secType == "mmodel") mesh_names.push_back(secName);
+  }
 }
 
+static struct gmi_model_ops ops;
+
 void gmi_register_cap(void)
 {
   ops.begin = begin;
@@ -368,70 +448,95 @@ void gmi_register_cap(void)
   ops.bbox = bbox;
   ops.is_discrete_ent = is_discrete_ent;
   ops.destroy = destroy;
-  /* gmi_register(create_cre, "cre"); */
+  gmi_register(create_cre, "cre");
   /* gmi_register(create_native, "xmt_txt"); */
   /* gmi_register(create_native, "x_t"); */
   /* gmi_register(create_native, "sat"); */
 }
 
-/* static gmi_model* owned_import(GeometryDatabaseInterface* gi) */
-/* { */
-/*   gmi_model* m = gmi_import_cap(gi); */
-/*   ((cap_model*)m)->owned = true; */
-/*   return m; */
-/* } */
-
-/* struct gmi_model* gmi_cap_load(const char* creFileName) */
-/* { */
-/*   if (!gmi_has_ext(creFileName, "cre")) */
-/*     gmi_fail("gmi_cap_load: cre file must have .cre extension"); */
-/*   const std::string gdbName("Geometry Database : SMLIB");// Switch Create with SMLIB for CAD */
-/*   const std::string mdbName("Mesh Database : Create"); */
-/*   const std::string adbName("Attribution Database : Create"); */
-
-/*   CapstoneModule  cs("test", gdbName.c_str(), mdbName.c_str(), adbName.c_str()); */
-
-/*   GeometryDatabaseInterface     *g = cs.get_geometry(); */
-/*   MeshDatabaseInterface         *m = cs.get_mesh(); */
-/*   AppContext                    *c = cs.get_context(); */
+static gmi_model* owned_import(GDBI* g, M_GModel gmodel) {
+  cap_model* m = reinterpret_cast<cap_model*>(gmi_import_cap(g, gmodel));
+  m->owned = true;
+  return reinterpret_cast<gmi_model*>(m);
+}
 
-/*   PCU_ALWAYS_ASSERT(g); */
-/*   PCU_ALWAYS_ASSERT(m); */
-/*   PCU_ALWAYS_ASSERT(c); */
+struct gmi_model* gmi_cap_load(const char* creFileName) {
+  if (!gmi_has_ext(creFileName, "cre"))
+    gmi_fail("gmi_cap_load: cre file must have .cre extension");
+  if (!cs_module) gmi_fail("gmi_cap_load called before gmi_cap_start");
 
-/*   v_string filenames; */
-/*   filenames.push_back(creFileName); */
+  std::vector<std::string> mesh_names;
+  gmi_cap_probe(creFileName, mesh_names);
+  return gmi_cap_load_selective(creFileName, mesh_names);
+}
 
-/*   M_GModel gmodel = cs.load_files(filenames); */
+struct gmi_model* gmi_cap_load_selective(
+  const char* creFileName, const std::vector<std::string>& mesh_names
+) {
+  if (!gmi_has_ext(creFileName, "cre"))
+    gmi_fail("gmi_cap_load_selective: CRE file must have .cre extension");
+  if (!cs_module) gmi_fail("gmi_cap_load_selective: called before gmi_cap_start");
+  static bool called = false;
+  if (!called) called = true;
+  else {
+    lion_eprint(1,
+      "WARNING: gmi_cap_load_selective called more than once. gmi_cap operations"
+      " may fail.\n"
+    );
+  }
+  auto ctx = cs_module->get_context();
+  FunctionPtr fn(get_function(ctx, "LoadCreateData"));
+  set_input(fn, "FileName", creFileName);
+  set_input(fn, "Meshes", mesh_names);
+  auto proc = get_context_processor(ctx);
+  if (proc->execute(fn) != STATUS_OK)
+    gmi_fail("gmi_cap_load_selective: failed to read CRE file");
+  M_GModel gmodel;
+  get_output(fn, "Model", gmodel);
+  MG_API_CALL(cs_module->get_geometry(), set_current_model(gmodel));
+  return owned_import(cs_module->get_geometry(), gmodel);
+}
 
-/*   int numBreps; */
-/*   g->get_num_breps(numBreps); */
-/*   PCU_ALWAYS_ASSERT(numBreps == 1); */
+void gmi_cap_write(struct gmi_model* model, const char* creFileName) {
+  cap_model* cm = reinterpret_cast<cap_model*>(model);
+  auto ctx = cm->geomInterface->get_context();
+  FunctionPtr fn(get_function(ctx, "SaveCreateData"));
+  set_input(fn, "Model", cm->gmodel);
+  set_input(fn, "FileName", creFileName);
+  auto proc = get_context_processor(ctx);
+  if (proc->execute(fn) != STATUS_OK)
+    gmi_fail("gmi_cap_write: failed to write the model");
+}
 
-/*   return owned_import(g); */
-/* } */
+gmi_model* gmi_import_cap(GDBI* gi) {
+  M_GModel gmodel;
+  MG_API_CALL(gi, get_current_model(gmodel));
+  return gmi_import_cap(gi, gmodel);
+}
 
-gmi_model* gmi_import_cap(GeometryDatabaseInterface* gi)
+gmi_model* gmi_import_cap(GDBI* gi, M_GModel gmodel)
 {
   cap_model* m;
   m = (cap_model*)malloc(sizeof(*m));
   m->model.ops = &ops;
   m->geomInterface = gi;
+  m->gmodel = gmodel;
   M_GBRep brep;
   int numBreps;
 
   m->geomInterface->get_num_breps(numBreps);
-  PCU_ALWAYS_ASSERT(numBreps == 1);
+  if (numBreps != 1)
+    gmi_fail("gmi_import_cap: loaded CRE with numBreps != 1");
   m->geomInterface->get_brep_by_index(0, brep);
-  m->geomInterface->get_num_topos(brep, VERTEX, m->model.n[0]);
-  m->geomInterface->get_num_topos(brep, EDGE,	m->model.n[1]);
-  m->geomInterface->get_num_topos(brep, FACE,	m->model.n[2]);
-  m->geomInterface->get_num_topos(brep, REGION, m->model.n[3]);
+  m->geomInterface->get_num_topos(brep, Geometry::VERTEX, m->model.n[0]);
+  m->geomInterface->get_num_topos(brep, Geometry::EDGE, m->model.n[1]);
+  m->geomInterface->get_num_topos(brep, Geometry::FACE, m->model.n[2]);
+  m->geomInterface->get_num_topos(brep, Geometry::REGION, m->model.n[3]);
   m->owned = false;
   return &m->model;
 }
 
-GeometryDatabaseInterface* gmi_export_cap(gmi_model* m)
+GDBI* gmi_export_cap(gmi_model* m)
 {
   cap_model* cm = (cap_model*)m;
   return cm->geomInterface;
diff --git a/gmi_cap/gmi_cap.h b/gmi_cap/gmi_cap.h
index f2e331d5a..6ee457198 100644
--- a/gmi_cap/gmi_cap.h
+++ b/gmi_cap/gmi_cap.h
@@ -9,29 +9,182 @@
 *******************************************************************************/
 #ifndef GMI_CAP_H
 #define GMI_CAP_H
+/**
+ * \file gmi_cap.h
+ * \brief Concrete Capstone GMI implementation.
+ *
+ * This interface can be used in two ways:
+ *
+ * 1. to import an existing Capstone model; and
+ * 2. to load the model and (optionally) meshes from a CRE file by wrapping
+ *    Capstone calls.
+ *
+ * The first option results in a non-owned model (gmi_cap_start does not need
+ * to be called). The second option results in an owned model (gmi_cap_start
+ * MUST be called) for which destruction is handled by gmi_destroy. Both
+ * options require gmi_register_cap to be called. It is recommended to call
+ * gmi_cap_start (if desired) before gmi_register_cap.
+ */
 
+#ifdef __cplusplus
+#include <string>
+#include <vector>
+#include <CreateMG_CoreDefs.h>
+// Forward declaration
+namespace CreateMG {
+  namespace Geometry { class GeometryDatabaseInterface; }
+  typedef Geometry::GeometryDatabaseInterface GDBI;
+}
+#else
+extern "C" {
+#endif
+// Forward declarations
+struct gmi_model;
+struct gmi_ent;
 
-#include "gmi.h"
-#include "CapstoneModule.h"
-#include "CreateMG_Framework_Geometry.h"
-#include "CreateMG_Framework_Mesh.h" 
-
-using namespace CreateMG;
-using namespace CreateMG::Attribution;
-using namespace CreateMG::Mesh;
-using namespace CreateMG::Geometry;
-
-gmi_ent* toGmiEntity(M_GTopo topo);
-M_GTopo fromGmiEntity(gmi_ent* g);
+/**
+ * \defgroup gmi_cap Capstone geometric model interface
+ *
+ * gmi_cap implements gmi. The interface in gmi_cap.h includes additional
+ * functions to initialize the library and simplify loading of Capstone models.
+ *
+ * \{
+ */
 
+/**
+ * \brief Initialize gmi_cap library.
+ *
+ * \note This call is required before calling any other gmi_cap functions.
+ */
 void gmi_cap_start(void);
+/**
+ * \brief Finalize gmi_cap library.
+ *
+ * \note No gmi_cap functions except gmi_cap_import and gmi_cap_export should
+ * bec called after this function.
+ */
 void gmi_cap_stop(void);
+/**
+ * \brief Register gmi_cap in GMI's global list.
+ *
+ * \note This call is required before every gmi_cap function except
+ * gmi_cap_start.
+ *
+ * This call registers gmi_cap with GMI. It also registers the ".cre" extension
+ * with the Create loader so that gmi_load can be used with ".cre" files.
+ */
 void gmi_register_cap(void);
-
+/**
+ * \brief Load a gmi_model from a CRE file.
+ *
+ * This method is called by gmi_load when the input file has a ".cre"
+ * extension. It loads the file with all associated mesh models. Use
+ * gmi_cap_load_selective to only load certain meshes.
+ *
+ * \param creFileName the name of the CRE file
+ * \pre gmi_cap_stop must be called before this method.
+ * \pre creFileName ends in ".cre"
+ */
 struct gmi_model* gmi_cap_load(const char* creFileName);
-struct gmi_model* gmi_import_cap(GeometryDatabaseInterface* gi);
-GeometryDatabaseInterface* gmi_export_cap(struct gmi_model* m);
+/**
+ * \brief Write a Capstone gmi_model and associated mesh to a CRE file.
+ *
+ * \param model A previously loaded gmi_model from gmi_cap.
+ * \param creFileName filename for the new CRE file.
+ */
+void gmi_cap_write(struct gmi_model* model, const char* creFileName);
+
+#ifdef __cplusplus
+
+/**
+ * \brief Probe a CRE file for associated mesh names.
+ *
+ * This function calls Capstone APIs to probe the CRE file header for meshes.
+ * It is intended to be used in conjunction with gmi_cap_load_selective.
+ *
+ * \param[in]  creFileName CRE file to probe
+ * \param[out] mesh_names A vector to be filled with mesh names
+ */
+void gmi_cap_probe(
+  const char* creFileName, std::vector<std::string>& mesh_names
+);
+/**
+ * \brief Probe a CRE file for associated mesh names and contents.
+ *
+ * This function calls Capstone APIs to probe the CRE file header for model
+ * contents, mesh names, and mesh contents. It is intended to be used in
+ * conjunction with gmi_cap_load_selective, and/or to assist when there are
+ * multiple loadable meshes (i.e. by examining contents or displaying an
+ * interactive menu).
+ *
+ * \warning The format of model_content and mesh_contents is governed by the
+ * Capstone library and we cannot make any claims about its stability.
+ *
+ * \param[in]  creFileName CRE file to probe
+ * \param[out] model_content A string indicating mesh contents
+ * \param[out] mesh_names A vector to be filled with mesh names
+ * \param[out] mesh_contents A vector to be filled with strings indicating mesh
+ *                           contents
+ */
+void gmi_cap_probe(
+  const char* creFileName, std::string& model_content,
+  std::vector<std::string>& mesh_names, std::vector<std::string>& mesh_contents
+);
+/**
+ * \brief Load a gmi_model from a CRE file and select the associated loaded
+ * meshes.
+ *
+ * Unlike other data formats, CRE files contain both geometries and meshes. As
+ * such, the meshes must be selected at load time. The normal gmi_cap_load
+ * (which is registered to gmi_load with the ".cre" extension) loads all
+ * meshes by default. In certain scenarios this is undesirable (e.g. when you
+ * only want the first mesh in a CRE file multiple meshes, or when you only
+ * want to load the mesh on one rank but need the geometry on all ranks).
+ *
+ * \param creFileName The CRE file to load
+ * \param mesh_names A vector of mesh names to load, which may be empty.
+ * \return A new gmi_model linked to a Capstone geometry model
+ */
+struct gmi_model* gmi_cap_load_selective(
+  const char* creFileName, const std::vector<std::string>& mesh_names
+);
 
+gmi_ent* toGmiEntity(CreateMG::M_GTopo topo);
+CreateMG::M_GTopo fromGmiEntity(gmi_ent* g);
+
+/**
+ * \brief Import a non-owned Capstone model.
+ *
+ * Loads the current geometry model on the interface.
+ *
+ * \param gi Geometry interface to load.
+ * \pre gi has a current model.
+ */
+struct gmi_model* gmi_import_cap(CreateMG::GDBI* gi);
+/**
+ * \brief Import a non-owned, specific Capstone model.
+ *
+ * \param gi Geometry interface to load from.
+ * \param gmodel The existing geometry model to load.
+ */
+struct gmi_model* gmi_import_cap(
+  CreateMG::GDBI* gi, CreateMG::M_GModel gmodel
+);
+/**
+ * \brief Expose the underlying Capstone geometry interface from a gmi_model.
+ *
+ * \param m Loaded or imported gmi_model.
+ */
+CreateMG::GDBI* gmi_export_cap(struct gmi_model* m);
+
+#endif // __cplusplus
+
+
+#ifndef __cplusplus
+} // extern "C"
 #endif
 
+/** \} */
+
+#endif // GMI_CAP_H
 
diff --git a/ma/maSnap.cc b/ma/maSnap.cc
index e754af553..4fea96021 100644
--- a/ma/maSnap.cc
+++ b/ma/maSnap.cc
@@ -68,7 +68,7 @@ static size_t isSurfUnderlyingFaceDegenerate(
     m->getFirstDerivative(g, p, uTan, vTan);
     double uTanSize = uTan.getLength();
     double vTanSize = vTan.getLength();
-#ifdef HAVE_CAPSTONE
+#ifdef PUMI_HAS_CAPSTONE
     uTanSize = uTan * uTan;
     vTanSize = vTan * vTan;
 #endif
@@ -150,7 +150,7 @@ static void interpolateParametricCoordinateOnEdge(
   p[1] = 0.0;
   p[2] = 0.0;
 
-#ifdef HAVE_CAPSTONE
+#ifdef PUMI_HAS_CAPSTONE
   // account for non-uniform parameterization of model-edge
   Vector X[3];
   Vector para[2] = {a, b};
@@ -491,7 +491,7 @@ static void interpolateParametricCoordinatesOnRegularFace(
    * 2) we only check for faces that are periodic
    */
 
-#ifndef HAVE_CAPSTONE
+#ifndef PUMI_HAS_CAPSTONE
   // this need to be done for faces, only
   if (dim != 2)
     return;
@@ -526,7 +526,7 @@ static void interpolateParametricCoordinatesOnFace(
   size_t num = isSurfUnderlyingFaceDegenerate(m, g, axes, vals);
 
   if (num > 0) { // the underlying surface is degenerate
-#ifndef HAVE_CAPSTONE
+#ifndef PUMI_HAS_CAPSTONE
     interpolateParametricCoordinatesOnDegenerateFace(m, g, t, a, b, axes, vals, p);
 #else
     // account for non-uniform parameterization of model-edge
diff --git a/pumi-meshes b/pumi-meshes
index 7be0e261d..27dca7523 160000
--- a/pumi-meshes
+++ b/pumi-meshes
@@ -1 +1 @@
-Subproject commit 7be0e261d69d45a6f0614c76c66eff2bba81f4ba
+Subproject commit 27dca75239ed386594b7eb35d69c2d58fb28cfe3
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 85b3db547..f6d2ec6aa 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -216,19 +216,25 @@ if(ENABLE_SIMMETRIX)
   test_exe_func(highOrderSizeFields highOrderSizeFields.cc)
 endif()
 
-if(ENABLE_CAPSTONE)
+if(PUMI_ENABLE_CAPSTONE)
   util_exe_func(capVol capVol.cc)
-  target_include_directories(capVol PRIVATE "${PROJECT_SOURCE_DIR}/capstone_clis")
   test_exe_func(cap_inClosureOf cap_inClosureOf.cc)
   test_exe_func(cap_closestPoint cap_closestPoint.cc)
-  if(HAVE_CAPSTONE_SIZINGMETRICTOOL)
+  if(PUMI_HAS_CAPSTONE_SIZINGMETRICTOOL)
     util_exe_func(cap_smooth cap_smooth.cc)
   endif()
   util_exe_func(cap2vtk cap2vtk.cc)
   util_exe_func(capGeomTest capGeomTest.cc)
   util_exe_func(capCheckParam capCheckParam.cc)
   util_exe_func(capAdapt capAdapt.cc)
-  target_link_libraries(capAdapt framework_meshing)
+  util_exe_func(capProbe capProbe.cc)
+  util_exe_func(capLoadSome capLoadSome.cc)
+  if(PUMI_HAS_CAPSTONE_SIZINGMETRICTOOL)
+    util_exe_func(capNativeAdapt capNativeAdapt.cc)
+    target_link_libraries(capNativeAdapt
+      framework_application framework_mesh framework_meshing
+    )
+  endif()
 endif()
 
 # send all the newly added utility executable targets
diff --git a/test/cap2vtk.cc b/test/cap2vtk.cc
index 1b2c696e7..954b77771 100644
--- a/test/cap2vtk.cc
+++ b/test/cap2vtk.cc
@@ -1,143 +1,79 @@
+#include <iostream>
+#include <stdexcept>
+
 #include <PCU.h>
+#include <apf.h>
 #include <apfCAP.h>
-#include <apfMDS.h>
-#include <samSz.h>
-#include <queue>
-#include <ma.h>
+#include <apfMesh2.h>
 #include <gmi.h>
 #include <gmi_cap.h>
-#include <gmi_null.h>
-#include <gmi_mesh.h>
-#include <apf.h>
-#include <apfConvert.h>
-#include <apfMesh2.h>
-#include <apfNumbering.h>
-#include <apfShape.h>
-#include <ma.h>
+#include <lionPrint.h>
 #include <pcu_util.h>
-#include <cstdlib>
-#include <cstring>
-#include <iostream>
-#include <iomanip>
-#include <vector>
-#include <math.h>
-
-
-#include "CapstoneModule.h"
-#include "CreateMG_Framework_Core.h"
-#include "CreateMG_Framework_Analysis.h"
-#include "CreateMG_Framework_Application.h"
-#include "CreateMG_Framework_Attributes.h"
-#include "CreateMG_Framework_Core.h"
-#include "CreateMG_Framework_Geometry.h"
-#include "CreateMG_Framework_Mesh.h"
-
-using namespace CreateMG;
-using namespace CreateMG::Attribution;
-using namespace CreateMG::Mesh;
-using namespace CreateMG::Geometry;
-
 
 int main(int argc, char** argv)
 {
+  int retval = 0;
   pcu::Init(&argc, &argv);
   {
   pcu::PCU PCUObj;
-
-  if (argc != 3) {
-    if(0==PCUObj.Self())
-      std::cerr << "usage: " << argv[0]
-        << " <cre file .cre> <output folder name>\n";
-    return EXIT_FAILURE;
-  }
-
-
-  gmi_register_mesh();
-  gmi_register_null();
-
-  const char* creFileName = argv[1];
-  const char* folderName = argv[2];
-
-  // load capstone mesh
-  // create an instance of the Capstone Module activating CREATE/CREATE/CREATE
-  // for the Geometry/Mesh/Attribution databases
-  /* const std::string gdbName("Geometry Database : Create");// Switch Create with SMLIB for CAD */
-  const std::string gdbName("Geometry Database : SMLIB");// Switch Create with SMLIB for CAD
-  const std::string mdbName("Mesh Database : Create");
-  const std::string adbName("Attribution Database : Create");
-
-  CapstoneModule  cs("the_module", gdbName.c_str(), mdbName.c_str(), adbName.c_str());
-
-  GeometryDatabaseInterface     *g = cs.get_geometry();
-  MeshDatabaseInterface         *m = cs.get_mesh();
-  AppContext                    *c = cs.get_context();
-
-
-  PCU_ALWAYS_ASSERT(g);
-  PCU_ALWAYS_ASSERT(m);
-  PCU_ALWAYS_ASSERT(c);
-
-  v_string filenames;
-  filenames.push_back(creFileName);
-
-  M_GModel gmodel = cs.load_files(filenames);
-
-  int numbreps = 0;
-  MG_CALL(g->get_num_breps(numbreps));
-  std::cout << "number of b reps is " << numbreps << std::endl;
-  if(numbreps == 0)
-      error(HERE, ERR_INVALID_INPUT, "Model is empty");
-
-  M_MModel mmodel;
-  // Pick the volume mesh model from associated mesh models to this geom model
-  std::vector<M_MModel> mmodels;
-  MG_API_CALL(m, get_associated_mesh_models(gmodel, mmodels));
-  for(std::size_t i = 0; i < mmodels.size(); ++i)
-  {
-      M_MModel ammodel = mmodels[i];
-      std::size_t numregs = 0;
-      std::size_t numfaces = 0;
-      std::size_t numedges = 0;
-      std::size_t numverts = 0;
-      MG_API_CALL(m, set_current_model(ammodel));
-      MG_API_CALL(m, get_num_topos(TOPO_REGION, numregs));
-      MG_API_CALL(m, get_num_topos(TOPO_FACE, numfaces));
-      MG_API_CALL(m, get_num_topos(TOPO_EDGE, numedges));
-      MG_API_CALL(m, get_num_topos(TOPO_VERTEX, numverts));
-      std::cout << "num regions is " << numregs << std::endl;
-      std::cout << "num faces   is " << numfaces << std::endl;
-      std::cout << "num edges   is " << numedges << std::endl;
-      std::cout << "num verts   is " << numverts << std::endl;
-      std::cout << "-----------" << std::endl;
-      if(numregs > 0)
-      {
-	  mmodel = ammodel;
-	  break;
+  lion_set_verbosity(1);
+  try {
+    if (PCUObj.Peers() > 1)
+      throw std::runtime_error("Capstone meshes are (currently) only serial");
+
+    if (argc != 3 && argc != 4) {
+      if (PCUObj.Self() == 0) {
+        std::cerr << "USAGE: " << argv[0]
+          << " <input.cre> [meshname] <output.vtk>" << std::endl;
       }
+      throw std::invalid_argument("invalid arguments");
+    }
+
+    const char* creFileName = argv[1];
+    const char* meshName = argc == 4 ? argv[2] : nullptr;
+    const char* folderName = argc == 4 ? argv[3] : argv[2];
+
+    gmi_cap_start();
+    gmi_register_cap();
+    try {
+      std::string model_content;
+      std::vector<std::string> mesh_names, mesh_contents;
+      gmi_cap_probe(creFileName, model_content, mesh_names, mesh_contents);
+      if (meshName) {
+        if (std::count(mesh_names.begin(), mesh_names.end(), meshName) == 0) {
+          throw std::invalid_argument(
+            "mesh named " + std::string(meshName) + " not found"
+          );
+        }
+      } else if (mesh_names.size() > 1) {
+        std::cout << "More than one mesh in CRE file."
+          " Run again with meshname argument after CRE filename." << std::endl;
+        for (size_t i = 0; i < mesh_names.size(); ++i) {
+          std::cout << "Mesh " << i << " `" << mesh_names[i] << "`: "
+            << mesh_contents[i] << std::endl;
+        }
+        throw std::runtime_error(
+          "meshname not provided for file with multiple meshes"
+        );
+      } else meshName = mesh_names.front().c_str();
+      PCU_DEBUG_ASSERT(meshName);
+      gmi_model* model = gmi_cap_load_selective(creFileName, {meshName});
+      apf::Mesh2* mesh = apf::createCapMesh(model, meshName, &PCUObj);
+      apf::writeVtkFiles(folderName, mesh);
+      apf::destroyMesh(mesh);
+    } catch (...) {
+      gmi_cap_stop();
+      std::rethrow_exception(std::current_exception());
+    }
+    gmi_cap_stop();
+  } catch (const std::exception& e) {
+    if (PCUObj.Self() == 0) std::cerr << "ERROR: " << e.what() << std::endl;
+    retval = 1;
+  } catch (...) {
+    if (PCUObj.Self() == 0) std::cerr << "UNKNOWN ERROR" << std::endl;
+    retval = 1;
   }
-
-  /* SET THE ADJACENCIES */
-  MG_API_CALL(m, set_adjacency_state(REGION2FACE|
-                                     REGION2EDGE|
-                                     REGION2VERTEX|
-                                     FACE2EDGE|
-                                     FACE2VERTEX));
-  MG_API_CALL(m, set_reverse_states());
-  MG_API_CALL(m, set_adjacency_scope(TOPO_EDGE, SCOPE_FULL));
-  MG_API_CALL(m, set_adjacency_scope(TOPO_FACE, SCOPE_FULL));
-  MG_API_CALL(m, compute_adjacency());
-
-
-  gmi_cap_start();
-  gmi_register_cap();
-
-  // convert the mesh to apf/mds mesh
-
-  apf::Mesh2* mesh = apf::createMesh(m,g,&PCUObj);
-
-  apf::writeVtkFiles(folderName, mesh);
-
-  gmi_cap_stop();
-  }
+  } // PCU Object scope
   pcu::Finalize();
+  return retval;
 }
diff --git a/test/capAdapt.cc b/test/capAdapt.cc
index b5c9b6aab..f194691bc 100644
--- a/test/capAdapt.cc
+++ b/test/capAdapt.cc
@@ -9,12 +9,9 @@
  * \author Eric Mestreau
  */
 #include <iostream>
-#include <chrono>
 #include <exception>
 #include <memory>
 #include <stdexcept>
-#include <thread>
-#include <type_traits>
 
 #include <PCU.h>
 #include <apf.h>
@@ -33,113 +30,47 @@
 #include <apfMETIS.h>
 #include <apfZoltan.h>
 
-#include "CapstoneModule.h"
-#include "CreateMG_SizingMetricTool.h"
-
-using namespace CreateMG;
-using namespace CreateMG::Attribution;
-using namespace CreateMG::Mesh;
-using namespace CreateMG::Geometry;
-
 namespace {
 
 /** \brief Print nested exceptions. */
 void print_exception(const std::exception& e, int level = 0);
 
-/** \brief Load bulk file of one type into std::vector. */
-template <class T>
-typename std::enable_if<std::is_trivially_copyable<T>::value, void>::type
-load_file(const std::string& fname, std::vector<T>& data) {
-  try {
-    std::ifstream file(fname, std::ios::in | std::ios::binary);
-    T val;
-    while (file) {
-      file.read(reinterpret_cast<char*>(std::addressof(val)), sizeof(T));
-      if (file.gcount() == sizeof(T)) data.push_back(val);
-    }
-  } catch(...) {
-    std::throw_with_nested(
-      std::runtime_error("failed to load the file `" + fname + "`")
-    );
-  }
-}
-
 /** \brief Return the filename without the extension. */
 std::string get_no_extension(const std::string& fname) {
   return fname.substr(0, fname.find_last_of("."));
 }
 
-/** \brief Decompose a 6 component metric tensor into frames and scales. */
-void decompose_metric(const v_double& metric, apf::Matrix3x3& frame,
-  apf::Vector3& scalar) {
-  apf::Matrix3x3 m(metric[0], metric[1], metric[2],
-    metric[1], metric[3], metric[4],
-    metric[2], metric[4], metric[5]);
-  int n = apf::eigen(m, &frame[0], &scalar[0]);
-  PCU_DEBUG_ASSERT(n == 3);
-}
-
-/**
- * \brief Load sizing metric in mesh vertex iteration order.
- * 
- * \param m Capstone mesh database.
- * \param mapFileName path to vmap file with vertex mappings (for ordering).
- * \param sizFielName path to bulk sizing file.
- * */
-std::vector<Metric6> loadSizing(
-  MeshDatabaseInterface *m, std::string mapFileName, std::string sizFileName
-);
-
-/** \brief Load Metric6 vector into APF frame and scale fields. */
-void convertMetric6toAnisoFields(
-  MeshDatabaseInterface *m, const std::vector<Metric6>& sizing6,
-  apf::Field* frameField, apf::Field* scaleField
-);
-
 /** \brief Print edge cuts. */
 void printEdgeCuts(apf::Mesh2 *m);
 
-/** \brief Create an apf::Mesh2* from CapstoneModule and bulk sizing file. */
-ma::Mesh* makeApfInterfaceWithSizing(
-  CapstoneModule& cs, std::unique_ptr<pcu::PCU>& pcu, bool original,
-  const std::string& vmapFile, const std::string& sizingFile, bool smooth
-);
-
 /** \brief Migrate all elements back to rank 0. */
 void migrateHome(apf::Mesh2* mesh);
 
-/** \brief Take a serial MDS mesh, partition, adapt, then localize. */
-void parallelAdapt(ma::Mesh* mesh, gmi_model* mod, pcu::PCU& PCU, int maxiter);
-
 /** \brief Command line argument processing class. */
-class Args {
+struct Args {
 public:
   Args() {}
   Args(int argc, char* argv[]) { parse(argc, argv); }
   void parse(int argc, char* argv[]);
   void print_usage(const char* argv0) const;
   /**
-   * \brief Get maximum iterations argument.
+   * \brief maximum iterations argument.
    *
    * Pass -1 to trust MeshAdapt or a number <= 10 to specify max iterations.
    */
-  int maxiter() const noexcept { return maxiter_; }
+  int maxiter{-1};
   /**
    * \brief Get smoothing argument. If true, request smoothing if supported.
    */
-  bool smooth() const noexcept { return smooth_; }
-
-  const std::string& creFilename() const noexcept { return creFilename_; }
-  const std::string& sizingFilename() const noexcept {
-    return sizingFilename_;
-  }
-private:
-  int maxiter_{-1};
-  bool smooth_{false};
-  std::string creFilename_;
-  std::string sizingFilename_;
+  bool smooth{false};
+  std::string creFilename, sizingFilename, before, after;
 };
 
+/** \brief Take a serial MDS mesh, partition, adapt, then localize. */
+void parallelAdapt(
+  ma::Mesh* mesh, gmi_model* mod, pcu::PCU& PCU, const Args& args
+);
+
 } // namespace
 
 int main(int argc, char** argv) {
@@ -167,112 +98,73 @@ int main(int argc, char** argv) {
     /* LOAD CAPSTONE MESH */
     if (PCU.Self() == 0)
       std::cout << "STATUS: Init Capstone" << std::endl;
-    const std::string analysis("Kestrel");
-    CapstoneModule cs("capAdapt",
-      "Geometry Database : SMLIB",
-      "Mesh Database : Create",
-      "Attribution Database : Create"
-    );
-    MeshDatabaseInterface     *m = cs.get_mesh();
-    PCU_ALWAYS_ASSERT(m);
 
     if (PCU.Self() == 0)
-      std::cout << "STATUS: Loading CRE file : " << args.creFilename()
-        << std::endl;
-    M_GModel gmodel;
+      std::cout << "STATUS: Loading CRE file : " << args.creFilename;
+    gmi_model* model = nullptr;
     // Load CRE file sequentially to avoid NFS errors.
-    for (int p = 0; p < PCU.Peers(); ++p) {
-      if (PCU.Self() == p) {
-        gmodel = cs.load_files(v_string(1, args.creFilename()));
+    for (int i = 0; i < PCU.Peers(); ++i) {
+      if (PCU.Self() == i) {
+        if (PCU.Self() == 0) model = gmi_load(args.creFilename.c_str());
+        else model = gmi_cap_load_selective(args.creFilename.c_str(), {});
       }
       PCU.Barrier();
     }
 
-    // Load mesh.
-    M_MModel mmodel;
-    std::vector<M_MModel> mmodels;
-    MG_API_CALL(m, get_associated_mesh_models(gmodel, mmodels));
-    PCU_ALWAYS_ASSERT(mmodels.size() == 1);
-    MG_API_CALL(m, set_current_model(mmodels[0]));
-    mmodel = mmodels[0];
-    cs.set_analysis_name(analysis);
-
-    /* SET THE ADJACENCIES */
-    if (PCU.Self() == 0) {
-      std::cout << "STATUS: Compute required adjacencies" << std::endl;
-      MG_API_CALL(m, set_adjacency_state(REGION2FACE | REGION2EDGE |
-        REGION2VERTEX | FACE2EDGE | FACE2VERTEX));
-      MG_API_CALL(m, set_reverse_states());
-      MG_API_CALL(m, compute_adjacency());
-    }
-
     // Create apf::Mesh2 wrapper.
     if (PCU.Self() == 0)
       std::cout << "STATUS: Create apf::Mesh interface" << std::endl;
 
-    ma::Mesh* apfCapMesh = nullptr, *apfMesh = nullptr;
+    ma::Mesh* apfMesh = nullptr;
     bool original = PCU.Self() == 0;
-    int parts = PCU.Peers();
     auto soloPCU = PCU.Split(PCU.Self(), 0);
 
-    apfCapMesh = makeApfInterfaceWithSizing(
-      cs, soloPCU, original,
-      get_no_extension(args.creFilename()) + ".avm.vmap",
-      args.sizingFilename(), args.smooth()
-    );
     if (original) {
+      ma::Mesh* apfCapMesh = apf::createCapMesh(model, soloPCU.get());
+      apf::disownCapModel(apfCapMesh);
+      if (!apf::loadCapSizingFile(
+        apfCapMesh, args.sizingFilename,
+        get_no_extension(args.creFilename) + ".avm.vmap",
+        "adapt_scales", "adapt_frames", true, "Kestrel"
+      )) {
+        throw std::runtime_error("failed to load sizing file");
+      }
       /* CONVERT TO APF MDS MESH AND WRITE VTK */
       std::cout << "STATUS: Convert to APF MDS Mesh" << std::endl;
       double t_start = pcu::Time();
-      apfMesh = apf::createMdsMesh(apfCapMesh->getModel(), apfCapMesh, true);
+      apfMesh = apf::createMdsMesh(model, apfCapMesh, true);
       double t_end = pcu::Time();
       std::cout << "INFO: Mesh converted in " << t_end - t_start << " seconds."
         << std::endl;
+      apf::disownMdsModel(apfMesh);
+      apf::destroyMesh(apfCapMesh);
       apfMesh->verify();
       apf::printStats(apfMesh);
-      apf::writeVtkFiles("before-mds", apfMesh);
     }
 
     PCU.Barrier();
 
-    parallelAdapt(apfMesh, apfCapMesh->getModel(), PCU, args.maxiter());
+    parallelAdapt(apfMesh, model, PCU, args);
 
     if (original) {
-      /* COPY THE MESH MODEL TO KEEP */
-      MG_API_CALL(m, create_associated_model(mmodel, gmodel, "MeshAdapt"));
-
-      /* SETUP NEW MESH ADJACENCIES */
-      MG_API_CALL(m, set_adjacency_state(REGION2FACE|REGION2EDGE|
-        REGION2VERTEX|FACE2EDGE|FACE2VERTEX));
-      MG_API_CALL(m, set_reverse_states());
-      MG_API_CALL(m, compute_adjacency());
-
-      // Destroy old fields.
-      apf::destroyField(apfCapMesh->findField("adapt_scales"));
-      apf::destroyField(apfCapMesh->findField("adapt_frames"));
-
+      apf::Mesh2* mdsMesh = apfMesh;
+      apfMesh = apf::makeEmptyCapMesh(model, "MeshAdapt", soloPCU.get());
+      apf::disownCapModel(apfMesh);
       auto t0 = pcu::Time();
       /* WRITE MDS MESH TO CRE */
-      apf::convert(apfMesh, apfCapMesh);
+      apf::convert(mdsMesh, apfMesh);
       auto t1 = pcu::Time();
       std::cout << "Converting MDS to CRE: " << t1 - t0 << " seconds"
         << std::endl;
+      apf::destroyMesh(mdsMesh);
       
       /* SAVE FINAL CRE FILE */
-      cs.save_mesh_file("after_cap.vtk", mmodel);
-      auto creOFileName = get_no_extension(args.creFilename())
+      auto creOFileName = get_no_extension(args.creFilename)
         + "_adapted.cre";
-      cs.save_file(creOFileName,gmodel);
-      
-      /* PRINT ADAPTED MESH INFO */
-      std::string info;
-      m->print_info(mmodel, info);
-      std::cout << info << std::endl;
-
-      apf::destroyMesh(apfCapMesh);
+      gmi_cap_write(model, creOFileName.c_str());
       apf::destroyMesh(apfMesh);
     }
-    PCU.Barrier();
+    gmi_destroy(model);
   } catch (const std::exception& e) {
     std::cerr << "FATAL: ";
     print_exception(e);
@@ -297,125 +189,6 @@ void print_exception(const std::exception& e, int level) {
   } catch (...) {}
 }
 
-std::vector<Metric6> loadSizing(
-  MeshDatabaseInterface *m, std::string mapFileName, std::string sizFileName
-) {
-  std::vector<std::size_t> vmap;
-  std::vector<double>      sizing;
-  std::cout << "STATUS: Loading VMap file: " << mapFileName << std::endl;
-  load_file<std::size_t>(mapFileName, vmap);
-  std::cout << "INFO: VMap number of vertices = " << vmap.size() - 1
-    << "    vmap[0] = " << vmap[0]  << std::endl;
-  std::cout << "STATUS: Loading Sizing File: " << sizFileName
-    << std::endl;
-  load_file<double>(sizFileName, sizing);
-  std::cout << "INFO: Sizing File number of vertices = "
-    << sizing.size()/6 << std::endl;
-  std::size_t numVertices;
-  MG_API_CALL(m, get_num_topos(Mesh::TOPO_VERTEX, numVertices));
-  if (sizing.size() / 6 != numVertices) {
-    throw std::runtime_error(
-      "sizing field does not match mesh: " + std::to_string(sizing.size() / 6)
-      + " / " + std::to_string(numVertices)
-    );
-  }
-
-  /* COPY SIZEFIELD TO METRIC6 */
-  std::vector<Metric6> sizing6(numVertices);
-  for (size_t si = 0, i = 1; i < vmap.size(); ++i) {
-    size_t tid = vmap[i];
-    PCU_DEBUG_ASSERT(tid != 0);
-    M_MTopo mtopo = m->get_topo_by_id(Mesh::TOPO_VERTEX, tid);
-    for (size_t j = 0; j < 6; ++j, ++si) {
-      sizing6[tid - 1][j] = sizing[si];
-    }
-  }
-  return sizing6;
-}
-
-void convertMetric6toAnisoFields(
-  MeshDatabaseInterface *m, const std::vector<Metric6>& sizing6,
-  apf::Field* frameField, apf::Field* scaleField
-) {
-  v_double       sz(6);
-  apf::Vector3   scalar;
-  apf::Matrix3x3 frame;
-  for (std::size_t i=0; i<sizing6.size(); i++) {
-    M_MTopo mtopo = m->get_topo_by_id(Mesh::TOPO_VERTEX,i + 1);
-    for (std::size_t j=0; j<6; j++) sz[j] = sizing6[i][j];
-    
-    // field eigen values
-    decompose_metric(sz,frame,scalar);
-    for (int i=0; i<3; i++) {
-      scalar[i] = 1.0/sqrt(scalar[i]);
-    }
-    // MeshAdapt wants frames on the columns.
-    frame = apf::transpose(frame);
-    apf::MeshEntity* ent = (apf::MeshEntity*)mtopo.get();
-    apf::setVector(scaleField, ent, 0, scalar);
-    apf::setMatrix(frameField, ent, 0, frame);
-  }
-}
-
-ma::Mesh* makeApfInterfaceWithSizing(
-  CapstoneModule& cs, std::unique_ptr<pcu::PCU>& pcu, bool original,
-  const std::string& vmapFile, const std::string& sizingFile, bool smooth
-) {
-  // All ranks need apfCapMesh->getModel().
-  ma::Mesh* apfCapMesh = apf::createMesh(
-    cs.get_mesh(), cs.get_geometry(), pcu.get()
-  );
-
-  if (original) {
-    /* LOAD THE SIZING FILE AND VMAP FILE                        */
-    /* CHECK VALIDITY OF THE SIZING FIELD AGAINST THE MESH MODEL */
-    std::cout << "STATUS: Loading sizing" << std::endl;
-    std::vector<Metric6> sizing6 = loadSizing(
-      cs.get_mesh(), vmapFile, sizingFile
-    );
-
-    if (smooth) {
-      /* SMOOTH SIZING FIELD */
-      std::cout << "STATUS: Smoothing sizing field." << std::endl;
-      auto smooth_tool = get_sizing_metric_tool(
-        cs.get_context(), "CreateSmoothingBase"
-      );
-      if (smooth_tool == nullptr) {
-        throw std::runtime_error("unable to find \"CreateSmoothingBase\"");
-      } else {
-        smooth_tool->set_context(cs.get_context());
-        M_MModel mmodel;
-        MG_API_CALL(cs.get_mesh(), get_current_model(mmodel));
-        smooth_tool->set_metric(mmodel, "sizing6", sizing6);
-        std::vector<Metric6> ometric;
-        auto analysis = cs.get_analysis_name();
-        smooth_tool->smooth_metric(mmodel, analysis, "sizing6", ometric);
-        sizing6 = ometric;
-      }
-    } else {
-      std::cout << "INFO: not smoothing sizing field." << std::endl;
-    }
-
-    /* CREATE THE ANISOTROPIC FIELD */
-    std::cout << "STATUS: Creating Anisotropic size-field from the file: "
-      << sizingFile << std::endl;
-    apf::Field* frameField = apf::createFieldOn(
-      apfCapMesh, "adapt_frames", apf::MATRIX
-    );
-    apf::Field* scaleField = apf::createFieldOn(
-      apfCapMesh, "adapt_scales", apf::VECTOR
-    );
-    convertMetric6toAnisoFields(
-      cs.get_mesh(), sizing6, frameField, scaleField
-    );
-
-    /* WRITE THE BEFORE ADAPT MESH TO VTK USING APF VTK WRITER */
-    std::cout << "STATUS: Writing VTK file (before)" << std::endl;
-    apf::writeVtkFiles("before", apfCapMesh);
-  }
-  return apfCapMesh;
-}
-
 void migrateHome(apf::Mesh2* mesh) {
   auto t0 = pcu::Time();
   apf::Migration* plan = new apf::Migration(mesh);
@@ -436,7 +209,9 @@ void migrateHome(apf::Mesh2* mesh) {
       << " seconds" << std::endl;
 }
 
-void parallelAdapt(ma::Mesh* mesh, gmi_model* model, pcu::PCU& PCU, int maxiter) {
+void parallelAdapt(
+  ma::Mesh* mesh, gmi_model* model, pcu::PCU& PCU, const Args& args
+) {
   auto t0 = pcu::Time();
   bool original = false;
   pcu::PCU* oldPCU = nullptr;
@@ -474,10 +249,12 @@ void parallelAdapt(ma::Mesh* mesh, gmi_model* model, pcu::PCU& PCU, int maxiter)
     );
     plan = nullptr; // plan is freed by apf::repeatMdsMesh.
     apf::printStats(mesh);
-    apf::writeVtkFiles("before-mds-split", mesh);
   }
   apf::disownMdsModel(mesh);
   printEdgeCuts(mesh);
+  if (!args.before.empty()) {
+    apf::writeVtkFiles(args.before.c_str(), mesh);
+  }
 
   apf::Field* mdsScaleField = mesh->findField("adapt_scales");
   apf::Field* mdsFrameField = mesh->findField("adapt_frames");
@@ -491,14 +268,16 @@ void parallelAdapt(ma::Mesh* mesh, gmi_model* model, pcu::PCU& PCU, int maxiter)
     ma::configure(mesh, mdsScaleField, mdsFrameField)
   );
   in->shouldForceAdaptation = true;
-  if (maxiter != -1) in->maximumIterations = maxiter;
+  if (args.maxiter != -1) in->maximumIterations = args.maxiter;
 
   if (original) std::cout << "STATUS: Adapting" << std::endl;
   ma::adapt(in);
 
   /* WRITE THE AFTER ADAPT MESH TO VTK USING APF VTK WRITER */
-  if (original) std::cout << "STATUS: Writing VTK file (after)" << std::endl;
-  apf::writeVtkFiles("after-mds", mesh);
+  if (!args.after.empty()) {
+    if (original) std::cout << "STATUS: Writing VTK file (after)" << std::endl;
+    apf::writeVtkFiles(args.after.c_str(), mesh);
+  }
 
   if (parts > 1) migrateHome(mesh);
   if (original) {
@@ -542,20 +321,30 @@ void Args::parse(int argc, char* argv[]) {
   for (int i = 1; i < argc; ++i) {
     if (argv[i][0] == '-') {
       switch (argv[i][1]) {
+      case 'B':
+        if (argv[i][2]) before = &argv[i][2];
+        else if (++i < argc) before = argv[i];
+        else throw std::invalid_argument("missing argument to -B");
+        break;
+      case 'A':
+        if (argv[i][2]) after = &argv[i][2];
+        else if (++i < argc) after = argv[i];
+        else throw std::invalid_argument("missing argument to -A");
+        break;
       case 'i':
-        if (argv[i][2]) maxiter_ = std::stoi(&argv[i][2]);
-        else if (++i < argc) maxiter_ = std::stoi(argv[i]);
+        if (argv[i][2]) maxiter = std::stoi(&argv[i][2]);
+        else if (++i < argc) maxiter = std::stoi(argv[i]);
         else throw std::invalid_argument("missing argument to -i");
         break;
       case 's':
-        smooth_ = true;
+        smooth = true;
         break;
       default:
         throw std::invalid_argument("invalid option -" + argv[i][1]);
       }
     } else {
-      if (creFilename_.empty()) creFilename_ = argv[i];
-      else if (sizingFilename_.empty()) sizingFilename_ = argv[i];
+      if (creFilename.empty()) creFilename = argv[i];
+      else if (sizingFilename.empty()) sizingFilename = argv[i];
       else {
         throw std::invalid_argument(
           std::string("invalid argument `") + argv[i] + "`"
@@ -566,8 +355,8 @@ void Args::parse(int argc, char* argv[]) {
 }
 
 void Args::print_usage(const char* argv0) const {
-  std::cout << "USAGE: " << argv0 << " [-i MAXITER] [-s] INPUT.CRE SIZING.DAT"
-    << std::endl;
+  std::cout << "USAGE: " << argv0 << " [-B BEFORE.VTK] [-A AFTER.VTK]"
+    " [-i MAXITER] [-s] INPUT.CRE SIZING.DAT" << std::endl;
 }
 
 } // namsepace
diff --git a/test/capCheckParam.cc b/test/capCheckParam.cc
index a349bff2c..d008f3dbe 100644
--- a/test/capCheckParam.cc
+++ b/test/capCheckParam.cc
@@ -1,29 +1,15 @@
-#include <PCU.h>
-#include <queue>
-#include <apf.h>
-#include <pcu_util.h>
-#include <cstdlib>
-#include <cstring>
 #include <iostream>
 #include <iomanip>
-#include <vector>
 
+#include <PCU.h>
+#include <apf.h>
+#include <apfCAP.h>
+#include <apfMesh2.h>
+#include <gmi_cap.h>
+#include <lionPrint.h>
+#include <pcu_util.h>
 
-#include "CapstoneModule.h"
-#include "CreateMG_Framework_Core.h"
-#include "CreateMG_Framework_Analysis.h"
-#include "CreateMG_Framework_Application.h"
-#include "CreateMG_Framework_Attributes.h"
-#include "CreateMG_Framework_Core.h"
-#include "CreateMG_Framework_Geometry.h"
-#include "CreateMG_Framework_Mesh.h"
-
-using namespace CreateMG;
-using namespace CreateMG::Attribution;
-using namespace CreateMG::Mesh;
-using namespace CreateMG::Geometry;
-
-void checkParametrization(MeshDatabaseInterface* mdb, GeometryDatabaseInterface* gdb);
+void checkParametrization(apf::Mesh* mesh);
 
 int main(int argc, char** argv)
 {
@@ -39,117 +25,67 @@ int main(int argc, char** argv)
 
   const char* creFileName = argv[1];
 
-  // load capstone mesh
-  // create an instance of the Capstone Module activating CREATE/CREATE/CREATE
-  // for the Geometry/Mesh/Attribution databases
-  /* const std::string gdbName("Geometry Database : Create");// Switch Create with SMLIB for CAD */
-  const std::string gdbName("Geometry Database : SMLIB");// Switch Create with SMLIB for CAD
-  const std::string mdbName("Mesh Database : Create");
-  const std::string adbName("Attribution Database : Create");
-
-  CapstoneModule  cs("the_module", gdbName.c_str(), mdbName.c_str(), adbName.c_str());
-
-  GeometryDatabaseInterface     *g = cs.get_geometry();
-  MeshDatabaseInterface         *m = cs.get_mesh();
-  AppContext                    *c = cs.get_context();
-
-
-  PCU_ALWAYS_ASSERT(g);
-  PCU_ALWAYS_ASSERT(m);
-  PCU_ALWAYS_ASSERT(c);
-
-  v_string filenames;
-  filenames.push_back(creFileName);
-
-  M_GModel gmodel = cs.load_files(filenames);
-
-  int numbreps = 0;
-  MG_CALL(g->get_num_breps(numbreps));
-  std::cout << "number of b reps is " << numbreps << std::endl;
-  if(numbreps == 0)
-      error(HERE, ERR_INVALID_INPUT, "Model is empty");
-
-  M_MModel mmodel;
-  // Pick the volume mesh model from associated mesh models to this geom model
-  std::vector<M_MModel> mmodels;
-  MG_API_CALL(m, get_associated_mesh_models(gmodel, mmodels));
-  for(std::size_t i = 0; i < mmodels.size(); ++i)
-  {
-      M_MModel ammodel = mmodels[i];
-      std::size_t numregs = 0;
-      std::size_t numfaces = 0;
-      std::size_t numedges = 0;
-      std::size_t numverts = 0;
-      MG_API_CALL(m, set_current_model(ammodel));
-      MG_API_CALL(m, get_num_topos(TOPO_REGION, numregs));
-      MG_API_CALL(m, get_num_topos(TOPO_FACE, numfaces));
-      MG_API_CALL(m, get_num_topos(TOPO_EDGE, numedges));
-      MG_API_CALL(m, get_num_topos(TOPO_VERTEX, numverts));
-      std::cout << "num regions is " << numregs << std::endl;
-      std::cout << "num faces   is " << numfaces << std::endl;
-      std::cout << "num edges   is " << numedges << std::endl;
-      std::cout << "num verts   is " << numverts << std::endl;
-      std::cout << "-----------" << std::endl;
-      if(numregs > 0)
-      {
-	  mmodel = ammodel;
-	  break;
-      }
-  }
-
-  /* SET THE ADJACENCIES */
-  MG_API_CALL(m, set_adjacency_state(REGION2FACE|
-                                     REGION2EDGE|
-                                     REGION2VERTEX|
-                                     FACE2EDGE|
-                                     FACE2VERTEX));
-  MG_API_CALL(m, set_reverse_states());
-  MG_API_CALL(m, set_adjacency_scope(TOPO_EDGE, SCOPE_FULL));
-  MG_API_CALL(m, set_adjacency_scope(TOPO_FACE, SCOPE_FULL));
-  MG_API_CALL(m, compute_adjacency());
+  lion_set_verbosity(1);
+  gmi_cap_start();
+  gmi_register_cap();
 
+  gmi_model* model = gmi_cap_load(creFileName);
+  apf::Mesh2* mesh = apf::createCapMesh(model, &PCUobj);
 
   // check parametrization using capstone apis
-  checkParametrization(m, g);
-
+  checkParametrization(mesh);
+  apf::destroyMesh(mesh);
+  gmi_cap_stop();
   } // pcu object scope
   pcu::Finalize();
 }
 
-void checkParametrization(MeshDatabaseInterface* mdb, GeometryDatabaseInterface* gdb)
-{
-   MeshSmartIterator miter(mdb);
-   mdb->get_topo_iterator(TOPO_VERTEX, miter);
-   int count = 0;
-   double sum = 0.0;
-   for(mdb->iterator_begin(miter); !mdb->iterator_end(miter); mdb->iterator_next(miter)) {
-     M_MTopo vert = mdb->iterator_value(miter);
-     M_GTopo geom;
-     GeometryTopoType gtype;
-     mdb->get_geom_entity(vert, gtype, geom);
-     if (!gdb->is_face(geom)) continue;
-     double range_u[2];
-     double range_v[2];
-     gdb->get_parametrization_range(geom, 0, range_u[0], range_u[1]);
-     gdb->get_parametrization_range(geom, 1, range_v[0], range_v[1]);
-     GeometryTopoType gtype1;
-     double u,v;
-     mdb->get_vertex_uv_parameters(vert, u, v, gtype1);
-     PCU_ALWAYS_ASSERT(gtype1 == gtype);
-
-     // coordinate from mesh
-     apf::Vector3 coord;
-     mdb->get_vertex_coord(vert, &(coord[0]));
-
-     // coordinate from surface
-     vec3d x;
-     gdb->get_point(geom, vec3d(u, v, 0.0), x);
-     apf::Vector3 pcoord(x[0], x[1], x[2]);
-
-     if (count < 50)
-       printf("%d, %e, %e, %e, %e, %e, %e, %e\n", count, u, v, range_u[0], range_u[1], range_v[0], range_v[1], (coord-pcoord).getLength());
-     sum += (coord-pcoord) * (coord-pcoord);
-     count++;
-   }
-   printf("norm of the difference vector is %e\n", std::sqrt(sum));
+void checkParametrization(apf::Mesh* mesh) {
+  int count[2] = {0, 0};
+  double sum = 0.0;
+  std::cout <<
+    "ct, (        u,         v), (     umin,      umax), "
+    "(     vmin,      vmax), diff\n";
+  apf::MeshIterator* it = mesh->begin(0);
+  for (apf::MeshEntity* vtx; (vtx = mesh->iterate(it));) {
+    apf::ModelEntity* me = mesh->toModel(vtx);
+    int dim = mesh->getModelType(me);
+    if (dim != 1 && dim != 2) continue;
+    double range_u[2], range_v[2];
+    mesh->getPeriodicRange(me, 0, range_u);
+    mesh->getPeriodicRange(me, 1, range_v);
+    apf::Vector3 xi;
+    mesh->getParam(vtx, xi);
+    // coordinate from mesh
+    apf::Vector3 coord;
+    mesh->getPoint(vtx, 0, coord);
+    // coordinate from surface
+    apf::Vector3 pcoord;
+    mesh->snapToModel(me, xi, pcoord);
+    apf::Vector3 diff = coord - pcoord;
+
+    constexpr int dim_printmax = 25;
+    if (count[dim - 1] < dim_printmax) {
+      std::cout << std::setw(2) << count[dim - 1] << ", "
+        << std::setprecision(3) << std::scientific
+        << '(' << std::setw(8) << xi.x() << ", "
+        << std::setw(8) << xi.y() << "), "
+        << '(' << std::setw(4) << range_u[0] << ", "
+        << std::setw(4) << range_u[1] << "), ";
+      if (dim == 2)
+        std::cout << '(' << std::setw(4) << range_v[0] << ", "
+          << std::setw(4) << range_v[1] << "), ";
+      else
+        std::cout << "----------------------, ";
+      std::cout << std::setw(4) << std::defaultfloat << diff.getLength()
+        << std::endl;
+    } else if (count[dim - 1] == dim_printmax) {
+      std::cout << "Skipping printing remaining "
+        << (dim == 1 ? "edges" : "faces") << std::endl;
+    }
+    sum += diff * diff;
+    count[dim - 1]++;
+  }
+  std::cout << "norm of the difference vector is " << std::sqrt(sum)
+    << std::endl;
 }
diff --git a/test/capGeomTest.cc b/test/capGeomTest.cc
index 25f681fa2..589f47fc4 100644
--- a/test/capGeomTest.cc
+++ b/test/capGeomTest.cc
@@ -1,36 +1,17 @@
+#include <iostream>
+#include <iomanip>
+
 #include <PCU.h>
+#include <apf.h>
 #include <apfCAP.h>
 #include <apfMDS.h>
+#include <apfMesh2.h>
 #include <gmi.h>
 #include <gmi_cap.h>
-#include <gmi_null.h>
 #include <gmi_mesh.h>
-#include <apf.h>
-#include <apfConvert.h>
-#include <apfMesh2.h>
-#include <apfNumbering.h>
-#include <ma.h>
+#include <gmi_null.h>
+#include <lionPrint.h>
 #include <pcu_util.h>
-#include <cstdlib>
-#include <cstring>
-#include <iostream>
-#include <iomanip>
-
-
-#include "CapstoneModule.h"
-#include "CreateMG_Framework_Core.h"
-#include "CreateMG_Framework_Analysis.h"
-#include "CreateMG_Framework_Application.h"
-#include "CreateMG_Framework_Attributes.h"
-#include "CreateMG_Framework_Core.h"
-#include "CreateMG_Framework_Geometry.h"
-#include "CreateMG_Framework_Mesh.h"
-
-using namespace CreateMG;
-using namespace CreateMG::Attribution;
-using namespace CreateMG::Mesh;
-using namespace CreateMG::Geometry;
-
 
 char const* const typeName[4] =
 {"vertex",
@@ -44,85 +25,24 @@ void visualizeFace(gmi_model* model, gmi_ent* entity, int n, int m, const char*
 void visualizeEdge(gmi_model* model, gmi_ent* entity, int n, const char* fileName, pcu::PCU *PCUObj);
 void visualizeEdges(gmi_model* model, int n, const char* fileName, pcu::PCU *PCUObj);
 
-
-
 int main(int argc, char** argv)
 {
   pcu::Init(&argc, &argv);
   { // pcu object scope
   pcu::PCU PCUObj;
 
+  PCU_ALWAYS_ASSERT(argc == 2);
+
   gmi_register_mesh();
   gmi_register_null();
-
-  // load capstone mesh
-  // create an instance of the Capstone Module activating CREATE/CREATE/CREATE
-  // for the Geometry/Mesh/Attribution databases
-  /* const std::string gdbName("Geometry Database : Create");// Switch Create with SMLIB for CAD */
-  const std::string gdbName("Geometry Database : SMLIB");// Switch Create with SMLIB for CAD
-  const std::string mdbName("Mesh Database : Create");
-  const std::string adbName("Attribution Database : Create");
-
-  CapstoneModule  cs("the_module", gdbName.c_str(), mdbName.c_str(), adbName.c_str());
-
-  GeometryDatabaseInterface     *g = cs.get_geometry();
-  MeshDatabaseInterface         *m = cs.get_mesh();
-  AppContext                    *c = cs.get_context();
-
-
-  PCU_ALWAYS_ASSERT(g);
-  PCU_ALWAYS_ASSERT(m);
-  PCU_ALWAYS_ASSERT(c);
-
-  v_string filenames;
-  filenames.push_back(argv[1]);
-
-  M_GModel gmodel = cs.load_files(filenames);
-
-  int numbreps = 0;
-  MG_CALL(g->get_num_breps(numbreps));
-  std::cout << "number of b reps is " << numbreps << std::endl;
-  if(numbreps == 0)
-      error(HERE, ERR_INVALID_INPUT, "Model is empty");
-
-  M_MModel mmodel;
-  // Pick the volume mesh model from associated mesh models to this geom model
-  std::vector<M_MModel> mmodels;
-  MG_API_CALL(m, get_associated_mesh_models(gmodel, mmodels));
-  for(std::size_t i = 0; i < mmodels.size(); ++i)
-  {
-      M_MModel ammodel = mmodels[i];
-      std::size_t numregs = 0;
-      std::size_t numfaces = 0;
-      std::size_t numedges = 0;
-      std::size_t numverts = 0;
-      MG_API_CALL(m, set_current_model(ammodel));
-      MG_API_CALL(m, get_num_topos(TOPO_REGION, numregs));
-      MG_API_CALL(m, get_num_topos(TOPO_FACE, numfaces));
-      MG_API_CALL(m, get_num_topos(TOPO_EDGE, numedges));
-      MG_API_CALL(m, get_num_topos(TOPO_VERTEX, numverts));
-      std::cout << "num regions is " << numregs << std::endl;
-      std::cout << "num faces   is " << numfaces << std::endl;
-      std::cout << "num edges   is " << numedges << std::endl;
-      std::cout << "num verts   is " << numverts << std::endl;
-      std::cout << "-----------" << std::endl;
-      if(numregs > 0)
-      {
-	  mmodel = ammodel;
-	  break;
-      }
-  }
-
-
   gmi_cap_start();
   gmi_register_cap();
 
-
-  apf::Mesh2* mesh0 = apf::createMesh(m,g,&PCUObj);
+  gmi_model* model = gmi_load(argv[1]);
+  apf::Mesh2* mesh0 = apf::createCapMesh(model, &PCUObj);
+  apf::disownCapModel(mesh0);
   apf::writeVtkFiles("mesh_no_param", mesh0);
 
-  gmi_model* model = gmi_import_cap(g);
-
   printInfo(model, 0);
   printf("\n");
   printInfo(model, 1);
@@ -147,7 +67,7 @@ int main(int argc, char** argv)
   printf("creating mesh with param field\n");
 
 
-  apf::Mesh2* mesh = apf::createMesh(m,g,&PCUObj);
+  apf::Mesh2* mesh = apf::createCapMesh(model, &PCUObj);
   apf::Field* pf  = apf::createFieldOn(mesh, "param_field", apf::VECTOR);
   apf::Field* idf  = apf::createFieldOn(mesh, "id", apf::SCALAR);
   apf::MeshEntity* e;
@@ -255,8 +175,8 @@ void visualizeFace(gmi_model* model, gmi_ent* entity, int n, int m, const char*
   double dv = (v_range[1] - v_range[0]) / (m-1);
 
   // make the array of vertex coordinates in the physical space
-  std::vector<ma::Vector> ps;
-  std::vector<ma::Vector> uvs;
+  std::vector<apf::Vector3> ps;
+  std::vector<apf::Vector3> uvs;
   for (int j = 0; j < m; j++) {
     for (int i = 0; i < n; i++) {
       double params[2];
@@ -264,18 +184,18 @@ void visualizeFace(gmi_model* model, gmi_ent* entity, int n, int m, const char*
       params[1] = v_range[0] + j * dv;
       double position[3];
       gmi_eval(model, entity, &params[0], &position[0]);
-      ma::Vector p(position[0], position[1], position[2]);
+      apf::Vector3 p(position[0], position[1], position[2]);
       ps.push_back(p);
-      ma::Vector uv(params[0], params[1], 0.);
+      apf::Vector3 uv(params[0], params[1], 0.);
       uvs.push_back(uv);
     }
   }
 
   // make the vertexes and set the coordinates using the array
-  std::vector<ma::Entity*> vs;
+  std::vector<apf::MeshEntity*> vs;
   apf::Mesh2* mesh = apf::makeEmptyMdsMesh(gmi_load(".null"), 2, false, PCUObj);
   for (size_t i = 0; i < ps.size(); i++) {
-    ma::Entity* vert = mesh->createVert(0);
+    apf::MeshEntity* vert = mesh->createVert(0);
     mesh->setPoint(vert, 0, ps[i]);
     vs.push_back(vert);
   }
@@ -283,7 +203,7 @@ void visualizeFace(gmi_model* model, gmi_ent* entity, int n, int m, const char*
   assert(vs.size() == ps.size());
 
 
-  ma::Entity* v[3];
+  apf::MeshEntity* v[3];
   // make the lower/upper t elems
   for (int i = 0; i < n-1; i++) {
     for (int j = 0; j < m-1; j++) {
@@ -307,8 +227,8 @@ void visualizeFace(gmi_model* model, gmi_ent* entity, int n, int m, const char*
 
   // parametric coordinates
   apf::Field* f = apf::createFieldOn(mesh, "param_coords", apf::VECTOR);
-  ma::Entity* e;
-  ma::Iterator* it;
+  apf::MeshEntity* e;
+  apf::MeshIterator* it;
   it = mesh->begin(0);
   int count = 0;
   while ( (e = mesh->iterate(it)) ) {
@@ -330,9 +250,9 @@ void visualizeFace(gmi_model* model, gmi_ent* entity, int n, int m, const char*
     double param[2] = {uvs[count].x(), uvs[count].y()};
     gmi_first_derivative(model, entity, param, uTangent, vTangent);
     gmi_normal(model, entity, param, normal);
-    apf::setVector(ut, e, 0, ma::Vector(uTangent[0], uTangent[1], uTangent[2]));
-    apf::setVector(vt, e, 0, ma::Vector(vTangent[0], vTangent[1], vTangent[2]));
-    apf::setVector(nv, e, 0, ma::Vector(normal[0], normal[1], normal[2]));
+    apf::setVector(ut, e, 0, apf::Vector3(uTangent[0], uTangent[1], uTangent[2]));
+    apf::setVector(vt, e, 0, apf::Vector3(vTangent[0], vTangent[1], vTangent[2]));
+    apf::setVector(nv, e, 0, apf::Vector3(normal[0], normal[1], normal[2]));
     count++;
   }
 
@@ -353,24 +273,24 @@ void visualizeEdge(gmi_model* model, gmi_ent* entity, int n, const char* fileNam
   double du = (u_range[1] - u_range[0]) / (n-1);
 
   // make the array of vertex coordinates in the physical space
-  std::vector<ma::Vector> ps;
-  std::vector<ma::Vector> us;
+  std::vector<apf::Vector3> ps;
+  std::vector<apf::Vector3> us;
   for (int i = 0; i < n; i++) {
     double params[2];
     params[0] = u_range[0] + i * du;
     double position[3];
     gmi_eval(model, entity, &params[0], &position[0]);
-    ma::Vector p(position[0], position[1], position[2]);
+    apf::Vector3 p(position[0], position[1], position[2]);
     ps.push_back(p);
-    ma::Vector uv(params[0], params[1], 0.);
+    apf::Vector3 uv(params[0], params[1], 0.);
     us.push_back(uv);
   }
 
   // make the vertexes and set the coordinates using the array
-  std::vector<ma::Entity*> vs;
+  std::vector<apf::MeshEntity*> vs;
   apf::Mesh2* mesh = apf::makeEmptyMdsMesh(gmi_load(".null"), 1, false, PCUObj);
   for (size_t i = 0; i < ps.size(); i++) {
-    ma::Entity* vert = mesh->createVert(0);
+    apf::MeshEntity* vert = mesh->createVert(0);
     mesh->setPoint(vert, 0, ps[i]);
     vs.push_back(vert);
   }
@@ -378,7 +298,7 @@ void visualizeEdge(gmi_model* model, gmi_ent* entity, int n, const char* fileNam
   assert(vs.size() == ps.size());
 
 
-  ma::Entity* v[2];
+  apf::MeshEntity* v[2];
   // make the lower/upper t elems
   for (int i = 0; i < n-1; i++) {
     v[0] = vs[i];
@@ -392,8 +312,8 @@ void visualizeEdge(gmi_model* model, gmi_ent* entity, int n, const char* fileNam
   apf::printStats(mesh);
 
   apf::Field* f = apf::createFieldOn(mesh, "param_coords", apf::VECTOR);
-  ma::Entity* e;
-  ma::Iterator* it;
+  apf::MeshEntity* e;
+  apf::MeshIterator* it;
   it = mesh->begin(0);
   int count = 0;
   while ( (e = mesh->iterate(it)) ) {
@@ -420,30 +340,30 @@ void visualizeEdges(gmi_model* model, int n, const char* fileName, pcu::PCU *PCU
     double du = (u_range[1] - u_range[0]) / (n-1);
 
     // make the array of vertex coordinates in the physical space
-    std::vector<ma::Vector> ps;
-    std::vector<ma::Vector> us;
+    std::vector<apf::Vector3> ps;
+    std::vector<apf::Vector3> us;
     for (int i = 0; i < n; i++) {
       double params[2];
       params[0] = u_range[0] + i * du;
       double position[3];
       gmi_eval(model, entity, &params[0], &position[0]);
-      ma::Vector p(position[0], position[1], position[2]);
+      apf::Vector3 p(position[0], position[1], position[2]);
       ps.push_back(p);
-      ma::Vector uv(params[0], params[1], 0.);
+      apf::Vector3 uv(params[0], params[1], 0.);
       us.push_back(uv);
     }
 
     // make the vertexes and set the coordinates using the array
-    std::vector<ma::Entity*> vs;
+    std::vector<apf::MeshEntity*> vs;
     for (size_t i = 0; i < ps.size(); i++) {
-      ma::Entity* vert = mesh->createVert(0);
+      apf::MeshEntity* vert = mesh->createVert(0);
       mesh->setPoint(vert, 0, ps[i]);
       vs.push_back(vert);
     }
 
     assert(vs.size() == ps.size());
 
-    ma::Entity* v[2];
+    apf::MeshEntity* v[2];
     // make the lower/upper t elems
     for (int i = 0; i < n-1; i++) {
       v[0] = vs[i];
diff --git a/test/capLoadSome.cc b/test/capLoadSome.cc
new file mode 100644
index 000000000..11ae8637a
--- /dev/null
+++ b/test/capLoadSome.cc
@@ -0,0 +1,26 @@
+#include <iostream>
+#include <string>
+#include <vector>
+
+#include <PCU.h>
+#include <gmi.h>
+#include <gmi_cap.h>
+
+int main(int argc, char* argv[]) {
+  if (argc < 2) {
+    std::cerr << "USAGE: " << argv[0] << " <input.cre> [meshname...]"
+      << std::endl;
+    return 1;
+  }
+  pcu::Init(&argc, &argv);
+  gmi_cap_start();
+  gmi_register_cap();
+  std::string model_content;
+  std::vector<std::string> mesh_names;
+  for (int i = 2; i < argc; ++i) mesh_names.push_back(argv[i]);
+  gmi_model* model = gmi_cap_load_selective(argv[1], mesh_names);
+  gmi_destroy(model);
+  gmi_cap_stop();
+  pcu::Finalize();
+  return 0;
+}
diff --git a/test/capNativeAdapt.cc b/test/capNativeAdapt.cc
new file mode 100644
index 000000000..a9afc9bad
--- /dev/null
+++ b/test/capNativeAdapt.cc
@@ -0,0 +1,114 @@
+#include <exception>
+#include <string>
+#include <vector>
+
+#include <PCU.h>
+#include <apf.h>
+#include <apfCAP.h>
+#include <apfMesh2.h>
+#include <lionPrint.h>
+#include <gmi_cap.h>
+#include <pcu_util.h>
+
+#include <CreateMG_AppProcessor.h>
+#include <CreateMG_Framework_Mesh.h>
+#include <CreateMG_Function.h>
+#include <CreateMG_SizingMetricTool.h>
+
+namespace {
+
+/** \brief Print nested exceptions. */
+void print_exception(const std::exception& e, int level = 0);
+
+} // namespace
+
+int main(int argc, char* argv[]) {
+  lion_set_verbosity(1);
+  pcu::Init(&argc, &argv);
+  try {
+    pcu::PCU PCU;
+    if (argc < 4) {
+      std::cerr << "USAGE: " << argv[0] << " <in.cre> <out.cre> <sizing.dat>"
+        " <in.vmap>" << std::endl;
+      throw 1;
+    }
+
+    const char *creFile = argv[1], *outFile = argv[2], *sizingFile = argv[3],
+      *vmapFile = argv[4];
+
+    gmi_cap_start();
+    gmi_register_cap();
+    try {
+      gmi_model* model = gmi_cap_load(creFile);
+      apf::Mesh2* mesh = apf::createCapMesh(model, &PCU);
+      try {
+        auto mdbi = apf::exportCapNative(mesh);
+        auto ctx = mdbi->get_context();
+        auto proc = CreateMG::get_context_processor(ctx);
+        auto fn = CreateMG::get_function(ctx, "AdaptMesh");
+        CreateMG::M_MModel mmodel;
+        MG_API_CALL(mdbi, get_current_model(mmodel));
+        CreateMG::set_input(fn, "MeshModel", mmodel);
+        // CreateMG::set_input(fn, "Analysis", mmodel);
+        auto sTool = CreateMG::get_sizing_metric_tool(
+          ctx, "CreateSmoothingBase"
+        );
+        PCU_ALWAYS_ASSERT(sTool);
+        sTool->set_context(ctx);
+        std::vector<CreateMG::Metric6> sizing6;
+        apf::loadCapSizingFileMetrics(mesh, sizingFile, vmapFile, sizing6);
+        sTool->set_metric(mmodel, "sizing6", sizing6);
+        CreateMG::set_input(fn, "TensorData", "sizing6");
+        // DiscreteCurvature, Complexity, MaxComplexity
+        auto t0 = pcu::Time();
+        if (proc->execute(fn) != CreateMG::STATUS_OK) {
+          std::cerr << "ERROR: failed to adapt mesh" << std::endl;
+          throw 1;
+        }
+        auto t1 = pcu::Time();
+        std::cout << "Capstone AdaptMesh ran in " << t1 - t0 << std::endl;
+        CreateMG::get_output(fn, "Output", mmodel);
+        std::string oname;
+        MG_API_CALL(mdbi, get_model_name(mmodel, oname));
+        apf::Mesh2* omesh = apf::createCapMesh(model, oname.c_str(), &PCU);
+        apf::disownCapModel(omesh);
+        apf::writeVtkFiles((std::string(outFile) + ".vtk").c_str(), omesh);
+        gmi_cap_write(model, outFile);
+        apf::destroyMesh(omesh);
+      } catch (...) {
+        apf::destroyMesh(mesh);
+        std::rethrow_exception(std::current_exception());
+      }
+      apf::destroyMesh(mesh);
+    } catch (const std::exception& e) {
+      if (PCU.Self() == 0) {
+        std::cerr << "ERROR: ";
+        print_exception(e);
+      }
+      gmi_cap_stop();
+      throw 1;
+    } catch (...) {
+      gmi_cap_stop();
+      throw 1;
+    }
+    gmi_cap_stop();
+  } catch (...) {
+    pcu::Finalize();
+    return 1;
+  }
+  pcu::Finalize();
+  return 0;
+}
+
+namespace {
+
+void print_exception(const std::exception& e, int level) {
+  std::cerr << std::string(level * 2, ' ') << e.what() << '\n';
+  try {
+    std::rethrow_if_nested(e);
+  } catch (const std::exception& nestedE) {
+    print_exception(nestedE, level + 1);
+  } catch (...) {}
+}
+
+} // namespace
diff --git a/test/capProbe.cc b/test/capProbe.cc
new file mode 100644
index 000000000..fb8e731a7
--- /dev/null
+++ b/test/capProbe.cc
@@ -0,0 +1,31 @@
+#include <iostream>
+#include <string>
+#include <vector>
+
+#include <apf.h>
+#include <apfMesh.h>
+#include <apfCAP.h>
+#include <gmi.h>
+#include <gmi_cap.h>
+#include <PCU.h>
+
+int main(int argc, char* argv[]) {
+  if (argc != 2) {
+    std::cerr << "USAGE: " << argv[0] << " <input.cre>" << std::endl;
+    return 1;
+  }
+  pcu::Init(&argc, &argv);
+  gmi_cap_start();
+  gmi_register_cap();
+  std::string model_content;
+  std::vector<std::string> mesh_names, mesh_contents;
+  gmi_cap_probe(argv[1], model_content, mesh_names, mesh_contents);
+  std::cout << "model: " << model_content << std::endl;
+  for (std::size_t i = 0; i < mesh_names.size(); ++i) {
+    std::cout << "mesh (" << i << ") `" << mesh_names[i] << "`: "
+      << mesh_contents[i] << std::endl;
+  }
+  gmi_cap_stop();
+  pcu::Finalize();
+  return 0;
+}
diff --git a/test/capVol.cc b/test/capVol.cc
index 7e44c2e36..a6d0fc897 100644
--- a/test/capVol.cc
+++ b/test/capVol.cc
@@ -1,276 +1,353 @@
-#include <cstring>
-#include <cstdlib>
+#include <exception>
+#include <stdexcept>
 
-// Output
-#include <lionPrint.h>
-
-// Parallelism
 #include <PCU.h>
-#include <pcu_util.h>
-
-// Mesh interfaces
 #include <apf.h>
 #include <apfCAP.h>
-
-// Geometry interfaces
+#include <apfConvert.h>
+#include <apfMDS.h>
+#include <apfMETIS.h>
+#include <apfMesh2.h>
+#include <apfPartition.h>
+#include <apfZoltan.h>
 #include <gmi.h>
 #include <gmi_cap.h>
-
-// Mesh adapt
+#include <lionPrint.h>
 #include <ma.h>
-
-using namespace CreateMG;
-using namespace CreateMG::Attribution;
-using namespace CreateMG::Mesh;
-using namespace CreateMG::Geometry;
+#include <parma.h>
+#include <pcu_util.h>
 
 #include "capVolSizeFields.h"
 
 namespace {
 
-void myExit(int exit_code = EXIT_SUCCESS) {
-  gmi_cap_stop();
-  pcu::Finalize();
-  exit(exit_code);
+/** \brief Print nested exceptions. */
+void print_exception(
+  const pcu::PCU& PCU, const std::exception& e, int level = 0
+);
+
+struct Args {
+  Args() {}
+  void parse(int argc, char* argv[]);
+  static void print_usage(const char* argv0);
+  std::string before, after, in, out;
+  bool analytic{false}, volume{false}, mds{false};
+  int sf{0};
+  enum Partitioner { Default, Zoltan, METIS, Parma } splitter, balancer;
+};
+
+void Args::print_usage(const char *argv0) {
+  std::cout << "USAGE: " << argv0 << " [-agm] [-B BEFORE.VTK] [-A AFTER.VTK] "
+    "[-s SPLITTER] [-b BALANCER] <size-field> <IN.CRE> <OUT.CRE>" << std::endl;
+  std::cout << R"help(
+Flags:
+-B BEFORE.VTK  Write the file BEFORE.VTK before adaptation with the
+               adapt_frames and adapt_scales fields.
+-A AFTER.VTK   Write the file AFTER.VTK after adaptation.
+-a             Evaluate size-field analytically during adaptation. The default
+               is to evaluate once, write the frames/scales, then transfer and
+               interpolate during adaptation.
+-g             Force mesh volume generation.
+-m             Convert mesh to MDS during adaptation (required for parallel
+               adaptation).
+-s SPLITTER    Force the selected splitter to be used. Possible values are:
+               Zoltan, METIS, Parma. The default is to use whatever is
+               available, in the order listed previously.
+-b BALANCER    Force the selected balancer to be used during adaptation.
+               Possible values and default priorities are the same as for
+               splitter selection. Using this option forces load-balancing
+               after each step, whereas the default only balances when the
+               estimated imbalance exceeds a threshold (10%%).
+
+SIZE-FIELDS:
+1, for uniform anisotropic size-field
+2, for wing-shock size-field
+3, for cube-shock size-field
+4, for cylinder boundary-layer size-field
+)help";
 }
 
-void writeCre(CapstoneModule& cs, const std::string& filename) {
-  GeometryDatabaseInterface    *gdbi = cs.get_geometry();
-  MeshDatabaseInterface        *mdbi = cs.get_mesh();
-  AppContext		       *ctx = cs.get_context();
-
-  // Get the CRE writer.
-  Writer *creWriter = get_writer(ctx, "Create Native Writer");
-  if (!creWriter) {
-    lion_eprint(1, "FATAL: Could not find the CRE writer.\n");
-    myExit(EXIT_FAILURE);
-  }
+ma::Mesh* loadAdaptMesh(
+  pcu::PCU* PCU, gmi_model* model, bool volume, bool mds
+);
 
-  IdMapping idmapping;
-  std::vector<M_MModel> mmodels;
-  M_GModel gmodel;
-  M_MModel mmodel;
-  gdbi->get_current_model(gmodel);
-  mdbi->get_current_model(mmodel);
-  mmodels.clear();
-  mmodels.push_back(mmodel);
-  creWriter->write(ctx, gmodel, mmodels, filename.c_str(), idmapping);
-}
-
-void printUsage(char *argv0) {
-  printf("USAGE: %s [-agwv] <size-field> <create_file.cre>\n", argv0);
-  printf("Flags:\n"
-  "-a\tEvaluate size-field analytically.\n"
-  "-g\tForce mesh generation.\n"
-  "-v\tEnable verbose output.\n"
-  "-w\tWrite before.vtk, after.vtk, and after.cre.\n"
-  "SIZE-FIELDS:\n"
-  "%d, for uniform anisotropic size-field\n"
-  "%d, for wing-shock size-field\n"
-  "%d, for cube-shock size-field\n"
-  "%d, for cylinder boundary-layer size-field\n", 1, 2, 3, 4);
-}
+void parallelAdapt(
+  pcu::PCU* PCU, gmi_model* model, apf::Mesh2* mesh, const Args& args
+);
 
 } // namespace
 
 int main(int argc, char** argv) {
+  lion_set_verbosity(1); // Initialize logging.
+  int retval = 0;
   // Initialize parallelism.
   pcu::Init(&argc, &argv);
   {
   pcu::PCU PCUObj;
+  try {
+    // Check arguments or print usage.
+    Args args;
+    try {
+      args.parse(argc, argv);
+    } catch (const std::exception& e) {
+      if (PCUObj.Self() == 0) args.print_usage(argv[0]);
+      std::throw_with_nested(
+        std::runtime_error("invalid command line arguments")
+      );
+    }
 
-  // Initialize logging.
-  lion_set_stdout(stdout);
-  lion_set_stderr(stderr);
-
-  // Check arguments or print usage.
-  if (argc < 3) {
-    if (PCUObj.Self() == 0) {
-      printUsage(argv[0]);
+    if (PCUObj.Peers() > 1 && !args.mds)
+      throw std::runtime_error("parallel run without -m flag");
+
+    // Initialize GMI.
+    gmi_cap_start();
+    gmi_register_cap();
+
+    try {
+      gmi_model* capGeomModel = nullptr;
+      if (PCUObj.Self() == 0) capGeomModel = gmi_cap_load(args.in.c_str());
+      else capGeomModel = gmi_cap_load_selective(args.in.c_str(), {});
+      auto soloPCU = PCUObj.Split(PCUObj.Self(), 0);
+      ma::Mesh* mesh = nullptr;
+      if (PCUObj.Self() == 0) {
+        mesh = loadAdaptMesh(
+          soloPCU.get(), capGeomModel, args.volume, args.mds
+        );
+        // APF default routine will typically fail to verify surface meshes.
+        if (args.volume) mesh->verify();
+      }
+      parallelAdapt(&PCUObj, capGeomModel, mesh, args);
+      if (PCUObj.Self() == 0) {
+        if (args.volume) mesh->verify();
+        if (args.mds) {
+          apf::Mesh2* mdsMesh = mesh;
+          mesh = apf::makeEmptyCapMesh(
+            capGeomModel, "MeshAdapt", soloPCU.get()
+          );
+          apf::disownCapModel(mesh);
+          apf::convert(mdsMesh, mesh);
+          apf::destroyMesh(mdsMesh);
+        }
+        gmi_cap_write(capGeomModel, args.out.c_str());
+        apf::destroyMesh(mesh);
+      }
+      gmi_destroy(capGeomModel);
+      gmi_cap_stop();
+    } catch(...) {
+      gmi_cap_stop();
+      std::rethrow_exception(std::current_exception());
     }
-    myExit(EXIT_FAILURE);
+  } catch (const std::exception& e) {
+    if (PCUObj.Self() == 0) std::cerr << "ERROR: ";
+    print_exception(PCUObj, e);
+    retval = 1;
+  } catch(...) {
+    if (PCUObj.Self() == 0)
+      std::cerr << "ERROR: unspecified error" << std::endl;
+    retval = 1;
   }
+  } // PCUObj scope
+  pcu::Finalize();
+  return retval;
+}
 
-  // Parse arguments.
-  bool volume_flag = false, write_flag = false, analytic_flag = false,
-       verbose_flag = false;
-  for (int i = 1; i < argc - 2; ++i) {
+namespace {
+
+void print_exception(const pcu::PCU& PCU, const std::exception& e, int level) {
+  if (PCU.Self() == 0)
+    std::cerr << std::string(level * 2, ' ') << e.what() << '\n';
+  try {
+    std::rethrow_if_nested(e);
+  } catch (const std::exception& nestedE) {
+    print_exception(PCU, nestedE, level + 1);
+  } catch (...) {}
+}
+
+void Args::parse(int argc, char* argv[]) {
+  constexpr int positional_args = 3;
+  if (argc < 1 + positional_args)
+    throw std::runtime_error("missing positional argument(s)");
+  auto strarg = [argc, argv](int& i, int& j) -> std::string {
+    if (argv[i][j + 1]) {
+      int oldj = j;
+      while (argv[i][j + 1] != '\0') ++j;
+      return &argv[i][oldj + 1];
+    } else if (i + 1 < argc - positional_args) {
+      ++i;
+      for (j = 0; argv[i][j + 1] != '\0'; ++j);
+      return argv[i];
+    } else throw std::runtime_error("missing argument for -" + argv[i][j]);
+  };
+  auto strlower = [](std::string str) {
+    for (auto& c : str) c = std::tolower(c);
+    return str;
+  };
+  auto str2part = [strlower](std::string str) -> Partitioner {
+    str = strlower(str);
+    if (str == "metis") return METIS;
+    if (str == "parma") return Parma;
+    if (str == "zoltan") return Zoltan;
+    throw std::invalid_argument("invalid partitioner name: " + str);
+  };
+  int i;
+  for (i = 1; i < argc - positional_args; ++i) {
     if (*argv[i] == '-') {
       for (int j = 1; argv[i][j] != '\0'; ++j) {
         switch(argv[i][j]) {
-        case 'a':
-          analytic_flag = true;
-          break;
-        case 'g':
-          volume_flag = true;
-          break;
-        case 'v':
-          verbose_flag = true;
-          lion_set_verbosity(1);
-          break;
-        case 'w':
-          write_flag = true;
-          break;
-        default:
-          printf("Error: invalid flag.\n");
-          printUsage(argv[0]);
-          myExit(EXIT_FAILURE);
+        case 'A': after = strarg(i, j); break;
+        case 'B': before = strarg(i, j); break;
+        case 'a': analytic = true; break;
+        case 'g': volume = true; break;
+        case 'm': mds = true; break;
+        case 's': splitter = str2part(strarg(i, j)); break;
+        case 'b': balancer = str2part(strarg(i, j)); break;
+        default: throw std::runtime_error("unrecognized flag: -" + argv[i][j]);
         }
       }
-    }
+    } else break;
   }
+  sf = std::atoi(argv[i]);
+  in = argv[i + 1];
+  out = argv[i + 2];
+}
 
-  const char* createFileName = argv[argc - 1];
-  int mode = atoi(argv[argc - 2]);
-
-  // Initialize GMI.
-  gmi_cap_start();
-  gmi_register_cap();
-  // create an instance of the Capstone Module activating SMLIB/CREATE/CREATE
-  // for the Geometry/Mesh/Attribution databases
-  const std::string gdbName("Geometry Database : SMLIB");
-  const std::string mdbName("Mesh Database : Create");
-  const std::string adbName("Attribution Database : Create");
-
-  CapstoneModule  cs("capTest", gdbName.c_str(), mdbName.c_str(), adbName.c_str());
-
-  GeometryDatabaseInterface     *g = cs.get_geometry();
-  MeshDatabaseInterface         *m = cs.get_mesh();
-  AppContext                    *c = cs.get_context();
-
-  PCU_ALWAYS_ASSERT(g);
-  PCU_ALWAYS_ASSERT(m);
-  PCU_ALWAYS_ASSERT(c);
-
-  // Load Capstone mesh.
-  v_string filenames;
-  filenames.push_back(createFileName);
-  M_GModel gmodel = cs.load_files(filenames);
+ma::Mesh* loadAdaptMesh(
+  pcu::PCU* PCU, gmi_model* model, bool volume, bool mds
+) {
+  // Load Capstone mesh (with optional volume generation).
+  ma::Mesh* capMesh = nullptr;
+  if (volume) {
+    int dim = 3;
+    capMesh = apf::generateCapMesh(model, dim, PCU);
+    // FIXME: maybe create copy of mesh model to work on (preserve original).
+  } else capMesh = apf::createCapMesh(model, PCU);
+  apf::disownCapModel(capMesh);
+  // Optionally convert to MDS
+  ma::Mesh* adaptMesh = nullptr;
+  if (mds) {
+    adaptMesh = apf::createMdsMesh(model, capMesh, true);
+    apf::disownMdsModel(adaptMesh); // Model is managed in main.
+    apf::destroyMesh(capMesh);
+  } else adaptMesh = capMesh;
+  return adaptMesh;
+}
 
-  if (volume_flag) {
-    M_MModel mmodel = cs.generate_mesh();
-    if (mmodel.is_invalid()) {
-      lion_eprint(1, "FATAL: Failed to mesh the model.\n");
-      myExit(EXIT_FAILURE);
-    }
-    MG_API_CALL(m, set_current_model(mmodel));
-  } else {
-    // Use the first existing mesh model.
-    std::vector<M_MModel> mmodels;
-    MG_API_CALL(m, get_associated_mesh_models(gmodel, mmodels));
-    PCU_ALWAYS_ASSERT(mmodels.size() == 1);
-    MG_API_CALL(m, set_current_model(mmodels[0]));
+apf::Splitter* makeSplitter(Args::Partitioner ptnr, apf::Mesh2* mesh) {
+  switch (ptnr) {
+  case Args::Zoltan:
+    return apf::makeZoltanSplitter(mesh, apf::GRAPH, apf::PARTITION);
+  case Args::METIS: return apf::makeMETISsplitter(mesh);
+  case Args::Parma: return Parma_MakeRibSplitter(mesh);
+  default:
+    #if defined(PUMI_HAS_ZOLTAN)
+    return apf::makeZoltanSplitter(
+      adaptMesh, apf::GRAPH, apf::PARTITION
+    );
+    #elif defined(PUMI_HAS_METIS)
+    return apf::makeMETISsplitter(mesh);
+    #else
+    return Parma_MakeRibSplitter(mesh);
+    #endif
   }
+}
 
-  if (write_flag) {
-    writeCre(cs, "core_capVol_before.cre");
+ma::AnisotropicFunction* makeUDF(int mode, ma::Mesh* mesh) {
+  switch (mode) {
+    case 1: return new UniformAniso(mesh);
+    case 2: return new WingShock(mesh, 50);
+    case 3: return new Shock2(mesh);
+    case 4: return new CylBoundaryLayer(mesh);
+    default: throw std::runtime_error("invalid size-field");
   }
+}
 
-  // Calculate adjacencies.
-  MG_API_CALL(m, set_adjacency_state(REGION2FACE|REGION2EDGE|REGION2VERTEX|
-				     FACE2EDGE|FACE2VERTEX));
-  MG_API_CALL(m, set_reverse_states());
-  MG_API_CALL(m, compute_adjacency());
-
-  // Make APF adapter over Capstone mesh.
-  ma::Mesh* apfCapMesh = apf::createMesh(m, g, &PCUObj);
+void migrateHome(apf::Mesh2* mesh) {
+  auto t0 = pcu::Time();
+  apf::Migration* plan = new apf::Migration(mesh);
+  apf::MeshIterator* it = mesh->begin(mesh->getDimension());
+  for (apf::MeshEntity* e; (e = mesh->iterate(it));) plan->send(e, 0);
+  mesh->end(it);
+  mesh->migrate(plan); // destroys plan
+  auto map0 = apf::Multiply(0);
+  apf::remapPartition(mesh, map0);
+  auto t1 = pcu::Time();
+  if (mesh->getPCU()->Self() == 0)
+    std::cout << "INFO: Migrated home in " << t1 - t0 << " seconds"
+      << std::endl;
+}
 
-  // Choose appropriate size-field.
-  ma::AnisotropicFunction* sf = nullptr;
-  switch (mode) {
-    case 1:
-      sf = new UniformAniso(apfCapMesh);
-      break;
-    case 2:
-      sf = new WingShock(apfCapMesh, 50);
-      break;
-    case 3:
-      sf = new Shock(apfCapMesh);
-      break;
-    case 4:
-      sf = new CylBoundaryLayer(apfCapMesh);
-      break;
-    default:
-      lion_eprint(1, "FATAL: Invalid size-field.\n");
-      myExit(EXIT_FAILURE);
+void parallelAdapt(
+  pcu::PCU* PCU, gmi_model* model, apf::Mesh2* mesh, const Args& args
+) {
+  pcu::PCU* oldPCU = nullptr;
+  if (mesh) oldPCU = mesh->getPCU();
+  if (PCU->Peers() > 1) {
+    apf::Migration* plan = nullptr;
+    if (PCU->Self() == 0) {
+      apf::Splitter* splitter = makeSplitter(args.splitter, mesh);
+      apf::MeshTag* weights = Parma_WeighByMemory(mesh);
+      plan = splitter->split(weights, 1.10, PCU->Peers());
+      apf::removeTagFromDimension(
+        mesh, weights, mesh->getDimension()
+      );
+      mesh->destroyTag(weights);
+      delete splitter;
+      mesh->switchPCU(PCU);
+    }
+    mesh = apf::repeatMdsMesh(mesh, model, plan, PCU->Peers(), PCU);
+    apf::disownMdsModel(mesh);
   }
+  // Choose appropriate size-field.
+  std::unique_ptr<ma::AnisotropicFunction> sf(makeUDF(args.sf, mesh));
 
   // Make pumi fields for the frames and scales for anisotropic size-fields.
   apf::Field* frameField = nullptr;
   apf::Field* scaleField = nullptr;
-  ma::Input *in = nullptr;
-  if (!analytic_flag || write_flag) {
-    frameField = apf::createFieldOn(apfCapMesh, "adapt_frames", apf::MATRIX);
-    scaleField = apf::createFieldOn(apfCapMesh, "adapt_scales", apf::VECTOR);
-
-    ma::Entity *v;
-    ma::Iterator* it = apfCapMesh->begin(0);
-    while( (v = apfCapMesh->iterate(it)) ) {
+  if (!args.analytic || !args.before.empty()) {
+    frameField = apf::createFieldOn(mesh, "adapt_frames", apf::MATRIX);
+    scaleField = apf::createFieldOn(mesh, "adapt_scales", apf::VECTOR);
+    apf::MeshIterator* it = mesh->begin(0);
+    for (apf::MeshEntity* v; (v = mesh->iterate(it));) {
       ma::Vector s;
       ma::Matrix f;
       sf->getValue(v, f, s);
       apf::setVector(scaleField, v, 0, s);
       apf::setMatrix(frameField, v, 0, f);
     }
-    apfCapMesh->end(it);
+    mesh->end(it);
 
-    if (write_flag) {
-      apf::writeVtkFiles("core_capVol_before", apfCapMesh);
+    if (!args.before.empty()) apf::writeVtkFiles(args.before.c_str(), mesh);
+    if (args.analytic) { // Cleanup if fields were only for before.vtk
+      apf::destroyField(frameField);
+      apf::destroyField(scaleField);
+      frameField = scaleField = nullptr;
     }
   }
 
-  if (!analytic_flag) {
-    // Pass the field data.
-    in = ma::makeAdvanced(ma::configure(apfCapMesh, scaleField, frameField));
-  } else {
-    // Pass the function.
-    in = ma::makeAdvanced(ma::configure(apfCapMesh, sf));
-  }
-
-  in->shouldSnap = true;
-  in->shouldTransferParametric = true;
-  in->shouldFixShape = true;
-  in->shouldForceAdaptation = true;
-  if (apfCapMesh->getDimension() == 2)
-    in->goodQuality = 0.04; // this is mean-ratio squared
-  else // 3D meshes
-    in->goodQuality = 0.027; // this is the mean-ratio cubed
-  in->maximumIterations = 10;
-
-  if (verbose_flag) {
-    // Adapt with verbose logging but without intermediate VTKs.
-    ma::adaptVerbose(in, false);
-  } else {
-    ma::adapt(in);
-  }
-
-  if (volume_flag) {
-    // We can't verify surface meshes.
-    apfCapMesh->verify();
-  }
-
-  if (write_flag) {
-    apf::writeVtkFiles("core_capVol_after", apfCapMesh);
-    writeCre(cs, "core_capVol_after.cre");
-  }
-
-  /* PRINT ADAPTED MESH INFO */
-  if (verbose_flag) {
-    M_MModel mmdl;
-    m->get_current_model(mmdl);
-    std::string info;
-    m->print_info(mmdl, info);
-    lion_oprint(1, "%s", info.c_str());
+  ma::Input *in = nullptr;
+  if (args.analytic) in = ma::makeAdvanced(ma::configure(mesh, sf.get()));
+  else in = ma::makeAdvanced(ma::configure(mesh, scaleField, frameField));
+  switch (args.balancer) {
+  case Args::Zoltan:
+    in->shouldRunPreZoltan = in->shouldRunMidZoltan = true;
+    in->shouldRunPostZoltan = true;
+    break;
+  case Args::METIS:
+    in->shouldRunPreMetis = in->shouldRunMidMetis = true;
+    in->shouldRunPostMetis = true;
+    break;
+  case Args::Parma:
+    in->shouldRunPreParma = in->shouldRunMidParma = true;
+    in->shouldRunPostParma = true;
+    break;
+  default:;
   }
 
-  // Clean up.
-  if (frameField) apf::destroyField(frameField);
-  if (scaleField) apf::destroyField(scaleField);
-  apf::destroyMesh(apfCapMesh);
-  delete sf;
+  ma::adapt(in);
+  if (!args.after.empty()) apf::writeVtkFiles(args.after.c_str(), mesh);
+  if (PCU->Peers() > 1) migrateHome(mesh);
+  if (PCU->Self() != 0) apf::destroyMesh(mesh);
+  else mesh->switchPCU(oldPCU);
+}
 
-  // Exit calls.
-  gmi_cap_stop();
-  }
-  pcu::Finalize();
-}
\ No newline at end of file
+} // namespace
diff --git a/test/capVolSizeFields.h b/test/capVolSizeFields.h
index eeb5b13f1..cae266fb9 100644
--- a/test/capVolSizeFields.h
+++ b/test/capVolSizeFields.h
@@ -105,6 +105,41 @@ class Shock : public ma::AnisotropicFunction
     ma::Mesh* mesh;
 };
 
+/**
+ * \brief Cube shock version 2.
+ */
+class Shock2 : public ma::AnisotropicFunction
+{
+  public:
+    Shock2(ma::Mesh* m) {
+      mesh = m;
+      ma::Vector bmin, bmax;
+      ma::getBoundingBox(mesh, bmin, bmax);
+      x0 = (bmax.x() + bmin.x()) / 2;
+      double apothem = (bmax.x() - bmin.x()) / 2;
+      hmin = apothem / 128, hmax = apothem / 2;
+    }
+    virtual void getValue(ma::Entity* v, ma::Matrix& R, ma::Vector& H)
+    {
+      ma::Vector p = ma::getPosition(mesh, v);
+      double x = p[0];
+      double s = hmax * (1 - exp(-4*std::abs(x - x0)));
+      double hx = std::max(hmin, s);
+      double hy = hmax;
+      double hz = hmax;
+      /* // principal directions */
+      R = ma::Matrix(
+        1.0, 0.0, 0.0,
+        0.0, 1.0, 0.0,
+		    0.0, 0.0, 1.0
+      );
+      H = ma::Vector(hx, hy, hz);
+    }
+  private:
+    ma::Mesh* mesh;
+    double x0, hmin, hmax;
+};
+
 // Shock along a cylinder like the one in
 // https://www.scorec.rpi.edu/~cwsmith/SCOREC//SimModSuite_2023.1-230428dev/MeshSimAdapt/group__MSAEx__GenAnisoMesh.html
 class CylBoundaryLayer : public ma::AnisotropicFunction {
@@ -116,35 +151,22 @@ class CylBoundaryLayer : public ma::AnisotropicFunction {
   CylBoundaryLayer(ma::Mesh* mesh): m(mesh) {
     ma::Vector bmin(HUGE_VAL, HUGE_VAL, HUGE_VAL),
                bmax(-HUGE_VAL, -HUGE_VAL, -HUGE_VAL);
-    ma::Iterator *it = m->begin(0);
-    ma::Entity *e;
-    while ((e = m->iterate(it))) {
-      ma::Vector p = ma::getPosition(m, e);
-
-      if (p.x() < bmin.x()) bmin.x() = p.x();
-      if (p.x() > bmax.x()) bmax.x() = p.x();
-      if (p.y() < bmin.y()) bmin.y() = p.y();
-      if (p.y() > bmax.y()) bmax.y() = p.y();
-      if (p.z() < bmin.z()) bmin.z() = p.z();
-      if (p.z() > bmax.z()) bmax.z() = p.z();
-    }
-    m->end(it);
-
+    ma::getBoundingBox(mesh, bmin, bmax);
     ma::Vector len = bmax - bmin;
     if (len.x() > len.y() && len.x() > len.z()) {
-      lion_oprint(1, "Shock2 along X axis.\n");
+      lion_oprint(1, "CylBoundaryLayer along X axis.\n");
       length = len.x();
       axis = Axis::X;
       // Radial lengths should be the same, but we find average diameter and
       // then halve that.
       radius = (len.y() / 2 + len.z() / 2) / 2;
     } else if (len.y() > len.x() && len.y() > len.z()) {
-      lion_oprint(1, "Shock2 along Y axis.\n");
+      lion_oprint(1, "CylBoundaryLayer along Y axis.\n");
       length = len.y();
       axis = Axis::Y;
       radius = (len.x() / 2 + len.z() / 2) / 2;
     } else {
-      lion_oprint(1, "Shock2 along Z axis.\n");
+      lion_oprint(1, "CylBoundaryLayer along Z axis.\n");
       length = len.z();
       axis = Axis::Z;
       radius = (len.x() / 2 + len.y() / 2) / 2;
diff --git a/test/cap_closestPoint.cc b/test/cap_closestPoint.cc
index 4db2b59c2..0902af2f7 100644
--- a/test/cap_closestPoint.cc
+++ b/test/cap_closestPoint.cc
@@ -1,9 +1,9 @@
 #include <apf.h>
+#include <apfMesh2.h>
 #include <apfCAP.h>
 #include <gmi_cap.h>
 #include <lionPrint.h>
 #include <pcu_util.h>
-#include <CapstoneModule.h>
 #include <PCU.h>
 
 int main (int argc, char* argv[]) {
@@ -11,16 +11,15 @@ int main (int argc, char* argv[]) {
   { // pcu object scope
   pcu::PCU PCUObj;
   lion_set_verbosity(1);
+  gmi_cap_start();
   gmi_register_cap();
 
   PCU_ALWAYS_ASSERT(argc == 2);
-  std::string creFile(argv[1]);
+  const char* creFile(argv[1]);
   // 1. Load model.
-  CapstoneModule cs("cap_inClosureOf", "Geometry Database : SMLIB",
-    "Mesh Database : Create", "Attribution Database : Create");
-  cs.load_files(v_string(1, creFile));
+  gmi_model* model = gmi_cap_load(creFile);
   // 2. CreateMesh.
-  apf::Mesh2* m = apf::createMesh(cs.get_mesh(), cs.get_geometry(), &PCUObj);
+  apf::Mesh2* m = apf::createCapMesh(model, &PCUObj);
 
   PCU_ALWAYS_ASSERT(m->canGetClosestPoint());
 
@@ -35,6 +34,7 @@ int main (int argc, char* argv[]) {
   PCU_ALWAYS_ASSERT((to - apf::Vector3(0.5, 0.34, 0.5)).getLength() < 0.0001);
 
   apf::destroyMesh(m);
+  gmi_cap_stop();
   } // pcu object scope
   pcu::Finalize();
 }
diff --git a/test/cap_inClosureOf.cc b/test/cap_inClosureOf.cc
index 63ebb828b..87e7a4a6f 100644
--- a/test/cap_inClosureOf.cc
+++ b/test/cap_inClosureOf.cc
@@ -1,26 +1,26 @@
+#include <PCU.h>
 #include <apf.h>
 #include <apfCAP.h>
+#include <apfMesh2.h>
+#include <gmi.h>
 #include <gmi_cap.h>
 #include <lionPrint.h>
 #include <pcu_util.h>
-#include <CapstoneModule.h>
-#include <PCU.h>
 
 int main (int argc, char* argv[]) {
   pcu::Init(&argc, &argv);
   { // pcu object scope
   pcu::PCU PCUObj;
   lion_set_verbosity(1);
+  gmi_cap_start();
   gmi_register_cap();
 
   PCU_ALWAYS_ASSERT(argc == 2);
-  std::string creFile(argv[1]);
+  const char* creFile(argv[1]);
   // 1. Load model.
-  CapstoneModule cs("cap_inClosureOf", "Geometry Database : SMLIB",
-    "Mesh Database : Create", "Attribution Database : Create");
-  cs.load_files(v_string(1, creFile));
+  gmi_model* model = gmi_load(creFile);
   // 2. CreateMesh.
-  apf::Mesh2* m = apf::createMesh(cs.get_mesh(), cs.get_geometry(), &PCUObj);
+  apf::Mesh2* m = apf::createCapMesh(model, &PCUObj);
   // 3. Get region 1 ModelEntity*.
   apf::ModelEntity* rgn = m->findModelEntity(3, 1);
   PCU_ALWAYS_ASSERT(rgn);
@@ -45,6 +45,7 @@ int main (int argc, char* argv[]) {
   PCU_ALWAYS_ASSERT(m->isInClosureOf(m->findModelEntity(0, 8), f1));
 
   apf::destroyMesh(m);
+  gmi_cap_stop();
   } // pcu object scope
   pcu::Finalize();
 }
diff --git a/test/cap_smooth.cc b/test/cap_smooth.cc
index c66cb6115..b84b742b9 100644
--- a/test/cap_smooth.cc
+++ b/test/cap_smooth.cc
@@ -2,7 +2,7 @@
 #include <apfCAP.h>
 
 int main (void) {
-  PCU_ALWAYS_ASSERT(apf::has_smoothCAPAnisoSizes());
+  PCU_ALWAYS_ASSERT(apf::has_smoothCapAnisoSizes());
   // FIXME: Test apf::smoothCAPAnisoSizes.
   return 0;
 }
diff --git a/test/testing.cmake b/test/testing.cmake
index 9f93a1d66..507420cbe 100644
--- a/test/testing.cmake
+++ b/test/testing.cmake
@@ -978,18 +978,45 @@ if (PCU_COMPRESS)
     WORKING_DIRECTORY ${MDIR}/run)
 endif()
 
-if(ENABLE_CAPSTONE)
-  mpi_test(capCyl 1 ./capVol -v 1 ${MESHES}/cap/cyl_surf_only.cre)
-  mpi_test(capWing 1 ./capVol -v 2 ${MESHES}/cap/wing_surf_only.cre)
-  mpi_test(capCube 1 ./capVol -v 3 ${MESHES}/cap/cube_surf_only.cre)
-  mpi_test(capCyl2 1 ./capVol -v 4 ${MESHES}/cap/cyl_surf_only.cre)
-  mpi_test(capVolCyl 1 ./capVol -vg 1 ${MESHES}/cap/cyl_surf_only.cre)
-  mpi_test(capVolWing 1 ./capVol -vg 2 ${MESHES}/cap/wing_surf_only.cre)
-  mpi_test(capVolCube 1 ./capVol -vg 3 ${MESHES}/cap/cube_surf_only.cre)
-  mpi_test(capVolCyl2 1 ./capVol -vg 4 ${MESHES}/cap/cyl_surf_only.cre)
+if(PUMI_ENABLE_CAPSTONE)
+  set(MDIR "${MESHES}/cap")
+  mpi_test(capVol_Wing 1 ./capVol 2 ${MDIR}/wing_surf_only.cre capVol_Wing.cre
+  )
+  mpi_test(capVol_WingMds 1
+    ./capVol -m 2 ${MDIR}/wing_surf_only.cre capVol_WingMds.cre
+  )
+  mpi_test(capVol_BLCylMds3D 1
+    ./capVol -gm 4 ${MDIR}/cyl_surf_only.cre capVol_BLCylMds3D.cre
+  )
+  mpi_test(capVol_BLCylMds3D-4 4
+    ./capVol -gm 4 ${MDIR}/cyl_surf_only.cre capVol_BLCylMds3D.cre
+  )
+  mpi_test(capVol_CubeMds3D 1
+    ./capVol -gm 3 ${MDIR}/cube_surf_only.cre capVol_CubeMds3D.cre
+  )
+  mpi_test(capVol_CubeMds3D-4 4
+    ./capVol -gm 3 ${MDIR}/cube_surf_only.cre capVol_CubeMds3D.cre
+  )
+  if(PUMI_TEST_CAPVOL_EXTRA)
+    mpi_test(capVol_BLCyl3D 1
+      ./capVol -g 4 ${MDIR}/cyl_surf_only.cre capVol_BLCyl3D.cre
+    )
+    mpi_test(capVol_Cube3D 1
+      ./capVol -g 3 ${MDIR}/cube_surf_only.cre capVol_Cube3D.cre
+    )
+    mpi_test(capVol_Cyl3D 1 ./capVol -g 1 ${MDIR}/cyl_surf_only.cre out.cre)
+    mpi_test(capVol_Cube3D 1 ./capVol -g 3 ${MDIR}/cube_surf_only.cre out.cre)
+    mpi_test(capVol_BLCyl3D 1 ./capVol -g 4 ${MDIR}/cyl_surf_only.cre out.cre)
+    mpi_test(capVol_CylMds 1 ./capVol -m 1 ${MDIR}/cyl_surf_only.cre out.cre)
+    mpi_test(capVol_CubeMds 1 ./capVol -m 3 ${MDIR}/cube_surf_only.cre out.cre)
+    mpi_test(capVol_BLCylMds 1 ./capVol -m 4 ${MDIR}/cyl_surf_only.cre out.cre)
+    mpi_test(capVol_CylMds3D 1 ./capVol -gm 1 ${MDIR}/cyl_surf_only.cre out.cre)
+    mpi_test(capVol_Wing3D 1 ./capVol -g 2 ${MDIR}/wing_surf_only.cre out.cre)
+    mpi_test(capVol_WingMds3D 1 ./capVol -gm 2 ${MDIR}/wing_surf_only.cre out.cre)
+  endif()
   mpi_test(cap_inClosureOf 1 ./cap_inClosureOf ${MESHES}/cap/cube_surf_only.cre)
   mpi_test(cap_closestPoint 1 ./cap_closestPoint ${MESHES}/cap/cube_surf_only.cre)
-  if(HAVE_CAPSTONE_SIZINGMETRICTOOL)
+  if(PUMI_HAS_CAPSTONE_SIZINGMETRICTOOL)
     mpi_test(cap_smooth 1 ./cap_smooth)
   endif()
 endif()