diff --git a/pyproject.toml b/pyproject.toml
index 7189fb20..faf0aeb7 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name            = "sympde"
-version         = "0.19.2"
+version         = "0.20.0"
 description     = "Symbolic calculus for partial differential equations (and variational forms)"
 readme          = "README.rst"
 requires-python = ">= 3.9"
diff --git a/sympde/expr/evaluation.py b/sympde/expr/evaluation.py
index 9edda70a..931c57e3 100644
--- a/sympde/expr/evaluation.py
+++ b/sympde/expr/evaluation.py
@@ -1,13 +1,5 @@
-# coding: utf-8
-
-from itertools import product
-
-
-import numpy as np
-
-from sympy import Abs, S, cacheit
+from sympy import Abs, S, cacheit, Expr
 from sympy import Indexed, Matrix, ImmutableDenseMatrix
-from sympy import expand
 from sympy.core import Basic, Symbol
 from sympy.core import Add, Mul, Pow
 from sympy.core.expr import AtomicExpr
@@ -19,19 +11,20 @@
 from sympde.core.algebra import (Dot_1d,
                                  Dot_2d, Inner_2d, Cross_2d,
                                  Dot_3d, Inner_3d, Cross_3d)
-from sympde.core.utils import random_string
+#from sympde.core.utils import random_string
 
 from sympde.calculus import jump, avg, minus, plus
 from sympde.calculus import Jump, is_zero
 from sympde.calculus.core import _generic_ops, _diff_ops
+from sympde.calculus.core import MinusInterfaceOperator, PlusInterfaceOperator
 from sympde.calculus.matrices import SymbolicDeterminant, Inverse, Transpose
-from sympde.calculus.matrices import MatSymbolicPow, MatrixElement, SymbolicTrace
+from sympde.calculus.matrices import MatrixElement, SymbolicTrace
 
 from sympde.topology.basic   import BasicDomain, Union, Interval
 from sympde.topology.basic   import Boundary, Interface
 from sympde.topology.basic   import InteriorDomain
 from sympde.topology.domain  import NormalVector, TangentVector, NCube, NCubeInterior
-from sympde.topology.mapping import JacobianSymbol, InterfaceMapping, MultiPatchMapping, JacobianInverseSymbol
+from sympde.topology.mapping import JacobianSymbol, InterfaceMapping, JacobianInverseSymbol
 from sympde.topology.mapping import LogicalExpr, PullBack
 
 # TODO fix circular dependency between sympde.expr.evaluation and sympde.topology.mapping
@@ -40,7 +33,6 @@
 from sympde.topology.space import VectorFunction
 from sympde.topology.space import IndexedVectorFunction
 from sympde.topology.space import Trace
-from sympde.topology.space import element_of
 from sympde.topology.space import ScalarFunctionSpace
 
 from sympde.topology.derivatives import _partial_derivatives
@@ -132,10 +124,7 @@ def _get_trials_tests(expr, *, flatten=False):
     tests : tuple
         All scalar test functions in the given expression, with vector
         functions decomposed into their scalar components.
-
-
     """
-
     if not isinstance(expr, (BasicForm, BasicExpr)):
         raise TypeError("Expression must be of type BasicForm or BasicExpr, got '{}' instead".format(type(expr)))
 
@@ -187,7 +176,7 @@ def _to_matrix_form(expr, *, trials=None, tests=None, domain=None):
     tests : iterable
         List of all scalar test functions (after unpacking vector functions).
 
-    domain : Domain
+    domain : BasicDomain
         domain of expr
 
     Returns
@@ -434,10 +423,29 @@ def _nullify(expr, u, us):
 
 #==============================================================================
 class KernelExpression(Basic):
+    """
+    A generic SymPDE expression `expr` defined over the `target` domain. This
+    is a lightweight object which stores just the input values. In addition, if
+    the expression is a 1x1 matrix, the content of its only entry is taken.
+
+    Parameters
+    ----------
+    target : BasicDomain (from sympde.topology.basic)
+        The domain in which the expression is defined. This may be a Domain, a
+        Boundary, an Interface, or an object of any subclasses.
+
+    expr : Expr (from sympy)
+        The mathematical expression.
+    """
     def __new__(cls, target, expr):
+
+        assert isinstance(target, BasicDomain)
+        assert isinstance(expr, Expr)
+
         if isinstance(expr, (Matrix, ImmutableDenseMatrix)):
-            n,m = expr.shape
-            if n*m == 1: expr = expr[0,0]
+            n, m = expr.shape
+            if n * m == 1:
+                expr = expr[0,0]
 
         return Basic.__new__(cls, target, expr)
 
@@ -451,18 +459,31 @@ def expr(self):
 
 #==============================================================================
 class DomainExpression(KernelExpression):
-    pass
+    def __new__(cls, target, expr):
+        assert isinstance(target, InteriorDomain)
+        return KernelExpression.__new__(cls, target, expr)
 
 #==============================================================================
 class BoundaryExpression(KernelExpression):
-    pass
+    def __new__(cls, target, expr):
+        assert isinstance(target, Boundary)
+        return KernelExpression.__new__(cls, target, expr)
 
 #==============================================================================
+function_types = (ScalarFunction, VectorFunction, MinusInterfaceOperator, PlusInterfaceOperator)
+
 class InterfaceExpression(KernelExpression):
+
     def __new__(cls, target, u, v, expr):
+
+        assert isinstance(target, Interface)
+        assert isinstance(u, function_types)
+        assert isinstance(v, function_types)
+
         obj = KernelExpression.__new__(cls, target, expr)
         obj._trial = u
         obj._test  = v
+
         return obj
 
     @property
@@ -476,33 +497,45 @@ def trial(self):
 #==============================================================================
 class TerminalExpr(CalculusFunction):
     """
-     This class takes a SymPDE expression and transforms its vector operators into atomic ones
-     for a specified domain.
+    This class takes a SymPDE expression and transforms its vector operators
+    into atomic ones for a specified domain.
+
+    Parameters
+    ----------
+    expr : Expr | Matrix | ImmutableDenseMatrix (from sympy)
+                | LogicalExpr (from sympde.topology.mapping)
+        The mathematical expression.
 
-     Parameters
-     ----------
-     expr: sympy.Expr
-      The mathematical expression.
+    domain : BasicDomain (from sympde.topology.basic)
+        The domain in which the expression is defined. This may be a Domain, a
+        Boundary, an Interface, or an object of any subclasses.
 
-     domain: Domain
-      The domain in which the expression is defined.
+    evaluate : bool, optional
+        Whether the input expression should be evaluated (default = True).
 
-     Returns
-     -------
-      r: sympy.Expr or tuple of sympy.Expr
-       The atomized expression.
+    Returns
+    -------
+    Expr | tuple[Expr] | TerminalExpr
+        The atomized expression. A tuple of expressions is returned in the case
+        of `expr` being a BasicForm (i.e. a BilinearForm, LinearForm, or
+        Functional), as the domain could contain multiple patches. An
+        unevaluated `TerminalExpr` object (which only stores its arguments) is
+        returned if `evaluate is False`.
     """
-    def __new__(cls, expr, domain, **options):
+    def __new__(cls, expr: Expr, domain: BasicDomain, *, evaluate: bool=True):
 
-        if options.pop('evaluate', True):
-            r = cls.eval(expr, domain)
-        else:
-            r = None
+        expr_types = (Expr, Matrix, ImmutableDenseMatrix, LogicalExpr)
 
-        if r is None:
-            return Basic.__new__(cls, expr, domain)
+        assert isinstance(expr, expr_types)
+        assert isinstance(domain, BasicDomain)
+        assert isinstance(evaluate, bool)
+
+        if evaluate:
+            obj = cls.eval(expr, domain)
         else:
-            return r
+            obj = Basic.__new__(cls, expr, domain)
+
+        return obj
 
     def __getitem__(self, indices, **kw_args):
         if is_sequence(indices):
@@ -522,9 +555,33 @@ def domain(self):
     @classmethod
     @cacheit
     def eval(cls, expr, domain):
-        """."""
+        """
+        Process the whole AST of `expr` to generate a new expression where the
+        vector operators have been replaced by atomic operators. The input
+        `domain` is used to determine essential information like the number
+        of dimensions, the coordinate variables, the boundary orientation, etc.
+
+        Parameters
+        ----------
+        expr : Expr (from sympy)
+            The mathematical expression.
+
+        domain : BasicDomain (from sympde.topology.basic)
+            The domain in which the expression is defined. This may be a Domain, a
+            Boundary, an Interface, or an object of any subclasses.
+
+        Returns
+        -------
+        Expr | tuple[Expr] | TerminalExpr
+            The atomized expression. A tuple of expressions is returned in the case
+            of `expr` being a BasicForm (i.e. a BilinearForm, LinearForm, or
+            Functional), as the domain could contain multiple patches. An
+            unevaluated `TerminalExpr` object (which only stores its arguments) is
+            returned if `evaluate is False`.
+        """
 
         dim = domain.dim
+
         if isinstance(expr, Add):
             args = [cls.eval(a, domain=domain) for a in expr.args]
             o = args[0]
@@ -734,8 +791,8 @@ def eval(cls, expr, domain):
 
             ls = []
             for key in d_int:
-                domain, u,v = key
-                expr    = d_int[domain, u, v]
+                domain, u, v = key
+                expr = d_int[domain, u, v]
 
                 newexpr = cls.eval(expr, domain=domain)
                 if newexpr != 0:
@@ -787,7 +844,7 @@ def eval(cls, expr, domain):
             return new(*args)
 
         elif isinstance(expr, Trace):
-            # TODO treate different spaces
+            # TODO treat different spaces
             if expr.order == 0:
                 return cls.eval(expr.expr, domain=domain)
 
@@ -824,8 +881,8 @@ def eval(cls, expr, domain):
             return ImmutableDenseMatrix(lines)
 
         elif isinstance(expr, LogicalExpr):
-            domain    = expr.domain
-            expr      = cls(expr.expr, domain=domain)
+            domain = expr.domain
+            expr   = cls(expr.expr, domain=domain)
             return LogicalExpr(expr, domain=domain)
 
         return expr
diff --git a/sympde/expr/tests/test_expr.py b/sympde/expr/tests/test_expr.py
index c14fe9be..d29bd077 100644
--- a/sympde/expr/tests/test_expr.py
+++ b/sympde/expr/tests/test_expr.py
@@ -1,5 +1,3 @@
-# coding: utf-8
-
 # TODO: - add assert to every test
 
 import pytest
@@ -35,21 +33,14 @@
 from sympde.expr.evaluation import TerminalExpr
 
 #==============================================================================
-
 def test_linear_expr_2d_1():
 
     domain = Domain('Omega', dim=2)
-    x,y = domain.coordinates
-
-    kappa = Constant('kappa', is_real=True)
-    mu    = Constant('mu'   , is_real=True)
-
+    x, y = domain.coordinates
     V = ScalarFunctionSpace('V', domain)
+    u, u1, u2 = elements_of(V, names='u, u1, u2')
+    v, v1, v2 = elements_of(V, names='v, v1, v2')
 
-    u,u1,u2 = [element_of(V, name=i) for i in ['u', 'u1', 'u2']]
-    v,v1,v2 = [element_of(V, name=i) for i in ['v', 'v1', 'v2']]
-
-    # ...
     l = LinearExpr(v, x*y*v)
     print(l)
     print(l.expr)
@@ -57,9 +48,7 @@ def test_linear_expr_2d_1():
     # TODO
 #    print(l(v1+v2))
     print('')
-    # ...
 
-    # ...
     l = LinearExpr((v1,v2), x*v1 + y*v2)
     print(l)
     print(l.expr)
@@ -67,21 +56,15 @@ def test_linear_expr_2d_1():
     # TODO
 #    print(l(u1+v1, u2+v2))
     print('')
-    # ...
 
 #==============================================================================
 def test_linear_expr_2d_2():
 
     domain = Domain('Omega', dim=2)
-    x,y = domain.coordinates
-
-    kappa = Constant('kappa', is_real=True)
-    mu    = Constant('mu'   , is_real=True)
-
+    x, y = domain.coordinates
     V = VectorFunctionSpace('V', domain)
-
-    u,u1,u2 = [element_of(V, name=i) for i in ['u', 'u1', 'u2']]
-    v,v1,v2 = [element_of(V, name=i) for i in ['v', 'v1', 'v2']]
+    u, u1, u2 = elements_of(V, names='u, u1, u2')
+    v, v1, v2 = elements_of(V, names='v, v1, v2')
 
     g = Tuple(x,y)
     l = LinearExpr(v, dot(g, v))
@@ -91,9 +74,7 @@ def test_linear_expr_2d_2():
     # TODO
 #    print(l(v1+v2))
     print('')
-    # ...
 
-    # ...
     g1 = Tuple(x,0)
     g2 = Tuple(0,y)
     l = LinearExpr((v1,v2), dot(g1, v1) + dot(g2, v2))
@@ -103,7 +84,6 @@ def test_linear_expr_2d_2():
     # TODO
 #    print(l(u1+v1, u2+v2))
     print('')
-    # ...
 
 #==============================================================================
 def test_linear_form_2d_1():
@@ -111,136 +91,86 @@ def test_linear_form_2d_1():
     domain = Domain('Omega', dim=2)
     B1 = Boundary(r'\Gamma_1', domain)
 
-    x,y = domain.coordinates
-
-    kappa = Constant('kappa', is_real=True)
-    mu    = Constant('mu'   , is_real=True)
-    nn    = NormalVector('nn')
-
+    x, y = domain.coordinates
+    nn = NormalVector('nn')
     V = ScalarFunctionSpace('V', domain)
-    W = VectorFunctionSpace('W', domain)
 
-    u,u1,u2 = [element_of(V, name=i) for i in ['u', 'u1', 'u2']]
-    v,v1,v2 = [element_of(V, name=i) for i in ['v', 'v1', 'v2']]
-    # ...
+    u, u1, u2 = elements_of(V, names='u, u1, u2')
+    v, v1, v2 = elements_of(V, names='v, v1, v2')
+
     int_0 = lambda expr: integral(domain , expr)
     int_1 = lambda expr: integral(B1, expr)
 
     l = LinearForm(v, int_0(x*y*v))
+    assert l.domain == domain.interior
+    assert l(v1) == int_0(x*y*v1)
 
-    assert(l.domain == domain.interior)
-    assert(l(v1) == int_0(x*y*v1))
-    # ...
-
-    # ...
-    g  = Tuple(x**2, y**2)
-    l  = LinearForm(v, int_1(v*dot(g, nn)))
+    g = Tuple(x**2, y**2)
+    l = LinearForm(v, int_1(v*dot(g, nn)))
     print(l)
+    assert l.domain == B1
+#    assert l(v1) == int_1(v1*trace_1(g, B1)) # TODO [YG, 27.01.2021]: fix trace
+    assert l(v1) == int_1(v1*dot(g, nn))
 
-    assert(l.domain == B1)
-#    assert(l(v1) == int_1(v1*trace_1(g, B1))) # TODO [YG, 27.01.2021]: fix trace
-    assert(l(v1) == int_1(v1*dot(g, nn)))
-    # ...
-
-    # ...
     g = Tuple(x**2, y**2)
     l = LinearForm(v, int_1(v*dot(g, nn)) + int_0(x*y*v))
-
-    assert(len(l.domain.args) == 2)
+    assert len(l.domain.args) == 2
     for i in l.domain.args:
-        assert(i in [domain.interior, B1])
+        assert i in [domain.interior, B1]
+#    assert l(v1) == int_1(v1*trace_1(g, B1)) + int_0(x*y*v1) # TODO [YG, 27.01.2021]: fix trace
+    assert l(v1) == int_1(v1*dot(g, nn)) + int_0(x*y*v1)
 
-#    assert(l(v1) == int_1(v1*trace_1(g, B1)) + int_0(x*y*v1)) # TODO [YG, 27.01.2021]: fix trace
-    assert(l(v1) == int_1(v1*dot(g, nn)) + int_0(x*y*v1))
-    # ...
-
-    # ...
     l1 = LinearForm(v1, int_0(x*y*v1))
-    l = LinearForm(v, l1(v))
-
-    assert(l.domain == domain.interior)
-    assert(l(u1) == int_0(x*y*u1))
-    # ...
+    l  = LinearForm(v, l1(v))
+    assert l.domain == domain.interior
+    assert l(u1) == int_0(x*y*u1)
 
-    # ...
     g = Tuple(x,y)
     l1 = LinearForm(u1, int_0(x*y*u1))
     l2 = LinearForm(u2, int_0(dot(grad(u2), g)))
+    l  = LinearForm(v, l1(v) + l2(v))
+    assert l.domain == domain.interior
+    assert l(v1) == int_0(x*y*v1) + int_0(dot(grad(v1), g))
 
-    l = LinearForm(v, l1(v) + l2(v))
-
-    assert(l.domain == domain.interior)
-    assert(l(v1) == int_0(x*y*v1) + int_0(dot(grad(v1), g)))
-    # ...
-
-    # ...
-    pn, wn = [element_of(V, name=i) for i in ['pn', 'wn']]
-
-    tau   = element_of(V, name='tau')
-    sigma = element_of(V, name='sigma')
-
-    Re    = Constant('Re', real=True)
-    dt    = Constant('dt', real=True)
-    alpha = Constant('alpha', real=True)
-
+    pn, wn, tau, sigma = elements_of(V, names='pn, wn, tau, sigma')
+    Re = Constant('Re', real=True)
+    dt = Constant('dt', real=True)
     l1 = LinearForm(tau, int_0(bracket(pn, wn)*tau - 1./Re * dot(grad(tau), grad(wn))))
-
-
-
-    l = LinearForm((tau, sigma), dt*l1(tau))
-
-    assert(l.domain == domain.interior)
-    assert(l(u1,u2).expand() == int_0(-1.0*dt*dot(grad(u1), grad(wn))/Re) +
-                       int_0(dt*u1*bracket(pn, wn)))
-    # ...
+    l  = LinearForm((tau, sigma), dt*l1(tau))
+    assert l.domain == domain.interior
+    assert l(u1,u2).expand() == int_0(-1.0*dt*dot(grad(u1), grad(wn))/Re) + \
+                       int_0(dt*u1*bracket(pn, wn))
 
 #==============================================================================
 def test_linear_form_2d_2():
 
     domain = Domain('Omega', dim=2)
-    B1 = Boundary(r'\Gamma_1', domain)
-
-    x,y = domain.coordinates
-
-    kappa = Constant('kappa', is_real=True)
-    mu    = Constant('mu'   , is_real=True)
+    x, y = domain.coordinates
 
     V = VectorFunctionSpace('V', domain)
+    u, u1, u2 = elements_of(V, names='u, u1, u2')
+    v, v1, v2 = elements_of(V, names='v, v1, v2')
 
-    u,u1,u2 = [element_of(V, name=i) for i in ['u', 'u1', 'u2']]
-    v,v1,v2 = [element_of(V, name=i) for i in ['v', 'v1', 'v2']]
-
-    # ...
-    int_0 = lambda expr: integral(domain , expr)
-    int_1 = lambda expr: integral(B1, expr)
+    int_0 = lambda expr: integral(domain, expr)
 
-    g = Matrix((x,y))
+    g = Matrix((x, y))
     l = LinearForm(v, int_0(dot(g, v)))
+    assert l.domain == domain.interior
+    assert l(v1) == int_0(dot(g, v1))
 
-    assert(l.domain == domain.interior)
-    assert(l(v1) == int_0(dot(g, v1)))
-    # ...
-
-    # ...
-    g = Matrix((x,y))
+    g = Matrix((x, y))
     l1 = LinearForm(v1, int_0(dot(g, v1)))
     l = LinearForm(v, l1(v))
+    assert l.domain == domain.interior
+    assert l(u1) == int_0(dot(g, u1))
 
-    assert(l.domain == domain.interior)
-    assert(l(u1) == int_0(dot(g, u1)))
-    # ...
-
-    # ...
-    g1 = Matrix((x,0))
-    g2 = Matrix((0,y))
+    g1 = Matrix((x, 0))
+    g2 = Matrix((0, y))
     l1 = LinearForm(v1, int_0(dot(v1, g1)))
     l2 = LinearForm(v2, int_0(dot(v2, g2)))
-
-    l = LinearForm(v, l1(v) + l2(v))
-
-    assert(l.domain == domain.interior)
-    assert(l(u) == int_0(dot(u, g1)) + int_0(dot(u, g2)))
-    # ...
+    l  = LinearForm(v, l1(v) + l2(v))
+    assert l.domain == domain.interior
+    assert l(u) == int_0(dot(u, g1)) + int_0(dot(u, g2))
 
 #==============================================================================
 def test_bilinear_form_2d_1():
@@ -248,120 +178,78 @@ def test_bilinear_form_2d_1():
     domain = Domain('Omega', dim=2)
     B1 = Boundary(r'\Gamma_1', domain)
 
-    x,y = domain.coordinates
-
     kappa = Constant('kappa', is_real=True)
-    mu    = Constant('mu'   , is_real=True)
     nn    = NormalVector('nn')
 
     V = ScalarFunctionSpace('V', domain)
+    u, u1, u2 = elements_of(V, names='u, u1, u2')
+    v, v1, v2 = elements_of(V, names='v, v1, v2')
 
-    u,u1,u2 = [element_of(V, name=i) for i in ['u', 'u1', 'u2']]
-    v,v1,v2 = [element_of(V, name=i) for i in ['v', 'v1', 'v2']]
-
-    # ...
     int_0 = lambda expr: integral(domain , expr)
     int_1 = lambda expr: integral(B1, expr)
 
     a = BilinearForm((u,v), int_0(u*v))
+    assert a.domain == domain.interior
+    assert a(u1,v1) == int_0(u1*v1)
 
-    assert(a.domain == domain.interior)
-    assert(a(u1,v1) == int_0(u1*v1))
-    # ...
-
-    # ...
     a = BilinearForm((u,v), int_0(u*v + dot(grad(u), grad(v))))
+    assert a.domain == domain.interior
+    assert a(u1,v1) == int_0(u1*v1) + int_0(dot(grad(u1), grad(v1)))
 
-    assert(a.domain == domain.interior)
-    assert(a(u1,v1) == int_0(u1*v1) + int_0(dot(grad(u1), grad(v1))))
-    # ...
-
-    # ...
     a = BilinearForm((u,v), int_1(v*dot(grad(u), nn)))
-
-    assert(a.domain == B1)
-#    assert(a(u1,v1) == int_1(v1*trace_1(grad(u1), B1))) # TODO [YG, 27.01.2021]: fix trace
+    assert a.domain == B1
+#    assert a(u1,v1) == int_1(v1*trace_1(grad(u1), B1)) # TODO [YG, 27.01.2021]: fix trace
     assert a(u1,v1) == int_1(v1*dot(grad(u1), nn))
-    # ...
 
-    # ...
     a = BilinearForm((u,v), int_0(u*v) + int_1(v*dot(grad(u), nn)))
-
     # TODO a.domain are not ordered
-    assert(len(a.domain.args) == 2)
+    assert len(a.domain.args) == 2
     for i in a.domain.args:
-        assert(i in [domain.interior, B1])
-
-#    assert(a(u1,v1) == int_0(u1*v1) + int_1(v1*trace_1(grad(u1), B1))) # TODO [YG, 27.01.2021]: fix trace
+        assert i in [domain.interior, B1]
+#    assert a(u1,v1) == int_0(u1*v1) + int_1(v1*trace_1(grad(u1), B1)) # TODO [YG, 27.01.2021]: fix trace
     assert a(u1,v1) == int_0(u1*v1) + int_1(v1*dot(grad(u1), nn))
-    # ...
 
-    # ...
     a1 = BilinearForm((u1,v1), int_0(u1*v1))
     a = BilinearForm((u,v), a1(u,v))
-
     assert(a.domain == domain.interior)
     assert(a(u2,v2) == int_0(u2*v2))
-    # ...
 
-    # ...
     a1 = BilinearForm((u1,v1), int_0(u1*v1))
     a2 = BilinearForm((u2,v2), int_0(dot(grad(u2), grad(v2))))
     a = BilinearForm((u,v), a1(u,v) + kappa*a2(u,v))
-
     assert(a.domain == domain.interior)
     assert(a(u,v) == int_0(u*v) + int_0(kappa*dot(grad(u), grad(v))))
-    # ...
 
 #==============================================================================
 def test_bilinear_form_2d_2():
 
     domain = Domain('Omega', dim=2)
-    B1 = Boundary(r'\Gamma_1', domain)
-
-    x,y = domain.coordinates
-
     kappa = Constant('kappa', is_real=True)
-    mu    = Constant('mu'   , is_real=True)
-
     V = VectorFunctionSpace('V', domain)
+    u, u1, u2 = elements_of(V, names='u, u1, u2')
+    v, v1, v2 = elements_of(V, names='v, v1, v2')
 
-    u,u1,u2 = [element_of(V, name=i) for i in ['u', 'u1', 'u2']]
-    v,v1,v2 = [element_of(V, name=i) for i in ['v', 'v1', 'v2']]
-
-    # ...
     int_0 = lambda expr: integral(domain , expr)
-    int_1 = lambda expr: integral(B1, expr)
 
     a = BilinearForm((u,v), int_0(dot(u,v)))
+    assert a.domain == domain.interior
+    assert a(u1,v1) == int_0(dot(u1,v1))
 
-    assert(a.domain == domain.interior)
-    assert(a(u1,v1) == int_0(dot(u1,v1)))
-    # ...
-
-    # ...
     a = BilinearForm((u,v), int_0(dot(u,v) + inner(grad(u), grad(v))))
+    assert a.domain == domain.interior
+    assert a(u1,v1) == int_0(dot(u1,v1)) + int_0(inner(grad(u1), grad(v1)))
 
-    assert(a.domain == domain.interior)
-    assert(a(u1,v1) == int_0(dot(u1,v1)) + int_0(inner(grad(u1), grad(v1))))
-    # ...
-
-    # ...
     a1 = BilinearForm((u1,v1), int_0(dot(u1,v1)))
     a = BilinearForm((u,v), a1(u,v))
+    assert a.domain == domain.interior
+    assert a(u2,v2) == int_0(dot(u2,v2))
 
-    assert(a.domain == domain.interior)
-    assert(a(u2,v2) == int_0(dot(u2,v2)))
-    # ...
-
-    # ...
     a1 = BilinearForm((u1,v1), int_0(dot(u1,v1)))
     a2 = BilinearForm((u2,v2), int_0(inner(grad(u2), grad(v2))))
     a = BilinearForm((u,v), a1(u,v) + kappa*a2(u,v))
 
-    assert(a.domain == domain.interior)
-    assert(a(u,v) == int_0(dot(u,v)) + int_0(kappa*inner(grad(u), grad(v))))
-    # ...
+    assert a.domain == domain.interior
+    assert a(u,v) == int_0(dot(u,v)) + int_0(kappa*inner(grad(u), grad(v)))
 
 #==============================================================================
 def test_terminal_expr_linear_2d_1():
@@ -369,55 +257,41 @@ def test_terminal_expr_linear_2d_1():
     domain = Domain('Omega', dim=2)
     B1 = Boundary(r'\Gamma_1', domain)
 
-    x,y = domain.coordinates
+    x, y = domain.coordinates
 
     kappa = Constant('kappa', is_real=True)
     mu    = Constant('mu'   , is_real=True)
     nn    = NormalVector('nn')
 
     V = ScalarFunctionSpace('V', domain)
+    v, v1, v2 = elements_of(V, names='v, v1, v2')
 
-    u,u1,u2 = [element_of(V, name=i) for i in ['u', 'u1', 'u2']]
-    v,v1,v2 = [element_of(V, name=i) for i in ['v', 'v1', 'v2']]
-
-    # ...
-    int_0 = lambda expr: integral(domain , expr)
+    int_0 = lambda expr: integral(domain, expr)
     int_1 = lambda expr: integral(B1, expr)
 
     l = LinearForm(v, int_0(x*y*v))
-
     print(TerminalExpr(l, domain))
     print('')
-    # ...
 
-    # ...
     l = LinearForm(v, int_0(x*y*v + v))
     print(TerminalExpr(l, domain))
     print('')
-    # ...
 
-    # ...
     g = Matrix((x**2, y**2))
     l = LinearForm(v, int_1(v*dot(g, nn)))
     print(TerminalExpr(l, domain))
     print('')
-    # ...
 
-    # ...
     g = Matrix((x**2, y**2))
     l = LinearForm(v, int_1(v*dot(g, nn)) + int_0(x*y*v))
     print(TerminalExpr(l, domain))
     print('')
-    # ...
 
-    # ...
     l1 = LinearForm(v1, int_0(x*y*v1))
     l = LinearForm(v, l1(v))
     print(TerminalExpr(l, domain))
     print('')
-    # ...
 
-    # ...
     g = Matrix((x,y))
     l1 = LinearForm(v1, int_0(x*y*v1))
     l2 = LinearForm(v2, int_0(dot(grad(v2), g)))
@@ -425,17 +299,13 @@ def test_terminal_expr_linear_2d_1():
     l = LinearForm(v, l1(v) + l2(v))
     print(TerminalExpr(l, domain))
     print('')
-    # ...
 
-    # ...
     l1 = LinearForm(v1, int_0(x*y*v1))
     l2 = LinearForm(v1, int_0(v1))
     l = LinearForm(v, l1(v) + kappa*l2(v))
     print(TerminalExpr(l, domain))
     print('')
-    # ...
 
-    # ...
     g = Matrix((x**2, y**2))
     l1 = LinearForm(v1, int_0(x*y*v1))
     l2 = LinearForm(v1, int_0(v1))
@@ -443,7 +313,6 @@ def test_terminal_expr_linear_2d_1():
     l = LinearForm(v, l1(v) + kappa*l2(v) + mu*l3(v))
     print(TerminalExpr(l, domain))
     print('')
-    # ...
 
 #==============================================================================
 def test_terminal_expr_linear_2d_2():
@@ -451,74 +320,52 @@ def test_terminal_expr_linear_2d_2():
     domain = Domain('Omega', dim=2)
     B1 = Boundary(r'\Gamma_1', domain)
 
-    x,y = domain.coordinates
-
-    kappa = Constant('kappa', is_real=True)
-    mu    = Constant('mu'   , is_real=True)
-
+    x, y = domain.coordinates
     V = VectorFunctionSpace('V', domain)
+    v, v1, v2 = elements_of(V, names='v, v1, v2')
 
-    u,u1,u2 = [element_of(V, name=i) for i in ['u', 'u1', 'u2']]
-    v,v1,v2 = [element_of(V, name=i) for i in ['v', 'v1', 'v2']]
-
-    # ...
-    int_0 = lambda expr: integral(domain , expr)
+    int_0 = lambda expr: integral(domain, expr)
     int_1 = lambda expr: integral(B1, expr)
 
     g = Matrix((x,y))
     l = LinearForm(v, int_0(dot(g, v)))
     print(TerminalExpr(l, domain))
     print('')
-    # ...
 
-    # ...
     g = Matrix((x,y))
     l = LinearForm(v, int_0(dot(g, v) + div(v)))
     print(TerminalExpr(l, domain))
     print('')
-    # ...
 
     # TODO
-#    # ...
 #    g = Tuple(x**2, y**2)
 #    l = LinearForm(v, v*trace_1(g, B1))
 #    print(TerminalExpr(l))
 #    print('')
-#    # ...
 #
-#    # ...
 #    g = Tuple(x**2, y**2)
 #    l = LinearForm(v, v*trace_1(g, B1) + x*y*v)
 #    print(TerminalExpr(l))
 #    print('')
-#    # ...
 #
-#    # ...
 #    l1 = LinearForm(v1, x*y*v1)
 #    l = LinearForm(v, l1(v))
 #    print(TerminalExpr(l))
 #    print('')
-#    # ...
 #
-#    # ...
 #    g = Tuple(x,y)
 #    l1 = LinearForm(v1, x*y*v1)
 #    l2 = LinearForm(v2, dot(grad(v2), g))
-#
 #    l = LinearForm(v, l1(v) + l2(v))
 #    print(TerminalExpr(l))
 #    print('')
-#    # ...
 #
-#    # ...
 #    l1 = LinearForm(v1, x*y*v1)
 #    l2 = LinearForm(v1, v1)
 #    l = LinearForm(v, l1(v) + kappa*l2(v))
 #    print(TerminalExpr(l))
 #    print('')
-#    # ...
 #
-#    # ...
 #    g = Tuple(x**2, y**2)
 #    l1 = LinearForm(v1, x*y*v1)
 #    l2 = LinearForm(v1, v1)
@@ -526,33 +373,22 @@ def test_terminal_expr_linear_2d_2():
 #    l = LinearForm(v, l1(v) + kappa*l2(v) + mu*l3(v))
 #    print(TerminalExpr(l))
 #    print('')
-#    # ...
 
 #==============================================================================
 def test_terminal_expr_linear_2d_3():
 
     domain = Square()
-    B      = domain.boundary
-
-    x,y = domain.coordinates
-
-    kappa = Constant('kappa', is_real=True)
-    mu    = Constant('mu'   , is_real=True)
-    nn    = NormalVector('nn')
-
+    B  = domain.boundary
+    nn = NormalVector('nn')
     V = ScalarFunctionSpace('V', domain)
+    v = element_of(V, name='v')
 
-    u,u1,u2 = [element_of(V, name=i) for i in ['u', 'u1', 'u2']]
-    v,v1,v2 = [element_of(V, name=i) for i in ['v', 'v1', 'v2']]
-
-    # ...
-    int_0 = lambda expr: integral(domain , expr)
+    int_0 = lambda expr: integral(domain, expr)
     int_1 = lambda expr: integral(B, expr)
 
     l = LinearForm(v, int_1(dot(grad(v), nn)))
     print(TerminalExpr(l, domain))
     print('')
-    # ...
 
 #==============================================================================
 def test_terminal_expr_linear_2d_4():
@@ -560,24 +396,16 @@ def test_terminal_expr_linear_2d_4():
     D1 = InteriorDomain('D1', dim=2)
     D2 = InteriorDomain('D2', dim=2)
     domain = Union(D1, D2)
-
-    x,y = domain.coordinates
-
-    kappa = Constant('kappa', is_real=True)
-    mu    = Constant('mu'   , is_real=True)
-
+    x, y = domain.coordinates
     V = ScalarFunctionSpace('V', domain)
+    v = element_of(V, name='v')
 
-    u,u1,u2 = [element_of(V, name=i) for i in ['u', 'u1', 'u2']]
-    v,v1,v2 = [element_of(V, name=i) for i in ['v', 'v1', 'v2']]
-
-    # ...
     int_0 = lambda expr: integral(domain , expr)
 
     l = LinearForm(v, int_0(x*y*v))
     print(TerminalExpr(l, domain))
     print('')
-    # ...
+
 #==============================================================================
 def test_terminal_expr_linear_2d_5():
 
@@ -623,7 +451,6 @@ def test_terminal_expr_linear_2d_5():
 
     print(TerminalExpr(l, domain))
     print('')
-    # ...
 
 #==============================================================================
 def test_terminal_expr_bilinear_2d_1():
@@ -631,8 +458,6 @@ def test_terminal_expr_bilinear_2d_1():
     domain = Domain('Omega', dim=2)
     B1     = Boundary(r'\Gamma_1', domain)
 
-    x,y = domain.coordinates
-
     kappa = Constant('kappa', is_real=True)
     mu    = Constant('mu'   , is_real=True)
     eps   = Constant('eps', real=True)
@@ -640,75 +465,57 @@ def test_terminal_expr_bilinear_2d_1():
 
     V = ScalarFunctionSpace('V', domain)
 
-    u,u1,u2 = [element_of(V, name=i) for i in ['u', 'u1', 'u2']]
-    v,v1,v2 = [element_of(V, name=i) for i in ['v', 'v1', 'v2']]
-
-    # ...
+    u, u1, u2 = elements_of(V, names='u, u1, u2')
+    v, v1, v2 = elements_of(V, names='v, v1, v2')
 
-    int_0 = lambda expr: integral(domain , expr)
+    int_0 = lambda expr: integral(domain, expr)
     int_1 = lambda expr: integral(B1, expr)
 
     a = BilinearForm((u,v), int_0(u*v))
     print(a)
     print(TerminalExpr(a, domain))
     print('')
-    # ...
 
-    # ...
     a = BilinearForm((u,v), int_0(dot(grad(u),grad(v))))
     print(a)
     print(TerminalExpr(a, domain))
     print('')
-    # ...
 
-    # ...
     a = BilinearForm((u,v), int_0(u*v + dot(grad(u),grad(v))))
     print(a)
     print(TerminalExpr(a, domain))
     print('')
-    # ...
 
-    # ...
     a = BilinearForm((u,v), int_0(u*v + dot(grad(u),grad(v))) + int_1(v*dot(grad(u), nn)) )
     print(a)
     print(TerminalExpr(a, domain))
     print('')
-    # ...
 
-    # ...
     a = BilinearForm(((u1,u2),(v1,v2)), int_0(u1*v1 + u2*v2))
     print(a)
     print(TerminalExpr(a, domain))
     print('')
-    # ...
 
-    # ...
     a1 = BilinearForm((u1,v1), int_0(u1*v1))
     a = BilinearForm((u,v), a1(u,v))
     print(a)
     print(TerminalExpr(a, domain))
     print('')
-    # ...
 
-    # ...
     a1 = BilinearForm((u1,v1), int_0(u1*v1))
     a2 = BilinearForm((u2,v2), int_0(dot(grad(u2), grad(v2))))
     a = BilinearForm((u,v), a1(u,v) + a2(u,v))
     print(a)
     print(TerminalExpr(a, domain))
     print('')
-    # ...
 
-    # ...
     a1 = BilinearForm((u1,v1), int_0(u1*v1))
     a2 = BilinearForm((u2,v2), int_0(dot(grad(u2), grad(v2))))
     a = BilinearForm((u,v), a1(u,v) + kappa*a2(u,v))
     print(a)
     print(TerminalExpr(a, domain))
     print('')
-    # ...
 
-    # ...
     a1 = BilinearForm((u1,v1), int_0(u1*v1))
     a2 = BilinearForm((u2,v2), int_0(dot(grad(u2), grad(v2))))
     a3 = BilinearForm((u,v), int_1(v*dot(grad(u), nn)))
@@ -716,7 +523,6 @@ def test_terminal_expr_bilinear_2d_1():
     print(a)
     print(TerminalExpr(a, domain))
     print('')
-    # ...
 
     # ... Poisson with Nitsch method
     a0 = BilinearForm((u,v), int_0(dot(grad(u),grad(v))))
@@ -727,55 +533,35 @@ def test_terminal_expr_bilinear_2d_1():
     print(a)
     print(TerminalExpr(a, domain))
     print('')
-    # ...
 
 #==============================================================================
 def test_terminal_expr_bilinear_2d_2():
 
     domain = Domain('Omega', dim=2)
-    B1     = Boundary(r'\Gamma_1', domain)
-
-    x,y = domain.coordinates
-
-    kappa = Constant('kappa', is_real=True)
-    mu    = Constant('mu'   , is_real=True)
-    nn    = NormalVector('nn')
-
     V = VectorFunctionSpace('V', domain)
+    u, v = elements_of(V, names='u, v')
 
-    u,u1,u2 = [element_of(V, name=i) for i in ['u', 'u1', 'u2']]
-    v,v1,v2 = [element_of(V, name=i) for i in ['v', 'v1', 'v2']]
-
-    # ...
-    int_0 = lambda expr: integral(domain , expr)
-    int_1 = lambda expr: integral(B1, expr)
+    int_0 = lambda expr: integral(domain, expr)
 
     a = BilinearForm((u,v), int_0(dot(u,v)))
     print(TerminalExpr(a, domain))
     print('')
 
-    # ...
     a = BilinearForm((u,v), int_0(inner(grad(u),grad(v))))
     print(TerminalExpr(a, domain))
     print('')
-    # ...
 
-    # ...
     a = BilinearForm((u,v), int_0(dot(u,v) + inner(grad(u),grad(v))))
     print(TerminalExpr(a, domain))
     print('')
-    # ...
 
+#==============================================================================
 def test_terminal_expr_bilinear_2d_3():
 
     domain = Square()
-
     V = ScalarFunctionSpace('V', domain)
-
     B = domain.boundary
-
-    v = element_of(V, name='v')
-    u = element_of(V, name='u')
+    u, v = elements_of(V, names='u, v')
 
     kappa = Constant('kappa', is_real=True)
     mu    = Constant('mu'   , is_real=True)
@@ -786,66 +572,50 @@ def test_terminal_expr_bilinear_2d_3():
 
     # nitsche
     a0  = BilinearForm((u,v), int_0(dot(grad(v),grad(u))))
-
     a_B = BilinearForm((u,v), int_1(-u*dot(grad(v), nn) \
                               -v*dot(grad(u), nn) \
                               +kappa*u*v))
-
     a = BilinearForm((u,v), a0(u,v) + a_B(u,v))
-
     print(TerminalExpr(a, domain))
     print('')
 
     a = BilinearForm((u,v), int_0(u*v + dot(grad(u),grad(v))) + int_1(v*dot(grad(u), nn)) )
     print(TerminalExpr(a, domain))
     print('')
-    # ...
 
-    # ...
     a = BilinearForm((u,v), int_0(u*v))
     print(TerminalExpr(a, domain))
     print('')
-    # ...
 
-    # ...
     a1 = BilinearForm((u,v), int_0(u*v))
     a = BilinearForm((u,v), a1(u,v))
     print(TerminalExpr(a, domain))
     print('')
-    # ...
 
-    # ...
     a1 = BilinearForm((u,v), int_0(u*v))
     a2 = BilinearForm((u,v), int_0(dot(grad(u), grad(v))))
     a = BilinearForm((u,v), a1(u,v) + a2(u,v))
     print(TerminalExpr(a, domain))
     print('')
-    # ...
 
-    # ...
     a1 = BilinearForm((u,v), int_0(u*v))
     a2 = BilinearForm((u,v), int_0(dot(grad(u), grad(v))))
     a = BilinearForm((u,v), a1(u,v) + kappa*a2(u,v))
     print(TerminalExpr(a, domain))
     print('')
-    # ...
 
-    # ...
     a1 = BilinearForm((u,v), int_0(u*v))
     a2 = BilinearForm((u,v), int_0(dot(grad(u), grad(v))))
     a3 = BilinearForm((u,v), int_1(v*dot(grad(u), nn)))
     a = BilinearForm((u,v), a1(u,v) + kappa*a2(u,v) + mu*a3(u,v))
     print(TerminalExpr(a, domain))
     print('')
-#    # ...
 
 #==============================================================================
 def test_terminal_expr_bilinear_2d_4():
 
     domain = Domain('Omega', dim=2)
 
-    x,y = domain.coordinates
-
     V = VectorFunctionSpace('V', domain)
     W = ScalarFunctionSpace('W', domain)
 
@@ -854,7 +624,7 @@ def test_terminal_expr_bilinear_2d_4():
     p = element_of(W, name='p')
     q = element_of(W, name='q')
 
-    int_0 = lambda expr: integral(domain , expr)
+    int_0 = lambda expr: integral(domain, expr)
 
     # stokes
     a = BilinearForm((u,v), int_0(inner(grad(v), grad(u))))
@@ -878,14 +648,13 @@ def test_terminal_expr_bilinear_3d_1():
     u,v   = elements_of(V, names='u,v')
     um,vm = elements_of(VM, names='u,v')
 
-    int_0 = lambda expr: integral(domain , expr)
-    int_1 = lambda expr: integral(mapped_domain , expr)
+    int_0 = lambda expr: integral(domain, expr)
+    int_1 = lambda expr: integral(mapped_domain, expr)
 
     J   = M.det_jacobian
     det = dx1(M[0])*dx2(M[1])*dx3(M[2]) - dx1(M[0])*dx2(M[2])*dx3(M[1]) - dx1(M[1])*dx2(M[0])*dx3(M[2])\
         + dx1(M[1])*dx2(M[2])*dx3(M[0]) + dx1(M[2])*dx2(M[0])*dx3(M[1]) - dx1(M[2])*dx2(M[1])*dx3(M[0])
 
-
     a1 = BilinearForm((u,v), int_0(dot(grad(u),grad(v))))
     a2 = BilinearForm((um,vm), int_1(dot(grad(um),grad(vm))))
     a3 = BilinearForm((u,v), int_0(J*dot(grad(u),grad(v))))
@@ -926,7 +695,6 @@ def test_terminal_expr_3d_1():
 @pytest.mark.skip(reason="New linearize() function does not accept 'LinearExpr' objects")
 def test_linearize_expr_2d_1():
     domain = Domain('Omega', dim=2)
-    x,y = domain.coordinates
 
     V1 = ScalarFunctionSpace('V1', domain)
     W1 = VectorFunctionSpace('W1', domain)
@@ -934,12 +702,9 @@ def test_linearize_expr_2d_1():
     v1 = element_of(V1, name='v1')
     w1 = element_of(W1, name='w1')
 
-    alpha = Constant('alpha')
-
     F = element_of(V1, name='F')
     G = element_of(W1, 'G')
 
-
     # ...
     l = LinearExpr(v1, F**2*v1)
     a = linearize(l, F, trials='u1')
@@ -986,16 +751,10 @@ def test_linearize_expr_2d_1():
 @pytest.mark.skip(reason="New linearize() function does not accept 'LinearExpr' objects")
 def test_linearize_expr_2d_2():
     domain = Domain('Omega', dim=2)
-    x,y = domain.coordinates
 
     V1 = ScalarFunctionSpace('V1', domain)
-
     v1 = element_of(V1, name='v1')
-
-    alpha = Constant('alpha')
-
-    F = element_of(V1, name='F')
-    G = element_of(V1, name='G')
+    F  = element_of(V1, name='F')
 
     # ...
     l1 = LinearExpr(v1, F**2*v1)
@@ -1005,7 +764,7 @@ def test_linearize_expr_2d_2():
     print(a)
 
     expected = linearize(l1, F, trials='u1')
-    assert( linearize(l, F, trials='u1') == expected )
+    assert linearize(l, F, trials='u1') == expected
     # ...
 
 #==============================================================================
@@ -1149,16 +908,12 @@ def test_linearize_form_2d_4():
     domain  = Domain('Omega', dim=2)
     Gamma_N = Boundary(r'\Gamma_N', domain)
 
-    x,y = domain.coordinates
-
     V = ScalarFunctionSpace('V', domain)
-
     v  = element_of(V, name='v')
     u  = element_of(V, name='u')
     du = element_of(V, name='du')
 
     int_0 = lambda expr: integral(domain , expr)
-    int_1 = lambda expr: integral(Gamma_N, expr)
 
 #    g = Matrix((cos(pi*x)*sin(pi*y),
 #              sin(pi*x)*cos(pi*y)))
@@ -1176,18 +931,15 @@ def test_linearize_form_2d_4():
 def test_area_2d_1():
 
     domain = Domain('Omega', dim=2)
-    x,y = domain.coordinates
-
-    mu    = Constant('mu'   , is_real=True)
 
     e = ElementDomain()
     area = Area(e)
 
     V = ScalarFunctionSpace('V', domain)
 
-    u,v = [element_of(V, name=i) for i in ['u', 'v']]
+    u, v = elements_of(V, names=['u', 'v'])
 
-    int_0 = lambda expr: integral(domain , expr)
+    int_0 = lambda expr: integral(domain, expr)
 
     # ...
     a = BilinearForm((v,u), int_0(area * u * v))
@@ -1286,9 +1038,6 @@ def test_user_function_2d_1():
     domain = Domain('Omega', dim=2)
     x,y = domain.coordinates
 
-    kappa = Constant('kappa', is_real=True)
-    mu    = Constant('mu'   , is_real=True)
-
     # right hand side
     f = Function('f')
 
@@ -1322,13 +1071,10 @@ def test_functional_2d_1():
     domain = Domain('Omega', dim=2)
     x,y = domain.coordinates
 
-    kappa = Constant('kappa', is_real=True)
-    mu    = Constant('mu'   , is_real=True)
-
     V = ScalarFunctionSpace('V', domain)
     F = element_of(V, name='F')
 
-    int_0 = lambda expr: integral(domain , expr)
+    int_0 = lambda expr: integral(domain, expr)
 
     # ...
     expr = x*y
@@ -1443,29 +1189,20 @@ def test_norm_2d_2():
 def test_bilinear_form_2d_3():
 
     domain = Domain('Omega', dim=2)
-    B1 = Boundary(r'\Gamma_1', domain)
-
-    x,y = domain.coordinates
-
-    kappa = Constant('kappa', is_real=True)
-    mu    = Constant('mu'   , is_real=True)
-
     V = ScalarFunctionSpace('V', domain)
+    u = element_of(V, name='u')
+    v = element_of(V, name='v')
 
-    u,u1,u2 = [element_of(V, name=i) for i in ['u', 'u1', 'u2']]
-    v,v1,v2 = [element_of(V, name=i) for i in ['v', 'v1', 'v2']]
-
-    int_0 = lambda expr: integral(domain , expr)
-    int_1 = lambda expr: integral(B1, expr)
+    int_0 = lambda expr: integral(domain, expr)
 
     # ...
     a = BilinearForm((u,v), int_0(u*v))
-    assert(a.is_symmetric)
+    assert a.is_symmetric
     # ...
 
     # ...
     a = BilinearForm((u,v), int_0(dot(grad(u), grad(v))))
-    assert(a.is_symmetric)
+    assert a.is_symmetric
     # ...
 
 #==============================================================================
@@ -1473,28 +1210,20 @@ def test_bilinear_form_2d_3():
 def test_bilinear_form_2d_4():
 
     domain = Domain('Omega', dim=2)
-    B1 = Boundary(r'\Gamma_1', domain)
-
-    x,y = domain.coordinates
-
-    kappa = Constant('kappa', is_real=True)
-    mu    = Constant('mu'   , is_real=True)
-
     V = VectorFunctionSpace('V', domain)
+    u = element_of(V, name='u')
+    v = element_of(V, name='v')
 
-    u,u1,u2 = [element_of(V, name=i) for i in ['u', 'u1', 'u2']]
-    v,v1,v2 = [element_of(V, name=i) for i in ['v', 'v1', 'v2']]
+    int_0 = lambda expr: integral(domain, expr)
 
-    int_0 = lambda expr: integral(domain , expr)
-    int_1 = lambda expr: integral(B1, expr)
     # ...
-    a = BilinearForm((u,v), int_0(dot(u,v)))
-    assert(a.is_symmetric)
+    a = BilinearForm((u, v), int_0(dot(u, v)))
+    assert a.is_symmetric
     # ...
 
     # ...
-    a = BilinearForm((u,v), int_0(inner(grad(u), grad(v))))
-    assert(a.is_symmetric)
+    a = BilinearForm((u, v), int_0(inner(grad(u), grad(v))))
+    assert a.is_symmetric
     # ...
 
 #==============================================================================
@@ -1656,53 +1385,74 @@ def test_linearity_bilinear_form_2d_1():
 #==============================================================================
 def test_interface_2d_1():
 
+#    # ...
+#    def two_patches():
+#
+#        from sympde.topology import Connectivity, Interface
+#
+#        A = Square('A')
+#        B = Square('B')
+#
+#        A = A.interior
+#        B = B.interior
+#
+#        bnd_A_1 = Boundary(r'\Gamma_1', A, axis=0, ext=-1)
+#        bnd_A_2 = Boundary(r'\Gamma_2', A, axis=0, ext=1)
+#        bnd_A_3 = Boundary(r'\Gamma_3', A, axis=1, ext=-1)
+#        bnd_A_4 = Boundary(r'\Gamma_4', A, axis=1, ext=1)
+#
+#        bnd_B_1 = Boundary(r'\Gamma_1', B, axis=0, ext=-1)
+#        bnd_B_2 = Boundary(r'\Gamma_2', B, axis=0, ext=1)
+#        bnd_B_3 = Boundary(r'\Gamma_3', B, axis=1, ext=-1)
+#        bnd_B_4 = Boundary(r'\Gamma_4', B, axis=1, ext=1)
+#
+#        connectivity = Connectivity()
+#        connectivity['I'] = Interface('I', bnd_A_2, bnd_B_1, ornt=1)
+#        Omega = Domain('Omega',
+#                       interiors=[A, B],
+#                       boundaries=[bnd_A_1, bnd_A_2, bnd_A_3, bnd_A_4, bnd_B_1, bnd_B_2, bnd_B_3, bnd_B_4],
+#                       connectivity=connectivity)
+#
+#        return Omega
+#    # ...
+
     # ...
     def two_patches():
 
-        from sympde.topology import InteriorDomain
-        from sympde.topology import Connectivity, Interface
-
         A = Square('A')
         B = Square('B')
 
-        A = A.interior
-        B = B.interior
-
-        connectivity = Connectivity()
-
-        bnd_A_1 = Boundary(r'\Gamma_1', A, axis=0, ext=-1)
-        bnd_A_2 = Boundary(r'\Gamma_2', A, axis=0, ext=1)
-        bnd_A_3 = Boundary(r'\Gamma_3', A, axis=1, ext=-1)
-        bnd_A_4 = Boundary(r'\Gamma_4', A, axis=1, ext=1)
-
-        bnd_B_1 = Boundary(r'\Gamma_1', B, axis=0, ext=-1)
-        bnd_B_2 = Boundary(r'\Gamma_2', B, axis=0, ext=1)
-        bnd_B_3 = Boundary(r'\Gamma_3', B, axis=1, ext=-1)
-        bnd_B_4 = Boundary(r'\Gamma_4', B, axis=1, ext=1)
-
-        connectivity['I'] = Interface('I', bnd_A_2, bnd_B_1)
-
-        Omega = Domain('Omega',
-                       interiors=[A, B],
-                       boundaries=[bnd_A_1, bnd_A_2, bnd_A_3, bnd_A_4, bnd_B_1, bnd_B_2, bnd_B_3, bnd_B_4],
-                       connectivity=connectivity)
-
-        return Omega
+        patches = [A, B]
+        connectivity = [((A, 0, 1), (B, 0, -1), 1)]
+        return Domain.join(patches, connectivity, 'Omega')
     # ...
 
     # create a domain with an interface
     domain = two_patches()
-    interfaces = domain.interfaces
+    interface = domain.interfaces
+
+    # ...
+    # Make as many checks as possible on the Interface object
+    from sympde.topology.basic import Interface
+    assert isinstance(interface, Interface)
+    assert interface.minus is domain.subdomains[0].get_boundary(axis=0, ext= 1)
+    assert interface.plus  is domain.subdomains[1].get_boundary(axis=0, ext=-1)
+    assert interface.ornt == 1
+    assert interface.dim  == 2
+    assert interface.mapping is None
+    assert interface.logical_domain is None
+    # ...
 
     V = ScalarFunctionSpace('V', domain)
+    assert V.is_broken
 
-    u,v = elements_of(V, names='u, v')
+    u, v = elements_of(V, names='u, v')
 
-    print(integral(interfaces, u*v))
+    print(integral(interface, u*v))
 
     expr  = integral(domain, dot(grad(v),grad(u)))
-    expr += integral(interfaces, - avg(Dn(u)) * jump(v)
-                                 + avg(Dn(v)) * jump(u))
+    expr += integral(interface, - avg(Dn(u)) * jump(v)
+                                + avg(Dn(v)) * jump(u))
     a  = BilinearForm((u,v), expr)
     print(a)
 
@@ -1714,7 +1464,7 @@ def test_interface_integral_1():
     B = Square('B')
 
     domains = [A, B]
-    connectivity = [((0, 0, 1),(1, 0, -1))]
+    connectivity = [((0, 0, 1), (1, 0, -1), 1)]
     domain = Domain.join(domains, connectivity, 'domain')
     # ...
 
@@ -1778,15 +1528,15 @@ def test_interface_integral_2():
     B = Square('B')
 
 
-    domains = [A, B]
-    connectivity = [((0, 0, 1),(1, 0, -1))]
-    domain = Domain.join(domains, connectivity, 'domain')
+    patches = [A, B]
+    connectivity = [((0, 0, 1), (1, 0, -1), 1)]
+    domain = Domain.join(patches, connectivity, 'domain')
     # ...
 
     x,y = domain.coordinates
 
     V = ScalarFunctionSpace('V', domain, kind=None)
-    assert(V.is_broken)
+    assert V.is_broken
 
     u, u1, u2, u3 = elements_of(V, names='u, u1, u2, u3')
     v, v1, v2, v3 = elements_of(V, names='v, v1, v2, v3')
@@ -1821,9 +1571,10 @@ def test_interface_integral_3():
     C = Square('C')
 
 
-    domains = [A, B, C]
-    connectivity = [((0, 0, 1),(1, 0, -1)),((1,0,1),(2,0,-1))]
-    domain = Domain.join(domains, connectivity, 'domain')
+    patches = [A, B, C]
+    connectivity = [((0, 0, 1), (1, 0, -1), 1),
+                    ((1, 0, 1), (2, 0, -1), 1)]
+    domain = Domain.join(patches, connectivity, 'domain')
 
     x,y = domain.coordinates
 
@@ -1863,16 +1614,16 @@ def test_interface_integral_4():
     A = Square('A')
     B = Square('B')
 
-    domains = [A, B]
-    connectivity = [((0, 0, 1),(1, 0, -1))]
-    domain = Domain.join(domains, connectivity, 'AB')
+    patches = [A, B]
+    connectivity = [((0, 0, 1), (1, 0, -1), 1)]
+    domain = Domain.join(patches, connectivity, 'AB')
 
     x,y = domain.coordinates
     assert all([x.name == 'x1', y.name == 'x2'])
 
     V1 = ScalarFunctionSpace('V1', domain, kind=None)
     V2 = VectorFunctionSpace('V2', domain, kind=None)
-    assert(V1.is_broken)
+    assert V1.is_broken
 
     u1, v1 = elements_of(V1, names='u1, v1')
     u2, v2 = elements_of(V2, names='u2, v2')
diff --git a/sympde/topology/basic.py b/sympde/topology/basic.py
index bce622f3..c572ea45 100644
--- a/sympde/topology/basic.py
+++ b/sympde/topology/basic.py
@@ -439,26 +439,85 @@ class Interface(BasicDomain):
     """
     Represents an interface between two subdomains through two boundaries.
 
-    Examples
-
+    Parameters
+    ----------
+    name : str
+        Name of the interface.
+
+    bnd_minus : Boundary
+        Boundary on the "minus" side of the interface.
+
+    bnd_plus : Boundary
+        Boundary on the "plus" side of the interface.
+
+    mapping : Mapping, optional
+        Mapping from the logical domain to the physical domain, if available.
+
+    logical_domain : BasicDomain, optional
+        Logical domain associated with the interface, if available. It should
+        be consistent with the mapping if provided.
+
+    ornt : int | Iterable[int], optional
+        Orientation of the interface. For 1D interfaces, this is not needed and
+        should be set to None. For 2D interfaces, this should be either -1 or 1.
+        For 3D interfaces, this should be a tuple of three integers, each being
+        either -1 or 1.
+
+    Notes
+    -----
+    The orientations are specified in the same manner as in GeoPDES, see e.g.
+    <https://github.com/rafavzqz/geopdes/blob/master/geopdes/doc/geo_specs_mp_v21.txt#L193-L237>
+    and
+    T. Dokken, E. Quak, V. Skytt. Requirements from Isogeometric Analysis for changes in product design ontologies, 2010.
     """
-    def __new__(cls, name, bnd_minus, bnd_plus, mapping=None, logical_domain=None, ornt=None):
+    def __new__(cls, name, bnd_minus, bnd_plus, *, mapping=None, logical_domain=None, ornt=None):
 
-        if not isinstance(name     , str    ): raise TypeError(name)
+        if not isinstance(name     , str     ): raise TypeError(name)
         if not isinstance(bnd_minus, Boundary): raise TypeError(bnd_minus)
         if not isinstance(bnd_plus , Boundary): raise TypeError(bnd_plus)
 
+        # Check that the dimensions of the two boundaries are the same
         if bnd_minus.dim != bnd_plus.dim:
-            raise TypeError('Dimension mismatch: {} != {}'.format(
-                bnd_minus.dim, bnd_plus.dim))
+            raise ValueError(f'Dimension mismatch between boundaries: {bnd_minus.dim} != {bnd_plus.dim}')
+        else:
+            # Number of logical dimensions of the interface
+            ldim = bnd_minus.dim
 
+        # Mapping and logical domain must be provided together or not at all
         assert mapping is None and logical_domain is None or\
                mapping is not None and logical_domain is not None
 
-        if bnd_minus.dim == 3:
+        # If provided, check that mapping is consistent with the boundaries
+        if mapping is not None:
+            from sympde.topology.mapping import Mapping
+            if not isinstance(mapping, Mapping):
+                raise TypeError(f'mapping must be of type Mapping, got {type(mapping)} instead')
+            if mapping.ldim != ldim:
+                raise ValueError(f'mapping should have logical dimension = {ldim}, got {mapping.ldim} instead')
+
+        # If provided, check that logical domain is consistent with the boundaries
+        if logical_domain is not None:
+            if not isinstance(logical_domain, BasicDomain):
+                raise TypeError(f'logical_domain must be of type BasicDomain, got {type(logical_domain)} instead')
+            if logical_domain.dim != ldim:
+                raise ValueError(f'Logical domain should have dimension = {ldim}, got {logical_domain.dim} instead')
+
+        # Check that orientation is provided in the correct format depending on the dimension
+        if ldim == 1:
+            assert ornt is None, 'ornt is not needed for 1D interfaces'
+        elif ldim == 2:
+            assert ornt in (-1, 1), 'ornt must be either -1 or 1 for 2D interfaces'
+        elif ldim == 3:
             ornt = tuple(ornt)
+            assert len(ornt) == 3, 'ornt must be a tuple of length 3 for 3D interfaces'
+            assert all(o in (-1, 1) for o in ornt), 'each element of ornt must be either -1 or 1 for 3D interfaces'
+        else:
+            raise ValueError(f'Unsupported dimension: {ldim}')
 
+        # Strong requirement: the two boundaries must be defined over the same axis
+        # TODO [YG 10.02.2026]: relax this requirement ASAP
         assert bnd_minus.axis == bnd_plus.axis
+
         obj = Basic.__new__(cls, name, bnd_minus, bnd_plus, ornt)
         obj._mapping        = mapping
         obj._logical_domain = logical_domain
@@ -520,10 +579,10 @@ def __init__(self, data=None):
         if data is None:
             data = {}
         else:
-            assert isinstance( data, dict )
-            for k,v in data.items():
-                assert( isinstance( k, str ) )
-                assert( isinstance(v, Interface) )
+            assert isinstance(data, dict)
+            for k, v in data.items():
+                assert isinstance(k, str)
+                assert isinstance(v, Interface)
         self._data = data
 
     @property
@@ -543,7 +602,7 @@ def todict(self):
         connectivity = {}
         data = dict(sorted(self._data.items()))
         for name, v in data.items():
-            connectivity[name] = [v.minus.todict(), v.plus.todict()]
+            connectivity[name] = [v.minus.todict(), v.plus.todict(), v.ornt]
         connectivity = dict(sorted(connectivity.items()))
         # ...
 
diff --git a/sympde/topology/domain.py b/sympde/topology/domain.py
index 12627ac8..2dde676e 100644
--- a/sympde/topology/domain.py
+++ b/sympde/topology/domain.py
@@ -351,7 +351,7 @@ def from_file( cls, filename ):
         if not(ext == '.h5'):
             raise ValueError('> Only h5 files are supported')
         # ...
-        from sympde.topology.mapping import Mapping, MultiPatchMapping
+        from sympde.topology.mapping import Mapping
 
         h5  = h5py.File( filename, mode='r' )
         yml = yaml.load( h5['topology.yml'][()], Loader=yaml.SafeLoader )
@@ -387,7 +387,7 @@ def from_file( cls, filename ):
             boundaries.append(bd)
 
         connectivity = []
-        for _,(minus, plus) in d_connectivity.items():
+        for _,(minus, plus, ornt) in d_connectivity.items():
             minus_name = minus['patch']
             minus_axis = int(minus['axis'])
             minus_ext  = int(minus['ext'])
@@ -397,7 +397,8 @@ def from_file( cls, filename ):
             plus_axis = int(plus['axis'])
             plus_ext  = int(plus['ext'])
             plus_patch_i = patch_index[plus_name]
-            interface = ((minus_patch_i, minus_axis, minus_ext),(plus_patch_i, plus_axis, plus_ext))
+            interface = ((minus_patch_i, minus_axis, minus_ext),
+                         ( plus_patch_i,  plus_axis,  plus_ext), ornt)
 
             connectivity.append(interface)
 
@@ -488,9 +489,7 @@ def join(cls, patches, connectivity, name):
 
         # the multi-patch domain
         Omega = Domain.join(patches=patches, connectivity=connectivity, name='Omega')
-        
         """
-
         assert isinstance(patches, (tuple, list))
         assert isinstance(connectivity, (tuple, list))
         assert isinstance(name, str)
@@ -501,7 +500,7 @@ def join(cls, patches, connectivity, name):
             return patches[0]
 
         assert all(p.dim==patches[0].dim for p in patches)
-        dim = int(patches[0].dim)
+        ldim = int(patches[0].dim)
 
         patch_given_by_indices = (len(connectivity) > 0 and isinstance(connectivity[0][0][0], int))
 
@@ -518,14 +517,21 @@ def join(cls, patches, connectivity, name):
                 patch_plus  = cn[1][0]
             bnd_minus = patch_minus.get_boundary(axis=cn[0][1], ext=cn[0][2])
             bnd_plus  =  patch_plus.get_boundary(axis=cn[1][1], ext=cn[1][2])
-            if   dim == 1: ornt = None
-            elif dim == 2: ornt = 1 if len(cn) == 2 else cn[2]
-            elif dim == 3:
-                flag = 1 if len(cn) == 2 else cn[2][0]
-                ornt1 = 1 if len(cn) == 2 else cn[2][1]
-                ornt2 = 1 if len(cn) == 2 else cn[2][2]
-                ornt = (flag, ornt1, ornt2)
 
+            # Check that orientation is provided in the correct format depending on the dimension
+            ornt = cn[2] if len(cn) == 3 else None
+            if ldim == 1:
+                assert ornt is None, 'ornt is not needed for 1D interfaces'
+            elif ldim == 2:
+                assert ornt in (-1, 1), 'ornt must be either -1 or 1 for 2D interfaces'
+            elif ldim == 3:
+                ornt = tuple(ornt)
+                assert len(ornt) == 3, 'ornt must be a tuple of length 3 for 3D interfaces'
+                assert all(o in (-1, 1) for o in ornt), 'each element of ornt must be either -1 or 1 for 3D interfaces'
+            else:
+                raise ValueError(f'Unsupported dimension: {ldim}')
+
+            # Create a new Interface object using the join method of the Boundary objects
             interface = bnd_minus.join(bnd_plus, ornt=ornt)
             if interface.name in interfaces:
                 interface = bnd_plus.join(bnd_minus, ornt=ornt)
diff --git a/sympde/topology/mapping.py b/sympde/topology/mapping.py
index 6efa5c90..c4844316 100644
--- a/sympde/topology/mapping.py
+++ b/sympde/topology/mapping.py
@@ -54,6 +54,7 @@
     'JacobianSymbol',
     'LogicalExpr',
     'MappedDomain',
+    'Mapping',
     'MappingApplication',
     'MultiPatchMapping',
     'PullBack',
@@ -318,8 +319,11 @@ def logical_coordinates( self ):
         else:
             return self._logical_coordinates
 
+    # Applying the mapping to a logical domain returns a mapped domain
     def __call__(self, domain):
-        assert(isinstance(domain, BasicDomain))
+        assert isinstance(domain, BasicDomain)
+        assert domain.logical_domain is None
+        assert domain.dim == self.ldim
         return MappedDomain(self, domain)
 
     @property
diff --git a/sympde/topology/tests/test_logical_expr.py b/sympde/topology/tests/test_logical_expr.py
index 59702e84..b66d3d86 100644
--- a/sympde/topology/tests/test_logical_expr.py
+++ b/sympde/topology/tests/test_logical_expr.py
@@ -221,9 +221,9 @@ def test_logical_expr_2d_2():
     D1 = M1(A)
     D2 = M2(B)
 
-    domains = [D1, D2]
-    connectivity = [((0, 0, 1),(1, 0, -1))]
-    domain = Domain.join(domains, connectivity, 'domain')
+    patches = [D1, D2]
+    connectivity = [((0, 0, 1), (1, 0, -1), 1)]
+    domain = Domain.join(patches, connectivity, 'domain')
 
 
     V1 = ScalarFunctionSpace('V1', domain, kind='h1')
@@ -275,11 +275,12 @@ def test_logical_expr_2d_3():
     D2 = M2(B)
 
 
-    domains = [D1, D2]
-    connectivity = [((0, 0, 1),(1, 0, -1))]
-    domain = Domain.join(domains, connectivity, 'domain')
+    patches = [D1, D2]
+    connectivity = [((0, 0, 1), (1, 0, -1), 1)]
+    domain = Domain.join(patches, connectivity, 'domain')
 
     V = VectorFunctionSpace('V', domain, kind='hcurl')
+    assert V.is_broken
 
     u,v = [element_of(V, name=i) for i in ['u', 'v']]
 
diff --git a/sympde/topology/tests/test_topology.py b/sympde/topology/tests/test_topology.py
index adf2cd7b..aa3e3dda 100644
--- a/sympde/topology/tests/test_topology.py
+++ b/sympde/topology/tests/test_topology.py
@@ -1,39 +1,32 @@
-# coding: utf-8
-
-from sympy.tensor import Indexed
+import os
 
 from sympde.topology import InteriorDomain, Union
-from sympde.topology import Boundary, NormalVector, TangentVector
-from sympde.topology import Connectivity, Edge
+from sympde.topology import Boundary
 from sympde.topology import Domain, ElementDomain
 from sympde.topology import Area, Mapping
 from sympde.topology import Interface
 from sympde.topology import Line, Square, Cube
 from sympde.topology import IdentityMapping
 
-import os
-
 base_dir = os.path.dirname(os.path.realpath(__file__))
 topo_dir = os.path.join(base_dir, 'data')
 
-
 #==============================================================================
 def test_interior_domain():
     D1 = InteriorDomain('D1', dim=2)
     D2 = InteriorDomain('D2', dim=2)
 
-    assert( D1.todict() == {'name': 'D1', 'mapping':'None'} )
-    assert( D2.todict() == {'name': 'D2', 'mapping':'None'} )
+    assert D1.todict() == {'name': 'D1', 'mapping':'None'}
+    assert D2.todict() == {'name': 'D2', 'mapping':'None'}
 
-    assert( Union(D2, D1) == Union(D1, D2) )
+    assert Union(D2, D1) == Union(D1, D2)
 
     D = Union(D1, D2)
 
-    assert(D.dim == 2)
-    assert(len(D) == 2)
-    assert( D.todict() == [{'name': 'D1', 'mapping':'None'},
-                           {'name': 'D2', 'mapping':'None'}] )
-
+    assert D.dim == 2
+    assert len(D) == 2
+    assert D.todict() == [{'name': 'D1', 'mapping':'None'},
+                          {'name': 'D2', 'mapping':'None'}]
 
 #==============================================================================
 def test_topology_1():
@@ -49,12 +42,12 @@ def test_topology_1():
     D1 = M1(A)
     D2 = M2(B)
 
-    domains = [D1, D2]
-    connectivity = [((0, 0, 1), (1, 0, -1))]
-    Omega = Domain.join(domains, connectivity, 'domain')
+    patches = [D1, D2]
+    connectivity = [((0, 0, 1), (1, 0, -1), 1)]
+    Omega = Domain.join(patches, connectivity, 'domain')
 
     interfaces = Omega.interfaces
-    assert(isinstance(interfaces, Interface))
+    assert isinstance(interfaces, Interface)
 
     # export
     Omega.export('omega.h5')
@@ -63,7 +56,7 @@ def test_topology_1():
     # read it again and check that it has the same description as Omega
     D = Domain.from_file('omega.h5')
 
-    assert( D.todict() == Omega.todict() )
+    assert D.todict() == Omega.todict()
 
 #==============================================================================
 def test_domain_1():
@@ -78,9 +71,9 @@ def test_domain_1():
                    interiors=[Omega_1, Omega_2],
                    boundaries=[Gamma_1, Gamma_2, Gamma_3])
 
-    assert( Omega.dim == 2 )
-    assert( len(Omega.interior) == 2 )
-    assert( len(Omega.boundary) == 3 )
+    assert Omega.dim == 2
+    assert len(Omega.interior) == 2
+    assert len(Omega.boundary) == 3
 
 #==============================================================================
 def test_boundary_1():
@@ -93,9 +86,9 @@ def test_boundary_1():
                    interiors=[Omega_1],
                    boundaries=[Gamma_1, Gamma_2])
 
-    assert(Omega.boundary == Union(Gamma_1, Gamma_2))
-    assert(Omega.boundary.complement(Gamma_1) == Gamma_2)
-    assert(Omega.boundary - Gamma_1 == Gamma_2)
+    assert Omega.boundary == Union(Gamma_1, Gamma_2)
+    assert Omega.boundary.complement(Gamma_1) == Gamma_2
+    assert Omega.boundary - Gamma_1 == Gamma_2
 
 #==============================================================================
 def test_boundary_2():
@@ -109,23 +102,23 @@ def test_boundary_2():
                    interiors=[Omega_1],
                    boundaries=[Gamma_1, Gamma_2, Gamma_3])
 
-    assert(Omega.boundary == Union(Gamma_1, Gamma_2, Gamma_3))
-    assert(Omega.boundary.complement(Gamma_1) == Union(Gamma_2, Gamma_3))
-    assert(Omega.boundary - Gamma_1 == Union(Gamma_2, Gamma_3))
+    assert Omega.boundary == Union(Gamma_1, Gamma_2, Gamma_3)
+    assert Omega.boundary.complement(Gamma_1) == Union(Gamma_2, Gamma_3)
+    assert Omega.boundary - Gamma_1 == Union(Gamma_2, Gamma_3)
 
 #==============================================================================
 def test_boundary_3():
     Omega_1 = InteriorDomain('Omega_1', dim=2)
 
     Gamma_1 = Boundary(r'\Gamma_1', Omega_1, axis=0, ext=-1)
-    Gamma_4 = Boundary(r'\Gamma_4', Omega_1, axis=1, ext=1)
+    Gamma_4 = Boundary(r'\Gamma_4', Omega_1, axis=1, ext= 1)
 
     Omega = Domain('Omega',
                    interiors=[Omega_1],
                    boundaries=[Gamma_1, Gamma_4])
 
-    assert(Omega.get_boundary(axis=0, ext=-1) == Gamma_1)
-    assert(Omega.get_boundary(axis=1, ext=1) == Gamma_4)
+    assert Omega.get_boundary(axis=0, ext=-1) == Gamma_1
+    assert Omega.get_boundary(axis=1, ext= 1) == Gamma_4
 
 #==============================================================================
 def test_element():
@@ -142,7 +135,7 @@ def test_element():
     a = Area(D1)
     print(a)
 
-    assert(Area(D) ==  Area(D1) + Area(D2))
+    assert Area(D) ==  Area(D1) + Area(D2)
 
 #==============================================================================
 def test_domain_join_line():
@@ -178,7 +171,8 @@ def test_domain_join_line():
     ABC = Domain.join(domains, connectivity, 'ABC')
 
     assert ABC.interior == Union(A.interior, B.interior, C.interior)
-    assert ABC.interfaces == Union(Interface('A|B', AB_bnd_minus, AB_bnd_plus),Interface('B|C', BC_bnd_minus, BC_bnd_plus))
+    assert ABC.interfaces == Union(Interface('A|B', AB_bnd_minus, AB_bnd_plus),
+                                   Interface('B|C', BC_bnd_minus, BC_bnd_plus))
     print(list(ABC.connectivity.items()))
     print('')
     # ...
@@ -193,14 +187,14 @@ def test_domain_join_square():
     # ...
 
 
-    domains = [A, B]
-    connectivity = [((0, 0, 1), (1, 0, -1))]
-    AB = Domain.join(domains, connectivity, 'AB')
+    patches = [A, B]
+    connectivity = [((0, 0, 1), (1, 0, -1), 1)]
+    AB = Domain.join(patches, connectivity, 'AB')
 
-    AB_bnd_minus = A.get_boundary(axis=0, ext=1)
+    AB_bnd_minus = A.get_boundary(axis=0, ext= 1)
     AB_bnd_plus  = B.get_boundary(axis=0, ext=-1)
 
-    BC_bnd_minus = B.get_boundary(axis=0, ext=1)
+    BC_bnd_minus = B.get_boundary(axis=0, ext= 1)
     BC_bnd_plus  = C.get_boundary(axis=0, ext=-1)
 
     print(AB)
@@ -209,9 +203,10 @@ def test_domain_join_square():
     print(AB.connectivity)
     # ...
 
-    domains = [A, B, C]
-    connectivity = [((0, 0, 1),(1, 0, -1)), ((1, 0, 1), (2, 0, -1))]
-    ABC = Domain.join(domains, connectivity, 'ABC')
+    patches = [A, B, C]
+    connectivity = [((0, 0, 1), (1, 0, -1), 1),
+                    ((1, 0, 1), (2, 0, -1), 1)]
+    ABC = Domain.join(patches, connectivity, 'ABC')
 
 
     print(ABC)
@@ -229,15 +224,15 @@ def test_get_subdomain():
     C = Square('C')
     # ...
 
-    domains = [A, B]
-    connectivity = [((0, 0, 1), (1, 0, -1))]
-    AB = Domain.join(domains, connectivity, 'AB')
+    patches = [A, B]
+    connectivity = [((0, 0, 1), (1, 0, -1), 1)]
+    AB = Domain.join(patches, connectivity, 'AB')
 
     # ...
 
-    domains = [A, B, C]
-    connectivity = [((0, 0, 1), (1, 0, -1)), ((1, 0, 1), (2, 0, -1))]
-    ABC = Domain.join(domains, connectivity, 'ABC')
+    patches = [A, B, C]
+    connectivity = [((0, 0, 1), (1, 0, -1), 1), ((1, 0, 1), (2, 0, -1), 1)]
+    ABC = Domain.join(patches, connectivity, 'ABC')
 
     A_1 = AB.get_subdomain('A')
     A_2 = ABC.get_subdomain('A')
@@ -261,62 +256,62 @@ def test_get_subdomain():
 def test_2d_domain_without_bnd():
 
     OmegaLog1 = Square('OmegaLog1', bounds1 = (0,.5), bounds2 = (0,.5))
-    mapping_1 = IdentityMapping('M1',2)
-    domain_1     = mapping_1(OmegaLog1)
+    mapping_1 = IdentityMapping('M1', 2)
+    domain_1  = mapping_1(OmegaLog1)
     OmegaLog2 = Square('OmegaLog2', bounds1 = (0,.5), bounds2 = (.5,1.))
-    mapping_2 = IdentityMapping('M2',2)
-    domain_2     = mapping_2(OmegaLog2)
+    mapping_2 = IdentityMapping('M2', 2)
+    domain_2  = mapping_2(OmegaLog2)
     OmegaLog3 = Square('OmegaLog3', bounds1 = (.5,1.), bounds2 = (0,.5))
-    mapping_3 = IdentityMapping('M3',2)
-    domain_3     = mapping_3(OmegaLog3)
+    mapping_3 = IdentityMapping('M3', 2)
+    domain_3  = mapping_3(OmegaLog3)
     OmegaLog4 = Square('OmegaLog4', bounds1 = (.5,1.), bounds2 = (.5,1.))
-    mapping_4 = IdentityMapping('M4',2)
-    domain_4     = mapping_4(OmegaLog4)
-
-    domains=[domain_1,domain_2,domain_3,domain_4]
-
-    connectivity = [((0, 0, 1), (2, 0,-1)),
-                    ((1, 0, 1), (3, 0,-1)),
-                    ((2, 0, 1), (0, 0,-1)),
-                    ((3, 0, 1), (1, 0,-1)),
-                    ((0, 1, 1), (1, 1,-1)),
-                    ((2, 1, 1), (3, 1,-1)),
-                    ((1, 1, 1), (0, 1,-1)),
-                    ((3, 1, 1), (2, 1,-1))]
-    domain = Domain.join(domains, connectivity, 'domain')
+    mapping_4 = IdentityMapping('M4', 2)
+    domain_4  = mapping_4(OmegaLog4)
+
+    patches =  [domain_1, domain_2, domain_3, domain_4]
+    connectivity = [((0, 0, 1), (2, 0,-1), 1),
+                    ((1, 0, 1), (3, 0,-1), 1),
+                    ((2, 0, 1), (0, 0,-1), 1),
+                    ((3, 0, 1), (1, 0,-1), 1),
+                    ((0, 1, 1), (1, 1,-1), 1),
+                    ((2, 1, 1), (3, 1,-1), 1),
+                    ((1, 1, 1), (0, 1,-1), 1),
+                    ((3, 1, 1), (2, 1,-1), 1)]
+    domain = Domain.join(patches, connectivity, 'domain')
 
     assert len(domain.interior) == 4
     assert len(domain.interfaces) == 8
 
+#==============================================================================
 def test_3d_domain_without_bnd():
 
     OmegaLog1 = Cube('OmegaLog1', bounds1 = (0,.5), bounds2 = (0,.5), bounds3 = (0,1))
-    mapping_1 = IdentityMapping('M1',2)
-    domain_1     = mapping_1(OmegaLog1)
+    mapping_1 = IdentityMapping('M1', 3)
+    domain_1  = mapping_1(OmegaLog1)
     OmegaLog2 = Cube('OmegaLog2', bounds1 = (0,.5), bounds2 = (.5,1.), bounds3 = (0,1))
-    mapping_2 = IdentityMapping('M2',2)
-    domain_2     = mapping_2(OmegaLog2)
+    mapping_2 = IdentityMapping('M2', 3)
+    domain_2  = mapping_2(OmegaLog2)
     OmegaLog3 = Cube('OmegaLog3', bounds1 = (.5,1.), bounds2 = (0,.5), bounds3 = (0,1))
-    mapping_3 = IdentityMapping('M3',2)
-    domain_3     = mapping_3(OmegaLog3)
+    mapping_3 = IdentityMapping('M3', 3)
+    domain_3  = mapping_3(OmegaLog3)
     OmegaLog4 = Cube('OmegaLog4', bounds1 = (.5,1.), bounds2 = (.5,1.), bounds3 = (0,1))
-    mapping_4 = IdentityMapping('M4',2)
-    domain_4     = mapping_4(OmegaLog4)
-
-    domains=[domain_1,domain_2,domain_3,domain_4]
-
-    connectivity = [((0, 0, 1),(2, 0, -1),(1,1,1)),
-                    ((1,0,1),(3,0,-1),(1,1,1)),
-                    ((2,0,1),(0,0,-1),(1,1,1)),
-                    ((3,0,1),(1,0,-1),(1,1,1)),
-                    ((0,1,1),(1,1,-1),(1,1,1)),
-                    ((2,1,1),(3,1,-1),(1,1,1)),
-                    ((1,1,1),(0,1,-1),(1,1,1)),
-                    ((3,1,1),(2,1,-1),(1,1,1))]
-    domain = Domain.join(domains, connectivity, 'domain')
+    mapping_4 = IdentityMapping('M4', 3)
+    domain_4  = mapping_4(OmegaLog4)
+
+    patches = [domain_1, domain_2, domain_3, domain_4]
+    connectivity = [((0, 0, 1), (2, 0,-1), (1, 1, 1)),
+                    ((1, 0, 1), (3, 0,-1), (1, 1, 1)),
+                    ((2, 0, 1), (0, 0,-1), (1, 1, 1)),
+                    ((3, 0, 1), (1, 0,-1), (1, 1, 1)),
+                    ((0, 1, 1), (1, 1,-1), (1, 1, 1)),
+                    ((2, 1, 1), (3, 1,-1), (1, 1, 1)),
+                    ((1, 1, 1), (0, 1,-1), (1, 1, 1)),
+                    ((3, 1, 1), (2, 1,-1), (1, 1, 1))]
+    domain = Domain.join(patches, connectivity, 'domain')
 
     assert len(domain.interior) == 4
     assert len(domain.interfaces) == 8
+
 #==============================================================================
 def test_hash():
     A = Square('A', bounds1=(0, 1), bounds2=(0, 1))