diff --git a/src/Cgismo.h.in b/src/Cgismo.h.in
index 5e5ce3b..0005188 100644
--- a/src/Cgismo.h.in
+++ b/src/Cgismo.h.in
@@ -37,6 +37,10 @@
 // gsModelling
 #include <gsCInterface/gsCFitting.h>
 
+// gsUtils
+#include <gsCInterface/gsCQuasiInterpolate.h>
+#include <gsCInterface/gsCL2Projection.h>
+
 @gsOptional_includes@
 
 //
diff --git a/src/gsCBasis.cpp b/src/gsCBasis.cpp
index 6616a7e..5af6edf 100644
--- a/src/gsCBasis.cpp
+++ b/src/gsCBasis.cpp
@@ -5,22 +5,115 @@
 #include <gsCInterface/gsMacros.h>
 #include <gsCInterface/gsCBasis.h>
 
+using namespace gismo;
+
+template <int dim>
+void gsHTensorBasis_elements_into_impl(gsCBasis * b, bool getKnotBoxes,
+                                                     bool getIndexBoxes,
+                                                     bool getLevels,
+                                                     gsCMatrix*    knotBoxes,
+                                                     gsCMatrixInt* indexBoxes,
+                                                     gsCVectorInt* levels)
+{
+    int N   = RICAST_B(b)->numElements();
+
+    auto * el = RICAST_M(knotBoxes);
+    auto * bx = RICAST_Mi(indexBoxes);
+    auto * lv = RICAST_Vi(levels);
+
+    if (getKnotBoxes) el->resize(dim,2*N);
+    if (getIndexBoxes) bx->resize(2*dim+1,N);
+    if (getLevels) lv->resize(N);
+
+    auto domain = RICAST_B(b)->domain();
+    auto domIt  = domain->beginAll();
+    auto domEnd = domain->endAll();
+
+    GISMO_ENSURE((dynamic_cast<gsHDomainIterator<double,dim> *>(domIt.get())),"The domain iterator is not a hierarchical domain iterator");
+    gsHTensorBasis<dim,double> * basis = dynamic_cast<gsHTensorBasis<dim,double> *>(RICAST_B(b));
+
+    int id=0;
+    gsVector<double,dim> low, upp;
+    for (; domIt<domEnd; ++domIt, ++id)
+    {
+        gsHDomainIterator<double,dim> * domItH = dynamic_cast<gsHDomainIterator<double,dim> *>(domIt.get());
+        low = domIt.lowerCorner();
+        upp = domIt.upperCorner();
+        if (getKnotBoxes)
+        {
+            el->col(2*id) = low;
+            el->col(2*id+1) = upp;
+        }
+        if (getLevels)
+        {
+            lv->at(id) = domItH->getLevel();
+        }
+        if (getIndexBoxes)
+        {
+            for(int j = 0; j < dim;j++)
+            {
+                // Convert the parameter coordinates to (unique) knot indices
+                const gsKnotVector<real_t> & kv = basis->tensorLevel(domItH->getLevel()).knots(j);
+                int k1 = (std::upper_bound(kv.domainUBegin(), kv.domainUEnd(),
+                                        low[j] ) - 1).uIndex();
+                int k2 = (std::upper_bound(kv.domainUBegin(), kv.domainUEnd()+1,
+                                        upp[j] ) - 1).uIndex();
+
+                // Trivial cells trigger some refinement
+                if ( k1 == k2)
+                {
+                    if (0!=k1) {--k1;}
+                    ++k2;
+                }
+
+                // Store the data...
+                (*bx)(0,id) = domItH->getLevel();
+                (*bx)(1+j,id) = k1;
+                (*bx)(1+j+dim,id) = k2;
+            }
+        }
+    }
+}
+
+
 #ifdef __cplusplus
 extern "C"
 {
 #endif
 
+GISMO_EXPORT gsCBasis * gsBasis_read(char* filename)
+{
+    gsFileData<> data(filename);
+    if (data.hasAny< gsBasis<> >())
+    {
+        gsBasis<>::uPtr ptr = data.getAnyFirst< gsBasis<> >();
+        return RICAST_CB(ptr.release());
+    }
+    else
+    {
+        gsWarn<<"[G+Smo] No gsBasis found in file "<<filename<<"\n";
+        return NULL;
+    }
+}
+
+GISMO_EXPORT void gsBasis_write(gsCBasis * obj, char* filename)
+{
+    gsFileData<> data;
+    data.add(*RICAST_B(obj));
+    data.save(filename);
+}
+
 GISMO_EXPORT gsCBasis * gsBSplineBasis_create(gsCKnotVector * KV)
 {
     auto * KV_ptr = RICAST_KV (KV);
-    return RICAST_CB (new gismo::gsBSplineBasis<double>(*KV_ptr) );
+    return RICAST_CB (new gsBSplineBasis<double>(*KV_ptr) );
 }
 
 GISMO_EXPORT gsCBasis* gsTensorBSplineBasis2_create(gsCKnotVector* KV1, gsCKnotVector* KV2)
 {
     auto * KV1_ptr = RICAST_KV (KV1);
     auto * KV2_ptr = RICAST_KV (KV2);
-    return RICAST_CB(new gismo::gsTensorBSplineBasis<2,double>(*KV1_ptr,*KV2_ptr) );
+    return RICAST_CB(new gsTensorBSplineBasis<2,double>(*KV1_ptr,*KV2_ptr) );
 }
 
 GISMO_EXPORT gsCBasis* gsTensorBSplineBasis3_create(gsCKnotVector* KV1, gsCKnotVector* KV2,
@@ -29,7 +122,7 @@ GISMO_EXPORT gsCBasis* gsTensorBSplineBasis3_create(gsCKnotVector* KV1, gsCKnotV
     auto * KV1_ptr = RICAST_KV (KV1);
     auto * KV2_ptr = RICAST_KV (KV2);
     auto * KV3_ptr = RICAST_KV (KV3);
-    return RICAST_CB(new gismo::gsTensorBSplineBasis<3,double>(*KV1_ptr,*KV2_ptr,*KV3_ptr));
+    return RICAST_CB(new gsTensorBSplineBasis<3,double>(*KV1_ptr,*KV2_ptr,*KV3_ptr));
 }
 
 GISMO_EXPORT gsCBasis* gsTensorBSplineBasis4_create(gsCKnotVector* KV1, gsCKnotVector* KV2,
@@ -39,90 +132,95 @@ GISMO_EXPORT gsCBasis* gsTensorBSplineBasis4_create(gsCKnotVector* KV1, gsCKnotV
     auto * KV2_ptr = RICAST_KV (KV2);
     auto * KV3_ptr = RICAST_KV (KV3);
     auto * KV4_ptr = RICAST_KV (KV4);
-    return RICAST_CB(new gismo::gsTensorBSplineBasis<4,double>(*KV1_ptr,*KV2_ptr,*KV3_ptr,*KV4_ptr));
+    return RICAST_CB(new gsTensorBSplineBasis<4,double>(*KV1_ptr,*KV2_ptr,*KV3_ptr,*KV4_ptr));
 }
 
 GISMO_EXPORT gsCBasis * gsNurbsBasis_create(gsCBasis * b, gsCMatrix * weights)
 {
-    auto * basis_ptr = reinterpret_cast< gismo::gsBSplineBasis<double>* >(b);
+    auto * basis_ptr = reinterpret_cast< gsBSplineBasis<double>* >(b);
     auto * w = RICAST_M(weights);
-    return RICAST_CB(new gismo::gsNurbsBasis<double>(basis_ptr,*w));
+    return RICAST_CB(new gsNurbsBasis<double>(basis_ptr,*w));
 }
 
 GISMO_EXPORT gsCBasis* gsTensorNurbsBasis2_create(gsCBasis* b, gsCMatrix * weights)
-{ 
-    auto * basis_ptr = reinterpret_cast< gismo::gsTensorBSplineBasis<2,double>* >(b);
+{
+    auto * basis_ptr = reinterpret_cast< gsTensorBSplineBasis<2,double>* >(b);
     auto * w = RICAST_M(weights);
-    return RICAST_CB(new  gismo::gsTensorNurbsBasis<2,double>(basis_ptr,*w));
+    return RICAST_CB(new  gsTensorNurbsBasis<2,double>(basis_ptr,*w));
 }
 
 GISMO_EXPORT gsCBasis* gsTensorNurbsBasis3_create(gsCBasis* b, gsCMatrix * weights)
 {
-    auto * basis_ptr = reinterpret_cast< gismo::gsTensorBSplineBasis<3,double>* >(b);
+    auto * basis_ptr = reinterpret_cast< gsTensorBSplineBasis<3,double>* >(b);
     auto * w = RICAST_M(weights);
-    return RICAST_CB(new  gismo::gsTensorNurbsBasis<3,double>(basis_ptr,*w));
+    return RICAST_CB(new  gsTensorNurbsBasis<3,double>(basis_ptr,*w));
 }
 
 GISMO_EXPORT gsCBasis* gsTensorNurbsBasis4_create(gsCBasis* b, gsCMatrix * weights)
 {
-    auto * basis_ptr = reinterpret_cast< gismo::gsTensorBSplineBasis<4,double>* >(b);
+    auto * basis_ptr = reinterpret_cast< gsTensorBSplineBasis<4,double>* >(b);
     auto * w = RICAST_M(weights);
-    return RICAST_CB(new  gismo::gsTensorNurbsBasis<4,double>(basis_ptr,*w));
+    return RICAST_CB(new  gsTensorNurbsBasis<4,double>(basis_ptr,*w));
 }
 
 
 GISMO_EXPORT gsCBasis* gsTHBSplineBasis1_create(gsCBasis* b)
 {
-    auto * basis_ptr = reinterpret_cast< gismo::gsTensorBSplineBasis<1,double>* >(b);
-    return RICAST_CB(new gismo::gsTHBSplineBasis<1,double>(*basis_ptr,false));
+    auto * basis_ptr = reinterpret_cast< gsTensorBSplineBasis<1,double>* >(b);
+    return RICAST_CB(new gsTHBSplineBasis<1,double>(*basis_ptr,false));
 }
 
 GISMO_EXPORT gsCBasis* gsTHBSplineBasis2_create(gsCBasis* b)
-{ 
-    auto * basis_ptr = reinterpret_cast< gismo::gsTensorBSplineBasis<2,double>* >(b);
-    return RICAST_CB(new  gismo::gsTHBSplineBasis<2,double>(*basis_ptr,false));
+{
+    auto * basis_ptr = reinterpret_cast< gsTensorBSplineBasis<2,double>* >(b);
+    return RICAST_CB(new  gsTHBSplineBasis<2,double>(*basis_ptr,false));
 }
 
 GISMO_EXPORT gsCBasis* gsTHBSplineBasis3_create(gsCBasis* b)
 {
-    auto * basis_ptr = reinterpret_cast< gismo::gsTensorBSplineBasis<3,double>* >(b);
-    return RICAST_CB(new  gismo::gsTHBSplineBasis<3,double>(*basis_ptr,false));
+    auto * basis_ptr = reinterpret_cast< gsTensorBSplineBasis<3,double>* >(b);
+    return RICAST_CB(new  gsTHBSplineBasis<3,double>(*basis_ptr,false));
 }
 
 GISMO_EXPORT gsCBasis* gsTHBSplineBasis4_create(gsCBasis* b)
 {
-    auto * basis_ptr = reinterpret_cast< gismo::gsTensorBSplineBasis<4,double>* >(b);
-    return RICAST_CB(new  gismo::gsTHBSplineBasis<4,double>(*basis_ptr,false));
+    auto * basis_ptr = reinterpret_cast< gsTensorBSplineBasis<4,double>* >(b);
+    return RICAST_CB(new  gsTHBSplineBasis<4,double>(*basis_ptr,false));
 }
 
 GISMO_EXPORT gsCBasis* gsHBSplineBasis1_create(gsCBasis* b)
 {
-    auto * basis_ptr = reinterpret_cast< gismo::gsTensorBSplineBasis<1,double>* >(b);
-    return RICAST_CB(new gismo::gsHBSplineBasis<1,double>(*basis_ptr,false));
+    auto * basis_ptr = reinterpret_cast< gsTensorBSplineBasis<1,double>* >(b);
+    return RICAST_CB(new gsHBSplineBasis<1,double>(*basis_ptr,false));
 }
 
 GISMO_EXPORT gsCBasis* gsHBSplineBasis2_create(gsCBasis* b)
 {
-    auto * basis_ptr = reinterpret_cast< gismo::gsTensorBSplineBasis<2,double>* >(b);
-    return RICAST_CB(new  gismo::gsHBSplineBasis<2,double>(*basis_ptr,false));
+    auto * basis_ptr = reinterpret_cast< gsTensorBSplineBasis<2,double>* >(b);
+    return RICAST_CB(new  gsHBSplineBasis<2,double>(*basis_ptr,false));
 }
 
 GISMO_EXPORT gsCBasis* gsHBSplineBasis3_create(gsCBasis* b)
 {
-    auto * basis_ptr = reinterpret_cast< gismo::gsTensorBSplineBasis<3,double>* >(b);
-    return RICAST_CB(new  gismo::gsHBSplineBasis<3,double>(*basis_ptr,false));
+    auto * basis_ptr = reinterpret_cast< gsTensorBSplineBasis<3,double>* >(b);
+    return RICAST_CB(new  gsHBSplineBasis<3,double>(*basis_ptr,false));
 }
 
 GISMO_EXPORT gsCBasis* gsHBSplineBasis4_create(gsCBasis* b)
 {
-    auto * basis_ptr = reinterpret_cast< gismo::gsTensorBSplineBasis<4,double>* >(b);
-    return RICAST_CB(new  gismo::gsHBSplineBasis<4,double>(*basis_ptr,false));
+    auto * basis_ptr = reinterpret_cast< gsTensorBSplineBasis<4,double>* >(b);
+    return RICAST_CB(new  gsHBSplineBasis<4,double>(*basis_ptr,false));
 }
 
 //
 // Methods, gsBasis
 //
 
+GISMO_EXPORT gsCBasis* gsBasis_clone(gsCBasis * b)
+{
+    return RICAST_CB(RICAST_B(b)->clone().release());
+}
+
 GISMO_EXPORT void gsBasis_active_into(gsCBasis * b,
                               gsCMatrix * u,
                               gsCMatrixInt * result)
@@ -148,7 +246,7 @@ GISMO_EXPORT void gsBasis_deriv2Single_into(gsCBasis * b,
 
 GISMO_EXPORT gsCBasis * gsBasis_component(gsCBasis * b, int dir)
 {
-    gismo::gsBasis<double> * c = & RICAST_B(b)->component(dir);
+    gsBasis<double> * c = & RICAST_B(b)->component(dir);
     return reinterpret_cast<gsCBasis*>(c);
 }
 
@@ -165,7 +263,7 @@ GISMO_EXPORT int gsBasis_size(gsCBasis * b)
 { return RICAST_B(b)->size(); }
 
 GISMO_EXPORT gsCMatrix* gsBasis_support(gsCBasis * b, int i)
-{ return reinterpret_cast<gsCMatrix*>( new gismo::gsMatrix<double>(RICAST_B(b)->support(i)) ); }
+{ return reinterpret_cast<gsCMatrix*>( new gsMatrix<double>(RICAST_B(b)->support(i)) ); }
 
 GISMO_EXPORT void gsBasis_uniformRefine(gsCBasis * b, int numKnots, int mul, int dir)
 { RICAST_B(b)->uniformRefine(numKnots, mul, dir); }
@@ -179,28 +277,93 @@ GISMO_EXPORT void gsBasis_refineElements(gsCBasis * b, int * boxData, int boxSiz
 GISMO_EXPORT void gsBasis_refine(gsCBasis * b, gsCMatrix * boxes, int refExt)
 { RICAST_B(b)->refine(*RICAST_M(boxes),refExt); }
 
+GISMO_EXPORT void gsBasis_degreeElevate(gsCBasis * b, int i, int dir)
+{ RICAST_B(b)->degreeElevate(i,dir); }
+
+GISMO_EXPORT void gsBasis_boundary_into(gsCBasis * b, int side, gsCMatrixInt * result)
+{ *RICAST_Mi(result) = RICAST_B(b)->boundary(side); }
+
+GISMO_EXPORT void gsBasis_boundaryOffset_into(gsCBasis * b, int side, int offset, gsCMatrixInt * result)
+{ *RICAST_Mi(result) = RICAST_B(b)->boundaryOffset(side,offset); }
+
+GISMO_EXPORT void gsBasis_elements_into(gsCBasis * b, gsCMatrix* elements)
+{
+    auto * el = RICAST_M(elements);
+    el->resize(RICAST_B(b)->domainDim(),2*RICAST_B(b)->numElements());
+    auto domain = RICAST_B(b)->domain();
+    auto domIt  = domain->beginAll();
+    auto domEnd = domain->endAll();
+    int id=0;
+    for (; domIt<domEnd; ++domIt, ++id)
+    {
+        el->col(2*id) = domIt.lowerCorner();
+        el->col(2*id+1) = domIt.upperCorner();
+    }
+}
+
+GISMO_EXPORT void gsBasis_elementsBdr_into(gsCBasis * b, int side, gsCMatrix* elements)
+{
+    auto * el = RICAST_M(elements);
+    el->resize(RICAST_B(b)->domainDim(),2*RICAST_B(b)->numElements());
+    auto domain = RICAST_B(b)->domain();
+    auto domIt  = domain->beginBdr(side);
+    auto domEnd = domain->endBdr(side);
+    int id=0;
+    for (; domIt<domEnd; ++domIt, ++id)
+    {
+        el->col(2*id) = domIt.lowerCorner();
+        el->col(2*id+1) = domIt.upperCorner();
+    }
+}
+
+GISMO_EXPORT void gsHTensorBasis_elements_into(gsCBasis * b, bool getKnotBoxes,
+                                                             bool getIndexBoxes,
+                                                             bool getLevels,
+                                                             gsCMatrix*    knotBoxes,
+                                                             gsCMatrixInt* indexBoxes,
+                                                             gsCVectorInt* levels)
+{
+    switch (RICAST_B(b)->domainDim())
+    {
+        case 1:
+            gsHTensorBasis_elements_into_impl<1>(b,getKnotBoxes,getIndexBoxes,getLevels,knotBoxes,indexBoxes,levels);
+            break;
+        case 2:
+            gsHTensorBasis_elements_into_impl<2>(b,getKnotBoxes,getIndexBoxes,getLevels,knotBoxes,indexBoxes,levels);
+            break;
+        case 3:
+            gsHTensorBasis_elements_into_impl<3>(b,getKnotBoxes,getIndexBoxes,getLevels,knotBoxes,indexBoxes,levels);
+            break;
+        case 4:
+            gsHTensorBasis_elements_into_impl<4>(b,getKnotBoxes,getIndexBoxes,getLevels,knotBoxes,indexBoxes,levels);
+            break;
+        default:
+            GISMO_ERROR("gsHTensorBasis_elements_into: Dimension not supported");
+    }
+}
+
 //
 // Methods, Other
 //
 
 GISMO_EXPORT gsCKnotVector * gsBSplineBasis_knots(gsCBasis * b)
 {
-    gismo::gsKnotVector<double> * KV= &reinterpret_cast< gismo::gsBSplineBasis<double>* >(b)->knots();
+    gsKnotVector<double> * KV= &reinterpret_cast< gsBSplineBasis<double>* >(b)->knots();
     return reinterpret_cast<gsCKnotVector*>(KV);
 }
 
 GISMO_EXPORT gsCKnotVector * gsTensorBSplineBasis_knots(gsCBasis * b, int dir)
 {
-    gismo::gsKnotVector<double> * KV=NULL;
+    gsKnotVector<double> * KV=NULL;
     GISMO_ASSERT(RICAST_B(b)->domainDim()>=dir,"gsTensorBSplineBasis_knots: dir out of range");
     switch (RICAST_B(b)->domainDim())
     {
         case 2:
-            KV = &reinterpret_cast< gismo::gsTensorBSplineBasis<2,double>* >(b)->knots(dir);
+            KV = &reinterpret_cast< gsTensorBSplineBasis<2,double>* >(b)->knots(dir);
         case 3:
-            KV = &reinterpret_cast< gismo::gsTensorBSplineBasis<3,double>* >(b)->knots(dir);
+            KV = &reinterpret_cast< gsTensorBSplineBasis<3,double>* >(b)->knots(dir);
         case 4:
-            KV = &reinterpret_cast< gismo::gsTensorBSplineBasis<4,double>* >(b)->knots(dir);
+            KV = &reinterpret_cast< gsTensorBSplineBasis<4,double>* >(b)->knots(dir);
     }
 
     return reinterpret_cast<gsCKnotVector*>(KV);
@@ -211,13 +374,13 @@ GISMO_EXPORT int gsHTensorBasis_numLevels(gsCBasis * b)
     switch (RICAST_B(b)->domainDim())
     {
         case 1:
-            return reinterpret_cast< gismo::gsHTensorBasis<1,double>* >(b)->numLevels();
+            return reinterpret_cast< gsHTensorBasis<1,double>* >(b)->numLevels();
         case 2:
-            return reinterpret_cast< gismo::gsHTensorBasis<2,double>* >(b)->numLevels();
+            return reinterpret_cast< gsHTensorBasis<2,double>* >(b)->numLevels();
         case 3:
-            return reinterpret_cast< gismo::gsHTensorBasis<3,double>* >(b)->numLevels();
+            return reinterpret_cast< gsHTensorBasis<3,double>* >(b)->numLevels();
         case 4:
-            return reinterpret_cast< gismo::gsHTensorBasis<4,double>* >(b)->numLevels();
+            return reinterpret_cast< gsHTensorBasis<4,double>* >(b)->numLevels();
         default:
             GISMO_ERROR("gsHTensorBasis_numLevels: domainDim not supported");
     }
@@ -228,13 +391,13 @@ GISMO_EXPORT int gsHTensorBasis_maxLevel(gsCBasis * b)
     switch (RICAST_B(b)->domainDim())
     {
         case 1:
-            return reinterpret_cast< gismo::gsHTensorBasis<1,double>* >(b)->maxLevel();
+            return reinterpret_cast< gsHTensorBasis<1,double>* >(b)->maxLevel();
         case 2:
-            return reinterpret_cast< gismo::gsHTensorBasis<2,double>* >(b)->maxLevel();
+            return reinterpret_cast< gsHTensorBasis<2,double>* >(b)->maxLevel();
         case 3:
-            return reinterpret_cast< gismo::gsHTensorBasis<3,double>* >(b)->maxLevel();
+            return reinterpret_cast< gsHTensorBasis<3,double>* >(b)->maxLevel();
         case 4:
-            return reinterpret_cast< gismo::gsHTensorBasis<4,double>* >(b)->maxLevel();
+            return reinterpret_cast< gsHTensorBasis<4,double>* >(b)->maxLevel();
         default:
             GISMO_ERROR("gsHTensorBasis_maxLevel: domainDim not supported");
     }
@@ -245,13 +408,13 @@ GISMO_EXPORT int gsHTensorBasis_levelOf(gsCBasis * b, int i)
     switch (RICAST_B(b)->domainDim())
     {
         case 1:
-            return reinterpret_cast< gismo::gsHTensorBasis<1,double>* >(b)->levelOf(i);
+            return reinterpret_cast< gsHTensorBasis<1,double>* >(b)->levelOf(i);
         case 2:
-            return reinterpret_cast< gismo::gsHTensorBasis<2,double>* >(b)->levelOf(i);
+            return reinterpret_cast< gsHTensorBasis<2,double>* >(b)->levelOf(i);
         case 3:
-            return reinterpret_cast< gismo::gsHTensorBasis<3,double>* >(b)->levelOf(i);
+            return reinterpret_cast< gsHTensorBasis<3,double>* >(b)->levelOf(i);
         case 4:
-            return reinterpret_cast< gismo::gsHTensorBasis<4,double>* >(b)->levelOf(i);
+            return reinterpret_cast< gsHTensorBasis<4,double>* >(b)->levelOf(i);
         default:
             GISMO_ERROR("gsHTensorBasis_levelOf: domainDim not supported");
     }
@@ -263,13 +426,13 @@ GISMO_EXPORT int gsHTensorBasis_getLevelAtPoint(gsCBasis * b, gsCMatrix * Pt)
     switch (RICAST_B(b)->domainDim())
     {
         case 1:
-            return reinterpret_cast< gismo::gsHTensorBasis<1,double>* >(b)->getLevelAtPoint(*m);
+            return reinterpret_cast< gsHTensorBasis<1,double>* >(b)->getLevelAtPoint(*m);
         case 2:
-            return reinterpret_cast< gismo::gsHTensorBasis<2,double>* >(b)->getLevelAtPoint(*m);
+            return reinterpret_cast< gsHTensorBasis<2,double>* >(b)->getLevelAtPoint(*m);
         case 3:
-            return reinterpret_cast< gismo::gsHTensorBasis<3,double>* >(b)->getLevelAtPoint(*m);
+            return reinterpret_cast< gsHTensorBasis<3,double>* >(b)->getLevelAtPoint(*m);
         case 4:
-            return reinterpret_cast< gismo::gsHTensorBasis<4,double>* >(b)->getLevelAtPoint(*m);
+            return reinterpret_cast< gsHTensorBasis<4,double>* >(b)->getLevelAtPoint(*m);
         default:
             GISMO_ERROR("gsHTensorBasis_getLevelAtPoint: domainDim not supported");
     }
@@ -280,13 +443,13 @@ GISMO_EXPORT gsCBasis * gsHTensorBasis_tensorLevel(gsCBasis * b, int l)
     switch (RICAST_B(b)->domainDim())
     {
         case 1:
-            return RICAST_CB(new  gismo::gsBSplineBasis<double>(reinterpret_cast< gismo::gsHTensorBasis<1,double>* >(b)->tensorLevel(l)));
+            return RICAST_CB(new  gsBSplineBasis<double>(reinterpret_cast< gsHTensorBasis<1,double>* >(b)->tensorLevel(l)));
         case 2:
-            return RICAST_CB(new  gismo::gsTensorBSplineBasis<2,double>(reinterpret_cast< gismo::gsHTensorBasis<2,double>* >(b)->tensorLevel(l)));
+            return RICAST_CB(new  gsTensorBSplineBasis<2,double>(reinterpret_cast< gsHTensorBasis<2,double>* >(b)->tensorLevel(l)));
         case 3:
-            return RICAST_CB(new  gismo::gsTensorBSplineBasis<3,double>(reinterpret_cast< gismo::gsHTensorBasis<3,double>* >(b)->tensorLevel(l)));
+            return RICAST_CB(new  gsTensorBSplineBasis<3,double>(reinterpret_cast< gsHTensorBasis<3,double>* >(b)->tensorLevel(l)));
         case 4:
-            return RICAST_CB(new  gismo::gsTensorBSplineBasis<4,double>(reinterpret_cast< gismo::gsHTensorBasis<4,double>* >(b)->tensorLevel(l)));
+            return RICAST_CB(new  gsTensorBSplineBasis<4,double>(reinterpret_cast< gsHTensorBasis<4,double>* >(b)->tensorLevel(l)));
         default:
             GISMO_ERROR("gsHTensorBasis_tensorLevel: domainDim not supported");
     }
@@ -297,13 +460,13 @@ GISMO_EXPORT void gsHTensorBasis_treeLeafSize(gsCBasis * b)
     switch (RICAST_B(b)->domainDim())
     {
         case 1:
-            reinterpret_cast< gismo::gsHTensorBasis<1,double>* >(b)->tree().leafSize();
+            reinterpret_cast< gsHTensorBasis<1,double>* >(b)->tree().leafSize();
         case 2:
-            reinterpret_cast< gismo::gsHTensorBasis<2,double>* >(b)->tree().leafSize();
+            reinterpret_cast< gsHTensorBasis<2,double>* >(b)->tree().leafSize();
         case 3:
-            reinterpret_cast< gismo::gsHTensorBasis<3,double>* >(b)->tree().leafSize();
+            reinterpret_cast< gsHTensorBasis<3,double>* >(b)->tree().leafSize();
         case 4:
-            reinterpret_cast< gismo::gsHTensorBasis<4,double>* >(b)->tree().leafSize();
+            reinterpret_cast< gsHTensorBasis<4,double>* >(b)->tree().leafSize();
         default:
             GISMO_ERROR("gsHTensorBasis_treeLeaveSize: domainDim not supported");
     }
@@ -314,18 +477,20 @@ GISMO_EXPORT void gsHTensorBasis_treePrintLeaves(gsCBasis * b)
     switch (RICAST_B(b)->domainDim())
     {
         case 1:
-            reinterpret_cast< gismo::gsHTensorBasis<1,double>* >(b)->tree().printLeaves();
+            reinterpret_cast< gsHTensorBasis<1,double>* >(b)->tree().printLeaves();
         case 2:
-            reinterpret_cast< gismo::gsHTensorBasis<2,double>* >(b)->tree().printLeaves();
+            reinterpret_cast< gsHTensorBasis<2,double>* >(b)->tree().printLeaves();
         case 3:
-            reinterpret_cast< gismo::gsHTensorBasis<3,double>* >(b)->tree().printLeaves();
+            reinterpret_cast< gsHTensorBasis<3,double>* >(b)->tree().printLeaves();
         case 4:
-            reinterpret_cast< gismo::gsHTensorBasis<4,double>* >(b)->tree().printLeaves();
+            reinterpret_cast< gsHTensorBasis<4,double>* >(b)->tree().printLeaves();
         default:
             GISMO_ERROR("gsHTensorBasis_treePrintLeaves: domainDim not supported");
     }
 }
 
+
+
 #ifdef __cplusplus
 }
 #endif
\ No newline at end of file
diff --git a/src/gsCBasis.h b/src/gsCBasis.h
index ab92602..131a5e4 100644
--- a/src/gsCBasis.h
+++ b/src/gsCBasis.h
@@ -4,6 +4,11 @@ extern "C"
 {
 #endif
 
+    #include<stdbool.h>
+
+    GISMO_EXPORT gsCBasis * gsBasis_read(char* filename);
+    GISMO_EXPORT void       gsBasis_write(gsCBasis * obj, char* filename);
+
 #   define gsBasis_print gsFunctionSet_print
 #   define gsBasis_delete gsFunctionSet_delete
 
@@ -26,6 +31,8 @@ extern "C"
     //
     // Methods, gsBasis
     //
+    GISMO_EXPORT gsCBasis* gsBasis_clone(gsCBasis * b);
+
     GISMO_EXPORT void gsBasis_active_into(gsCBasis * b, gsCMatrix * u, gsCMatrixInt * result);
 
     GISMO_EXPORT void gsBasis_evalSingle_into(gsCBasis * b, int i, gsCMatrix * u, gsCMatrix * result);
@@ -38,11 +45,26 @@ extern "C"
     GISMO_EXPORT int gsBasis_dim(gsCBasis * b);
     GISMO_EXPORT int gsBasis_size(gsCBasis * b);
     GISMO_EXPORT gsCMatrix* gsBasis_support(gsCBasis * b, int i);
-    
+
     GISMO_EXPORT void gsBasis_uniformRefine(gsCBasis * b, int numKnots, int mul, int dir);
     GISMO_EXPORT void gsBasis_refineElements(gsCBasis * b, int * boxData, int boxSize);
     GISMO_EXPORT void gsBasis_refine(gsCBasis * b, gsCMatrix * boxes, int refExt);
 
+    GISMO_EXPORT void gsBasis_degreeElevate(gsCBasis * b, int i, int dir);
+
+    GISMO_EXPORT void gsBasis_boundary_into(gsCBasis * b, int side, gsCMatrixInt * result);
+    GISMO_EXPORT void gsBasis_boundaryOffset_into(gsCBasis * b, int side, int offset, gsCMatrixInt * result);
+
+    GISMO_EXPORT void gsBasis_elements_into(gsCBasis * b, gsCMatrix*);
+    GISMO_EXPORT void gsBasis_elementsBdr_into(gsCBasis * b, int, gsCMatrix*);
+
+    GISMO_EXPORT void gsHTensorBasis_elements_into(gsCBasis * b, bool getKnotBoxes,
+                                                                 bool getIndexBoxes,
+                                                                 bool getLevels,
+                                                                 gsCMatrix*    knotBoxes,
+                                                                 gsCMatrixInt* indexBoxes,
+                                                                 gsCVectorInt* levels);
+
     // TODO:
     // - DegreeElevate
 
diff --git a/src/gsCBoundaryConditions.cpp b/src/gsCBoundaryConditions.cpp
index bb55877..7f41932 100644
--- a/src/gsCBoundaryConditions.cpp
+++ b/src/gsCBoundaryConditions.cpp
@@ -8,9 +8,11 @@ extern "C"
 {
 #endif
 
+using namespace gismo;
+
 GISMO_EXPORT gsCBoundaryConditions * gsBoundaryConditions_create()
 {
-    return RICAST_CBC(new gismo::gsBoundaryConditions<double>());
+    return RICAST_CBC(new gsBoundaryConditions<double>());
 }
 
 GISMO_EXPORT void gsBoundaryConditions_addCondition(gsCBoundaryConditions * bc,
@@ -19,18 +21,18 @@ GISMO_EXPORT void gsBoundaryConditions_addCondition(gsCBoundaryConditions * bc,
                                                     int ctype,
                                                     gsCFunctionSet * fun,
                                                     int unknown,
-                                                    int component,
-                                                    bool parametric)
+                                                    bool parametric,
+                                                    int component)
 {
-    gismo::boxSide bside(side);
-    gismo::gsFunctionSet<double> * f_ptr = RICAST_F(fun);
+    boxSide bside(side);
+    gsFunctionSet<double> * f_ptr = RICAST_F(fun);
     RICAST_BC(bc)->addCondition(patch,
                                 bside,
-                                (gismo::condition_type::type)ctype,
+                                (condition_type::type)ctype,
                                 f_ptr,
                                 unknown,
-                                component,
-                                parametric);
+                                parametric,
+                                component);
 }
 
 GISMO_EXPORT void gsBoundaryConditions_addCornerValue(gsCBoundaryConditions * bc,
@@ -40,7 +42,7 @@ GISMO_EXPORT void gsBoundaryConditions_addCornerValue(gsCBoundaryConditions * bc
                                                       int unknown,
                                                       int component)
 {
-    gismo::boxCorner bcorner(corner);
+    boxCorner bcorner(corner);
     RICAST_BC(bc)->addCornerValue(bcorner, value, patch, unknown, component);
 }
 
diff --git a/src/gsCFitting.cpp b/src/gsCFitting.cpp
index 078d0c6..4bb7a10 100644
--- a/src/gsCFitting.cpp
+++ b/src/gsCFitting.cpp
@@ -9,58 +9,60 @@ extern "C"
 {
 #endif
 
+using namespace gismo;
+
 GISMO_EXPORT gsCFitting * gsFitting_create(gsCMatrix * param_values, gsCMatrix * points, gsCBasis * basis)
 {
     auto * param_values_ptr = RICAST_M(param_values);
     auto * points_ptr = RICAST_M(points);
     auto * basis_ptr = RICAST_B(basis);
-    return reinterpret_cast<gsCFitting*>(new gismo::gsFitting<double>(*param_values_ptr, *points_ptr, *basis_ptr));
+    return reinterpret_cast<gsCFitting*>(new gsFitting<double>(*param_values_ptr, *points_ptr, *basis_ptr));
 }
 
 GISMO_EXPORT void gsFitting_delete(gsCFitting * fitter)
 {
-    delete reinterpret_cast<gismo::gsFitting<double>*>(fitter);
+    delete reinterpret_cast<gsFitting<double>*>(fitter);
 }
 
 GISMO_EXPORT void gsFitting_compute(gsCFitting * fitter, double lambda)
 {
-    reinterpret_cast<gismo::gsFitting<double>*>(fitter)->compute(lambda);
+    reinterpret_cast<gsFitting<double>*>(fitter)->compute(lambda);
 }
 
 GISMO_EXPORT void gsFitting_parameterCorrection(gsCFitting * fitter, double accuracy, int maxIter, double tolOrth)
 {
-    reinterpret_cast<gismo::gsFitting<double>*>(fitter)->parameterCorrection(accuracy, maxIter, tolOrth);
+    reinterpret_cast<gsFitting<double>*>(fitter)->parameterCorrection(accuracy, maxIter, tolOrth);
 }
 
 GISMO_EXPORT void gsFitting_computeErrors(gsCFitting * fitter)
 {
-    reinterpret_cast<gismo::gsFitting<double>*>(fitter)->computeErrors();
+    reinterpret_cast<gsFitting<double>*>(fitter)->computeErrors();
 }
 
 GISMO_EXPORT double gsFitting_minPointError(gsCFitting * fitter)
 {
-    return reinterpret_cast<gismo::gsFitting<double>*>(fitter)->minPointError();
+    return reinterpret_cast<gsFitting<double>*>(fitter)->minPointError();
 }
 
 GISMO_EXPORT double gsFitting_maxPointError(gsCFitting * fitter)
 {
-    return reinterpret_cast<gismo::gsFitting<double>*>(fitter)->maxPointError();
+    return reinterpret_cast<gsFitting<double>*>(fitter)->maxPointError();
 }
 
 GISMO_EXPORT double* gsFitting_pointWiseErrors(gsCFitting * fitter)
 {
-    const double * errors = reinterpret_cast< const double* >(reinterpret_cast<gismo::gsFitting<double>*>(fitter)->pointWiseErrors().data());
+    const double * errors = reinterpret_cast< const double* >(reinterpret_cast<gsFitting<double>*>(fitter)->pointWiseErrors().data());
     return const_cast<double *>(errors);
 }
 
 GISMO_EXPORT int gsFitting_numPointsBelow(gsCFitting * fitter, double threshold)
 {
-    return reinterpret_cast<gismo::gsFitting<double>*>(fitter)->numPointsBelow(threshold);
+    return reinterpret_cast<gsFitting<double>*>(fitter)->numPointsBelow(threshold);
 }
 
 GISMO_EXPORT gsCGeometry* gsFitting_result(gsCFitting * fitter)
 {
-    auto * result = reinterpret_cast<gismo::gsFitting<double>*>(fitter)->result();
+    auto * result = reinterpret_cast<gsFitting<double>*>(fitter)->result();
     return RICAST_CG(result->clone().release());
 }
 
diff --git a/src/gsCFunctionExpr.cpp b/src/gsCFunctionExpr.cpp
index 963e1c5..7143c2d 100644
--- a/src/gsCFunctionExpr.cpp
+++ b/src/gsCFunctionExpr.cpp
@@ -3,70 +3,72 @@
 #include <gsCInterface/gsCMemory.h>
 #include <gsCInterface/gsMacros.h>
 
+using namespace gismo;
+
 #ifdef __cplusplus
 extern "C"
 {
 #endif
 
-GISMO_EXPORT gsCFunctionExpr * gsFunctionExpr1_create(const char * expression_string, 
+GISMO_EXPORT gsCFunctionExpr * gsFunctionExpr1_create(const char * expression_string,
                                                       short_t ddim)
 {
-    return RICAST_CF(new gismo::gsFunctionExpr<double>(expression_string, 
+    return RICAST_CF(new gsFunctionExpr<double>(expression_string,
                                                        ddim));
 }
 
-GISMO_EXPORT gsCFunctionExpr * gsFunctionExpr2_create(const char * expression_string1, 
-                                                      const char * expression_string2, 
+GISMO_EXPORT gsCFunctionExpr * gsFunctionExpr2_create(const char * expression_string1,
+                                                      const char * expression_string2,
                                                       short_t ddim)
 {
-    return RICAST_CF(new gismo::gsFunctionExpr<double>(expression_string1, 
+    return RICAST_CF(new gsFunctionExpr<double>(expression_string1,
                                                        expression_string2, ddim));
 }
 
-GISMO_EXPORT gsCFunctionExpr * gsFunctionExpr3_create(const char * expression_string1, 
-                                                      const char * expression_string2, 
-                                                      const char * expression_string3, 
+GISMO_EXPORT gsCFunctionExpr * gsFunctionExpr3_create(const char * expression_string1,
+                                                      const char * expression_string2,
+                                                      const char * expression_string3,
                                                       short_t ddim)
 {
-    return RICAST_CF(new gismo::gsFunctionExpr<double>(expression_string1, 
-                                                       expression_string2, 
-                                                       expression_string3, 
+    return RICAST_CF(new gsFunctionExpr<double>(expression_string1,
+                                                       expression_string2,
+                                                       expression_string3,
                                                        ddim));
 }
 
-GISMO_EXPORT gsCFunctionExpr * gsFunctionExpr4_create(const char * expression_string1, 
-                                                      const char * expression_string2, 
-                                                      const char * expression_string3, 
-                                                      const char * expression_string4, 
+GISMO_EXPORT gsCFunctionExpr * gsFunctionExpr4_create(const char * expression_string1,
+                                                      const char * expression_string2,
+                                                      const char * expression_string3,
+                                                      const char * expression_string4,
                                                       short_t ddim)
 {
-    return RICAST_CF(new gismo::gsFunctionExpr<double>(expression_string1, 
-                                                       expression_string2, 
-                                                       expression_string3, 
-                                                       expression_string4, 
+    return RICAST_CF(new gsFunctionExpr<double>(expression_string1,
+                                                       expression_string2,
+                                                       expression_string3,
+                                                       expression_string4,
                                                        ddim));
 }
 
-GISMO_EXPORT gsCFunctionExpr * gsFunctionExpr9_create(const char * expression_string1, 
-                                                      const char * expression_string2, 
-                                                      const char * expression_string3, 
-                                                      const char * expression_string4, 
-                                                      const char * expression_string5, 
-                                                      const char * expression_string6, 
-                                                      const char * expression_string7, 
-                                                      const char * expression_string8, 
-                                                      const char * expression_string9, 
+GISMO_EXPORT gsCFunctionExpr * gsFunctionExpr9_create(const char * expression_string1,
+                                                      const char * expression_string2,
+                                                      const char * expression_string3,
+                                                      const char * expression_string4,
+                                                      const char * expression_string5,
+                                                      const char * expression_string6,
+                                                      const char * expression_string7,
+                                                      const char * expression_string8,
+                                                      const char * expression_string9,
                                                       short_t ddim)
 {
-    return RICAST_CF(new gismo::gsFunctionExpr<double>(expression_string1, 
-                                                       expression_string2, 
-                                                       expression_string3, 
-                                                       expression_string4, 
-                                                       expression_string5, 
-                                                       expression_string6, 
-                                                       expression_string7, 
-                                                       expression_string8, 
-                                                       expression_string9, 
+    return RICAST_CF(new gsFunctionExpr<double>(expression_string1,
+                                                       expression_string2,
+                                                       expression_string3,
+                                                       expression_string4,
+                                                       expression_string5,
+                                                       expression_string6,
+                                                       expression_string7,
+                                                       expression_string8,
+                                                       expression_string9,
                                                        ddim));
 }
 
diff --git a/src/gsCFunctionSet.cpp b/src/gsCFunctionSet.cpp
index 4b69651..59c8ad5 100644
--- a/src/gsCFunctionSet.cpp
+++ b/src/gsCFunctionSet.cpp
@@ -5,13 +5,15 @@
 #include <gsCInterface/gsCMemory.h>
 #include <gsCInterface/gsMacros.h>
 
+using namespace gismo;
+
 #ifdef __cplusplus
 extern "C"
 {
 #endif
 
 GISMO_EXPORT void gsFunctionSet_print(gsCFunctionSet * fs)
-{ gsInfo<<RICAST_F(fs); }
+{ gsInfo<<*RICAST_F(fs)<<"\n"; }
 
 GISMO_EXPORT void gsFunctionSet_delete(gsCFunctionSet * ptr)
 { delete RICAST_F(ptr); }
@@ -26,7 +28,7 @@ GISMO_EXPORT int gsFunctionSet_nPieces(gsCFunctionSet * fs)
 { return RICAST_F(fs)->nPieces(); }
 
 GISMO_EXPORT gsCMatrix* gsFunctionSet_support(gsCFunctionSet * fs)
-{ return RICAST_CM( new gismo::gsMatrix<double>(RICAST_F(fs)->support()) ); }
+{ return RICAST_CM( new gsMatrix<double>(RICAST_F(fs)->support()) ); }
 
 GISMO_EXPORT void gsFunctionSet_eval_into(gsCFunctionSet * fs,
                             gsCMatrix * u,
@@ -37,6 +39,36 @@ GISMO_EXPORT void gsFunctionSet_deriv_into(gsCFunctionSet * fs,
                              gsCMatrix * u,
                              gsCMatrix * result)
 { RICAST_F(fs)->deriv_into(*RICAST_M(u), *RICAST_M(result) ); }
+
+GISMO_EXPORT void gsFunctionSet_deriv2_into(gsCFunctionSet * fs,
+                                            gsCMatrix * u,
+                                            gsCMatrix * result)
+{ RICAST_F(fs)->deriv2_into(*RICAST_M(u), *RICAST_M(result) ); }
+
+GISMO_EXPORT void gsFunctionSet_evalAllDers1_into(gsCFunctionSet * fs,
+                                                  gsCMatrix * u,
+                                                  gsCMatrix * values,
+                                                  gsCMatrix * deriv)
+{
+    std::vector<gsMatrix<> > result(2);
+    RICAST_F(fs)->evalAllDers_into(*RICAST_M(u), 1, result, false);
+    *RICAST_M(values) = give(result[0]);
+    *RICAST_M(deriv) = give(result[1]);
+}
+
+GISMO_EXPORT void gsFunctionSet_evalAllDers2_into(gsCFunctionSet * fs,
+                                                  gsCMatrix * u,
+                                                  gsCMatrix * values,
+                                                  gsCMatrix * deriv,
+                                                  gsCMatrix * deriv2)
+{
+    std::vector<gsMatrix<> > result(3);
+    RICAST_F(fs)->evalAllDers_into(*RICAST_M(u), 2, result, false);
+    *RICAST_M(values) = give(result[0]);
+    *RICAST_M(deriv) = give(result[1]);
+    *RICAST_M(deriv2) = give(result[2]);
+}
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/src/gsCFunctionSet.h b/src/gsCFunctionSet.h
index 3dfefd8..31a8ab9 100644
--- a/src/gsCFunctionSet.h
+++ b/src/gsCFunctionSet.h
@@ -21,6 +21,23 @@ extern "C"
                                                gsCMatrix * u,
                                                gsCMatrix * result);
 
+    GISMO_EXPORT void gsFunctionSet_deriv2_into(gsCFunctionSet * fs,
+                                                gsCMatrix * u,
+                                                gsCMatrix * result);
+
+    // gsFunctionSet::evalAllDers_into(u, 1, result,false);
+    GISMO_EXPORT void gsFunctionSet_evalAllDers1_into(gsCFunctionSet * fs,
+                                                      gsCMatrix * u,
+                                                      gsCMatrix * values,
+                                                      gsCMatrix * deriv);
+
+    // gsFunctionSet::evalAllDers_into(u, 2, result,false);
+    GISMO_EXPORT void gsFunctionSet_evalAllDers2_into(gsCFunctionSet * fs,
+                                                      gsCMatrix * u,
+                                                      gsCMatrix * values,
+                                                      gsCMatrix * deriv,
+                                                      gsCMatrix * deriv2);
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/src/gsCGeometry.cpp b/src/gsCGeometry.cpp
index b21ec18..bb3426f 100644
--- a/src/gsCGeometry.cpp
+++ b/src/gsCGeometry.cpp
@@ -12,137 +12,172 @@ extern "C"
 {
 #endif
 
+using namespace gismo;
+
+GISMO_EXPORT gsCGeometry * gsGeometry_read(char* filename)
+{
+    gsFileData<> data(filename);
+    if (data.hasAny< gsGeometry<> >())
+    {
+        gsGeometry<>::uPtr ptr = data.getAnyFirst< gsGeometry<> >();
+        return RICAST_CG(ptr.release());
+    }
+    else
+    {
+        gsWarn<<"[G+Smo] No gsGeometry found in file "<<filename<<"\n";
+        return NULL;
+    }
+}
+
+GISMO_EXPORT void gsGeometry_write(gsCGeometry * obj, char* filename)
+{
+    gsFileData<> data;
+    data.add(*RICAST_G(obj));
+    data.save(filename);
+}
+
 GISMO_EXPORT gsCGeometry* gsBSpline_create(gsCBasis* b, gsCMatrix * coefs)
 {
-    auto * basis_ptr = reinterpret_cast< gismo::gsBSplineBasis<double>* >(b);
+    auto * basis_ptr = reinterpret_cast< gsBSplineBasis<double>* >(b);
     auto * m = RICAST_M(coefs);
-    return RICAST_CG(new gismo::gsBSpline<double>(*basis_ptr,*m));
+    return RICAST_CG(new gsBSpline<double>(*basis_ptr,*m));
 }
 
 GISMO_EXPORT gsCGeometry* gsTensorBSpline2_create(gsCBasis* b, gsCMatrix * coefs)
 {
-    gismo::gsTensorBSplineBasis<2,double>* basis_ptr = reinterpret_cast< gismo::gsTensorBSplineBasis<2,double>* >(b);
+    gsTensorBSplineBasis<2,double>* basis_ptr = reinterpret_cast< gsTensorBSplineBasis<2,double>* >(b);
     auto * m = RICAST_M(coefs);
-    return RICAST_CG(new gismo::gsTensorBSpline<2,double>(*basis_ptr,*m));
+    return RICAST_CG(new gsTensorBSpline<2,double>(*basis_ptr,*m));
 }
 
 GISMO_EXPORT gsCGeometry* gsTensorBSpline3_create(gsCBasis* b, gsCMatrix * coefs)
 {
-    gismo::gsTensorBSplineBasis<3,double>* basis_ptr = reinterpret_cast< gismo::gsTensorBSplineBasis<3,double>* >(b);
+    gsTensorBSplineBasis<3,double>* basis_ptr = reinterpret_cast< gsTensorBSplineBasis<3,double>* >(b);
     auto * m = RICAST_M(coefs);
-    return RICAST_CG(new gismo::gsTensorBSpline<3,double>(*basis_ptr,*m));
+    return RICAST_CG(new gsTensorBSpline<3,double>(*basis_ptr,*m));
 }
 
 GISMO_EXPORT gsCGeometry* gsTensorBSpline4_create(gsCBasis* b, gsCMatrix * coefs)
 {
-    auto * basis_ptr = reinterpret_cast< gismo::gsTensorBSplineBasis<4,double>* >(b);
+    auto * basis_ptr = reinterpret_cast< gsTensorBSplineBasis<4,double>* >(b);
     auto * m = RICAST_M(coefs);
-    return RICAST_CG(new gismo::gsTensorBSpline<4,double>(*basis_ptr,*m));
+    return RICAST_CG(new gsTensorBSpline<4,double>(*basis_ptr,*m));
 }
 
 GISMO_EXPORT gsCGeometry* gsNurbs_create(gsCBasis* b, gsCMatrix * coefs)
 {
-    auto * basis_ptr = reinterpret_cast< gismo::gsNurbsBasis<double>* >(b);
+    auto * basis_ptr = reinterpret_cast< gsNurbsBasis<double>* >(b);
     auto * m = RICAST_M(coefs);
-    return RICAST_CG(new gismo::gsNurbs<double>(*basis_ptr,*m));
+    return RICAST_CG(new gsNurbs<double>(*basis_ptr,*m));
 }
 
 GISMO_EXPORT gsCGeometry* gsTensorNurbs2_create(gsCBasis* b, gsCMatrix * coefs)
 {
-    gismo::gsTensorNurbsBasis<2,double>* basis_ptr = reinterpret_cast< gismo::gsTensorNurbsBasis<2,double>* >(b);
+    gsTensorNurbsBasis<2,double>* basis_ptr = reinterpret_cast< gsTensorNurbsBasis<2,double>* >(b);
     auto * m = RICAST_M(coefs);
-    return RICAST_CG(new gismo::gsTensorNurbs<2,double>(*basis_ptr,*m));
+    return RICAST_CG(new gsTensorNurbs<2,double>(*basis_ptr,*m));
 }
 
 GISMO_EXPORT gsCGeometry* gsTensorNurbs3_create(gsCBasis* b, gsCMatrix * coefs)
 {
-    gismo::gsTensorNurbsBasis<3,double>* basis_ptr = reinterpret_cast< gismo::gsTensorNurbsBasis<3,double>* >(b);
+    gsTensorNurbsBasis<3,double>* basis_ptr = reinterpret_cast< gsTensorNurbsBasis<3,double>* >(b);
     auto * m = RICAST_M(coefs);
-    return RICAST_CG(new gismo::gsTensorNurbs<3,double>(*basis_ptr,*m));
+    return RICAST_CG(new gsTensorNurbs<3,double>(*basis_ptr,*m));
 }
 
 GISMO_EXPORT gsCGeometry* gsTensorNurbs4_create(gsCBasis* b, gsCMatrix * coefs)
 {
-    auto * basis_ptr = reinterpret_cast< gismo::gsTensorNurbsBasis<4,double>* >(b);
+    auto * basis_ptr = reinterpret_cast< gsTensorNurbsBasis<4,double>* >(b);
     auto * m = RICAST_M(coefs);
-    return RICAST_CG(new gismo::gsTensorNurbs<4,double>(*basis_ptr,*m));
+    return RICAST_CG(new gsTensorNurbs<4,double>(*basis_ptr,*m));
 }
 
 
 GISMO_EXPORT gsCGeometry* gsTHBSpline1_create(gsCBasis* b, gsCMatrix * coefs)
 {
-    auto * basis_ptr = reinterpret_cast< gismo::gsTHBSplineBasis<1,double>* >(b);
+    auto * basis_ptr = reinterpret_cast< gsTHBSplineBasis<1,double>* >(b);
     auto * m = RICAST_M(coefs);
-    return RICAST_CG(new gismo::gsTHBSpline<1,double>(*basis_ptr,*m));
+    return RICAST_CG(new gsTHBSpline<1,double>(*basis_ptr,*m));
 }
 
 GISMO_EXPORT gsCGeometry* gsTHBSpline2_create(gsCBasis* b, gsCMatrix * coefs)
 {
-    auto * basis_ptr = reinterpret_cast< gismo::gsTHBSplineBasis<2,double>* >(b);
+    auto * basis_ptr = reinterpret_cast< gsTHBSplineBasis<2,double>* >(b);
     auto * m = RICAST_M(coefs);
-    return RICAST_CG(new gismo::gsTHBSpline<2,double>(*basis_ptr,*m));
+    return RICAST_CG(new gsTHBSpline<2,double>(*basis_ptr,*m));
 }
 
 GISMO_EXPORT gsCGeometry* gsTHBSpline3_create(gsCBasis* b, gsCMatrix * coefs)
 {
-    auto * basis_ptr = reinterpret_cast< gismo::gsTHBSplineBasis<3,double>* >(b);
+    auto * basis_ptr = reinterpret_cast< gsTHBSplineBasis<3,double>* >(b);
     auto * m = RICAST_M(coefs);
-    return RICAST_CG(new gismo::gsTHBSpline<3,double>(*basis_ptr,*m));
+    return RICAST_CG(new gsTHBSpline<3,double>(*basis_ptr,*m));
 }
 
 GISMO_EXPORT gsCGeometry* gsTHBSpline4_create(gsCBasis* b, gsCMatrix * coefs)
 {
-    auto * basis_ptr = reinterpret_cast< gismo::gsTHBSplineBasis<4,double>* >(b);
+    auto * basis_ptr = reinterpret_cast< gsTHBSplineBasis<4,double>* >(b);
     auto * m = RICAST_M(coefs);
-    return RICAST_CG(new gismo::gsTHBSpline<4,double>(*basis_ptr,*m));
+    return RICAST_CG(new gsTHBSpline<4,double>(*basis_ptr,*m));
 }
 
 GISMO_EXPORT gsCGeometry* gsHBSpline1_create(gsCBasis* b, gsCMatrix * coefs)
 {
-    auto * basis_ptr = reinterpret_cast< gismo::gsHBSplineBasis<1,double>* >(b);
+    auto * basis_ptr = reinterpret_cast< gsHBSplineBasis<1,double>* >(b);
     auto * m = RICAST_M(coefs);
-    return RICAST_CG(new gismo::gsHBSpline<1,double>(*basis_ptr,*m));
+    return RICAST_CG(new gsHBSpline<1,double>(*basis_ptr,*m));
 }
 
 GISMO_EXPORT gsCGeometry* gsHBSpline2_create(gsCBasis* b, gsCMatrix * coefs)
 {
-    auto * basis_ptr = reinterpret_cast< gismo::gsHBSplineBasis<2,double>* >(b);
+    auto * basis_ptr = reinterpret_cast< gsHBSplineBasis<2,double>* >(b);
     auto * m = RICAST_M(coefs);
-    return RICAST_CG(new gismo::gsHBSpline<2,double>(*basis_ptr,*m));
+    return RICAST_CG(new gsHBSpline<2,double>(*basis_ptr,*m));
 }
 
 GISMO_EXPORT gsCGeometry* gsHBSpline3_create(gsCBasis* b, gsCMatrix * coefs)
 {
-    auto * basis_ptr = reinterpret_cast< gismo::gsHBSplineBasis<3,double>* >(b);
+    auto * basis_ptr = reinterpret_cast< gsHBSplineBasis<3,double>* >(b);
     auto * m = RICAST_M(coefs);
-    return RICAST_CG(new gismo::gsHBSpline<3,double>(*basis_ptr,*m));
+    return RICAST_CG(new gsHBSpline<3,double>(*basis_ptr,*m));
 }
 
 GISMO_EXPORT gsCGeometry* gsHBSpline4_create(gsCBasis* b, gsCMatrix * coefs)
 {
-    auto * basis_ptr = reinterpret_cast< gismo::gsHBSplineBasis<4,double>* >(b);
+    auto * basis_ptr = reinterpret_cast< gsHBSplineBasis<4,double>* >(b);
     auto * m = RICAST_M(coefs);
-    return RICAST_CG(new gismo::gsHBSpline<4,double>(*basis_ptr,*m));
+    return RICAST_CG(new gsHBSpline<4,double>(*basis_ptr,*m));
 }
 
 GISMO_EXPORT gsCGeometry* gsTensorBSpline2_slice(gsCGeometry * g, int direction, double parameter)
 {
-    auto * g_ptr = reinterpret_cast< gismo::gsTensorBSpline<2,double>* >(g);
-    typedef typename gismo::gsTensorBSpline<2,double>::BoundaryGeometryType GeometryBdr;
+    auto * g_ptr = reinterpret_cast< gsTensorBSpline<2,double>* >(g);
+    typedef typename gsTensorBSpline<2,double>::BoundaryGeometryType GeometryBdr;
     GeometryBdr * bdr = new GeometryBdr();
     g_ptr->slice(direction, parameter, *bdr);
     return RICAST_CG(bdr);
 }
 
+GISMO_EXPORT gsCGeometry* gsGeometry_clone(gsCGeometry * g)
+{
+    return RICAST_CG(RICAST_G(g)->clone().release());
+}
 
 GISMO_EXPORT gsCBasis* gsGeometry_basis(gsCGeometry * g)
-{ return RICAST_CB(&RICAST_G(g)->basis()); }
+{
+    return RICAST_CB(RICAST_G(g)->basis().clone().release());
+}
 
 GISMO_EXPORT void gsGeometry_coefs_into(gsCGeometry * g, gsCMatrix * coefs)
 {
     *RICAST_M(coefs) = RICAST_G(g)->coefs();
 }
 
+GISMO_EXPORT void gsGeometry_setCoefs(gsCGeometry * g, gsCMatrix * coefs)
+{
+    RICAST_G(g)->coefs() = *RICAST_M(coefs);
+}
+
 GISMO_EXPORT void gsGeometry_uniformRefine(gsCGeometry * g, int numKnots, int mul, int dir)
 { RICAST_G(g)->uniformRefine(numKnots, mul, dir); }
 
@@ -152,6 +187,12 @@ GISMO_EXPORT void gsGeometry_refineElements(gsCGeometry * g, int * boxData, int
     RICAST_G(g)->refineElements(boxes);
 }
 
+GISMO_EXPORT void gsGeometry_refine(gsCGeometry * g, gsCMatrix * boxes, int refExt)
+{
+    std::vector<int> boxData = RICAST_G(g)->basis().asElements(*RICAST_M(boxes),refExt);
+    RICAST_G(g)->refineElements(boxData);
+}
+
 GISMO_EXPORT void gsGeometry_degreeElevate(gsCGeometry * g, int i, int dir)
 { RICAST_G(g)->degreeElevate(i, dir); }
 
diff --git a/src/gsCGeometry.h b/src/gsCGeometry.h
index d7ca940..f97dc4c 100644
--- a/src/gsCGeometry.h
+++ b/src/gsCGeometry.h
@@ -4,6 +4,9 @@ extern "C"
 {
 #endif
 
+    GISMO_EXPORT gsCGeometry * gsGeometry_read(char* filename);
+    GISMO_EXPORT void          gsGeometry_write(gsCGeometry * obj, char* filename);
+
 #   define gsGeometry_print gsFunctionSet_print
 #   define gsGeometry_delete gsFunctionSet_delete
 #   define gsBSpline_delete gsFunctionSet_delete
@@ -30,13 +33,17 @@ extern "C"
 
     GISMO_EXPORT gsCGeometry* gsTensorBSpline2_slice(gsCGeometry * g, int direction, double parameter);
 
+    GISMO_EXPORT gsCGeometry* gsGeometry_clone(gsCGeometry * g);
+
     GISMO_EXPORT gsCBasis* gsGeometry_basis(gsCGeometry * g);
 
     GISMO_EXPORT void gsGeometry_coefs_into(gsCGeometry * g, gsCMatrix * coef);
+    GISMO_EXPORT void gsGeometry_setCoefs(gsCGeometry * g, gsCMatrix * coef);
 
 
     GISMO_EXPORT void gsGeometry_uniformRefine(gsCGeometry * b, int numKnots, int mul, int dir);
     GISMO_EXPORT void gsGeometry_refineElements(gsCGeometry * b, int * boxData, int boxSize);
+    GISMO_EXPORT void gsGeometry_refine(gsCGeometry * b, gsCMatrix * boxes, int refExt);
 
     GISMO_EXPORT void gsGeometry_degreeElevate(gsCGeometry * b, int i, int dir);
 
diff --git a/src/gsCGeometryTransform.cpp b/src/gsCGeometryTransform.cpp
index 14ef390..2161ec7 100644
--- a/src/gsCGeometryTransform.cpp
+++ b/src/gsCGeometryTransform.cpp
@@ -8,6 +8,8 @@
 #include <gsCInterface/gsCBasis.h>
 #include <gsCInterface/gsMacros.h>
 
+using namespace gismo;
+
 #ifdef __cplusplus
 extern "C"
 {
@@ -18,11 +20,11 @@ extern "C"
     {
         auto * g_ptr = RICAST_G(g);
         auto * mm = RICAST_M(m);
-        auto * vv = RICAST_V(v);  
-        return RICAST_CG(new gismo::gsGeometryTransform<double>(g_ptr,*mm, *vv));
+        auto * vv = RICAST_V(v);
+        return RICAST_CG(new gsGeometryTransform<double>(g_ptr,*mm, *vv));
     }
 
-    
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/src/gsCKnotVector.cpp b/src/gsCKnotVector.cpp
index bc87442..08eb1b8 100644
--- a/src/gsCKnotVector.cpp
+++ b/src/gsCKnotVector.cpp
@@ -3,48 +3,58 @@
 #include <gsCInterface/gsCTypes.h>
 #include <gsCInterface/gsCKnotVector.h>
 
+using namespace gismo;
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif
+
 GISMO_EXPORT gsCKnotVector * gsKnotVector_create(double* knots, int size)
 {
-    return reinterpret_cast<gsCKnotVector*>(new gismo::gsKnotVector<double>(-1,knots,knots+size));
+    return reinterpret_cast<gsCKnotVector*>(new gsKnotVector<double>(-1,knots,knots+size));
 }
 
 GISMO_EXPORT void gsKnotVector_delete(gsCKnotVector * kv)
 {
-    delete reinterpret_cast<gismo::gsKnotVector<double>*>(kv);
+    delete reinterpret_cast<gsKnotVector<double>*>(kv);
 }
 
 GISMO_EXPORT void gsKnotVector_print(gsCKnotVector * kv)
 {
-    reinterpret_cast<gismo::gsKnotVector<double>*>(kv)->print(gsInfo);
+    reinterpret_cast<gsKnotVector<double>*>(kv)->print(gsInfo);
 }
 
 GISMO_EXPORT double* gsKnotVector_data(gsCKnotVector * kv)
 {
-    return const_cast<double *>(reinterpret_cast<gismo::gsKnotVector<double>*>(kv)->get().data());
+    return const_cast<double *>(reinterpret_cast<gsKnotVector<double>*>(kv)->get().data());
 }
 
 GISMO_EXPORT int gsKnotVector_size(gsCKnotVector * kv)
 {
-    return reinterpret_cast<gismo::gsKnotVector<double>*>(kv)->size();
+    return reinterpret_cast<gsKnotVector<double>*>(kv)->size();
 }
 
 GISMO_EXPORT int gsKnotVector_uSize(gsCKnotVector * kv)
 {
-    return reinterpret_cast<gismo::gsKnotVector<double>*>(kv)->uSize();
+    return reinterpret_cast<gsKnotVector<double>*>(kv)->uSize();
 }
 
 GISMO_EXPORT int gsKnotVector_numElements(gsCKnotVector * kv)
 {
-    return reinterpret_cast<gismo::gsKnotVector<double>*>(kv)->numElements();
+    return reinterpret_cast<gsKnotVector<double>*>(kv)->numElements();
 }
 
 // GISMO_EXPORT double* gsKnotVector_unique(gsCKnotVector * kv)
 // {
-//     return reinterpret_cast<gismo::gsKnotVector<double>*>(kv)->unique().data();
+//     return reinterpret_cast<gsKnotVector<double>*>(kv)->unique().data();
 // }
 
 //GISMO_EXPORT double * greville(void * kv)
 //{
-//    gismo::gsKnotVector<double>* _kv = static_cast<gismo::gsKnotVector<double>*>(kv);
+//    gsKnotVector<double>* _kv = static_cast<gsKnotVector<double>*>(kv);
 //}
 
+#ifdef __cplusplus
+}
+#endif
\ No newline at end of file
diff --git a/src/gsCL2Projection.cpp b/src/gsCL2Projection.cpp
new file mode 100644
index 0000000..36fee07
--- /dev/null
+++ b/src/gsCL2Projection.cpp
@@ -0,0 +1,33 @@
+
+#include <gismo.h>
+#include <gsCInterface/gsCTypes.h>
+#include <gsCInterface/gsMacros.h>
+#include <gsCInterface/gsCQuadRule.h>
+
+using namespace gismo;
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif
+
+GISMO_EXPORT double gsL2Projection_into( gsCFunctionSet * projectionBasis,
+                                         gsCMultiBasis * integrationBasis,
+                                         gsCMultiPatch * geometryMap,
+                                         gsCFunctionSet * sourceFunction,
+                                         gsCMatrix * coefs,
+                                         gsCOptionList * options)
+{
+    auto * projBasis_ptr = RICAST_F(projectionBasis);
+    auto * intBasis_ptr = RICAST_MB(integrationBasis);
+    auto * geomMap_ptr = RICAST_MP(geometryMap);
+    auto * sourceFunc_ptr = RICAST_F(sourceFunction);
+    auto * coefs_ptr = RICAST_M(coefs);
+    auto * options_ptr = reinterpret_cast<gsOptionList*>(options);
+    double error= gsL2Projection<double>::project(*projBasis_ptr, *intBasis_ptr, *geomMap_ptr, *sourceFunc_ptr, *coefs_ptr, *options_ptr);
+    return error;
+}
+
+#ifdef __cplusplus
+}
+#endif
\ No newline at end of file
diff --git a/src/gsCL2Projection.h b/src/gsCL2Projection.h
new file mode 100644
index 0000000..e9d71e2
--- /dev/null
+++ b/src/gsCL2Projection.h
@@ -0,0 +1,16 @@
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif
+
+    GISMO_EXPORT double gsL2Projection_into( gsCFunctionSet * projectionBasis,
+                                             gsCMultiBasis * integrationBasis,
+                                             gsCMultiPatch * geometryMap,
+                                             gsCFunctionSet * sourceFunction,
+                                             gsCMatrix * coefs,
+                                             gsCOptionList * options);
+
+#ifdef __cplusplus
+}
+#endif
diff --git a/src/gsCMatrix.cpp b/src/gsCMatrix.cpp
index 294864b..a97f162 100644
--- a/src/gsCMatrix.cpp
+++ b/src/gsCMatrix.cpp
@@ -3,14 +3,21 @@
 #include <gsCInterface/gsMacros.h>
 #include <gsCInterface/gsCMatrix.h>
 
+using namespace gismo;
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif
+
 GISMO_EXPORT gsCMatrix * gsMatrix_create(void)
-{ return RICAST_CM(new gismo::gsMatrix<double>()); }
+{ return RICAST_CM(new gsMatrix<double>()); }
 
 GISMO_EXPORT gsCMatrix * gsMatrix_create_rc(int rows, int cols)
-{ return RICAST_CM(new gismo::gsMatrix<double>(rows,cols)); }
+{ return RICAST_CM(new gsMatrix<double>(rows,cols)); }
 
 GISMO_EXPORT gsCMatrix * gsMatrix_create_rcd(int rows, int cols, double * data)
-{ return RICAST_CM(new gismo::gsMatrix<double>(gismo::gsAsMatrix<double>(data,rows,cols))); }
+{ return RICAST_CM(new gsMatrix<double>(gsAsMatrix<double>(data,rows,cols))); }
 
 GISMO_EXPORT void gsMatrix_delete(gsCMatrix * m)
 { delete RICAST_M(m); }
@@ -32,3 +39,7 @@ GISMO_EXPORT int gsMatrix_cols(gsCMatrix * m)
 
 GISMO_EXPORT void gsMatrix_setZero(gsCMatrix * m)
 { RICAST_M(m)->setZero(); }
+
+#ifdef __cplusplus
+}
+#endif
\ No newline at end of file
diff --git a/src/gsCMatrixInt.cpp b/src/gsCMatrixInt.cpp
index 479225b..8dca053 100644
--- a/src/gsCMatrixInt.cpp
+++ b/src/gsCMatrixInt.cpp
@@ -3,14 +3,21 @@
 #include <gsCInterface/gsMacros.h>
 #include <gsCInterface/gsCMatrixInt.h>
 
+using namespace gismo;
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif
+
 GISMO_EXPORT gsCMatrixInt * gsMatrixInt_create(void)
-{ return reinterpret_cast<gsCMatrixInt*>(new gismo::gsMatrix<int>()); }
+{ return reinterpret_cast<gsCMatrixInt*>(new gsMatrix<int>()); }
 
 GISMO_EXPORT gsCMatrixInt * gsMatrixInt_create_rc(int rows, int cols)
-{ return reinterpret_cast<gsCMatrixInt*>(new gismo::gsMatrix<int>(rows,cols)); }
+{ return reinterpret_cast<gsCMatrixInt*>(new gsMatrix<int>(rows,cols)); }
 
 GISMO_EXPORT gsCMatrixInt * gsMatrixInt_create_rcd(int rows, int cols, int * data)
-{ return reinterpret_cast<gsCMatrixInt*>(new gismo::gsMatrix<int>(gismo::gsAsMatrix<int>(data,rows,cols))); }
+{ return reinterpret_cast<gsCMatrixInt*>(new gsMatrix<int>(gsAsMatrix<int>(data,rows,cols))); }
 
 GISMO_EXPORT void gsMatrixInt_delete(gsCMatrixInt * m)
 { delete RICAST_Mi(m); }
@@ -32,3 +39,7 @@ GISMO_EXPORT int gsMatrixInt_cols(gsCMatrixInt * m)
 
 GISMO_EXPORT void gsMatrixInt_setZero(gsCMatrixInt * m)
 { RICAST_Mi(m)->setZero(); }
+
+#ifdef __cplusplus
+}
+#endif
\ No newline at end of file
diff --git a/src/gsCMultiBasis.cpp b/src/gsCMultiBasis.cpp
index 5148096..7d1f604 100644
--- a/src/gsCMultiBasis.cpp
+++ b/src/gsCMultiBasis.cpp
@@ -5,15 +5,28 @@
 #include <gsCInterface/gsCMemory.h>
 #include <gsCInterface/gsMacros.h>
 
+using namespace gismo;
+
 #ifdef __cplusplus
 extern "C"
 {
 #endif
 
-GISMO_EXPORT gsCMultiBasis * gsMultiBasis_create(gsCMultiPatch * mp)
+GISMO_EXPORT gsCMultiBasis * gsMultiBasis_create()
+{
+    return RICAST_CMP(new gsMultiBasis<double>());
+}
+
+GISMO_EXPORT gsCMultiBasis * gsMultiBasis_create_basis(gsCBasis * basis)
+{
+    auto * basis_ptr = RICAST_B(basis);
+    return RICAST_CMP(new gsMultiBasis<double>(*basis_ptr));
+}
+
+GISMO_EXPORT gsCMultiBasis * gsMultiBasis_create_patches(gsCMultiPatch * mp)
 {
     auto * mp_ptr = RICAST_MP(mp);
-    return RICAST_CMP(new gismo::gsMultiBasis<double>(*mp_ptr));
+    return RICAST_CMP(new gsMultiBasis<double>(*mp_ptr));
 }
 
 GISMO_EXPORT void gsMultiBasis_addBasis(gsCMultiBasis* mb, gsCBasis* basis)
@@ -29,6 +42,29 @@ GISMO_EXPORT gsCBasis * gsMultiBasis_basis(gsCMultiBasis * mb, int i)
     return RICAST_CB(&mb_ptr->basis(i));
 }
 
+GISMO_EXPORT gsCMultiBasis * gsMultiBasis_read(char* filename)
+{
+    gsFileData<> data(filename);
+    if (data.hasAny< gsMultiBasis<> >())
+    {
+        gsMultiBasis<>::uPtr ptr = data.getAnyFirst< gsMultiBasis<> >();
+        return RICAST_CMB(ptr.release());
+    }
+    else
+    {
+        gsWarn<<"[G+Smo] No gsMultiBasis found in file "<<filename<<"\n";
+        return NULL;
+    }
+}
+
+GISMO_EXPORT void gsMultiBasis_write(gsCMultiBasis * obj, char* filename)
+{
+    gsFileData<> data;
+    data.add(*RICAST_MB(obj));
+    data.save(filename);
+}
+
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/src/gsCMultiBasis.h b/src/gsCMultiBasis.h
index 8a610f8..c346a04 100644
--- a/src/gsCMultiBasis.h
+++ b/src/gsCMultiBasis.h
@@ -3,12 +3,18 @@
 extern "C"
 {
 #endif
-    GISMO_EXPORT gsCMultiBasis * gsMultiBasis_create(gsCMultiPatch* mp);
+
+    GISMO_EXPORT gsCMultiBasis * gsMultiBasis_create();
+    GISMO_EXPORT gsCMultiBasis * gsMultiBasis_create_basis(gsCBasis * basis);
+    GISMO_EXPORT gsCMultiBasis * gsMultiBasis_create_patches(gsCMultiPatch* mp);
 
     GISMO_EXPORT void gsMultiBasis_addBasis(gsCMultiBasis* mp, gsCGeometry* geo);
 
     GISMO_EXPORT gsCBasis * gsMultiBasis_basis(gsCMultiBasis * mp, int i);
 
+    GISMO_EXPORT gsCMultiBasis * gsMultiBasis_read(char* filename);
+    GISMO_EXPORT void            gsMultiBasis_write(gsCMultiBasis * obj, char* filename);
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/src/gsCMultiPatch.cpp b/src/gsCMultiPatch.cpp
index 392764d..22b5061 100644
--- a/src/gsCMultiPatch.cpp
+++ b/src/gsCMultiPatch.cpp
@@ -5,6 +5,8 @@
 #include <gsCInterface/gsCMemory.h>
 #include <gsCInterface/gsMacros.h>
 
+using namespace gismo;
+
 #ifdef __cplusplus
 extern "C"
 {
@@ -12,7 +14,13 @@ extern "C"
 
 GISMO_EXPORT gsCMultiPatch* gsMultiPatch_create()
 {
-    return RICAST_CMP(new gismo::gsMultiPatch<double>());
+    return RICAST_CMP(new gsMultiPatch<double>());
+}
+
+GISMO_EXPORT gsCMultiPatch* gsMultiPatch_create_geometry(gsCGeometry * geo)
+{
+    auto * geo_ptr = RICAST_G(geo);
+    return RICAST_CMP(new gsMultiPatch<double>(*geo_ptr));
 }
 
 GISMO_EXPORT void gsMultiPatch_addPatch(gsCMultiPatch* mp, gsCGeometry* geo)
@@ -40,6 +48,46 @@ GISMO_EXPORT void gsMultiPatch_computeTopology(gsCMultiPatch * mp)
     mp_ptr->computeTopology();
 }
 
+GISMO_EXPORT void gsMultiPatch_embed(gsCMultiPatch * mp, int dim)
+{
+    auto * mp_ptr = RICAST_MP(mp);
+    mp_ptr->embed(dim);
+}
+
+GISMO_EXPORT void gsMultiPatch_uniformRefine(gsCMultiPatch * mp, int numKnots, int mul, int dir)
+{
+    auto * mp_ptr = RICAST_MP(mp);
+    mp_ptr->uniformRefine(numKnots, mul, dir);
+}
+
+GISMO_EXPORT void gsMultiPatch_degreeElevate(gsCMultiPatch * mp, int i, int dir)
+{
+    auto * mp_ptr = RICAST_MP(mp);
+    mp_ptr->degreeElevate(i, dir);
+}
+
+GISMO_EXPORT gsCMultiPatch * gsMultiPatch_read(char* filename)
+{
+    gsFileData<> data(filename);
+    if (data.hasAny< gsMultiPatch<> >())
+    {
+        gsMultiPatch<>::uPtr ptr = data.getAnyFirst< gsMultiPatch<> >();
+        return RICAST_CMP(ptr.release());
+    }
+    else
+    {
+        gsWarn<<"[G+Smo] No gsMultiPatch found in file "<<filename<<"\n";
+        return NULL;
+    }
+}
+
+GISMO_EXPORT void gsMultiPatch_write(gsCMultiPatch * obj, char* filename)
+{
+    gsFileData<> data;
+    data.add(*RICAST_MP(obj));
+    data.save(filename);
+}
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/src/gsCMultiPatch.h b/src/gsCMultiPatch.h
index 1928d2c..b6aa916 100644
--- a/src/gsCMultiPatch.h
+++ b/src/gsCMultiPatch.h
@@ -4,6 +4,7 @@ extern "C"
 {
 #endif
     GISMO_EXPORT gsCMultiPatch* gsMultiPatch_create();
+    GISMO_EXPORT gsCMultiPatch* gsMultiPatch_create_geometry(gsCGeometry * geo);
 
     GISMO_EXPORT void gsMultiPatch_addPatch(gsCMultiPatch* mp, gsCGeometry* geo);
 
@@ -12,6 +13,13 @@ extern "C"
     GISMO_EXPORT gsCGeometry * gsMultiPatch_patch(gsCMultiPatch * mp, int i);
 
     GISMO_EXPORT void gsMultiPatch_computeTopology(gsCMultiPatch * mp);
+    GISMO_EXPORT void gsMultiPatch_embed(gsCMultiPatch * mp, int dim);
+
+    GISMO_EXPORT void gsMultiPatch_uniformRefine(gsCMultiPatch * mp, int numKnots, int mul, int dir);
+    GISMO_EXPORT void gsMultiPatch_degreeElevate(gsCMultiPatch * mp, int i, int dir);
+
+    GISMO_EXPORT gsCMultiPatch * gsMultiPatch_read(char* filename);
+    GISMO_EXPORT void            gsMultiPatch_write(gsCMultiPatch * obj, char* filename);
 
 #ifdef __cplusplus
 }
diff --git a/src/gsCNurbsCreator.cpp b/src/gsCNurbsCreator.cpp
index 9495169..17b8261 100644
--- a/src/gsCNurbsCreator.cpp
+++ b/src/gsCNurbsCreator.cpp
@@ -7,6 +7,8 @@
 #include <gsCInterface/gsCBasis.h>
 #include <gsCInterface/gsMacros.h>
 
+using namespace gismo;
+
 #ifdef __cplusplus
 extern "C"
 {
@@ -14,72 +16,72 @@ extern "C"
 
 GISMO_EXPORT gsCGeometry* gsNurbsCreator_BSplineUnitInterval(int deg)
 {
-    return RICAST_CG(gismo::gsNurbsCreator<double>::BSplineUnitInterval(deg).release());
+    return RICAST_CG(gsNurbsCreator<double>::BSplineUnitInterval(deg).release());
 }
 
 GISMO_EXPORT gsCGeometry* gsNurbsCreator_BSplineRectangle(double low_x, double low_y, double upp_x, double upp_y, double turndeg)
 {
-    return RICAST_CG(gismo::gsNurbsCreator<double>::BSplineRectangle(low_x, low_y, upp_x, upp_y, turndeg).release());
+    return RICAST_CG(gsNurbsCreator<double>::BSplineRectangle(low_x, low_y, upp_x, upp_y, turndeg).release());
 }
 
 GISMO_EXPORT gsCGeometry* gsNurbsCreator_BSplineTrapezium(double Lbot, double Ltop, double H, double d, double turndeg)
 {
-    return RICAST_CG(gismo::gsNurbsCreator<double>::BSplineTrapezium(Lbot, Ltop, H, d, turndeg).release());
+    return RICAST_CG(gsNurbsCreator<double>::BSplineTrapezium(Lbot, Ltop, H, d, turndeg).release());
 }
 
 GISMO_EXPORT gsCGeometry* gsNurbsCreator_BSplineSquare(double r, double x, double y)
 {
-    return RICAST_CG(gismo::gsNurbsCreator<double>::BSplineSquare(r, x, y).release());
+    return RICAST_CG(gsNurbsCreator<double>::BSplineSquare(r, x, y).release());
 }
 
 GISMO_EXPORT gsCMultiPatch* gsNurbsCreator_BSplineSquareGrid(int n, int m, double r, double lx, double ly)
 {
-    return RICAST_CMP(new gismo::gsMultiPatch<double>(gismo::gsNurbsCreator<double>::BSplineSquareGrid(n, m, r, lx, ly)));
+    return RICAST_CMP(new gsMultiPatch<double>(gsNurbsCreator<double>::BSplineSquareGrid(n, m, r, lx, ly)));
 }
 
 GISMO_EXPORT gsCGeometry* gsNurbsCreator_BSplineCube(double r, double x, double y, double z)
 {
-    return RICAST_CG(gismo::gsNurbsCreator<double>::BSplineCube(r, x, y, z).release());
+    return RICAST_CG(gsNurbsCreator<double>::BSplineCube(r, x, y, z).release());
 }
 
 GISMO_EXPORT gsCMultiPatch* gsNurbsCreator_BSplineCubeGrid(int n, int m, int p, double r, double lx, double ly, double lz)
 {
-    return RICAST_CMP(new gismo::gsMultiPatch<double>(gismo::gsNurbsCreator<double>::BSplineCubeGrid(n, m, p, r, lx, ly, lz)));
+    return RICAST_CMP(new gsMultiPatch<double>(gsNurbsCreator<double>::BSplineCubeGrid(n, m, p, r, lx, ly, lz)));
 }
 
 GISMO_EXPORT gsCGeometry* gsNurbsCreator_NurbsQuarterAnnulus(double r0, double r1)
 {
-    return RICAST_CG(gismo::gsNurbsCreator<double>::NurbsQuarterAnnulus(r0, r1).release());
+    return RICAST_CG(gsNurbsCreator<double>::NurbsQuarterAnnulus(r0, r1).release());
 }
 
 GISMO_EXPORT gsCGeometry* gsNurbsCreator_NurbsAnnulus(double r0, double r1)
 {
-    return RICAST_CG(gismo::gsNurbsCreator<double>::NurbsAnnulus(r0, r1).release());
+    return RICAST_CG(gsNurbsCreator<double>::NurbsAnnulus(r0, r1).release());
 }
 
 GISMO_EXPORT gsCGeometry* gsNurbsCreator_BSplineSaddle()
 {
-    return RICAST_CG(gismo::gsNurbsCreator<double>::BSplineSaddle().release());
+    return RICAST_CG(gsNurbsCreator<double>::BSplineSaddle().release());
 }
 
 GISMO_EXPORT gsCGeometry* gsNurbsCreator_NurbsSphere(double r, double x, double y, double z)
 {
-    return RICAST_CG(gismo::gsNurbsCreator<double>::NurbsSphere(r, x, y, z).release());
+    return RICAST_CG(gsNurbsCreator<double>::NurbsSphere(r, x, y, z).release());
 }
 
 GISMO_EXPORT gsCGeometry* gsNurbsCreator_NurbsCircle(double r, double x, double y)
 {
-    return RICAST_CG(gismo::gsNurbsCreator<double>::NurbsCircle(r, x, y).release());
+    return RICAST_CG(gsNurbsCreator<double>::NurbsCircle(r, x, y).release());
 }
 
 GISMO_EXPORT gsCGeometry* gsNurbsCreator_BSplineTriangle(double H, double W)
 {
-    return RICAST_CG(gismo::gsNurbsCreator<double>::BSplineTriangle(H, W).release());
+    return RICAST_CG(gsNurbsCreator<double>::BSplineTriangle(H, W).release());
 }
 
 GISMO_EXPORT gsCMultiPatch* gsNurbsCreator_BSplineStar(int N, double R0, double R1)
 {
-    return RICAST_CMP(new gismo::gsMultiPatch<double>(gismo::gsNurbsCreator<double>::BSplineStar(N, R0, R1)));
+    return RICAST_CMP(new gsMultiPatch<double>(gsNurbsCreator<double>::BSplineStar(N, R0, R1)));
 }
 
 #ifdef __cplusplus
diff --git a/src/gsCOptionList.cpp b/src/gsCOptionList.cpp
index 8ce3234..c1dc9d5 100644
--- a/src/gsCOptionList.cpp
+++ b/src/gsCOptionList.cpp
@@ -4,6 +4,8 @@
 #include <gsCInterface/gsMacros.h>
 #include <gsCInterface/gsCOptionList.h>
 
+using namespace gismo;
+
 #ifdef __cplusplus
 extern "C"
 {
@@ -11,77 +13,77 @@ extern "C"
 
 GISMO_EXPORT gsCOptionList * gsOptionList_create()
 {
-    return reinterpret_cast<gsCOptionList*>(new gismo::gsOptionList());
+    return reinterpret_cast<gsCOptionList*>(new gsOptionList());
 }
 
 GISMO_EXPORT void gsOptionList_delete(gsCOptionList * ol)
 {
-    delete reinterpret_cast<gismo::gsOptionList*>(ol);
+    delete reinterpret_cast<gsOptionList*>(ol);
 }
 
 GISMO_EXPORT void gsOptionList_print(gsCOptionList * ol)
 {
-    reinterpret_cast<gismo::gsOptionList*>(ol)->print(gsInfo);
+    reinterpret_cast<gsOptionList*>(ol)->print(gsInfo);
 }
 
 GISMO_EXPORT void gsOptionList_addString(gsCOptionList * ol, const char * label, const char * description, const char * value)
 {
-    reinterpret_cast<gismo::gsOptionList*>(ol)->addString(label, description, value);
+    reinterpret_cast<gsOptionList*>(ol)->addString(label, description, value);
 }
 
 GISMO_EXPORT void gsOptionList_addInt(gsCOptionList * ol, const char * label, const char * description, int value)
 {
-    reinterpret_cast<gismo::gsOptionList*>(ol)->addInt(label, description, value);
+    reinterpret_cast<gsOptionList*>(ol)->addInt(label, description, value);
 }
 
 GISMO_EXPORT void gsOptionList_addReal(gsCOptionList * ol, const char * label, const char * description, double value)
 {
-    reinterpret_cast<gismo::gsOptionList*>(ol)->addReal(label, description, value);
+    reinterpret_cast<gsOptionList*>(ol)->addReal(label, description, value);
 }
 
 GISMO_EXPORT void gsOptionList_addSwitch(gsCOptionList * ol, const char * label, const char * description, int value)
 {
-    reinterpret_cast<gismo::gsOptionList*>(ol)->addSwitch(label, description, value);
+    reinterpret_cast<gsOptionList*>(ol)->addSwitch(label, description, value);
 }
 
 GISMO_EXPORT void gsOptionList_setString(gsCOptionList * ol, const char * label, const char * value)
 {
-    reinterpret_cast<gismo::gsOptionList*>(ol)->setString(label, value);
+    reinterpret_cast<gsOptionList*>(ol)->setString(label, value);
 }
 
 GISMO_EXPORT void gsOptionList_setInt(gsCOptionList * ol, const char * label, int value)
 {
-    reinterpret_cast<gismo::gsOptionList*>(ol)->setInt(label, value);
+    reinterpret_cast<gsOptionList*>(ol)->setInt(label, value);
 }
 
 GISMO_EXPORT void gsOptionList_setReal(gsCOptionList * ol, const char * label, double value)
 {
-    reinterpret_cast<gismo::gsOptionList*>(ol)->setReal(label, value);
+    reinterpret_cast<gsOptionList*>(ol)->setReal(label, value);
 }
 
 GISMO_EXPORT void gsOptionList_setSwitch(gsCOptionList * ol, const char * label, int value)
 {
-    reinterpret_cast<gismo::gsOptionList*>(ol)->setSwitch(label, value);
+    reinterpret_cast<gsOptionList*>(ol)->setSwitch(label, value);
 }
 
 GISMO_EXPORT char* gsOptionList_getString(gsCOptionList * ol, const char * label)
 {
-    return const_cast<char*>(reinterpret_cast<gismo::gsOptionList*>(ol)->getString(label).c_str());
+    return const_cast<char*>(reinterpret_cast<gsOptionList*>(ol)->getString(label).c_str());
 }
 
 GISMO_EXPORT int gsOptionList_getInt(gsCOptionList * ol, const char * label)
 {
-    return reinterpret_cast<gismo::gsOptionList*>(ol)->getInt(label);
+    return reinterpret_cast<gsOptionList*>(ol)->getInt(label);
 }
 
 GISMO_EXPORT double gsOptionList_getReal(gsCOptionList * ol, const char * label)
 {
-    return reinterpret_cast<gismo::gsOptionList*>(ol)->getReal(label);
+    return reinterpret_cast<gsOptionList*>(ol)->getReal(label);
 }
 
 GISMO_EXPORT int gsOptionList_getSwitch(gsCOptionList * ol, const char * label)
 {
-    return reinterpret_cast<gismo::gsOptionList*>(ol)->getSwitch(label);
+    return reinterpret_cast<gsOptionList*>(ol)->getSwitch(label);
 }
 
 #ifdef __cplusplus
diff --git a/src/gsCQuadRule.cpp b/src/gsCQuadRule.cpp
index 00836ae..afd3359 100644
--- a/src/gsCQuadRule.cpp
+++ b/src/gsCQuadRule.cpp
@@ -4,6 +4,8 @@
 #include <gsCInterface/gsMacros.h>
 #include <gsCInterface/gsCQuadRule.h>
 
+using namespace gismo;
+
 #ifdef __cplusplus
 extern "C"
 {
@@ -12,45 +14,45 @@ extern "C"
 GISMO_EXPORT gsCQuadRule* gsQuadrature_get(gsCBasis * basis, gsCOptionList* options, int fixDir)
 {
     auto * basis_ptr = RICAST_B(basis);
-    auto * options_ptr = reinterpret_cast<gismo::gsOptionList*>(options);
-    return reinterpret_cast<gsCQuadRule*>(gismo::gsQuadrature::getPtr(*basis_ptr, *options_ptr, fixDir).release());
+    auto * options_ptr = reinterpret_cast<gsOptionList*>(options);
+    return reinterpret_cast<gsCQuadRule*>(gsQuadrature::getPtr(*basis_ptr, *options_ptr, fixDir).release());
 }
 
 GISMO_EXPORT gsCQuadRule* gsGaussRule_create(int d, int* numNodes, int digits)
 {
     std::vector<int> nodes(numNodes, numNodes + d);
-    return reinterpret_cast<gsCQuadRule*>(new gismo::gsGaussRule<double>(gismo::gsAsVector<int>(nodes), digits));
+    return reinterpret_cast<gsCQuadRule*>(new gsGaussRule<double>(gsAsVector<int>(nodes), digits));
 }
 
 GISMO_EXPORT gsCQuadRule* gsLobattoRule_create(int d, int* numNodes, int digits)
 {
     std::vector<int> nodes(numNodes, numNodes + d);
-    return reinterpret_cast<gsCQuadRule*>(new gismo::gsLobattoRule<double>(gismo::gsAsVector<int>(nodes), digits));
+    return reinterpret_cast<gsCQuadRule*>(new gsLobattoRule<double>(gsAsVector<int>(nodes), digits));
 }
 
 GISMO_EXPORT void gsQuadRule_delete(gsCQuadRule * qr)
 {
-    delete reinterpret_cast<gismo::gsQuadRule<double>*>(qr);
+    delete reinterpret_cast<gsQuadRule<double>*>(qr);
 }
 
 GISMO_EXPORT int gsQuadRule_numNodes(gsCQuadRule * qr)
 {
-    return reinterpret_cast<gismo::gsQuadRule<double>*>(qr)->numNodes();
+    return reinterpret_cast<gsQuadRule<double>*>(qr)->numNodes();
 }
 
 GISMO_EXPORT int gsQuadRule_dim(gsCQuadRule * qr)
 {
-    return reinterpret_cast<gismo::gsQuadRule<double>*>(qr)->dim();
+    return reinterpret_cast<gsQuadRule<double>*>(qr)->dim();
 }
 
 GISMO_EXPORT void gsQuadRule_mapTo(gsCQuadRule* qr, gsCVector* lower, gsCVector* upper, gsCMatrix* nodes, gsCVector* weights)
 {
-    reinterpret_cast<gismo::gsQuadRule<double>*>(qr)->mapTo(*RICAST_V(lower), *RICAST_V(upper), *RICAST_M(nodes), *RICAST_V(weights));
+    reinterpret_cast<gsQuadRule<double>*>(qr)->mapTo(*RICAST_V(lower), *RICAST_V(upper), *RICAST_M(nodes), *RICAST_V(weights));
 }
 
 GISMO_EXPORT void gsQuadRule_mapToScalar(gsCQuadRule* qr, double startVal, double endVal, gsCMatrix* nodes, gsCVector* weights)
 {
-    reinterpret_cast<gismo::gsQuadRule<double>*>(qr)->mapTo(startVal, endVal, *RICAST_M(nodes), *RICAST_V(weights));
+    reinterpret_cast<gsQuadRule<double>*>(qr)->mapTo(startVal, endVal, *RICAST_M(nodes), *RICAST_V(weights));
 }
 
 
diff --git a/src/gsCQuasiInterpolate.cpp b/src/gsCQuasiInterpolate.cpp
new file mode 100644
index 0000000..b14768c
--- /dev/null
+++ b/src/gsCQuasiInterpolate.cpp
@@ -0,0 +1,48 @@
+
+#include <gismo.h>
+#include <gsCInterface/gsCTypes.h>
+#include <gsCInterface/gsMacros.h>
+#include <gsCInterface/gsCQuadRule.h>
+
+using namespace gismo;
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif
+
+
+GISMO_EXPORT void gsQuasiInterpolate_localIntpl_into( gsCBasis * basis,
+                                                      gsCFunctionSet * fun,
+                                                      gsCMatrix * result)
+{
+    auto basis_ptr = RICAST_B(basis);
+    auto fun_ptr = RICAST_F(fun);
+    auto result_ptr = RICAST_M(result);
+    gsQuasiInterpolate<double>::localIntpl(*basis_ptr, fun_ptr->function(0), *result_ptr);
+}
+
+GISMO_EXPORT void gsQuasiInterpolate_Taylor_into( gsCBasis * basis,
+                                                  gsCFunctionSet * fun,
+                                                  int deg,
+                                                  gsCMatrix * result)
+{
+    auto basis_ptr = RICAST_B(basis);
+    auto fun_ptr = RICAST_F(fun);
+    auto result_ptr = RICAST_M(result);
+    gsQuasiInterpolate<double>::Taylor(*basis_ptr, fun_ptr->function(0), deg, *result_ptr);
+}
+
+GISMO_EXPORT void gsQuasiInterpolate_Schoenberg_into( gsCBasis * basis,
+                                                      gsCFunctionSet * fun,
+                                                      gsCMatrix * result)
+{
+    auto basis_ptr = RICAST_B(basis);
+    auto fun_ptr = RICAST_F(fun);
+    auto result_ptr = RICAST_M(result);
+    gsQuasiInterpolate<double>::Schoenberg(*basis_ptr, fun_ptr->function(0), *result_ptr);
+}
+
+#ifdef __cplusplus
+}
+#endif
\ No newline at end of file
diff --git a/src/gsCQuasiInterpolate.h b/src/gsCQuasiInterpolate.h
new file mode 100644
index 0000000..c9c3a10
--- /dev/null
+++ b/src/gsCQuasiInterpolate.h
@@ -0,0 +1,22 @@
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif
+
+    GISMO_EXPORT void gsQuasiInterpolate_localIntpl_into( gsCBasis * basis,
+                                                          gsCFunctionSet * fun,
+                                                          gsCMatrix * result);
+
+    GISMO_EXPORT void gsQuasiInterpolate_Taylor_into( gsCBasis * basis,
+                                                      gsCFunctionSet * fun,
+                                                      int deg,
+                                                      gsCMatrix * result);
+
+    GISMO_EXPORT void gsQuasiInterpolate_Schoenberg_into( gsCBasis * basis,
+                                                          gsCFunctionSet * fun,
+                                                          gsCMatrix * result);
+
+#ifdef __cplusplus
+}
+#endif
diff --git a/src/gsCReadFile.cpp b/src/gsCReadFile.cpp
index 69f0fed..4151646 100644
--- a/src/gsCReadFile.cpp
+++ b/src/gsCReadFile.cpp
@@ -1,31 +1,39 @@
 #include <gismo.h>
 #include <gsCInterface/gsCReadFile.h>
+#include <gsCInterface/gsCTypes.h>
+#include <gsCInterface/gsMacros.h>
 
+using namespace gismo;
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif
 GISMO_EXPORT void* gsCReadFile(char* filename)
 {
-    gsInfo << "Loading file: " << filename << std::endl;
+    gsInfo << "[G+Smo] Loading file: " << filename << std::endl;
 
     void* result = NULL;
 
-    gismo::gsFileData<> data(filename);
-    if      ( data.hasAny< gismo::gsGeometry<>  >() )
+    gsFileData<> data(filename);
+    if      ( data.hasAny< gsGeometry<>  >() )
     {
-        gismo::gsGeometry<>::uPtr ptr = data.getAnyFirst< gismo::gsGeometry<>  >();
+        gsGeometry<>::uPtr ptr = data.getAnyFirst< gsGeometry<>  >();
         result = ptr.release();
     }
-    else if ( data.hasAny< gismo::gsBasis<>  >() )
+    else if ( data.hasAny< gsBasis<>  >() )
     {
-        gismo::gsBasis<>::uPtr ptr = data.getAnyFirst< gismo::gsBasis<>  >();
+        gsBasis<>::uPtr ptr = data.getAnyFirst< gsBasis<>  >();
         result = ptr.release();
     }
-    else if ( data.hasAny< gismo::gsMultiPatch<>  >() )
+    else if ( data.hasAny< gsMultiPatch<>  >() )
     {
-        gismo::gsMultiPatch<>::uPtr ptr = data.getAnyFirst< gismo::gsMultiPatch<>  >();
+        gsMultiPatch<>::uPtr ptr = data.getAnyFirst< gsMultiPatch<>  >();
         result = ptr.release();
     }
-    else if ( data.hasAny< gismo::gsMultiBasis<>  >() )
+    else if ( data.hasAny< gsMultiBasis<>  >() )
     {
-        gismo::gsMultiBasis<>::uPtr ptr = data.getAnyFirst< gismo::gsMultiBasis<>  >();
+        gsMultiBasis<>::uPtr ptr = data.getAnyFirst< gsMultiBasis<>  >();
         result = ptr.release();
     }
     else
@@ -35,3 +43,7 @@ GISMO_EXPORT void* gsCReadFile(char* filename)
     }
     return result;
 }
+
+#ifdef __cplusplus
+}
+#endif
\ No newline at end of file
diff --git a/src/gsCReadFile.h b/src/gsCReadFile.h
index 46ea51c..960645d 100644
--- a/src/gsCReadFile.h
+++ b/src/gsCReadFile.h
@@ -1,9 +1,12 @@
+#include <gsCInterface/gsCTypes.h>
+
 #ifdef __cplusplus
 extern "C"
 {
 #endif
 
-    GISMO_EXPORT void* gsCReadFile(char* filename);
+    GISMO_EXPORT void* gsCReadFile (char* filename);
+    GISMO_EXPORT void  gsCWriteFile(gsCFunctionSet * obj, char * filename);
 
 #ifdef __cplusplus
 }
diff --git a/src/gsCSparseMatrix.cpp b/src/gsCSparseMatrix.cpp
index eb274c0..240dcd8 100644
--- a/src/gsCSparseMatrix.cpp
+++ b/src/gsCSparseMatrix.cpp
@@ -3,8 +3,15 @@
 #include <gsCInterface/gsMacros.h>
 #include <gsCInterface/gsCSparseMatrix.h>
 
+using namespace gismo;
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif
+
 GISMO_EXPORT gsCSparseMatrix * gsSparseMatrix_create(void)
-{ return RICAST_CSM(new gismo::gsMatrix<double>()); }
+{ return RICAST_CSM(new gsSparseMatrix<double>()); }
 
 GISMO_EXPORT void gsSparseMatrix_delete(gsCSparseMatrix * m)
 { delete RICAST_SM(m); }
@@ -15,37 +22,65 @@ GISMO_EXPORT void gsSparseMatrix_print(gsCSparseMatrix * m)
 // GISMO_EXPORT double * gsSparseMatrix_data(gsCSparseMatrix * m)
 // { return RICAST_SM(m)->data(); }
 
-GISMO_EXPORT void gsSparseMatrix_setFromTriplets(gsCSparseMatrix * m, int * rows, int * cols, double * values, int nnz)
+GISMO_EXPORT double* gsSparseMatrix_valuePtr(gsCSparseMatrix * m)
+{ return RICAST_SM(m)->valuePtr(); } // get pointer to matrix values
+GISMO_EXPORT int*    gsSparseMatrix_innerIndexPtr(gsCSparseMatrix * m)
+{ return RICAST_SM(m)->innerIndexPtr(); } // get pointer to matrix rows
+GISMO_EXPORT int*    gsSparseMatrix_outerIndexPtr(gsCSparseMatrix * m)
+{ return RICAST_SM(m)->outerIndexPtr(); } // get pointer to matrix columns
+
+GISMO_EXPORT int gsSparseMatrix_rows(gsCSparseMatrix * m)
+{ return RICAST_SM(m)->rows(); }
+
+GISMO_EXPORT int gsSparseMatrix_cols(gsCSparseMatrix * m)
+{ return RICAST_SM(m)->cols(); }
+
+GISMO_EXPORT int gsSparseMatrix_nnz(gsCSparseMatrix * m)
+{ return RICAST_SM(m)->nonZeros(); }
+
+GISMO_EXPORT void gsSparseMatrix_setFromTriplets(gsCSparseMatrix * m, gsCVectorInt * rows, gsCVectorInt * cols, gsCVector * values)
 {
-    std::vector<int> rows_vec(rows, rows + nnz);
-    std::vector<int> cols_vec(cols, cols + nnz);
-    std::vector<double> values_vec(values, values + nnz);
+    auto * R = RICAST_Vi(rows);
+    auto * C = RICAST_Vi(cols);
+    auto * V = RICAST_V(values);
+
+    GISMO_ENSURE(R->size() == C->size() && R->size() == V->size(), "Input vectors must have the same size.");
 
-    gismo::gsSparseEntries<double> entries;
-    entries.reserve(nnz);
-    for (int i = 0; i < nnz; i++)
-        entries.add(rows_vec[i], cols_vec[i], values_vec[i]);
+    gsSparseEntries<double> entries;
+    entries.reserve(R->size());
+    for (int i = 0; i < R->size(); i++)
+        entries.add((*R)[i], (*C)[i], (*V)[i]);
 
+    RICAST_SM(m)->resize(R->size(), C->size());
     RICAST_SM(m)->setFrom(entries);
 }
 
-GISMO_EXPORT void gsSparseMatrix_intoTriplets(gsCSparseMatrix * m, int * rows, int * cols, double * values)
+GISMO_EXPORT void gsSparseMatrix_intoTriplets(gsCSparseMatrix * m, gsCVectorInt * rows, gsCVectorInt * cols, gsCVector * vals)
 {
     auto * sm = RICAST_SM(m);
+    auto * R = RICAST_Vi(rows);
+    auto * C = RICAST_Vi(cols);
+    auto * V = RICAST_V(vals);
+
+    R->resize(sm->nonZeros());
+    C->resize(sm->nonZeros());
+    V->resize(sm->nonZeros());
+    auto R_it = R->begin();
+    auto C_it = C->begin();
+    auto V_it = V->begin();
+
     for (int i = 0; i!=sm->outerSize(); i++)
     {
-        for (typename gismo::gsSparseMatrix<double>::InnerIterator it(*sm,i); it; ++it)
+        for (typename gsSparseMatrix<double>::InnerIterator it(*sm,i); it;
+                    ++it, ++R_it, ++C_it, ++V_it)
         {
-            *rows = it.row();
-            *cols = it.col();
-            *values = it.value();
-            ++rows; ++cols; ++values;
+            *R_it = it.row();
+            *C_it = it.col();
+            *V_it = it.value();
         }
     }
 }
 
-GISMO_EXPORT int gsSparseMatrix_rows(gsCSparseMatrix * m)
-{ return RICAST_SM(m)->rows(); }
-
-GISMO_EXPORT int gsSparseMatrix_cols(gsCSparseMatrix * m)
-{ return RICAST_SM(m)->cols(); }
+#ifdef __cplusplus
+}
+#endif
\ No newline at end of file
diff --git a/src/gsCSparseMatrix.h b/src/gsCSparseMatrix.h
index 20324ed..143821b 100644
--- a/src/gsCSparseMatrix.h
+++ b/src/gsCSparseMatrix.h
@@ -10,11 +10,17 @@ extern "C"
     GISMO_EXPORT void gsSparseMatrix_print(gsCSparseMatrix * m);
     // GISMO_EXPORT double* gsSparseMatrix_data(gsCSparseMatrix * m);   // get pointer to matrix data
 
-    GISMO_EXPORT void gsSparseMatrix_setFromTriplets(gsCSparseMatrix * m, int * rows, int * cols, double * values, int nnz);
-    GISMO_EXPORT void gsSparseMatrix_intoTriplets(gsCSparseMatrix * m, int * rows, int * cols, double * values);
+    GISMO_EXPORT double* gsSparseMatrix_valuePtr(gsCSparseMatrix * m); // get pointer to matrix values
+    GISMO_EXPORT int*    gsSparseMatrix_innerIndexPtr(gsCSparseMatrix * m);    // get pointer to matrix rows
+    GISMO_EXPORT int*    gsSparseMatrix_outerIndexPtr(gsCSparseMatrix * m);    // get pointer to matrix cols
+    GISMO_EXPORT int     gsSparseMatrix_rows(gsCSparseMatrix * m);
+    GISMO_EXPORT int     gsSparseMatrix_cols(gsCSparseMatrix * m);
+    GISMO_EXPORT int     gsSparseMatrix_nnz(gsCSparseMatrix * m);
+
+
+    GISMO_EXPORT void gsSparseMatrix_setFromTriplets(gsCSparseMatrix * m, gsCVectorInt * rows, gsCVectorInt * cols, gsCVector * values);
+    GISMO_EXPORT void gsSparseMatrix_intoTriplets(gsCSparseMatrix * m, gsCVectorInt * rows, gsCVectorInt * cols, gsCVector * vals);
 
-    GISMO_EXPORT int gsSparseMatrix_rows(gsCSparseMatrix * m);
-    GISMO_EXPORT int gsSparseMatrix_cols(gsCSparseMatrix * m);
 
 #ifdef __cplusplus
 }
diff --git a/src/gsCVector.cpp b/src/gsCVector.cpp
index b498bd4..c85d9ee 100644
--- a/src/gsCVector.cpp
+++ b/src/gsCVector.cpp
@@ -3,14 +3,21 @@
 #include <gsCInterface/gsMacros.h>
 #include <gsCInterface/gsCVector.h>
 
+using namespace gismo;
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif
+
 GISMO_EXPORT gsCVector * gsVector_create(void)
-{ return reinterpret_cast<gsCVector*>(new gismo::gsVector<double>()); }
+{ return reinterpret_cast<gsCVector*>(new gsVector<double>()); }
 
 GISMO_EXPORT gsCVector * gsVector_create_r(int rows)
-{ return reinterpret_cast<gsCVector*>(new gismo::gsVector<double>(rows)); }
+{ return reinterpret_cast<gsCVector*>(new gsVector<double>(rows)); }
 
 GISMO_EXPORT gsCVector * gsVector_create_rd(int rows, double * data)
-{ return reinterpret_cast<gsCVector*>(new gismo::gsVector<double>(gismo::gsAsVector<double>(data,rows))); }
+{ return reinterpret_cast<gsCVector*>(new gsVector<double>(gsAsVector<double>(data,rows))); }
 
 GISMO_EXPORT void gsVector_delete(gsCVector * m)
 { delete RICAST_V(m); }
@@ -32,3 +39,7 @@ GISMO_EXPORT int gsVector_cols(gsCVector * m)
 
 GISMO_EXPORT void gsVector_setZero(gsCVector * m)
 { RICAST_V(m)->setZero(); }
+
+#ifdef __cplusplus
+}
+#endif
\ No newline at end of file
diff --git a/src/gsCVectorInt.cpp b/src/gsCVectorInt.cpp
index f9e05c4..1d1174b 100644
--- a/src/gsCVectorInt.cpp
+++ b/src/gsCVectorInt.cpp
@@ -3,14 +3,21 @@
 #include <gsCInterface/gsCVectorInt.h>
 #include <gsCInterface/gsMacros.h>
 
+using namespace gismo;
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif
+
 GISMO_EXPORT gsCVectorInt * gsVectorInt_create(void)
-{ return reinterpret_cast<gsCVectorInt*>(new gismo::gsVector<int>()); }
+{ return reinterpret_cast<gsCVectorInt*>(new gsVector<int>()); }
 
 GISMO_EXPORT gsCVectorInt * gsVectorInt_create_r(int rows)
-{ return reinterpret_cast<gsCVectorInt*>(new gismo::gsVector<int>(rows)); }
+{ return reinterpret_cast<gsCVectorInt*>(new gsVector<int>(rows)); }
 
 GISMO_EXPORT gsCVectorInt * gsVectorInt_create_rd(int rows, int * data)
-{ return reinterpret_cast<gsCVectorInt*>(new gismo::gsVector<int>(gismo::gsAsVector<int>(data,rows))); }
+{ return reinterpret_cast<gsCVectorInt*>(new gsVector<int>(gsAsVector<int>(data,rows))); }
 
 GISMO_EXPORT void gsVectorInt_delete(gsCVectorInt * m)
 { delete RICAST_Vi(m); }
@@ -32,3 +39,7 @@ GISMO_EXPORT int gsVectorInt_cols(gsCVectorInt * m)
 
 GISMO_EXPORT void gsVectorInt_setZero(gsCVectorInt * m)
 { RICAST_Vi(m)->setZero(); }
+
+#ifdef __cplusplus
+}
+#endif
\ No newline at end of file