diff --git a/src/bvpy/domains/abstract.py b/src/bvpy/domains/abstract.py
index c284b621833bd94f30fba53336ff3ab75b3ccf3e..2aa00aca77bc40f636248e39be97f2640e9218c0 100644
--- a/src/bvpy/domains/abstract.py
+++ b/src/bvpy/domains/abstract.py
@@ -391,6 +391,11 @@ class AbstractDomain(ABC):
         if self.mesh is not None:
             self.model.mesh.clear()
 
+        if not self.volumes:
+            self.model.addPhysicalGroup(dim, [self.surfaces[0]])
+        else:
+            self.model.addPhysicalGroup(dim, [self.volumes[0]])
+            
         self.model.mesh.generate(dim)
         self.mesh, cell_markers, facet_markers = gmshio.model_to_mesh(self.model, MPI.COMM_SELF, 0)
         # self.bdata = _generate_default_boundary_data(self.mesh)
diff --git a/src/bvpy/domains/primitives.py b/src/bvpy/domains/primitives.py
index 52e751c225627247c6ae8a91b068a6d161fb3d5b..b8c73901870f24aa1d0beb202b4b6a2c9e393e64 100644
--- a/src/bvpy/domains/primitives.py
+++ b/src/bvpy/domains/primitives.py
@@ -150,10 +150,8 @@ class Rectangle(AbstractDomain, OccModel):
         None
 
         """
-        _dim = 2
         self.surfaces[0] = self.factory.addRectangle(x, y, 0, length, width)
         self.synchronize()
-        self.model.addPhysicalGroup(_dim, [self.surfaces[0]])
 
 class Disk(AbstractDomain, OccModel):
     """Instanciates a circular domain.
@@ -277,14 +275,10 @@ class Disk(AbstractDomain, OccModel):
         None
 
         """
-        _dim = 2
-        _tag = [1]
-
         circle = self.factory.addCircle(*center, radius)
         curve_loop = self.factory.addCurveLoop([circle], circle)
         self.surfaces[0] = self.factory.addPlaneSurface([curve_loop])
         self.synchronize()
-        self.model.addPhysicalGroup(_dim, [self.surfaces[0]])
 
 
 class HemiSphere(AbstractDomain, BuiltInModel):
@@ -486,6 +480,7 @@ class Cube(AbstractDomain, OccModel):
     def __init__(self, x=0, y=0, z=0, dx=1, dy=1, dz=1, **kwargs):
         super(Cube, self).__init__(**kwargs)
         self._characteristic_size = {}
+        self._cell_type = 'tetra'
         self.geometry(x, y, z, dx, dy, dz)
 
     def geometry(self, x, y, z, dx, dy, dz):
@@ -498,6 +493,7 @@ class Cylinder(AbstractDomain, OccModel):
                  dx=0, dy=0, dz=1, r=0.5, open=False, **kwargs):
         super(Cylinder, self).__init__(**kwargs)
         self._characteristic_size = {}
+        self._cell_type = 'tetra'
         self.geometry(x, y, z, dx, dy, dz, r, open)
 
     def geometry(self, x, y, z, dx, dy, dz, r, open):
@@ -604,6 +600,7 @@ class Sphere(AbstractDomain, OccModel):
         """
         super(Sphere, self).__init__(**kwargs)
         self._characteristic_size = {'Radius': radius}
+        self._cell_type = 'tetra'
         self.geometry(radius, center)
 
     def geometry(self, radius, center):
@@ -622,17 +619,9 @@ class Sphere(AbstractDomain, OccModel):
         None
 
         """
-        # self.volumes[0] = self.factory.addSphere(*center, radius)
-        # self.synchronize()
-
-        _dim = 2
-        _tag = [1]
-        self.surfaces[0] = self.factory.addSphere(*center, radius)
+        self.volumes[0] = self.factory.addSphere(*center, radius)
         self.synchronize()
-        self.model.addPhysicalGroup(_dim, [self.surfaces[0]])
-        # self.model.mesh.generate(_dim)
-
-
+        
 class Ellipsoid(AbstractDomain, OccModel):
     """Instanciates an Ellispoid domain.
 
@@ -731,6 +720,7 @@ class Ellipsoid(AbstractDomain, OccModel):
         super(Ellipsoid, self).__init__(**kwargs)
         self._characteristic_size = {'Oxy radius': radius_xy,
                                      'Oz radius': radius_z}
+        self._cell_type = 'tetra'
         self.geometry(radius_xy, radius_z, center)
 
     def geometry(self, radius_xy, radius_z, center=[0, 0, 0]):
@@ -756,7 +746,6 @@ class Ellipsoid(AbstractDomain, OccModel):
         self.volumes[0] = self.factory.addSphere(*center, radius_z)
         self.factory.dilate([(3, self.volumes[0])], *center, ratio, ratio, 1)
         self.synchronize()
-        self.model.addPhysicalGroup(2, [self.volumes[0]])
 
 
 class Torus(AbstractDomain, OccModel):
@@ -883,4 +872,3 @@ class Torus(AbstractDomain, OccModel):
         """
         self.surfaces[0] = self.factory.addTorus(*center, main_radius, tube_radius)
         self.synchronize()
-        self.model.addPhysicalGroup(2, [self.surfaces[0]])
diff --git a/tutorials/bvpy_tutorial_2_domains.ipynb b/tutorials/bvpy_tutorial_2_domains.ipynb
index e99e1299a5acce48012b6257b5903b32d1843d34..21223f26862f5f82ca3222de60c257dade1f0782 100644
--- a/tutorials/bvpy_tutorial_2_domains.ipynb
+++ b/tutorials/bvpy_tutorial_2_domains.ipynb
@@ -46,15 +46,7 @@
    "cell_type": "code",
    "execution_count": 1,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "####### Rectangle\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "from bvpy.domains import Rectangle\n",
     "domain = Rectangle()"
@@ -77,7 +69,6 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "####### Rectangle 2\n",
       "Number of triangles within the domain AFTER discretization: 4\n"
      ]
     }
@@ -153,7 +144,7 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "72bb02fedfcd464faa86c5f07de28589",
+       "model_id": "962f1bb475cb45eca0a1228b7b6c41ae",
        "version_major": 2,
        "version_minor": 0
       },
@@ -191,21 +182,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 5,
    "metadata": {},
    "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Info    : Clearing all models and views...\n",
-      "Info    : Done clearing all models and views\n"
-     ]
-    },
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "7face6d24f314334a03f15805958abf8",
+       "model_id": "e561ccff5399482eb3d4928de75f32ee",
        "version_major": 2,
        "version_minor": 0
       },
@@ -226,8 +209,7 @@
     "\n",
     "d1 = Rectangle(length=2, width=3, clear=True)\n",
     "\n",
-    "# TODO methode pas propre ??\n",
-    "d1.set_cell_size(.1)\n",
+    "d1.set_cell_size(.05)\n",
     "\n",
     "plot(d1)\n"
    ]
@@ -249,13 +231,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 6,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "4a3a113f374547749ae6a5a01cdeff4c",
+       "model_id": "b5db0aa22c1945e0b82b57bafae8c249",
        "version_major": 2,
        "version_minor": 0
       },
@@ -281,13 +263,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 7,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "298b8023f3014b1f811f3a36b17a7765",
+       "model_id": "3d3f14eed8fa4bf4b371dd7395cd9005",
        "version_major": 2,
        "version_minor": 0
       },
@@ -314,7 +296,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 1,
    "metadata": {},
    "outputs": [
     {
@@ -336,6 +318,8 @@
    ],
    "source": [
     "from bvpy.domains import HemiSphere\n",
+    "from bvpy.utils.visu import plot, set_renderer\n",
+    "\n",
     "hsph = HemiSphere(cell_size=0.1, clear=True)\n",
     "plot(hsph)"
    ]
@@ -349,21 +333,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 9,
    "metadata": {},
    "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Info    : Clearing all models and views...\n",
-      "Info    : Done clearing all models and views\n"
-     ]
-    },
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "1013193846f44b1787cc7f3d66354a06",
+       "model_id": "f4c9b9021ef04184b2a14e75c4f45d32",
        "version_major": 2,
        "version_minor": 0
       },
@@ -391,13 +367,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 10,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "f2c103eabca941508ed07a529681d6c5",
+       "model_id": "f282dc8ea0864d71b137db207b31f06f",
        "version_major": 2,
        "version_minor": 0
       },
@@ -424,13 +400,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 11,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "67b2c90594e34945a4e117801cd1ae51",
+       "model_id": "47d8a09373f244449677019fb66eafd7",
        "version_major": 2,
        "version_minor": 0
       },
@@ -471,13 +447,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 12,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "2984515db25541c1b1888b3a3a5ab4fe",
+       "model_id": "890f8e1bdad14aae8149deda32465d1d",
        "version_major": 2,
        "version_minor": 0
       },
@@ -497,14 +473,12 @@
     " \n",
     "set_renderer('notebook')\n",
     "\n",
-    "#TODO : how these operations are done\n",
     "dsk1 = Disk(center=[0, -.25, 0], clear=True, cell_size=.1)\n",
     "dsk2 = Disk(center=[0, .25, 0])\n",
-    "# dsk3 = Disk(radius=.3)\n",
-    "rec1 = Rectangle(y=-.05, x=-.5, length=10, width=.1)\n",
+    "dsk3 = Disk(radius=.3)\n",
+    "rec1 = Rectangle(y=-.05, x=-.5, length=1, width=.1)\n",
     "\n",
-    "# cmplx_dom = dsk1 & dsk2 - (dsk3 + rec1)\n",
-    "cmplx_dom = rec1 \n",
+    "cmplx_dom = dsk1 & dsk2 - (dsk3 + rec1)\n",
     "\n",
     "plot(cmplx_dom)"
    ]
@@ -525,10 +499,28 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 13,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "9acdb0a75a5441a79109755fdb731a71",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "ViewInteractiveWidget(height=768, layout=Layout(height='auto', width='100%'), width=1024)"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
    "source": [
+    "from bvpy.domains import Sphere\n",
+    "from bvpy.utils.visu import plot, set_renderer\n",
+    "\n",
     "sph1 = Sphere(clear=True)\n",
     "sph2 = Sphere(center=[0,1,0], radius=.5)\n",
     "\n",
@@ -547,9 +539,24 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 14,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "132ea445eae74e7f836119b4b35da5b6",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "ViewInteractiveWidget(height=768, layout=Layout(height='auto', width='100%'), width=1024)"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
    "source": [
     "big = Sphere(radius=10, cell_type='tetra', clear=True)\n",
     "sml = Sphere(radius=9, cell_type='tetra')\n",
@@ -646,7 +653,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3.10.6 ('bvpyx-fenicsx-dev')",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },