diff --git a/docs/src/functions.md b/docs/src/functions.md
index 84af1fd..d2d57fb 100644
--- a/docs/src/functions.md
+++ b/docs/src/functions.md
@@ -211,6 +211,66 @@ Comodo.geosphere
 Comodo.hexbox
 ```
 
+## `con_face_edge`
+```@docs
+Comodo.con_face_edge
+```
+
+## `con_edge_face`
+```@docs
+Comodo.con_edge_face
+```
+
+## `con_face_face`
+```@docs
+Comodo.con_face_face
+```
+
+## `con_face_face_v`
+```@docs
+Comodo.con_face_face_v
+```
+
+## `con_vertex_simplex`
+```@docs
+Comodo.con_vertex_simplex
+```
+
+## `con_vertex_face`
+```@docs
+Comodo.con_vertex_face
+```
+
+## `con_vertex_edge`
+```@docs
+Comodo.con_vertex_edge
+```
+
+## `con_edge_edge`
+```@docs
+Comodo.con_edge_edge
+```
+
+## `con_vertex_vertex_f`
+```@docs
+Comodo.con_vertex_vertex_f
+```
+
+## `con_vertex_vertex`
+```@docs
+Comodo.con_vertex_vertex
+```
+
+## `meshconnectivity`
+```@docs
+Comodo.meshconnectivity
+```
+
+## `mergevertices`
+```@docs
+Comodo.mergevertices
+```
+
 ## `smoothmesh_laplacian`
 ```@docs
 Comodo.smoothmesh_laplacian
@@ -221,6 +281,16 @@ Comodo.smoothmesh_laplacian
 Comodo.smoothmesh_hc
 ```
 
+## `quadplate`
+```@docs
+Comodo.quadplate
+```
+
+## `quadsphere`
+```@docs
+Comodo.quadsphere
+```
+
 ## `loflinear`
 ```@docs
 Comodo.loftlinear
diff --git a/examples/demo_dualclad.jl b/examples/demo_dualclad.jl
index 7457bf6..97ea394 100644
--- a/examples/demo_dualclad.jl
+++ b/examples/demo_dualclad.jl
@@ -3,12 +3,17 @@ using GLMakie
 using GeometryBasics
 using FileIO
 using StaticArrays
+using Makie.Colors
+
+c1 = RGBf(1.0, 0.30196078431372547, 0.023529411764705882)
+c2 = RGBf(0.2235294117647059, 1.0, 0.0784313725490196)
+c3 = RGBf(0.5803921568627451, 0.3411764705882353, 0.9215686274509803)
 
 #=
 This demo shows the use of the dualclag function. 
 =#
 
-testCase = 2
+testCase = 1
 if testCase == 1 
     F,V = geosphere(2,1.0)
 elseif testCase == 2 
@@ -61,17 +66,20 @@ end
 
 
 s = 0.5
-con_type = :edge
+con_type = :face
 Fs,Fq,Vs = dualclad(F,V,s; connectivity=con_type)
 
-# Rotate the coordinates
-fig = Figure(size = (800,800))
-ax = Axis3(fig[1, 1], aspect = :data)
+# Visualisation
+strokewidth =1
 
-# hp1 = poly!(ax, GeometryBasics.Mesh(V,F), color=:blue,transparency=true,strokewidth=0.15,strokecolor=:black, shading = FastShading)
+Vn = V-0.01*vertexnormal(F,V)
+
+fig = Figure(size = (1200,1200))
+ax = Axis3(fig[1, 1], aspect = :data)
 
-hp2 = poly!(ax, GeometryBasics.Mesh(Vs,Fs), color=:white,transparency=false,strokewidth=0.15,strokecolor=:black, shading = FastShading)
-hp3 = poly!(ax, GeometryBasics.Mesh(Vs,Fq), color=:white,transparency=false,strokewidth=0.15,strokecolor=:black, shading = FastShading)
+hp1 = mesh!(ax, GeometryBasics.Mesh(Vn,F), color=c3,transparency=false, shading = FastShading)
+hp2 = poly!(ax, GeometryBasics.Mesh(Vs,Fs), color=c1,transparency=false,strokewidth=strokewidth,strokecolor=:white, shading = FastShading)
+hp3 = poly!(ax, GeometryBasics.Mesh(Vs,Fq), color=c2,transparency=false,strokewidth=strokewidth,strokecolor=:white, shading = FastShading)
 
 # hp2 = mesh!(ax, GeometryBasics.Mesh(Vs,Fs), color=:white,transparency=false, shading = FastShading)
 # hp3 = mesh!(ax, GeometryBasics.Mesh(Vs,Fq), color=:white,transparency=false, shading = FastShading)
@@ -90,8 +98,8 @@ end
 
 set_close_to!(hSlider,0.5)
 
-# fileName = comododir()*"/assets/temp/dualclad_anim_05.mp4"
-# slider2anim(fig,hSlider,fileName; backforth=true, duration=3)
+fileName = comododir()*"/assets/temp/dualclad_anim_06.mp4"
+slider2anim(fig,hSlider,fileName; backforth=true, duration=6)
 
 fig
 
diff --git a/examples/demo_mergevertices.jl b/examples/demo_mergevertices.jl
index 15179e4..de28038 100644
--- a/examples/demo_mergevertices.jl
+++ b/examples/demo_mergevertices.jl
@@ -50,7 +50,7 @@ poly!(ax2,M2,strokewidth=2,color=:white, shading = FastShading)
 
 fig[2, :]=hSlider
 
-fileName = comododir()*"/assets/img/mergevertices_anim.gif"
-slider2anim(fig,hSlider,fileName; backforth=true, duration=2)
+# fileName = comododir()*"/assets/temp/mergevertices_anim.mp4"
+# slider2anim(fig,hSlider,fileName; backforth=true, duration=2)
 
 fig
diff --git a/examples/demo_regiontrimesh.jl b/examples/demo_regiontrimesh.jl
index d16c75e..aa50fd1 100644
--- a/examples/demo_regiontrimesh.jl
+++ b/examples/demo_regiontrimesh.jl
@@ -104,8 +104,8 @@ ax1 = Axis3(fig[1, 1],aspect = :data,title="Multi-region meshing",azimuth=-pi/2,
 # hp3 = scatter!(ax1, Vg,markersize=8,color=d,colormap=:Spectral)
 
 
-Fp,Vp = separate_vertices(F,V)
-Cp = simplex2vertexdata(Fp,C)
+Fp,Vp = separate_vertices(F,V) # Give each face its own point set 
+Cp = simplex2vertexdata(Fp,C) # Convert face color data to vertex color data 
 hp4 = poly!(ax1,GeometryBasics.Mesh(Vp,Fp), strokewidth=1,color=Cp, strokecolor=:black, shading = FastShading, transparency=false,colormap=Makie.Categorical(Makie.Reverse(:Spectral)))
 
 
diff --git a/examples/demo_scalesimplex.jl b/examples/demo_scalesimplex.jl
index 01bd619..b60f99b 100644
--- a/examples/demo_scalesimplex.jl
+++ b/examples/demo_scalesimplex.jl
@@ -21,27 +21,36 @@ F = tofaces(faces(M))
 V = topoints(coordinates(M))
 F,V = mergevertices(F,V)
 
-n = 1
-F,V = subtri(F,V,n; method=:Loop)
+# n = 1
+# F,V = subtri(F,V,n; method=:Loop)
 
 s = 0.5
 Fs,Vs = scalesimplex(F,V,s)
 
 
 ## Visualisation
-fig = Figure(size = (800,800))
+
+strokewidth = 0.1
+strokecolor = :white
+# Create slightly inward offset version so the visualisation does not co-inside
+Vn = V-0.01*vertexnormal(F,V)
+
+fig = Figure(size = (1400,800))
 
 ax1 = Axis3(fig[1, 1], aspect = :data,title="Single scale factor")
-mesh!(ax1, GeometryBasics.Mesh(V,F), color=:white,transparency=false,shading = FastShading)
-hp1 = poly!(ax1, GeometryBasics.Mesh(Vs,Fs), color=:green,transparency=false,strokewidth=1,strokecolor=:black,shading = FastShading)
+mesh!(ax1, GeometryBasics.Mesh(Vn,F), color=:white,transparency=false,shading = FastShading)
+hp1 = poly!(ax1, GeometryBasics.Mesh(Vs,Fs), color=:green,transparency=false,strokewidth=strokewidth,strokecolor=strokecolor,shading = FastShading)
+# hp1 = mesh!(ax1, GeometryBasics.Mesh(Vs,Fs), color=:green,transparency=false,shading = FastShading)
 
 ax2 = Axis3(fig[1, 2], aspect = :data,title="Per vertex scale factor")
-mesh!(ax2, GeometryBasics.Mesh(V,F), color=:white,transparency=false,shading = FastShading)
-hp2 = poly!(ax2, GeometryBasics.Mesh(Vs,Fs), color=:green,transparency=false,strokewidth=1,strokecolor=:black,shading = FastShading)
+mesh!(ax2, GeometryBasics.Mesh(Vn,F), color=:white,transparency=false,shading = FastShading)
+hp2 = poly!(ax2, GeometryBasics.Mesh(Vs,Fs), color=:green,transparency=false,strokewidth=strokewidth,strokecolor=strokecolor,shading = FastShading)
+# hp2 = mesh!(ax1, GeometryBasics.Mesh(Vs,Fs), color=:green,transparency=false,shading = FastShading)
 
-ax2 = Axis3(fig[1, 3], aspect = :data,title="Per face scale factor")
-mesh!(ax2, GeometryBasics.Mesh(V,F), color=:white,transparency=false,shading = FastShading)
-hp3 = poly!(ax2, GeometryBasics.Mesh(Vs,Fs), color=:green,transparency=false,strokewidth=1,strokecolor=:black,shading = FastShading)
+ax3 = Axis3(fig[1, 3], aspect = :data,title="Per face scale factor")
+mesh!(ax3, GeometryBasics.Mesh(Vn,F), color=:white,transparency=false,shading = FastShading)
+hp3 = poly!(ax3, GeometryBasics.Mesh(Vs,Fs), color=:green,transparency=false,strokewidth=strokewidth,strokecolor=strokecolor,shading = FastShading)
+# hp3 = mesh!(ax3, GeometryBasics.Mesh(Vs,Fs), color=:green,transparency=false,shading = FastShading)
 
 
 stepRange = range(1.0,0.0,50)
@@ -75,6 +84,8 @@ set_close_to!(hSlider,0.5)
 # fileName = comododir()*"/assets/temp/surface_mesh_smoothing_anim.mp4"
 # slider2anim(fig,hSlider,fileName; backforth=true, duration=2)
 
+
+
 fig
 
 
diff --git a/examples/demo_subquad.jl b/examples/demo_subquad.jl
index 1fadd8e..2d81d27 100644
--- a/examples/demo_subquad.jl
+++ b/examples/demo_subquad.jl
@@ -1,33 +1,90 @@
 using Comodo
 using GLMakie
 using GeometryBasics
+using LinearAlgebra
 
-## Define example input
-r = 0.5 #radius
-M = platonicsolid(2,r) # Get an example quadrilateral mesh (for a cube in this case)
-V = coordinates(M)
-F = faces(M)
+testCase = 2
 
+## Define example input
+if testCase == 1 # cube
+    r = 1.0 #radius
+    M = platonicsolid(2,r) # Get an example quadrilateral mesh (for a cube in this case)
+    V = coordinates(M)
+    F = faces(M)
+elseif testCase == 2 # Extruded prism/cylinder with nc points
+    r = 1.0
+    nc = 3
+    Vc = circlepoints(r,nc;dir=:cw)    
+    d = norm(Vc[1]-Vc[2])
+    n = normalizevector(Vec{3, Float64}(0.0,0.0,1.0))
+    s = 1    
+    F,V = extrudecurve(Vc,d;s=s, n=n, num_steps=2, close_loop=true,face_type=:quad)
+    M = GeometryBasics.Mesh(V,F)
+elseif testCase == 3 # Quadrangulated hemi-sphere
+    n = 1
+    r = 1.0
+    F,V = quadsphere(n,r)
+    VC = simplexcenter(F,V) # Finding triangle centre coordiantes
+    F = [F[i] for i in findall(map(v-> v[3]>0,VC))] # Remove some faces using z of central coordinates
+    F,V = remove_unused_vertices(F,V) # Cleanup/remove unused vertices after faces were removed
+    M = GeometryBasics.Mesh(V,F)
+end
 ## Refine mesh using `subquad` and the default "linear" method
 Fn1,Vn1=subquad(F,V,1) # Split once 
 Fn2,Vn2=subquad(F,V,2) # Split twice
 Fn3,Vn3=subquad(F,V,3) # Split 3 times
 
+Fn4,Vn4=subquad(F,V,1;method=:Catmull_Clark) # Split once 
+Fn5,Vn5=subquad(F,V,2;method=:Catmull_Clark) # Split twice
+Fn6,Vn6=subquad(F,V,3;method=:Catmull_Clark) # Split 3 times
+
+Fn7,Vn7=subquad(F,V,1; method=:Catmull_Clark, constrain_boundary=true) # Split once 
+Fn8,Vn8=subquad(F,V,2; method=:Catmull_Clark, constrain_boundary=true) # Split twice
+Fn9,Vn9=subquad(F,V,3; method=:Catmull_Clark, constrain_boundary=true) # Split 3 times
+
+
 ## Visualization
+strokewidth1 = 1
+lineWidth = 4
 fig = Figure(size=(1600,800))
 
-ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Refined n=1")
-wireframe!(ax1,M, linewidth=8,color=:red, overdraw=false)
-poly!(ax1,GeometryBasics.Mesh(Vn1,Fn1), strokewidth=3,color=:white,shading=FastShading,transparency=false)
-
-ax2 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Refined n=2")
-wireframe!(ax2,M, linewidth=8,color=:red, overdraw=false)
-poly!(ax2,GeometryBasics.Mesh(Vn2,Fn2), strokewidth=3,color=:white,shading=FastShading,transparency=false)
+ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Linear, n=1")
+wireframe!(ax1,M, linewidth=lineWidth,color=:red, overdraw=false)
+poly!(ax1,GeometryBasics.Mesh(Vn1,Fn1), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false)
 
-ax3 = Axis3(fig[1, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Refined n=3")
-hp1 = wireframe!(ax3,M, linewidth=8,color=:red, overdraw=false)
-hp2 = poly!(ax3,GeometryBasics.Mesh(Vn3,Fn3), strokewidth=3,color=:white,shading=FastShading,transparency=false)
+ax2 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Linearly, n=2")
+wireframe!(ax2,M, linewidth=lineWidth,color=:red, overdraw=false)
+poly!(ax2,GeometryBasics.Mesh(Vn2,Fn2), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false)
 
+ax3 = Axis3(fig[1, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Linearly, n=3")
+hp1 = wireframe!(ax3,M, linewidth=lineWidth,color=:red, overdraw=false)
+hp2 = poly!(ax3,GeometryBasics.Mesh(Vn3,Fn3), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false)
 Legend(fig[1, 4],[hp1,hp2],["Initial","Refined"])
 
+ax4 = Axis3(fig[2, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Catmull_Clark, n=1")
+wireframe!(ax4,M, linewidth=lineWidth,color=:red, overdraw=false)
+poly!(ax4,GeometryBasics.Mesh(Vn4,Fn4), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false)
+
+ax5 = Axis3(fig[2, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Catmull_Clark, n=2")
+wireframe!(ax5,M, linewidth=lineWidth,color=:red, overdraw=false)
+poly!(ax5,GeometryBasics.Mesh(Vn5,Fn5), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false)
+
+ax6 = Axis3(fig[2, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Catmull_Clark, n=3")
+hp1 = wireframe!(ax6,M, linewidth=lineWidth,color=:red, overdraw=false)
+hp2 = poly!(ax6,GeometryBasics.Mesh(Vn6,Fn6), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false)
+Legend(fig[2, 4],[hp1,hp2],["Initial","Refined"])
+
+ax4 = Axis3(fig[3, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Catmull_Clark, constrained boundary, n=1")
+wireframe!(ax4,M, linewidth=lineWidth,color=:red, overdraw=false)
+poly!(ax4,GeometryBasics.Mesh(Vn7,Fn7), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false)
+
+ax5 = Axis3(fig[3, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Catmull_Clark, constrained boundary, n=2")
+wireframe!(ax5,M, linewidth=lineWidth,color=:red, overdraw=false)
+poly!(ax5,GeometryBasics.Mesh(Vn8,Fn8), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false)
+
+ax6 = Axis3(fig[3, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Catmull_Clark, constrained boundary, n=3")
+hp1 = wireframe!(ax6,M, linewidth=lineWidth,color=:red, overdraw=false)
+hp2 = poly!(ax6,GeometryBasics.Mesh(Vn9,Fn9), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false)
+Legend(fig[3, 4],[hp1,hp2],["Initial","Refined"])
+
 fig
\ No newline at end of file
diff --git a/examples/demo_subtri.jl b/examples/demo_subtri.jl
index 7d2c969..393fc28 100644
--- a/examples/demo_subtri.jl
+++ b/examples/demo_subtri.jl
@@ -20,10 +20,21 @@ one, and 3 at each corner). The following refinement methods are implemented:
 =#
 
 ## Define example input
-r = 0.5 #radius
-M = platonicsolid(4,r) # Get example triangulated mesh (e.g. icosahedron)
-V = coordinates(M)
-F = faces(M)
+testCase = 2
+if testCase == 1 # icosahedron
+    r = 0.5 #radius
+    M = platonicsolid(4,r) 
+    V = coordinates(M)
+    F = faces(M)
+elseif testCase == 2 # Extruded prism/cylinder with nc points
+    r = 1.0
+    nc = 3
+    Vc = circlepoints(r,nc;dir=:cw)    
+    d = norm(Vc[1]-Vc[2])    
+    s = 1    
+    F,V = extrudecurve(Vc,d;s=s, num_steps=2, close_loop=true,face_type=:tri_slash)
+    M = GeometryBasics.Mesh(V,F)
+end
 
 ## Refine triangulation using `subtri` and the default :linear method
 
@@ -35,33 +46,53 @@ Fn4,Vn4=subtri(F,V,1; method=:Loop) # Split once
 Fn5,Vn5=subtri(F,V,2; method=:Loop) # Split twice
 Fn6,Vn6=subtri(F,V,3; method=:Loop) # Split 3 times
 
+Fn7,Vn7=subtri(F,V,1; method=:Loop, constrain_boundary=true) # Split once 
+Fn8,Vn8=subtri(F,V,2; method=:Loop, constrain_boundary=true) # Split twice
+Fn9,Vn9=subtri(F,V,3; method=:Loop, constrain_boundary=true) # Split 3 times
+
+
 ## Visualization
+strokewidth1 = 1
+lineWidth = 4
 fig = Figure(size=(1600,800))
 
-ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Linearly refined n=1")
-wireframe!(ax1,M, linewidth=8,color=:red, overdraw=false)
-poly!(ax1,GeometryBasics.Mesh(Vn1,Fn1), strokewidth=3,color=:white,shading=FastShading,transparency=false)
+ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Linear, n=1")
+wireframe!(ax1,M, linewidth=lineWidth,color=:red, overdraw=false)
+poly!(ax1,GeometryBasics.Mesh(Vn1,Fn1), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false)
 
-ax2 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Linearly refined n=2")
-wireframe!(ax2,M, linewidth=8,color=:red, overdraw=false)
-poly!(ax2,GeometryBasics.Mesh(Vn2,Fn2), strokewidth=3,color=:white,shading=FastShading,transparency=false)
+ax2 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Linear, n=2")
+wireframe!(ax2,M, linewidth=lineWidth,color=:red, overdraw=false)
+poly!(ax2,GeometryBasics.Mesh(Vn2,Fn2), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false)
 
-ax3 = Axis3(fig[1, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Linearly refined n=3")
-hp1 = wireframe!(ax3,M, linewidth=8,color=:red, overdraw=false)
-hp2 = poly!(ax3,GeometryBasics.Mesh(Vn3,Fn3), strokewidth=3,color=:white,shading=FastShading,transparency=false)
+ax3 = Axis3(fig[1, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Linear, n=3")
+hp1 = wireframe!(ax3,M, linewidth=lineWidth,color=:red, overdraw=false)
+hp2 = poly!(ax3,GeometryBasics.Mesh(Vn3,Fn3), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false)
 Legend(fig[1, 4],[hp1,hp2],["Initial","Refined"])
 
-ax4 = Axis3(fig[2, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Loop-method refined n=1")
-wireframe!(ax4,M, linewidth=8,color=:red, overdraw=false)
-poly!(ax4,GeometryBasics.Mesh(Vn4,Fn4), strokewidth=3,color=:white,shading=FastShading,transparency=false)
+ax4 = Axis3(fig[2, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Loop, n=1")
+wireframe!(ax4,M, linewidth=lineWidth,color=:red, overdraw=false)
+poly!(ax4,GeometryBasics.Mesh(Vn4,Fn4), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false)
 
-ax5 = Axis3(fig[2, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Loop-method refined n=2")
-wireframe!(ax5,M, linewidth=8,color=:red, overdraw=false)
-poly!(ax5,GeometryBasics.Mesh(Vn5,Fn5), strokewidth=3,color=:white,shading=FastShading,transparency=false)
+ax5 = Axis3(fig[2, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Loop, n=2")
+wireframe!(ax5,M, linewidth=lineWidth,color=:red, overdraw=false)
+poly!(ax5,GeometryBasics.Mesh(Vn5,Fn5), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false)
 
-ax6 = Axis3(fig[2, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Loop-method refined n=3")
-hp1 = wireframe!(ax6,M, linewidth=8,color=:red, overdraw=false)
-hp2 = poly!(ax6,GeometryBasics.Mesh(Vn6,Fn6), strokewidth=3,color=:white,shading=FastShading,transparency=false)
+ax6 = Axis3(fig[2, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Loop, n=3")
+hp1 = wireframe!(ax6,M, linewidth=lineWidth,color=:red, overdraw=false)
+hp2 = poly!(ax6,GeometryBasics.Mesh(Vn6,Fn6), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false)
 Legend(fig[2, 4],[hp1,hp2],["Initial","Refined"])
 
+ax4 = Axis3(fig[3, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Loop, constrained boundary, n=1")
+wireframe!(ax4,M, linewidth=lineWidth,color=:red, overdraw=false)
+poly!(ax4,GeometryBasics.Mesh(Vn7,Fn7), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false)
+
+ax5 = Axis3(fig[3, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Loop, constrained boundary, n=2")
+wireframe!(ax5,M, linewidth=lineWidth,color=:red, overdraw=false)
+poly!(ax5,GeometryBasics.Mesh(Vn8,Fn8), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false)
+
+ax6 = Axis3(fig[3, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Loop, constrained boundary, n=3")
+hp1 = wireframe!(ax6,M, linewidth=lineWidth,color=:red, overdraw=false)
+hp2 = poly!(ax6,GeometryBasics.Mesh(Vn9,Fn9), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false)
+Legend(fig[3, 4],[hp1,hp2],["Initial","Refined"])
+
 fig
diff --git a/src/functions.jl b/src/functions.jl
index f587197..849fdf5 100644
--- a/src/functions.jl
+++ b/src/functions.jl
@@ -1575,14 +1575,25 @@ method both refines and smoothes the geometry through spline approximation.
 1. [Charles Loop, _Smooth Subdivision Surfaces Based on Triangles_, M.S. Mathematics Thesis, University of Utah. 1987.](https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/thesis-10.pdf)
 2. [Jos Stam, Charles Loop, _Quad/Triangle Subdivision_, doi: 10.1111/1467-8659.t01-2-00647](https://doi.org/10.1111/1467-8659.t01-2-00647)
 """
-function subtri(F::Vector{NgonFace{3,TF}},V::Vector{Point{ND,TV}},n::Int64; method = :linear) where TF<:Integer where ND where TV <: Real
+function subtri(F::Vector{NgonFace{3,TF}},V::Vector{Point{ND,TV}},n::Int64; method = :linear, constrain_boundary=false) where TF<:Integer where ND where TV <: Real
     
     if iszero(n)
         return F,V
     elseif isone(n)
         E = meshedges(F)
         Eu,_,indReverse = gunique(E; return_unique=true, return_index=true, return_inverse=true, sort_entries=true)
-        
+        # Check for boundary edges        
+        count_E2F = count_edge_face(F,Eu,indReverse)
+        B_boundary = isone.(count_E2F)
+        if any(B_boundary)
+            treatBoundary = true
+            Eb = Eu[B_boundary]
+            indB = unique(reduce(vcat,Eb))
+            con_V2V = con_vertex_vertex(Eb,V)
+        else
+            treatBoundary = false
+        end
+
         Fm1 = [TriangleFace{TF}(a.+length(V)) for a in eachrow(reshape(indReverse,length(F),length(F[1])))] 
         Fm2 = Vector{TriangleFace{TF}}(undef,length(Fm1))
         Fm3 = Vector{TriangleFace{TF}}(undef,length(Fm1))
@@ -1603,39 +1614,48 @@ function subtri(F::Vector{NgonFace{3,TF}},V::Vector{Point{ND,TV}},n::Int64; meth
             # Create complete point set
             Vn = [V; simplexcenter(Eu,V)]  # Old and new mid-edge points          
         elseif method == :Loop #Loop subdivision 
-    
             # New mid-edge like vertices
             Vm = Vector{Point{ND,TV}}(undef,length(Eu)) 
-            for q in eachindex(Eu) # For each edge index                        
-                F_touch = F[con_E2F[q]] # Faces sharing current edge, mostly 2 but 1 for a boundary edge
-                indVerticesTouch = Vector{TF}() 
-                for f in F_touch        
-                    b = f.!=Eu[q][1] .&& f.!=Eu[q][2]      
-                    if any(b)  
-                        append!(indVerticesTouch,f[b])           
-                    end
-                end        
-                Vm[q]=3/8 .*(V[Eu[q][1]] .+ V[Eu[q][2]])  .+ 1/8 .* (V[indVerticesTouch[1]] .+ V[indVerticesTouch[2]])
+            for q in eachindex(Eu) # For each edge index       
+                if treatBoundary && B_boundary[q]              
+                    Vm[q] = 1/2 .*(V[Eu[q][1]] .+ V[Eu[q][2]]) # Normal mid-edge point
+                else
+                    F_touch = F[con_E2F[q]] 
+                    indVerticesTouch = Vector{TF}() 
+                    for f in F_touch        
+                        b = f.!=Eu[q][1] .&& f.!=Eu[q][2]      
+                        if any(b)  
+                            append!(indVerticesTouch,f[b])           
+                        end
+                    end                   
+                    Vm[q] = 3/8 .*(V[Eu[q][1]] .+ V[Eu[q][2]])  .+ 1/8 .* (V[indVerticesTouch[1]] .+ V[indVerticesTouch[2]])
+                end
             end
     
             # Modified vertices for original vertices
-            Vv = Vector{Point{ND,TV}}(undef,length(V))
+            Vv = deepcopy(V)
             for q in eachindex(V)            
-                B_vert_face = [any(f.==q) for f in F]
-                F_touch = F[B_vert_face] # Faces mostly 2 but 1 for a boundary edge
-                indVerticesTouch = Vector{TF}()
-                for f in F_touch                
-                    indTouch = f[f.!=q]        
-                    for i in indTouch 
-                        if i ∉ indVerticesTouch 
-                            push!(indVerticesTouch,i)
+                if treatBoundary && in(q,indB)
+                    if !constrain_boundary
+                        Vv[q] = 6/8*V[q] + 1/8*(V[con_V2V[q][1]]+V[con_V2V[q][2]]) 
+                    end
+                else
+                    B_vert_face = [any(f.==q) for f in F]
+                    F_touch = F[B_vert_face] 
+                    indVerticesTouch = Vector{TF}()
+                    for f in F_touch                
+                        indTouch = f[f.!=q]        
+                        for i in indTouch 
+                            if i ∉ indVerticesTouch 
+                                push!(indVerticesTouch,i)
+                            end
                         end
                     end
+                    N = length(indVerticesTouch)                
+                    v_sum = sum(V[indVerticesTouch],dims=1)[1]                
+                    β = 1/N * (5/8-(3/8 +1/4*cos((2*π)/N))^2)        
+                    Vv[q] = (1-N*β) .* V[q] .+ β*v_sum                
                 end
-                N = length(indVerticesTouch)                
-                v_sum = sum(V[indVerticesTouch],dims=1)[1]                
-                β = 1/N * (5/8-(3/8 +1/4*cos((2*π)/N))^2)        
-                Vv[q] = (1-N*β) .* V[q] .+ β*v_sum                
             end    
             # Create complete point set
             Vn = [Vv;Vm] # Updated orignals and new "mid-edge-ish" points
@@ -1646,8 +1666,8 @@ function subtri(F::Vector{NgonFace{3,TF}},V::Vector{Point{ND,TV}},n::Int64; meth
         return Fn,Vn    
 
     elseif n>1
-        for _ =1:n
-            F,V = subtri(F,V,1; method=method)
+        for _ = 1:n
+            F,V = subtri(F,V,1; method=method, constrain_boundary=constrain_boundary)
         end
         return F,V
     else
@@ -1681,16 +1701,26 @@ spline approximation.
 # References
 1. [E. Catmull and J. Clark, _Recursively generated B-spline surfaces on arbitrary topological meshes_, Computer-Aided Design, vol. 10, no. 6, pp. 350-355, Nov. 1978, doi: 10.1016/0010-4485(78)90110-0](https://doi.org/10.1016/0010-4485(78)90110-0).
 """
-function subquad(F::Vector{NgonFace{4,TF}},V::Vector{Point{ND,TV}},n::Int64; method=:linear) where TF<:Integer where ND where TV <: Real
+function subquad(F::Vector{NgonFace{4,TF}},V::Vector{Point{ND,TV}},n::Int64; method=:linear, constrain_boundary=false) where TF<:Integer where ND where TV <: Real
     if iszero(n)
         return F,V
     elseif isone(n)
 
         # Get edges
         E = meshedges(F) # Non-unique edges
-
         Eu,_,indReverse = gunique(E; return_unique=true, return_index=true, return_inverse=true, sort_entries=true)
-        
+        # Check for boundary edges        
+        count_E2F = count_edge_face(F,Eu,indReverse)
+        B_boundary = isone.(count_E2F)
+        if any(B_boundary)
+            treatBoundary = true
+            Eb = Eu[B_boundary]
+            indB = unique(reduce(vcat,Eb))
+            con_V2V = con_vertex_vertex(Eb,V)
+        else
+            treatBoundary = false
+        end
+
         con_F2E = con_face_edge(F,Eu,indReverse)
         con_E2F = con_edge_face(F,Eu,indReverse)
         con_V2E = con_vertex_edge(Eu,V)
@@ -1703,22 +1733,32 @@ function subquad(F::Vector{NgonFace{4,TF}},V::Vector{Point{ND,TV}},n::Int64; met
             Vn = [V;Ve;Vf] # Joined point set
         elseif method ==:Catmull_Clark
             # Mid face points
-            Vf = simplexcenter(F,V)  
+            Vf = simplexcenter(F,V) # Mid face points
             Ve_mid = simplexcenter(Eu,V) # Mid edge points
 
             # Edge points 
             Ve = Vector{Point{ND,TV}}(undef,length(Eu)) # Initialize edge points
-            for q in eachindex(Eu)                         
-                Ve[q] = (mean(Vf[con_E2F[q]],dims=1)[1] .+ Ve_mid[q])./2.0
+            for q in eachindex(Eu)      
+                if treatBoundary && B_boundary[q]              
+                    Ve[q] = Ve_mid[q] # Normal mid-edge point
+                else                   
+                    Ve[q] = (mean(Vf[con_E2F[q]],dims=1)[1] .+ Ve_mid[q])./2.0
+                end
             end
 
             # Vertex points 
-            Vv = Vector{Point{ND,TV}}(undef,length(V)) # Initialize vertex points
+            Vv = deepcopy(V) # Initialize vertex points
             for q in eachindex(V) # Loop over all vertices
-                indF = con_V2F[q]
-                indE = con_V2E[q]
-                N = length(indF) # Number of faces (or edges) touching this vertex                    
-                Vv[q] = (mean(Vf[indF],dims=1)[1] .+ 2.0.*mean(Ve_mid[indE],dims=1)[1] .+ (N-3.0).*V[q])./N
+                if treatBoundary && in(q,indB)
+                    if !constrain_boundary
+                        Vv[q] = 6/8*V[q] + 1/8*(V[con_V2V[q][1]]+V[con_V2V[q][2]]) 
+                    end
+                else
+                    indF = con_V2F[q]
+                    indE = con_V2E[q]
+                    N = length(indF) # Number of faces (or edges) touching this vertex                    
+                    Vv[q] = (mean(Vf[indF],dims=1)[1] .+ 2.0.*mean(Ve_mid[indE],dims=1)[1] .+ (N-3.0).*V[q])./N
+                end
             end
 
             Vn = [Vv;Ve;Vf] # Joined point set
@@ -1739,7 +1779,7 @@ function subquad(F::Vector{NgonFace{4,TF}},V::Vector{Point{ND,TV}},n::Int64; met
         return Fn,Vn
     elseif n>1
         for _ =1:n
-            F,V = subquad(F,V,1;method=method)
+            F,V = subquad(F,V,1;method=method, constrain_boundary=constrain_boundary)
         end
         return F,V
     else
@@ -1859,7 +1899,20 @@ function hexbox(boxDim,boxEl)
     return E,V,F,Fb,CFb_type
 end
 
+"""
+    con_face_edge(F,E_uni=nothing,indReverse=nothing)
+
+Returns the edges connected to each face.
+
+# Description
 
+This function computes the face-edge connectivity. The input faces `F` (and 
+optionally also the unique edges `E_uni` and reverse indices `indReverse` to map
+to the non-unique edges, see also `gunique`) are used to create a list of edges 
+connected to each face. If `F` contains N faces then the output contains N such 
+lists. For triangles the output contains 3 edges per faces, for quads 4 per face
+and so on.   
+"""
 function con_face_edge(F,E_uni=nothing,indReverse=nothing)
     if isnothing(E_uni) | isnothing(indReverse)
         E = meshedges(F)
@@ -1869,6 +1922,20 @@ function con_face_edge(F,E_uni=nothing,indReverse=nothing)
 end
 
 
+"""
+con_edge_face(F,E_uni=nothing,indReverse=nothing)
+
+Returns the faces connected to each edge.
+
+# Description
+
+This function computes the edge-face connectivity. The input faces `F` (and 
+optionally also the unique edges `E_uni` and reverse indices `indReverse` to map
+to the non-unique edges, see also `gunique`) are used to create a list of faces 
+connected to each edges. If `E_uni` contains N edges then the output contains 
+N such lists. For non-boundary edges each edge should connect to 2 faces. 
+Boundary edges connect to just 1 face.  
+"""
 function con_edge_face(F,E_uni=nothing,indReverse=nothing)
     if isnothing(E_uni) || isnothing(indReverse)
         E = meshedges(F)
@@ -1886,6 +1953,23 @@ function con_edge_face(F,E_uni=nothing,indReverse=nothing)
 end
 
 
+"""
+    con_face_face(F,E_uni=nothing,indReverse=nothing,con_E2F=nothing,con_F2E=nothing)
+
+Returns the edge-connected faces for each face.
+
+# Description
+
+This function computes the face-face connectivity for each face. The input faces
+`F` are used to create a list of faces connected to each face by a shared edge.
+For non-boundary triangles for instance the output contains 3 edges per faces 
+(which may be less for boundary triangles), and similarly non-boundary quads 
+would each have 4 edge-connected faces. Additional optional inputs include: the 
+unique edges `E_uni`, the reverse indices `indReverse` to map to the non-unique 
+edges (see also `gunique`), as well as the edge-face `con_E2F` and face-edge 
+`con_F2E` connectivity. These are all needed for computing the face-face 
+connectivity and supplying them if already computed therefore save time.  
+"""
 function con_face_face(F,E_uni=nothing,indReverse=nothing,con_E2F=nothing,con_F2E=nothing)
     if length(F)>1 # More than one face so compute connectivity
         if isnothing(E_uni)| isnothing(indReverse)
@@ -1912,7 +1996,22 @@ function con_face_face(F,E_uni=nothing,indReverse=nothing,con_E2F=nothing,con_F2
     end
 end
 
+"""
+    con_face_face_v(F,V=nothing,con_V2F=nothing)
+
+Returns the vertex-connected faces for each face.
 
+# Description
+
+This function computes the face-face connectivity for each face. The input faces
+`F` are used to create a list of faces connected to each face by a shared vertex.
+Additional optional inputs include: the vertices `V`, and the vertex-face 
+connectivity `con_V2F`. In terms of vertices only the number of vertices, i.e. 
+`length(V)` is neede, if `V` is not provided it is assumed that `length(V)` 
+corresponds to the largest index in `F`. The vertex-face connectivity if not
+supplied, will be computed by this function, hence computational time may be saved 
+if it was already computed. 
+"""
 function con_face_face_v(F,V=nothing,con_V2F=nothing)
     if length(F)>1 # More than one face so compute connectivity
         if isnothing(con_V2F) 
@@ -1932,7 +2031,20 @@ function con_face_face_v(F,V=nothing,con_V2F=nothing)
     end
 end
 
+"""
+    con_vertex_simplex(F,V=nothing)
+
+Returns how vertices connect to simplices
 
+# Description
+
+This function computes the vertex-simplex connectivity for each vertex. The input
+simplices `F` are used to create a list of simplices connected to each vertex.
+Additional optional inputs include: the vertices `V`.
+In terms of vertices only the number of vertices, i.e. `length(V)` is needed, 
+if `V` is not provided it is assumed that `length(V)` corresponds to the largest
+index in `F`.  
+"""
 function con_vertex_simplex(F,V=nothing)
     if isnothing(V)
         n = maximum(reduce(vcat,F))
@@ -1948,17 +2060,53 @@ function con_vertex_simplex(F,V=nothing)
     return con_V2F
 end
 
+"""
+    con_vertex_face(F,V=nothing)
+
+Returns how vertices connect to faces
 
+# Description
+This function is an alias of `con_vertex_simplex`, and computes the vertex-face 
+connectivity for each vertex. The input faces `F` are used to create a list of 
+faces connected to each vertex. Additional optional inputs include: the vertices `V`.
+In terms of vertices only the number of vertices, i.e. `length(V)` is needed, 
+if `V` is not provided it is assumed that `length(V)` corresponds to the largest
+index in `F`.  
+"""
 function con_vertex_face(F,V=nothing)
     return con_vertex_simplex(F,V)
 end
 
 
+"""
+    con_vertex_edge(F,V=nothing)
+
+Returns how vertices connect to edges
+
+# Description
+This function is an alias of `con_vertex_simplex`, and computes the vertex-edge 
+connectivity for each vertex. The input edges `E` are used to create a list of 
+edges connected to each vertex. Additional optional inputs include: the vertices `V`.
+In terms of vertices only the number of vertices, i.e. `length(V)` is needed, 
+if `V` is not provided it is assumed that `length(V)` corresponds to the largest
+index in `E`.  
+"""
 function con_vertex_edge(E,V=nothing)
     return con_vertex_simplex(E,V)
 end
 
+"""
+    con_edge_edge(E_uni,con_V2E=nothing)
+
+Returns the vertex-connected edges for each edge.
+
+# Description
 
+This function computes the edge-edge connectivity for each edge. The input edges
+`F` are used to create a list of edges connected to each edge by a shared vertex.
+Additional optional inputs include: `con_V2E` (the vertex-edge connectivity), which
+is instead computed when not provided.  
+"""
 function con_edge_edge(E_uni,con_V2E=nothing)
     if isnothing(con_V2E)
         con_V2E = con_vertex_edge(E_uni) 
@@ -1974,7 +2122,21 @@ function con_edge_edge(E_uni,con_V2E=nothing)
     return con_E2E
 end
 
+"""
+    con_vertex_vertex_f(F,V=nothing,con_V2F=nothing)
+
+Returns the face-connected vertices for each vertex.
 
+# Description
+
+This function computes the vertex-vertex connectivity for each vertex using the 
+vertex connected faces. The input faces `F` are used to create a list of vertices
+connected to each vertex by a shared face. Additional optional inputs include: 
+the vertices `V` and `con_V2F` (the vertex-face connectivity). In terms of vertices
+only the number of vertices, i.e. `length(V)` is needed, if `V` is not provided 
+it is assumed that `length(V)` corresponds to the largest index in `F`. The 
+vertex-face connectivity `con_V2F` is needed, hence is computed when not provided.  
+"""
 function con_vertex_vertex_f(F,V=nothing,con_V2F=nothing)
     if isnothing(V)
         n = maximum(reduce(vcat,F))
@@ -2000,7 +2162,21 @@ function con_vertex_vertex_f(F,V=nothing,con_V2F=nothing)
 
 end
 
+"""
+    con_vertex_vertex(E,V=nothing,con_V2E=nothing)
+
+Returns the edge-connected vertices for each vertex.
 
+# Description
+
+This function computes the vertex-vertex connectivity for each vertex using the 
+vertex connected edges. The input edges `E` are used to create a list of vertices
+connected to each vertex by a shared edge. Additional optional inputs include: 
+the vertices `V` and `con_V2E` (the vertex-edge connectivity). In terms of vertices
+only the number of vertices, i.e. `length(V)` is needed, if `V` is not provided 
+it is assumed that `length(V)` corresponds to the largest index in `E`. The 
+vertex-edge connectivity `con_V2E` is needed, hence is computed when not provided.  
+"""
 function con_vertex_vertex(E,V=nothing,con_V2E=nothing)
     if isnothing(V)
         n = maximum(reduce(vcat,E))
@@ -2025,6 +2201,25 @@ function con_vertex_vertex(E,V=nothing,con_V2E=nothing)
     return con_V2V
 end
 
+"""
+    meshconnectivity(F::Vector{NgonFace{N,TF}},V::Vector{Point{ND,TV}}) where N where TF<:Integer where ND where TV<:Real
+
+Returns all mesh connectivity data
+
+# Description
+This function returns the `ConnectivitySet`, i.e. all mesh connectivity data for 
+the input mesh defined by the faces `F` and the vertices `V`. 
+The `ConnectivitySet` contains the following connectivity descriptions: 
+* face-edge
+* edge-face
+* face-face
+* face-face (wrt vertices)
+* vertex-face
+* vertex-edge
+* edge-edge
+* vertex-vertex
+* vertex-vertex (wrt faces)
+"""
 function meshconnectivity(F::Vector{NgonFace{N,TF}},V::Vector{Point{ND,TV}}) where N where TF<:Integer where ND where TV<:Real
 
     # EDGE-VERTEX connectivity
@@ -2061,6 +2256,16 @@ function meshconnectivity(F::Vector{NgonFace{N,TF}},V::Vector{Point{ND,TV}}) whe
     return ConnectivitySet(E_uni, con_E2F, con_E2E, F,  con_F2E, con_F2F, con_V2E, con_V2F, con_V2V, con_V2V_f, con_F2F_v) 
 end 
 
+"""
+    mergevertices(F::Vector{NgonFace{N,TF}},V::Vector{Point{ND,TV}}; roundVertices = true, numDigitsMerge=nothing) where N where TF<:Integer where ND where TV<:Real
+
+Merges points that coincide
+
+# Description
+This function take the faces `F` and vertices `V` and merges points that are sufficiently 
+similar. Once points are merged the indices in `F` are corrected for the new reduced
+point set. 
+"""
 function mergevertices(F::Vector{NgonFace{N,TF}},V::Vector{Point{ND,TV}}; roundVertices = true, numDigitsMerge=nothing) where N where TF<:Integer where ND where TV<:Real
 
     m = length(V)
@@ -2212,6 +2417,18 @@ function smoothmesh_hc(F::Vector{NgonFace{N,TF}},V::Vector{Point{ND,TV}}, n=1, 
     end
 end
 
+
+"""
+    quadplate(plateDim,plateElem)
+
+Returns a quad mesh for a plate
+
+# Description
+This function creates a quadrilateral mesh (faces `F` and vertices `V`) for a 
+plate. The dimensions in the x-, and y-direction are specified in the input vector 
+`plateDim`, and the number of elements to use in each direction in the input 
+vector `plateElem`. 
+"""
 function quadplate(plateDim,plateElem)
     num_x = plateElem[1]+1
     num_y = plateElem[2]+1
@@ -2235,6 +2452,18 @@ function quadplate(plateDim,plateElem)
     return F, V
 end
 
+"""
+    quadsphere(n::Int64,r::T) where T <: Real
+
+Returns a quadrangulated sphere
+
+# Description
+This function creates a quadrilateral mesh (faces `F` and vertices `V`) for a sphere
+with a radius defined by the input `r`. The input `n` defines the density of sphere
+mesh. The quad mesh is constructed using `subquad` subdivision of a regular cube, 
+whereby `n` sets the number of splitting iterations to use. Using `n=0` therefore
+returns a cube.
+"""
 function quadsphere(n::Int64,r::T) where T <: Real
     M = platonicsolid(2,r)
     F = faces(M)
diff --git a/test/runtests.jl b/test/runtests.jl
index 4441f91..82ff2cc 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1397,40 +1397,72 @@ end
 
 @testset "subtri" verbose = true begin
     eps_level = 1e-4
-    M = platonicsolid(4, 1.0)
+    M = platonicsolid(4, 1.0) # icosahedron with radius 1.0 
     V = coordinates(M)
     F = faces(M)
-    n = 3
+    n = 2
 
     @testset "Errors" begin
-        @test_throws ArgumentError subtri(F, V, n; method=:wrong)
+        @test_throws ArgumentError subtri(F, V, n; method=:wrong)        
         @test_throws ArgumentError subtri(F, V, -1)
     end
 
     @testset "linear" begin
-        Fn, Vn = subtri(F, V, n; method=:linear)
-
-        @test Fn isa Vector{TriangleFace{Int64}}
-        @test length(Fn) == 1280
-        @test Fn[1] == TriangleFace(163, 323, 243)
-        @test Vn isa Vector{Point3{Float64}}
-        @test length(Vn) == 642
-        @test isapprox(Vn[1], [0.0, -0.5257311121191336, -0.85065080835204], atol=eps_level)
-
         Fn, Vn = subtri(F, V, 0)
         @test Fn == F
         @test Vn == V
+
+        Fn, Vn = subtri(F, V, n; method=:linear)
+        ind = round.(Int64,range(1,length(Vn),6)) # Sample indices
+
+        @test eltype(F) == eltype(Fn)
+        @test length(Fn) == 4^n*length(F)        
+        @test eltype(V) == eltype(Vn)        
+        @test isapprox(Vn[ind], Point{3, Float64}[[0.0, -0.5257311121191336, -0.85065080835204], 
+        [-0.2628655560595668, 0.6881909602355868, 0.42532540417602], [0.0, -0.6881909602355868, -0.42532540417602], 
+        [0.1314327780297834, -0.7694208842938134, 0.21266270208801], [-0.21266270208801, 0.3942983340893502, 0.7694208842938134], 
+        [-0.63798810626403, 0.1314327780297834, -0.6069610361773602]], atol=eps_level)        
     end
 
     @testset "Loop" begin
+        Fn, Vn = subtri(F, V, 0; method=:Loop)
+        @test Fn == F
+        @test Vn == V
+
         Fn, Vn = subtri(F, V, n; method=:Loop)
+        ind = round.(Int64,range(1,length(Vn),6)) # Sample indices
 
-        @test Fn isa Vector{TriangleFace{Int64}}
-        @test length(Fn) == 1280
-        @test Fn[1] == TriangleFace(163, 323, 243)
-        @test Vn isa Vector{Point3{Float64}}
-        @test length(Vn) == 642
-        @test isapprox(Vn[1], [0.0, -0.37343167032888536, -0.6042251350677821], atol=eps_level)
+        @test eltype(F) == eltype(Fn)
+        @test length(Fn) == 4^n*length(F)        
+        @test eltype(V) == eltype(Vn)        
+        @test isapprox(Vn[ind], Point{3, Float64}[[-4.668111393312408e-18, -0.37854357368761526, -0.6124963684494119], 
+        [-0.22191241796362113, 0.5809742527544328, 0.3590618347908117], [0.0, -0.6134756030005282, -0.37014980570299627], 
+        [0.10988309631672716, -0.681387091313344, 0.19235522107345335], [-0.17398693193931355, 0.3182970619500977, 0.6225453022912995], 
+        [-0.5150154647544891, 0.10752983753681045, -0.4922839938894113]], atol=eps_level)
+    end
+
+    @testset "Loop, constrained boundary" begin
+        # Example with boundary edges (extruded prism)
+        r = 1.0
+        nc = 3
+        Vc = circlepoints(r,nc;dir=:cw)    
+        d = norm(Vc[1]-Vc[2])        
+        s = 1    
+        F,V = extrudecurve(Vc,d;s=s, num_steps=2, close_loop=true,face_type=:tri_slash)
+        
+        Fn, Vn = subtri(F, V, 0; method=:Loop, constrain_boundary=true)
+        @test Fn == F
+        @test Vn == V
+
+        Fn, Vn = subtri(F, V, n; method=:Loop, constrain_boundary=true)
+        ind = round.(Int64,range(1,length(Vn),6)) # Sample indices
+
+        @test eltype(F) == eltype(Fn)
+        @test length(Fn) == 4^n*length(F)        
+        @test eltype(V) == eltype(Vn)        
+        @test isapprox(Vn[ind], Point{3, Float64}[[1.0, 0.0, 0.0], [-0.4999999999999999, 1.1102230246251565e-16, 1.7320508075688776], 
+        [-0.054687500000000125, -0.5277342304311423, 0.8660254037844388], [0.2734375000000002, 0.39241776108982385, 1.2990381056766582], 
+        [-0.3359374999999997, 0.5818608181676699, 1.2990381056766582], [-0.1249999999999997, 0.6495190528383291, 1.7320508075688776]], atol=eps_level)
     end
 
 end
@@ -1449,45 +1481,77 @@ end
     end
 
     @testset "linear" begin
-        Fn, Vn = subquad(F, V, n; method=:linear)
-
-        @test Fn isa Vector{QuadFace{Int64}}
-        @test length(Fn) == 384
-        @test Fn[1] == QuadFace(1, 99, 291, 187)
-        @test Vn isa Vector{Point3{Float64}}
-        @test length(Vn) == 386
-        @test isapprox(Vn[1], [-0.5773502691896258, -0.5773502691896258, -0.5773502691896258], atol=eps_level)
-        
-        Fn, Vn = subquad(F, V, 0)
+        Fn, Vn = subquad(F, V, 0; method=:Linear)
         @test Fn == F
         @test Vn == V
+
+        Fn, Vn = subquad(F, V, n; method=:linear)        
+        ind = round.(Int64,range(1,length(Vn),6)) # Sample indices
+        @test eltype(F) == eltype(Fn)
+        @test length(Fn) == 4^n*length(F)        
+        @test eltype(V) == eltype(Vn)  
+        @test isapprox(Vn[ind], Point{3, Float64}[[-0.5773502691896258, -0.5773502691896258, -0.5773502691896258], 
+        [0.2886751345948129, -0.2886751345948129, -0.5773502691896258], [0.4330127018922194, 0.5773502691896258, -0.5773502691896258], 
+        [-0.5773502691896258, 0.2886751345948129, 0.14433756729740646], [0.14433756729740646, -0.14433756729740646, 0.5773502691896258], 
+        [0.14433756729740646, -0.5773502691896258, -0.43301270189221935]], atol=eps_level)        
     end
 
     @testset "Catmull_Clark" begin
+        Fn, Vn = subquad(F, V, 0; method=:Catmull_Clark)
+        @test Fn == F
+        @test Vn == V
+
         Fn, Vn = subquad(F, V, n; method=:Catmull_Clark)
+        ind = round.(Int64,range(1,length(Vn),6)) # Sample indices
+        @test eltype(F) == eltype(Fn)
+        @test length(Fn) == 4^n*length(F)        
+        @test eltype(V) == eltype(Vn)  
+        @test isapprox(Vn[ind], Point{3, Float64}[[-0.2895661072324513, -0.2895661072324513, -0.2895661072324513], 
+        [0.18414681640860342, -0.18414681640860342, -0.4257509268592806], [0.24244200756986245, 0.31179169810728824, -0.31179169810728824], 
+        [-0.44889359308574917, 0.19141642225574024, 0.09422035643025145], [0.0977285611909523, -0.0977285611909523, 0.47360764269461497], 
+        [0.0907121516695506, -0.407828803431474, -0.2770228830681994]], atol=eps_level)        
+    end
+
+    @testset "Catmull_Clark, constrained boundary" begin
+        # Example with boundary edges (extruded prism)
+        r = 1.0
+        nc = 3
+        Vc = circlepoints(r,nc;dir=:cw)    
+        d = norm(Vc[1]-Vc[2])        
+        s = 1    
+        F,V = extrudecurve(Vc,d;s=s, num_steps=2, close_loop=true,face_type=:quad)
 
-        @test Fn isa Vector{QuadFace{Int64}}
-        @test length(Fn) == 384
-        @test Fn[1] == QuadFace(1, 99, 291, 187)
-        @test Vn isa Vector{Point3{Float64}}
-        @test length(Vn) == 386
-        @test isapprox(Vn[1], [-0.2895661072324513, -0.2895661072324513, -0.2895661072324513], atol=eps_level)
-    end
+        Fn, Vn = subquad(F, V, 0; method=:Catmull_Clark, constrain_boundary=true)
+        @test Fn == F
+        @test Vn == V
 
+        Fn, Vn = subquad(F, V, n; method=:Catmull_Clark, constrain_boundary=true)
+        ind = round.(Int64,range(1,length(Vn),6)) # Sample indices
+        @test eltype(F) == eltype(Fn)
+        @test length(Fn) == 4^n*length(F)        
+        @test eltype(V) == eltype(Vn)  
+        @test isapprox(Vn[ind], Point{3, Float64}[[1.0, 0.0, 0.0], [0.625, -0.21650635094610965, 0.0], 
+        [-0.4516601562499999, 0.16068830734281586, 0.8660254037844388], [-0.4296874999999998, 0.31757083898540706, 1.0825317547305486], 
+        [-0.25976562500000017, -0.6664336115059938, 1.515544456622768], [0.07128906250000025, 0.515894039363777, 0.2165063509461097]], atol=eps_level)        
+    end
 end
 
-@testset "geosphere" begin
+
+@testset "geosphere" begin    
+    eps_level = 1e-4
     r = 1.0
+    
     n = 3
-    F, V = geosphere(n, r)
-    eps_level = 1e-4
-
+    F, V = geosphere(n, r)    
+    ind = round.(Int64,range(1,length(Vn),6)) # Sample indices
     @test F isa Vector{TriangleFace{Int64}}
-    @test length(F) == 1280
-    @test F[1] == [163,323,243]
-    @test V isa Vector{Point3{Float64}}
-    @test length(V) == 642
-    @test isapprox(V[1], [0.0, -0.5257311121191336 , -0.85065080835204], atol=eps_level)
+    @test length(F) == 4^n*20    
+    @test V isa Vector{Point3{Float64}}    
+    @test isapprox(V[ind],Point{3, Float64}[[0.0, -0.5257311121191336, -0.85065080835204], 
+    [-0.5877852522924731, -0.6881909602355868, -0.42532540417602], [0.26286555605956685, -0.16245984811645317, -0.9510565162951538], 
+    [-0.25989191300775444, 0.43388856455269487, 0.8626684804161864], [0.9129824929322992, 0.3996070517018512, 0.08232358003196016], 
+    [0.13279247682790243, -0.2201170274732924, 0.9663925974024391]], atol=eps_level)
+    @test isapprox(norm.(V),fill(r,length(V))) # Test radii
 end