diff --git a/setup.py b/setup.py
index b2189932..2baa270a 100644
--- a/setup.py
+++ b/setup.py
@@ -38,7 +38,7 @@
 
 # Dependencies
 install_requires = [
-    'sympy>=1.2',
+    'sympy>=1.2,<1.5',
     'h5py',
     'pytest',
     'pyyaml',
diff --git a/sympde/topology/basic.py b/sympde/topology/basic.py
index d6cf3f57..839bd63d 100644
--- a/sympde/topology/basic.py
+++ b/sympde/topology/basic.py
@@ -62,7 +62,7 @@ def __new__(cls, name, dim=None):
 
     @property
     def name(self):
-        return self._args[0]
+        return self.args[0]
 
     @property
     def target(self):
@@ -122,14 +122,14 @@ def __new__(cls, *args):
 
     @property
     def dim(self):
-        return self._args[0].dim
+        return self.args[0].dim
 
     def __len__(self):
-        return len(self._args)
+        return len(self.args)
 
     def complement(self, arg):
         if isinstance(arg, Union):
-            arg = arg._args
+            arg = arg.args
         elif isinstance(arg, BasicDomain):
             arg = [arg]
         elif arg is None:
@@ -137,16 +137,16 @@ def complement(self, arg):
         else:
             TypeError('Invalid argument {}'.format(arg))
 
-        return Union(*[i for i in self._args if (i not in arg)])
+        return Union(*[i for i in self.args if (i not in arg)])
 
     def __sub__(self, other):
         return self.complement(other)
 
     def todict(self):
-        return [i.todict() for i in self._args]
+        return [i.todict() for i in self.args]
 
     def as_tuple(self):
-        ls = [i for i in self._args]
+        ls = [i for i in self.args]
         return tuple(ls)
 
 #==============================================================================
@@ -161,11 +161,12 @@ def __new__(cls, *args, name=None):
         obj = Basic.__new__(cls, *args)
         obj._dim = sum(i.dim for i in args)
         obj._name = name
+
         return obj
 
     @property
     def domains(self):
-        return self._args
+        return self.args
 
 #==============================================================================
 class Interval(InteriorDomain):
@@ -193,7 +194,7 @@ def __new__(cls, name=None, coordinate=None, bounds=None):
 
     @property
     def name(self):
-        return self._args[0]
+        return self.args[0]
 
     @property
     def bounds(self):
@@ -209,37 +210,35 @@ class Boundary(BasicDomain):
     """
     def __new__(cls, name, domain, axis=None, ext=None):
 
-        if not( axis is None ):
-            assert(isinstance(axis, int))
+        if axis is not None:
+            assert isinstance(axis, int)
 
-        if not( ext is None ):
-            assert(isinstance(ext, int))
+        if ext is not None:
+            assert isinstance(ext, int)
+
+        obj = Basic.__new__(cls, name, domain, axis, ext)
 
-        obj = Basic.__new__(cls, name)
-        obj._domain = domain
-        obj._axis = axis
-        obj._ext = ext
         return obj
 
     @property
     def name(self):
-        return self._args[0]
+        return self.args[0]
 
     @property
     def domain(self):
-        return self._domain
-
-    @property
-    def dim(self):
-        return self.domain.dim
+        return self.args[1]
 
     @property
     def axis(self):
-        return self._axis
+        return self.args[2]
 
     @property
     def ext(self):
-        return self._ext
+        return self.args[3]
+
+    @property
+    def dim(self):
+        return self.domain.dim
 
     def _sympystr(self, printer):
         sstr = printer.doprint
@@ -264,19 +263,6 @@ def todict(self):
 
         return OrderedDict(sorted(d.items()))
 
-    def __hash__(self):
-        return hash(self.name)
-
-    def __lt__(self, other):
-        #add this method to avoid sympy error in Basic.compare
-        return 0
-
-    def __eq__(self, other):
-        if not isinstance(other, Boundary):
-            return False
-
-        return ( self.name == other.name ) and ( self.domain is other.domain )
-
 #==============================================================================
 class Interface(BasicDomain):
     """
@@ -303,15 +289,15 @@ def dim(self):
 
     @property
     def name(self):
-        return self._args[0]
+        return self.args[0]
 
     @property
     def minus(self):
-        return self._args[1]
+        return self.args[1]
 
     @property
     def plus(self):
-        return self._args[2]
+        return self.args[2]
 
     def _sympystr(self, printer):
         sstr = printer.doprint
diff --git a/sympde/topology/domain.py b/sympde/topology/domain.py
index 0fe1f1b1..0f0ddd71 100644
--- a/sympde/topology/domain.py
+++ b/sympde/topology/domain.py
@@ -9,11 +9,12 @@
 from collections import OrderedDict
 from collections import abc
 
+from sympy import Integer
 from sympy.core.singleton import Singleton
-from sympy.core.compatibility import with_metaclass
-from sympy.core import Basic, Symbol
+from sympy.core.compatibility import with_metaclass, is_sequence
+from sympy.core import Basic, symbols
 from sympy.core.containers import Tuple
-from sympy.tensor import IndexedBase
+from sympy.tensor import IndexedBase, Indexed
 from sympy.core import Add, Mul, Pow
 from sympy.core.expr import AtomicExpr
 
@@ -127,15 +128,15 @@ def __new__(cls, name, interiors=None, boundaries=None, dim=None,
 
     @property
     def name(self):
-        return self._args[0]
+        return self.args[0]
 
     @property
     def interior(self):
-        return self._args[1]
+        return self.args[1]
 
     @property
     def boundary(self):
-        return self._args[2]
+        return self.args[2]
 
     @property
     def connectivity(self):
@@ -166,7 +167,7 @@ def interior_names(self):
             return [self.interior.name]
 
         elif isinstance(self.interior, Union):
-            return [i.name for i in self.interior._args]
+            return [i.name for i in self.interior.args]
 
 
     def _sympystr(self, printer):
@@ -206,7 +207,7 @@ def get_boundary(self, name=None, axis=None, ext=None):
         # ...
 
         if isinstance(self.boundary, Union):
-            x = [i for i in self.boundary._args if i.name == name]
+            x = [i for i in self.boundary.args if i.name == name]
             if len(x) == 0:
                 raise ValueError('> could not find boundary {}'.format(name))
 
@@ -222,7 +223,7 @@ def get_boundary(self, name=None, axis=None, ext=None):
     def get_interior(self, name):
         """return interior by name."""
         if isinstance(self.interior, Union):
-            x = [i for i in self.interior._args if i.name == name]
+            x = [i for i in self.interior.args if i.name == name]
             if len(x) == 0:
                 raise ValueError('> could not find interior {}'.format(name))
 
@@ -293,8 +294,8 @@ def from_file( cls, filename ):
 
         if dtype == 'None': dtype = None
 
-        if not( dtype is None ):
-            constructor =  eval(dtype['type'])
+        if dtype is not None:
+            constructor = globals()[dtype['type']]
             return constructor(domain_name, **dtype['parameters'])
 
         # ... create sympde InteriorDomain (s)
@@ -421,40 +422,46 @@ def __hash__(self):
         return self.domain.__hash__() + self._periods.__hash__()
 
 #==============================================================================
-class Line(Domain):
-    def __new__(cls, name=None, bounds=None):
-        if name is None:
-            name = 'Line'
+# Ncube's properties (in addition to Domain's properties):
+#   . min_coords (default value is tuple of zeros)
+#   . max_coords (default value is tuple of ones)
+#
+class NCube(Domain):
 
-        x  = Symbol('x')
-        Ix = Interval(name, coordinate=x, bounds=bounds)
+    def __new__(cls, name, dim, min_coords, max_coords):
 
-        Gamma_1 = Boundary(r'\Gamma_1', Ix, axis=0, ext=-1)
-        Gamma_2 = Boundary(r'\Gamma_2', Ix, axis=0, ext=1)
+        assert isinstance(name, str)
+        assert isinstance(dim, (int, Integer))
+        assert isinstance(min_coords, (tuple, list, Tuple))
+        assert isinstance(max_coords, (tuple, list, Tuple))
 
-        interior   = Ix
-        boundaries = [Gamma_1, Gamma_2]
+        if not name:
+            raise ValueError("Name must be provided")
 
-        dtype = {'type': 'Line', 'parameters': {}}
+        if dim < 1:
+            raise ValueError("Number of dimensions must be at least 1")
 
-        return Domain(name, interiors=[interior],
-                      boundaries=boundaries, dtype=dtype)
+        if not (dim == len(min_coords) == len(max_coords)):
+            raise ValueError("Input arguments must have 'dim' components")
 
+        if not all(xmin < xmax for xmin, xmax in zip(min_coords, max_coords)):
+            raise ValueError("Min coordinates must be smaller than max")
 
-#==============================================================================
-# TODO: verify is the boundaries' domain should be a 1D Interval instead
-class Square(Domain):
-    def __new__(cls, name=None):
-        if name is None:
-            name = 'Square'
-
-        x  = Symbol('x')
-        Ix = Interval('Ix', coordinate=x)
+        # Choose coordinate names. TODO: use unique convention
+        if dim <= 3:
+            coord_names = ('x', 'y', 'z')[:dim]
+        else:
+            coord_names = 'x1:{}'.format(dim + 1)
 
-        y  = Symbol('y')
-        Iy = Interval('Iy', coordinate=y)
+        coordinates = symbols(coord_names)
+        intervals   = [Interval('I_{}'.format(c.name), coordinate=c, bounds=(xmin, xmax))
+                       for c, xmin, xmax in zip(coordinates, min_coords, max_coords)]
 
-        interior = ProductDomain(Ix, Iy, name=name)
+        if len(intervals) == 1:
+            interior = intervals[0]
+        else:
+            interior = ProductDomain(*intervals, name=name)
+            interior = InteriorDomain(interior)
 
         boundaries = []
         i = 1
@@ -463,51 +470,116 @@ def __new__(cls, name=None):
                 bnd_name = r'\Gamma_{}'.format(i)
                 Gamma = Boundary(bnd_name, interior, axis=axis, ext=ext)
                 boundaries += [Gamma]
-
                 i += 1
 
-        interior = InteriorDomain(interior)
+        # Choose which type to use:
+        #   a) if dim <= 3, use Line, Square or Cube;
+        #   b) if dim <= 4, use a generic 'NCube' type.
+        #
+        # Moreover, store all initialization parameters in a 'dtype' dictionary.
+        # This dictionary will be written to file when exporting the geometry,
+        # and it must contain all information necessary for building a new object
+        # by calling the appropriate constructor:
+        #
+        #   cls = globals()[dtype['type']]
+        #   domain = cls(name, **dtype['parameters'])
+        #
+        if dim == 1:
+            cls = Line
+            dtype = {'type': 'Line',
+                     'parameters': {'bounds': [min_coords[0], max_coords[0]]}}
+        elif dim == 2:
+            cls = Square
+            dtype = {'type': 'Square',
+                     'parameters': {'bounds1': [min_coords[0], max_coords[0]],
+                                    'bounds2': [min_coords[1], max_coords[1]]}}
+        elif dim == 3:
+            cls = Cube
+            dtype = {'type': 'Cube',
+                     'parameters': {'bounds1': [min_coords[0], max_coords[0]],
+                                    'bounds2': [min_coords[1], max_coords[1]],
+                                    'bounds3': [min_coords[2], max_coords[2]]}}
+        else:
+            dtype = {'type': 'NCube',
+                     'parameters': {'dim'       : dim,
+                                    'min_coords': [*min_coords],
+                                    'max_coords': [*max_coords]}}
 
-        dtype = {'type': 'Square', 'parameters': {}}
+        # Create instance of given type
+        obj = super().__new__(cls, name, interiors=[interior],
+                boundaries=boundaries, dtype=dtype)
 
-        return Domain(name, interiors=[interior],
-                      boundaries=boundaries, dtype=dtype)
+        # Store attributes in object
+        obj._coordinates = tuple(coordinates)
+        obj._min_coords  = tuple(min_coords)
+        obj._max_coords  = tuple(max_coords)
+
+        # Return object
+        return obj
+
+    @classmethod
+    def from_file(cls, filename):
+        msg = "Class method 'from_file' must be called on 'Domain' base class"
+        raise TypeError(msg)
+
+    @property
+    def min_coords(self):
+        return self._min_coords
+
+    @property
+    def max_coords(self):
+        return self._max_coords
 
 #==============================================================================
-# TODO: verify is the boundaries' domain should be a 2D ProductDomain instead
-class Cube(Domain):
-    def __new__(cls, name=None):
-        if name is None:
-            name = 'Cube'
+class Line(NCube):
 
-        x  = Symbol('x')
-        Ix = Interval('Ix', coordinate=x)
+    def __new__(cls, name='Line', bounds=(0, 1)):
+        dim = 1
+        min_coords = (bounds[0],)
+        max_coords = (bounds[1],)
+        return super().__new__(cls, name, dim, min_coords, max_coords)
 
-        y  = Symbol('y')
-        Iy = Interval('Iy', coordinate=y)
+    @property
+    def bounds(self):
+        return (self.min_coords[0], self.max_coords[0])
 
-        z  = Symbol('z')
-        Iz = Interval('Iz', coordinate=z)
+#==============================================================================
+class Square(NCube):
 
-        interior = ProductDomain(Ix, Iy, Iz, name=name)
+    def __new__(cls, name='Square', bounds1=(0, 1), bounds2=(0, 1)):
+        dim = 2
+        min_coords = (bounds1[0], bounds2[0])
+        max_coords = (bounds1[1], bounds2[1])
+        return super().__new__(cls, name, dim, min_coords, max_coords)
 
-        boundaries = []
-        i = 1
-        for axis in range(interior.dim):
-            for ext in [-1, 1]:
-                bnd_name = r'\Gamma_{}'.format(i)
-                Gamma = Boundary(bnd_name, interior, axis=axis, ext=ext)
-                boundaries += [Gamma]
+    @property
+    def bounds1(self):
+        return (self.min_coords[0], self.max_coords[0])
 
-                i += 1
+    @property
+    def bounds2(self):
+        return (self.min_coords[1], self.max_coords[1])
+
+#==============================================================================
+class Cube(NCube):
 
-        interior = InteriorDomain(interior)
+    def __new__(cls, name='Cube', bounds1=(0, 1), bounds2=(0, 1), bounds3=(0, 1)):
+        dim = 3
+        min_coords = (bounds1[0], bounds2[0], bounds3[0])
+        max_coords = (bounds1[1], bounds2[1], bounds3[1])
+        return super().__new__(cls, name, dim, min_coords, max_coords)
 
-        dtype = {'type': 'Cube', 'parameters': {}}
+    @property
+    def bounds1(self):
+        return (self.min_coords[0], self.max_coords[0])
 
-        return Domain(name, interiors=[interior],
-                      boundaries=boundaries, dtype=dtype)
+    @property
+    def bounds2(self):
+        return (self.min_coords[1], self.max_coords[1])
 
+    @property
+    def bounds3(self):
+        return (self.min_coords[2], self.max_coords[2])
 
 #==============================================================================
 class BoundaryVector(IndexedBase):
@@ -540,7 +612,7 @@ def __new__(cls, domain):
 
     @property
     def domain(self):
-        return self._args[0]
+        return self.args[0]
 
 class DomainArea(BasicArea):
     pass
@@ -575,20 +647,18 @@ def __new__(cls, *args, **options):
             return r
 
     @classmethod
-    def eval(cls, *_args):
+    def eval(cls, *args):
         """."""
 
-        if not _args:
+        if not args:
             return
 
-        if not len(_args) == 1:
+        if not len(args) == 1:
             raise ValueError('Expecting one argument')
 
-        expr = _args[0]
+        expr = args[0]
         if isinstance(expr, Union):
-            args = expr._args
-            args = [cls.eval(a) for a in expr.args]
-            return Add(*args)
+            return Add(*[cls.eval(a) for a in expr.args])
 
         elif isinstance(expr, ElementDomain):
             return ElementArea(expr)
diff --git a/sympde/topology/gallery.py b/sympde/topology/gallery.py
deleted file mode 100644
index ccb934f4..00000000
--- a/sympde/topology/gallery.py
+++ /dev/null
@@ -1,4 +0,0 @@
-# coding: utf-8
-
-
-
diff --git a/sympde/topology/tests/test_gallery.py b/sympde/topology/tests/test_gallery.py
index 47382203..b75ef824 100644
--- a/sympde/topology/tests/test_gallery.py
+++ b/sympde/topology/tests/test_gallery.py
@@ -3,11 +3,11 @@
 from sympy.abc import x,y,z
 from sympy import Tuple
 
-from sympde.topology import Interval, ProductDomain
-from sympde.topology import Domain, Line, Square, Cube
+from sympde.topology import Interval, ProductDomain, InteriorDomain, Domain
+from sympde.topology import Line, Square, Cube
 
 #==============================================================================
-def test_interval_1():
+def test_interval():
     Ix = Interval('Ix', coordinate=x)
     Iy = Interval('Iy', coordinate=y)
     Iz = Interval('Iz', coordinate=z)
@@ -19,35 +19,214 @@ def test_interval_1():
     assert(D_xy.dim == 2)
 
 #==============================================================================
-def test_gallery():
-    # ...
-    line = Line()
-    assert(line.dim == 1)
-    assert(len(line.boundary) == 2)
-    # ...
-
-    # ...
-    cube = Cube()
-    assert(cube.dim == 3)
-    assert(len(cube.boundary) == 6)
-    # ...
+def test_unit_line():
+
+    # Create 1D domain (Line) from interval [0, 1]
+    domain = Line('line')
+
+    assert isinstance(domain, Line)
+
+    # BasicDomain's attributes
+    assert domain.dim  == 1
+    assert domain.name == 'line'
+    assert domain.coordinates == x
+
+    # Domain's attributes
+    assert isinstance(domain.interior, Interval)
+    assert len(domain.boundary) == 2
+    assert domain.dtype == {'type': 'Line',
+                            'parameters': {'bounds': [0, 1]}}
+
+    # NCube's attributes
+    assert domain.min_coords == (0,)
+    assert domain.max_coords == (1,)
+
+    # Line's attributes
+    assert domain.bounds == (0, 1)
+
+    # Export to file, read it again and compare with original domain
+    domain.export('domain.h5')
+    D = Domain.from_file('domain.h5')
+    assert D == domain
+
+#==============================================================================
+def test_generic_line():
+
+    # Create 1D domain (Line) from interval [-3, 4]
+    domain = Line('line', bounds=(-3, 4))
+
+    assert isinstance(domain, Line)
+
+    # BasicDomain's attributes
+    assert domain.dim  == 1
+    assert domain.name == 'line'
+    assert domain.coordinates == x
+
+    # Domain's attributes
+    assert isinstance(domain.interior, Interval)
+    assert len(domain.boundary) == 2
+    assert domain.dtype == {'type': 'Line',
+                            'parameters': {'bounds': [-3, 4]}}
+
+    # NCube's attributes
+    assert domain.min_coords == (-3,)
+    assert domain.max_coords == ( 4,)
+
+    # Line's attributes
+    assert domain.bounds == (-3, 4)
+
+    # Export to file, read it again and compare with original domain
+    domain.export('domain.h5')
+    D = Domain.from_file('domain.h5')
+    assert D == domain
 
 #==============================================================================
-def test_square():
-    square = Square('square')
+def test_unit_square():
+
+    # Create 2D square domain [0, 1]^2
+    domain = Square('square')
+
+    assert isinstance(domain, Square)
+
+    # BasicDomain's attributes
+    assert domain.dim  == 2
+    assert domain.name == 'square'
+    assert domain.coordinates == (x, y)
 
-    assert(square.dim == 2)
-    assert(len(square.boundary) == 4)
+    # Domain's attributes
+    assert isinstance(domain.interior, InteriorDomain)
+    assert isinstance(domain.interior.target, ProductDomain)
+    assert all(isinstance(i, Interval) for i in domain.interior.target.domains)
+    assert len(domain.interior.target.domains) == 2
+    assert len(domain.boundary) == 4
+    assert domain.dtype == {'type': 'Square',
+                            'parameters': {'bounds1': [0, 1],
+                                           'bounds2': [0, 1]}}
 
-    square.export('square.h5')
+    # NCube's attributes
+    assert domain.min_coords == (0, 0)
+    assert domain.max_coords == (1, 1)
+
+    # Square's attributes
+    assert domain.bounds1 == (0, 1)
+    assert domain.bounds2 == (0, 1)
+
+    # Export to file, read it again and compare with original domain
+    domain.export('domain.h5')
+    D = Domain.from_file('domain.h5')
+    assert D == domain
+
+#==============================================================================
+def test_rectangle():
 
-    # read it again
-    D = Domain.from_file('square.h5')
+    # Create 2D rectangular domain [1, 5] X [3, 7]
+    domain = Square('rectangle', bounds1=(1, 5), bounds2=(3, 7))
 
-    # TODO BUG not working after naming boundaries \Gamma_1 etc
-#    assert( D == square )
+    assert isinstance(domain, Square)
 
+    # BasicDomain's attributes
+    assert domain.dim  == 2
+    assert domain.name == 'rectangle'
+    assert domain.coordinates == (x, y)
 
+    # Domain's attributes
+    assert isinstance(domain.interior, InteriorDomain)
+    assert isinstance(domain.interior.target, ProductDomain)
+    assert all(isinstance(i, Interval) for i in domain.interior.target.domains)
+    assert len(domain.interior.target.domains) == 2
+    assert len(domain.boundary) == 4
+    assert domain.dtype == {'type': 'Square',
+                            'parameters': {'bounds1': [1, 5],
+                                           'bounds2': [3, 7]}}
+
+    # NCube's attributes
+    assert domain.min_coords == (1, 3)
+    assert domain.max_coords == (5, 7)
+
+    # Square's attributes
+    assert domain.bounds1 == (1, 5)
+    assert domain.bounds2 == (3, 7)
+
+    # Export to file, read it again and compare with original domain
+    domain.export('domain.h5')
+    D = Domain.from_file('domain.h5')
+    assert D == domain
+
+#==============================================================================
+def test_unit_cube():
+
+    # Create 3D cube domain [0, 1]^3
+    domain = Cube('cube')
+
+    assert isinstance(domain, Cube)
+
+    # Check object attributes
+    assert domain.dim  == 3
+    assert domain.name == 'cube'
+    assert domain.coordinates == (x, y, z)
+
+    # Domain's attributes
+    assert isinstance(domain.interior, InteriorDomain)
+    assert isinstance(domain.interior.target, ProductDomain)
+    assert all(isinstance(i, Interval) for i in domain.interior.target.domains)
+    assert len(domain.interior.target.domains) == 3
+    assert len(domain.boundary) == 6
+    assert domain.dtype == {'type': 'Cube',
+                            'parameters': {'bounds1': [0, 1],
+                                           'bounds2': [0, 1],
+                                           'bounds3': [0, 1]}}
+
+    # NCube's attributes
+    assert domain.min_coords == (0, 0, 0)
+    assert domain.max_coords == (1, 1, 1)
+
+    # Cube's attributes
+    assert domain.bounds1 == (0, 1)
+    assert domain.bounds2 == (0, 1)
+    assert domain.bounds3 == (0, 1)
+
+    # Export to file, read it again and compare with original domain
+    domain.export('domain.h5')
+    D = Domain.from_file('domain.h5')
+    assert D == domain
+
+#==============================================================================
+def test_orthogonal_hexahedron():
+
+    # Create 3D orthogonal hexahedron [-1, 1] X [0, 10] X [0, 2]
+    domain = Cube('hexahedron', bounds1=(-1, 1), bounds2=(0, 10), bounds3=(0, 2))
+
+    assert isinstance(domain, Cube)
+
+    # Check object attributes
+    assert domain.dim  == 3
+    assert domain.name == 'hexahedron'
+    assert domain.coordinates == (x, y, z)
+
+    # Domain's attributes
+    assert isinstance(domain.interior, InteriorDomain)
+    assert isinstance(domain.interior.target, ProductDomain)
+    assert all(isinstance(i, Interval) for i in domain.interior.target.domains)
+    assert len(domain.interior.target.domains) == 3
+    assert len(domain.boundary) == 6
+    assert domain.dtype == {'type': 'Cube',
+                            'parameters': {'bounds1': [-1, 1],
+                                           'bounds2': [0, 10],
+                                           'bounds3': [0,  2]}}
+
+    # NCube's attributes
+    assert domain.min_coords == (-1, 0, 0)
+    assert domain.max_coords == (1, 10, 2)
+
+    # Cube's attributes
+    assert domain.bounds1 == (-1, 1)
+    assert domain.bounds2 == (0, 10)
+    assert domain.bounds3 == (0,  2)
+
+    # Export to file, read it again and compare with original domain
+    domain.export('domain.h5')
+    D = Domain.from_file('domain.h5')
+    assert D == domain
 
 #==============================================================================
 # CLEAN UP SYMPY NAMESPACE
@@ -57,9 +236,9 @@ def teardown_module():
     from sympy import cache
     cache.clear_cache()
 
-    # Remove output file generated by test_square()
+    # Remove output file generated by tests
     import os
-    fname = 'square.h5'
+    fname = 'domain.h5'
     if os.path.exists(fname):
         os.remove(fname)