diff --git a/sympde/core/utils.py b/sympde/core/utils.py
index 7a734dd..ae0fbca 100644
--- a/sympde/core/utils.py
+++ b/sympde/core/utils.py
@@ -1,9 +1,163 @@
 # coding: utf-8
 
+import re
 import string
 import random
+from sympy.utilities.iterables import cartes
 
+__all__ = ('random_string', 'expand_name_patterns')
+
+#==============================================================================
 def random_string( n ):
     chars    = string.ascii_uppercase + string.ascii_lowercase
     selector = random.SystemRandom()
     return ''.join( selector.choice( chars ) for _ in range( n ) )
+
+#==============================================================================
+_range = re.compile('([0-9]*:[0-9]+|[a-zA-Z]?:[a-zA-Z])')
+
+def expand_name_patterns(names, *, seq=None):
+    """
+    Expand names according to the given pattern or list of patterns, following
+    exactly the same rules of Sympy's "symbols" function.
+
+    Parameters
+    ----------
+    names: str or iterable
+        Variable name, or name pattern, or list of name patterns, or list of
+        [list of ...] name patterns.
+
+    seq: bool
+        Set to True if an iterable container is needed for a single name.
+
+    Results
+    -------
+    res: str or iterable
+        Single name string, or list of [list of ...] name strings (having the
+        same type as the input parameter 'names').
+
+    See also
+    --------
+    https://docs.sympy.org/latest/modules/core.html#sympy.core.symbol.symbols
+
+    Notes
+    -----
+    This code is extracted from Sympy's "symbols" function, but it only
+    generates variables' names without creating new objects, which can be
+    instantiated at a later stage. This creates a very clear logical
+    separation between 1) names' generation, and 2) objects' instantiation,
+    plus it allows one to reuse the same name expansion rules in different
+    locations of the library, for generating different objects. Further,
+    since separate functions take on conceptually different tasks, it is much
+    easier to test them individually.
+
+    This approach would also be beneficial to Sympy's "symbols", as it would
+    allow for other Python libraries based on Sympy (like SymPDE) to be
+    compliant with Sympy's name expansion rules without duplicating its code.
+    Hence, in the future we would like to propose the integration of this
+    function into Sympy.
+
+    """
+    result = []
+
+    if isinstance(names, str):
+        marker = 0
+        literals = [r'\,', r'\:', r'\ ']
+        for i in range(len(literals)):
+            lit = literals.pop(0)
+            if lit in names:
+                while chr(marker) in names:
+                    marker += 1
+                lit_char = chr(marker)
+                marker += 1
+                names = names.replace(lit, lit_char)
+                literals.append((lit_char, lit[1:]))
+        def literal(s):
+            if literals:
+                for c, l in literals:
+                    s = s.replace(c, l)
+            return s
+
+        names = names.strip()
+        as_seq = names.endswith(',')
+        if as_seq:
+            names = names[:-1].rstrip()
+        if not names:
+            raise ValueError('no symbols given')
+
+        # split on commas
+        names = [n.strip() for n in names.split(',')]
+        if not all(n for n in names):
+            raise ValueError('missing symbol between commas')
+
+        # split on spaces
+        for i in range(len(names) - 1, -1, -1):
+            names[i: i + 1] = names[i].split()
+
+        # if provided, get 'seq' value from parameter list
+        if seq is None:
+            seq = as_seq
+        else:
+            if not isinstance(seq, bool):
+                msg = 'type(seq) must be bool, got {} instead'.format(type(seq))
+                raise TypeError(msg)
+
+        # expand ranges
+        for name in names:
+            if not name:
+                raise ValueError('missing symbol')
+
+            if ':' not in name:
+                result.append(literal(name))
+                continue
+
+            split = _range.split(name)
+            # remove 1 layer of bounding parentheses around ranges
+            for i in range(len(split) - 1):
+                if i and ':' in split[i] and split[i] != ':' and \
+                        split[i - 1].endswith('(') and \
+                        split[i + 1].startswith(')'):
+                    split[i - 1] = split[i - 1][:-1]
+                    split[i + 1] = split[i + 1][1:]
+            for i, s in enumerate(split):
+                if ':' in s:
+                    if s[-1].endswith(':'):
+                        raise ValueError('missing end range')
+                    a, b = s.split(':')
+                    if b[-1] in string.digits:
+                        a = 0 if not a else int(a)
+                        b = int(b)
+                        split[i] = [str(c) for c in range(a, b)]
+                    else:
+                        a = a or 'a'
+                        split[i] = [string.ascii_letters[c] for c in range(
+                            string.ascii_letters.index(a),
+                            string.ascii_letters.index(b) + 1)]  # inclusive
+                    if not split[i]:
+                        break
+                else:
+                    split[i] = [s]
+            else:
+                seq = True
+                if len(split) == 1:
+                    names = split[0]
+                else:
+                    names = [''.join(s) for s in cartes(*split)]
+                if literals:
+                    result.extend([literal(s) for s in names])
+                else:
+                    result.extend([s for s in names])
+
+        if not seq and len(result) <= 1:
+            if not result:
+                return ()
+            return result[0]
+
+        return tuple(result)
+
+    else:
+        # recursive call on list of names, seq is ignored
+        for name in names:
+            result.append(expand_name_patterns(name))
+
+        return type(names)(result)
diff --git a/sympde/expr/expr.py b/sympde/expr/expr.py
index 1d6fe1f..5b7cd10 100644
--- a/sympde/expr/expr.py
+++ b/sympde/expr/expr.py
@@ -1,71 +1,35 @@
 # coding: utf-8
 
-from itertools import groupby
-
-from sympy.core import Basic
-from sympy.core import Symbol
-from sympy.core import Function
-from sympy.simplify.simplify import simplify
-from sympy import collect
-from sympy.series.order import Order
-from sympy.core import Expr, Add, Mul, Pow
-from sympy import S
 from sympy import Dummy
+from sympy import Matrix
+from sympy.core import Basic
+from sympy.core import Expr, Add, Mul
+from sympy.core.expr import AtomicExpr
 from sympy.core.numbers import Zero as sy_Zero
 from sympy.core.containers import Tuple
-from sympy import Indexed, IndexedBase, Matrix, ImmutableDenseMatrix
-from sympy import Integer, Float
-from sympy.core.expr import AtomicExpr
-from sympy.physics.quantum import TensorProduct
-from sympy.series.series import series
 from sympy.core.compatibility import is_sequence
 
-from sympde.core.basic import _coeffs_registery
 from sympde.core.basic import CalculusFunction
 from sympde.core.basic import Constant
-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.calculus import Dot, Inner, Cross
-from sympde.calculus import Grad, Rot, Curl, Div, Hessian
-from sympde.calculus import Bracket
-from sympde.calculus import Laplace
-from sympde.calculus.core import _generic_ops
-
-from sympde.topology import BasicDomain, Domain, MappedDomain, Union, Interval
-from sympde.topology import BoundaryVector, NormalVector, TangentVector
-from sympde.topology import Boundary, Connectivity, Interface
-from sympde.topology.derivatives import _partial_derivatives
-from sympde.topology.derivatives import partial_derivative_as_symbol
-from sympde.topology.derivatives import sort_partial_derivatives
-from sympde.topology.derivatives import get_atom_derivatives
-from sympde.topology.derivatives import dx, dy, dz
-from sympde.topology.derivatives import (Grad_1d, Div_1d,
-                                         Grad_2d, Curl_2d, Rot_2d, Div_2d,
-                                         Grad_3d, Curl_3d, Div_3d)
-from sympde.topology.derivatives import Bracket_2d
-from sympde.topology.derivatives import Laplace_1d, Laplace_2d, Laplace_3d
-from sympde.topology.derivatives import Hessian_1d, Hessian_2d, Hessian_3d
-from sympde.topology.space import BasicFunctionSpace
+from sympde.calculus import Dot, Inner
+from sympde.calculus import Grad, Hessian
+from sympde.topology import BasicDomain, Union
+from sympde.topology import NormalVector
+from sympde.topology import Boundary, Interface
 from sympde.topology.space import ScalarFunctionSpace
 from sympde.topology.space import ProductSpace
 from sympde.topology.space import ScalarTestFunction
 from sympde.topology.space import VectorTestFunction
-from sympde.topology.space import IndexedTestTrial
-from sympde.topology.space import Unknown, VectorUnknown
 from sympde.topology.space import Trace, trace_0, trace_1
-from sympde.topology.space import ScalarField, VectorField, IndexedVectorField
-from sympde.topology.measure import CanonicalMeasure
-from sympde.topology.measure import CartesianMeasure
-from sympde.topology.measure import Measure
+from sympde.topology.space import ScalarField, VectorField
 
 from .errors import UnconsistentLinearExpressionError
-from .basic import BasicForm
+from .basic  import BasicForm
 from .basic  import BasicExpr
 from .basic  import is_linear_form, _sanitize_arguments
 
+#==============================================================================
 def expand(expr):
     from sympy import expand as _expand
 
@@ -82,6 +46,7 @@ def expand(expr):
         args[i] = c
 
     return Add(*args)
+
 #==============================================================================
 def _get_domain(expr):
     # expr is an integral of BasicExpr or Add of Integral of BasicExpr
@@ -187,7 +152,6 @@ class BilinearExpr(BasicExpr):
         variables = Tuple(*self.variables[0], *self.variables[1])
         return self.expr.xreplace(dict(list(zip(self.variables, args))))
 
-
 #==============================================================================
 class Integral(CalculusFunction):
 
@@ -225,10 +189,7 @@ class Integral(CalculusFunction):
         else:
             return DomainIntegral(expr, domain)
 
-
-
 #==============================================================================
-
 class DomainIntegral(AtomicExpr):
     _op_priority = 20
     @property
@@ -253,6 +214,7 @@ class DomainIntegral(AtomicExpr):
 
     def __hash__(self):
         return hash(self.expr) + hash(self.domain)
+
 #==============================================================================
 class BoundaryIntegral(AtomicExpr):
     _op_priority = 20
@@ -294,7 +256,6 @@ class BoundaryIntegral(AtomicExpr):
 
         return Basic.__new__(cls, expr, domain)
 
-
     @property
     def expr(self):
         return self._args[0]
@@ -346,7 +307,6 @@ class InterfaceIntegral(AtomicExpr):
     def __hash__(self):
         return hash(self.expr) + hash(self.domain)
 
-
 #==============================================================================
 class Functional(BasicForm):
     is_functional = True
@@ -394,7 +354,6 @@ class Functional(BasicForm):
 #    def _eval_nseries(self, x, n, logx):
 #        return self.expr._eval_nseries(x, n, logx)
 
-
 #==============================================================================
 class LinearForm(BasicForm):
     is_linear = True
@@ -476,7 +435,6 @@ class LinearForm(BasicForm):
         # ...
         return expr
 
-
 #==============================================================================
 class BilinearForm(BasicForm):
     is_bilinear = True
@@ -669,103 +627,71 @@ class Norm(Functional):
     def exponent(self):
         return self._exponent
 
-
 #==============================================================================
 def linearize(form, fields, trials=None):
-    """linearize a LinearForm around the fields."""
-    # ...
-    #TODO add LinearForm
-    if not isinstance(form, (LinearForm, LinearExpr)):
+    """
+    Parameters
+    ----------
+    form : LinearForm
+        Linear form F(v; u) to be linearized. F is nonlinear in the parameter u
+
+    fields : [Scalar|Vector]TestFunction or tuple
+        Function u with respect to which F should be differentiated.
+        Tuple if element of ProductSpace.
+
+    trials : [Scalar|Vector]TestFunction or tuple
+        Trial function to be used in resulting bilinear form.
+
+    Results
+    -------
+    a : BilinearForm
+        a(u, v) is bilinear form obtained from linearization of F(v; u).
+
+    """
+    if not isinstance(form, LinearForm):
         raise TypeError('> Expecting a LinearForm')
 
+    # ...
     if not isinstance(fields, (list, tuple, Tuple)):
         fields = [fields]
 
-    fields = [f.space.field(str(f.name)) for f in fields]
-    form   = form.annotate()
-
-    for f in fields:
-
-        if not isinstance(f, (ScalarField, VectorField)):
-            raise TypeError('{} is not ScalarField/VectorField'.format(f))
+    for field in fields:
+        if not isinstance(field, (ScalarTestFunction, VectorTestFunction)):
+            raise TypeError('> Expecting a TestFunction', field)
 
-    if not(trials is None):
+    # ...
+    if trials:
         if not isinstance(trials, (list, tuple, Tuple)):
             trials = [trials]
 
-        assert( all([isinstance(i, (str, ScalarTestFunction, VectorTestFunction)) for i in trials]) )
-        assert( len(fields) == len(trials) )
-
-        newtrials = []
-        for i in trials:
-            if isinstance(i, (ScalarTestFunction, VectorTestFunction)):
-                newtrials += [i.name]
-
-            else:
-                newtrials += [i]
-
-        trials = newtrials
-    # ...
-
-    if isinstance(form, LinearForm):
-        is_form = True
-        expr = form.body
-        domain = form.domain
-
+        for trial in trials:
+            if not isinstance(trial, (ScalarTestFunction, VectorTestFunction)):
+                raise TypeError('> Expecting a TestFunction', trial)
     else:
-        is_form = False
-        expr = form.expr
-
-    test_functions = form.variables
-    fields         = Tuple(*fields)
-
-    # ...
-    trial_functions = []
-    newargs         = []
-    eps  = Constant('eps_' + random_string( 4 ))
-    for i,x in enumerate(fields):
-        tag  = random_string( 4 )
-
-        if trials is None:
-            name = x.name + '_' + tag
-        else:
-            name = trials[i]
-
-        if isinstance(x, ScalarField):
-            trial  = ScalarTestFunction(x.space, name=name)
-
-        elif isinstance(x, VectorField):
-            trial  = VectorTestFunction(x.space, name=name)
+        trials = [x.space.element(x.name + '_trial_' + random_string(4)) for x in fields]
 
-        else:
-            raise TypeError('Only ScalarTestFunction , VectorTestFunction  are available')
-
-        newargs         += [x + eps*trial]
-        trial_functions += [trial]
     # ...
 
-    # ...
-    newexpr = expr
-    for k,v in zip(fields, newargs):
-        newexpr = newexpr.subs(k,v)
-    # ...
+    eps        = Constant('eps_' + random_string(4))
+    new_fields = [field + eps*trial for field, trial in zip(fields, trials)]
 
-    newexpr = expand(newexpr)
+    integrals     = form.expr.args if isinstance(form.expr, Add) else [form.expr]
+    new_integrals = []
 
-    e = newexpr.series(eps, 0, 2)
-    d = collect(e, eps, evaluate=False)
-    expr = d[eps]
+    for I in integrals:
 
-#    print('> linearize = ', expr)
-#    import sys; sys.exit(0)
+        g0 = I.expr
+        g1 = expand(I.expr.subs(zip(fields, new_fields)))
+        dg_du = (g1 - g0).as_leading_term(eps).as_coefficient(eps)
 
-    test_trial = (trial_functions, test_functions)
+        if dg_du:
+            new_I = integral(I.domain, dg_du)
+            new_integrals.append(new_I)
 
-    if is_form:
-        return BilinearForm(test_trial, integral(domain, expr))
+    bilinear_expr = Add(*new_integrals)
+    tests = form.variables
 
-    else:
-        return BilinearExpr(test_trial, expr)
+    return BilinearForm((trials, tests), bilinear_expr)
 
 #==============================================================================
 def is_bilinear_form(expr):
@@ -895,7 +821,7 @@ def is_linear_expression(expr, args, debug=True):
 
     return True
 
-
+#==============================================================================
 def integral(domain, expr):
     """."""
     return Integral(expr, domain)
diff --git a/sympde/topology/domain.py b/sympde/topology/domain.py
index 70d8498..f30c31e 100644
--- a/sympde/topology/domain.py
+++ b/sympde/topology/domain.py
@@ -282,7 +282,7 @@ class Domain(BasicDomain):
         # ...
 
         h5  = h5py.File( filename, mode='r' )
-        yml = yaml.load( h5['topology.yml'][()] )
+        yml = yaml.load( h5['topology.yml'][()], Loader=yaml.SafeLoader )
 
         domain_name    = yml['name']
         dim            = yml['dim']
diff --git a/sympde/topology/mapping.py b/sympde/topology/mapping.py
index 8e4c9fa..4e6af3b 100644
--- a/sympde/topology/mapping.py
+++ b/sympde/topology/mapping.py
@@ -2,46 +2,32 @@
 
 from collections  import OrderedDict
 
-from sympy.core import Basic
 from sympy import Indexed, IndexedBase, Matrix, ImmutableDenseMatrix
-from sympy.core import Symbol
-from sympy.core.containers import Tuple
 from sympy import Function
-from sympy.core import Add, Mul, Pow
-from sympy.core.singleton import S
-from sympy.core.expr import AtomicExpr
-from sympy import Rational
-from sympy.core.numbers import ImaginaryUnit
 from sympy import sympify
-from sympy.core.function import AppliedUndef
+from sympy.core import Basic
+from sympy.core import Symbol
+from sympy.core import Add, Mul, Pow
+from sympy.core.expr       import AtomicExpr
+from sympy.core.numbers    import ImaginaryUnit
+from sympy.core.containers import Tuple
 
+from sympde.core       import Constant
 from sympde.core.basic import BasicMapping
-from sympde.core import Constant
-from sympde.core.algebra import (Dot_1d,
-                                 Dot_2d, Inner_2d, Cross_2d,
-                                 Dot_3d, Inner_3d, Cross_3d)
-
 from sympde.core.basic import CalculusFunction
 from sympde.core.basic import _coeffs_registery
-from .basic import BasicDomain
-from .domain import Domain
 
+from .basic       import BasicDomain
+from .space       import ScalarTestFunction, VectorTestFunction, IndexedTestTrial
+from .space       import ScalarField, VectorField, IndexedVectorField
+from .datatype    import HcurlSpaceType
 from .derivatives import dx, dy, dz
 from .derivatives import _partial_derivatives
 from .derivatives import get_atom_derivatives, get_index_derivatives_atom
-from .derivatives import (Grad_1d, Div_1d,
-                          Grad_2d, Curl_2d, Rot_2d, Div_2d,
-                          Grad_3d, Curl_3d, Div_3d)
-
-from .derivatives import dx1, dx2, dx3
 from .derivatives import _logical_partial_derivatives
 from .derivatives import get_atom_logical_derivatives, get_index_logical_derivatives_atom
 from .derivatives import LogicalGrad_1d, LogicalGrad_2d, LogicalGrad_3d
 
-from .space import ScalarTestFunction, VectorTestFunction, IndexedTestTrial
-from .space import ScalarField, VectorField, IndexedVectorField
-from .datatype import H1SpaceType, HcurlSpaceType, HdivSpaceType, L2SpaceType, UndefinedSpaceType
-
 #==============================================================================
 class Mapping(BasicMapping):
     """
diff --git a/sympde/topology/space.py b/sympde/topology/space.py
index bdef15b..42e4e5a 100644
--- a/sympde/topology/space.py
+++ b/sympde/topology/space.py
@@ -13,6 +13,7 @@ from sympy import Function
 from sympy.core.numbers import Zero as sy_Zero
 from sympy.core.singleton import S
 
+from sympde.core.utils import expand_name_patterns
 from sympde.core.utils import random_string
 from sympde.core.basic import CalculusFunction
 from sympde.core.basic import _coeffs_registery
@@ -20,18 +21,88 @@ from .basic import BasicDomain, Union, Interval
 from .datatype import SpaceType, dtype_space_registry
 from .datatype import RegularityType, dtype_regularity_registry
 
-
+#==============================================================================
 def element_of(space, name):
-    assert isinstance(space, BasicFunctionSpace)
-    if isinstance(name, (list,tuple,Tuple)):
-        return [space.element(nm) for nm in name]
-    return space.element(name)
+    """ Create a single element of a given space (possibly a ProductSpace).
+
+    Parameters
+    ----------
+    space : ScalarFunctionSpace | VectorFunctionSpace | ProductSpace
+        Function space from which a single element should be created.
+
+    name : str | iterable
+        If space is ProductSpace, 'name' must either be an explicit list of
+        function names, or a single pattern string that will be expanded into
+        such a list. Otherwise, 'name' must be a simple string.
+
+    Results
+    -------
+    res : ScalarTestFunction | VectorTestFunction | iterable
+        Single element taken from the given space. If space is ProductSpace,
+        an element is a list of functions; otherwise, it is a single function.
+
+    """
+    if not isinstance(space, BasicFunctionSpace):
+        raise TypeError(space)
+    names = expand_name_patterns(name)
+    return _recursive_element_of(space, names)
+
+#----------------------------------------
+def _recursive_element_of(space, names):
+
+    if isinstance(names, str):
+        name = names
+        return space.element(name)
+
+    elif isinstance(space, ProductSpace):
+        spaces = space.spaces
+        result = [_recursive_element_of(s, n) for s, n in zip(spaces, names)]
+        return type(names)(result)
+
+    else:
+        msg = "To create multiple elements of same space, use 'elements_of'."
+        raise ValueError(msg)
 
+#==============================================================================
 def elements_of(space, names):
-    assert(isinstance(names, str))
-    names = names.split(',')
-    names = [n.replace(' ', '') for n in names]
-    return [element_of(space, n) for n in names]
+    """ Create multiple elements of same space (possibly a ProductSpace).
+
+    Parameters
+    ----------
+    space : ScalarFunctionSpace | VectorFunctionSpace | ProductSpace
+
+    names : str | iterable
+        Pattern or list of patterns from which a list of function names is
+        produced.
+
+    Results
+    -------
+    res : iterable
+        Multiple elements taken from the given space. If space is ProductSpace,
+        each element is a list of functions; otherwise, each element is a
+        single function.
+
+    """
+    if not isinstance(space, BasicFunctionSpace):
+        raise TypeError(space)
+    names = expand_name_patterns(names, seq=True)
+    return _recursive_elements_of(space, names)
+
+#----------------------------------------
+def _recursive_elements_of(space, names):
+
+    if isinstance(names, str):
+        name = names
+        return space.element(name)
+
+    elif isinstance(space, ProductSpace):
+        spaces = space.spaces
+        result = [_recursive_elements_of(s, n) for s, n in zip(spaces, names)]
+        return type(names)(result)
+
+    else:
+        result = [_recursive_elements_of(space, n) for n in names]
+        return type(names)(result)
 
 #==============================================================================
 class BasicFunctionSpace(Basic):

use element_of instead of ScalarTestFunction and VectorTestFunction in our unit tests
diff --git a/sympde/calculus/tests/test_calculus.py b/sympde/calculus/tests/test_calculus.py
index 415a9ab..4fe3ae1 100644
--- a/sympde/calculus/tests/test_calculus.py
+++ b/sympde/calculus/tests/test_calculus.py
@@ -2,20 +2,11 @@
 
 import pytest
 
-from sympy import Symbol
-from sympy.core.containers import Tuple
-from sympy import symbols
-from sympy import IndexedBase
-from sympy import Matrix
-from sympy import Function
-from sympy import pi, cos, sin
-from sympy import srepr
-from sympy import expand
 from sympy import Function
 from sympy import Integer, Float, Rational
-from sympy.physics.quantum import TensorProduct
+from sympy import expand
 
-from sympde.core import Constant
+from sympde.core     import Constant
 from sympde.calculus import grad, dot, inner, outer, cross, rot, curl, div
 from sympde.calculus import laplace, hessian, bracket, convect, D, conv
 from sympde.calculus import ArgumentTypeError
@@ -23,12 +14,8 @@ from sympde.calculus import jump, avg, Dn, minus, plus
 from sympde.topology import Domain
 from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
 from sympde.topology import ProductSpace
-from sympde.topology import H1Space, HcurlSpace, HdivSpace, L2Space, UndefinedSpace
-from sympde.topology import TestFunction, ScalarTestFunction, VectorTestFunction
-from sympde.topology import ScalarField, VectorField
 from sympde.topology import element_of, elements_of
 
-
 #==============================================================================
 def test_zero_derivative():
 
@@ -59,7 +46,7 @@ def test_zero_derivative():
     # 2D convection of constant scalar field
     domain = Domain('Omega', dim=2)
     W = VectorFunctionSpace('W', domain)
-    F = VectorField(W, name='F')
+    F = element_of(W, name='F')
 
     assert convect(F, 1)              == 0  # native int
     assert convect(F, 2.3)            == 0  # native float
@@ -72,7 +59,7 @@ def test_zero_derivative():
     # 3D convection of constant scalar field
     domain = Domain('Omega', dim=3)
     Z = VectorFunctionSpace('Z', domain)
-    G = VectorField(Z, name='G')
+    G = element_of(Z, name='G')
 
     assert convect(G, 1)              == 0  # native int
     assert convect(G, 2.3)            == 0  # native float
@@ -91,12 +78,12 @@ def test_calculus_2d_1():
 
     alpha, beta, gamma = [Constant(i) for i in ['alpha','beta','gamma']]
 
-    f,g,h = [ScalarField(V, name=i) for i in ['f','g','h']]
-    F,G,H = [VectorField(W, i) for i in ['F','G','H']]
+    f, g, h = elements_of(V, names='f, g, h')
+    F, G, H = elements_of(W, names='F, G, H')
 
     # ... scalar gradient properties
     assert( grad(f+g) == grad(f) + grad(g) )
-    assert( grad(alpha*f) == alpha*grad(f) )
+    assert( grad(alpha*h) == alpha*grad(h) )
     assert( grad(alpha*f + beta*g) == alpha*grad(f) + beta*grad(g)  )
 
     assert( grad(f*g) == f*grad(g) + g*grad(f) )
@@ -107,7 +94,7 @@ def test_calculus_2d_1():
 
     # ... vector gradient properties
     assert( grad(F+G) == grad(F) + grad(G) )
-    assert( grad(alpha*F) == alpha*grad(F) )
+    assert( grad(alpha*H) == alpha*grad(H) )
     assert( grad(alpha*F + beta*G) == alpha*grad(F) + beta*grad(G)  )
 
     assert( grad(dot(F,G)) == convect(F, G) + convect(G, F) + cross(F, curl(G)) - cross(curl(F), G) )
@@ -115,19 +102,19 @@ def test_calculus_2d_1():
 
     # ... curl properties
     assert( curl(f+g) == curl(f) + curl(g) )
-    assert( curl(alpha*f) == alpha*curl(f) )
+    assert( curl(alpha*h) == alpha*curl(h) )
     assert( curl(alpha*f + beta*g) == alpha*curl(f) + beta*curl(g)  )
     # ...
 
     # ... laplace properties
     assert( laplace(f+g) == laplace(f) + laplace(g) )
-    assert( laplace(alpha*f) == alpha*laplace(f) )
+    assert( laplace(alpha*h) == alpha*laplace(h) )
     assert( laplace(alpha*f + beta*g) == alpha*laplace(f) + beta*laplace(g)  )
     # ...
 
     # ... divergence properties
     assert( div(F+G) == div(F) + div(G) )
-    assert( div(alpha*F) == alpha*div(F) )
+    assert( div(alpha*H) == alpha*div(H) )
     assert( div(alpha*F + beta*G) == alpha*div(F) + beta*div(G)  )
 
     assert( div(cross(F,G)) == -dot(F, curl(G)) + dot(G, curl(F)) )
@@ -135,24 +122,19 @@ def test_calculus_2d_1():
 
     # ... rot properties
     assert( rot(F+G) == rot(F) + rot(G) )
-    assert( rot(alpha*F) == alpha*rot(F) )
+    assert( rot(alpha*H) == alpha*rot(H) )
     assert( rot(alpha*F + beta*G) == alpha*rot(F) + beta*rot(G)  )
     # ...
 
-    expr = grad(dot(F,G))
-    print(expr)
-
 #==============================================================================
 def test_calculus_2d_2():
     domain = Domain('Omega', dim=2)
 
-    V = ScalarFunctionSpace('V', domain)
     W = VectorFunctionSpace('W', domain)
 
     alpha, beta, gamma = [Constant(i) for i in ['alpha','beta','gamma']]
 
-    f,g,h = [ScalarField(V, name=i) for i in ['f','g','h']]
-    F,G,H = [VectorField(W, i) for i in ['F','G','H']]
+    F, G, H = elements_of(W, names='F, G, H')
 
     # ... vector gradient properties
     assert( convect(F+G, H) == convect(F,H) + convect(G,H) )
@@ -169,7 +151,7 @@ def test_calculus_2d_3():
 
     alpha, beta, gamma = [Constant(i) for i in ['alpha','beta','gamma']]
 
-    f,g,h = [ScalarField(V, name=i) for i in ['f','g','h']]
+    f, g = elements_of(V, names='f, g')
     K = Function('K')(x,y)
     S = Function('S')(x,y)
 
@@ -187,12 +169,12 @@ def test_calculus_3d():
 
     alpha, beta, gamma = [Constant(i) for i in ['alpha','beta','gamma']]
 
-    f,g,h = [ScalarField(V, name=i) for i in ['f','g','h']]
-    F,G,H = [VectorField(W, i) for i in ['F','G','H']]
+    f, g, h = elements_of(V, names='f, g, h')
+    F, G, H = elements_of(W, names='F, G, H')
 
     # ... gradient properties
     assert( grad(f+g) == grad(f) + grad(g) )
-    assert( grad(alpha*f) == alpha*grad(f) )
+    assert( grad(alpha*h) == alpha*grad(h) )
     assert( grad(alpha*f + beta*g) == alpha*grad(f) + beta*grad(g)  )
 
     assert( grad(f*g) == f*grad(g) + g*grad(f) )
@@ -203,7 +185,7 @@ def test_calculus_3d():
 
     # ... curl properties
     assert( curl(F+G) == curl(F) + curl(G) )
-    assert( curl(alpha*F) == alpha*curl(F) )
+    assert( curl(alpha*H) == alpha*curl(H) )
     assert( curl(alpha*F + beta*G) == alpha*curl(F) + beta*curl(G)  )
 
     assert( curl(cross(F,G)) == F*div(G) - G*div(F) - convect(F, G) + convect(G, F) )
@@ -212,14 +194,14 @@ def test_calculus_3d():
 
     # ... laplace properties
     assert( laplace(f+g) == laplace(f) + laplace(g) )
-    assert( laplace(alpha*f) == alpha*laplace(f) )
+    assert( laplace(alpha*h) == alpha*laplace(h) )
     assert( laplace(alpha*f + beta*g) == alpha*laplace(f) + beta*laplace(g)  )
     assert( laplace(f*g) == f*laplace(g) + g*laplace(f) + 2*dot(grad(f), grad(g)) )
     # ...
 
     # ... divergence properties
     assert( div(F+G) == div(F) + div(G) )
-    assert( div(alpha*F) == alpha*div(F) )
+    assert( div(alpha*H) == alpha*div(H) )
     assert( div(alpha*F + beta*G) == alpha*div(F) + beta*div(G)  )
 
     assert( div(cross(F,G)) == -dot(F, curl(G)) + dot(G, curl(F)) )
@@ -227,8 +209,8 @@ def test_calculus_3d():
     # ...
 
     # ...
-    assert( curl(grad(f)) == 0 )
-    assert( div(curl(F)) == 0 )
+    assert( curl(grad(h)) == 0 )
+    assert( div(curl(H)) == 0 )
     assert( div(cross(grad(f), grad(g))) == 0 )
     assert( curl(curl(F)) == grad(div(F)) - laplace(F))
     assert( curl(f*grad(g)) == cross(grad(f), grad(g)) )
@@ -246,14 +228,11 @@ def test_calculus_3d_3():
     V     = ScalarFunctionSpace('V', domain, kind=None)
     W     = VectorFunctionSpace('W', domain, kind=None)
 
-    f,g,h = [ScalarField(V, name=i) for i in ['f','g','h']]
-    F,G,H = [VectorField(W, i) for i in ['F','G','H']]
-
     X = ProductSpace(H1, Hcurl, Hdiv, L2)
-    v0, v1, v2, v3 = TestFunction(X, ['v0', 'v1', 'v2', 'v3'])
+    v0, v1, v2, v3 = element_of(X, ['v0', 'v1', 'v2', 'v3'])
 
-    v = TestFunction(V, 'v')
-    w = TestFunction(W, 'w')
+    v = element_of(V, 'v')
+    w = element_of(W, 'w')
 
     # ... consistency of grad operator
     # we can apply grad on an undefined space type or H1
@@ -307,12 +286,7 @@ def test_calculus_3d_4():
 
     alpha, beta, gamma = [Constant(i) for i in ['alpha','beta','gamma']]
 
-    F,G,H = [VectorField(W, i) for i in ['F','G','H']]
-
-    # ...
-    expr = inner(D(F), D(G))
-    print(expr)
-    # ...
+    F, G, H = elements_of(W, names='F, G, H')
 
     # ...
     expected = alpha*inner(D(F), D(G)) + beta*inner(D(F), D(H))
@@ -327,11 +301,7 @@ def test_calculus_3d_5():
 
     alpha, beta, gamma = [Constant(i) for i in ['alpha','beta','gamma']]
 
-    E,F,G,H = [VectorField(W, i) for i in ['E','F','G','H']]
-
-    # ...
-    expr = outer(F, G)
-    # ...
+    F, G, H = elements_of(W, names='F, G, H')
 
     # ...
     expected = alpha*outer(F, G) + beta*outer(F, H)
diff --git a/sympde/core/tests/test_utils.py b/sympde/core/tests/test_utils.py
new file mode 100644
index 0000000..53e581e
--- /dev/null
+++ b/sympde/core/tests/test_utils.py
@@ -0,0 +1,33 @@
+from sympde.core.utils import expand_name_patterns
+
+#==============================================================================
+def test_expand_name_patterns():
+
+    assert expand_name_patterns('x,y,z') == ('x', 'y', 'z')
+    assert expand_name_patterns('x y z') == ('x', 'y', 'z')
+
+    assert expand_name_patterns('x,') == ('x',)
+    assert expand_name_patterns('x', seq=True) == ('x',)
+
+    assert expand_name_patterns(('a', 'b', 'c')) == ('a', 'b', 'c')
+    assert expand_name_patterns(['a', 'b', 'c']) == ['a', 'b', 'c']
+    assert expand_name_patterns({'a', 'b', 'c'}) == {'a', 'b', 'c'}
+
+    assert expand_name_patterns('x:4') == ('x0', 'x1', 'x2', 'x3')
+    assert expand_name_patterns('x5:10') == ('x5', 'x6', 'x7', 'x8', 'x9')
+    assert expand_name_patterns('x7(1:3)') == ('x71', 'x72')
+    assert expand_name_patterns('x2:5, y:2') == ('x2', 'x3', 'x4', 'y0', 'y1')
+    assert expand_name_patterns(('x2:5', 'y:2')) == (('x2', 'x3', 'x4'), ('y0', 'y1'))
+
+    assert expand_name_patterns(':c') == ('a', 'b', 'c')
+    assert expand_name_patterns('x:c') == ()
+    assert expand_name_patterns('x:z') == ('x', 'y', 'z')
+    assert expand_name_patterns('x(:c)') == ('xa', 'xb', 'xc')
+    assert expand_name_patterns('a:b, x:z') == ('a', 'b', 'x', 'y', 'z')
+    assert expand_name_patterns(('a:b', 'x:z')) == (('a', 'b'), ('x', 'y', 'z'))
+
+    assert expand_name_patterns('x:2(1:3)') == ('x01', 'x02', 'x11', 'x12')
+    assert expand_name_patterns(':3:2') == ('00', '01', '10', '11', '20', '21')
+
+    assert expand_name_patterns('x((a:b))') == ('x(a)', 'x(b)')
+    assert expand_name_patterns(r'x(:1\,:2)') == ('x(0,0)', 'x(0,1)')
diff --git a/sympde/expr/tests/test_equation_2d.py b/sympde/expr/tests/test_equation_2d.py
index 962d213..4fe7abd 100644
--- a/sympde/expr/tests/test_equation_2d.py
+++ b/sympde/expr/tests/test_equation_2d.py
@@ -1,43 +1,26 @@
 # coding: utf-8
 
-import pytest
+#import pytest
 
-from sympy import Symbol
 from sympy.core.containers import Tuple
-from sympy import symbols
-from sympy import IndexedBase
-from sympy import Matrix
-from sympy import Function
-from sympy import pi, cos, sin, exp
-from sympy import srepr
-from sympy.physics.quantum import TensorProduct
-
-from sympde.core import Constant
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
-from sympde.calculus import laplace, hessian, bracket
-from sympde.topology import (dx, dy, dz)
-from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
-from sympde.topology import ScalarField, VectorField
-from sympde.topology import element_of
-from sympde.topology import ProductSpace
-from sympde.topology import ScalarTestFunction
-from sympde.topology import VectorTestFunction
-from sympde.topology import Unknown
-from sympde.topology import Domain, Boundary, NormalVector, TangentVector
-from sympde.topology import Trace, trace_0, trace_1
-from sympde.topology import Square
-from sympde.topology import ElementDomain
-from sympde.topology import Area
-
-from sympde.expr import BilinearForm, LinearForm, integral
-from sympde.expr import TerminalExpr
-
+from sympy import pi, cos, sin
+
+from sympde.core        import Constant
+from sympde.calculus    import grad, dot, inner
+from sympde.calculus    import laplace, bracket
+from sympde.topology    import ScalarFunctionSpace, VectorFunctionSpace
+from sympde.topology    import element_of
+from sympde.topology    import ProductSpace
+from sympde.topology    import Domain, Boundary, NormalVector
+from sympde.topology    import Square
+from sympde.topology    import ElementDomain
+from sympde.topology    import Area
+from sympde.expr        import BilinearForm, LinearForm, integral
+from sympde.expr        import Equation, EssentialBC
 from sympde.expr.errors import UnconsistentLhsError
 from sympde.expr.errors import UnconsistentRhsError
 from sympde.expr.errors import UnconsistentBCError
 
-from sympde.expr import Equation, EssentialBC, NewtonIteration
-
 DIM = 2
 domain = Domain('Omega', dim=DIM)
 
diff --git a/sympde/expr/tests/test_essential_bc.py b/sympde/expr/tests/test_essential_bc.py
index 0eed27f..b83ea9e 100644
--- a/sympde/expr/tests/test_essential_bc.py
+++ b/sympde/expr/tests/test_essential_bc.py
@@ -1,33 +1,11 @@
 # coding: utf-8
 
-import pytest
-
-from sympy import Symbol
-from sympy.core.containers import Tuple
-from sympy import symbols
-from sympy import IndexedBase
-from sympy import Matrix
-from sympy import Function
-from sympy import pi, cos, sin
-from sympy import srepr
-from sympy.physics.quantum import TensorProduct
-
-from sympde.core import Constant
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
-from sympde.calculus import laplace, hessian, bracket
-from sympde.topology import (dx, dy, dz)
+from sympde.core     import Constant
+from sympde.calculus import grad, dot
 from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
-from sympde.topology import ScalarField, VectorField
-from sympde.topology import ProductSpace
-from sympde.topology import ScalarTestFunction
-from sympde.topology import VectorTestFunction
-from sympde.topology import Unknown
-from sympde.topology import Domain, Boundary, NormalVector, TangentVector
-from sympde.topology import Trace, trace_0, trace_1
-from sympde.topology import Square
-
-from sympde.expr import Equation, EssentialBC
-
+from sympde.topology import element_of
+from sympde.topology import Domain, Boundary, NormalVector
+from sympde.expr     import EssentialBC
 
 #==============================================================================
 def test_essential_bc_1():
@@ -36,8 +14,8 @@ def test_essential_bc_1():
     V = ScalarFunctionSpace('V', domain)
     W = VectorFunctionSpace('W', domain)
 
-    v = ScalarTestFunction(V, name='v')
-    w = VectorTestFunction(W, name='w')
+    v = element_of(V, name='v')
+    w = element_of(W, name='w')
 
     B1 = Boundary(r'\Gamma_1', domain)
     nn = NormalVector('nn')
diff --git a/sympde/expr/tests/test_expr_2d.py b/sympde/expr/tests/test_expr_2d.py
index 479f739..5e8a2ae 100644
--- a/sympde/expr/tests/test_expr_2d.py
+++ b/sympde/expr/tests/test_expr_2d.py
@@ -4,30 +4,22 @@
 
 import pytest
 
-from sympy import Symbol
 from sympy.core.containers import Tuple
-from sympy import symbols
-from sympy import IndexedBase
 from sympy import Matrix
 from sympy import Function
 from sympy import pi, cos, sin, exp
-from sympy import srepr
-from sympy.physics.quantum import TensorProduct
 
 from sympde.core import Constant
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
-from sympde.calculus import laplace, hessian, bracket, convect
+from sympde.calculus import grad, dot, inner, rot, div
+from sympde.calculus import laplace, bracket, convect
 from sympde.calculus import jump, avg, Dn, minus, plus
-from sympde.topology import (dx, dy, dz)
+#from sympde.topology import dx
 from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
-from sympde.topology import ProductSpace
 from sympde.topology import element_of, elements_of
-from sympde.topology import Unknown
 from sympde.topology import InteriorDomain, Union
-from sympde.topology import Boundary, NormalVector, TangentVector
+from sympde.topology import Boundary, NormalVector
 from sympde.topology import Domain
-from sympde.topology import Trace, trace_0, trace_1
-from sympde.topology import Mapping
+from sympde.topology import trace_1
 from sympde.topology import Square
 from sympde.topology import ElementDomain
 from sympde.topology import Area
@@ -39,16 +31,6 @@ from sympde.expr.expr import Functional, Norm
 from sympde.expr.expr import linearize
 from sympde.expr.evaluation import TerminalExpr
 
-from sympde.expr.errors import UnconsistentError
-from sympde.expr.errors import UnconsistentLinearExpressionError
-from sympde.expr.errors import UnconsistentLhsError
-from sympde.expr.errors import UnconsistentRhsError
-from sympde.expr.errors import UnconsistentBCError
-
-DIM = 2
-VERBOSE = False
-VERBOSE = True
-
 #==============================================================================
 def test_linear_expr_2d_1():
 
@@ -938,8 +920,9 @@ def test_terminal_expr_bilinear_2d_4():
     print('')
 
 #==============================================================================
+@pytest.mark.skip(reason="New linearize() function does not accept 'LinearExpr' objects")
 def test_linearize_expr_2d_1():
-    domain = Domain('Omega', dim=DIM)
+    domain = Domain('Omega', dim=2)
     x,y = domain.coordinates
 
     V1 = ScalarFunctionSpace('V1', domain)
@@ -997,8 +980,9 @@ 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=DIM)
+    domain = Domain('Omega', dim=2)
     x,y = domain.coordinates
 
     V1 = ScalarFunctionSpace('V1', domain)
@@ -1023,140 +1007,143 @@ def test_linearize_expr_2d_2():
 
 #==============================================================================
 def test_linearize_form_2d_1():
-    domain = Domain('Omega', dim=DIM)
-    x,y = domain.coordinates
-
-    V1 = ScalarFunctionSpace('V1', domain)
-    W1 = VectorFunctionSpace('W1', domain)
-
-    v1 = element_of(V1, name='v1')
-    w1 = element_of(W1, name='w1')
+    domain = Domain('Omega', dim=2)
 
-    alpha = Constant('alpha')
+    V = ScalarFunctionSpace('V', domain)
+    W = VectorFunctionSpace('W', domain)
 
-    F = element_of(V1, name='F')
-    G = element_of(W1, 'G')
+    v, F, u = elements_of(V, names='v, F, u')
+    w, G, m = elements_of(W, names='w, G, m')
 
     int_0 = lambda expr: integral(domain , expr)
 
     # ...
-    l = LinearForm(v1, int_0(F**2*v1))
-    a = linearize(l, F, trials='u1')
-    print(a)
+    l  = LinearForm(v, int_0(F**2*v))
+    a  = linearize(l, F, trials=u)
+    assert a(u, v) == int_0(2 * F * u * v)
     # ...
 
     # ...
-    l = LinearForm(v1, int_0(dot(grad(F), grad(F))*v1))
-    a = linearize(l, F, trials='u1')
-    print(a)
+    l = LinearForm(v, int_0(dot(grad(F), grad(F))*v))
+    a = linearize(l, F, trials=u)
+    assert a(u, v) == int_0(2 * dot(grad(F), grad(u)) * v)
     # ...
 
     # ...
-    l = LinearForm(v1, int_0(exp(-F)*v1))
-    a = linearize(l, F, trials='u1')
-    print(a)
+    l = LinearForm(v, int_0(exp(-F)*v))
+    a = linearize(l, F, trials=u)
+    assert a(u, v) == int_0(-exp(-F) * u * v)
     # ...
 
     # ...
-    l = LinearForm(v1, int_0(cos(F)*v1))
-    a = linearize(l, F, trials='u1')
-    print(a)
+    l = LinearForm(v, int_0(cos(F)*v))
+    a = linearize(l, F, trials=u)
+    assert a(u, v) == int_0(-sin(F) * u * v)
     # ...
 
     # ...
-    l = LinearForm(v1, int_0(cos(F**2)*v1))
-    a = linearize(l, F, trials='u1')
-    print(a)
+    l = LinearForm(v, int_0(cos(F**2)*v))
+    a = linearize(l, F, trials=u)
+    assert a(u, v) == int_0(-2 * F *sin(F**2) * u * v)
     # ...
 
     # ...
-    l = LinearForm(v1, int_0(F**2*dot(grad(F), grad(v1))))
-    a = linearize(l, F, trials='u1')
-    print(a)
+    l = LinearForm(v, int_0(F**2*dot(grad(F), grad(v))))
+    a = linearize(l, F, trials=u)
+    assert a(u, v) == int_0(2 * F * u * dot(grad(F), grad(v)) + F**2 * dot(grad(u), grad(v)))
     # ...
 
     # ...
-    l = LinearForm(w1, int_0(dot(rot(G), grad(G))*w1))
-    a = linearize(l, G, trials='u1')
-    print(a)
+    l = LinearForm(w, int_0(dot(rot(G), grad(G))*w))
+    a = linearize(l, G, trials=m)
+    assert a(m, w) == int_0((dot(rot(m), grad(G)) + dot(rot(G), grad(m))) * w)
     # ...
 
 #==============================================================================
 def test_linearize_form_2d_2():
-    domain = Domain('Omega', dim=DIM)
-    x,y = domain.coordinates
-
-    V1 = ScalarFunctionSpace('V1', domain)
-
-    v1 = element_of(V1, name='v1')
+    domain = Domain('Omega', dim=2)
 
-    alpha = Constant('alpha')
+    V = ScalarFunctionSpace('V', domain)
 
-    F = element_of(V1, name='F')
-    G = element_of(V1, name='G')
+    v, F, u = elements_of(V, names='v, F, u')
 
     int_0 = lambda expr: integral(domain , expr)
 
     # ...
-    l1 = LinearForm(v1, int_0(F**2*v1))
-    l  = LinearForm(v1, l1(v1))
+    l1 = LinearForm(v, int_0(F**2*v))
+    l  = LinearForm(v, l1(v))
 
-    a = linearize(l, F, trials='u1')
-    print(a)
+    a = linearize(l, F, trials=u)
 
-    expected = linearize(l1, F, trials='u1')
-    assert( linearize(l, F, trials='u1') == expected )
+    expected = linearize(l1, F, trials=u)
+    assert a == expected
     # ...
 
 #==============================================================================
 def test_linearize_form_2d_3():
     """steady Euler equation."""
-    domain = Domain('Omega', dim=DIM)
-    x,y = domain.coordinates
+    domain = Domain('Omega', dim=2)
 
     U = VectorFunctionSpace('U', domain)
     W = ScalarFunctionSpace('W', domain)
 
+    # Test functions
     v   = element_of(U, name='v')
-    phi =       element_of(W, name='phi')
-    q   =       element_of(W, name='q')
+    phi = element_of(W, name='phi')
+    q   = element_of(W, name='q')
 
+    # Steady-state fields
     U_0   = element_of(U, name='U_0')
     Rho_0 = element_of(W, name='Rho_0')
     P_0   = element_of(W, name='P_0')
 
+    # Trial functions (displacements from steady-state)
+    d_u   = element_of(U, name='d_u')
+    d_rho = element_of(W, name='d_rho')
+    d_p   = element_of(W, name='d_p')
+
+    # Shortcut
     int_0 = lambda expr: integral(domain , expr)
 
-    # ...
-    expr = div(Rho_0*U_0) * phi
-    l1 = LinearForm(phi, int_0(expr))
+    # The Euler equations are a system of three non-linear equations; for each of
+    # them we create a linear form in the test functions (phi, v, q) respectively.
+    e1 = div(Rho_0 * U_0)
+    l1 = LinearForm(phi, int_0(e1 * phi))
 
-    expr = Rho_0*dot(convect(U_0, grad(U_0)), v) + dot(grad(P_0), v)
-    l2 = LinearForm(v, int_0(expr))
+    e2 = Rho_0 * convect(U_0, U_0) + grad(P_0)
+    l2 = LinearForm(v, int_0(dot(e2, v)))
 
-    expr = dot(U_0, grad(P_0)) * q + P_0 * div(U_0) * q
-    l3 = LinearForm(q, int_0(expr))
+    e3 = div(P_0 * U_0)
+    l3 = LinearForm(q, int_0(e3 * q))
     # ...
 
-    a1 = linearize(l1, [Rho_0, U_0], trials=['d_rho', 'd_u'])
-    print(a1)
-    print('')
+    # Linearize l1, l2 and l3 separately
+    a1 = linearize(l1, fields=[Rho_0, U_0     ], trials=[d_rho, d_u     ])
+    a2 = linearize(l2, fields=[Rho_0, U_0, P_0], trials=[d_rho, d_u, d_p])
+    a3 = linearize(l3, fields=[       U_0, P_0], trials=[       d_u, d_p])
 
-    a2 = linearize(l2, [Rho_0, U_0, P_0], trials=['d_rho', 'd_u', 'd_p'])
-    print(a2)
-    print('')
+    # Check individual bilinear forms
+    d_e1 = div(U_0 * d_rho + Rho_0 * d_u)
+    d_e2 = d_rho * convect(U_0, U_0) + \
+           Rho_0 * convect(d_u, U_0) + \
+           Rho_0 * convect(U_0, d_u) + grad(d_p)
+    d_e3 = div(d_p * U_0 + P_0 * d_u)
 
-    a3 = linearize(l3, [P_0, U_0], trials=['d_p', 'd_u'])
-    print(a3)
-    print('')
+    assert a1([d_rho, d_u     ], phi) == int_0(d_e1 * phi)
+    assert a2([d_rho, d_u, d_p], v  ) == int_0(dot(d_e2, v))
+    assert a3([       d_u, d_p], q  ) == int_0(d_e3 * q)
 
+    # Linearize linear form of system: l = l1 + l2 + l3
     l = LinearForm((phi, v, q), l1(phi) + l2(v) + l3(q))
-    a = linearize(l, [Rho_0, U_0, P_0], trials=['d_rho', 'd_u', 'd_p'])
-    print(a)
+    a = linearize(l, fields=[Rho_0, U_0, P_0], trials=[d_rho, d_u, d_p])
+
+    # Check composite linear form
+    assert a([d_rho, d_u, d_p], [phi, v, q]) == \
+            int_0(d_e1 * phi + dot(d_e2, v) + d_e3 * q)
 
 #==============================================================================
 def test_linearize_form_2d_4():
-    domain = Domain('Omega', dim=DIM)
+    domain  = Domain('Omega', dim=2)
     Gamma_N = Boundary(r'\Gamma_N', domain)
 
     x,y = domain.coordinates
@@ -1164,10 +1151,9 @@ def test_linearize_form_2d_4():
     V = ScalarFunctionSpace('V', domain)
 
     v  = element_of(V, name='v')
+    u  = element_of(V, name='u')
     du = element_of(V, name='du')
 
-    u = element_of(V, name='u')
-
     int_0 = lambda expr: integral(domain , expr)
     int_1 = lambda expr: integral(Gamma_N, expr)
 
@@ -1180,9 +1166,8 @@ def test_linearize_form_2d_4():
 
     # linearising l around u, using du
     a = linearize(l, u, trials=du)
-    # ...
-    print(a)
-    # ...
+
+    assert a(du, v) == int_0(dot(grad(v), grad(du)) + 4.*exp(-u) * du * v)
 
 #==============================================================================
 def test_area_2d_1():
diff --git a/sympde/expr/tests/test_find_2d.py b/sympde/expr/tests/test_find_2d.py
index 5f88a44..816e954 100644
--- a/sympde/expr/tests/test_find_2d.py
+++ b/sympde/expr/tests/test_find_2d.py
@@ -1,42 +1,19 @@
 # coding: utf-8
 
-import pytest
-
-from sympy import Symbol
-from sympy.core.containers import Tuple
-from sympy import symbols
-from sympy import IndexedBase
-from sympy import Matrix
-from sympy import Function
-from sympy import pi, cos, sin, exp
-from sympy import srepr
-from sympy.physics.quantum import TensorProduct
-
-from sympde.core import Constant
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
-from sympde.calculus import laplace, hessian, bracket
-from sympde.topology import (dx, dy, dz)
-from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
-from sympde.topology import ScalarField, VectorField
-from sympde.topology import ProductSpace
-from sympde.topology import ScalarTestFunction
-from sympde.topology import VectorTestFunction
-from sympde.topology import Unknown
-from sympde.topology import Domain, Boundary, NormalVector, TangentVector
-from sympde.topology import Trace, trace_0, trace_1
-from sympde.topology import Square
-from sympde.topology import ElementDomain
-from sympde.topology import Area
-
-from sympde.expr import BilinearForm, LinearForm,integral
-from sympde.expr import TerminalExpr
-from sympde.expr import find
-
-from sympde.expr.errors import UnconsistentLhsError
-from sympde.expr.errors import UnconsistentRhsError
-from sympde.expr.errors import UnconsistentBCError
-
-from sympde.expr import Equation, EssentialBC, NewtonIteration
+from sympy import cos
+#from sympy import exp
+
+from sympde.core     import Constant
+from sympde.calculus import grad, dot
+#from sympde.topology import dx, dy
+from sympde.topology import ScalarFunctionSpace
+from sympde.topology import element_of
+from sympde.topology import Domain, Boundary
+from sympde.expr     import BilinearForm, LinearForm, integral
+from sympde.expr     import TerminalExpr
+from sympde.expr     import find
+from sympde.expr     import EssentialBC
+#from sympde.expr     import NewtonIteration
 
 DIM = 2
 domain = Domain('Omega', dim=DIM)
@@ -47,8 +24,8 @@ def test_find_2d_1():
     V = ScalarFunctionSpace('V', domain)
     U = ScalarFunctionSpace('U', domain)
 
-    v = ScalarTestFunction(V, name='v')
-    u = ScalarTestFunction(U, name='u')
+    v = element_of(V, name='v')
+    u = element_of(U, name='u')
 
     x,y = domain.coordinates
 
@@ -124,9 +101,9 @@ def test_find_2d_1():
 #
 #    x,y = domain.coordinates
 #
-#    Un = ScalarField(V, name='Un')
+#    Un = element_of(V, name='Un')
 #
-#    v = ScalarTestFunction(V, name='v')
+#    v = element_of(V, name='v')
 #
 #    f  = -4.*exp(-Un)
 #    l = LinearForm(v, dot(grad(v), grad(Un)) - f*v )
@@ -134,7 +111,7 @@ def test_find_2d_1():
 #    eq = NewtonIteration(l, Un, trials='u')
 #    # ...
 #
-#    u = ScalarTestFunction(V, name='u')
+#    u = element_of(V, name='u')
 #
 #    # ...
 #    expected =  -4.0*u*v*exp(-Un) + dx(u)*dx(v) + dy(u)*dy(v)
diff --git a/sympde/expr/tests/test_newton_2d.py b/sympde/expr/tests/test_newton_2d.py
index 985f3fd..6a78f87 100644
--- a/sympde/expr/tests/test_newton_2d.py
+++ b/sympde/expr/tests/test_newton_2d.py
@@ -1,41 +1,13 @@
 # coding: utf-8
 
-import pytest
-
-from sympy import Symbol
-from sympy.core.containers import Tuple
-from sympy import symbols
-from sympy import IndexedBase
-from sympy import Matrix
-from sympy import Function
-from sympy import pi, cos, sin, exp
-from sympy import srepr
-from sympy.physics.quantum import TensorProduct
-
-from sympde.core import Constant
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
-from sympde.calculus import laplace, hessian, bracket
-from sympde.topology import (dx, dy, dz)
-from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
-from sympde.topology import ScalarField, VectorField
-from sympde.topology import ProductSpace
-from sympde.topology import ScalarTestFunction
-from sympde.topology import VectorTestFunction
-from sympde.topology import Unknown
-from sympde.topology import Domain, Boundary, NormalVector, TangentVector
-from sympde.topology import Trace, trace_0, trace_1
-from sympde.topology import Square
-from sympde.topology import ElementDomain
-from sympde.topology import Area
-
-from sympde.expr import BilinearForm, LinearForm,integral
-from sympde.expr import TerminalExpr
-
-from sympde.expr.errors import UnconsistentLhsError
-from sympde.expr.errors import UnconsistentRhsError
-from sympde.expr.errors import UnconsistentBCError
-
-from sympde.expr import Equation, EssentialBC, NewtonIteration
+from sympy import exp
+
+from sympde.calculus import grad, dot
+from sympde.topology import ScalarFunctionSpace
+from sympde.topology import element_of
+from sympde.topology import Domain, Boundary
+from sympde.expr     import LinearForm, integral
+from sympde.expr     import EssentialBC, NewtonIteration
 
 DIM = 2
 domain = Domain('Omega', dim=DIM)
@@ -50,32 +22,22 @@ def test_newton_2d_1():
 
     x,y = domain.coordinates
 
-    Un = ScalarField(V, name='Un')
-
-    v = ScalarTestFunction(V, name='v')
+    Un = element_of(V, name='Un')
+    v  = element_of(V, name='v')
     
     int_0 = lambda expr: integral(domain , expr)
 
     f  = -4.*exp(-Un)
     l = LinearForm(v, int_0(dot(grad(v), grad(Un)) - f*v) )
 
-    eq = NewtonIteration(l, Un, trials='u')
-    # ...
-
-    u = ScalarTestFunction(V, name='u')
-
-    # ...
-    expected =  -4.0*u*v*exp(-Un) + dx(u)*dx(v) + dy(u)*dy(v)
-
-    expr = TerminalExpr(eq.lhs)[0]
-    assert(expr.expr == expected)
+    u  = element_of(V, name='u')
+    eq = NewtonIteration(l, Un, trials=u)
     # ...
 
     # ...
-    expected = -4.0*v*exp(-Un) - dx(Un)*dx(v) - dy(Un)*dy(v)
+    expected = int_0(-4.0*u*v*exp(-Un) + dot(grad(u), grad(v)))
 
-    expr = TerminalExpr(eq.rhs)[0]
-    assert(expr.expr == expected)
+    assert(eq.lhs.expr == expected)
     # ...
 
     # ...
@@ -83,7 +45,6 @@ def test_newton_2d_1():
     eq = NewtonIteration(l, Un, bc=bc, trials=u)
     # ...
 
-
 #==============================================================================
 # CLEAN UP SYMPY NAMESPACE
 #==============================================================================
diff --git a/sympde/expr/tests/test_tensor_2d.py b/sympde/expr/tests/test_tensor_2d.py
index 0f755cc..71580de 100644
--- a/sympde/expr/tests/test_tensor_2d.py
+++ b/sympde/expr/tests/test_tensor_2d.py
@@ -2,42 +2,20 @@
 
 # TODO: - add assert to every test
 
-import pytest
-
-from sympy import Symbol
 from sympy.core.containers import Tuple
-from sympy import symbols
-from sympy import IndexedBase
-from sympy import Matrix
-from sympy import Function
-from sympy import pi, cos, sin, exp
-from sympy import srepr
-from sympy import expand
-from sympy.physics.quantum import TensorProduct
-
-from sympde.core import Constant
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
-from sympde.calculus import laplace, hessian, bracket, convect
-from sympde.topology import (dx, dy, dz)
+
+from sympde.core     import Constant
+from sympde.calculus import grad, dot, curl, div
+#from sympde.calculus import laplace
+#from sympde.topology import dx
 from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
-from sympde.topology import Field, TestFunction
-from sympde.topology import ProductSpace
-from sympde.topology import InteriorDomain, Union
-from sympde.topology import Boundary, NormalVector, TangentVector
+from sympde.topology import elements_of
+from sympde.topology import Boundary
 from sympde.topology import Domain
-from sympde.topology import Trace, trace_0, trace_1
 from sympde.topology import Mapping
-from sympde.topology import Square
-from sympde.topology import ElementDomain
-from sympde.topology import Area
-
-from sympde.expr.expr import LinearExpr, BilinearExpr
-from sympde.expr.expr import LinearForm, BilinearForm, integral
-from sympde.expr.expr import DomainIntegral, BoundaryIntegral
-from sympde.expr.expr import Functional, Norm
-from sympde.expr.expr import linearize
+
+from sympde.expr.expr import BilinearForm, integral
 from sympde.expr.evaluation import TensorExpr
-from sympde.topology import elements_of
 
 
 #==============================================================================
@@ -52,8 +30,7 @@ def test_tensorize_2d_1():
     mu    = Constant('mu'   , is_real=True)
 
     V = ScalarFunctionSpace('V', domain)
-
-    u,v = [TestFunction(V, name=i) for i in ['u', 'v']]
+    u, v = elements_of(V, names='u, v')
 
     int_0 = lambda expr: integral(domain , expr)
 
@@ -80,8 +57,7 @@ def test_tensorize_2d_2():
     mu    = Constant('mu'   , is_real=True)
 
     V = VectorFunctionSpace('V', domain)
-
-    u,v = [TestFunction(V, name=i) for i in ['u', 'v']]
+    u, v = elements_of(V, names='u, v')
 
     int_0 = lambda expr: integral(domain , expr)
     # ...
@@ -108,8 +84,7 @@ def test_tensorize_2d_1_mapping():
     mu    = Constant('mu'   , is_real=True)
 
     V = ScalarFunctionSpace('V', domain)
-
-    u,v = [TestFunction(V, name=i) for i in ['u', 'v']]
+    u, v = elements_of(V, names='u, v')
 
     int_0 = lambda expr: integral(domain , expr)
     # ...
@@ -131,9 +106,7 @@ def test_tensorize_2d_2_mapping():
     domain = Domain('Omega', dim=DIM)
 
     V = VectorFunctionSpace('V', domain)
-
-    v = TestFunction(V, name='v')
-    u = TestFunction(V, name='u')
+    u, v = elements_of(V, names='u, v')
 
     c = Constant('c')
 
@@ -149,7 +122,7 @@ def test_tensorize_2d_3():
     domain = Domain('Omega', dim=2)
 
     V = ScalarFunctionSpace('V', domain)
-    u,v = elements_of(V, names='u,v')
+    u, v = elements_of(V, names='u,v')
 
     bx = Constant('bx')
     by = Constant('by')
diff --git a/sympde/exterior/tests/test_compiler.py b/sympde/exterior/tests/test_compiler.py
index c52d60a..09e38a5 100644
--- a/sympde/exterior/tests/test_compiler.py
+++ b/sympde/exterior/tests/test_compiler.py
@@ -2,28 +2,18 @@
 
 import pytest
 
-
-from sympy import symbols
-from sympy import Tuple
-from sympy import Matrix
-from sympy import srepr
-
-from sympde import Constant
-
 from sympde.topology import Domain
 from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
-from sympde.topology import ProductSpace
-from sympde.topology import H1Space, HcurlSpace, HdivSpace, L2Space, UndefinedSpace
-from sympde.topology import TestFunction
+from sympde.topology import element_of
 from sympde.topology import Field
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
+from sympde.calculus import grad, dot, inner, cross, curl, div
 
-from sympde.exterior import d, wedge, ip, ld, delta, jp, Ld
+#from sympde.exterior import d, wedge, ip, delta, jp
+from sympde.exterior import ld
 from sympde.exterior import DifferentialForm
 from sympde.exterior import ExteriorCalculusExpr, augmented_expression
 from sympde.calculus.errors import ArgumentTypeError
 
-
 #==============================================================================
 def test_compiler_3d_1():
 
@@ -37,7 +27,7 @@ def test_compiler_3d_1():
 
     X = H1 * Hcurl * Hdiv * L2
 
-    v0, v1, v2, v3 = TestFunction(X, ['v0', 'v1', 'v2', 'v3'])
+    v0, v1, v2, v3 = element_of(X, name='v0, v1, v2, v3')
 
     beta = Field(V, 'beta')
 
@@ -182,8 +172,8 @@ def test_compiler_3d_2():
 
     X = H1 * Hcurl * Hdiv * L2
 
-    v0, v1, v2, v3 = TestFunction(X, ['v0', 'v1', 'v2', 'v3'])
-    u0, u1, u2, u3 = TestFunction(X, ['u0', 'u1', 'u2', 'u3'])
+    v = element_of(X, name='v0, v1, v2, v3')
+    u = element_of(X, name='u0, u1, u2, u3')
 
     beta = Field(V, 'beta')
 
@@ -208,14 +198,14 @@ def test_compiler_3d_2():
 #    # ...
 
     # ... Mul operator
-    expr = u0*v0
-    print(ExteriorCalculusExpr(expr, tests=[v0]))
+    expr = u[0] * v[0]
+    print(ExteriorCalculusExpr(expr, tests=[v[0]]))
 
-    expr = u0*div(v2)
-    print(ExteriorCalculusExpr(expr, tests=[v2]))
+    expr = u[0] * div(v[2])
+    print(ExteriorCalculusExpr(expr, tests=[v[2]]))
 
-    expr = v0*div(u2)
-    print(ExteriorCalculusExpr(expr, tests=[v0]))
+    expr = v[0] * div(u[2])
+    print(ExteriorCalculusExpr(expr, tests=[v[0]]))
     # ...
 
 #    # ... Add operator
@@ -243,8 +233,8 @@ def test_compiler_3d_poisson():
 
     X = Hdiv * L2
 
-    sigma, u = TestFunction(X, ['sigma', 'u'])
-    tau,   v = TestFunction(X, [  'tau', 'v'])
+    sigma, u = element_of(X, name='sigma, u')
+    tau,   v = element_of(X, name='tau,   v')
 
     expr = dot(sigma, tau) + div(tau)*u + div(sigma)*v
     print(ExteriorCalculusExpr(expr, tests=[tau,v]))
@@ -264,8 +254,8 @@ def test_compiler_3d_stokes():
 
     X = Hdiv * L2
 
-    u,p = TestFunction(X, ['u', 'p'])
-    v,q = TestFunction(X, ['v', 'q'])
+    u, p = element_of(X, name='u, p')
+    v, q = element_of(X, name='v, q')
 
     with pytest.raises(ArgumentTypeError):
         expr = inner(grad(u), grad(v)) - div(v)*p + q*div(u)
@@ -277,8 +267,8 @@ def test_compiler_3d_stokes():
 
     X = Hdiv * L2
 
-    u,p = TestFunction(X, ['u', 'p'])
-    v,q = TestFunction(X, ['v', 'q'])
+    u, p = element_of(X, name='u, p')
+    v, q = element_of(X, name='v, q')
 
     expr = inner(grad(u), grad(v)) - div(v)*p + q*div(u)
     atoms = {u: DifferentialForm('u', index=2, dim=domain.dim),
diff --git a/sympde/printing/tests/test_latex.py b/sympde/printing/tests/test_latex.py
index 6e251ba..7bc91de 100644
--- a/sympde/printing/tests/test_latex.py
+++ b/sympde/printing/tests/test_latex.py
@@ -2,24 +2,14 @@
 
 from sympy import sin, cos, pi
 
-from sympde.core import Constant
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
-
-from sympde.topology import (dx, dy, dz)
+from sympde.calculus import grad, dot, inner, rot, div
+#from sympde.topology import (dx, dy, dz)
+from sympde.topology import Domain, Boundary
 from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
-from sympde.topology import ScalarField, VectorField
-from sympde.topology import ProductSpace
-from sympde.topology import ScalarTestFunction
-from sympde.topology import VectorTestFunction
-from sympde.topology import TestFunction
-from sympde.topology import Field
-from sympde.topology import Unknown
-from sympde.topology import Domain, Boundary, NormalVector, TangentVector
-#from sympde.expr import tensorize
-from sympde.expr import BilinearForm, LinearForm, integral
+from sympde.topology import element_of
+from sympde.expr     import BilinearForm, LinearForm, integral
 from sympde.exterior import d, wedge, ip, jp, delta, hodge
 from sympde.exterior import DifferentialForm
-
 from sympde.printing.latex import latex
 
 
@@ -33,10 +23,9 @@ def test_latex_1d():
 
     x = V.coordinates
 
-    v = ScalarTestFunction(V, name='v')
-    u = ScalarTestFunction(V, name='u')
-
-    F = ScalarField(V, name='F')
+    v = element_of(V, name='v')
+    u = element_of(V, name='u')
+#    F = element_of(V, name='F')
 
     int_0 = lambda expr: integral(domain , expr)
     
@@ -65,10 +54,9 @@ def test_latex_2d_1():
 
     x,y = V.coordinates
 
-    v = ScalarTestFunction(V, name='v')
-    u = ScalarTestFunction(V, name='u')
-
-    F = ScalarField(V, name='F')
+    v = element_of(V, name='v')
+    u = element_of(V, name='u')
+#    F = element_of(V, name='F')
     
     int_0 = lambda expr: integral(domain , expr)
 
@@ -97,10 +85,10 @@ def test_latex_2d_2():
 
     x,y = V.coordinates
 
-    v = VectorTestFunction(V, name='v')
-    u = VectorTestFunction(V, name='u')
+    v = element_of(V, name='v')
+    u = element_of(V, name='u')
+#    F = element_of(V, name='F')
 
-    F = VectorField(V, name='F')
     int_0 = lambda expr: integral(domain , expr)
 
     assert(latex(v) == r'\mathbf{v}')
@@ -125,10 +113,9 @@ def test_latex_3d_1():
 
     x,y,z = V.coordinates
 
-    v = ScalarTestFunction(V, name='v')
-    u = ScalarTestFunction(V, name='u')
-
-    F = ScalarField(V, name='F')
+    v = element_of(V, name='v')
+    u = element_of(V, name='u')
+#    F = element_of(V, name='F')
 
     int_0 = lambda expr: integral(domain , expr)
     
@@ -157,10 +144,8 @@ def test_latex_3d_2():
 
     x,y,z = V.coordinates
 
-    v = VectorTestFunction(V, name='v')
-    u = VectorTestFunction(V, name='u')
-
-    F = VectorField(V, name='F')
+    v = element_of(V, name='v')
+    u = element_of(V, name='u')
 
     int_0 = lambda expr: integral(domain , expr)
     
@@ -191,8 +176,8 @@ def test_latex_2d_3():
 
     x = V.coordinates
 
-    v = ScalarTestFunction(V, name='v')
-    u = ScalarTestFunction(V, name='u')
+    v = element_of(V, name='v')
+    u = element_of(V, name='u')
 
     int_0 = lambda expr: integral(domain , expr)
     int_1 = lambda expr: integral(B1, expr)
@@ -223,10 +208,10 @@ def test_latex_2d_4():
     V = VectorFunctionSpace('V', domain)
     W = ScalarFunctionSpace('W', domain)
 
-    v = VectorTestFunction(V, name='v')
-    u = VectorTestFunction(V, name='u')
-    p = ScalarTestFunction(W, name='p')
-    q = ScalarTestFunction(W, name='q')
+    v = element_of(V, name='v')
+    u = element_of(V, name='u')
+    p = element_of(W, name='p')
+    q = element_of(W, name='q')
 
     int_0 = lambda expr: integral(domain , expr)
 
@@ -249,9 +234,8 @@ def test_latex_2d_5():
     # ... abstract model
     W1 = VectorFunctionSpace('W1', domain)
 
-    w1 = VectorTestFunction(W1, name='w1')
-
-    F = VectorField(W1, 'F')
+    w1 = element_of(W1, name='w1')
+    F  = element_of(W1, 'F')
 
     int_0 = lambda expr: integral(domain , expr)
     
@@ -292,7 +276,7 @@ def test_latex_ec_3d_1():
     domain = Domain('Omega', dim=3)
     V = VectorFunctionSpace('V', domain)
 
-    beta = Field(V, 'beta')
+    beta = element_of(V, 'beta')
     # ...
 
     print(latex(u_0))
diff --git a/sympde/topology/tests/test_field_2d.py b/sympde/topology/tests/test_field_2d.py
index 888814f..677b079 100644
--- a/sympde/topology/tests/test_field_2d.py
+++ b/sympde/topology/tests/test_field_2d.py
@@ -1,25 +1,14 @@
 # coding: utf-8
 
-import pytest
-
-from sympy import Symbol
-from sympy.core.containers import Tuple
-from sympy import symbols
-from sympy import IndexedBase
 from sympy import Matrix
-from sympy import Function
-from sympy import pi, cos, sin
-from sympy import srepr
-from sympy.physics.quantum import TensorProduct
-
-from sympde.core import Constant
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
-from sympde.calculus import laplace, hessian
+
+from sympde.calculus import grad, inner
 from sympde.topology import Domain
-from sympde.topology import dx, dy, dz
-from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
-from sympde.topology import VectorField
+from sympde.topology import dx
+from sympde.topology import VectorFunctionSpace
+from sympde.topology import element_of
 
+#==============================================================================
 
 DIM = 2
 domain = Domain('Omega', dim=DIM)
@@ -28,10 +17,10 @@ domain = Domain('Omega', dim=DIM)
 def test_field_2d_1():
     print('============ test_field_2d_1 =============')
 
-    W = VectorFunctionSpace('W', domain)
-    x,y = W.coordinates
+#    x, y = domain.coordinates
 
-    F = VectorField(W, 'F')
+    W = VectorFunctionSpace('W', domain)
+    F = element_of(W, 'F')
 
     assert( dx(F) == Matrix([[dx(F[0]), dx(F[1])]]) )
 
diff --git a/sympde/topology/tests/test_logical_expr.py b/sympde/topology/tests/test_logical_expr.py
index 4361c2d..3a5d486 100644
--- a/sympde/topology/tests/test_logical_expr.py
+++ b/sympde/topology/tests/test_logical_expr.py
@@ -1,26 +1,18 @@
 # coding: utf-8
 
-from sympy.core.containers import Tuple
+from sympy import Symbol, symbols
 from sympy import Matrix
-from sympy.tensor import IndexedBase
-from sympy import symbols, simplify, Symbol
 from sympy import expand
 from sympy import cos, sin, sqrt, pi
 
-from sympde.topology import Mapping, DetJacobian
-from sympde.topology import Domain
-
-from sympde.core import Constant
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
-from sympde.calculus import laplace, hessian, bracket, convect
-from sympde.topology import (dx, dy, dz)
-from sympde.topology import (dx1, dx2, dx3)
+from sympde.core     import Constant
+from sympde.topology import Domain, Mapping, DetJacobian
+from sympde.topology import dx, dy
+from sympde.topology import dx1, dx2, dx3
 from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
-from sympde.topology import ScalarTestFunction
-from sympde.topology import VectorTestFunction
+from sympde.topology import element_of
 from sympde.topology import LogicalExpr
 from sympde.topology import SymbolicExpr
-from sympde.topology import element_of
 from sympde.topology import IdentityMapping
 from sympde.topology import PolarMapping
 from sympde.topology import TargetMapping
@@ -40,7 +32,7 @@ def test_logical_expr_1d_1():
 
     V = ScalarFunctionSpace('V', domain)
 
-    u,v = [ScalarTestFunction(V, name=i) for i in ['u', 'v']]
+    u,v = [element_of(V, name=i) for i in ['u', 'v']]
 
     det_M = DetJacobian(M)
     #print('det = ', det_M)
@@ -87,7 +79,7 @@ def test_symbolic_expr_1d_1():
 
     V = ScalarFunctionSpace('V', domain)
 
-    u,v = [ScalarTestFunction(V, name=i) for i in ['u', 'v']]
+    u = element_of(V, name='u')
 
     det_M = DetJacobian(M)
     det_M = SymbolicExpr(det_M)
@@ -159,8 +151,8 @@ def test_logical_expr_2d_1():
     V = ScalarFunctionSpace('V', domain)
     W = VectorFunctionSpace('V', domain)
 
-    u,v = [ScalarTestFunction(V, name=i) for i in ['u', 'v']]
-    w = VectorTestFunction(W, name='w')
+    u,v = [element_of(V, name=i) for i in ['u', 'v']]
+    w = element_of(W, name='w')
 
     det_M = DetJacobian(M)
     #print('det = ', det_M)
@@ -221,7 +213,7 @@ def test_symbolic_expr_2d_1():
 
     V = ScalarFunctionSpace('V', domain)
 
-    u,v = [ScalarTestFunction(V, name=i) for i in ['u', 'v']]
+    u = element_of(V, name='u')
 
     det_M = DetJacobian(M)
     det_M = SymbolicExpr(det_M)
@@ -299,7 +291,7 @@ def test_logical_expr_3d_1():
 
     V = ScalarFunctionSpace('V', domain)
 
-    u,v = [ScalarTestFunction(V, name=i) for i in ['u', 'v']]
+    u,v = [element_of(V, name=i) for i in ['u', 'v']]
 
     det_M = DetJacobian(M)
     #print('det = ', det_M)
@@ -352,8 +344,7 @@ def test_symbolic_expr_3d_1():
     alpha = Constant('alpha')
 
     V = ScalarFunctionSpace('V', domain)
-
-    u,v = [ScalarTestFunction(V, name=i) for i in ['u', 'v']]
+    u = element_of(V, 'u')
 
     det_M = DetJacobian(M)
     det_M = SymbolicExpr(det_M)
diff --git a/sympde/topology/tests/test_space.py b/sympde/topology/tests/test_space.py
index 3bc5614..32e2242 100644
--- a/sympde/topology/tests/test_space.py
+++ b/sympde/topology/tests/test_space.py
@@ -1,12 +1,12 @@
 # coding: utf-8
 
-from sympde.core     import Constant
-from sympde.calculus import grad, curl, div
+from sympde.calculus import grad, div
 from sympde.topology import Domain
 from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
-from sympde.topology import ProductSpace
+#from sympde.topology import ProductSpace
+from sympde.topology import element_of
 from sympde.topology import H1Space, HcurlSpace, HdivSpace, L2Space, UndefinedSpace
-from sympde.topology import TestFunction, ScalarTestFunction, VectorTestFunction
+from sympde.topology import ScalarTestFunction, VectorTestFunction
 from sympde.topology import Field, ScalarField, VectorField
 from sympde.topology import Projector
 
@@ -65,18 +65,21 @@ def test_space_2d_1():
     # ...
 
     # ...
-    v1, v2, v3 = TestFunction(V, ['v1', 'v2', 'v3'])
+    v1, v2, v3 = element_of(V, 'v1, v2, v3')
 
     assert(isinstance(v1, ScalarTestFunction))
     assert(isinstance(v2, ScalarTestFunction))
     assert(isinstance(v3, ScalarTestFunction))
+    assert(v1.space is V1)
+    assert(v2.space is V2)
+    assert(v3.space is V3)
     assert(v1.name == 'v1')
     assert(v2.name == 'v2')
     assert(v3.name == 'v3')
     # ...
 
     # ...
-    u1, u2, u3 = TestFunction(U, ['u1', 'u2', 'u3'])
+    u1, u2, u3 = element_of(U, 'u1, u2, u3')
 
     assert(isinstance(u1, ScalarTestFunction))
     assert(isinstance(u2, ScalarTestFunction))

pip install -e .
parse_log_pytest
pytest
pytest
python -m pip install --upgrade pip python -m pip uninstall -y sympde
3.6
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
attrs @ file:///opt/conda/conda-bld/attrs_1642510447205/work
cached-property==1.5.2
certifi==2021.5.30
h5py==3.1.0
importlib-metadata @ file:///tmp/build/80754af9/importlib-metadata_1631916693255/work
iniconfig @ file:///home/linux1/recipes/ci/iniconfig_1610983019677/work
more-itertools @ file:///tmp/build/80754af9/more-itertools_1637733554872/work
mpmath==1.3.0
numpy==1.19.5
packaging @ file:///tmp/build/80754af9/packaging_1637314298585/work
pluggy @ file:///tmp/build/80754af9/pluggy_1615976315926/work
py @ file:///opt/conda/conda-bld/py_1644396412707/work
pyparsing @ file:///tmp/build/80754af9/pyparsing_1635766073266/work
pytest==6.2.4
PyYAML==6.0.1
-e git+https://github.com/pyccel/sympde.git@c7461ffa7de44c968d2fd9c9996a8fdcfa7debcf#egg=sympde
sympy==1.9
toml @ file:///tmp/build/80754af9/toml_1616166611790/work
typing_extensions @ file:///opt/conda/conda-bld/typing_extensions_1647553014482/work
yamlloader==1.2.2
zipp @ file:///tmp/build/80754af9/zipp_1633618647012/work

name: sympde
channels:
  - defaults
  - https://repo.anaconda.com/pkgs/main
  - https://repo.anaconda.com/pkgs/r
  - conda-forge
dependencies:
  - _libgcc_mutex=0.1=main
  - _openmp_mutex=5.1=1_gnu
  - attrs=21.4.0=pyhd3eb1b0_0
  - ca-certificates=2025.2.25=h06a4308_0
  - certifi=2021.5.30=py36h06a4308_0
  - importlib-metadata=4.8.1=py36h06a4308_0
  - importlib_metadata=4.8.1=hd3eb1b0_0
  - iniconfig=1.1.1=pyhd3eb1b0_0
  - ld_impl_linux-64=2.40=h12ee557_0
  - libffi=3.3=he6710b0_2
  - libgcc-ng=11.2.0=h1234567_1
  - libgomp=11.2.0=h1234567_1
  - libstdcxx-ng=11.2.0=h1234567_1
  - more-itertools=8.12.0=pyhd3eb1b0_0
  - ncurses=6.4=h6a678d5_0
  - openssl=1.1.1w=h7f8727e_0
  - packaging=21.3=pyhd3eb1b0_0
  - pluggy=0.13.1=py36h06a4308_0
  - py=1.11.0=pyhd3eb1b0_0
  - pyparsing=3.0.4=pyhd3eb1b0_0
  - pytest=6.2.4=py36h06a4308_2
  - python=3.6.13=h12debd9_1
  - readline=8.2=h5eee18b_0
  - setuptools=58.0.4=py36h06a4308_0
  - sqlite=3.45.3=h5eee18b_0
  - tk=8.6.14=h39e8969_0
  - toml=0.10.2=pyhd3eb1b0_0
  - typing_extensions=4.1.1=pyh06a4308_0
  - wheel=0.37.1=pyhd3eb1b0_0
  - xz=5.6.4=h5eee18b_1
  - zipp=3.6.0=pyhd3eb1b0_0
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - cached-property==1.5.2
      - h5py==3.1.0
      - mpmath==1.3.0
      - numpy==1.19.5
      - pip==21.3.1
      - pyyaml==6.0.1
      - sympy==1.9
      - yamlloader==1.2.2
prefix: /opt/conda/envs/sympde

sympde/core/tests/test_utils.py::test_expand_name_patterns
sympde/expr/tests/test_expr_2d.py::test_terminal_expr_linear_2d_1 sympde/expr/tests/test_expr_2d.py::test_terminal_expr_linear_2d_3 sympde/expr/tests/test_expr_2d.py::test_terminal_expr_bilinear_2d_1 sympde/expr/tests/test_expr_2d.py::test_terminal_expr_bilinear_2d_3 sympde/expr/tests/test_expr_2d.py::test_linearize_form_2d_1 sympde/expr/tests/test_expr_2d.py::test_linearize_form_2d_3 sympde/expr/tests/test_expr_2d.py::test_area_2d_1 sympde/expr/tests/test_expr_2d.py::test_stabilization_2d_1 sympde/expr/tests/test_expr_2d.py::test_user_function_2d_1 sympde/expr/tests/test_expr_2d.py::test_functional_2d_1 sympde/expr/tests/test_expr_2d.py::test_norm_2d_1 sympde/expr/tests/test_expr_2d.py::test_interface_integral_1 sympde/expr/tests/test_expr_2d.py::test_interface_integral_2 sympde/expr/tests/test_expr_2d.py::test_interface_integral_3 sympde/expr/tests/test_tensor_2d.py::test_tensorize_2d_1 sympde/expr/tests/test_tensor_2d.py::test_tensorize_2d_1_mapping sympde/expr/tests/test_tensor_2d.py::test_tensorize_2d_3 sympde/printing/tests/test_latex.py::test_latex_2d_2 sympde/printing/tests/test_latex.py::test_latex_3d_2
sympde/topology/mapping.py sympde/topology/space.py sympde/core/utils.py sympde/expr/expr.py sympde/topology/domain.py
sympde/topology/mapping.py sympde/topology/space.py sympde/core/utils.py sympde/expr/expr.py sympde/topology/domain.py