diff --git a/README.rst b/README.rst
index d24f650..7d05101 100644
--- a/README.rst
+++ b/README.rst
@@ -123,6 +123,113 @@ Here is the cell quality as computed according to the minimum scaled jacobian.
 .. image:: https://github.com/pyvista/tetgen/raw/main/doc/images/sphere_qual.png
 
 
+Using a Background Mesh
+-----------------------
+A background mesh in TetGen is used to define a mesh sizing function for
+adaptive mesh refinement. This function informs TetGen of the desired element
+size throughout the domain, allowing for detailed refinement in specific areas
+without unnecessary densification of the entire mesh. Here's how to utilize a
+background mesh in your TetGen workflow:
+
+1. **Generate the Background Mesh**: Create a tetrahedral mesh that spans the
+   entirety of your input piecewise linear complex (PLC) domain. This mesh will
+   serve as the basis for your sizing function.
+
+2. **Define the Sizing Function**: At the nodes of your background mesh, define
+   the desired mesh sizes. This can be based on geometric features, proximity
+   to areas of interest, or any criterion relevant to your simulation needs.
+
+3. **Optional: Export the Background Mesh and Sizing Function**: Save your
+   background mesh in the TetGen-readable `.node` and `.ele` formats, and the
+   sizing function values in a `.mtr` file. These files will be used by TetGen
+   to guide the mesh generation process.
+
+4. **Run TetGen with the Background Mesh**: Invoke TetGen, specifying the
+   background mesh. TetGen will adjust the mesh according to the provided
+   sizing function, refining the mesh where smaller elements are desired.
+
+**Full Example**
+
+To illustrate, consider a scenario where you want to refine a mesh around a
+specific region with increased detail. The following steps and code snippets
+demonstrate how to accomplish this with TetGen and PyVista:
+
+1. **Prepare Your PLC and Background Mesh**:
+
+   .. code-block:: python
+
+      import pyvista as pv
+      import tetgen
+      import numpy as np
+
+      # Load or create your PLC
+      sphere = pv.Sphere(theta_resolution=10, phi_resolution=10)
+
+      # Generate a background mesh with desired resolution
+      def generate_background_mesh(bounds, resolution=20, eps=1e-6):
+          x_min, x_max, y_min, y_max, z_min, z_max = bounds
+          grid_x, grid_y, grid_z = np.meshgrid(
+              np.linspace(xmin - eps, xmax + eps, resolution),
+              np.linspace(ymin - eps, ymax + eps, resolution),
+              np.linspace(zmin - eps, zmax + eps, resolution),
+              indexing="ij",
+          )
+          return pv.StructuredGrid(grid_x, grid_y, grid_z).triangulate()
+
+      bg_mesh = generate_background_mesh(sphere.bounds)
+
+2. **Define the Sizing Function and Write to Disk**:
+
+   .. code-block:: python
+
+      # Define sizing function based on proximity to a point of interest
+      def sizing_function(points, focus_point=np.array([0, 0, 0]), max_size=1.0, min_size=0.1):
+          distances = np.linalg.norm(points - focus_point, axis=1)
+          return np.clip(max_size - distances, min_size, max_size)
+
+      bg_mesh.point_data['target_size'] = sizing_function(bg_mesh.points)
+
+      # Optionally write out the background mesh
+      def write_background_mesh(background_mesh, out_stem):
+          """Write a background mesh to a file.
+
+          This writes the mesh in tetgen format (X.b.node, X.b.ele) and a X.b.mtr file
+          containing the target size for each node in the background mesh.
+          """
+          mtr_content = [f"{background_mesh.n_points} 1"]
+          target_size = background_mesh.point_data["target_size"]
+          for i in range(background_mesh.n_points):
+              mtr_content.append(f"{target_size[i]:.8f}")
+
+      write_background_mesh(bg_mesh, 'bgmesh.b')
+
+3. **Use TetGen with the Background Mesh**:
+
+
+   Directly pass the background mesh from PyVista to ``tetgen``:
+
+   .. code-block:: python
+
+      tet_kwargs = dict(order=1, mindihedral=20, minratio=1.5)
+      tet = tetgen.TetGen(mesh)
+      tet.tetrahedralize(bgmesh=bgmesh, **tet_kwargs)
+      refined_mesh = tet.grid
+
+   Alternatively, use the background mesh files.
+
+   .. code-block:: python
+
+      tet = tetgen.TetGen(sphere)
+      tet.tetrahedralize(bgmeshfilename='bgmesh.b', **tet_kwargs)
+      refined_mesh = tet.grid
+
+
+This example demonstrates generating a background mesh, defining a spatially
+varying sizing function, and using this background mesh to guide TetGen in
+refining a PLC. By following these steps, you can achieve adaptive mesh
+refinement tailored to your specific simulation requirements.
+
+
 Acknowledgments
 ---------------
 Software was originally created by Hang Si based on work published in
diff --git a/pyproject.toml b/pyproject.toml
index 707e1f8..cf9acea 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -19,12 +19,12 @@ filterwarnings = [
 [tool.cibuildwheel]
 archs = ["auto64"]  # 64-bit only
 skip = "pp* *musllinux* cp37-*"  # disable PyPy, musl-based wheels, and Python < 3.8
-test-requires = "pytest pymeshfix"
+before-test = "pip install -r requirements_test.txt"
 test-command = "pytest {project}/tests"
 
 [tool.cibuildwheel.macos]
 # https://cibuildwheel.readthedocs.io/en/stable/faq/#apple-silicon
-archs = ["x86_64", "universal2"]
+archs = ["universal2"]
 test-skip = ["*_arm64", "*_universal2:arm64"]
 
 [tool.codespell]
diff --git a/requirements_test.txt b/requirements_test.txt
index 6e3302b..4999a79 100644
--- a/requirements_test.txt
+++ b/requirements_test.txt
@@ -1,3 +1,4 @@
+meshio
+pymeshfix>=0.13.4
 pytest
 pytest-cov
-pymeshfix>=0.13.4
diff --git a/tests/test_background_mesh.py b/tests/test_background_mesh.py
new file mode 100644
index 0000000..8823809
--- /dev/null
+++ b/tests/test_background_mesh.py
@@ -0,0 +1,121 @@
+"""Test the mesh resizing feature of tetgen with sizing function."""
+from pathlib import Path
+import tempfile
+
+import numpy as np
+import pyvista as pv
+import tetgen
+
+
+def sizing_function(points):
+    """Return the target size at a given point.
+
+    This is just an example. You can use any function you want.
+    """
+    x, y, z = points.T
+    return np.where(x < 0, 0.5, 0.1)
+
+
+def generate_background_mesh(boundary_mesh, resolution=20, eps=1e-6):
+    """Generate a new mesh with the same bounds as the boundary meshy.
+
+    We will use this as a background mesh for the sizing function.
+    """
+    xmin, xmax, ymin, ymax, zmin, zmax = boundary_mesh.bounds
+    new_vertices = np.meshgrid(
+        np.linspace(xmin - eps, xmax + eps, resolution),
+        np.linspace(ymin - eps, ymax + eps, resolution),
+        np.linspace(zmin - eps, zmax + eps, resolution),
+        indexing="ij",
+    )
+
+    # tetgen supports only tetrahedral meshes
+    return pv.StructuredGrid(*new_vertices).triangulate()
+
+
+def write_background_mesh(background_mesh, out_stem):
+    """Write a background mesh to a file.
+
+    This writes the mesh in tetgen format (X.b.node, X.b.ele) and a X.b.mtr file
+    containing the target size for each node in the background mesh.
+    """
+    mtr_content = [f"{background_mesh.n_points} 1"]
+    target_size = background_mesh.point_data["target_size"]
+    for i in range(background_mesh.n_points):
+        mtr_content.append(f"{target_size[i]:.8f}")
+
+    pv.save_meshio(f"{out_stem}.node", background_mesh)
+    mtr_file = f"{out_stem}.mtr"
+
+    with open(mtr_file, "w") as f:
+        f.write("\n".join(mtr_content))
+
+
+def mesh_resizing_with_bgmeshfilename(mesh, bgmesh, **kwargs):
+    """Performs mesh resizing with a background mesh file."""
+    with tempfile.TemporaryDirectory() as tmpdir:
+        tmpfile = Path(tmpdir).joinpath("bgmesh.b")
+        write_background_mesh(bgmesh, tmpfile)
+
+        # Pass the background mesh file to tetgen
+        tet = tetgen.TetGen(mesh)
+        tet.tetrahedralize(bgmeshfilename=str(tmpfile), metric=1, **kwargs)
+        grid = tet.grid
+
+    # extract cells below the 0 xy plane
+    # cell_center = grid.cell_centers().points
+    # subgrid = grid.extract_cells(cell_center[:, 2] < 0)
+
+    # debug: plot this
+    # plotter = pv.Plotter()
+    # plotter.set_background("w")
+    # plotter.add_mesh(subgrid, "lightgrey", lighting=True)
+    # plotter.add_mesh(grid, "r", "wireframe")
+    # plotter.add_legend([[" Input Mesh ", "r"], [" Tessellated Mesh ", "black"]])
+    # plotter.show() # Uncomment for visualisation of resized mesh
+
+    return grid
+
+
+def mesh_resizing_with_pyvista_bgmesh(mesh, bgmesh, **kwargs):
+    """Performs mesh resizing with a pyvista bgmesh."""
+    # Pass the background mesh to tetgen
+    tet = tetgen.TetGen(mesh)
+    tet.tetrahedralize(bgmesh=bgmesh, metric=1, **kwargs)
+    grid = tet.grid
+
+    # Extract cells below the 0 xy plane
+    # cell_center = grid.cell_centers().points
+    # subgrid = grid.extract_cells(cell_center[:, 2] < 0)
+
+    # Debug: uncomment for visualisation of resized mesh
+    # plotter = pv.Plotter()
+    # plotter.set_background("w")
+    # plotter.add_mesh(subgrid, "lightgrey", lighting=True)
+    # plotter.add_mesh(grid, "r", "wireframe")
+    # plotter.add_legend([[" Input Mesh ", "r"], [" Tessellated Mesh ", "black"]])
+    # plotter.show()
+    return grid
+
+
+def test_mesh_resizing():
+    """Test the mesh resizing feature of tetgen with sizing function."""
+    sphere = pv.Sphere(theta_resolution=10, phi_resolution=10)
+    tet_kwargs = dict(order=1, mindihedral=20, minratio=1.5)
+
+    # Vanilla tetgen for reference
+    tet = tetgen.TetGen(sphere)
+    tet.tetrahedralize(**tet_kwargs)
+    grid = tet.grid
+
+    # Generate background mesh
+    bgmesh = generate_background_mesh(sphere)
+    bgmesh.point_data["target_size"] = sizing_function(bgmesh.points)
+
+    resized_grid_file = mesh_resizing_with_bgmeshfilename(sphere, bgmesh, **tet_kwargs)
+    assert resized_grid_file.n_points >= grid.n_points
+
+    resized_grid_direct = mesh_resizing_with_pyvista_bgmesh(sphere, bgmesh, **tet_kwargs)
+    assert resized_grid_direct.n_points > grid.n_points
+
+    assert resized_grid_file == resized_grid_direct
diff --git a/tetgen/cython/tetgen/_tetgen.pyx b/tetgen/cython/tetgen/_tetgen.pyx
index d7d12c7..fb9c66d 100644
--- a/tetgen/cython/tetgen/_tetgen.pyx
+++ b/tetgen/cython/tetgen/_tetgen.pyx
@@ -5,16 +5,18 @@
 # cython nonecheck=False
 
 import numpy as np
+
 cimport numpy as np
 
 import ctypes
+
 from cython cimport view
+from libc.string cimport strcpy
 
 # # Numpy must be initialized. When using numpy from C or Cython you must
 # # _always_ do that, or you will have segfaults
 # np.import_array()
 
-
 """ Wrapped tetgen class """
 cdef extern from "tetgen_wrap.h":
     cdef cppclass tetgenio_wrap:
@@ -31,17 +33,21 @@ cdef extern from "tetgen_wrap.h":
 
         # Loads Arrays directly to tetgenio object
         void LoadArray(int, double*, int, int*)
-        
+        # Loads MTR Arrays directly to tetgenio object
+        void LoadMTRArray(int, double*, int, int*, double*)
+        # Loads tetmesh from file
+        bint LoadTetMesh(char*, int)
 
-cdef extern from "tetgen.h":
 
+cdef extern from "tetgen.h":
     cdef cppclass tetrahedralize:
-        int tetrahedralize(char*, tetgenio_wrap*, tetgenio_wrap*)
+        int tetrahedralize(char*, tetgenio_wrap*, tetgenio_wrap*, tetgenio_wrap*, tetgenio_wrap*)
+
 
     cdef cppclass tetgenbehavior:
         void tetgenbehavior()
 
-        # Switches of TetGen. 
+        # Switches of TetGen.
         int plc
         int psc
         int refine
@@ -125,10 +131,15 @@ cdef extern from "tetgen.h":
         double coarsen_percent
         double elem_growth_ratio
         double refine_progress_ratio
+        char bgmeshfilename[1024]
 
     # Different calls depending on using settings input
-    cdef void tetrahedralize(tetgenbehavior*, tetgenio_wrap*, tetgenio_wrap*) except +
-    cdef void tetrahedralize(char*, tetgenio_wrap*, tetgenio_wrap*) except +
+    cdef void tetrahedralize(tetgenbehavior*, tetgenio_wrap*, tetgenio_wrap*, tetgenio_wrap*, tetgenio_wrap*) except +
+    cdef void tetrahedralize(char*, tetgenio_wrap*, tetgenio_wrap*, tetgenio_wrap*, tetgenio_wrap*) except +
+
+
+cdef extern from "tetgen.h" namespace "tetgenbehavior":
+    cdef enum objecttype: NODES, POLY, OFF, PLY, STL, MEDIT, VTK, MESH, NEU_MESH
 
 
 cdef class PyBehavior:
@@ -147,23 +158,23 @@ cdef class PyTetgenio:
 
 
     def ReturnNodes(self):
-        """ Returns nodes from tetgen """ 
-        
+        """ Returns nodes from tetgen """
+
         # Create python copy of array
         cdef int npoints  = self.c_tetio.numberofpoints*3
         cdef double [::1] nodes = np.empty(npoints)
-        
+
         cdef int i
         for i in range(npoints):
             nodes[i] = self.c_tetio.pointlist[i]
 
         return np.asarray(nodes).reshape((-1, 3))
-        
-        
+
+
     def ReturnTetrahedrals(self, order):
         """ Returns tetrahedrals from tetgen """
 
-        # Determine         
+        # Determine
         if order == 1:
             arrsz = self.c_tetio.numberoftetrahedra*4
         else:
@@ -171,7 +182,7 @@ cdef class PyTetgenio:
 
         # Create python copy of tetrahedral array
         cdef int [::1] tets = np.empty(arrsz, ctypes.c_int)
-        
+
         cdef int i
         cdef int arrsz_c = arrsz
         for i in range(arrsz_c):
@@ -188,7 +199,7 @@ cdef class PyTetgenio:
             """
             vtkQuadraticTetra
             The ordering of the ten points defining the cell is point ids
-            (0-3, 4-9) where ids 0-3 are the four tetra vertices, 
+            (0-3, 4-9) where ids 0-3 are the four tetra vertices,
             and point ids 4-9 are the midedge nodes between:
             (0,1), (1,2), (2,0), (0,3), (1,3), and (2,3)
             """
@@ -198,17 +209,29 @@ cdef class PyTetgenio:
             tetcopy[:, 7] = tetarr[:, 5]
             tetcopy[:, 8] = tetarr[:, 8]
             tetcopy[:, 9] = tetarr[:, 4]
-            
+
             return tetcopy
-            
-        
+
+
     def LoadMesh(self, double [::1] points, int [::1] faces):
         """ Loads points and faces into TetGen """
         npoints = points.size/3
         nfaces = faces.size/3
         self.c_tetio.LoadArray(npoints, &points[0], nfaces, &faces[0])
-    
-    
+
+
+    def LoadMTRMesh(self, double [::1] points, int [::1] tets, double [::1] mtr):
+        """ Loads points and tets into TetGen """
+        npoints = points.size/3
+        ntets = tets.size/4
+        self.c_tetio.LoadMTRArray(npoints, &points[0], ntets, &tets[0], &mtr[0])
+
+
+    def LoadTetMesh(self, char* filename, int order):
+        """ Loads tetmesh from file """
+        return self.c_tetio.LoadTetMesh(filename, order)
+
+
 def Tetrahedralize(
         v,
         f,
@@ -298,44 +321,57 @@ def Tetrahedralize(
         coarsen_percent=1.0,
         elem_growth_ratio=0.0,
         refine_progress_ratio=0.333,
+        bgmeshfilename='',
+        bgmesh_v=None,
+        bgmesh_tet=None,
+        bgmesh_mtr=None,
     ):
 
     """Tetgen function to interface with TetGen C++ program."""
     # convert switches to c object
     cdef char *cstring = switches
 
+    # convert bgmeshfilename to c object
+    cdef char bgmeshfilename_c[1024]
+    cdef bytes py_bgmeshfilename = bgmeshfilename.encode('utf-8')
+    cdef char* bgmeshfilename_py = py_bgmeshfilename
+    strcpy(bgmeshfilename_c, bgmeshfilename_py)
+
     # Check that inputs are valid
-    if not v.flags['C_CONTIGUOUS']:
-        if v.dtype != np.float64:
-            v = np.ascontiguousarray(v, dtype=np.float64)
-        else:
-            v = np.ascontiguousarray(v)
+    def cast_to_cint(x):
+        return np.ascontiguousarray(x, dtype=ctypes.c_int)
 
-    elif v.dtype != np.float64:
-        v = v.astype(np.float64)
+    def cast_to_cdouble(x):
+        return np.ascontiguousarray(x, dtype=np.float64)
 
-    # Ensure inputs are of the right type
-    if not f.flags['C_CONTIGUOUS']:
-        if f.dtype != ctypes.c_int:
-            f = np.ascontiguousarray(f, dtype=ctypes.c_int)
-        else:
-            f = np.ascontiguousarray(f)
+    v = cast_to_cdouble(v)
+    f = cast_to_cint(f)
+
+    if bgmesh_v is not None:
+        bgmesh_v = cast_to_cdouble(bgmesh_v)
+    if bgmesh_tet is not None:
+        bgmesh_tet = cast_to_cint(bgmesh_tet)
+    if bgmesh_mtr is not None:
+        bgmesh_mtr = cast_to_cdouble(bgmesh_mtr)
 
-    elif f.dtype != ctypes.c_int:
-        f = f.astype(ctypes.c_int)
-    
     # Create input class
     tetgenio_in = PyTetgenio()
     tetgenio_in.LoadMesh(v.ravel(), f.ravel())
-        
+
     # Create output class
-    tetgenio_out = PyTetgenio()        
-    
+    tetgenio_out = PyTetgenio()
+
+    tetgenio_bg = PyTetgenio()
+    if bgmesh_mtr is not None:
+        tetgenio_bg.LoadMTRMesh(bgmesh_v, bgmesh_tet, bgmesh_mtr)
+    if bgmeshfilename:
+        tetgenio_bg.LoadTetMesh(bgmeshfilename_c, <int>objecttype.NODES)
+
     if switches:
-        tetrahedralize(cstring, &tetgenio_in.c_tetio, &tetgenio_out.c_tetio)
+        tetrahedralize(cstring, &tetgenio_in.c_tetio, &tetgenio_out.c_tetio, NULL, &tetgenio_bg.c_tetio)
         if b'o2' in switches:
             order = 2
-    
+
     else: # set defaults or user input settings
         # Enable plc if checking for self intersections
         if diagnose:
@@ -345,7 +381,7 @@ def Tetrahedralize(
         if order != 1 and order != 2:
             raise Exception('Settings error: Order should be 1 or 2')
 
-        # Set behavior 
+        # Set behavior
         behavior = PyBehavior()
         behavior.c_behavior.plc = plc
         behavior.c_behavior.psc = psc
@@ -428,14 +464,14 @@ def Tetrahedralize(
         behavior.c_behavior.coarsen_percent = coarsen_percent
         behavior.c_behavior.elem_growth_ratio = elem_growth_ratio
         behavior.c_behavior.refine_progress_ratio = refine_progress_ratio
-    
+
         # Process from C++ side using behavior object
-        tetrahedralize(&behavior.c_behavior, &tetgenio_in.c_tetio, 
-                       &tetgenio_out.c_tetio)
+        tetrahedralize(&behavior.c_behavior, &tetgenio_in.c_tetio,
+                       &tetgenio_out.c_tetio, NULL, &tetgenio_bg.c_tetio)
+
 
     # Returns vertices and tetrahedrals of new mesh
     nodes = tetgenio_out.ReturnNodes()
     tets = tetgenio_out.ReturnTetrahedrals(order)
 
-    return nodes, tets  
-
+    return nodes, tets
diff --git a/tetgen/cython/tetgen/tetgen_wrap.cxx b/tetgen/cython/tetgen/tetgen_wrap.cxx
index 558e2bc..4c67cbb 100644
--- a/tetgen/cython/tetgen/tetgen_wrap.cxx
+++ b/tetgen/cython/tetgen/tetgen_wrap.cxx
@@ -30,7 +30,7 @@ void tetgenio_wrap::LoadArray(int npoints, double* points, int nfaces,
     // Initialize a face
     f = &facetlist[i];
     init(f);
-    
+
     // Each facet has one polygon, no hole, and each polygon has a three\
     //vertices
     f->numberofpolygons = 1;
@@ -47,3 +47,39 @@ void tetgenio_wrap::LoadArray(int npoints, double* points, int nfaces,
   }
 }
 
+void tetgenio_wrap::LoadMTRArray(int npoints, double* points, int ntets,
+                                 int* tetarr, double* mtrpoints)
+{
+  int i, j;
+  int count = 0;
+
+  // Allocate memory for points and store them
+  numberofpoints = npoints;
+  pointlist = new double[npoints*3];
+
+  for(i = 0; i < npoints*3; i++) {
+    pointlist[i] = points[i];
+  }
+
+  // Populate pointmtrlist
+  numberofpointmtrs = 1;
+  pointmtrlist = new double[npoints];
+  for (i = 0; i < npoints; i++) {
+    pointmtrlist[i] = mtrpoints[i];
+  }
+
+  // Load tets (assumes 4 nodes per tetrahedron)
+  numberoftetrahedra = ntets;
+  numberofcorners = 4;
+  tetrahedronlist = new int[ntets*numberofcorners];
+  numberoftetrahedronattributes = 0;
+  for (i = 0; i < ntets*numberofcorners; i++) {
+    tetrahedronlist[i] = tetarr[i];
+  }
+}
+
+
+// Wrapper around loadtetmesh function for filename support
+bool tetgenio_wrap::LoadTetMesh(char* filename, int object) {
+    return load_tetmesh(filename, object);
+}
diff --git a/tetgen/cython/tetgen/tetgen_wrap.h b/tetgen/cython/tetgen/tetgen_wrap.h
index 7e0581a..fab6a26 100644
--- a/tetgen/cython/tetgen/tetgen_wrap.h
+++ b/tetgen/cython/tetgen/tetgen_wrap.h
@@ -11,6 +11,8 @@ class tetgenio_wrap : public tetgenio
         facet *f;
         polygon *p;
         void LoadArray(int, double*, int, int*);
+        void LoadMTRArray(int, double*, int, int*, double*);
+        bool LoadTetMesh(char*, int);
 
         //destructor
 //        ~myRectangle();
diff --git a/tetgen/pytetgen.py b/tetgen/pytetgen.py
index a0bf9e4..871f2f9 100644
--- a/tetgen/pytetgen.py
+++ b/tetgen/pytetgen.py
@@ -9,6 +9,7 @@
 LOG = logging.getLogger(__name__)
 LOG.setLevel("CRITICAL")
 
+MTR_POINTDATA_KEY = "target_size"
 
 invalid_input = TypeError(
     "Invalid input. Must be either a pyvista.PolyData object or vertex and face arrays."
@@ -301,6 +302,8 @@ def tetrahedralize(
         elem_growth_ratio=0.0,
         refine_progress_ratio=0.333,
         switches=None,
+        bgmeshfilename="",
+        bgmesh=None,
     ):
         """Generate tetrahedrals interior to the surface mesh.
 
@@ -381,6 +384,13 @@ def tetrahedralize(
             quadradic tetrahedrals.  Set order to 2 to output
             quadradic tetrahedrals.
 
+        bgmeshfilename : str, optional
+            Filename of the background mesh.
+
+        bgmesh : pv.UnstructuredGrid
+            Background mesh to be processed. Must be composed of only linear
+            tetra. Cannot specify both ``bgmeshfilename`` and ``bgmesh``.
+
         Returns
         -------
         nodes : numpy.ndarray
@@ -579,8 +589,21 @@ def tetrahedralize(
         if verbose == 0:
             quiet = 1
 
+        # Validate background mesh parameters
+        if bgmesh and bgmeshfilename:
+            raise ValueError("Cannot specify both bgmesh and bgmeshfilename")
+        if bgmesh or bgmeshfilename:
+            # Passing a background mesh only makes sense with metric set to true
+            # (will be silently ignored otherwise)
+            metric = 1
+        if bgmesh:
+            bgmesh_v, bgmesh_tet, bgmesh_mtr = self._process_bgmesh(bgmesh)
+        else:
+            bgmesh_v, bgmesh_tet, bgmesh_mtr = None, None, None
+
         # Call library
         plc = True  # always true
+
         try:
             self.node, self.elem = _tetgen.Tetrahedralize(
                 self.v,
@@ -667,6 +690,10 @@ def tetrahedralize(
                 coarsen_percent,
                 elem_growth_ratio,
                 refine_progress_ratio,
+                bgmeshfilename,
+                bgmesh_v,
+                bgmesh_tet,
+                bgmesh_mtr,
             )
         except RuntimeError as e:
             raise RuntimeError(
@@ -678,7 +705,7 @@ def tetrahedralize(
         # check if a mesh was generated
         if not np.any(self.node):
             raise RuntimeError(
-                "Failed to tetrahedralize.\n" + "May need to repair surface by making it manifold"
+                "Failed to tetrahedralize.\nMay need to repair surface by making it manifold"
             )
 
         # Return nodes and elements
@@ -747,3 +774,33 @@ def write(self, filename, binary=False):
         to convert to other formats in order to import into FEA software.
         """
         self.grid.save(filename, binary)
+
+    @staticmethod
+    def _process_bgmesh(mesh):
+        """Process a background mesh.
+
+        Parameters
+        ----------
+        bgmesh : pyvista.UnstructuredGrid
+            Background mesh to be processed. Must be composed of only linear tetra,
+
+        Returns
+        -------
+        bgmesh_v: numpy.ndarray
+            Vertex array of the background mesh.
+        bgmesh_tet: numpy.ndarray
+            Tet array of the background mesh.
+        bgmesh_mtr: numpy.ndarray
+            Target size array of the background mesh.
+        """
+        if MTR_POINTDATA_KEY not in mesh.point_data:
+            raise ValueError("Background mesh does not have target size information")
+
+        # Celltype check
+        if not (mesh.celltypes == pv.CellType.TETRA).all():
+            raise ValueError("`mesh` must contain only tetrahedrons")
+
+        bgmesh_v = mesh.points.astype(np.float64, copy=True).ravel()
+        bgmesh_tet = mesh.cell_connectivity.astype(np.int32, copy=True)
+        bgmesh_mtr = mesh.point_data[MTR_POINTDATA_KEY].astype(np.float64, copy=True).ravel()
+        return bgmesh_v, bgmesh_tet, bgmesh_mtr