diff --git a/examples/Test_Plot.ipynb b/examples/Test_Plot.ipynb
index f843623..61bed89 100644
--- a/examples/Test_Plot.ipynb
+++ b/examples/Test_Plot.ipynb
@@ -44,14 +44,14 @@
     "fem = Solver(tria)\n",
     "evals, evecs = fem.eigs(k=3)\n",
     "evDict = dict()\n",
-    "evDict['Refine'] = 0\n",
-    "evDict['Degree'] = 1\n",
-    "evDict['Dimension'] = 2\n",
-    "evDict['Elements'] = len(tria.t)\n",
-    "evDict['DoF'] = len(tria.v)\n",
-    "evDict['NumEW'] = 3\n",
-    "evDict['Eigenvalues'] = evals\n",
-    "evDict['Eigenvectors'] = evecs\n",
+    "evDict[\"Refine\"] = 0\n",
+    "evDict[\"Degree\"] = 1\n",
+    "evDict[\"Dimension\"] = 2\n",
+    "evDict[\"Elements\"] = len(tria.t)\n",
+    "evDict[\"DoF\"] = len(tria.v)\n",
+    "evDict[\"NumEW\"] = 3\n",
+    "evDict[\"Eigenvalues\"] = evals\n",
+    "evDict[\"Eigenvectors\"] = evecs\n",
     "io.write_ev(\"../data/cubeTria.ev\", evDict)"
    ]
   },
@@ -68,7 +68,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plot.plot_tria_mesh(tria, vfunc=evecs[:,1], xrange=None, yrange=None, zrange=None, showcaxis=False, caxis=None)"
+    "plot.plot_tria_mesh(\n",
+    "    tria,\n",
+    "    vfunc=evecs[:, 1],\n",
+    "    xrange=None,\n",
+    "    yrange=None,\n",
+    "    zrange=None,\n",
+    "    showcaxis=False,\n",
+    "    caxis=None,\n",
+    ")"
    ]
   },
   {
@@ -84,7 +92,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plot.plot_tria_mesh(tria, vfunc=evecs[:,1], xrange=[-2, 2], yrange=[-2, 2], zrange=[-2, 2], showcaxis=True, caxis=[-0.3, 0.5])"
+    "plot.plot_tria_mesh(\n",
+    "    tria,\n",
+    "    vfunc=evecs[:, 1],\n",
+    "    xrange=[-2, 2],\n",
+    "    yrange=[-2, 2],\n",
+    "    zrange=[-2, 2],\n",
+    "    showcaxis=True,\n",
+    "    caxis=[-0.3, 0.5],\n",
+    ")"
    ]
   },
   {
@@ -111,14 +127,14 @@
     "fem = Solver(tetra)\n",
     "evals, evecs = fem.eigs(k=3)\n",
     "evDict = dict()\n",
-    "evDict['Refine'] = 0\n",
-    "evDict['Degree'] = 1\n",
-    "evDict['Dimension'] = 2\n",
-    "evDict['Elements'] = len(tetra.t)\n",
-    "evDict['DoF'] = len(tetra.v)\n",
-    "evDict['NumEW'] = 3\n",
-    "evDict['Eigenvalues'] = evals\n",
-    "evDict['Eigenvectors'] = evecs\n",
+    "evDict[\"Refine\"] = 0\n",
+    "evDict[\"Degree\"] = 1\n",
+    "evDict[\"Dimension\"] = 2\n",
+    "evDict[\"Elements\"] = len(tetra.t)\n",
+    "evDict[\"DoF\"] = len(tetra.v)\n",
+    "evDict[\"NumEW\"] = 3\n",
+    "evDict[\"Eigenvalues\"] = evals\n",
+    "evDict[\"Eigenvectors\"] = evecs\n",
     "io.write_ev(\"../data/cubeTetra.ev\", evDict)"
    ]
   },
@@ -135,7 +151,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plot.plot_tet_mesh(tetra, vfunc=evecs[:,1], xrange=None, yrange=None, zrange=None, showcaxis=False, caxis=None)"
+    "plot.plot_tet_mesh(\n",
+    "    tetra,\n",
+    "    vfunc=evecs[:, 1],\n",
+    "    xrange=None,\n",
+    "    yrange=None,\n",
+    "    zrange=None,\n",
+    "    showcaxis=False,\n",
+    "    caxis=None,\n",
+    ")"
    ]
   },
   {
@@ -151,7 +175,16 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plot.plot_tet_mesh(tetra, cutting=('f>0'), vfunc=evecs[:,1], xrange=[-2, 2], yrange=[-2, 2], zrange=[-2, 2], showcaxis=True, caxis=[-0.3, 0.5])"
+    "plot.plot_tet_mesh(\n",
+    "    tetra,\n",
+    "    cutting=(\"f>0\"),\n",
+    "    vfunc=evecs[:, 1],\n",
+    "    xrange=[-2, 2],\n",
+    "    yrange=[-2, 2],\n",
+    "    zrange=[-2, 2],\n",
+    "    showcaxis=True,\n",
+    "    caxis=[-0.3, 0.5],\n",
+    ")"
    ]
   }
  ],
diff --git a/examples/Test_ShapeDNA.ipynb b/examples/Test_ShapeDNA.ipynb
index e646978..26a14d1 100644
--- a/examples/Test_ShapeDNA.ipynb
+++ b/examples/Test_ShapeDNA.ipynb
@@ -89,8 +89,8 @@
    "source": [
     "# compute eigenvalues and eigenvectors for tria mesh\n",
     "ev = shapedna.compute_shapedna(tria, k=3)\n",
-    "ev['Eigenvectors']\n",
-    "ev['Eigenvalues']"
+    "ev[\"Eigenvectors\"]\n",
+    "ev[\"Eigenvalues\"]"
    ]
   },
   {
@@ -218,8 +218,8 @@
    "source": [
     "# compute eigenvalues and eigenvectors for tet mesh\n",
     "evTet = shapedna.compute_shapedna(tet, k=3)\n",
-    "evTet['Eigenvectors']\n",
-    "evTet['Eigenvalues']"
+    "evTet[\"Eigenvectors\"]\n",
+    "evTet[\"Eigenvalues\"]"
    ]
   },
   {
diff --git a/examples/Test_TetMesh.ipynb b/examples/Test_TetMesh.ipynb
index 7ee13f2..06d9d4e 100644
--- a/examples/Test_TetMesh.ipynb
+++ b/examples/Test_TetMesh.ipynb
@@ -30,12 +30,32 @@
    "outputs": [],
    "source": [
     "# cube9 (cube with center node)\n",
-    "points = [[0,0,0], [1,0,0], [1,1,0], [0,1,0], [0,0,1],\n",
-    "         [1,0,1], [1,1,1], [0,1,1], [0.5,0.5,0.5]]\n",
-    "tets = [[0,5,8,1], [0,4,5,8], [2,5,6,8], [1,5,2,8],\n",
-    "        [6,7,3,8], [6,3,2,8], [0,3,4,8], [3,7,4,8],\n",
-    "        [0,1,2,8], [0,2,3,8], [4,6,5,8], [4,7,6,8]]\n",
-    "T = TetMesh(points,tets)"
+    "points = [\n",
+    "    [0, 0, 0],\n",
+    "    [1, 0, 0],\n",
+    "    [1, 1, 0],\n",
+    "    [0, 1, 0],\n",
+    "    [0, 0, 1],\n",
+    "    [1, 0, 1],\n",
+    "    [1, 1, 1],\n",
+    "    [0, 1, 1],\n",
+    "    [0.5, 0.5, 0.5],\n",
+    "]\n",
+    "tets = [\n",
+    "    [0, 5, 8, 1],\n",
+    "    [0, 4, 5, 8],\n",
+    "    [2, 5, 6, 8],\n",
+    "    [1, 5, 2, 8],\n",
+    "    [6, 7, 3, 8],\n",
+    "    [6, 3, 2, 8],\n",
+    "    [0, 3, 4, 8],\n",
+    "    [3, 7, 4, 8],\n",
+    "    [0, 1, 2, 8],\n",
+    "    [0, 2, 3, 8],\n",
+    "    [4, 6, 5, 8],\n",
+    "    [4, 7, 6, 8],\n",
+    "]\n",
+    "T = TetMesh(points, tets)"
    ]
   },
   {
diff --git a/examples/Test_TetMesh_Geodesics.ipynb b/examples/Test_TetMesh_Geodesics.ipynb
index 5710b44..3050301 100644
--- a/examples/Test_TetMesh_Geodesics.ipynb
+++ b/examples/Test_TetMesh_Geodesics.ipynb
@@ -84,9 +84,9 @@
     }
    ],
    "source": [
-    "T = TetMesh.read_vtk('../data/cubeTetra.vtk')\n",
-    "#T.is_oriented()\n",
-    "T.orient_()\n"
+    "T = TetMesh.read_vtk(\"../data/cubeTetra.vtk\")\n",
+    "# T.is_oriented()\n",
+    "T.orient_()"
    ]
   },
   {
@@ -120,9 +120,9 @@
    "source": [
     "from lapy import Solver\n",
     "\n",
-    "fem = Solver(T,lump=True)\n",
+    "fem = Solver(T, lump=True)\n",
     "\n",
-    "evals, evec = fem.eigs(10)\n"
+    "evals, evec = fem.eigs(10)"
    ]
   },
   {
@@ -152,11 +152,20 @@
    },
    "outputs": [],
    "source": [
-    "evnum=1\n",
-    "cutting = ['x<0.5']\n",
+    "evnum = 1\n",
+    "cutting = [\"x<0.5\"]\n",
     "# also here we comment all plots to reduce file size\n",
     "# uncomment and take a look\n",
-    "plot_tet_mesh(T,vfunc=evals[evnum]*evec[:,evnum],plot_edges=None,plot_levels=False,cutting=cutting,edge_color='rgb(50,50,50)',html_output=False,flatshading=True)\n"
+    "plot_tet_mesh(\n",
+    "    T,\n",
+    "    vfunc=evals[evnum] * evec[:, evnum],\n",
+    "    plot_edges=None,\n",
+    "    plot_levels=False,\n",
+    "    cutting=cutting,\n",
+    "    edge_color=\"rgb(50,50,50)\",\n",
+    "    html_output=False,\n",
+    "    flatshading=True,\n",
+    ")"
    ]
   },
   {
@@ -174,9 +183,9 @@
    "source": [
     "from lapy.diffgeo import compute_divergence, compute_gradient\n",
     "\n",
-    "grad = compute_gradient(T,evec[:,evnum])\n",
-    "divx = -compute_divergence(T,grad)\n",
-    "vfunc = Bi*divx\n"
+    "grad = compute_gradient(T, evec[:, evnum])\n",
+    "divx = -compute_divergence(T, grad)\n",
+    "vfunc = Bi * divx"
    ]
   },
   {
@@ -185,8 +194,17 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "cutting = ['x<0.5']\n",
-    "plot_tet_mesh(T,vfunc=vfunc,plot_edges=None,plot_levels=False,cutting=cutting,edge_color='rgb(50,50,50)',html_output=False,flatshading=True)\n"
+    "cutting = [\"x<0.5\"]\n",
+    "plot_tet_mesh(\n",
+    "    T,\n",
+    "    vfunc=vfunc,\n",
+    "    plot_edges=None,\n",
+    "    plot_levels=False,\n",
+    "    cutting=cutting,\n",
+    "    edge_color=\"rgb(50,50,50)\",\n",
+    "    html_output=False,\n",
+    "    flatshading=True,\n",
+    ")"
    ]
   },
   {
@@ -213,7 +231,7 @@
     }
    ],
    "source": [
-    "np.max(np.abs(vfunc-(evals[evnum]*evec[:,evnum])))"
+    "np.max(np.abs(vfunc - (evals[evnum] * evec[:, evnum])))"
    ]
   },
   {
@@ -252,9 +270,18 @@
     "tria = T.boundary_tria()\n",
     "bvert = np.unique(tria.t)\n",
     "\n",
-    "u = heat.diffusion(T,bvert,m=1)\n",
-    "cutting = ['x<0.5']\n",
-    "plot_tet_mesh(T,vfunc=u,plot_edges=None,plot_levels=True,cutting=cutting,edge_color='rgb(50,50,50)',html_output=False,flatshading=True)"
+    "u = heat.diffusion(T, bvert, m=1)\n",
+    "cutting = [\"x<0.5\"]\n",
+    "plot_tet_mesh(\n",
+    "    T,\n",
+    "    vfunc=u,\n",
+    "    plot_edges=None,\n",
+    "    plot_levels=True,\n",
+    "    cutting=cutting,\n",
+    "    edge_color=\"rgb(50,50,50)\",\n",
+    "    html_output=False,\n",
+    "    flatshading=True,\n",
+    ")"
    ]
   },
   {
@@ -273,7 +300,7 @@
     "# get gradients\n",
     "tfunc = compute_gradient(T, u)\n",
     "# flip and normalize\n",
-    "X = -tfunc / np.sqrt((tfunc ** 2).sum(1))[:, np.newaxis]\n",
+    "X = -tfunc / np.sqrt((tfunc**2).sum(1))[:, np.newaxis]\n",
     "X = np.nan_to_num(X)\n",
     "# compute divergence\n",
     "divx = compute_divergence(T, X)"
@@ -310,13 +337,13 @@
     "except ImportError:\n",
     "    useCholmod = False\n",
     "\n",
-    "A, B = fem.stiffness, fem.mass # computed above when creating Solver\n",
+    "A, B = fem.stiffness, fem.mass  # computed above when creating Solver\n",
+    "\n",
+    "H = A\n",
+    "b0 = -divx\n",
     "\n",
-    "H=A\n",
-    "b0=-divx\n",
-    "        \n",
     "# solve H x = b0\n",
-    "print(\"Matrix Format now: \"+H.getformat())\n",
+    "print(\"Matrix Format now: \" + H.getformat())\n",
     "if useCholmod:\n",
     "    print(\"Solver: cholesky decomp - performance optimal ...\")\n",
     "    chol = cholesky(H)\n",
@@ -346,9 +373,18 @@
     }
    ],
    "source": [
-    "cutting = ['x<0.5']\n",
-    "plot_tet_mesh(T,vfunc=x,plot_edges=None,plot_levels=True,cutting=cutting,edge_color='rgb(50,50,50)',html_output=False,flatshading=True)\n",
-    "max(x), 0.5*np.sqrt(3.0)"
+    "cutting = [\"x<0.5\"]\n",
+    "plot_tet_mesh(\n",
+    "    T,\n",
+    "    vfunc=x,\n",
+    "    plot_edges=None,\n",
+    "    plot_levels=True,\n",
+    "    cutting=cutting,\n",
+    "    edge_color=\"rgb(50,50,50)\",\n",
+    "    html_output=False,\n",
+    "    flatshading=True,\n",
+    ")\n",
+    "max(x), 0.5 * np.sqrt(3.0)"
    ]
   },
   {
@@ -382,15 +418,24 @@
     "from lapy.diffgeo import compute_geodesic_f\n",
     "\n",
     "tria = T.boundary_tria()\n",
-    "bvert=np.unique(tria.t)\n",
+    "bvert = np.unique(tria.t)\n",
     "\n",
     "# get heat diffusion\n",
-    "u = heat.diffusion(T,bvert, m=1)\n",
+    "u = heat.diffusion(T, bvert, m=1)\n",
     "\n",
-    "gu = compute_geodesic_f(T,u)\n",
+    "gu = compute_geodesic_f(T, u)\n",
     "\n",
-    "cutting = ['x<0.5']\n",
-    "plot_tet_mesh(T,vfunc=gu,plot_edges=None,plot_levels=True,cutting=cutting,edge_color='rgb(50,50,50)',html_output=False,flatshading=True)"
+    "cutting = [\"x<0.5\"]\n",
+    "plot_tet_mesh(\n",
+    "    T,\n",
+    "    vfunc=gu,\n",
+    "    plot_edges=None,\n",
+    "    plot_levels=True,\n",
+    "    cutting=cutting,\n",
+    "    edge_color=\"rgb(50,50,50)\",\n",
+    "    html_output=False,\n",
+    "    flatshading=True,\n",
+    ")"
    ]
   },
   {
@@ -407,15 +452,24 @@
    "outputs": [],
    "source": [
     "# test function is squared distance to each vertex\n",
-    "v1func =  T.v[:,0]* T.v[:,0] + T.v[:,1]* T.v[:,1] + T.v[:,2]* T.v[:,2]\n",
+    "v1func = T.v[:, 0] * T.v[:, 0] + T.v[:, 1] * T.v[:, 1] + T.v[:, 2] * T.v[:, 2]\n",
     "\n",
-    "grad = compute_gradient(T,v1func)\n",
+    "grad = compute_gradient(T, v1func)\n",
     "# glength = np.sqrt(np.sum(grad * grad, axis=1))\n",
     "# fcols=glength\n",
-    "fcols=grad[:,2]\n",
+    "fcols = grad[:, 2]\n",
     "# cutting = ['x<0.5']\n",
     "cutting = None\n",
-    "plot_tet_mesh(T,vfunc=None,tfunc=fcols,plot_edges=None,plot_levels=False,cutting=cutting,edge_color='rgb(50,50,50)',html_output=False)"
+    "plot_tet_mesh(\n",
+    "    T,\n",
+    "    vfunc=None,\n",
+    "    tfunc=fcols,\n",
+    "    plot_edges=None,\n",
+    "    plot_levels=False,\n",
+    "    cutting=cutting,\n",
+    "    edge_color=\"rgb(50,50,50)\",\n",
+    "    html_output=False,\n",
+    ")"
    ]
   },
   {
@@ -433,8 +487,18 @@
    "source": [
     "divx = compute_divergence(T, grad)\n",
     "divx2 = Bi * divx\n",
-    "cutting = ['z<0.5']\n",
-    "plot_tet_mesh(T,vfunc=divx2,plot_edges=True,plot_levels=False,cutting=cutting,edge_color='rgb(50,50,50)',html_output=False,flatshading=True,caxis=[0, 8])"
+    "cutting = [\"z<0.5\"]\n",
+    "plot_tet_mesh(\n",
+    "    T,\n",
+    "    vfunc=divx2,\n",
+    "    plot_edges=True,\n",
+    "    plot_levels=False,\n",
+    "    cutting=cutting,\n",
+    "    edge_color=\"rgb(50,50,50)\",\n",
+    "    html_output=False,\n",
+    "    flatshading=True,\n",
+    "    caxis=[0, 8],\n",
+    ")"
    ]
   },
   {
diff --git a/examples/Test_TriaMesh.ipynb b/examples/Test_TriaMesh.ipynb
index b07f9af..faa5db3 100644
--- a/examples/Test_TriaMesh.ipynb
+++ b/examples/Test_TriaMesh.ipynb
@@ -31,12 +31,32 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "points = [[0.0,0.0,0.0], [0,1,0], [1,1,0], [1,0,0],\n",
-    "          [0,0,1], [0,1,1], [1,1,1], [1,0,1]]\n",
-    "trias = [[0,1,2], [2,3,0], [4,5,6], [6,7,4], [0,4,7], [7,3,0],\n",
-    "         [0,4,5], [5,1,0], [1,5,6], [6,2,1], [3,7,6], [6,2,3]]\n",
+    "points = [\n",
+    "    [0.0, 0.0, 0.0],\n",
+    "    [0, 1, 0],\n",
+    "    [1, 1, 0],\n",
+    "    [1, 0, 0],\n",
+    "    [0, 0, 1],\n",
+    "    [0, 1, 1],\n",
+    "    [1, 1, 1],\n",
+    "    [1, 0, 1],\n",
+    "]\n",
+    "trias = [\n",
+    "    [0, 1, 2],\n",
+    "    [2, 3, 0],\n",
+    "    [4, 5, 6],\n",
+    "    [6, 7, 4],\n",
+    "    [0, 4, 7],\n",
+    "    [7, 3, 0],\n",
+    "    [0, 4, 5],\n",
+    "    [5, 1, 0],\n",
+    "    [1, 5, 6],\n",
+    "    [6, 2, 1],\n",
+    "    [3, 7, 6],\n",
+    "    [6, 2, 3],\n",
+    "]\n",
     "\n",
-    "T = TriaMesh(points,trias)"
+    "T = TriaMesh(points, trias)"
    ]
   },
   {
@@ -527,7 +547,7 @@
    ],
    "source": [
     "# Compute Volume (works only for oriented meshes)\n",
-    "T.volume() "
+    "T.volume()"
    ]
   },
   {
@@ -578,7 +598,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "T.normal_offset_(0.2*T.avg_edge_length())"
+    "T.normal_offset_(0.2 * T.avg_edge_length())"
    ]
   },
   {
@@ -594,7 +614,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "Tbdr = TriaMesh(T.v,T.t[2:,:])"
+    "Tbdr = TriaMesh(T.v, T.t[2:, :])"
    ]
   },
   {
diff --git a/examples/Test_TriaMesh_Geodesics.ipynb b/examples/Test_TriaMesh_Geodesics.ipynb
index 744c543..cc8a2e9 100644
--- a/examples/Test_TriaMesh_Geodesics.ipynb
+++ b/examples/Test_TriaMesh_Geodesics.ipynb
@@ -107,7 +107,7 @@
    "outputs": [],
    "source": [
     "q = T.tria_qualities()\n",
-    "plot_tria_mesh(T,plot_edges=True, tfunc=q)"
+    "plot_tria_mesh(T, plot_edges=True, tfunc=q)"
    ]
   },
   {
@@ -142,9 +142,9 @@
     "# compute first eigenfunction\n",
     "from lapy import Solver\n",
     "\n",
-    "fem = Solver(T,lump=True)\n",
+    "fem = Solver(T, lump=True)\n",
     "eval, evec = fem.eigs()\n",
-    "vfunc = evec[:,1]\n",
+    "vfunc = evec[:, 1]\n",
     "\n",
     "# also get A,B (lumped), and inverse of B (easy as it is diagonal due to lumping)\n",
     "A, B = fem.stiffness, fem.mass\n",
@@ -192,7 +192,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plot_tria_mesh(T,Bi*(A*vfunc),plot_edges=True)"
+    "plot_tria_mesh(T, Bi * (A * vfunc), plot_edges=True)"
    ]
   },
   {
@@ -208,7 +208,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plot_tria_mesh(T,eval[1]*vfunc,plot_edges=True)"
+    "plot_tria_mesh(T, eval[1] * vfunc, plot_edges=True)"
    ]
   },
   {
@@ -226,9 +226,9 @@
    "source": [
     "from lapy.diffgeo import compute_divergence, compute_gradient\n",
     "\n",
-    "grad = compute_gradient(T,vfunc)\n",
-    "divx = -compute_divergence(T,grad)\n",
-    "plot_tria_mesh(T,Bi*divx,plot_edges=True)"
+    "grad = compute_gradient(T, vfunc)\n",
+    "divx = -compute_divergence(T, grad)\n",
+    "plot_tria_mesh(T, Bi * divx, plot_edges=True)"
    ]
   },
   {
@@ -259,7 +259,7 @@
     "from lapy import heat\n",
     "\n",
     "bvert = T.boundary_loops()\n",
-    "u = heat.diffusion(T,bvert,m=1)"
+    "u = heat.diffusion(T, bvert, m=1)"
    ]
   },
   {
@@ -275,7 +275,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plot_tria_mesh(T,u,plot_edges=True,plot_levels=True)"
+    "plot_tria_mesh(T, u, plot_edges=True, plot_levels=True)"
    ]
   },
   {
@@ -294,10 +294,10 @@
     "import numpy as np\n",
     "\n",
     "# compute gradient of heat diffusion\n",
-    "tfunc = compute_gradient(T,u)\n",
+    "tfunc = compute_gradient(T, u)\n",
     "\n",
     "# normalize gradient\n",
-    "X = -tfunc / np.sqrt((tfunc**2).sum(1))[:,np.newaxis]\n",
+    "X = -tfunc / np.sqrt((tfunc**2).sum(1))[:, np.newaxis]\n",
     "X = np.nan_to_num(X)"
    ]
   },
@@ -314,7 +314,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "divx = compute_divergence(T,X)"
+    "divx = compute_divergence(T, X)"
    ]
   },
   {
@@ -349,15 +349,15 @@
     "except ImportError:\n",
     "    useCholmod = False\n",
     "\n",
-    "fem = Solver(T,lump=True)\n",
+    "fem = Solver(T, lump=True)\n",
     "A, B = fem.stiffness, fem.mass\n",
     "\n",
-    "H=-A\n",
-    "b0=divx\n",
-    "        \n",
+    "H = -A\n",
+    "b0 = divx\n",
+    "\n",
     "# solve H x = b0\n",
-    "# we don't need the B matrix here, as divx is the integrated divergence \n",
-    "print(\"Matrix Format now: \"+H.getformat())\n",
+    "# we don't need the B matrix here, as divx is the integrated divergence\n",
+    "print(\"Matrix Format now: \" + H.getformat())\n",
     "if useCholmod:\n",
     "    print(\"Solver: cholesky decomp - performance optimal ...\")\n",
     "    chol = cholesky(H)\n",
@@ -368,7 +368,7 @@
     "    x = lu.solve(b0)\n",
     "\n",
     "# remove shift\n",
-    "x = x-min(x)"
+    "x = x - min(x)"
    ]
   },
   {
@@ -384,7 +384,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plot_tria_mesh(T,x,plot_edges=True,plot_levels=True)"
+    "plot_tria_mesh(T, x, plot_edges=True, plot_levels=True)"
    ]
   },
   {
@@ -412,7 +412,7 @@
    ],
    "source": [
     "# max distance (smoothed)\n",
-    "(max(x), np.sqrt(2)/2)"
+    "(max(x), np.sqrt(2) / 2)"
    ]
   },
   {
@@ -437,9 +437,9 @@
     }
    ],
    "source": [
-    "vf = fem.poisson(-Bi*divx)\n",
+    "vf = fem.poisson(-Bi * divx)\n",
     "vf = vf - min(vf)\n",
-    "plot_tria_mesh(T,vf,plot_edges=True,plot_levels=True)"
+    "plot_tria_mesh(T, vf, plot_edges=True, plot_levels=True)"
    ]
   },
   {
@@ -466,7 +466,7 @@
     }
    ],
    "source": [
-    "max(abs(vf-x))"
+    "max(abs(vf - x))"
    ]
   },
   {
@@ -494,8 +494,8 @@
    "source": [
     "from lapy.diffgeo import compute_geodesic_f\n",
     "\n",
-    "gf = compute_geodesic_f(T,u)\n",
-    "plot_tria_mesh(T,gf,plot_edges=True,plot_levels=True)"
+    "gf = compute_geodesic_f(T, u)\n",
+    "plot_tria_mesh(T, gf, plot_edges=True, plot_levels=True)"
    ]
   },
   {
@@ -523,7 +523,7 @@
    ],
    "source": [
     "# again should be the same\n",
-    "max(abs(gf-x))"
+    "max(abs(gf - x))"
    ]
   },
   {
@@ -552,8 +552,8 @@
     "# testing if we can rotate the function\n",
     "from lapy.diffgeo import compute_rotated_f\n",
     "\n",
-    "gf = compute_rotated_f(T,vf)\n",
-    "plot_tria_mesh(T,gf,plot_edges=True,plot_levels=True)"
+    "gf = compute_rotated_f(T, vf)\n",
+    "plot_tria_mesh(T, gf, plot_edges=True, plot_levels=True)"
    ]
   },
   {
@@ -572,11 +572,11 @@
    "outputs": [],
    "source": [
     "# Load your mesh here and uncomment. The mesh should have not too many flat regions (not a cube)\n",
-    "#from lapy.diffgeo import tria_mean_curvature_flow\n",
-    "#from lapy.plot import plot_tria_mesh\n",
-    "#T = TriaIO.import_off(\"../data/???\")\n",
-    "#T2 = tria_mean_curvature_flow(T)\n",
-    "#plot_tria_mesh(T2,plot_edges=True,plot_levels=True)"
+    "# from lapy.diffgeo import tria_mean_curvature_flow\n",
+    "# from lapy.plot import plot_tria_mesh\n",
+    "# T = TriaIO.import_off(\"../data/???\")\n",
+    "# T2 = tria_mean_curvature_flow(T)\n",
+    "# plot_tria_mesh(T2,plot_edges=True,plot_levels=True)"
    ]
   }
  ],
diff --git a/lapy/diffgeo.py b/lapy/diffgeo.py
index d9518e9..ef34c2b 100644
--- a/lapy/diffgeo.py
+++ b/lapy/diffgeo.py
@@ -1,7 +1,7 @@
 """Differential Geometry Functions for meshes.
 
 This module includes gradient, divergence, curvature, geodesics,
-mean curvature flow etc. 
+mean curvature flow etc.
 
 Note, the interface is not yet final, some functions are split into
 tet and tria versions.
@@ -369,8 +369,8 @@ def tria_compute_rotated_f(tria, vfunc):
     # since the solution is ill defined (addition of constant)
     # we specify an arbitrary boundary condition at a single vertex to
     # remove that degree of freedom
-    dtup = ( np.array([0]), np.array([0.0]))
-    vf = fem.poisson(divf,dtup)
+    dtup = (np.array([0]), np.array([0.0]))
+    vf = fem.poisson(divf, dtup)
     return vf
 
 
diff --git a/lapy/heat.py b/lapy/heat.py
index ce3d8de..9e6de09 100644
--- a/lapy/heat.py
+++ b/lapy/heat.py
@@ -1,7 +1,7 @@
 """Functions for computing heat kernel and diffusion.
 
 Inputs are eigenvalues and eigenvectors (for heat kernel) and the
-mesh geometries (tet or tria mesh) for heat diffusion. 
+mesh geometries (tet or tria mesh) for heat diffusion.
 """
 
 import importlib
diff --git a/lapy/shapedna.py b/lapy/shapedna.py
index f222100..973f9c6 100644
--- a/lapy/shapedna.py
+++ b/lapy/shapedna.py
@@ -1,8 +1,8 @@
 """Functions for computing and comparing Laplace spectra.
 
 Includes code for solving the anisotropic Laplace-Beltrami eigenvalue
-problem as well as functions for normalization and comparison of 
-Laplace spectra. 
+problem as well as functions for normalization and comparison of
+Laplace spectra.
 """
 
 import numpy as np
diff --git a/lapy/solver.py b/lapy/solver.py
index 5c17899..a93c9a8 100644
--- a/lapy/solver.py
+++ b/lapy/solver.py
@@ -706,7 +706,7 @@ def poisson(self, h=0.0, dtup=(), ntup=()):  # poissonSolver
                 "h should be either scalar or column vector with row num of A"
             )
         if h.ndim == 1:
-            h = h[:,np.newaxis]
+            h = h[:, np.newaxis]
         # create vector d
         didx = []
         dvec = []
diff --git a/lapy/tria_mesh.py b/lapy/tria_mesh.py
index 8871b7a..f8757d3 100644
--- a/lapy/tria_mesh.py
+++ b/lapy/tria_mesh.py
@@ -1120,3 +1120,281 @@ def smooth_(self, n=1):
         vfunc = self.smooth_vfunc(self.v, n)
         self.v = vfunc
         return
+
+    def level_length(self, vfunc, level):
+        """Compute the length of level sets.
+
+        Parameters
+        ----------
+        vfunc : array
+            Float vector of values at vertices (here only scalar function 1D).
+        level : float | array
+            Level set value or array of level values.
+
+        Returns
+        -------
+        length : float | array
+            Length of level set (or array of lengths).
+        """
+        if vfunc.ndim > 1:
+            raise ValueError(f"vfunc needs to be 1-dim, but is {vfunc.ndim}-dim!")
+        levels = np.atleast_1d(level)
+        ll = np.empty((levels.size, 0))
+        for lnext in levels:
+            # get intersecting triangles
+            intersect = vfunc[self.t] > lnext
+            t_idx = np.where(
+                np.logical_or(
+                    np.sum(intersect, axis=1) == 1, np.sum(intersect, axis=1) == 2
+                )
+            )[0]
+            # reduce to triangles that intersect with level:
+            t_level = self.t[t_idx, :]
+            intersect = intersect[t_idx, :]
+            # trias have one vertex on one side and two on the other side of the
+            # level set. Here we invert trias with two true values, so that single
+            # vertex is true
+            intersect[np.sum(intersect, axis=1) > 1, :] = np.logical_not(
+                intersect[np.sum(intersect, axis=1) > 1, :]
+            )
+            # get idx within tria with single vertex:
+            idx_single = np.argmax(intersect, axis=1)
+            idx_o1 = (idx_single + 1) % 3
+            idx_o2 = (idx_single + 2) % 3
+            # get global idx
+            gidx0 = t_level[np.arange(t_level.shape[0]), idx_single]
+            gidx1 = t_level[np.arange(t_level.shape[0]), idx_o1]
+            gidx2 = t_level[np.arange(t_level.shape[0]), idx_o2]
+            # determine fraction along edges (for each triangle)
+            xl1 = (lnext - vfunc[gidx0]) / (vfunc[gidx1] - vfunc[gidx0])
+            xl2 = (lnext - vfunc[gidx0]) / (vfunc[gidx2] - vfunc[gidx0])
+            # determine points on the two edges (for each triangle)
+            p1 = (1 - xl1)[:, np.newaxis] * self.v[gidx0] + xl1[:, np.newaxis] * self.v[
+                gidx1
+            ]
+            p2 = (1 - xl2)[:, np.newaxis] * self.v[gidx0] + xl2[:, np.newaxis] * self.v[
+                gidx2
+            ]
+            # compute edge length between the points
+            lls = np.sqrt(((p1 - p2) ** 2).sum(1))
+            ll = np.append(ll, lls.sum())
+        if ll.size > 1:
+            return ll
+        elif ll.size > 0:
+            return ll[0]
+        else:
+            raise ValueError("No lengths computed, should never get here.")
+
+    @staticmethod
+    def __reduce_edges_to_path(edges, start_idx=None, get_edge_idx=False):
+        """Reduce undirected unsorted edges (index pairs) to path (index array).
+
+        Parameters
+        ----------
+        edges : list of shape (n, 2)
+            Pairs of positive integer node indices representing a undirected edges.
+            All indices from 0 to max(edges) must be used and graph needs to be
+            connected. Nodes cannot have more than 2 neighbors. Also needs exactly
+            two endnodes (nodes with only one neighbor). Tip for closed loops, cut
+            it open before passing to this function by removing one edge.
+        start_idx : int, default: None
+            Node with only one neighbor to start path. Optional, if not passed the
+            node with the smaller index will be selected as start point.
+        get_edge_idx : bool, default: False
+            Also return index of edge into edges for each path segment. This list
+            has length n, while path has length n+1.
+
+        Returns
+        -------
+        path: array
+            Array of length n+1 containing indices as single path from start to
+            endpoint.
+        edge_idx: array
+            Array of length n containing corresponding edge idx into edges for each
+            path segment. Only passed if get_edges_idx is True.
+        """
+        from scipy.sparse.csgraph import shortest_path
+
+        # Extract node indices and create a sparse adjacency matrix
+        edges = np.array(edges)
+        i = np.column_stack((edges[:, 0], edges[:, 1])).reshape(-1)
+        j = np.column_stack((edges[:, 1], edges[:, 0])).reshape(-1)
+        dat = np.ones(i.shape)
+        n = edges.max() + 1
+        adj_matrix = sparse.csr_matrix((dat, (i, j)), shape=(n, n))
+        # Find the degree of each node, sum over rows to get degree
+        degrees = np.asarray(adj_matrix.sum(axis=1)).ravel()
+        endpoints = np.where(degrees == 1)[0]
+        if len(endpoints) != 2:
+            raise ValueError(
+                "The graph does not have exactly two endpoints; invalid input."
+            )
+        if not start_idx:
+            start_idx = endpoints[0]
+        else:
+            if not np.isin(start_idx, endpoints):
+                raise ValueError(
+                    f"start_idx {start_idx} must be one of the endpoints {endpoints}."
+                )
+        # Traverse the graph by computing shortest path
+        dist_matrix = shortest_path(
+            csgraph=adj_matrix,
+            directed=False,
+            indices=start_idx,
+            return_predecessors=False,
+        )
+        if np.isinf(dist_matrix).any():
+            raise ValueError(
+                "Ensure connected graph with indices from 0 to max_idx without gaps."
+            )
+        # sort indices according to distance form start
+        path = dist_matrix.argsort()
+        # get edge idx of each segment from original list
+        enum = edges.shape[0]
+        dat = np.arange(enum) + 1
+        dat = np.column_stack((dat, dat)).reshape(-1)
+        eidx_matrix = sparse.csr_matrix((dat, (i, j)), shape=(n, n))
+        ei = path[0:-1]
+        ej = path[(np.arange(path.size - 1) + 1)]
+        eidx = np.asarray(eidx_matrix[ei, ej] - 1).ravel()
+        if get_edge_idx:
+            return path, eidx
+        else:
+            return path
+
+    @staticmethod
+    def __resample_polygon(path3d: np.ndarray, n_points: int = 100) -> np.ndarray:
+        # Helper: Cumulative Euclidean distance between successive polygon points.
+        # This will be the "x" for interpolation
+        d = np.cumsum(np.r_[0, np.sqrt((np.diff(path3d, axis=0) ** 2).sum(axis=1))])
+        # get linearly spaced points along the cumulative Euclidean distance
+        d_sampled = np.linspace(0, d.max(), n_points)
+        # interpolate x and y coordinates
+        xy_interp = np.c_[
+            np.interp(d_sampled, d, path3d[:, 0]),
+            np.interp(d_sampled, d, path3d[:, 1]),
+            np.interp(d_sampled, d, path3d[:, 2]),
+        ]
+        return xy_interp
+
+    @staticmethod
+    def __iterative_resample_polygon(
+        path3d: np.ndarray, n_points: int = 100, n_iter: int = 3) -> np.ndarray:
+        # Helper: resample multiple times to numerically stabilize the result to be
+        # truly equidistant
+        path3d_resampled = TriaMesh.__resample_polygon(path3d, n_points)
+        for _ in range(n_iter - 1):
+            path3d_resampled = TriaMesh.__resample_polygon(path3d_resampled, n_points)
+        return path3d_resampled
+
+    def level_path(self, vfunc, level, get_tria_idx=False, n_points=None):
+        """Extract levelset as a path.
+
+        Note: Only works for level sets that represent a single path with
+        exactly one start and one endpoint.
+
+        Parameters
+        ----------
+        vfunc : array
+            Float vector of values at vertices (here only scalar function 1D).
+        level : float
+            Level set value.
+        get_tria_idx : bool, default: False
+            Also return a list of triangle indices for each edge, default False.
+        n_points : int
+            Resample level set into n equidistant points. Cannot be combined
+            with get_tria_idx=True.
+
+        Returns
+        -------
+        path: array
+            Array of shape (n,3) containing coordinates of vertices on level path.
+        length : float
+            Length of level set.
+        tria_idx : array
+            Array of triangle index for each edge segment on the path (length n-1
+            if path is length n). Will only be returned if get_tria_idx is True.
+        """
+        if vfunc.ndim > 1:
+            raise ValueError(f"vfunc needs to be 1-dim, but is {vfunc.ndim}-dim!")
+        # get intersecting triangles
+        intersect = vfunc[self.t] > level
+        t_idx = np.where(
+            np.logical_or(
+                np.sum(intersect, axis=1) == 1, np.sum(intersect, axis=1) == 2
+            )
+        )[0]
+        # reduce to triangles that intersect with level:
+        t_level = self.t[t_idx, :]
+        intersect = intersect[t_idx, :]
+        # trias have one vertex on one side and two on the other side of the level set
+        # invert trias with two true values, so that single vertex is true
+        intersect[np.sum(intersect, axis=1) > 1, :] = np.logical_not(
+            intersect[np.sum(intersect, axis=1) > 1, :]
+        )
+        # get idx within tria with single vertex:
+        idx_single = np.argmax(intersect, axis=1)
+        idx_o1 = (idx_single + 1) % 3
+        idx_o2 = (idx_single + 2) % 3
+        # get global idx
+        gidx0 = t_level[np.arange(t_level.shape[0]), idx_single]
+        gidx1 = t_level[np.arange(t_level.shape[0]), idx_o1]
+        gidx2 = t_level[np.arange(t_level.shape[0]), idx_o2]
+        # sort edge indices (rows are trias, cols are the two vertex ids)
+        gg1 = np.sort(
+            np.concatenate((gidx0[:, np.newaxis], gidx1[:, np.newaxis]), axis=1)
+        )
+        gg2 = np.sort(
+            np.concatenate((gidx0[:, np.newaxis], gidx2[:, np.newaxis]), axis=1)
+        )
+        # concatenate all and get unique ones
+        gg = np.concatenate((gg1, gg2), axis=0)
+        gg_unique = np.unique(gg, axis=0)
+        # generate level set intersection points for unique edges
+        xl = (level - vfunc[gg_unique[:, 0]]) / (
+            vfunc[gg_unique[:, 1]] - vfunc[gg_unique[:, 0]]
+        )
+        p = (1 - xl)[:, np.newaxis] * self.v[gg_unique[:, 0]] + xl[
+            :, np.newaxis
+        ] * self.v[gg_unique[:, 1]]
+        # fill sparse matrix with new point indices (+1 to distinguish from zero)
+        A = sparse.csc_matrix(
+            (np.arange(gg_unique.shape[0]) + 1, (gg_unique[:, 0], gg_unique[:, 1]))
+        )
+        # for each tria create one edge via lookup in matrix
+        edge_idxs = (
+            np.concatenate((A[gg1[:, 0], gg1[:, 1]], A[gg2[:, 0], gg2[:, 1]]), axis=0).T
+            - 1
+        )
+        # lengths computation
+        p1 = np.squeeze(p[edge_idxs[:, 0]])
+        p2 = np.squeeze(p[edge_idxs[:, 1]])
+        llength = np.sqrt(((p1 - p2) ** 2).sum(1)).sum()
+        # compute path from unordered, not-directed edge list
+        # and return path as list of points, and path length
+        if get_tria_idx:
+            path, edge_idx = TriaMesh.__reduce_edges_to_path(
+                edge_idxs, get_edge_idx=get_tria_idx
+            )
+            # translate local edge id to global tria id
+            tria_idx = t_idx[edge_idx]
+        else:
+            path = TriaMesh.__reduce_edges_to_path(edge_idxs, get_tria_idx)
+        # remove duplicate vertices (happens when levelset hits a vertex almost
+        # perfectly)
+        path3d = p[path, :]
+        dd = ((path3d[0:-1, :] - path3d[1:, :]) ** 2).sum(1)
+        # append 1 (never delete last node, if identical to the one before, we delete
+        # the one before)
+        dd = np.append(dd, 1)
+        eps = 0.000001
+        path3d = path3d[dd > eps, :]
+        if get_tria_idx:
+            if n_points:
+                raise ValueError("n_points cannot be combined with get_tria_idx=True.")
+            tria_idx = tria_idx[dd[:-1] > eps]
+            return path3d, llength, tria_idx
+        else:
+            if n_points:
+                path3d = TriaMesh.__iterative_resample_polygon(path3d, n_points)
+            return path3d, llength
diff --git a/pyproject.toml b/pyproject.toml
index 5b3f128..ab6a4c1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -7,7 +7,7 @@ build-backend = 'setuptools.build_meta'
 
 [project]
 name = 'lapy'
-version = '1.2.0-dev'
+version = '1.2.0'
 description = 'A package for differential geometry on meshes (Laplace, FEM)'
 readme = 'README.md'
 license = {file = 'LICENSE'}