diff --git a/.gitignore b/.gitignore
index abf76f98e..e632fe94d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -258,3 +258,5 @@ IPCToolkitOptions.cmake
 src/ipc/config.hpp
 
 docs/_doxygen
+
+notebooks/*.[ch]pp
diff --git a/notebooks/area_weights.ipynb b/notebooks/area_weights.ipynb
new file mode 100644
index 000000000..6cf24b40a
--- /dev/null
+++ b/notebooks/area_weights.ipynb
@@ -0,0 +1,166 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy\n",
+    "import sympy\n",
+    "from sympy.printing import ccode\n",
+    "\n",
+    "from generate_cpp_code import *\n",
+    "from utils import norm"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "funcs = []"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Edge Length"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def edge_len(e0, e1):\n",
+    "    return norm(e1 - e0)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "x_2D = numpy.array(sympy.symbols(\" \".join([f\"e{i}_{d}\" for i in range(2) for d in \"xy\"])))\n",
+    "x_3D = numpy.array(sympy.symbols(\" \".join([f\"e{i}_{d}\" for i in range(2) for d in \"xyz\"])))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "funcs.append(CXXGradientGenerator(\n",
+    "    edge_len(*numpy.split(x_2D, 2)), x_2D, \n",
+    "    \"edge_length_gradient_2D\", out_param_name=\"dA\"))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "funcs.append(CXXGradientGenerator(\n",
+    "    edge_len(*numpy.split(x_3D, 2)), x_3D, \n",
+    "    \"edge_length_gradient_3D\", out_param_name=\"dA\"))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Triangle Area"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def triangle_area(t0, t1, t2):\n",
+    "    n = numpy.cross(t1 - t0, t2 - t0)\n",
+    "    return norm(n) / 2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "x = numpy.array(sympy.symbols(\" \".join([f\"t{i}_{d}\" for i in range(3) for d in \"xyz\"])))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "funcs.append(CXXGradientGenerator(\n",
+    "    triangle_area(*numpy.split(x, 3)), x,\n",
+    "    \"triangle_area_gradient\", out_param_name=\"dA\"))"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Generate Code"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "generate_hpp_file(funcs, \"area_gradient.hpp\")\n",
+    "generate_cpp_file(funcs, \"area_gradient.cpp\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.1"
+  },
+  "orig_nbformat": 4,
+  "vscode": {
+   "interpreter": {
+    "hash": "5c7b89af1651d0b8571dde13640ecdccf7d5a6204171d6ab33e7c296e100e08a"
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/notebooks/closest_point_grad.ipynb b/notebooks/closest_point_grad.ipynb
new file mode 100644
index 000000000..ec8b5ecc6
--- /dev/null
+++ b/notebooks/closest_point_grad.ipynb
@@ -0,0 +1,255 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy\n",
+    "from sympy import *\n",
+    "from sympy.matrices import *\n",
+    "\n",
+    "from generate_cpp_code import *\n",
+    "from utils import norm, jacobian"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Function(\"F\")(MatrixSymbol(\"A\", 2, 2)).diff()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "generators = []"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Point-Edge"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for dim in [2, 3]:\n",
+    "    p = numpy.array(symbols(\" \".join(\"p_x p_y p_z\".split()[:dim]))).reshape(dim, 1)\n",
+    "    e0 = numpy.array(symbols(\" \".join(\"e0_x e0_y e0_z\".split()[:dim]))).reshape(dim, 1)\n",
+    "    e1 = numpy.array(symbols(\" \".join(\"e1_x e1_y e1_z\".split()[:dim]))).reshape(dim, 1)\n",
+    "    x = numpy.vstack([p, e0, e1]).reshape(-1)\n",
+    "        \n",
+    "    e = e1 - e0\n",
+    "    c = ((p - e0).T @ e) / (e.T @ e)\n",
+    "\n",
+    "    generators.append(CXXJacobianGenerator(c, x, f\"point_edge_closest_point_{dim}D_jacobian\"))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Edge-Edge"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ea0 = numpy.array(symbols(\"ea0_x ea0_y ea0_z\"))\n",
+    "ea1 = numpy.array(symbols(\"ea1_x ea1_y ea1_z\"))\n",
+    "eb0 = numpy.array(symbols(\"eb0_x eb0_y eb0_z\"))\n",
+    "eb1 = numpy.array(symbols(\"eb1_x eb1_y eb1_z\"))\n",
+    "x = numpy.vstack([ea0, ea1, eb0, eb1]).reshape(-1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "eb_to_ea = ea0 - eb0\n",
+    "ea = ea1 - ea0\n",
+    "eb = eb1 - eb0\n",
+    "\n",
+    "coefMtr = numpy.array([\n",
+    "    [ea.dot(ea), -eb.dot(ea)],\n",
+    "    [-eb.dot(ea), eb.dot(eb)]\n",
+    "])\n",
+    "\n",
+    "rhs = numpy.array([[-eb_to_ea.dot(ea)], [eb_to_ea.dot(eb)]])\n",
+    "c = Matrix(coefMtr).inv() @ rhs"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "generators.append(CXXJacobianGenerator(c, x, \"edge_edge_closest_point_jacobian\"))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Point-Triangle"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "p = numpy.array(symbols(\"p_x p_y p_z\")).reshape(3, 1)\n",
+    "t0 = numpy.array(symbols(\"t0_x t0_y t0_z\")).reshape(3, 1)\n",
+    "t1 = numpy.array(symbols(\"t1_x t1_y t1_z\")).reshape(3, 1)\n",
+    "t2 = numpy.array(symbols(\"t2_x t2_y t2_z\")).reshape(3, 1)\n",
+    "x = numpy.vstack([p, t0, t1, t2]).reshape(-1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle \\left[\\begin{matrix}\\left(- t_{0 x} + t_{1 x}\\right)^{2} + \\left(- t_{0 y} + t_{1 y}\\right)^{2} + \\left(- t_{0 z} + t_{1 z}\\right)^{2} & \\left(- t_{0 x} + t_{1 x}\\right) \\left(- t_{0 x} + t_{2 x}\\right) + \\left(- t_{0 y} + t_{1 y}\\right) \\left(- t_{0 y} + t_{2 y}\\right) + \\left(- t_{0 z} + t_{1 z}\\right) \\left(- t_{0 z} + t_{2 z}\\right)\\\\\\left(- t_{0 x} + t_{1 x}\\right) \\left(- t_{0 x} + t_{2 x}\\right) + \\left(- t_{0 y} + t_{1 y}\\right) \\left(- t_{0 y} + t_{2 y}\\right) + \\left(- t_{0 z} + t_{1 z}\\right) \\left(- t_{0 z} + t_{2 z}\\right) & \\left(- t_{0 x} + t_{2 x}\\right)^{2} + \\left(- t_{0 y} + t_{2 y}\\right)^{2} + \\left(- t_{0 z} + t_{2 z}\\right)^{2}\\end{matrix}\\right]$"
+      ],
+      "text/plain": [
+       "Matrix([\n",
+       "[                                    (-t0_x + t1_x)**2 + (-t0_y + t1_y)**2 + (-t0_z + t1_z)**2, (-t0_x + t1_x)*(-t0_x + t2_x) + (-t0_y + t1_y)*(-t0_y + t2_y) + (-t0_z + t1_z)*(-t0_z + t2_z)],\n",
+       "[(-t0_x + t1_x)*(-t0_x + t2_x) + (-t0_y + t1_y)*(-t0_y + t2_y) + (-t0_z + t1_z)*(-t0_z + t2_z),                                     (-t0_x + t2_x)**2 + (-t0_y + t2_y)**2 + (-t0_z + t2_z)**2]])"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle \\left[\\begin{matrix}\\left(p_{x} - t_{0 x}\\right) \\left(- t_{0 x} + t_{1 x}\\right) + \\left(p_{y} - t_{0 y}\\right) \\left(- t_{0 y} + t_{1 y}\\right) + \\left(p_{z} - t_{0 z}\\right) \\left(- t_{0 z} + t_{1 z}\\right)\\\\\\left(p_{x} - t_{0 x}\\right) \\left(- t_{0 x} + t_{2 x}\\right) + \\left(p_{y} - t_{0 y}\\right) \\left(- t_{0 y} + t_{2 y}\\right) + \\left(p_{z} - t_{0 z}\\right) \\left(- t_{0 z} + t_{2 z}\\right)\\end{matrix}\\right]$"
+      ],
+      "text/plain": [
+       "Matrix([\n",
+       "[(p_x - t0_x)*(-t0_x + t1_x) + (p_y - t0_y)*(-t0_y + t1_y) + (p_z - t0_z)*(-t0_z + t1_z)],\n",
+       "[(p_x - t0_x)*(-t0_x + t2_x) + (p_y - t0_y)*(-t0_y + t2_y) + (p_z - t0_z)*(-t0_z + t2_z)]])"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "basis = numpy.empty((2, 3), dtype=object)\n",
+    "basis[0] = (t1 - t0).T\n",
+    "basis[1] = (t2 - t0).T\n",
+    "A = basis @ basis.T\n",
+    "b = basis @ (p - t0)\n",
+    "display(Matrix(A))\n",
+    "display(Matrix(b))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# f.write(gen_jac_code(numpy.array(Matrix(A).diff(x)).reshape(24, 2), x, \n",
+    "#                      \"point_triangle_LSQ_matrix_jacobian\"))\n",
+    "# f.write(gen_jac_code(numpy.array(Matrix(b).diff(x)).reshape(12, 2).T, x, \n",
+    "#                      \"point_triangle_LSQ_RHS_jacobian\"))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "c = Matrix(A).inv() @ b"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "generators.append(CXXJacobianGenerator(c, x, \"point_triangle_closest_point_jacobian\"))"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Generate Code"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "generate_hpp_file(generators, \"closest_point_grad.hpp\")\n",
+    "generate_cpp_file(generators, \"closest_point_grad.cpp\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.1"
+  },
+  "orig_nbformat": 4,
+  "vscode": {
+   "interpreter": {
+    "hash": "5c7b89af1651d0b8571dde13640ecdccf7d5a6204171d6ab33e7c296e100e08a"
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/notebooks/distances.ipynb b/notebooks/distances.ipynb
new file mode 100644
index 000000000..8eebc12f5
--- /dev/null
+++ b/notebooks/distances.ipynb
@@ -0,0 +1,107 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import sympy\n",
+    "from generate_cpp_code import *"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def line_line_distance(ea0 : np.ndarray, ea1 : np.ndarray, eb0 : np.ndarray, eb1 : np.ndarray):\n",
+    "    normal = np.cross(ea1 - ea0, eb1 - eb0)\n",
+    "    line_to_line = (eb0 - ea0).dot(normal)\n",
+    "    return line_to_line * line_to_line / normal.dot(normal)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([ea0_x, ea0_y, ea0_z], dtype=object)"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "points = []\n",
+    "for i in 'ab':\n",
+    "    for j in range(2):\n",
+    "        points.append(np.array(sympy.symbols([f'e{i}{j}_{d}' for d in \"xyz\"])))\n",
+    "points = np.array(points).T\n",
+    "points[:, 0]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "funcs = [CXXFunctionGenerator(\n",
+    "    line_line_distance(points[:, 0], points[:, 1], points[:, 2], points[:, 3]), \n",
+    "    points.flatten(order=\"F\"), \"line_line_distance\")]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "generate_hpp_file(funcs, \"distances.hpp\")\n",
+    "generate_cpp_file(funcs, \"distances.cpp\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.1"
+  },
+  "orig_nbformat": 4,
+  "vscode": {
+   "interpreter": {
+    "hash": "5c7b89af1651d0b8571dde13640ecdccf7d5a6204171d6ab33e7c296e100e08a"
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/notebooks/generate_cpp_code.py b/notebooks/generate_cpp_code.py
new file mode 100644
index 000000000..1d5e47752
--- /dev/null
+++ b/notebooks/generate_cpp_code.py
@@ -0,0 +1,107 @@
+import numpy as np
+import sympy
+from sympy import Matrix, MatrixSymbol
+from sympy.printing import ccode
+import subprocess
+
+from utils import jacobian
+
+
+def generate_code(expr, out_var_name="J"):
+    CSE_results = sympy.cse(
+        expr, sympy.numbered_symbols("t"), optimizations='basic')
+    lines = []
+    for helper in CSE_results[0]:
+        if isinstance(helper[1], MatrixSymbol):
+            lines.append(f'const auto {helper[0]}[{helper[1].size}];')
+            lines.append(sympy.cxxcode(helper[1], helper[0]))
+        else:
+            lines.append(
+                f'const auto {sympy.cxxcode(helper[0])} = {sympy.cxxcode(helper[1])};')
+
+    if out_var_name != None:
+        for i, result in enumerate(CSE_results[1]):
+            lines.append(sympy.cxxcode(result, out_var_name))
+    else:
+        for i, result in enumerate(CSE_results[1]):
+            lines.append(f"return {sympy.cxxcode(result)};")
+
+    return '\n'.join(lines)
+
+
+class CXXFunctionGenerator:
+    def __init__(self, expr, params, name):
+        self.expr = expr
+        self.params = params
+        self.name = name
+        self.comment = ""
+        self.out_param = None
+
+    def signature(self):
+        params = ", ".join(f"double {ccode(var)}" for var in self.params)
+        if self.out_param != None:
+            params += f", double {self.out_param}"
+        ret_type = "double" if self.out_param == None else "void"
+        return f"{self.comment}\n{ret_type} {self.name}({params});"
+
+    def __call__(self):
+        signature = self.signature()[:-1]  # remove semicolon
+        return f"""
+{signature}{{
+{generate_code(self.expr, None if self.out_param ==
+               None else self.out_param.split("[")[0])}
+}}
+"""
+
+
+class CXXGradientGenerator(CXXFunctionGenerator):
+    def __init__(self, expr, params, name, out_param_name="grad"):
+        super().__init__(expr, params, name)
+        self.expr = self.expr.diff(Matrix(params))
+        self.out_param = f"{out_param_name}[{np.prod(self.expr.shape)}]"
+
+
+class CXXJacobianGenerator(CXXFunctionGenerator):
+    def __init__(self, expr, params, name, out_param_name="J"):
+        super().__init__(expr, params, name)
+        self.expr = jacobian(expr, params)
+        J_shape = self.expr.shape
+        self.expr = Matrix(np.array(self.expr).flatten(order="F"))
+        self.out_param = f"{out_param_name}[{np.prod(self.expr.shape)}]"
+        self.comment = f"// {out_param_name} is ({J_shape[0]}×{J_shape[1]}) flattened in column-major order"
+
+
+class CXXHessianGenerator(CXXJacobianGenerator):
+    def __init__(self, expr, params, name, out_param_name="hess"):
+        super().__init__(expr.diff(Matrix(params)), params, name, out_param_name)
+
+
+def generate_hpp_file(code_generators, file_name, transformer=None):
+    newline = "\n"
+    with open(file_name, 'w') as f:
+        f.write(f"""\
+#include <{file_name[:-4]}.hpp>
+
+namespace ipc{{
+namespace autogen{{
+    {newline.join(code_generator.signature()
+                  for code_generator in code_generators)}
+}}
+}}
+""")
+    subprocess.run(["clang-format", str(file_name), "-i"])
+
+
+def generate_cpp_file(code_generators, file_name):
+    newline = "\n"
+    with open(file_name, 'w') as f:
+        f.write(f"""\
+#include <{file_name[:-4]}.hpp>
+
+namespace ipc{{
+namespace autogen{{
+    {newline.join(code_generator() for code_generator in code_generators)}
+}}
+}}
+""")
+    subprocess.run(["clang-format", str(file_name), "-i"])
diff --git a/notebooks/tangent_basis_grad.ipynb b/notebooks/tangent_basis_grad.ipynb
new file mode 100644
index 000000000..3290a6580
--- /dev/null
+++ b/notebooks/tangent_basis_grad.ipynb
@@ -0,0 +1,370 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "35b911ac",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy\n",
+    "import sympy\n",
+    "from sympy import *\n",
+    "from sympy.matrices import *\n",
+    "from sympy.printing import ccode\n",
+    "\n",
+    "from generate_cpp_code import *\n",
+    "from utils import *"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4e2d285b",
+   "metadata": {},
+   "source": [
+    "## Point-Point"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "8f808648",
+   "metadata": {},
+   "source": [
+    "### 2D"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "0529cbcc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "p0 = numpy.array(sympy.symbols(\"p0_x p0_y\"))\n",
+    "p1 = numpy.array(sympy.symbols(\"p1_x p1_y\"))\n",
+    "x = numpy.concatenate([p0.T, p1.T])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "d5f30896",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "p0_to_p1 = normalize(p1 - p0)\n",
+    "\n",
+    "T = numpy.array([-p0_to_p1[1], p0_to_p1[0]]).reshape(-1, 1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "42db9ee8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "generators = [CXXJacobianGenerator(T, x, \"point_point_tangent_basis_2D_jacobian\")]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f601fd47",
+   "metadata": {},
+   "source": [
+    "### 3D"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "0f9af310",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "p0 = numpy.array(sympy.symbols(\"p0_x p0_y p0_z\"))\n",
+    "p1 = numpy.array(sympy.symbols(\"p1_x p1_y p1_z\"))\n",
+    "x = numpy.concatenate([p0.T, p1.T])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "c96391a1",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(6, 3, 2)"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "T = numpy.empty([3, 2], dtype=object)\n",
+    "\n",
+    "p0_to_p1 = p1 - p0\n",
+    "\n",
+    "cross_x = numpy.cross(numpy.array([1, 0, 0]), p0_to_p1)\n",
+    "cross_y = numpy.cross(numpy.array([0, 1, 0]), p0_to_p1)\n",
+    "\n",
+    "cross = numpy.array([\n",
+    "    sympy.Piecewise(\n",
+    "        (cross_x[i], sq_norm(cross_x) > sq_norm(cross_y)), \n",
+    "        (cross_y[i], True))\n",
+    "    for i in range(3)])\n",
+    "\n",
+    "T[:,0] = normalize(cross)\n",
+    "T[:,1] = normalize(numpy.cross(p0_to_p1, cross))\n",
+    "\n",
+    "sympy.Matrix(T).diff(x).shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "id": "ff0604c8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "generators.append(CXXJacobianGenerator(T, x, \"point_point_tangent_basis_3D_jacobian\"))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9c84b471",
+   "metadata": {},
+   "source": [
+    "## Point-Edge"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2e3c6fac",
+   "metadata": {},
+   "source": [
+    "### 2D"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "id": "55ce1ce3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "p = numpy.array(sympy.symbols(\"p_x p_y\"))\n",
+    "e0 = numpy.array(sympy.symbols(\"e0_x e0_y\"))\n",
+    "e1 = numpy.array(sympy.symbols(\"e1_x e1_y\"))\n",
+    "x = numpy.concatenate([p.T, e0.T, e1.T])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "id": "7b845c4a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "T = numpy.empty([2, 1], dtype=object)\n",
+    "\n",
+    "T[:, 0] = normalize(e1 - e0)\n",
+    "\n",
+    "generators.append(CXXJacobianGenerator(T, x, \"point_edge_tangent_basis_2D_jacobian\"))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "51c598b7",
+   "metadata": {},
+   "source": [
+    "### 3D"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "id": "4ff855d7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "p = numpy.array(sympy.symbols(\"p_x p_y p_z\"))\n",
+    "e0 = numpy.array(sympy.symbols(\"e0_x e0_y e0_z\"))\n",
+    "e1 = numpy.array(sympy.symbols(\"e1_x e1_y e1_z\"))\n",
+    "x = numpy.concatenate([p.T, e0.T, e1.T])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "id": "737aad1c",
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "T = numpy.empty([3, 2], dtype=object)\n",
+    "\n",
+    "e = e1 - e0\n",
+    "T[:, 0] = normalize(e)\n",
+    "T[:, 1] = normalize(numpy.cross(e, p - e0))\n",
+    "\n",
+    "generators.append(CXXJacobianGenerator(T, x, \"point_edge_tangent_basis_3D_jacobian\"))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "dc6e87ed",
+   "metadata": {},
+   "source": [
+    "## Edge-Edge"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "id": "6bcac693",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ea0 = numpy.array(sympy.symbols(\"ea0_x ea0_y ea0_z\"))\n",
+    "ea1 = numpy.array(sympy.symbols(\"ea1_x ea1_y ea1_z\"))\n",
+    "eb0 = numpy.array(sympy.symbols(\"eb0_x eb0_y eb0_z\"))\n",
+    "eb1 = numpy.array(sympy.symbols(\"eb1_x eb1_y eb1_z\"))\n",
+    "x = numpy.concatenate([ea0.T, ea1.T, eb0.T, eb1.T])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "id": "cd07cb2a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "T = numpy.empty([3, 2], dtype=object)\n",
+    "\n",
+    "ea = ea1 - ea0\n",
+    "T[:, 0] = normalize(ea)\n",
+    "\n",
+    "normal = numpy.cross(ea, eb1 - eb0);\n",
+    "T[:, 1] = normalize(numpy.cross(normal, ea));\n",
+    "\n",
+    "generators.append(CXXJacobianGenerator(T, x, \"edge_edge_tangent_basis_jacobian\"))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "bcd647a7",
+   "metadata": {},
+   "source": [
+    "## Point-Triangle"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "id": "3207482c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "p = numpy.array(sympy.symbols(\"p_x p_y p_z\"))\n",
+    "t0 = numpy.array(sympy.symbols(\"t0_x t0_y t0_z\"))\n",
+    "t1 = numpy.array(sympy.symbols(\"t1_x t1_y t1_z\"))\n",
+    "t2 = numpy.array(sympy.symbols(\"t2_x t2_y t2_z\"))\n",
+    "x = numpy.concatenate([p.T, t0.T, t1.T, t2.T])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "id": "9ee95a7a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "e0 = t1 - t0\n",
+    "normal = numpy.cross(e0, t2 - t0)\n",
+    "T = numpy.array([normalize(e0), normalize(numpy.cross(normal, e0))]).T"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "id": "0f6c947c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "generators.append(CXXJacobianGenerator(T, x, \"point_triangle_tangent_basis_jacobian\"))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "id": "9e95c122",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "generate_hpp_file(generators, \"tangent_basis_grad.hpp\")\n",
+    "generate_cpp_file(generators, \"tangent_basis_grad.cpp\")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.1"
+  },
+  "varInspector": {
+   "cols": {
+    "lenName": 16,
+    "lenType": 16,
+    "lenVar": 40
+   },
+   "kernels_config": {
+    "python": {
+     "delete_cmd_postfix": "",
+     "delete_cmd_prefix": "del ",
+     "library": "var_list.py",
+     "varRefreshCmd": "print(var_dic_list())"
+    },
+    "r": {
+     "delete_cmd_postfix": ") ",
+     "delete_cmd_prefix": "rm(",
+     "library": "var_list.r",
+     "varRefreshCmd": "cat(var_dic_list()) "
+    }
+   },
+   "types_to_exclude": [
+    "module",
+    "function",
+    "builtin_function_or_method",
+    "instance",
+    "_Feature"
+   ],
+   "window_display": false
+  },
+  "vscode": {
+   "interpreter": {
+    "hash": "d6ea818022cdff55a53271c8f07b73d2a935415fa70f5c4a81dac4a46ca819a5"
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/notebooks/utils.py b/notebooks/utils.py
new file mode 100644
index 000000000..42cc29077
--- /dev/null
+++ b/notebooks/utils.py
@@ -0,0 +1,26 @@
+import numpy as np
+import sympy
+
+
+def sq_norm(x):
+    """Compute the squared norm of a vector x."""
+    return x.dot(x)
+
+
+def norm(x):
+    """Compute the norm of a vector x."""
+    return sympy.sqrt(sq_norm(x))
+
+
+def normalize(x):
+    """Normalize a vector x."""
+    return x / norm(x)
+
+
+def jacobian(F, x):
+    J = np.empty((x.size * F.shape[0], F.shape[1]), dtype=object)
+    for xi in range(x.size):
+        for Fi in range(F.shape[0]):
+            for Fj in range(F.shape[1]):
+                J[xi * F.shape[0] + Fi, Fj] = F[Fi, Fj].diff(x[xi])
+    return J