diff --git a/src/SIM/SIM1D.C b/src/SIM/SIM1D.C
index eaefd167..dd637bea 100644
--- a/src/SIM/SIM1D.C
+++ b/src/SIM/SIM1D.C
@@ -15,6 +15,7 @@
 #include "IFEM.h"
 #include "ASMs1D.h"
 #include "Functions.h"
+#include "ModelGenerator.h"
 #include "Utilities.h"
 #include "Vec3Oper.h"
 #include "tinyxml.h"
@@ -535,54 +536,10 @@ bool SIM1D::createFEMmodel (char)
 }
 
 
-ASMbase* SIM1D::createDefaultGeometry (const TiXmlElement* geo) const
+ModelGenerator* SIM1D::createModelGenerator(const TiXmlElement* geo) const
 {
-  std::string g2("100 1 0 0\n");
-  g2.append(1,'0'+nsd);
-
-  bool rational=false;
-  utl::getAttribute(geo,"rational",rational);
-  if (rational)
-    IFEM::cout << "\t Rational basis\n";
-  g2.append(rational?" 1":" 0");
-  g2.append("\n2 2\n0 0 1 1\n");
-
-  unsigned char d;
-  std::string XYZ;
-  if (utl::getAttribute(geo,"X0",XYZ))
-  {
-    IFEM::cout <<"\tX0 = "<< XYZ << std::endl;
-    g2.append(XYZ);
-  }
-  else
-  {
-    g2.append("0.0");
-    for (d = 1; d < nsd; d++)
-      g2.append(" 0.0");
-  }
-  if (rational)
-    g2.append(" 1.0");
-  g2.append("\n");
-  if (utl::getAttribute(geo,"X1",XYZ))
-  {
-    IFEM::cout <<"\tX1 = "<< XYZ << std::endl;
-    g2.append(XYZ);
-  }
-  else
-  {
-    XYZ = "1.0";
-    if (utl::getAttribute(geo,"L",XYZ))
-      IFEM::cout <<"\tLength scale = "<< XYZ << std::endl;
-    g2.append(XYZ);
-    for (d = 1; d < nsd; d++)
-      g2.append(" 0.0");
-  }
-  if (rational)
-    g2.append(" 1.0");
-  g2.append("\n");
-
-  std::istringstream unitLine(g2);
-  return this->readPatch(unitLine,0,{nf});
+  IFEM::cout <<"  Using default linear geometry basis on unit domain [0,1]";
+  return new DefaultGeometry1D(geo);
 }
 
 
diff --git a/src/SIM/SIM1D.h b/src/SIM/SIM1D.h
index 774e4e99..95ea46da 100644
--- a/src/SIM/SIM1D.h
+++ b/src/SIM/SIM1D.h
@@ -82,6 +82,10 @@ protected:
   //! \brief Parses the twist angle description along the curve.
   bool parseTwist(const TiXmlElement* elem);
 
+  //! \brief Creates a default single-patch model generator.
+  //! \param[in] geo XML element containing geometry defintion
+  virtual ModelGenerator* createModelGenerator(const TiXmlElement* geo) const;
+
   //! \brief Parses a data section from an XML document.
   //! \param[in] elem The XML element to parse
   virtual bool parse(const TiXmlElement* elem);
@@ -105,9 +109,6 @@ protected:
   //! \details Reimplemented to account for twist angle in beam problems.
   virtual bool createFEMmodel(char = 'y');
 
-  //! \brief Creates a default single-patch geometry.
-  virtual ASMbase* createDefaultGeometry(const TiXmlElement* geo) const;
-
 protected:
   unsigned char nf; //!< Number of scalar fields
 
diff --git a/src/SIM/SIM2D.C b/src/SIM/SIM2D.C
index ca05d49f..0b9e1d05 100644
--- a/src/SIM/SIM2D.C
+++ b/src/SIM/SIM2D.C
@@ -16,6 +16,7 @@
 #include "ASMs2DC1.h"
 #include "ImmersedBoundaries.h"
 #include "Functions.h"
+#include "ModelGenerator.h"
 #include "Utilities.h"
 #include "Vec3Oper.h"
 #include "tinyxml.h"
@@ -814,60 +815,10 @@ void SIM2D::clonePatches (const PatchVec& patches,
 }
 
 
-ASMbase* SIM2D::createDefaultGeometry (const TiXmlElement* geo) const
+ModelGenerator* SIM2D::createModelGenerator(const TiXmlElement* geo) const
 {
-  std::string g2("200 1 0 0\n");
-  g2.append(nsd > 2 ? "3" : "2");
-  bool rational=false;
-  utl::getAttribute(geo,"rational",rational);
-  if (rational)
-    IFEM::cout << "\t Rational basis\n";
-  g2.append(rational?" 1":" 0");
-  g2.append("\n2 2\n0 0 1 1\n2 2\n0 0 1 1");
-
-  Vec3 X0;
-  std::string corner;
-  if (utl::getAttribute(geo,"X0",corner)) {
-    std::stringstream str(corner); str >> X0;
-    IFEM::cout <<"  Corner: "<< X0 << std::endl;
-  }
-
-  double scale = 1.0;
-  if (utl::getAttribute(geo,"scale",scale))
-    IFEM::cout <<"  Scale: "<< scale << std::endl;
-
-  double Lx = 1.0, Ly = 1.0;
-  if (utl::getAttribute(geo,"Lx",Lx))
-    IFEM::cout <<"  Length in X: "<< Lx << std::endl;
-  Lx *= scale;
-  if (utl::getAttribute(geo,"Ly",Ly))
-    IFEM::cout <<"  Length in Y: "<< Ly << std::endl;
-  Ly *= scale;
-
-  std::stringstream str;
-  str <<"\n"<< X0.x <<" "<< X0.y;
-  if (nsd > 2) str <<" 0.0";
-  if (rational) str << " 1.0";
-  g2.append(str.str());
-  str.str("");
-  str <<"\n"<< X0.x+Lx <<" "<< X0.y;
-  if (nsd > 2) str <<" 0.0";
-  if (rational) str << " 1.0";
-  g2.append(str.str());
-  str.str("");
-  str <<"\n"<< X0.x <<" "<< X0.y+Ly;
-  if (nsd > 2) str <<" 0.0";
-  if (rational) str << " 1.0";
-  g2.append(str.str());
-  str.str("");
-  str <<"\n"<< X0.x+Lx <<" "<< X0.y+Ly;
-  if (nsd > 2) str <<" 0.0";
-  if (rational) str << " 1.0";
-  g2.append(str.str());
-  g2.append("\n");
-
-  std::istringstream unitSquare(g2);
-  return this->readPatch(unitSquare,0,nf);
+  IFEM::cout <<"  Using default linear geometry basis on unit domain [0,1]^2";
+  return new DefaultGeometry2D(geo);
 }
 
 
diff --git a/src/SIM/SIM2D.h b/src/SIM/SIM2D.h
index c35b4de0..5da34aa9 100644
--- a/src/SIM/SIM2D.h
+++ b/src/SIM/SIM2D.h
@@ -131,8 +131,9 @@ protected:
   virtual bool addConstraint(int patch, int lndx, int ldim,
                              int dirs, int code, int& ngnod, char basis = 1);
 
-  //! \brief Creates a default single-patch geometry.
-  virtual ASMbase* createDefaultGeometry(const TiXmlElement* geo) const;
+  //! \brief Creates a default single-patch model generator.
+  //! \param[in] geo XML element containing geometry defintion
+  virtual ModelGenerator* createModelGenerator(const TiXmlElement* geo) const;
 
 protected:
   CharVec nf;         //!< Number of scalar fields
diff --git a/src/SIM/SIM3D.C b/src/SIM/SIM3D.C
index aa864385..6eef2e5c 100644
--- a/src/SIM/SIM3D.C
+++ b/src/SIM/SIM3D.C
@@ -14,6 +14,7 @@
 #include "SIM3D.h"
 #include "ASMs3D.h"
 #include "Functions.h"
+#include "ModelGenerator.h"
 #include "Utilities.h"
 #include "IFEM.h"
 #include "Vec3Oper.h"
@@ -815,77 +816,10 @@ bool SIM3D::readNodes (std::istream& isn, int pchInd, int basis, bool oneBased)
 }
 
 
-ASMbase* SIM3D::createDefaultGeometry (const TiXmlElement* geo) const
+ModelGenerator* SIM3D::createModelGenerator(const TiXmlElement* geo) const
 {
-  std::string g2("700 1 0 0\n3 ");
-
-  bool rational = false;
-  utl::getAttribute(geo,"rational",rational);
-  if (rational)
-    IFEM::cout <<"  Rational basis"<< std::endl;
-
-  g2.append(rational ? "1\n" : "0\n");
-  g2.append("2 2\n0 0 1 1\n"
-            "2 2\n0 0 1 1\n"
-            "2 2\n0 0 1 1\n");
-
-  std::array<double,24> nodes =
-    {{ 0.0, 0.0, 0.0,
-       1.0, 0.0, 0.0,
-       0.0, 1.0, 0.0,
-       1.0, 1.0, 0.0,
-       0.0, 0.0, 1.0,
-       1.0, 0.0, 1.0,
-       0.0, 1.0, 1.0,
-       1.0, 1.0, 1.0 }};
-
-  double scale = 1.0;
-  if (utl::getAttribute(geo,"scale",scale))
-    IFEM::cout <<"\tscale = "<< scale << std::endl;
-
-  double Lx = 1.0, Ly = 1.0, Lz = 1.0;
-  if (utl::getAttribute(geo,"Lx",Lx))
-    IFEM::cout <<"\tLength in X: "<< Lx << std::endl;
-  Lx *= scale;
-  if (utl::getAttribute(geo,"Ly",Ly))
-    IFEM::cout <<"\tLength in Y: "<< Ly << std::endl;
-  Ly *= scale;
-  if (utl::getAttribute(geo,"Lz",Lz))
-    IFEM::cout <<"\tLength in Z: "<< Lz << std::endl;
-  Lz *= scale;
-
-  if (Lx != 1.0)
-    nodes[3] = nodes[9] = nodes[15] = nodes[21] = Lx;
-  if (Ly != 1.0)
-    nodes[7] = nodes[10] = nodes[19] = nodes[22] = Ly;
-  if (Lz != 1.0)
-    nodes[14] = nodes[17] = nodes[20] = nodes[23] = Lz;
-
-  std::string corner;
-  if (utl::getAttribute(geo,"X0",corner)) {
-    std::stringstream str(corner);
-    Vec3 X0;
-    str >> X0;
-    IFEM::cout <<"\tCorner: "<< X0 << std::endl;
-    for (size_t i = 0; i < nodes.size(); i += 3)
-    {
-      nodes[i]   += X0.x;
-      nodes[i+1] += X0.y;
-      nodes[i+2] += X0.z;
-    }
-  }
-
-  for (size_t i = 0; i < nodes.size(); i += 3)
-  {
-    std::stringstream str;
-    for (size_t j = 0; j < 3; j++)
-      str << nodes[i+j] <<" ";
-    g2.append(str.str());
-    g2.append(rational ? "1.0\n" : "\n");
-  }
-
-  std::istringstream unitCube(g2);
-  return this->readPatch(unitCube,0,nf);
+  IFEM::cout <<"  Using default linear geometry basis on unit domain [0,1]^3";
+  return new DefaultGeometry3D(geo);
 }
 
 
diff --git a/src/SIM/SIM3D.h b/src/SIM/SIM3D.h
index 9cbebb8f..b9b742ad 100644
--- a/src/SIM/SIM3D.h
+++ b/src/SIM/SIM3D.h
@@ -137,8 +137,9 @@ protected:
   bool addConstraint(int patch, int lndx, int line, double xi,
                      int dirs, char basis = 1);
 
-  //! \brief Creates a default single-patch geometry.
-  virtual ASMbase* createDefaultGeometry(const TiXmlElement* geo) const;
+  //! \brief Creates a default single-patch model generator.
+  //! \param[in] geo XML element containing geometry defintion
+  virtual ModelGenerator* createModelGenerator(const TiXmlElement* geo) const;
 
 protected:
   CharVec nf;         //!< Number of scalar fields
diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C
index 5664a7ca..1c4378f2 100644
--- a/src/SIM/SIMbase.C
+++ b/src/SIM/SIMbase.C
@@ -29,6 +29,7 @@
 #include "AnaSol.h"
 #include "Vec3Oper.h"
 #include "Functions.h"
+#include "ModelGenerator.h"
 #include "Profiler.h"
 #include "Utilities.h"
 #include "HDF5Writer.h"
@@ -495,15 +496,24 @@ bool SIMbase::parse (const TiXmlElement* elem)
       result = false;
   }
 
+  ModelGenerator* gen = nullptr;
   // Create the default geometry of no patchfile is specified
   if (myModel.empty() && !strcasecmp(elem->Value(),"geometry"))
     if (this->getNoParamDim() > 0 && !elem->FirstChildElement("patchfile"))
     {
-      IFEM::cout <<"  Using default linear geometry basis on unit domain [0,1]";
-      if (this->getNoParamDim() > 1) IFEM::cout <<"^"<< this->getNoParamDim();
-      IFEM::cout << std::endl;
-      nGlPatches = 1;
-      myModel.resize(1,this->createDefaultGeometry(elem));
+      // parse partitioning first
+      for (const TiXmlElement* part = elem->FirstChildElement("partitioning");
+                         part; part = part->NextSiblingElement("partitioning"))
+        result &= this->parseGeometryTag(part);
+
+      gen = this->createModelGenerator(elem);
+      myModel = gen->createGeometry(*this);
+      if (myPatches.empty())
+        nGlPatches = myModel.size();
+
+      TopologySet set = gen->createTopologySets(*this);
+      for (auto& it : set)
+        myEntitys[it.first] = it.second;
     }
 
   if (!strcasecmp(elem->Value(),"linearsolver")) {
@@ -529,6 +539,11 @@ bool SIMbase::parse (const TiXmlElement* elem)
     else if (!strcasecmp(elem->Value(),"discretization"))
       result &= opt.parseDiscretizationTag(child);
 
+  if (gen)
+    gen->createTopology(*this);
+
+  delete gen;
+
   return result;
 }
 
diff --git a/src/SIM/SIMbase.h b/src/SIM/SIMbase.h
index c61f5723..11c57ddd 100644
--- a/src/SIM/SIMbase.h
+++ b/src/SIM/SIMbase.h
@@ -32,6 +32,7 @@ class LinSolParams;
 class TimeStep;
 class SystemVector;
 class Vec4;
+class ModelGenerator;
 namespace LR { struct RefineData; }
 
 //! Property code to integrand map
@@ -674,9 +675,9 @@ public:
   RealFunc* getSclFunc(int code) const;
 
 protected:
-  //! \brief Creates a default single-patch geometry.
+  //! \brief Instantiate a generator for the finite element model.
   //! \param[in] geo XML element containing geometry defintion
-  virtual ASMbase* createDefaultGeometry(const TiXmlElement* geo) const = 0;
+  virtual ModelGenerator* createModelGenerator(const TiXmlElement* geo) const = 0;
 
   //! \brief Initializes material properties for integration of interior terms.
   virtual bool initMaterial(size_t) { return true; }
diff --git a/src/SIM/SIMdummy.h b/src/SIM/SIMdummy.h
index 4ff80961..80b8cd85 100644
--- a/src/SIM/SIMdummy.h
+++ b/src/SIM/SIMdummy.h
@@ -48,6 +48,9 @@ protected:
   { return nullptr; }
   //! \brief Preprocesses the result sampling points.
   virtual void preprocessResultPoints() {}
+  //! \brief Creates a model generator.
+  virtual ModelGenerator* createModelGenerator(const TiXmlElement*) const
+  { return nullptr; }
 };
 
 #endif
diff --git a/src/SIM/SIMgeneric.C b/src/SIM/SIMgeneric.C
index 98f8aa33..921bb695 100644
--- a/src/SIM/SIMgeneric.C
+++ b/src/SIM/SIMgeneric.C
@@ -13,14 +13,17 @@
 
 #include "SIMgeneric.h"
 #include "ASMbase.h"
+#include "ModelGenerator.h"
 
 
 void SIMgeneric::createDefaultModel ()
 {
   if (!myModel.empty()) return;
 
-  nGlPatches = 1;
-  myModel.resize(1,this->createDefaultGeometry(nullptr));
+  ModelGenerator* gen = this->createModelGenerator(nullptr);
+  myModel = gen->createGeometry(*this);
+  nGlPatches = myModel.size();
+  delete gen;
 }