diff --git a/README.md b/README.md
index 9f22cae..cd0d7be 100644
--- a/README.md
+++ b/README.md
@@ -212,13 +212,13 @@ in the terminal. Remote execution on Berkeley HPC (Savio) requires running
 sbatch saviojob.sh
 ```
 
-The data will be saved as VTK files. The file named as ```<prefix>_<electrode name>.vtk``` stores potential, field (RF only), and field square (RF only). The file named as ```<prefix>_<electrode name>_mesh.vtk``` stores the charge distribution, and the final mesh after adaptive mesh during simulation.
+The data will be saved as VTK files. The file named as ```<prefix>_<electrode_name>.vtk``` stores potential, field (RF only), and field square (RF only). The file named as ```<prefix>_<electrode_name>_mesh.vtk``` stores the charge distribution, and the final mesh after adaptive mesh during simulation.
 
 ### 4. Result analysis <div id='id-6'/> 
 
 Lastly, after we get the results, we can visualize the result by simply calling 
 ```
-Result.view(prefix, '<electrode name>')
+Result.view(prefix, '<electrode_name>')
 ```
 For example, if we want to visualize the field square of the SimpleTrap's RF electrode, we can simply call ```Result.view(prefix, 'RF')``` in the same folder as the result VTK files. Here is the visualization example
 
diff --git a/bem/area_constraints.py b/bem/area_constraints.py
index a97dbb0..45ad0ff 100644
--- a/bem/area_constraints.py
+++ b/bem/area_constraints.py
@@ -34,14 +34,12 @@ def __init__(self, start, end, radius, *a, **k):
 
     def lookup(self, x):
         # http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html
-        print((self.start-x).shape)
         d = np.linalg.norm(np.cross(self.start-x, self.end-x), axis=1)/np.linalg.norm(self.end-self.start)
         dmax = np.sqrt(np.linalg.norm(self.start-self.end)**2 + self.radius**2)
         dx1 = np.linalg.norm(self.start - x, axis=1)
         dx2 = np.linalg.norm(self.end - x, axis=1)
         # return np.where((d<self.radius), self.inside, self.outside)
         return np.where((d<self.radius)&(dx1<dmax)&(dx2<dmax), self.inside, self.outside)
-        # FIXME: cylinder has infinite length
 
 class Box(Constraint):
     def __init__(self, start=None, center=None, diagonal=None, end=None,
diff --git a/examples/SimpleTrap/SimpleTrap.ipynb b/examples/SimpleTrap/SimpleTrap.ipynb
index af15aac..be2dab8 100644
--- a/examples/SimpleTrap/SimpleTrap.ipynb
+++ b/examples/SimpleTrap/SimpleTrap.ipynb
@@ -19,15 +19,27 @@
     "import ipyparallel as ipp\n",
     "\n",
     "from bem import Electrodes, Sphere, Cylinder, Mesh, Grid, Configuration, Result\n",
-    "from bem.formats import stl"
+    "from bem.formats import stl\n",
+    "from bem.bemColors_lib.bemColors import bemColors\n",
+    "from helper_functions import helper_functions"
    ]
   },
   {
    "cell_type": "markdown",
-   "id": "d1fdc86f-b723-4269-a3ac-0499e4ded3dd",
+   "id": "840b9036-d3af-401a-b3f7-30ce6acbfc75",
    "metadata": {},
    "source": [
-    "# Import mesh from stl"
+    "# `I. Read STL and rename electrodes` \n",
+    "This part reads a colored model file `./SimpleTrap.stl`, and assigns a name for each color(thus for each electrode). The result is stored in `./SimpleTrap_mesh.vtk`.\n",
+    "\n",
+    "Here we name different electrodes and refine the meshes. The size of mesh deternmines the accuracy in next step. You should set the mesh finer and finer until the result converges. Also, you can set the mesh finer and test the convergence.\n",
+    "\n",
+    "## (1) color your electrodes (in fusion 360) and export stl\n",
+    "A trap always consists of different electrodes. Here we use different color to identify different electrodes(STL file stores color for each mesh). In fusion 360/Inventor, you can assign different appearance for different parts, of which we only cares about the color(because stl file only stores color), not texture. We recommend to use ../../bemCol_lib/bemCol.adsklib for different electrodes if you are usin Fusion 360. Same color for two electrodes means they are shorted. Then export to an stl file and move the file in the same folder as the mesh processing file. \n",
+    "\n",
+    "Below, \"prefix\" should be the name of the STL file. \n",
+    "\"a0\" should be close to the characteristic distance.\n",
+    "\"unit\" is the unit for the STL file."
    ]
   },
   {
@@ -38,7 +50,32 @@
    "outputs": [],
    "source": [
     "# base file name for outputs and inputs is the script name\n",
-    "prefix = \"SimpleTrap\""
+    "prefix = \"SimpleTrap\"\n",
+    "\n",
+    "a0 = 40e-6    # Distance from ion to electrode is 40 µm.\n",
+    "unit = 1e-6   # unit of the stl file is µm."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fd0c71e1-d7c1-47d4-8325-c26cb9210a7d",
+   "metadata": {},
+   "source": [
+    "Some notes:\n",
+    "* Why use standard color? We define a set of standard color lib `../../../bemCol_lib/bemCol.adsklib` for two reasons: 1. when exporting to .stl file, Fusion 360 will compress 0-256 color(RGBA32) into 0-32 color(RGBA16), which means that some similar colors in fusion 360 will become same color after export to .stl file. Our standard color lib avoid this problem. 2.If you use .stl file exported from Fusion 360 DIRECTLY, naming the electrodes would be very convinient. If you use other apps such as Inventor or Meshlab to generate .stl file, you can ignore the second reason because color encoding are quite different in different apps\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c4e6aaa4-6709-4e60-8332-d5b20237033f",
+   "metadata": {
+    "tags": []
+   },
+   "source": [
+    "## (2) The formal rename of electrode. \n",
+    "### (2.1) rename directly with color code\n",
+    "First read stl and drop all color codes. Be sure to put something random in rename argument to drop all color codes."
    ]
   },
   {
@@ -51,7 +88,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Import stl: /Users/nem0x/Documents/Pyckages/bem-savio/examples/SimpleTrap/SimpleTrap.stl \n",
+      "Import stl: C:\\Users\\3d_printing\\Documents\\Pyckages\\bem\\examples\\SimpleTrap\\SimpleTrap.stl \n",
       "\n",
       "Electrode colors (numbers):\n",
       "\n",
@@ -67,9 +104,6 @@
     }
    ],
    "source": [
-    "a0 = 40e-6    # Distance from ion to electrode is 40 µm.\n",
-    "unit = 1e-6   # unit of the stl file is µm.\n",
-    "\n",
     "# load electrode faces from colored stl\n",
     "# s_nta is intermediate processed stl file.\n",
     "s_nta = stl.read_stl(open(\"%s.stl\" % prefix, \"rb\"))\n",
@@ -81,10 +115,9 @@
   },
   {
    "cell_type": "markdown",
-   "id": "c4e6aaa4-6709-4e60-8332-d5b20237033f",
+   "id": "08a39274-32dc-4430-bec0-44ed2b79a19f",
    "metadata": {},
    "source": [
-    "## The formal rename of electrode. \n",
     "Assign each electrode a string name instead of its color coding. Use the numbers you get above.  \n",
     "`stl.stl_to_mesh()` prints normal vectors (different faces) in each electrode."
    ]
@@ -146,12 +179,98 @@
   },
   {
    "cell_type": "markdown",
-   "id": "ab1b7163-5174-4deb-9a82-892cad05efb6",
+   "id": "d06f63b6-ccb3-465b-8123-986709021233",
    "metadata": {},
    "source": [
-    "# Generate triangle mesh with constraints\n",
+    "### (2.2) rename with BemColors\n",
+    "\n",
+    "Below shows how to rename electrodes with `BemColors`\n",
     "\n",
-    "The meshes are 2-dimensional triangles on the surface of electrodes. The region enclosed by constraint shape can have finer mesh. Triangulation is done by `triangle` C library."
+    "The code is fully commented out, because the SimpleTrap.stl file doesn't use Fusion 360 colors. But The basic\n",
+    "\n",
+    "Firstly, we print all the colors appear in the stl file. If a color is of standard colors in `bemCol.adsklib` or named manually by function `bemCol.set_my_color()`, then the color will appear in a straightforward name, such as `bem0`,`bem1` ...... . Otherwise, the color will simply be named as `'_unkCol0'` , `'_unkCol1'`, ... \n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "id": "f5a63b75-7ec4-4e3b-97d1-b51b3f9fbb81",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# ele_col = bemColors(np.array(list(set(s_nta[2]))),('fusion360','export_stl'))\n",
+    "# ele_col.set_my_color(value = (178,178,178),cl_format = ('fusion360','export_stl','RGBA64'),name = 'self_defined')\n",
+    "# ele_col.print_stl_colors()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "871b225d-15b6-4b7a-9819-4d7be33afbdc",
+   "metadata": {},
+   "source": [
+    "Next, you need to assign a name for each color appeared above. With standard colors defined in `bemCol.adsklib`, the correspondence between color and electrode is clear, which make it easy for this step. You can also comment out some `ele_col.set_color_name()`, run all the codes again and observe the missing part in the printed figure, which corresponds to the electrode you comment out. The second method also serves as a double check. Besides, parameter in `stl.stl_to_mesh()` can be set as `quiet = False` to check whether the program reads planes correctly."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "id": "6d9f67f6-8743-4b01-b2f9-3ab39a88decf",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# # assign a name for each color\n",
+    "# ele_col.color_electrode(color = 'bem1',name = 'DC1')\n",
+    "# ele_col.color_electrode(color = 'bem2',name = 'DC2')\n",
+    "# ele_col.color_electrode(color = 'bem3',name = 'DC3')\n",
+    "# ele_col.color_electrode(color = 'bem4',name = 'DC4')\n",
+    "# ele_col.color_electrode(color = 'bem5',name = 'DC5')\n",
+    "# ele_col.color_electrode(color = 'bem6',name = 'DC6')\n",
+    "# ele_col.color_electrode(color = 'bem7',name = 'DC7')\n",
+    "# ele_col.color_electrode(color = 'bem8',name = 'DC8')\n",
+    "# ele_col.color_electrode(color = 'bem9',name = 'DC9')\n",
+    "# ele_col.color_electrode(color = 'bem10',name = 'DC10')\n",
+    "# ele_col.color_electrode(color = 'bem11',name = 'DC11')\n",
+    "# ele_col.color_electrode(color = 'bem12',name = 'DC12')\n",
+    "# ele_col.color_electrode(color = 'bem13',name = 'DC13')\n",
+    "# ele_col.color_electrode(color = 'bem14',name = 'DC14')\n",
+    "# ele_col.color_electrode(color = 'bem15',name = 'DC15')\n",
+    "# ele_col.color_electrode(color = 'bem16',name = 'DC16')\n",
+    "# ele_col.color_electrode(color = 'bem17',name = 'DC17')\n",
+    "# ele_col.color_electrode(color = 'bem18',name = 'DC18')\n",
+    "# ele_col.color_electrode(color = 'bem19',name = 'DC19')\n",
+    "# ele_col.color_electrode(color = 'bem20',name = 'DC20')\n",
+    "# ele_col.color_electrode(color = 'bem21',name = 'DC21')\n",
+    "# ele_col.color_electrode(color = 'bem30',name = 'DC0')\n",
+    "# ele_col.color_electrode(color = 'bem25',name = 'RF')\n",
+    "\n",
+    "# # print colors still with no name. These meshes will be neglected in the following codes\n",
+    "# ele_col.drop_colors()\n",
+    "\n",
+    "# # read stl into mesh with electrode names\n",
+    "# # unnamed meshes will not be imported at all\n",
+    "# mesh = Mesh.from_mesh(stl.stl_to_mesh(*s_nta, scale=a0/unit,\n",
+    "#     rename=ele_col.electrode_colors, quiet=True))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ab1b7163-5174-4deb-9a82-892cad05efb6",
+   "metadata": {
+    "tags": []
+   },
+   "source": [
+    "# `II. Remesh with constraint` \n",
+    "\n",
+    "In this step, we generate triangle mesh with constraints. The meshes are 2-dimensional triangles on the surface of electrodes. The region enclosed by constraint shape can have finer mesh. Triangulation is done by `triangle` C library. Folowing variables are all in unit `mesh_unit` now.\n",
+    "\n",
+    "Our remesh strategy consists of two steps of triangulation: \n",
+    "1. global triangulation without constraint. This step eliminate some long and sharp triangles by combining and dividing, and obtains a coarse grain triangulated model.\n",
+    "2. local triangulation with constraint. This step refines each triangles in step 1, the triangle density is defined by `mesh.areas_from_constraints`\n",
+    "\n",
+    "It is better to update triangle density gradually to accurately position the constraint boundary. Parameters in the below code block should be tuned specificly for different trap geometries.\n",
+    "\n",
+    "Notes:\n",
+    "- mesh.triangulate function takes in a parameter opts, after some processing of the wrap function, opts changes a little and prints as final opts. These final opts finally passes into Fastlap triangulate package, of which parameters and usage can be found in this link (http://www.cs.cmu.edu/~quake/triangle.research.html)."
    ]
   },
   {
@@ -164,6 +283,42 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
+      "start triangulate DC1\n",
+      "('final opts', 'q20Qzr')\n",
+      "finish triangulate DC1\n",
+      "start triangulate DC3\n",
+      "('final opts', 'q20Qzr')\n",
+      "finish triangulate DC3\n",
+      "start triangulate DC5\n",
+      "('final opts', 'q20Qzr')\n",
+      "finish triangulate DC5\n",
+      "start triangulate DC2\n",
+      "('final opts', 'q20Qzr')\n",
+      "finish triangulate DC2\n",
+      "start triangulate RF\n",
+      "('final opts', 'q20Qzr')\n",
+      "finish triangulate RF\n",
+      "start triangulate DC4\n",
+      "('final opts', 'q20Qzr')\n",
+      "finish triangulate DC4\n",
+      "start triangulate DC1\n",
+      "('final opts', 'q20Qzra')\n",
+      "finish triangulate DC1\n",
+      "start triangulate DC3\n",
+      "('final opts', 'q20Qzra')\n",
+      "finish triangulate DC3\n",
+      "start triangulate DC5\n",
+      "('final opts', 'q20Qzra')\n",
+      "finish triangulate DC5\n",
+      "start triangulate DC2\n",
+      "('final opts', 'q20Qzra')\n",
+      "finish triangulate DC2\n",
+      "start triangulate RF\n",
+      "('final opts', 'q20Qzra')\n",
+      "finish triangulate RF\n",
+      "start triangulate DC4\n",
+      "('final opts', 'q20Qzra')\n",
+      "finish triangulate DC4\n",
       "start triangulate DC1\n",
       "('final opts', 'q20Qzra')\n",
       "finish triangulate DC1\n",
@@ -182,18 +337,30 @@
       "start triangulate DC4\n",
       "('final opts', 'q20Qzra')\n",
       "finish triangulate DC4\n",
-      "Output vtk: /Users/nem0x/Documents/Pyckages/bem-savio/examples/SimpleTrap/SimpleTrap.vtk\n"
+      "Output vtk: C:\\Users\\3d_printing\\Documents\\Pyckages\\bem\\examples\\SimpleTrap\\SimpleTrap.vtk\n"
      ]
     }
    ],
    "source": [
-    "# set .1 max area within 3\n",
-    "# areas_from_constraints specifies sphere with finer mesh inside it.\n",
-    "mesh.areas_from_constraints(Sphere(center=np.array([0, 0, 1.]),\n",
-    "           radius=2, inside=0.2, outside=10))    # \"inside\", \"outside\" set different mesh densities.\n",
-    "# retriangulate quality and quiet with areas\n",
-    "mesh.triangulate(opts=\"q20Q\", new=False)\n",
-    "# save base mesh to vtk\n",
+    "# 1. Global triangulation\n",
+    "mesh.triangulate(opts=\"q20Q\",new = False)\n",
+    "\n",
+    "# 2. Local triangulation\n",
+    "# areas_from_constraints specifies sphere/cylinder with finer mesh inside it.\n",
+    "# \"inside\", \"outside\" set different mesh densities.\n",
+    "\n",
+    "# first refine\n",
+    "mesh.areas_from_constraints(Sphere(center=np.array([0, 0, 1.]), radius=2, inside=2., outside=10))    # Sphere constraint\n",
+    "# mesh.areas_from_constraints(Cylinder(start = np.array([-1, 0, 0]), end = np.array([1, 0, 0]), radius = 2, inside=2, outside=10)) # Cyliner constraint\n",
+    "mesh.triangulate(opts=\"q20Q\", new=False) # retriangulate quality and quiet with areas\n",
+    "\n",
+    "# second refine\n",
+    "mesh.areas_from_constraints(Sphere(center=np.array([0, 0, 1.]), radius=2, inside=0.2, outside=10))    # Sphere constraint\n",
+    "# mesh.areas_from_constraints(Cylinder(start = np.array([-1, 0, 0]), end = np.array([1, 0, 0]), radius = 2, inside=0.2, outside=10)) # Cyliner constraint\n",
+    "mesh.triangulate(opts=\"q20Q\", new=False) # retriangulate quality and quiet with areas\n",
+    "\n",
+    "\n",
+    "# 3. save base mesh to vtk\n",
     "mesh.to_vtk(prefix)\n",
     "print(\"Output vtk:\",os.path.abspath(\"./\"+prefix+\".vtk\"))    # output path"
    ]
@@ -203,7 +370,8 @@
    "id": "d2c3c568-f382-4373-85b2-c125c77f3497",
    "metadata": {},
    "source": [
-    "## Plot mesh\n",
+    "## Visualize in 3D to check mesh\n",
+    "Can be used to check whether electrode faces are colored properly, and whether electrodes are named properly.\n",
     "\n",
     "Plot mesh using pyvista plotter in a separate window.\n",
     "On macOS perhaps Linux as well, the window doesn't close properly after clicking close button, but the code can continue to run."
@@ -225,7 +393,8 @@
    "id": "a48cf7ce-c550-4249-a660-4cdeff01916c",
    "metadata": {},
    "source": [
-    "# Main boundary element calculations\n",
+    "# `III. Main boundary element calculations`\n",
+    "## (1) Define job function\n",
     "\n",
     "In `run_job` function, `job` is `Configuration` instance and `grid` is discretirized spatial grid (not the mesh). The general workflow (also the routine of BEM method) are:  \n",
     "1. `solve_singularities()` solves charge distributions by iterative methods to make it consistent with one electrode at 1V and others at 0V (unit potentials). `adapt_mesh()` refines meshes adaptively to achieve certain precision while solving sigulartities.\n",
@@ -241,7 +410,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 7,
    "id": "8d51a66b-90fd-42ce-8860-380af9584859",
    "metadata": {},
    "outputs": [],
@@ -267,11 +436,9 @@
    "id": "a48439af-df27-4640-b860-f11d959c5f79",
    "metadata": {},
    "source": [
-    "## Define grid\n",
-    "\n",
-    "Create a grid in unit of scaled length `l`. Only choose the interested region (trap center) to save time.\n",
+    "## (2) Define grid\n",
     "\n",
-    "For reference, to compute Seidelin trap, grid shape = (60, 60, 60) takes 266 s, while shape = (150, 150, 150) takes 3369 s."
+    "Create a grid in unit of scaled length `l`. Only choose the interested region (trap center) to save time."
    ]
   },
   {
@@ -311,12 +478,12 @@
    "id": "c922dff7-6961-4d06-ba4c-915b868fe3ee",
    "metadata": {},
    "source": [
-    "# Parallel computation using ipyparallel, which is compatible with Savio BRC-HPC"
+    "## (3) Parallel computation using ipyparallel, which is compatible with Savio BRC-HPC"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 9,
    "id": "6cabdbfb-a5f6-4084-aa41-7b6ed6567e80",
    "metadata": {},
    "outputs": [
@@ -324,25 +491,25 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Starting 8 engines with <class 'ipyparallel.cluster.launcher.LocalEngineSetLauncher'>\n"
+      "Starting 16 engines with <class 'ipyparallel.cluster.launcher.LocalEngineSetLauncher'>\n"
      ]
     },
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "INFO:ipyparallel.cluster.cluster.1666855509-orme:Starting 8 engines with <class 'ipyparallel.cluster.launcher.LocalEngineSetLauncher'>\n"
+      "INFO:ipyparallel.cluster.cluster.1666995915-liwq:Starting 16 engines with <class 'ipyparallel.cluster.launcher.LocalEngineSetLauncher'>\n"
      ]
     },
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "5a24adaf80f940d6ab0c1bbd6ee51211",
+       "model_id": "e2c2a6cecbd8427c9f0c9746e6608f4d",
        "version_major": 2,
        "version_minor": 0
       },
       "text/plain": [
-       "  0%|          | 0/8 [00:00<?, ?engine/s]"
+       "  0%|          | 0/16 [00:00<?, ?engine/s]"
       ]
      },
      "metadata": {},
@@ -351,7 +518,7 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "d640f1738be64727838bd2cacf50278c",
+       "model_id": "d37e3776c45d49408402aa8783674283",
        "version_major": 2,
        "version_minor": 0
       },
@@ -366,13 +533,14 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Computing time: 11.055103 s\n"
+      "Computing time: 9.654386 s\n"
      ]
     }
    ],
    "source": [
     "# create job list\n",
     "jobs = list(Configuration.select(mesh,\"DC.*\",\"RF\"))\n",
+    "# evaluate electric potential for each configurations: one electrode at 1V, the rest 0V.\n",
     "\n",
     "# parallel computation\n",
     "mycluster = ipp.Cluster()\n",
@@ -395,12 +563,16 @@
    "id": "ae146844-83bb-4b22-941f-9417dd671d35",
    "metadata": {},
    "source": [
-    "# View result in 3D"
+    "## (4) View result in 3D\n",
+    "\n",
+    "This part can be done in other files as well. You only need to define prefix and make sure the two result files for each electrode, i.e. ```<prefix>_<electrode_name>.vtk``` and ```<prefix>_<electrode_name>_mesh.vtk```, are in the same directory as the code. Also don't remember to put the base mesh ```<prefix>_mesh.vtk``` in the same folder.\n",
+    "\n",
+    "Result view inlcudes charge distribution, potential/field_square contours."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": 10,
    "id": "8b0925f3-61b3-42b9-8028-c84abb6cb33d",
    "metadata": {},
    "outputs": [],
@@ -419,7 +591,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 11,
    "id": "c8e64059-24dc-4a9c-8d50-27b4f388ca88",
    "metadata": {},
    "outputs": [],
@@ -437,7 +609,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 12,
    "id": "9f837b51-4b0f-4ba0-b6e5-9b894c5ccadc",
    "metadata": {},
    "outputs": [],
@@ -450,17 +622,118 @@
    "id": "9275120d-072b-46da-b8c6-76043a1cf716",
    "metadata": {},
    "source": [
-    "# Load result"
+    "# `IV. Load result for analysis`\n",
+    "\n",
+    "All analysis are intended to be done outside of the simulation file. Here shows how to access the simulated data, such as potential (DC and RF), field (RF only), field square (RF only), for each electrode."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 13,
    "id": "4da591e3-db45-4090-98a3-25ce564e2034",
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[[0.29602049, 0.20908952, 0.15270445, ..., 0.0081563 ,\n",
+       "         0.0076957 , 0.00727314],\n",
+       "        [0.32986564, 0.22669958, 0.16191063, ..., 0.00824515,\n",
+       "         0.00779002, 0.00737097],\n",
+       "        [0.35066552, 0.23678475, 0.16634174, ..., 0.0083209 ,\n",
+       "         0.00787349, 0.0074599 ],\n",
+       "        ...,\n",
+       "        [0.33296267, 0.22118485, 0.15304576, ..., 0.01101531,\n",
+       "         0.01051962, 0.01004972],\n",
+       "        [0.33107613, 0.22497251, 0.15900469, ..., 0.01126803,\n",
+       "         0.01074391, 0.01025011],\n",
+       "        [0.31920493, 0.22244071, 0.16108435, ..., 0.01151763,\n",
+       "         0.01096633, 0.01044903]],\n",
+       "\n",
+       "       [[0.29501693, 0.20868764, 0.15255272, ..., 0.00815696,\n",
+       "         0.0076965 , 0.00727403],\n",
+       "        [0.32900035, 0.22638087, 0.16180271, ..., 0.00824572,\n",
+       "         0.00779075, 0.00737185],\n",
+       "        [0.35000655, 0.23655553, 0.16628435, ..., 0.00832137,\n",
+       "         0.00787412, 0.00746071],\n",
+       "        ...,\n",
+       "        [0.3331624 , 0.22127416, 0.15309373, ..., 0.0110143 ,\n",
+       "         0.01051858, 0.01004873],\n",
+       "        [0.33125041, 0.22505041, 0.15903744, ..., 0.0112667 ,\n",
+       "         0.0107429 , 0.01024916],\n",
+       "        [0.31427077, 0.22021981, 0.15966461, ..., 0.01151651,\n",
+       "         0.01096513, 0.01044803]],\n",
+       "\n",
+       "       [[0.29407397, 0.20839777, 0.15247447, ..., 0.00815764,\n",
+       "         0.00769728, 0.00727487],\n",
+       "        [0.32819762, 0.22616496, 0.16176069, ..., 0.00824634,\n",
+       "         0.00779146, 0.00737263],\n",
+       "        [0.34938998, 0.23641893, 0.16627903, ..., 0.00832193,\n",
+       "         0.00787476, 0.0074614 ],\n",
+       "        ...,\n",
+       "        [0.33332176, 0.22137719, 0.15314227, ..., 0.01101316,\n",
+       "         0.01051769, 0.01004783],\n",
+       "        [0.41564396, 0.25487628, 0.16711187, ..., 0.01126544,\n",
+       "         0.01074176, 0.01024829],\n",
+       "        [0.31438244, 0.22028103, 0.15969845, ..., 0.01151513,\n",
+       "         0.01096388, 0.01044699]],\n",
+       "\n",
+       "       ...,\n",
+       "\n",
+       "       [[0.29440461, 0.20857692, 0.15258303, ..., 0.00815963,\n",
+       "         0.00769907, 0.00727646],\n",
+       "        [0.32833701, 0.22624301, 0.16181158, ..., 0.00824828,\n",
+       "         0.00779298, 0.007374  ],\n",
+       "        [0.34931039, 0.23637803, 0.16626518, ..., 0.00832358,\n",
+       "         0.00787623, 0.00746282],\n",
+       "        ...,\n",
+       "        [0.33392138, 0.22171341, 0.15329234, ..., 0.01100632,\n",
+       "         0.0105118 , 0.01004296],\n",
+       "        [0.44516375, 0.28576752, 0.19368984, ..., 0.01125869,\n",
+       "         0.01073611, 0.0102435 ],\n",
+       "        [0.31578175, 0.22088898, 0.15995597, ..., 0.01150834,\n",
+       "         0.01095829, 0.01044228]],\n",
+       "\n",
+       "       [[0.29569508, 0.20907355, 0.15278628, ..., 0.00815929,\n",
+       "         0.00769859, 0.00727585],\n",
+       "        [0.32947095, 0.22667806, 0.16198554, ..., 0.00824799,\n",
+       "         0.00779259, 0.00737343],\n",
+       "        [0.35019621, 0.23671306, 0.16639166, ..., 0.00832335,\n",
+       "         0.00787583, 0.00746237],\n",
+       "        ...,\n",
+       "        [0.33379646, 0.22163655, 0.15325252, ..., 0.01100715,\n",
+       "         0.01051263, 0.01004353],\n",
+       "        [0.33234261, 0.22555919, 0.15924789, ..., 0.01125968,\n",
+       "         0.01073693, 0.01024417],\n",
+       "        [0.31573034, 0.22084966, 0.15993239, ..., 0.01150997,\n",
+       "         0.01095922, 0.01044303]],\n",
+       "\n",
+       "       [[0.29705451, 0.20967617, 0.15306292, ..., 0.00815894,\n",
+       "         0.00769799, 0.00727509],\n",
+       "        [0.33070113, 0.22720135, 0.16222113, ..., 0.00824761,\n",
+       "         0.00779215, 0.0073729 ],\n",
+       "        [0.35118854, 0.23713188, 0.16657457, ..., 0.00832315,\n",
+       "         0.00787545, 0.0074617 ],\n",
+       "        ...,\n",
+       "        [0.33356024, 0.22153048, 0.15318983, ..., 0.01100786,\n",
+       "         0.01051346, 0.01004436],\n",
+       "        [0.40250648, 0.26222798, 0.17994094, ..., 0.01126098,\n",
+       "         0.01073781, 0.01024487],\n",
+       "        [0.31553337, 0.22076307, 0.15989186, ..., 0.01151104,\n",
+       "         0.01096048, 0.0104439 ]]])"
+      ]
+     },
+     "execution_count": 13,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
    "source": [
-    "rf_result = Result.from_vtk(prefix, 'RF')"
+    "rf_result = Result.from_vtk(prefix, 'RF')\n",
+    "rf_result.grid.to_xyz() # access the (x, y, z) coordinates\n",
+    "rf_result.potential # access the simulated potential data\n",
+    "rf_result.field # access the simulated field data\n",
+    "rf_result.field_square # access the simulated field_square data"
    ]
   },
   {
diff --git a/examples/SimpleTrap/SimpleTrap_DC1.vtk b/examples/SimpleTrap/SimpleTrap_DC1.vtk
index e5d6acb..14b9c58 100644
Binary files a/examples/SimpleTrap/SimpleTrap_DC1.vtk and b/examples/SimpleTrap/SimpleTrap_DC1.vtk differ
diff --git a/examples/SimpleTrap/SimpleTrap_DC1_mesh.vtk b/examples/SimpleTrap/SimpleTrap_DC1_mesh.vtk
index e1ffcc5..211c62f 100644
Binary files a/examples/SimpleTrap/SimpleTrap_DC1_mesh.vtk and b/examples/SimpleTrap/SimpleTrap_DC1_mesh.vtk differ
diff --git a/examples/SimpleTrap/SimpleTrap_DC2.vtk b/examples/SimpleTrap/SimpleTrap_DC2.vtk
index 1a5dbab..2b80082 100644
Binary files a/examples/SimpleTrap/SimpleTrap_DC2.vtk and b/examples/SimpleTrap/SimpleTrap_DC2.vtk differ
diff --git a/examples/SimpleTrap/SimpleTrap_DC2_mesh.vtk b/examples/SimpleTrap/SimpleTrap_DC2_mesh.vtk
index 9123ae9..444391a 100644
Binary files a/examples/SimpleTrap/SimpleTrap_DC2_mesh.vtk and b/examples/SimpleTrap/SimpleTrap_DC2_mesh.vtk differ
diff --git a/examples/SimpleTrap/SimpleTrap_DC3.vtk b/examples/SimpleTrap/SimpleTrap_DC3.vtk
index 1d792d3..3b4b9f7 100644
Binary files a/examples/SimpleTrap/SimpleTrap_DC3.vtk and b/examples/SimpleTrap/SimpleTrap_DC3.vtk differ
diff --git a/examples/SimpleTrap/SimpleTrap_DC3_mesh.vtk b/examples/SimpleTrap/SimpleTrap_DC3_mesh.vtk
index be6b69f..16050f5 100644
Binary files a/examples/SimpleTrap/SimpleTrap_DC3_mesh.vtk and b/examples/SimpleTrap/SimpleTrap_DC3_mesh.vtk differ
diff --git a/examples/SimpleTrap/SimpleTrap_DC4.vtk b/examples/SimpleTrap/SimpleTrap_DC4.vtk
index 5b56ae4..9060919 100644
Binary files a/examples/SimpleTrap/SimpleTrap_DC4.vtk and b/examples/SimpleTrap/SimpleTrap_DC4.vtk differ
diff --git a/examples/SimpleTrap/SimpleTrap_DC4_mesh.vtk b/examples/SimpleTrap/SimpleTrap_DC4_mesh.vtk
index 14d221e..50f1f80 100644
Binary files a/examples/SimpleTrap/SimpleTrap_DC4_mesh.vtk and b/examples/SimpleTrap/SimpleTrap_DC4_mesh.vtk differ
diff --git a/examples/SimpleTrap/SimpleTrap_DC5.vtk b/examples/SimpleTrap/SimpleTrap_DC5.vtk
index c5e1478..61fc137 100644
Binary files a/examples/SimpleTrap/SimpleTrap_DC5.vtk and b/examples/SimpleTrap/SimpleTrap_DC5.vtk differ
diff --git a/examples/SimpleTrap/SimpleTrap_DC5_mesh.vtk b/examples/SimpleTrap/SimpleTrap_DC5_mesh.vtk
index 998132c..605fffc 100644
Binary files a/examples/SimpleTrap/SimpleTrap_DC5_mesh.vtk and b/examples/SimpleTrap/SimpleTrap_DC5_mesh.vtk differ
diff --git a/examples/SimpleTrap/SimpleTrap_RF.vtk b/examples/SimpleTrap/SimpleTrap_RF.vtk
index 9c09a53..04761e7 100644
Binary files a/examples/SimpleTrap/SimpleTrap_RF.vtk and b/examples/SimpleTrap/SimpleTrap_RF.vtk differ
diff --git a/examples/SimpleTrap/SimpleTrap_RF_mesh.vtk b/examples/SimpleTrap/SimpleTrap_RF_mesh.vtk
index 2607fc9..809112c 100644
Binary files a/examples/SimpleTrap/SimpleTrap_RF_mesh.vtk and b/examples/SimpleTrap/SimpleTrap_RF_mesh.vtk differ
diff --git a/examples/SimpleTrap/SimpleTrap_mesh.vtk b/examples/SimpleTrap/SimpleTrap_mesh.vtk
index 6f1f33c..4ded886 100644
Binary files a/examples/SimpleTrap/SimpleTrap_mesh.vtk and b/examples/SimpleTrap/SimpleTrap_mesh.vtk differ