diff --git a/src/fem/py.rs b/src/fem/py.rs
index 90f0edb..566362c 100644
--- a/src/fem/py.rs
+++ b/src/fem/py.rs
@@ -2,6 +2,7 @@ use super::{
     super::py::PyIntermediateError, finite_element_data_from_inp, write_finite_elements_metrics,
     write_finite_elements_to_abaqus, write_finite_elements_to_exodus,
     write_finite_elements_to_mesh, write_finite_elements_to_vtk, Blocks, Connectivity, Coordinates,
+    Smoothing,
 };
 use pyo3::prelude::*;
 
@@ -45,11 +46,40 @@ impl FiniteElements {
         ))
     }
     /// Smooths the nodal coordinates according to the provided smoothing method.
-    pub fn smooth(&mut self, method: String) -> Result<(), PyIntermediateError> {
-        Ok(Err(format!(
-            "Invalid smoothing method {} specified.",
-            method
-        ))?)
+    #[pyo3(signature = (method="Taubin", hierarchical=false, iterations=10, pass_band=0.1, scale=0.6307))]
+    pub fn smooth(
+        &mut self,
+        method: &str,
+        hierarchical: bool,
+        iterations: usize,
+        pass_band: f64,
+        scale: f64,
+    ) -> Result<(), PyIntermediateError> {
+        let mut finite_elements = super::FiniteElements::from_data(
+            self.element_blocks.clone(),
+            self.element_node_connectivity.clone(),
+            self.nodal_coordinates.clone(),
+        );
+        finite_elements.calculate_node_element_connectivity()?;
+        finite_elements.calculate_node_node_connectivity()?;
+        if hierarchical {
+            finite_elements.calculate_nodal_hierarchy()?;
+        }
+        finite_elements.calculate_nodal_influencers();
+        match method {
+            "Gauss" | "gauss" | "Gaussian" | "gaussian" | "Laplacian" | "Laplace" | "laplacian"
+            | "laplace" => {
+                finite_elements.smooth(Smoothing::Laplacian(iterations, scale))?;
+            }
+            "Taubin" | "taubin" => {
+                finite_elements.smooth(Smoothing::Taubin(iterations, pass_band, scale))?;
+            }
+            _ => return Err(format!("Invalid smoothing method {} specified.", method))?,
+        }
+        self.element_blocks = finite_elements.element_blocks;
+        self.element_node_connectivity = finite_elements.element_node_connectivity;
+        self.nodal_coordinates = finite_elements.nodal_coordinates;
+        Ok(())
     }
     /// Writes the finite elements data to a new Exodus file.
     pub fn write_exo(&self, file_path: &str) -> Result<(), PyIntermediateError> {
diff --git a/tests/fem.py b/tests/fem.py
index a5c0f77..bd1faf9 100644
--- a/tests/fem.py
+++ b/tests/fem.py
@@ -3,10 +3,20 @@
 remove = [0]
 scale = [1, 1, 1]
 translate = [0, 0, 0]
+voxels = Voxels.from_npy('tests/input/letter_f_3d.npy')
 
 
-def test_write_inp_letter_f_3d():
-    voxels = Voxels.from_spn('tests/input/letter_f_3d.spn', [4, 5, 3])
+def test_smooth_laplace():
+    fem = voxels.as_finite_elements(remove, scale, translate)
+    fem.smooth(method='Laplace')
+
+
+def test_smooth_taubin():
+    fem = voxels.as_finite_elements(remove, scale, translate)
+    fem.smooth(method='Taubin')
+
+
+def test_write_inp():
     fem = voxels.as_finite_elements(remove, scale, translate)
     inp = 'target/letter_f_3d.inp'
     fem.write_inp(inp)
@@ -24,6 +34,11 @@ def test_write_inp_letter_f_3d():
             line = file.readline()
 
 
+def test_write_exo():
+    fem = voxels.as_finite_elements(remove, scale, translate)
+    fem.write_exo('target/letter_f_3d.exo')
+
+
 def test_write_inp_sparse():
     voxels = Voxels.from_spn('tests/input/sparse.spn', [5, 5, 5])
     fem = voxels.as_finite_elements(remove, scale, translate)
@@ -41,3 +56,18 @@ def test_write_inp_sparse():
         while line != '':
             assert gold.readline() == line
             line = file.readline()
+
+
+def test_write_mesh():
+    fem = voxels.as_finite_elements(remove, scale, translate)
+    fem.write_mesh('target/letter_f_3d.mesh')
+
+
+def test_write_metrics():
+    fem = voxels.as_finite_elements(remove, scale, translate)
+    fem.write_metrics('target/letter_f_3d.csv')
+
+
+def test_write_vtk():
+    fem = voxels.as_finite_elements(remove, scale, translate)
+    fem.write_vtk('target/letter_f_3d.vtk')
diff --git a/tests/voxel.py b/tests/voxel.py
index dccb7e1..b5adfb0 100644
--- a/tests/voxel.py
+++ b/tests/voxel.py
@@ -7,3 +7,15 @@ def test_write_npy():
           .write_npy('target/letter_f_3d.npy')
     assert (np.load('tests/input/letter_f_3d.npy') ==
             np.load('target/letter_f_3d.npy')).all()
+
+
+def test_write_spn():
+    input = 'tests/input/letter_f_3d.spn'
+    output = 'target/letter_f_3d.spn'
+    Voxels.from_spn(input, [4, 5, 3]) \
+          .write_spn(output)
+    with open(input) as gold, open(output) as file:
+        line = file.readline()
+        while line != '':
+            assert gold.readline() == line
+            line = file.readline()