diff --git a/src/gsC1SurfBasisEdge.h b/src/gsC1SurfBasisEdge.h
index b2c21dc..e89a956 100644
--- a/src/gsC1SurfBasisEdge.h
+++ b/src/gsC1SurfBasisEdge.h
@@ -263,22 +263,24 @@ namespace gismo
             const gsGeometry<T> & patch = m_mp.patch(0);
 
             // Initialize domain element iterator
-            typename gsBasis<T>::domainIter domIt = m_geo.basis(0).makeDomainIterator(boundary::none);
-
-#ifdef _OPENMP
-            for ( domIt->next(tid); domIt->good(); domIt->next(nt) )
-#else
-            for (; domIt->good(); domIt->next() )
-#endif
+            typename gsBasis<T>::domainIter domIt    = m_geo.basis(0).domain()->beginBdr(boundary::none);
+            typename gsBasis<T>::domainIter domItEnd = m_geo.basis(0).domain()->endBdr(boundary::none);
+
+#           ifdef _OPENMP
+            domIt += tid;
+            for ( ; domIt<domItEnd; domIt+=nt )
+#           else
+            for (; domIt<domItEnd; ++domIt )
+#           endif
             {
                 // Map the Quadrature rule to the element
-                quRule.mapTo( domIt->lowerCorner(), domIt->upperCorner(), quNodes, quWeights );
+                quRule.mapTo( domIt.lowerCorner(), domIt.upperCorner(), quNodes, quWeights );
 
                 // Perform required evaluations on the quadrature nodes
                 visitor_.evaluate(bf_index, typeBf, basis_g1, basis_geo, basis_plus, basis_minus, patch, quNodes, m_uv, m_gD, m_isBoundary);
 
                 // Assemble on element
-                visitor_.assemble(*domIt, quWeights);
+                visitor_.assemble(domIt, quWeights);
 
                 // Push to global matrix and right-hand side vector
 #pragma omp critical(localToGlobal)
diff --git a/src/gsC1SurfBasisVertex.h b/src/gsC1SurfBasisVertex.h
index 6e10b3f..eef665d 100644
--- a/src/gsC1SurfBasisVertex.h
+++ b/src/gsC1SurfBasisVertex.h
@@ -230,22 +230,24 @@ namespace gismo
             const gsGeometry<T> & patch = m_mp.patch(0);
 
             // Initialize domain element iterator
-            typename gsBasis<T>::domainIter domIt = m_geo.basis(0).makeDomainIterator(boundary::none);
-
-#ifdef _OPENMP
-            for ( domIt->next(tid); domIt->good(); domIt->next(nt) )
-#else
-            for (; domIt->good(); domIt->next() )
-#endif
+            typename gsBasis<T>::domainIter domIt    = m_geo.basis(0).domain()->beginBdr(boundary::none);
+            typename gsBasis<T>::domainIter domItEnd = m_geo.basis(0).domain()->endBdr(boundary::none);
+
+#           ifdef _OPENMP
+            domIt += tid;
+            for ( ; domIt<domItEnd; domIt+=nt )
+#           else
+            for (; domIt<domItEnd; ++domIt )
+#           endif
             {
                 // Map the Quadrature rule to the element
-                quRule.mapTo( domIt->lowerCorner(), domIt->upperCorner(), quNodes, quWeights );
+                quRule.mapTo( domIt.lowerCorner(), domIt.upperCorner(), quNodes, quWeights );
 #pragma omp critical(evaluate)
                 // Perform required evaluations on the quadrature nodes
                 visitor_.evaluate(basis_g1, basis_geo, m_basis_plus, m_basis_minus, patch, quNodes, m_gD, m_isBoundary, m_Phi);
 
                 // Assemble on element
-                visitor_.assemble(*domIt, quWeights);
+                visitor_.assemble(domIt, quWeights);
 
                 // Push to global matrix and right-hand side vector
 #pragma omp critical(localToGlobal)
diff --git a/src/gsC1SurfGluingData.h b/src/gsC1SurfGluingData.h
index 6330bf3..764f57e 100644
--- a/src/gsC1SurfGluingData.h
+++ b/src/gsC1SurfGluingData.h
@@ -176,23 +176,25 @@ class gsC1SurfGluingData : public gsC1SurfGD<T>
             // Initialize reference quadrature rule and visitor data
             visitor_.initialize(basis,quRule);
 
-            // Initialize domain element iterator -- using unknown 0
-            typename gsBasis<T>::domainIter domIt = basis.makeDomainIterator(boundary::none);
-
-    #ifdef _OPENMP
-            for ( domIt->next(tid); domIt->good(); domIt->next(nt) )
-    #else
-            for (; domIt->good(); domIt->next() )
-    #endif
+            // Initialize domain element iterator
+            typename gsBasis<T>::domainIter domIt    = basis.domain()->beginBdr(boundary::none);
+            typename gsBasis<T>::domainIter domItEnd = basis.domain()->endBdr(boundary::none);
+
+#           ifdef _OPENMP
+            domIt += tid;
+            for ( ; domIt<domItEnd; domIt+=nt )
+#           else
+            for (; domIt<domItEnd; ++domIt )
+#           endif
             {
                 // Map the Quadrature rule to the element
-                quRule.mapTo( domIt->lowerCorner(), domIt->upperCorner(), quNodes, quWeights );
+                quRule.mapTo( domIt.lowerCorner(), domIt.upperCorner(), quNodes, quWeights );
 
                 // Perform required evaluations on the quadrature nodes
                 visitor_.evaluate(quNodes, this->m_mp);
 
                 // Assemble on element
-                visitor_.assemble(*domIt, quWeights);
+                visitor_.assemble(domIt, quWeights);
 
                 // Push to global matrix and right-hand side vector
     #pragma omp critical(localToGlobal)
@@ -223,23 +225,25 @@ class gsC1SurfGluingData : public gsC1SurfGD<T>
             // Initialize reference quadrature rule and visitor data
             visitor_Beta.initialize(basis,quRule);
 
-            // Initialize domain element iterator -- using unknown 0
-            typename gsBasis<T>::domainIter domIt = basis.makeDomainIterator(boundary::none);
+            // Initialize domain element iterator
+            typename gsBasis<T>::domainIter domIt    = basis.domain()->beginBdr(boundary::none);
+            typename gsBasis<T>::domainIter domItEnd = basis.domain()->endBdr(boundary::none);
 
-    #ifdef _OPENMP
-            for ( domIt->next(tid); domIt->good(); domIt->next(nt) )
-    #else
-            for (; domIt->good(); domIt->next() )
-    #endif
+#           ifdef _OPENMP
+            domIt += tid;
+            for ( ; domIt<domItEnd; domIt+=nt )
+#           else
+            for (; domIt<domItEnd; ++domIt )
+#           endif
             {
                 // Map the Quadrature rule to the element
-                quRule.mapTo( domIt->lowerCorner(), domIt->upperCorner(), quNodes, quWeights );
+                quRule.mapTo( domIt.lowerCorner(), domIt.upperCorner(), quNodes, quWeights );
 
                 // Perform required evaluations on the quadrature nodes
                 visitor_Beta.evaluateBeta(quNodes, this->m_mp, sol);
 
                 // Assemble on element
-                visitor_Beta.assembleBeta(*domIt, quWeights);
+                visitor_Beta.assembleBeta(domIt, quWeights);
 
                 // Push to global matrix and right-hand side vector
     #pragma omp critical(localToGlobal)
diff --git a/src/gsC1SurfGluingDataAssembler.h b/src/gsC1SurfGluingDataAssembler.h
index 715446e..9bccb4c 100644
--- a/src/gsC1SurfGluingDataAssembler.h
+++ b/src/gsC1SurfGluingDataAssembler.h
@@ -127,7 +127,7 @@ class gsC1SurfGluingDataAssembler : public gsAssembler<T>
         m_system_beta_S.matrix().makeCompressed();
 
     }
-    
+
     template <class T, class bhVisitor>
     inline void gsC1SurfGluingDataAssembler<T, bhVisitor>::apply(bhVisitor & visitor,
                                                                index_t patchIndex,
@@ -155,23 +155,25 @@ class gsC1SurfGluingDataAssembler : public gsAssembler<T>
 
             //const gsGeometry<T> & patch = m_geo.patch(patchIndex); // 0 = patchindex
 
-            // Initialize domain element iterator -- using unknown 0
-            typename gsBasis<T>::domainIter domIt = basis.makeDomainIterator(boundary::none);
+            // Initialize domain element iterator
+            typename gsBasis<T>::domainIter domIt    = basis.domain()->beginBdr(boundary::none);
+            typename gsBasis<T>::domainIter domItEnd = basis.domain()->endBdr(boundary::none);
 
-    #ifdef _OPENMP
-            for ( domIt->next(tid); domIt->good(); domIt->next(nt) )
-    #else
-            for (; domIt->good(); domIt->next() )
-    #endif
+#           ifdef _OPENMP
+            domIt += tid;
+            for ( ; domIt<domItEnd; domIt+=nt )
+#           else
+            for (; domIt<domItEnd; ++domIt )
+#           endif
             {
                 // Map the Quadrature rule to the element
-                quRule.mapTo( domIt->lowerCorner(), domIt->upperCorner(), quNodes, quWeights );
+                quRule.mapTo( domIt.lowerCorner(), domIt.upperCorner(), quNodes, quWeights );
 
                 // Perform required evaluations on the quadrature nodes
                 visitor_.evaluate(basis, quNodes, m_uv, m_mp, m_gamma, m_isBoundary);
 
                 // Assemble on element
-                visitor_.assemble(*domIt, quWeights);
+                visitor_.assemble(domIt, quWeights);
 
                 // Push to global matrix and right-hand side vector
     #pragma omp critical(localToGlobal)
diff --git a/src/gsC1SurfGluingDataVisitor.h b/src/gsC1SurfGluingDataVisitor.h
index 1bc12a9..befaeea 100644
--- a/src/gsC1SurfGluingDataVisitor.h
+++ b/src/gsC1SurfGluingDataVisitor.h
@@ -263,7 +263,7 @@ class gsC1SurfGluingDataVisitor
         localRhsBeta.setZero(numActiveBeta, 1 ); //multiple right-hand sides
     }
 
-    inline void assemble(gsDomainIterator<T>    & element,
+    inline void assemble(gsDomainIteratorWrapper<T>    & element,
                          const gsVector<T>      & quWeights)
     {
         GISMO_UNUSED(element);
@@ -285,7 +285,7 @@ class gsC1SurfGluingDataVisitor
         }
     }
 
-    inline void assembleBeta(gsDomainIterator<T>    & element,
+    inline void assembleBeta(gsDomainIteratorWrapper<T>    & element,
                              const gsVector<T>      & quWeights)
     {
         GISMO_UNUSED(element);
diff --git a/src/gsC1SurfVisitorBasisEdge.h b/src/gsC1SurfVisitorBasisEdge.h
index 6e0baa4..f516799 100644
--- a/src/gsC1SurfVisitorBasisEdge.h
+++ b/src/gsC1SurfVisitorBasisEdge.h
@@ -158,7 +158,7 @@ namespace gismo
             } // Patch 1
         } // evaluate1
 
-        inline void assemble(gsDomainIterator<T>    & element,
+        inline void assemble(gsDomainIteratorWrapper<T>    & element,
                              const gsVector<T>      & quWeights)
         {
             GISMO_UNUSED(element);
diff --git a/src/gsC1SurfVisitorBasisVertex.h b/src/gsC1SurfVisitorBasisVertex.h
index 09bc5cc..68a1d0d 100644
--- a/src/gsC1SurfVisitorBasisVertex.h
+++ b/src/gsC1SurfVisitorBasisVertex.h
@@ -419,7 +419,7 @@ namespace gismo
 
         } // evaluate
 
-        inline void assemble(gsDomainIterator<T>    & element,
+        inline void assemble(gsDomainIteratorWrapper<T>    & element,
                              const gsVector<T>      & quWeights)
         {
             GISMO_UNUSED(element);
diff --git a/src/gsContainerBasis.h b/src/gsContainerBasis.h
index 2f79a5b..49ea2b8 100644
--- a/src/gsContainerBasis.h
+++ b/src/gsContainerBasis.h
@@ -130,12 +130,19 @@ namespace gismo
 
         index_t nPieces() const {return basisContainer.size();}
 
+        virtual memory::shared_ptr<gsDomain<T> > domain() const
+        {
+            return basisContainer[0].domain();
+        }
+
+        GISMO_DEPRECATED
         typename gsBasis<T>::domainIter makeDomainIterator(const boxSide & side) const override
         {
             // Using the inner basis for iterating
             return basisContainer[0].makeDomainIterator(side);
         }
 
+        GISMO_DEPRECATED
         typename gsBasis<T>::domainIter makeDomainIterator() const override
         {
             // Using the inner basis for iterating