diff --git a/src/mesh.jl b/src/mesh.jl
index 0202180bd55a6462e0cdbe9e1a8f0403f21a0ab5..d400bfd9d55e5c2623452e58883fce25543a30bb 100644
--- a/src/mesh.jl
+++ b/src/mesh.jl
@@ -687,14 +687,28 @@ function MeshInterpolator(smpl_x::AbstractArray{<:Float64}, ref_mesh::Mesh1d)
     # find all sample indices of nodes which fall into the same reference cell
     same_idxs = findall(ref_idx -> ref_idx == cell_idx, ref_cell_indices)
     smpl_idxs = [ current_idx ]
-    if length(same_idxs) == 1
+    if length(same_idxs) != 1
       append!(smpl_idxs, same_idxs[2:end])
     end
     # extract the sample nodes
     smpl_xs = smpl_x[smpl_idxs]
+
+    # generate interpolation matrix
+    box = ref_mesh.boxes[cell_idx]
+    (xmin, xmax), = box.extends
+    ref_xs = @. (2*smpl_xs - (xmax+xmin))/(xmax-xmin)
     # generate interpolation matrix
-    push!(interpolation_matrices,
-          (cell_idx, dg1d.barycentric_interp_matrix(smpl_xs, ref_x[:,cell_idx])))
+    if ref_mesh.element.kind === :modal_bspline1
+      push!(interpolation_matrices,
+            (cell_idx, dg1d.Bspline.vandermonde_matrix(ref_xs, ref_mesh.element.bs1)))
+    elseif ref_mesh.element.kind === :modal_bsplinea2
+      push!(interpolation_matrices,
+            (cell_idx, dg1d.Bspline.vandermonde_matrix(ref_xs, ref_mesh.element.bs2)))
+    else
+      push!(interpolation_matrices,
+            (cell_idx, dg1d.barycentric_interp_matrix(smpl_xs, ref_x[:,cell_idx])))
+    end
+
 
     current_idx += length(smpl_idxs)
     if current_idx > length(ref_cell_indices)