diff --git a/cmake/Gismo.jl.in b/cmake/Gismo.jl.in
index 7e67846..e21f9eb 100644
--- a/cmake/Gismo.jl.in
+++ b/cmake/Gismo.jl.in
@@ -8,6 +8,9 @@ include("Declarations.jl")
 # Load gsCore
 include("gsCore.jl")
 
+# Load gsMatrix
+include("gsMatrix.jl")
+
 # Load gsAssembler
 include("gsAssembler.jl")
 
@@ -17,9 +20,6 @@ include("gsHSplines.jl")
 # Load gsIO
 include("gsIO.jl")
 
-# Load gsMatrix
-include("gsMatrix.jl")
-
 # Load gsMatrix
 include("gsModeling.jl")
 
diff --git a/src/Gismo.jl b/src/Gismo.jl
index 94c79cc..57a5dff 100644
--- a/src/Gismo.jl
+++ b/src/Gismo.jl
@@ -8,6 +8,9 @@ include("Declarations.jl")
 # Load gsCore
 include("gsCore.jl")
 
+# Load gsMatrix
+include("gsMatrix.jl")
+
 # Load gsAssembler
 include("gsAssembler.jl")
 
@@ -17,9 +20,6 @@ include("gsHSplines.jl")
 # Load gsIO
 include("gsIO.jl")
 
-# Load gsMatrix
-include("gsMatrix.jl")
-
 # Load gsMatrix
 include("gsModeling.jl")
 
diff --git a/src/gsHSplines.jl b/src/gsHSplines.jl
index 070e64b..ab224b8 100644
--- a/src/gsHSplines.jl
+++ b/src/gsHSplines.jl
@@ -1,5 +1,8 @@
 export
-THBSplineBasis
+    THBSplineBasis,
+    HBSplineBasis,
+    THBSpline,
+    HBSpline
 ########################################################################
 # gsTHBSplineBasis
 ########################################################################
@@ -11,14 +14,11 @@ THBSplineBasis
 - `basis::Basis`: a basis object
 
 # Examples
-```jldoctest
-a = 1
-b = 2
-a + b
-
+```jldoctest output=(false)
+kv = KnotVector([0.,0.,0.,1.,1.,1.])
+b = BSplineBasis(kv)
+thb = THBSplineBasis(b)
 # output
-
-3
 ```
 """
 function THBSplineBasis(basis::Basis)::Basis
@@ -35,3 +35,100 @@ function THBSplineBasis(basis::Basis)::Basis
     end
     return Basis(b)
 end
+
+"""
+    HBSplineBasis(basis::Basis)
+
+# Arguments
+- `basis::Basis`: a basis object
+
+# Examples
+```jldoctest output=(false)
+kv = KnotVector([0.,0.,0.,1.,1.,1.])
+b = BSplineBasis(kv)
+thb = HBSplineBasis(b)
+# output
+```
+"""
+function HBSplineBasis(basis::Basis)::Basis
+    if (domainDim(basis)==1)
+        b = ccall((:gsHBSplineBasis1_create,libgismo),Ptr{gsCBasis},(Ptr{gsCBasis},),basis.ptr)
+    elseif (domainDim(basis)==2)
+        b = ccall((:gsHBSplineBasis2_create,libgismo),Ptr{gsCBasis},(Ptr{gsCBasis},),basis.ptr)
+    elseif (domainDim(basis)==3)
+        b = ccall((:gsHBSplineBasis3_create,libgismo),Ptr{gsCBasis},(Ptr{gsCBasis},),basis.ptr)
+    elseif (domainDim(basis)==4)
+        b = ccall((:gsHBSplineBasis4_create,libgismo),Ptr{gsCBasis},(Ptr{gsCBasis},),basis.ptr)
+    else
+        error("HBSplineBasis not implemented for this dimension")
+    end
+    return Basis(b)
+end
+
+"""
+    THBSpline(basis::Basis, coefs::Matrix{Cdouble})
+
+# Arguments
+- `basis::Basis`: a basis object
+- `coefs::Matrix{Cdouble}`: a matrix of coefficients
+
+# Examples
+```jldoctest output=(false)
+kv = KnotVector([0.,0.,0.,1.,1.,1.])
+b = BSplineBasis(kv)
+thb = THBSplineBasis(b)
+coefs = rand(3,2)
+g = THBSpline(thb,coefs)
+# output
+```
+"""
+function THBSpline(basis::Basis, coefs::Matrix{Cdouble})::Geometry
+    @assert Base.size(coefs,1) == Gismo.size(basis) "THBSpline: coefs must have the same number of rows as the number of degrees of freedom"
+    cc = EigenMatrix(Base.size(coefs,1),Base.size(coefs,2),pointer(coefs))
+    if (domainDim(basis)==1)
+        g = ccall((:gsTHBSpline1_create,libgismo),Ptr{gsCGeometry},(Ptr{gsCBasis},Ptr{EigenMatrix}),basis.ptr,cc.ptr)
+    elseif (domainDim(basis)==2)
+        g = ccall((:gsTHBSpline2_create,libgismo),Ptr{gsCGeometry},(Ptr{gsCBasis},Ptr{EigenMatrix}),basis.ptr,cc.ptr)
+    elseif (domainDim(basis)==3)
+        g = ccall((:gsTHBSpline3_create,libgismo),Ptr{gsCGeometry},(Ptr{gsCBasis},Ptr{EigenMatrix}),basis.ptr,cc.ptr)
+    elseif (domainDim(basis)==4)
+        g = ccall((:gsTHBSpline4_create,libgismo),Ptr{gsCGeometry},(Ptr{gsCBasis},Ptr{EigenMatrix}),basis.ptr,cc.ptr)
+    else
+        error("THBSpline not implemented for this dimension")
+    end
+    return Geometry(g)
+end
+
+"""
+    HBSpline(basis::Basis, coefs::Matrix{Cdouble})
+
+# Arguments
+- `basis::Basis`: a basis object
+- `coefs::Matrix{Cdouble}`: a matrix of coefficients
+
+# Examples
+```jldoctest output=(false)
+kv = KnotVector([0.,0.,0.,1.,1.,1.])
+b = BSplineBasis(kv)
+hb = HBSplineBasis(b)
+coefs = rand(3,2)
+g = HBSpline(hb,coefs)
+# output
+```
+"""
+function HBSpline(basis::Basis, coefs::Matrix{Cdouble})::Geometry
+    @assert Base.size(coefs,1) == Gismo.size(basis) "HBSpline: coefs must have the same number of rows as the number of degrees of freedom"
+    cc = EigenMatrix(Base.size(coefs,1),Base.size(coefs,2),pointer(coefs))
+    if (domainDim(basis)==1)
+        g = ccall((:gsHBSpline1_create,libgismo),Ptr{gsCGeometry},(Ptr{gsCBasis},Ptr{EigenMatrix}),basis.ptr,cc.ptr)
+    elseif (domainDim(basis)==2)
+        g = ccall((:gsHBSpline2_create,libgismo),Ptr{gsCGeometry},(Ptr{gsCBasis},Ptr{EigenMatrix}),basis.ptr,cc.ptr)
+    elseif (domainDim(basis)==3)
+        g = ccall((:gsHBSpline3_create,libgismo),Ptr{gsCGeometry},(Ptr{gsCBasis},Ptr{EigenMatrix}),basis.ptr,cc.ptr)
+    elseif (domainDim(basis)==4)
+        g = ccall((:gsHBSpline4_create,libgismo),Ptr{gsCGeometry},(Ptr{gsCBasis},Ptr{EigenMatrix}),basis.ptr,cc.ptr)
+    else
+        error("HBSpline not implemented for this dimension")
+    end
+    return Geometry(g)
+end
\ No newline at end of file
diff --git a/src/gsIO.jl b/src/gsIO.jl
index 0e8cabf..d5a275b 100644
--- a/src/gsIO.jl
+++ b/src/gsIO.jl
@@ -180,13 +180,8 @@ Get a string from the option list
 # Return
 - `string::Cstring`: the value
 
-# Examples
-#```jldoctest
-opt = OptionList(Dict("key1"=>"value1"))
-println(getString(opt,"key1"))
-# output
-value1
-```
+!!! warning
+    The returned value is a Cstring, and its conversion to a Julia string is not trivial.
 """
 function getString(opt::OptionList,key::String)::Cstring
     return ccall((:gsOptionList_getString,libgismo),Cstring,(Ptr{gsCOptionList},Cstring),opt.ptr,key)
@@ -208,6 +203,7 @@ opt = OptionList(Dict("key1"=>1))
 println(getInt(opt,"key1"))
 # output
 1
+```
 """
 function getInt(opt::OptionList,key::String)::Int
     return ccall((:gsOptionList_getInt,libgismo),Cint,(Ptr{gsCOptionList},Cstring),opt.ptr,key)
@@ -251,6 +247,7 @@ opt = OptionList(Dict("key1"=>true))
 println(getSwitch(opt,"key1"))
 # output
 true
+```
 """
 function getSwitch(opt::OptionList,key::String)::Bool
     return ccall((:gsOptionList_getSwitch,libgismo),Cint,(Ptr{gsCOptionList},Cstring),opt.ptr,key)