diff --git a/dependencies/conda_env.yml b/dependencies/conda_env.yml
index efeb7214..560938bc 100644
--- a/dependencies/conda_env.yml
+++ b/dependencies/conda_env.yml
@@ -8,9 +8,7 @@ channels:
 dependencies:
   - python>=3.9
   - cmake
-  - k3d
-  # - numpy=1.23 #Temporary error with numba calling numpy calling ufunc with no exceptions in latest numpy; see https://github.com/numba/numba/issues/8615
-  - vtk
+  - numpy
   - jupyterlab
   - yt
   - scipy<=1.13
@@ -18,7 +16,6 @@ dependencies:
   - plotly
   - pyyaml
   - ipywidgets
-  - meshio=5.2
   - tqdm
   - palettable
   - h5py
@@ -31,8 +28,11 @@ dependencies:
   - pandas
   - nodejs
   - twine
-  - sphinx=4.4.0
-  - breathe
+  - sphinx #=4.4.0
+  - pip
+  - pip:
+    - breathe
+  # - breathe #install with pip, otherwise it will install the wrong version, not compatible with sphinx
   - sphinx-copybutton
   - nbsphinx
   - sphinx_rtd_theme
diff --git a/docs/src/0_getting_started/0_prerequisites.rst b/docs/src/0_getting_started/0_prerequisites.rst
index 9f4a8673..e035e9a1 100644
--- a/docs/src/0_getting_started/0_prerequisites.rst
+++ b/docs/src/0_getting_started/0_prerequisites.rst
@@ -11,11 +11,11 @@ Below is a list of the software required to be able to compile and run Magritte.
 * `GCC <https://gcc.gnu.org/>`_ (version :literal:`5.0.0` or later): to compile the C++ part of Magritte.
 * `CMake <https://cmake.org/>`_ (version :literal:`3.18.0` or later): for building Magritte, organising compilation and linking;
 * `Anaconda <https://www.anaconda.com/blog/individual-edition-2020-11>`_ (optional): for managing the required Python packages;
-* `Gmsh <https://gmsh.info/>`_ (optional): for constructing geometric meshes for the models;
+.. * `Gmsh <https://gmsh.info/>`_ (optional): for constructing geometric meshes for the models;
 
 Although Anaconda is optional and you could, in priciple, use any other way of installing the required Python dependencies, we recommend and assume the use of Anaconda.
 The list of required Python packages can be found in the `conda environment file <https://github.com/Magritte-code/Magritte/blob/stable/dependencies/conda_env.yml>`_.
-Also Gmsh is optional, however, it will be used in most examples in this documentation to create models.
+.. Also Gmsh is optional, however, it will be used in most examples in this documentation to create models.
 
 Once these prerequisites are in place, Magritte can be installed following the :ref:`quickstart <link-quickstart>`.
 
@@ -25,9 +25,8 @@ Once these prerequisites are in place, Magritte can be installed following the :
 
 
 .. Hint::
-    Although we recommend to first try to install the required software in your preferred usual way, the folder `Magritte/dependencies/ <https://github.com/Magritte-code/Magritte/tree/stable/dependencies>`_ contains a set of bash scripts: :literal:`get_conda.sh`, :literal:`get_cmake.sh`,
-    and :literal:`get_gmsh.sh`, to quickly obtain the appropriate versions of miniconda3 (for Anaconda), CMake, and
-    Gmsh respectively. However, these scripts don't properly install the software, but rather download an executable binary, which is useful, e.g. if you don't have root access.
+    Although we recommend to first try to install the required software in your preferred usual way, the folder `Magritte/dependencies/ <https://github.com/Magritte-code/Magritte/tree/stable/dependencies>`_ contains a set of bash scripts: :literal:`get_conda.sh` and :literal:`get_cmake.sh`,
+    to quickly obtain the appropriate versions of miniconda3 (for Anaconda) and CMake respectively. However, these scripts don't properly install the software, but rather download an executable binary, which is useful, e.g. if you don't have root access.
     
     To use these, first clone the Magritte `GitHub <https://github.com/Magritte-code/Magritte>`_ repository
     
@@ -44,10 +43,6 @@ Once these prerequisites are in place, Magritte can be installed following the :
     .. code-block:: shell
     
         bash dependencies/get_cmake.sh
-        
-    .. code-block:: shell
-    
-        bash dependencies/get_gmsh.sh
 
     ⚠️  Don't forget to include the relevant paths to your shell's :literal:`$PATH` variable!
     (For convenience, the scripts will print the relevant path.)
diff --git a/docs/src/0_getting_started/2_installation.rst b/docs/src/0_getting_started/2_installation.rst
index d4991485..ac006b0d 100644
--- a/docs/src/0_getting_started/2_installation.rst
+++ b/docs/src/0_getting_started/2_installation.rst
@@ -49,7 +49,6 @@ Magritte has several dependencies, some of which are optional.
 **Optional**
 
 * `Anaconda <https://www.anaconda.com/blog/individual-edition-2020-11>`_, for managing the required Python packages;
-* `Gmsh <https://gmsh.info/>`_, version :literal:`4.6.0` or later, for meshing model geometries.
 
 
 Please note that :literal:`Paracabs` might have further dependencies depending
@@ -71,8 +70,6 @@ the model files.
 * :mod:`scipy`, for interpolation and spatial functions such as nearest neighbour calculations;
 * :mod:`healpy`, to sample directions from a discretized unit sphere;
 * :mod:`astropy`, for unit conversions and physical constants;
-* :mod:`meshio`, for reading and writing several types of mesh data structures;
-* :mod:`vtk`, for reading and writing vtk files;
 * :mod:`pyyaml`, for reading and writing yaml files;
 * :mod:`mpi4py`, for MPI (Message Passing Interface) functionality in Python;
 * :mod:`tqdm`, for progress bars;
diff --git a/magritte/mesher.py b/magritte/mesher.py
index a9f34b2b..74059e8a 100644
--- a/magritte/mesher.py
+++ b/magritte/mesher.py
@@ -1,5 +1,5 @@
 import sys, os
-import meshio
+# import meshio
 import numpy     as np
 import itertools as itt
 import numba
@@ -25,600 +25,600 @@ class Mesh:
     """
     Magritte custom mesh class (that ignores unconnected points).
     """
-    def __init__(self, meshFile):
-        """
-        Read the meshFile using meshio and remove unconnected points.
-
-        Parameters
-        ----------
-        meshFile : str
-            File name of the mesh file.
-        """
-        self.mesh      = meshio.read(meshFile)
-        self.points    = self.mesh.points
-        self.tetras    = self.mesh.cells_dict['tetra']
-        self.edges     = self.get_edges()
-        self.neighbors = self.get_neighbors()
-        self.boundary  = self.get_boundary()
-
-        # Remove non-connected points
-        non_point = self.get_non_point()
-        while not (non_point is None):
-            print(non_point)
-            self.del_non_point(non_point)
-            non_point = self.get_non_point()
-
-        return
-
-
-    def get_edges(self):
-        edges = set([])
-        for tetra in self.tetras:
-            for (A,B) in itt.combinations(tetra,2):
-                edges.add((A,B))
-        return np.array(list(edges))
-
-    def get_non_point(self):
-        for (p,nl) in enumerate(self.neighbors):
-            if (len(nl) == 0):
-                return p
-
-    def del_non_point(self, p):
-        # remove the non-connected point
-        self.points    = np.delete(self.points,    p, 0)
-        self.neighbors.pop(p)
-        # Change all indices, since a point got removed
-        self.tetras    = relocate_indices(self.tetras,    p)
-        self.edges     = relocate_indices(self.edges,     p)
-        # self.neighbors = relocate_indices(self.neighbors, p)
-
-    # Does not work properly, crashes Magritte.
-    def get_neighbors(self):
-        neighbors = [[] for _ in range(len(self.points))]
-        for edge in self.edges:
-            neighbors[edge[0]].append(edge[1])
-            neighbors[edge[1]].append(edge[0])
-        return neighbors
-
-    def get_tetra_volume(self, tetra):
-        a = self.points[tetra[0]]
-        b = self.points[tetra[1]]
-        c = self.points[tetra[2]]
-        d = self.points[tetra[3]]
-        return np.abs(np.dot(a-d,np.cross(b-d,c-d)))/6
-
-    def get_tetra_volumes(self):
-        return [self.get_tetra_volume(tetra) for tetra in self.tetras]
-
-    def get_tetra_position(self, tetra):
-        a = self.points[tetra[0]]
-        b = self.points[tetra[1]]
-        c = self.points[tetra[2]]
-        d = self.points[tetra[3]]
-        return 0.25*(a+b+c+d)
-
-    def get_tetra_positions(self):
-        return [self.get_tetra_position(tetra) for tetra in self.tetras]
-
-    def get_edge_length(self, edge):
-        a = self.points[edge[0]]
-        b = self.points[edge[1]]
-        r = b - a
-        return np.sqrt(np.dot(r,r))
-
-    def get_edge_lengths(self):
-        return [self.get_edge_length(edge) for edge in self.edges]
-
-    def get_edge_position(self, edge):
-        a = self.points[edge[0]]
-        b = self.points[edge[1]]
-        return 0.5*(a+b)
-
-    def get_edge_positions(self):
-        return [self.get_edge_position(edge) for edge in self.edges]
-
-    def get_boundary(self):
-        boundary = set([])
-        for elem in self.mesh.cells_dict['triangle']:
-            for p in elem:
-                boundary.add(p)
-        for elem in self.mesh.cells_dict['line']:
-            for p in elem:
-                boundary.add(p)
-        for elem in self.mesh.cells_dict['vertex']:
-            for p in elem:
-                boundary.add(p)
-        boundary = list(boundary)
-        return boundary
-
-
-def run(command):
-    """
-    Run command in shell and continuously print its output.
-
-    Parameters
-    ----------
-    command : str
-        The command to run in the shell.
-    """
-    # Run command pipe output
-    process = Popen(command, stdout=PIPE, shell=True)
-
-    # Continuously read output and print
-    while True:
-        # Extract output
-        line = process.stdout.readline().rstrip()
-        # Break if there is no more line, print otherwise
-        if not line:
-            break
-        else:
-            print(line.decode("utf-8"))
-
-    return
-
-
-def convert_msh_to_pos(meshName, replace:bool=False):
-    """
-    Convert a .msh file to a .pos file.
-
-    Parameters
-    ----------
-    modelName : str
-        Path to the .msh file to convert.
-    replace : bool
-        Whether or not to remove (and thus replace) the original .msh file.
-    """
-    # Remove extension from meshName
-    meshName, extension = os.path.splitext(meshName)
-
-    # create the converision gmsh script file
-    conversion_script = f'{meshName}_convert_to_pos.geo'
-
-    # Get the path to this folder
-    thisFolder = os.path.dirname(os.path.abspath(__file__))
-
-    # Extract the template
-    with open(f'{thisFolder}/templates/convert_to_pos.template', 'r') as file:
-        template = Template(file.read())
-
-    # Fill out the template and write the conversion script
-    with open(conversion_script,                                 'w') as file:
-        file.write(template.substitute(FILE_NAME=meshName))
-
-    # Run gmsh in a subprocess to convert the background mesh to the .pos format
-    run(f'gmsh -0 {conversion_script}')
-
-    # Remove the auxiliary file that is created (geo_unrolled) and the script file
-    os.remove(f'{meshName}_convert_to_pos.geo_unrolled')
-    os.remove(conversion_script)
-    if replace:
-        os.remove(f"{meshName}.msh")
-
-    return
-
-
-def boundary_cuboid (minVec, maxVec):
-    """
-    Retuns the gmsh script for a cuboid element that can be used as boundary.
-
-    Parameters
-    ----------
-    minVec : array_like
-        (x_min, y_min, z_min) vector defining the cuboid.
-    maxVec : array_like
-        (x_max, y_max, z_max) vector defining the cuboid.
-
-    Returns
-    -------
-    out : str
-        A string containing the gmsh script defining the cuboid.
-    """
-    # Get the path to this folder
-    thisFolder = os.path.dirname(os.path.abspath(__file__))
-
-    # Extract the template
-    with open(f'{thisFolder}/templates/cuboid.template', 'r') as file:
-        cuboid = Template(file.read())
-
-    # Return the filled out template
-    return cuboid.substitute(
-               I     = 1,
-               X_MIN = minVec[0],
-               X_MAX = maxVec[0],
-               Y_MIN = minVec[1],
-               Y_MAX = maxVec[1],
-               Z_MIN = minVec[2],
-               Z_MAX = maxVec[2] )
-
-
-def boundary_sphere (centre=np.zeros(3), radius=1.0):
-    """
-    Returns the gmsh script for a sphere element that can be used as boundary.
-
-    Parameters
-    ----------
-    centre : array_like
-        (x, y, z) coordinates of the centre of the sphere.
-    radius : float
-        Radius of the sphere.
-
-    Returns
-    -------
-    out : str
-        A string containing the gmsh script defining the sphere.
-    """
-    # Get the path to this folder
-    thisFolder = os.path.dirname(os.path.abspath(__file__))
-
-    # Extract the template
-    with open(f'{thisFolder}/templates/sphere.template', 'r') as file:
-        sphere = Template(file.read())
-
-    # Return the filled out template
-    return sphere.substitute(
-               I      = 1,
-               CX     = centre[0],
-               CY     = centre[1],
-               CZ     = centre[2],
-               RADIUS = radius    )
-
-
-def boundary_sphere_in_sphere (centre_in =np.zeros(3), radius_in =1.0,
-                               centre_out=np.zeros(3), radius_out=1.0 ):
-    """
-    Returns the gmsh script for a sphere inside a sphere that can be used as boundary.
-
-    Parameters
-    ----------
-    centre_in : array_like
-        (x, y, z) coordinates of the centre of the inner sphere.
-    radius_in : float
-        Radius of the inner sphere.
-    centre_out : array_like
-        (x, y, z) coordinates of the centre of the outer sphere.
-    radius_out : float
-        Radius of the outer sphere.
-
-    Returns
-    -------
-    out : str
-        A string containing the gmsh script defining the sphere in sphere.
-    """
-    # Get the path to this folder
-    thisFolder = os.path.dirname(os.path.abspath(__file__))
-
-    # Extract the template
-    with open(f'{thisFolder}/templates/sphere_in_sphere.template', 'r') as file:
-        sphere = Template(file.read())
-
-    # Return the filled out template
-    return sphere.substitute(
-               CX_IN      = centre_in[0],
-               CY_IN      = centre_in[1],
-               CZ_IN      = centre_in[2],
-               RADIUS_IN  = radius_in,
-               CX_OUT     = centre_out[0],
-               CY_OUT     = centre_out[1],
-               CZ_OUT     = centre_out[2],
-               RADIUS_OUT = radius_out    )
-
-
-def boundary_sphere_in_cuboid (centre_in=np.zeros(3), radius_in=1.0,
-                               minVec   =np.zeros(3), maxVec   =np.ones(3)):
-    """
-    Returns the gmsh script for a sphere inside a cuboid that can be used as boundary.
-
-    Parameters
-    ----------
-    centre_in : array_like
-        (x, y, z) coordinates of the centre of the sphere.
-    radius_in : float
-        Radius of the sphere.
-    minVec : array_like
-        (x_min, y_min, z_min) vector defining the cuboid.
-    maxVec : array_like
-        (x_max, y_max, z_max) vector defining the cuboid.
-
-    Returns
-    -------
-    out : str
-        A string containing the gmsh script defining the sphere in cuboid.
-    """
-    # Get the path to this folder
-    thisFolder = os.path.dirname(os.path.abspath(__file__))
-
-    # Extrace the template
-    with open(f'{thisFolder}/templates/sphere_in_cuboid.template', 'r') as file:
-        sphere = Template(file.read())
-
-    # Return the filled out template
-    return sphere.substitute(
-        CX_IN      = centre_in[0],
-        CY_IN      = centre_in[1],
-        CZ_IN      = centre_in[2],
-        RADIUS_IN  = radius_in,
-        X_MIN      = minVec[0],
-        X_MAX      = maxVec[0],
-        Y_MIN      = minVec[1],
-        Y_MAX      = maxVec[1],
-        Z_MIN      = minVec[2],
-        Z_MAX      = maxVec[2]   )
-
-
-def create_mesh_from_background(meshName, boundary, scale_min, scale_max):
-    """
-    Creates a mesh with element sizes defined on a background mesh.
-
-    Parameters
-    ----------
-    meshName : str
-        File name of the mesh file.
-    boundary : str
-        Gmsh script containing the definition of the boundary.
-    scale_min : float
-        Minimum desired element size of the mesh.
-    scale_max : float
-        Maximum desired element size of the mesh.
-    """
-    # Remove extension from meshName
-    meshName, extension = os.path.splitext(meshName)
-
-    # create the mesh generating gmsh script file
-    meshing_script = f'{meshName}.geo'
-    background     = f'{meshName}.pos'
-    resulting_mesh = f'{meshName}.vtk'
-
-    # Get the path to this folder
-    thisFolder = os.path.dirname(os.path.abspath(__file__))
-
-    # Extract the template
-    with open(f'{thisFolder}/templates/mesh_from_background.template', 'r') as file:
-        template = Template(file.read())
-
-    # Write the gmsh script
-    with open(meshing_script, 'w') as file:
-        file.write(template.substitute(
-            BOUNDARY   = boundary,
-            SCALE_MIN  = scale_min,
-            SCALE_MAX  = scale_max,
-            BACKGROUND = background   ))
-
-    # run gmsh in a subprocess to generate the mesh from the background
-    run(f'gmsh {meshing_script} -3 -saveall -o {resulting_mesh}')
-
-    # Remove the script file
-    os.remove(meshing_script)
-
-    return
-
-
-def create_mesh_from_function(meshName, boundary, scale_min, scale_max, scale_function):
-    """
-    Creates a mesh with element sizes defined by a scale function.
-
-    Parameters
-    ----------
-    meshName : str
-        File name of the mesh file.
-    boundary : str
-        Gmsh script containing the definition of the boundary.
-    scale_min : float
-        Minimum desired element size of the mesh.
-    scale_max : float
-        Maximum desired element size of the mesh.
-    scale_function : str
-        Formula for the function describing the desired mesh element sizes.
-    """
-    # Remove extension from meshName
-    meshName, extension = os.path.splitext(meshName)
-
-    # create the mesh generating gmsh script file
-    meshing_script = f'{meshName}.geo'
-    resulting_mesh = f'{meshName}.vtk'
-
-    # Get the path to this folder
-    thisFolder = os.path.dirname(os.path.abspath(__file__))
-
-    # Extract the template
-    with open(f'{thisFolder}/templates/mesh_from_function.template', 'r') as file:
-        template = Template(file.read())
-
-    # Write the gmsh script
-    with open(meshing_script, 'w') as file:
-        file.write(template.substitute(
-            BOUNDARY       = boundary,
-            SCALE_MIN      = scale_min,
-            SCALE_MAX      = scale_max,
-            SCALE_FUNCTION = scale_function ))
-
-    # Run gmsh in a subprocess to generate the mesh from the background
-    run(f'gmsh {meshing_script} -3 -saveall -o {resulting_mesh}')
-
-    # Remove the script file
-    os.remove(meshing_script)
-
-    return
-
-
-def generate_background_from_1D_data(meshName, R, data):
-    """
-    Generate a background mesh from 1D (radial) data.
-
-    Parameters
-    ----------
-    meshName : str
-        File name of the mesh file.
-    R : array_type
-        Array containing the radial positions at which the data is given.
-    data : array_type
-        Array containing the data to put a background.
-    """
-    # Remove extension from meshName
-    meshName, extension = os.path.splitext(meshName)
-
-    # Map the data to a set of concentric HEALPix spheres
-    nsides = 8
-    sphere = np.array(pixelfunc.pix2vec(nsides, range(12*nsides**2))).transpose()
-    points = []
-    data_s = []
-    for (i,r) in enumerate(R):
-        for s in Rotation.random().apply(sphere):
-            points.append(r*s)
-            data_s.append(data[i])
-
-    # Delaunay tetrahedralise the point set
-    delaunay = Delaunay(points)
-
-    # Write the mesh
-    meshio.write_points_cells(
-        filename   = f'{meshName}.msh',
-        points     =             delaunay.points,
-        cells      = {'tetra'  : delaunay.simplices},
-        point_data = {'weights': np.array(data_s)}   )
-
-    return
-
-
-@njit
-def get_pfs(nns):
-    pfs = np.zeros(nns.shape, dtype=np.int64)
-    tot = 0
-    for i, n in enumerate(nns):
-        tot    += n
-        pfs[i]  = tot
-    return pfs
-
-
-@njit
-def get_nbs(nbs, pfs, i):
-    if (i == 0):
-        return nbs[0:pfs[0]]
-    else:
-        return nbs[pfs[i-1]:pfs[i]]
-
-
-@njit
-def norm(vec_arr):
-    norms = np.zeros(vec_arr.shape[0])
-    for i in range(vec_arr.shape[0]):
-        for j in range(vec_arr.shape[1]):
-            norms[i] = norms[i] + vec_arr[i,j]**2
-    return np.sqrt(norms)
-
-
-@njit
-def get_Gs(tracer, nbs, pfs):
-    """
-    Get the maximum relative difference in the tracer function for each point.
-    """
-    Gs = np.zeros(tracer.shape[0])
-    for i in range(tracer.shape[0]):
-        tracer_i = tracer[i]
-        tracer_n = tracer[get_nbs(nbs, pfs, i)]
-        Gs[i] = np.max(np.abs(tracer_i - tracer_n) / (tracer_i + tracer_n))
-    return Gs
-
-
-@njit
-def get_Ls(pos, nbs, pfs):
-    """
-    Get the element size distributions for each point.
-    """
-    Ls = np.zeros(pos.shape[0])
-    for i in range(pos.shape[0]):
-        Ls[i] = np.mean(norm(pos[i] - pos[get_nbs(nbs, pfs, i)]))
-    return Ls
-
-
-@njit
-def get_weights(pos, nbs, nns, tracer, threshold=0.21, fmin=1.0, ftarget=2.15):
-    """
-    Get the weights (desired element sizes) for each point.
-    """
-    pfs = get_pfs (nns)
-    GGs = get_Gs  (tracer, nbs, pfs)
-    LLs = get_Ls  (pos,    nbs, pfs)
-
-    L_min    = fmin    * LLs
-    L_target = ftarget * LLs
-
-    weights = L_target
-    weights[GGs > threshold] = L_min[GGs > threshold]
-
-    return weights
-
-
-def reduce(
-        meshName,
-        delaunay,
-        tracer,
-        boundary,
-        scale_max = 1.0e+99,  # Maximum value of the scale parameter
-        scale_min = 0.0e+00,  # Minimum value of the scale parameter
-        threshold = 0.21,     # Threashold relative difference for coarsening
-        fmin      = 1.0,      # Don't allow refinenment
-        ftarget   = 2.15      # Approx 10^(-1/3) for approx 10 times fewer points
-    ):
-    """
-    Reduce a model mesh for a given tracer function.
-
-    Parameters
-    ----------
-    meshName : str
-        File name of the mesh file.
-    delaunay :
-        scipy.spatial.Delaunay object for the given mesh.
-    tracer : array_like
-        1D array containing the tabulated tracer function.
-    boundary : str
-        Gmsh script containing the definition of the boundary.
-    scale_max : float
-        Maximum desired element size of the mesh.
-    scale_min : float
-        Minimum desired element size of the mesh.
-    threshold : float
-        Threashold value for relative difference below which we coarsen.
-    fmin : float
-        Multiplication factor of the mesh element size where the relative difference is below the threashold.
-    ftarget : float
-        Multiplication factor of the mesh element size where the relative difference is below the threashold.
-    """
-    # Remove extension from meshName
-    meshName, extension = os.path.splitext(meshName)
-
-    # Extract mesh structure from delaunay object
-    (indptr, indices) = delaunay.vertex_neighbor_vertices
-    neighbors = [indices[indptr[k]:indptr[k+1]] for k in range(len(delaunay.points))]
-    nbs       = np.array([n for sublist in neighbors for n in sublist])
-    nns       = np.array([len(sublist) for sublist in neighbors])
-    pos       = delaunay.points
-
-    # Create a background mesh (in .msh format)
-    meshio.write_points_cells(
-        filename   = f'{meshName}.msh',
-        points     = delaunay.points,
-        cells      = {'tetra'  : delaunay.simplices},
-        point_data = {'weights': get_weights(pos=pos,
-                                             nbs=nbs,
-                                             nns=nns,
-                                             tracer=tracer,
-                                             threshold=threshold,
-                                             fmin=fmin,
-                                             ftarget=ftarget)}
-    )
-
-    # Convert .msh to .pos mesh for Gmsh
-    convert_msh_to_pos (meshName=meshName, replace=True)
-
-    # Create a new mesh from the background mesh
-    create_mesh_from_background(
-        meshName  = meshName,
-        boundary  = boundary,
-        scale_min = scale_min,
-        scale_max = scale_max
-    )
-
-    return
+#     def __init__(self, meshFile):
+#         """
+#         Read the meshFile using meshio and remove unconnected points.
+
+#         Parameters
+#         ----------
+#         meshFile : str
+#             File name of the mesh file.
+#         """
+#         self.mesh      = meshio.read(meshFile)
+#         self.points    = self.mesh.points
+#         self.tetras    = self.mesh.cells_dict['tetra']
+#         self.edges     = self.get_edges()
+#         self.neighbors = self.get_neighbors()
+#         self.boundary  = self.get_boundary()
+
+#         # Remove non-connected points
+#         non_point = self.get_non_point()
+#         while not (non_point is None):
+#             print(non_point)
+#             self.del_non_point(non_point)
+#             non_point = self.get_non_point()
+
+#         return
+
+
+#     def get_edges(self):
+#         edges = set([])
+#         for tetra in self.tetras:
+#             for (A,B) in itt.combinations(tetra,2):
+#                 edges.add((A,B))
+#         return np.array(list(edges))
+
+#     def get_non_point(self):
+#         for (p,nl) in enumerate(self.neighbors):
+#             if (len(nl) == 0):
+#                 return p
+
+#     def del_non_point(self, p):
+#         # remove the non-connected point
+#         self.points    = np.delete(self.points,    p, 0)
+#         self.neighbors.pop(p)
+#         # Change all indices, since a point got removed
+#         self.tetras    = relocate_indices(self.tetras,    p)
+#         self.edges     = relocate_indices(self.edges,     p)
+#         # self.neighbors = relocate_indices(self.neighbors, p)
+
+#     # Does not work properly, crashes Magritte.
+#     def get_neighbors(self):
+#         neighbors = [[] for _ in range(len(self.points))]
+#         for edge in self.edges:
+#             neighbors[edge[0]].append(edge[1])
+#             neighbors[edge[1]].append(edge[0])
+#         return neighbors
+
+#     def get_tetra_volume(self, tetra):
+#         a = self.points[tetra[0]]
+#         b = self.points[tetra[1]]
+#         c = self.points[tetra[2]]
+#         d = self.points[tetra[3]]
+#         return np.abs(np.dot(a-d,np.cross(b-d,c-d)))/6
+
+#     def get_tetra_volumes(self):
+#         return [self.get_tetra_volume(tetra) for tetra in self.tetras]
+
+#     def get_tetra_position(self, tetra):
+#         a = self.points[tetra[0]]
+#         b = self.points[tetra[1]]
+#         c = self.points[tetra[2]]
+#         d = self.points[tetra[3]]
+#         return 0.25*(a+b+c+d)
+
+#     def get_tetra_positions(self):
+#         return [self.get_tetra_position(tetra) for tetra in self.tetras]
+
+#     def get_edge_length(self, edge):
+#         a = self.points[edge[0]]
+#         b = self.points[edge[1]]
+#         r = b - a
+#         return np.sqrt(np.dot(r,r))
+
+#     def get_edge_lengths(self):
+#         return [self.get_edge_length(edge) for edge in self.edges]
+
+#     def get_edge_position(self, edge):
+#         a = self.points[edge[0]]
+#         b = self.points[edge[1]]
+#         return 0.5*(a+b)
+
+#     def get_edge_positions(self):
+#         return [self.get_edge_position(edge) for edge in self.edges]
+
+#     def get_boundary(self):
+#         boundary = set([])
+#         for elem in self.mesh.cells_dict['triangle']:
+#             for p in elem:
+#                 boundary.add(p)
+#         for elem in self.mesh.cells_dict['line']:
+#             for p in elem:
+#                 boundary.add(p)
+#         for elem in self.mesh.cells_dict['vertex']:
+#             for p in elem:
+#                 boundary.add(p)
+#         boundary = list(boundary)
+#         return boundary
+
+
+# def run(command):
+#     """
+#     Run command in shell and continuously print its output.
+
+#     Parameters
+#     ----------
+#     command : str
+#         The command to run in the shell.
+#     """
+#     # Run command pipe output
+#     process = Popen(command, stdout=PIPE, shell=True)
+
+#     # Continuously read output and print
+#     while True:
+#         # Extract output
+#         line = process.stdout.readline().rstrip()
+#         # Break if there is no more line, print otherwise
+#         if not line:
+#             break
+#         else:
+#             print(line.decode("utf-8"))
+
+#     return
+
+
+# def convert_msh_to_pos(meshName, replace:bool=False):
+#     """
+#     Convert a .msh file to a .pos file.
+
+#     Parameters
+#     ----------
+#     modelName : str
+#         Path to the .msh file to convert.
+#     replace : bool
+#         Whether or not to remove (and thus replace) the original .msh file.
+#     """
+#     # Remove extension from meshName
+#     meshName, extension = os.path.splitext(meshName)
+
+#     # create the converision gmsh script file
+#     conversion_script = f'{meshName}_convert_to_pos.geo'
+
+#     # Get the path to this folder
+#     thisFolder = os.path.dirname(os.path.abspath(__file__))
+
+#     # Extract the template
+#     with open(f'{thisFolder}/templates/convert_to_pos.template', 'r') as file:
+#         template = Template(file.read())
+
+#     # Fill out the template and write the conversion script
+#     with open(conversion_script,                                 'w') as file:
+#         file.write(template.substitute(FILE_NAME=meshName))
+
+#     # Run gmsh in a subprocess to convert the background mesh to the .pos format
+#     run(f'gmsh -0 {conversion_script}')
+
+#     # Remove the auxiliary file that is created (geo_unrolled) and the script file
+#     os.remove(f'{meshName}_convert_to_pos.geo_unrolled')
+#     os.remove(conversion_script)
+#     if replace:
+#         os.remove(f"{meshName}.msh")
+
+#     return
+
+
+# def boundary_cuboid (minVec, maxVec):
+#     """
+#     Retuns the gmsh script for a cuboid element that can be used as boundary.
+
+#     Parameters
+#     ----------
+#     minVec : array_like
+#         (x_min, y_min, z_min) vector defining the cuboid.
+#     maxVec : array_like
+#         (x_max, y_max, z_max) vector defining the cuboid.
+
+#     Returns
+#     -------
+#     out : str
+#         A string containing the gmsh script defining the cuboid.
+#     """
+#     # Get the path to this folder
+#     thisFolder = os.path.dirname(os.path.abspath(__file__))
+
+#     # Extract the template
+#     with open(f'{thisFolder}/templates/cuboid.template', 'r') as file:
+#         cuboid = Template(file.read())
+
+#     # Return the filled out template
+#     return cuboid.substitute(
+#                I     = 1,
+#                X_MIN = minVec[0],
+#                X_MAX = maxVec[0],
+#                Y_MIN = minVec[1],
+#                Y_MAX = maxVec[1],
+#                Z_MIN = minVec[2],
+#                Z_MAX = maxVec[2] )
+
+
+# def boundary_sphere (centre=np.zeros(3), radius=1.0):
+#     """
+#     Returns the gmsh script for a sphere element that can be used as boundary.
+
+#     Parameters
+#     ----------
+#     centre : array_like
+#         (x, y, z) coordinates of the centre of the sphere.
+#     radius : float
+#         Radius of the sphere.
+
+#     Returns
+#     -------
+#     out : str
+#         A string containing the gmsh script defining the sphere.
+#     """
+#     # Get the path to this folder
+#     thisFolder = os.path.dirname(os.path.abspath(__file__))
+
+#     # Extract the template
+#     with open(f'{thisFolder}/templates/sphere.template', 'r') as file:
+#         sphere = Template(file.read())
+
+#     # Return the filled out template
+#     return sphere.substitute(
+#                I      = 1,
+#                CX     = centre[0],
+#                CY     = centre[1],
+#                CZ     = centre[2],
+#                RADIUS = radius    )
+
+
+# def boundary_sphere_in_sphere (centre_in =np.zeros(3), radius_in =1.0,
+#                                centre_out=np.zeros(3), radius_out=1.0 ):
+#     """
+#     Returns the gmsh script for a sphere inside a sphere that can be used as boundary.
+
+#     Parameters
+#     ----------
+#     centre_in : array_like
+#         (x, y, z) coordinates of the centre of the inner sphere.
+#     radius_in : float
+#         Radius of the inner sphere.
+#     centre_out : array_like
+#         (x, y, z) coordinates of the centre of the outer sphere.
+#     radius_out : float
+#         Radius of the outer sphere.
+
+#     Returns
+#     -------
+#     out : str
+#         A string containing the gmsh script defining the sphere in sphere.
+#     """
+#     # Get the path to this folder
+#     thisFolder = os.path.dirname(os.path.abspath(__file__))
+
+#     # Extract the template
+#     with open(f'{thisFolder}/templates/sphere_in_sphere.template', 'r') as file:
+#         sphere = Template(file.read())
+
+#     # Return the filled out template
+#     return sphere.substitute(
+#                CX_IN      = centre_in[0],
+#                CY_IN      = centre_in[1],
+#                CZ_IN      = centre_in[2],
+#                RADIUS_IN  = radius_in,
+#                CX_OUT     = centre_out[0],
+#                CY_OUT     = centre_out[1],
+#                CZ_OUT     = centre_out[2],
+#                RADIUS_OUT = radius_out    )
+
+
+# def boundary_sphere_in_cuboid (centre_in=np.zeros(3), radius_in=1.0,
+#                                minVec   =np.zeros(3), maxVec   =np.ones(3)):
+#     """
+#     Returns the gmsh script for a sphere inside a cuboid that can be used as boundary.
+
+#     Parameters
+#     ----------
+#     centre_in : array_like
+#         (x, y, z) coordinates of the centre of the sphere.
+#     radius_in : float
+#         Radius of the sphere.
+#     minVec : array_like
+#         (x_min, y_min, z_min) vector defining the cuboid.
+#     maxVec : array_like
+#         (x_max, y_max, z_max) vector defining the cuboid.
+
+#     Returns
+#     -------
+#     out : str
+#         A string containing the gmsh script defining the sphere in cuboid.
+#     """
+#     # Get the path to this folder
+#     thisFolder = os.path.dirname(os.path.abspath(__file__))
+
+#     # Extrace the template
+#     with open(f'{thisFolder}/templates/sphere_in_cuboid.template', 'r') as file:
+#         sphere = Template(file.read())
+
+#     # Return the filled out template
+#     return sphere.substitute(
+#         CX_IN      = centre_in[0],
+#         CY_IN      = centre_in[1],
+#         CZ_IN      = centre_in[2],
+#         RADIUS_IN  = radius_in,
+#         X_MIN      = minVec[0],
+#         X_MAX      = maxVec[0],
+#         Y_MIN      = minVec[1],
+#         Y_MAX      = maxVec[1],
+#         Z_MIN      = minVec[2],
+#         Z_MAX      = maxVec[2]   )
+
+
+# def create_mesh_from_background(meshName, boundary, scale_min, scale_max):
+#     """
+#     Creates a mesh with element sizes defined on a background mesh.
+
+#     Parameters
+#     ----------
+#     meshName : str
+#         File name of the mesh file.
+#     boundary : str
+#         Gmsh script containing the definition of the boundary.
+#     scale_min : float
+#         Minimum desired element size of the mesh.
+#     scale_max : float
+#         Maximum desired element size of the mesh.
+#     """
+#     # Remove extension from meshName
+#     meshName, extension = os.path.splitext(meshName)
+
+#     # create the mesh generating gmsh script file
+#     meshing_script = f'{meshName}.geo'
+#     background     = f'{meshName}.pos'
+#     resulting_mesh = f'{meshName}.vtk'
+
+#     # Get the path to this folder
+#     thisFolder = os.path.dirname(os.path.abspath(__file__))
+
+#     # Extract the template
+#     with open(f'{thisFolder}/templates/mesh_from_background.template', 'r') as file:
+#         template = Template(file.read())
+
+#     # Write the gmsh script
+#     with open(meshing_script, 'w') as file:
+#         file.write(template.substitute(
+#             BOUNDARY   = boundary,
+#             SCALE_MIN  = scale_min,
+#             SCALE_MAX  = scale_max,
+#             BACKGROUND = background   ))
+
+#     # run gmsh in a subprocess to generate the mesh from the background
+#     run(f'gmsh {meshing_script} -3 -saveall -o {resulting_mesh}')
+
+#     # Remove the script file
+#     os.remove(meshing_script)
+
+#     return
+
+
+# def create_mesh_from_function(meshName, boundary, scale_min, scale_max, scale_function):
+#     """
+#     Creates a mesh with element sizes defined by a scale function.
+
+#     Parameters
+#     ----------
+#     meshName : str
+#         File name of the mesh file.
+#     boundary : str
+#         Gmsh script containing the definition of the boundary.
+#     scale_min : float
+#         Minimum desired element size of the mesh.
+#     scale_max : float
+#         Maximum desired element size of the mesh.
+#     scale_function : str
+#         Formula for the function describing the desired mesh element sizes.
+#     """
+#     # Remove extension from meshName
+#     meshName, extension = os.path.splitext(meshName)
+
+#     # create the mesh generating gmsh script file
+#     meshing_script = f'{meshName}.geo'
+#     resulting_mesh = f'{meshName}.vtk'
+
+#     # Get the path to this folder
+#     thisFolder = os.path.dirname(os.path.abspath(__file__))
+
+#     # Extract the template
+#     with open(f'{thisFolder}/templates/mesh_from_function.template', 'r') as file:
+#         template = Template(file.read())
+
+#     # Write the gmsh script
+#     with open(meshing_script, 'w') as file:
+#         file.write(template.substitute(
+#             BOUNDARY       = boundary,
+#             SCALE_MIN      = scale_min,
+#             SCALE_MAX      = scale_max,
+#             SCALE_FUNCTION = scale_function ))
+
+#     # Run gmsh in a subprocess to generate the mesh from the background
+#     run(f'gmsh {meshing_script} -3 -saveall -o {resulting_mesh}')
+
+#     # Remove the script file
+#     os.remove(meshing_script)
+
+#     return
+
+
+# def generate_background_from_1D_data(meshName, R, data):
+#     """
+#     Generate a background mesh from 1D (radial) data.
+
+#     Parameters
+#     ----------
+#     meshName : str
+#         File name of the mesh file.
+#     R : array_type
+#         Array containing the radial positions at which the data is given.
+#     data : array_type
+#         Array containing the data to put a background.
+#     """
+#     # Remove extension from meshName
+#     meshName, extension = os.path.splitext(meshName)
+
+#     # Map the data to a set of concentric HEALPix spheres
+#     nsides = 8
+#     sphere = np.array(pixelfunc.pix2vec(nsides, range(12*nsides**2))).transpose()
+#     points = []
+#     data_s = []
+#     for (i,r) in enumerate(R):
+#         for s in Rotation.random().apply(sphere):
+#             points.append(r*s)
+#             data_s.append(data[i])
+
+#     # Delaunay tetrahedralise the point set
+#     delaunay = Delaunay(points)
+
+#     # Write the mesh
+#     meshio.write_points_cells(
+#         filename   = f'{meshName}.msh',
+#         points     =             delaunay.points,
+#         cells      = {'tetra'  : delaunay.simplices},
+#         point_data = {'weights': np.array(data_s)}   )
+
+#     return
+
+
+# @njit
+# def get_pfs(nns):
+#     pfs = np.zeros(nns.shape, dtype=np.int64)
+#     tot = 0
+#     for i, n in enumerate(nns):
+#         tot    += n
+#         pfs[i]  = tot
+#     return pfs
+
+
+# @njit
+# def get_nbs(nbs, pfs, i):
+#     if (i == 0):
+#         return nbs[0:pfs[0]]
+#     else:
+#         return nbs[pfs[i-1]:pfs[i]]
+
+
+# @njit
+# def norm(vec_arr):
+#     norms = np.zeros(vec_arr.shape[0])
+#     for i in range(vec_arr.shape[0]):
+#         for j in range(vec_arr.shape[1]):
+#             norms[i] = norms[i] + vec_arr[i,j]**2
+#     return np.sqrt(norms)
+
+
+# @njit
+# def get_Gs(tracer, nbs, pfs):
+#     """
+#     Get the maximum relative difference in the tracer function for each point.
+#     """
+#     Gs = np.zeros(tracer.shape[0])
+#     for i in range(tracer.shape[0]):
+#         tracer_i = tracer[i]
+#         tracer_n = tracer[get_nbs(nbs, pfs, i)]
+#         Gs[i] = np.max(np.abs(tracer_i - tracer_n) / (tracer_i + tracer_n))
+#     return Gs
+
+
+# @njit
+# def get_Ls(pos, nbs, pfs):
+#     """
+#     Get the element size distributions for each point.
+#     """
+#     Ls = np.zeros(pos.shape[0])
+#     for i in range(pos.shape[0]):
+#         Ls[i] = np.mean(norm(pos[i] - pos[get_nbs(nbs, pfs, i)]))
+#     return Ls
+
+
+# @njit
+# def get_weights(pos, nbs, nns, tracer, threshold=0.21, fmin=1.0, ftarget=2.15):
+#     """
+#     Get the weights (desired element sizes) for each point.
+#     """
+#     pfs = get_pfs (nns)
+#     GGs = get_Gs  (tracer, nbs, pfs)
+#     LLs = get_Ls  (pos,    nbs, pfs)
+
+#     L_min    = fmin    * LLs
+#     L_target = ftarget * LLs
+
+#     weights = L_target
+#     weights[GGs > threshold] = L_min[GGs > threshold]
+
+#     return weights
+
+
+# def reduce(
+#         meshName,
+#         delaunay,
+#         tracer,
+#         boundary,
+#         scale_max = 1.0e+99,  # Maximum value of the scale parameter
+#         scale_min = 0.0e+00,  # Minimum value of the scale parameter
+#         threshold = 0.21,     # Threashold relative difference for coarsening
+#         fmin      = 1.0,      # Don't allow refinenment
+#         ftarget   = 2.15      # Approx 10^(-1/3) for approx 10 times fewer points
+#     ):
+#     """
+#     Reduce a model mesh for a given tracer function.
+
+#     Parameters
+#     ----------
+#     meshName : str
+#         File name of the mesh file.
+#     delaunay :
+#         scipy.spatial.Delaunay object for the given mesh.
+#     tracer : array_like
+#         1D array containing the tabulated tracer function.
+#     boundary : str
+#         Gmsh script containing the definition of the boundary.
+#     scale_max : float
+#         Maximum desired element size of the mesh.
+#     scale_min : float
+#         Minimum desired element size of the mesh.
+#     threshold : float
+#         Threashold value for relative difference below which we coarsen.
+#     fmin : float
+#         Multiplication factor of the mesh element size where the relative difference is below the threashold.
+#     ftarget : float
+#         Multiplication factor of the mesh element size where the relative difference is below the threashold.
+#     """
+#     # Remove extension from meshName
+#     meshName, extension = os.path.splitext(meshName)
+
+#     # Extract mesh structure from delaunay object
+#     (indptr, indices) = delaunay.vertex_neighbor_vertices
+#     neighbors = [indices[indptr[k]:indptr[k+1]] for k in range(len(delaunay.points))]
+#     nbs       = np.array([n for sublist in neighbors for n in sublist])
+#     nns       = np.array([len(sublist) for sublist in neighbors])
+#     pos       = delaunay.points
+
+#     # Create a background mesh (in .msh format)
+#     meshio.write_points_cells(
+#         filename   = f'{meshName}.msh',
+#         points     = delaunay.points,
+#         cells      = {'tetra'  : delaunay.simplices},
+#         point_data = {'weights': get_weights(pos=pos,
+#                                              nbs=nbs,
+#                                              nns=nns,
+#                                              tracer=tracer,
+#                                              threshold=threshold,
+#                                              fmin=fmin,
+#                                              ftarget=ftarget)}
+#     )
+
+#     # Convert .msh to .pos mesh for Gmsh
+#     convert_msh_to_pos (meshName=meshName, replace=True)
+
+#     # Create a new mesh from the background mesh
+#     create_mesh_from_background(
+#         meshName  = meshName,
+#         boundary  = boundary,
+#         scale_min = scale_min,
+#         scale_max = scale_max
+#     )
+
+#     return
 
 def remesh_point_cloud(positions, data, max_depth=9, threshold= 5e-2, hullorder = 3):
     '''
diff --git a/setup.py b/setup.py
index 3e5188fb..ec6b0eb3 100644
--- a/setup.py
+++ b/setup.py
@@ -50,7 +50,6 @@ def build_extension(self, ext):
         'scipy',
         'astropy',
         'numba',
-        'meshio',
         'healpy',
         'yt',
         'mpi4py',
diff --git a/tests/benchmarks/numeric/reduce_phantom_3D_model.py b/tests/benchmarks/numeric/reduce_phantom_3D_model.py
index 8b0f475b..becdb067 100644
--- a/tests/benchmarks/numeric/reduce_phantom_3D_model.py
+++ b/tests/benchmarks/numeric/reduce_phantom_3D_model.py
@@ -2,7 +2,6 @@
 import magritte.core   as magritte                     # Core functionality
 import magritte.mesher as mesher                       # Mesher
 import numpy           as np                           # Data structures
-import vtk                                             # Reading the model
 import warnings                                        # Hide warnings
 warnings.filterwarnings('ignore')                      # especially for yt
 import yt                                              # 3D plotting
@@ -10,7 +9,6 @@
 
 from tqdm                   import tqdm                # Progress bars
 from astropy                import constants           # Unit conversions
-from vtk.util.numpy_support import vtk_to_numpy        # Converting data
 from scipy.spatial          import Delaunay, cKDTree   # Finding neighbors
 from yt.funcs               import mylog               # To avoid yt output
 mylog.setLevel(40)                                     # as error messages
diff --git a/tests/benchmarks/numeric/vanZadelhoff_1_3D_mesher.py b/tests/benchmarks/numeric/vanZadelhoff_1_3D_mesher.py
index 18ab8b46..69e4edfa 100644
--- a/tests/benchmarks/numeric/vanZadelhoff_1_3D_mesher.py
+++ b/tests/benchmarks/numeric/vanZadelhoff_1_3D_mesher.py
@@ -35,29 +35,19 @@
     'b' : 1.0E-6
 }
 
-scale_max = 0.05 * r_out
-scale_min = 0.05 * r_in
-scale_cte = 0.05 * r_in
-scale_fun = f'{scale_cte / r_in**2} * (x*x + y*y + z*z)'
+#generate log grid for remeshing
+posoneDvals = np.logspace(np.log10(r_in), np.log10(r_out), 100)
+oneDvals = np.concatenate((-posoneDvals[::-1], posoneDvals))
+XYZ = np.meshgrid(oneDvals, oneDvals, oneDvals)
+positions = np.array([XYZ[0].flatten(), XYZ[1].flatten(), XYZ[2].flatten()]).T
 
-meshName = f'{moddir}/vanZadelhoff_1_3D_mesher.vtk'
-
-mesher.create_mesh_from_function(
-    meshName       = meshName,
-    boundary       = mesher.boundary_sphere_in_sphere(
-                         radius_in  = r_in,
-                         radius_out = r_out),
-    scale_min      = scale_min,
-    scale_max      = scale_max,
-    scale_function = scale_fun )
-
-mesh = mesher.Mesh(meshName)
-
-npoints = len(mesh.points)
-# nbs     = [n for sublist in mesh.neighbors for n in sublist]
-# n_nbs   = [len(sublist) for sublist in mesh.neighbors]
-
-rs = np.linalg.norm(mesh.points, axis=1)
+#remeshed the log grid #fairly low resolution
+r = np.sqrt(positions[:,0]**2 + positions[:,1]**2 + positions[:,2]**2)
+remesh_fun = 1/(r**2)#approximately how the density drops
+remesh_points, nbound = mesher.remesh_point_cloud(positions, remesh_fun, max_depth=12, threshold=4e-1, hullorder=3)
+remesh_points, nbound =mesher.point_cloud_add_spherical_outer_boundary(remesh_points, nbound, r_out)
+npoints = len(remesh_points)
+rs = np.linalg.norm(remesh_points, axis=1)
 
 
 def create_model (a_or_b):
@@ -87,7 +77,7 @@ def nTT (r):
     model.parameters.set_nlspecs           (nlspecs)
     model.parameters.set_nquads            (nquads)
 
-    model.geometry.points.position.set(mesh.points)
+    model.geometry.points.position.set(remesh_points)
     model.geometry.points.velocity.set(np.zeros((npoints, 3)))
 
 #     model.geometry.points.  neighbors.set(  nbs)
@@ -101,8 +91,8 @@ def nTT (r):
 
     model = setup.set_Delaunay_neighbor_lists (model)
 
-    model.parameters.set_nboundary(len(mesh.boundary))
-    model.geometry.boundary.boundary2point.set(mesh.boundary)
+    model.parameters.set_nboundary(nbound)
+    model.geometry.boundary.boundary2point.set(np.arange(nbound))
 
     model = setup.set_boundary_condition_CMB  (model)
     model = setup.set_uniform_rays            (model)
diff --git a/tests/benchmarks/numeric/vanZadelhoff_2_3D_mesher.py b/tests/benchmarks/numeric/vanZadelhoff_2_3D_mesher.py
index 2fa0f234..f2ab977e 100644
--- a/tests/benchmarks/numeric/vanZadelhoff_2_3D_mesher.py
+++ b/tests/benchmarks/numeric/vanZadelhoff_2_3D_mesher.py
@@ -33,7 +33,7 @@ def create_model (a_or_b):
     # Parameters
     dimension = 3
     nshells   = len(r_shell)
-    nrays     = 12 * 3**2
+    nrays     = 12 * 2**2
     nspecs    = 3
     nlspecs   = 1
     nquads    = 7
@@ -61,34 +61,21 @@ def create_model (a_or_b):
     r_in  = min(r_shell)
     r_out = max(r_shell)
 
-    scale_max = 0.1 * r_out
-    scale_min = 0.1 * r_in
-    scale_cte = 0.1 * r_in
-    scale_fun = f'{scale_cte / r_in**2} * (x*x + y*y + z*z)'
+    #generate log grid for remeshing
+    posoneDvals = np.logspace(np.log10(r_in/2), np.log10(r_out), 100)#slightly smaller r to allow for better refinement
+    oneDvals = np.concatenate((-posoneDvals[::-1], posoneDvals))
+    XYZ = np.meshgrid(oneDvals, oneDvals, oneDvals)
+    positions = np.array([XYZ[0].flatten(), XYZ[1].flatten(), XYZ[2].flatten()]).T
 
-    meshName = f'{moddir}/vanZadelhoff_2_3D_mesher.vtk'
+    #remeshed the log grid #fairly low resolution
+    r = np.sqrt(positions[:,0]**2 + positions[:,1]**2 + positions[:,2]**2)
+    remesh_fun = nH2_int(r)#approximately how the density drops
+    remesh_points, nbound = mesher.remesh_point_cloud(positions, remesh_fun, max_depth=10, threshold=3.5e-1, hullorder=3)
+    remesh_points, nbound = mesher.point_cloud_add_spherical_outer_boundary(remesh_points, nbound, r_out)
+    npoints = len(remesh_points)
+    rs = np.linalg.norm(remesh_points, axis=1)
 
-    mesher.create_mesh_from_function(
-        meshName       = meshName,
-        boundary       = mesher.boundary_sphere_in_sphere(
-                             radius_in  = r_in,
-                             radius_out = r_out),
-        scale_min      = scale_min,
-        scale_max      = scale_max,
-        scale_function = scale_fun )
-
-    mesh = mesher.Mesh(meshName)
-
-    position = mesh.points
-
-    npoints = len(mesh.points)
-    print(npoints)
-#     nbs     = [n for sublist in mesh.neighbors for n in sublist]
-#     n_nbs   = [len(sublist) for sublist in mesh.neighbors]
-
-    rs = np.linalg.norm(mesh.points, axis=1)
-
-    velocity = (vs_int(rs) / rs * position.T).T
+    velocity = (vs_int(rs) / rs * remesh_points.T).T
 
     model = magritte.Model ()
     model.parameters.set_spherical_symmetry(False)
@@ -100,7 +87,7 @@ def create_model (a_or_b):
     model.parameters.set_nlspecs           (nlspecs)
     model.parameters.set_nquads            (nquads)
 
-    model.geometry.points.position.set(position)
+    model.geometry.points.position.set(remesh_points)
     model.geometry.points.velocity.set(velocity)
 
 #     model.geometry.points.  neighbors.set(  nbs)
@@ -114,8 +101,8 @@ def create_model (a_or_b):
 
     model = setup.set_Delaunay_neighbor_lists (model)
 
-    model.parameters.set_nboundary(len(mesh.boundary))
-    model.geometry.boundary.boundary2point.set(mesh.boundary)
+    model.parameters.set_nboundary(nbound)
+    model.geometry.boundary.boundary2point.set(np.arange(nbound))
 
     model = setup.set_boundary_condition_CMB  (model)
     model = setup.set_uniform_rays            (model)
@@ -152,7 +139,7 @@ def run_model (a_or_b, nosave=False):
     model.compute_level_populations (True, 100)
     timer3.stop()
 
-    pops = np.array(model.lines.lineProducingSpecies[0].population).reshape((model.parameters.npoints(), 2))
+    pops = np.array(model.lines.lineProducingSpecies[0].population).reshape((model.parameters.npoints(), 21))
     abun = np.array(model.chemistry.species.abundance)[:,0]
     rs   = np.linalg.norm(np.array(model.geometry.points.position), axis=1)