diff --git a/.travis.yml b/.travis.yml
index 066eb8b57..11db5fef7 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -24,24 +24,16 @@ before_install:
   - h5pcc -showconfig -echo || true
   - export CC="mpicc" HDF5_MPI="ON" HDF5_DIR=/usr/lib/x86_64-linux-gnu/hdf5/openmpi
   - python -m pip install --upgrade pip
-  - python -m pip install -r requirements.txt
   - python -m pip uninstall -y psydac
   - python -m pip uninstall -y sympde
-  - python -m pip uninstall -y gelato
   - python -m pip uninstall -y pyccel
-  - wget https://raw.githubusercontent.com/pyccel/pyccel/master/requirements.txt -O requirements_pyccel.txt
-  - python -m pip install -r requirements_pyccel.txt
-  - python -m pip install git+https://github.com/pyccel/pyccel@master
-  - wget https://raw.githubusercontent.com/pyccel/sympde/master/requirements.txt -O requirements_sympde.txt
-  - python -m pip install -r requirements_sympde.txt
-  - python -m pip install git+https://github.com/pyccel/sympde@master
-  - wget https://raw.githubusercontent.com/pyccel/gelato/master/requirements.txt -O requirements_gelato.txt
-  - python -m pip install -r requirements_gelato.txt
-  - python -m pip install git+https://github.com/pyccel/gelato@master
+  - python -m pip uninstall -y gelato
+  - python -m pip install -r requirements.txt
 
 # command to install project
 install:
   - python -m pip install .
+  - python -m pip freeze
 
 before_script:
   - mkdir pytest
diff --git a/TODO.rst b/TODO.rst
index 497033ab3..7119d5764 100644
--- a/TODO.rst
+++ b/TODO.rst
@@ -50,7 +50,7 @@ API
 
 - add other solvers to the solver_driver (only cg is available now)
 
-- remove psydac/run_tests.sh. however, we need to clear the cache of sympy after 1d, 2d and 3d tests, otherwise pytest will crash.
+- however, we need to clear the cache of sympy after 1d, 2d and 3d tests, otherwise pytest will crash.
 
 - add sympde and pyccel install procedure; maybe wget to download the requierements files as requirements_pyccel.txt etc then call pip3
 
diff --git a/psydac/ddm/plot_cart_poisson.py b/examples/poisson_2d_jacobi_iteration/plot_cart_poisson.py
similarity index 100%
rename from psydac/ddm/plot_cart_poisson.py
rename to examples/poisson_2d_jacobi_iteration/plot_cart_poisson.py
diff --git a/psydac/ddm/test_cart_poisson.py b/examples/poisson_2d_jacobi_iteration/test_cart_poisson.py
similarity index 100%
rename from psydac/ddm/test_cart_poisson.py
rename to examples/poisson_2d_jacobi_iteration/test_cart_poisson.py
diff --git a/psydac/__init__.py b/psydac/__init__.py
index c498e6668..36d218598 100644
--- a/psydac/__init__.py
+++ b/psydac/__init__.py
@@ -1,8 +1,9 @@
 # -*- coding: UTF-8 -*-
-__version__ = "0.1"
-__all__     = ['api', 'cad', 'core', 'ddm', 'feec', 'fem',
+__all__     = ['__version__', 'api', 'cad', 'core', 'ddm', 'feec', 'fem',
                'linalg', 'mapping', 'utilities']
 
+from psydac.version import __version__
+
 from psydac import api
 from psydac import cad
 from psydac import core
diff --git a/psydac/api/ast/evaluation.py b/psydac/api/ast/evaluation.py
index e70012448..3a82d20a0 100644
--- a/psydac/api/ast/evaluation.py
+++ b/psydac/api/ast/evaluation.py
@@ -7,7 +7,10 @@
 from sympde.topology             import ScalarField
 from sympde.topology             import VectorField
 from sympde.topology             import SymbolicExpr
-from sympde.topology.space       import ScalarTestFunction
+from sympde.topology.space       import element_of
+from sympde.topology.space       import ScalarFunctionSpace
+from sympde.topology.space       import VectorFunctionSpace
+from sympde.topology.datatype    import dtype_space_registry
 from sympde.topology.derivatives import _logical_partial_derivatives
 
 from pyccel.ast.core      import IndexedVariable
@@ -35,6 +38,21 @@ def __new__(cls, space, mapping, discrete_boundary=None, name=None,
         if not isinstance(mapping, Mapping):
             raise TypeError('> Expecting a Mapping object')
 
+        #.....................................................................
+        # If vector space is of undefined type, we assume that each component
+        # lives in H1; otherwise we raise an error. TODO: improve!
+        if isinstance(space, VectorFunctionSpace):
+            if space.kind == dtype_space_registry['undefined']:
+                space = ScalarFunctionSpace(
+                    name   = space.name,
+                    domain = space.domain,
+                    kind   = 'H1'
+                )
+            else:
+                msg = 'Cannot evaluate vector spaces of kind {}'.format(space.kind)
+                raise NotImplementedError(msg)
+        #.....................................................................
+
         obj = SplBasic.__new__(cls, mapping, name=name,
                                prefix='eval_mapping', mapping=mapping,
                                is_rational_mapping=is_rational_mapping)
@@ -173,7 +191,7 @@ def _initialize(self):
         weights_elements = []
         if self.is_rational_mapping:
             # TODO check if 'w' exist already
-            weights_pts = ScalarField(self.space, name='w')
+            weights_pts = element_of(self.space, name='w')
 
             weights_elements = [weights_pts]
 
@@ -217,7 +235,7 @@ def _initialize(self):
         # ...
 
         # ...
-        Nj           = ScalarTestFunction(space, name='Nj')
+        Nj           = element_of(space, name='Nj')
         body         = []
         init_basis   = OrderedDict()
         updates      = []
@@ -364,7 +382,7 @@ def _initialize(self):
         # ...
 
         # ...
-        Nj = ScalarTestFunction(space, name='Nj')
+        Nj = element_of(space, name='Nj')
         body = []
         init_basis = OrderedDict()
         init_map   = OrderedDict()
@@ -503,7 +521,7 @@ def _initialize(self):
         # ...
 
         # ...
-        Nj = VectorField(space, name='Nj')
+        Nj = VectorField(space, name='Nj')  # TODO: use 'element_of'
         body = []
         init_basis = OrderedDict()
         init_map   = OrderedDict()
@@ -666,7 +684,7 @@ def _initialize(self):
         # ...
 
         # ...
-        Nj = ScalarTestFunction(space, name='Nj')
+        Nj = element_of(space, name='Nj')
         init_basis = OrderedDict()
         init_map   = OrderedDict()
 
@@ -860,7 +878,7 @@ def _initialize(self):
         weights_elements = []
         if self.is_rational_mapping:
             # TODO check if 'w' exist already
-            weights_pts = ScalarField(self.space, name='w')
+            weights_pts = element_of(self.space, name='w')
 
             weights_elements = [weights_pts]
 
@@ -905,7 +923,7 @@ def _initialize(self):
         # ...
 
         # ...
-        Nj             = ScalarTestFunction(space, name='Nj')
+        Nj             = element_of(space, name='Nj')
         body           = []
         init_basis     = OrderedDict()
         atomic_exprs   = self.elements + weights_elements
diff --git a/psydac/api/ast/fem.py b/psydac/api/ast/fem.py
index 56540ac5a..03bb31326 100644
--- a/psydac/api/ast/fem.py
+++ b/psydac/api/ast/fem.py
@@ -39,6 +39,7 @@
 from sympde.topology.space       import ProductSpace
 from sympde.topology.space       import ScalarTestFunction
 from sympde.topology.space       import VectorTestFunction
+from sympde.topology.space       import element_of
 from sympde.topology.space       import IndexedTestTrial
 from sympde.topology.derivatives import _partial_derivatives
 from sympde.topology.derivatives import _logical_partial_derivatives
@@ -617,13 +618,11 @@ def _initialize(self):
         if is_bilinear or is_linear:
             test_function = self.weak_form.test_functions
             if not isinstance(test_function, (tuple, Tuple)):
-                test_function = [test_function]
-                test_function = Tuple(*test_function)
+                test_function = Tuple(test_function)
 
         elif is_function:
-            test_function = ScalarTestFunction(self.weak_form.space, name='Nj')
-            test_function = [test_function]
-            test_function = Tuple(*test_function)
+            test_function = element_of(self.weak_form.space, name='Nj')
+            test_function = Tuple(test_function)
 
         # creation of symbolic vars
         if is_bilinear:
diff --git a/psydac/api/ast/utilities.py b/psydac/api/ast/utilities.py
index b0575685d..e6af6f014 100644
--- a/psydac/api/ast/utilities.py
+++ b/psydac/api/ast/utilities.py
@@ -12,6 +12,7 @@
 from sympde.topology.space       import ScalarTestFunction
 from sympde.topology.space       import VectorTestFunction
 from sympde.topology.space       import IndexedTestTrial
+from sympde.topology.space       import element_of
 from sympde.topology             import ScalarField
 from sympde.topology             import VectorField, IndexedVectorField
 from sympde.topology             import Mapping
@@ -462,7 +463,7 @@ def rationalize_eval_mapping(mapping, nderiv, space, indices_quad):
 
     # ... weights and their derivatives
     # TODO check if 'w' exist already
-    weights = ScalarField(space, name='w')
+    weights = element_of(space, name='w')
 
     weights_elements = [weights]
     if nderiv > 0:
diff --git a/psydac/api/tests/test_api_1d_compatible_spaces.py b/psydac/api/tests/test_api_1d_compatible_spaces.py
index 9c8d83ab2..3c199fa3c 100644
--- a/psydac/api/tests/test_api_1d_compatible_spaces.py
+++ b/psydac/api/tests/test_api_1d_compatible_spaces.py
@@ -1,37 +1,21 @@
 # -*- coding: UTF-8 -*-
 
-from sympy import pi, cos, sin, Tuple, Matrix
+from sympy               import pi, sin
+from scipy.sparse.linalg import spsolve
 
-from sympde.core import Constant
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
-from sympde.calculus import laplace, hessian
-from sympde.topology import (dx, dy, dz)
-from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
+from sympde.calculus import dot, div
+from sympde.topology import ScalarFunctionSpace
 from sympde.topology import ProductSpace
 from sympde.topology import element_of
-from sympde.topology import Boundary, NormalVector, TangentVector
-from sympde.topology import Domain, Line, Square, Cube
-from sympde.topology import Trace, trace_0, trace_1
-from sympde.topology import Union
-from sympde.expr import BilinearForm, LinearForm
-from sympde.expr import integral
-from sympde.expr import Norm, TerminalExpr
-from sympde.expr import find, EssentialBC
-
-
-from psydac.fem.basic   import FemField
-from psydac.fem.vector   import VectorFemField
-from psydac.api.discretization import discretize
-
-from numpy import linspace, zeros, allclose
-from mpi4py import MPI
-import pytest
+from sympde.topology import Line
+from sympde.expr     import BilinearForm, LinearForm
+from sympde.expr     import integral
+from sympde.expr     import Norm
+from sympde.expr     import find
 
-from scipy.sparse.linalg import cg, gmres
-from scipy import linalg
-from scipy.sparse.linalg import spsolve
+from psydac.fem.basic          import FemField
+from psydac.api.discretization import discretize
 
-import matplotlib.pyplot as plt
 #==============================================================================
 
 def run_system_1_1d_dir(f0, sol, ncells, degree):
@@ -76,7 +60,7 @@ def run_system_1_1d_dir(f0, sol, ncells, degree):
 
     ah.assemble()
 
-    M   = ah.linear_system.lhs.tosparse()
+    M   = ah.linear_system.lhs.tosparse().tocsc()
     rhs = ah.linear_system.rhs.toarray()
     sol = spsolve(M, rhs)
 
diff --git a/psydac/api/tests/test_api_2d_compatible_spaces.py b/psydac/api/tests/test_api_2d_compatible_spaces.py
index 154745f41..9c8607e04 100644
--- a/psydac/api/tests/test_api_2d_compatible_spaces.py
+++ b/psydac/api/tests/test_api_2d_compatible_spaces.py
@@ -1,39 +1,23 @@
 # -*- coding: UTF-8 -*-
 
-from sympy import pi, cos, sin, Tuple, Matrix
+import numpy as np
+from sympy import pi, cos, sin, Matrix
+from scipy import linalg
+from scipy.sparse.linalg import spsolve
 
-from sympde.core import Constant
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
-from sympde.calculus import laplace, hessian
-from sympde.topology import (dx, dy, dz)
+from sympde.calculus import grad, dot, inner, div
 from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
 from sympde.topology import ProductSpace
 from sympde.topology import element_of
-from sympde.topology import Boundary, NormalVector, TangentVector
-from sympde.topology import Domain, Line, Square, Cube
-from sympde.topology import Trace, trace_0, trace_1
-from sympde.topology import Union
+from sympde.topology import Square
 from sympde.expr import BilinearForm, LinearForm, integral
-from sympde.expr import Norm, TerminalExpr
+from sympde.expr import Norm
 from sympde.expr import find, EssentialBC
 
-
-from psydac.fem.basic   import FemField
-from psydac.fem.vector  import VectorFemField
+from psydac.fem.basic          import FemField
+from psydac.fem.vector         import VectorFemField
 from psydac.api.discretization import discretize
 
-from numpy import linspace, zeros, allclose
-import numpy as np
-from mpi4py import MPI
-import pytest
-
-from scipy.sparse.linalg import cg, gmres
-from scipy.sparse.linalg import spsolve
-from scipy import linalg
-
-import matplotlib.pyplot as plt
-from mpl_toolkits.mplot3d import Axes3D
-from matplotlib import animation
 #==============================================================================
 
 def run_system_1_2d_dir(f0, sol, ncells, degree):
@@ -83,7 +67,7 @@ def run_system_1_2d_dir(f0, sol, ncells, degree):
 
     # ...
     ah.assemble()
-    M   = ah.linear_system.lhs.tosparse()
+    M   = ah.linear_system.lhs.tosparse().tocsc()
     rhs = ah.linear_system.rhs.toarray()
     x   = spsolve(M, rhs)
     # ...
@@ -194,11 +178,17 @@ def run_system_2_2d_dir(f1, f2,u1, u2, ncells, degree):
 #            SERIAL TESTS
 ###############################################################################
 
+def test_api_system_1_2d_dir_1():
+    from sympy.abc import x,y
+
+    f0 =  -2*(2*pi)**2*sin(2*pi*x)*sin(2*pi*y)
+    u  = sin(2*pi*x)*sin(2*pi*y)
+    x = run_system_1_2d_dir(f0,u, ncells=[10,10], degree=[2,2])
+
 #==============================================================================
 def test_api_system_2_2d_dir_1():
     from sympy.abc import x,y
-    from sympy import cos
-    
+
     f1 = -x**2*(x - 1)**2*(24*y - 12) - 4*y*(x**2 + 4*x*(x - 1) + (x - 1)**2)*(2*y**2 - 3*y + 1) - 2*pi*cos(2*pi*x)
     f2 = 4*x*(2*x**2 - 3*x + 1)*(y**2 + 4*y*(y - 1) + (y - 1)**2) + y**2*(24*x - 12)*(y - 1)**2 + 2*pi*cos(2*pi*y)
     u1 = x**2*(-x + 1)**2*(4*y**3 - 6*y**2 + 2*y)
@@ -206,11 +196,3 @@ def test_api_system_2_2d_dir_1():
     p  = sin(2*pi*x) - sin(2*pi*y)
     
     x = run_system_2_2d_dir(f1, f2, u1, u2, ncells=[2**3,2**3], degree=[2,2])
-            
-def test_api_system_1_2d_dir_1():
-    from sympy.abc import x,y
-
-    f0 =  -2*(2*pi)**2*sin(2*pi*x)*sin(2*pi*y)
-    u  = sin(2*pi*x)*sin(2*pi*y)
-    x = run_system_1_2d_dir(f0,u, ncells=[10,10], degree=[2,2])
-
diff --git a/psydac/api/tests/test_api_2d_scalar.py b/psydac/api/tests/test_api_2d_scalar.py
index 95cb3b718..01ea0e0bb 100644
--- a/psydac/api/tests/test_api_2d_scalar.py
+++ b/psydac/api/tests/test_api_2d_scalar.py
@@ -1,31 +1,24 @@
 # -*- coding: UTF-8 -*-
 
+from mpi4py import MPI
 from sympy import pi, cos, sin
 from sympy.utilities.lambdify import implemented_function
+import pytest
 
-from sympde.core import Constant
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
-from sympde.calculus import laplace, hessian
-from sympde.topology import (dx, dy, dz)
-from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
-from sympde.topology import ProductSpace
+from sympde.calculus import grad, dot
+from sympde.calculus import laplace
+from sympde.topology import ScalarFunctionSpace
 from sympde.topology import element_of
-from sympde.topology import Boundary, NormalVector, TangentVector
-from sympde.topology import Domain, Line, Square, Cube
-from sympde.topology import Trace, trace_0, trace_1
+from sympde.topology import NormalVector
+from sympde.topology import Square
 from sympde.topology import Union
-from sympde.expr import BilinearForm, LinearForm, integral
-from sympde.expr import Norm
-from sympde.expr import find, EssentialBC
+from sympde.expr     import BilinearForm, LinearForm, integral
+from sympde.expr     import Norm
+from sympde.expr     import find, EssentialBC
 
-from psydac.fem.basic   import FemField
+from psydac.fem.basic          import FemField
 from psydac.api.discretization import discretize
 
-from numpy import linspace, zeros, allclose
-from mpi4py import MPI
-import pytest
-
-
 #==============================================================================
 def run_poisson_2d_dir(solution, f, ncells, degree, comm=None):
 
diff --git a/psydac/api/tests/test_api_2d_scalar_mapping.py b/psydac/api/tests/test_api_2d_scalar_mapping.py
index 6c73888b5..69ba30c5d 100644
--- a/psydac/api/tests/test_api_2d_scalar_mapping.py
+++ b/psydac/api/tests/test_api_2d_scalar_mapping.py
@@ -1,40 +1,24 @@
 # -*- coding: UTF-8 -*-
 
+from mpi4py import MPI
 from sympy import pi, cos, sin
-from sympy import S
-from sympy import Tuple
-from sympy import Matrix
-
-from sympde.core import Constant
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
-from sympde.calculus import laplace, hessian
-from sympde.topology import (dx, dy, dz)
-from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
-from sympde.topology import ProductSpace
+import pytest
+import os
+
+from sympde.calculus import grad, dot
+from sympde.calculus import laplace
+from sympde.topology import ScalarFunctionSpace
 from sympde.topology import element_of
-from sympde.topology import Unknown
-from sympde.topology import InteriorDomain, Union
-from sympde.topology import Boundary, NormalVector, TangentVector
-from sympde.topology import Domain, Line, Square, Cube
-from sympde.topology import Trace, trace_0, trace_1
+from sympde.topology import NormalVector
+from sympde.topology import Domain
 from sympde.topology import Union
-from sympde.topology import Mapping
 from sympde.expr import BilinearForm, LinearForm, integral
 from sympde.expr import Norm
 from sympde.expr import find, EssentialBC
 
-from psydac.fem.basic   import FemField
-from psydac.fem.vector  import ProductFemSpace, VectorFemField
+from psydac.fem.basic          import FemField
 from psydac.api.discretization import discretize
 
-from psydac.mapping.discrete import SplineMapping
-
-from numpy import linspace, zeros, allclose
-from mpi4py import MPI
-import pytest
-
-import os
-
 # ... get the mesh directory
 try:
     mesh_dir = os.environ['PSYDAC_MESH_DIR']
diff --git a/psydac/api/tests/test_api_2d_system.py b/psydac/api/tests/test_api_2d_system.py
index 880537398..29cdf30be 100644
--- a/psydac/api/tests/test_api_2d_system.py
+++ b/psydac/api/tests/test_api_2d_system.py
@@ -1,31 +1,23 @@
 # -*- coding: UTF-8 -*-
 
+from mpi4py import MPI
 from sympy import pi, cos, sin, Tuple, Matrix
+import numpy as np
+import pytest
 
-from sympde.core import Constant
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
-from sympde.calculus import laplace, hessian
-from sympde.topology import (dx, dy, dz)
+from sympde.calculus import grad, dot, inner
 from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
 from sympde.topology import ProductSpace
 from sympde.topology import element_of
-from sympde.topology import Boundary, NormalVector, TangentVector
-from sympde.topology import Domain, Line, Square, Cube
-from sympde.topology import Trace, trace_0, trace_1
-from sympde.topology import Union
+from sympde.topology import Square
 from sympde.expr import BilinearForm, LinearForm, integral
 from sympde.expr import Norm
 from sympde.expr import find, EssentialBC
 
-from psydac.fem.basic   import FemField
-from psydac.fem.vector   import VectorFemField
+from psydac.fem.basic          import FemField
+from psydac.fem.vector         import VectorFemField
 from psydac.api.discretization import discretize
 
-from numpy import linspace, zeros, allclose
-import numpy as np
-from mpi4py import MPI
-import pytest
-
 #==============================================================================
 
 def run_system_1_2d_dir(Fe, Ge, f0, f1, ncells, degree):
diff --git a/psydac/api/tests/test_api_2d_vector.py b/psydac/api/tests/test_api_2d_vector.py
index 21dd23690..efc785a91 100644
--- a/psydac/api/tests/test_api_2d_vector.py
+++ b/psydac/api/tests/test_api_2d_vector.py
@@ -1,28 +1,18 @@
 # -*- coding: UTF-8 -*-
 
 from sympy import Tuple, Matrix
-from sympy import pi, cos, sin
-
-from sympde.core import Constant
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
-from sympde.calculus import laplace, hessian
-from sympde.topology import (dx, dy, dz)
-from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
-from sympde.topology import ProductSpace
-from sympde.topology import element_of
-from sympde.topology import Boundary, NormalVector, TangentVector
-from sympde.topology import Domain, Line, Square, Cube
-from sympde.topology import Trace, trace_0, trace_1
-from sympde.topology import Union
-from sympde.expr import BilinearForm, LinearForm, integral
-from sympde.expr import Norm
-from sympde.expr import find, EssentialBC
-
-from psydac.fem.vector  import VectorFemField
-from psydac.api.discretization import discretize
+from sympy import pi, sin
 
-from numpy import linspace, zeros, allclose
+from sympde.calculus import grad, dot, inner
+from sympde.topology import VectorFunctionSpace
+from sympde.topology import element_of
+from sympde.topology import Square
+from sympde.expr     import BilinearForm, LinearForm, integral
+from sympde.expr     import Norm
+from sympde.expr     import find, EssentialBC
 
+from psydac.fem.vector         import VectorFemField
+from psydac.api.discretization import discretize
 
 #==============================================================================
 def run_vector_poisson_2d_dir(solution, f, ncells, degree):
diff --git a/psydac/api/tests/test_api_3d_scalar.py b/psydac/api/tests/test_api_3d_scalar.py
index 2f4a629f9..b4f017a98 100644
--- a/psydac/api/tests/test_api_3d_scalar.py
+++ b/psydac/api/tests/test_api_3d_scalar.py
@@ -1,29 +1,22 @@
 # -*- coding: UTF-8 -*-
 
+from mpi4py import MPI
 from sympy import pi, cos, sin
+import pytest
 
-from sympde.core import Constant
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
-from sympde.calculus import laplace, hessian
-from sympde.topology import (dx, dy, dz)
-from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
-from sympde.topology import ProductSpace
+from sympde.calculus import grad, dot
+from sympde.topology import ScalarFunctionSpace
 from sympde.topology import element_of
-from sympde.topology import Boundary, NormalVector, TangentVector
-from sympde.topology import Domain, Line, Square, Cube
-from sympde.topology import Trace, trace_0, trace_1
+from sympde.topology import NormalVector
+from sympde.topology import Cube
 from sympde.topology import Union
-from sympde.expr import BilinearForm, LinearForm, integral
-from sympde.expr import Norm
-from sympde.expr import find, EssentialBC
+from sympde.expr     import BilinearForm, LinearForm, integral
+from sympde.expr     import Norm
+from sympde.expr     import find, EssentialBC
 
-from psydac.fem.basic   import FemField
+from psydac.fem.basic          import FemField
 from psydac.api.discretization import discretize
 
-from numpy import linspace, zeros, allclose
-from mpi4py import MPI
-import pytest
-
 #==============================================================================
 def run_poisson_3d_dir(solution, f, ncells, degree, comm=None):
 
diff --git a/psydac/api/tests/test_api_3d_scalar_mapping.py b/psydac/api/tests/test_api_3d_scalar_mapping.py
index 48a780c1a..5ea919656 100644
--- a/psydac/api/tests/test_api_3d_scalar_mapping.py
+++ b/psydac/api/tests/test_api_3d_scalar_mapping.py
@@ -1,40 +1,24 @@
 # -*- coding: UTF-8 -*-
 
+from mpi4py import MPI
 from sympy import pi, cos, sin
-from sympy import S
-from sympy import Tuple
-from sympy import Matrix
-
-from sympde.core import Constant
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
-from sympde.calculus import laplace, hessian
-from sympde.topology import (dx, dy, dz)
+import pytest
+import os
+
+from sympde.calculus import grad, dot
 from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
 from sympde.topology import ProductSpace
 from sympde.topology import element_of
-from sympde.topology import Unknown
-from sympde.topology import InteriorDomain, Union
-from sympde.topology import Boundary, NormalVector, TangentVector
-from sympde.topology import Domain, Line, Square, Cube
-from sympde.topology import Trace, trace_0, trace_1
+from sympde.topology import NormalVector
 from sympde.topology import Union
-from sympde.topology import Mapping
 from sympde.topology import Domain
-from sympde.expr import BilinearForm, LinearForm, integral
-from sympde.expr import Norm
-from sympde.expr import find, EssentialBC
+from sympde.expr     import BilinearForm, LinearForm, integral
+from sympde.expr     import Norm
+from sympde.expr     import find, EssentialBC
 
-from psydac.fem.basic   import FemField
+from psydac.fem.basic          import FemField
 from psydac.api.discretization import discretize
 
-from psydac.mapping.discrete import SplineMapping
-
-from numpy import linspace, zeros, allclose
-from mpi4py import MPI
-import pytest
-
-import os
-
 # ... get the mesh directory
 try:
     mesh_dir = os.environ['PSYDAC_MESH_DIR']
diff --git a/psydac/api/tests/test_api_3d_vector.py b/psydac/api/tests/test_api_3d_vector.py
index 878f511a1..aa14bc10e 100644
--- a/psydac/api/tests/test_api_3d_vector.py
+++ b/psydac/api/tests/test_api_3d_vector.py
@@ -1,28 +1,19 @@
 # -*- coding: UTF-8 -*-
 
 from sympy import Tuple, Matrix
-from sympy import pi, cos, sin
+from sympy import pi, sin
 
-from sympde.core import Constant
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
-from sympde.calculus import laplace, hessian
-from sympde.topology import (dx, dy, dz)
-from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
+from sympde.calculus import grad, dot, inner
+from sympde.topology import VectorFunctionSpace
 from sympde.topology import ProductSpace
 from sympde.topology import element_of
-from sympde.topology import Boundary, NormalVector, TangentVector
-from sympde.topology import Domain, Line, Square, Cube
-from sympde.topology import Trace, trace_0, trace_1
-from sympde.topology import Union
-from sympde.expr import BilinearForm, LinearForm, integral
-from sympde.expr import Norm
-from sympde.expr import find, EssentialBC
-
-from psydac.fem.vector  import VectorFemField
-from psydac.api.discretization import discretize
-
-from numpy import linspace, zeros, allclose
+from sympde.topology import Cube
+from sympde.expr     import BilinearForm, LinearForm, integral
+from sympde.expr     import Norm
+from sympde.expr     import find, EssentialBC
 
+from psydac.fem.vector         import VectorFemField
+from psydac.api.discretization import discretize
 
 #==============================================================================
 def run_vector_poisson_3d_dir(solution, f, ncells, degree):
diff --git a/psydac/api/tests/test_api_3d_vector_mapping.py b/psydac/api/tests/test_api_3d_vector_mapping.py
index eeda248fb..4880e985c 100644
--- a/psydac/api/tests/test_api_3d_vector_mapping.py
+++ b/psydac/api/tests/test_api_3d_vector_mapping.py
@@ -1,30 +1,20 @@
 # -*- coding: UTF-8 -*-
 
+import os
 from sympy import Tuple, Matrix
-from sympy import pi, cos, sin
-
-from sympde.core import Constant
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
-from sympde.calculus import laplace, hessian
-from sympde.topology import (dx, dy, dz)
-from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
-from sympde.topology import ProductSpace
+from sympy import pi, sin
+
+from sympde.calculus import grad, dot, inner
+from sympde.topology import VectorFunctionSpace
 from sympde.topology import element_of
-from sympde.topology import Boundary, NormalVector, TangentVector
-from sympde.topology import Domain, Line, Square, Cube
-from sympde.topology import Trace, trace_0, trace_1
-from sympde.topology import Union
+from sympde.topology import Domain
 from sympde.expr import BilinearForm, LinearForm, integral
 from sympde.expr import Norm
 from sympde.expr import find, EssentialBC
 
-from psydac.fem.vector  import VectorFemField
+from psydac.fem.vector         import VectorFemField
 from psydac.api.discretization import discretize
 
-from numpy import linspace, zeros, allclose
-
-import os
-
 # ... get the mesh directory
 try:
     mesh_dir = os.environ['PSYDAC_MESH_DIR']
diff --git a/psydac/api/tests/test_api_expr_2d_scalar.py b/psydac/api/tests/test_api_expr_2d_scalar.py
index b88da5720..48f6180b9 100644
--- a/psydac/api/tests/test_api_expr_2d_scalar.py
+++ b/psydac/api/tests/test_api_expr_2d_scalar.py
@@ -1,36 +1,17 @@
 # -*- coding: UTF-8 -*-
 
-from sympy import pi, cos, sin
-from sympy.utilities.lambdify import implemented_function
+from mpi4py import MPI
+from sympy import pi
+import numpy as np
 
-from sympde.core import Constant
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
-from sympde.calculus import laplace, hessian
-from sympde.topology import (dx, dy, dz)
-from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
-from sympde.topology import ProductSpace
+from sympde.calculus import dot, div
+from sympde.topology import VectorFunctionSpace
 from sympde.topology import element_of
+from sympde.topology import Square
 
-from sympde.topology import Boundary, NormalVector, TangentVector
-from sympde.topology import Domain, Line, Square, Cube
-from sympde.topology import Trace, trace_0, trace_1
-from sympde.topology import Union
-from sympde.expr import BilinearForm, LinearForm, integral
-from sympde.expr import Norm
-from sympde.expr import find, EssentialBC
-
-from gelato.expr import GltExpr
-
-from psydac.fem.basic   import FemField
-from psydac.fem.vector  import VectorFemField
+from psydac.fem.vector         import VectorFemField
 from psydac.api.discretization import discretize
 
-import numpy as np
-from scipy.linalg import eig as eig_solver
-from mpi4py import MPI
-import pytest
-
-
 #==============================================================================
 def run_poisson_2d_dir(ncells, degree, comm=None):
 
diff --git a/psydac/api/tests/test_api_feec_2d.py b/psydac/api/tests/test_api_feec_2d.py
index c426b0568..1f5ecf963 100644
--- a/psydac/api/tests/test_api_feec_2d.py
+++ b/psydac/api/tests/test_api_feec_2d.py
@@ -1,43 +1,23 @@
 # -*- coding: UTF-8 -*-
 
-from sympy import Tuple, Matrix
-
-from sympde.core import Constant
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
-from sympde.calculus import laplace, hessian
-from sympde.topology import (dx, dy, dz)
-from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace, Derham
-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 Boundary, NormalVector, TangentVector
-from sympde.topology import Domain, Line, Square, Cube
-from sympde.topology import Trace, trace_0, trace_1
-from sympde.topology import Union
-from sympde.expr import BilinearForm, LinearForm, integral
-from sympde.expr import Norm, TerminalExpr
-from sympde.expr import find, EssentialBC
-
-
-from psydac.fem.basic   import FemField
-from psydac.fem.vector  import VectorFemField
-from psydac.api.discretization import discretize
-from psydac.fem.vector         import ProductFemSpace
-from psydac.feec.utilities     import Interpolation, interpolation_matrices
-from psydac.feec.derivatives         import Grad, Curl, Div
-from numpy import linspace, zeros, allclose
-import numpy as np
 from mpi4py import MPI
-import pytest
-
-from scipy.sparse.linalg import cg, gmres
 from scipy.sparse.linalg import spsolve
-from scipy import linalg
+import numpy as np
+#import matplotlib.pyplot as plt
+#from mpl_toolkits.mplot3d import Axes3D
+
+from sympde.calculus import dot, div
+from sympde.topology import Derham
+from sympde.topology import ProductSpace
+from sympde.topology import element_of, elements_of
+from sympde.topology import Square
+from sympde.expr     import BilinearForm, LinearForm, integral
+from sympde.expr     import Norm
+from sympde.expr     import find
+
+from psydac.fem.basic          import FemField
+from psydac.api.discretization import discretize
 
-import matplotlib.pyplot as plt
-from mpl_toolkits.mplot3d import Axes3D
-from matplotlib import animation
 #==============================================================================
 
 def run_system_1_2d_dir(f0, sol, ncells, degree):
@@ -50,10 +30,10 @@ def run_system_1_2d_dir(f0, sol, ncells, degree):
     L2     = derham.V2  
     X      = ProductSpace(Hdiv, L2)
 
-    F = ScalarField(L2, name='F')
+    p, q = elements_of(Hdiv, names='p, q')
+    u, v = elements_of(  L2, names='u, v')
 
-    p,q = [VectorTestFunction(Hdiv, name=i) for i in ['p', 'q']]
-    u,v = [ScalarTestFunction(L2, name=i) for i in ['u', 'v']]
+    F = element_of(L2, name='F')
 
     int_0 = lambda expr: integral(domain , expr)
 
diff --git a/psydac/api/tests/test_api_glt_2d_scalar.py b/psydac/api/tests/test_api_glt_2d_scalar.py
index 1a6d5caf5..668f7a737 100644
--- a/psydac/api/tests/test_api_glt_2d_scalar.py
+++ b/psydac/api/tests/test_api_glt_2d_scalar.py
@@ -72,7 +72,7 @@ def run_poisson_2d_dir(ncells, degree, comm=None):
     # ...
 
     # ... use eigenvalue solver
-    M = ah.assemble().tosparse().todense()
+    M = ah.assemble().toarray()
     w, v = eig_solver(M)
     eig = w.real
     eig.sort()
@@ -193,7 +193,7 @@ def run_variable_coeff_2d_dir(ncells, degree, comm=None):
     # ...
 
     # ... use eigenvalue solver
-    M = ah.assemble(c=0.2).tosparse().todense()
+    M = ah.assemble(c=0.2).toarray()
     w, v = eig_solver(M)
     eig = w.real
     eig.sort()
diff --git a/psydac/api/tests/test_api_glt_2d_scalar_mapping.py b/psydac/api/tests/test_api_glt_2d_scalar_mapping.py
index 8c071f5ca..f66c6a2b3 100644
--- a/psydac/api/tests/test_api_glt_2d_scalar_mapping.py
+++ b/psydac/api/tests/test_api_glt_2d_scalar_mapping.py
@@ -1,43 +1,19 @@
 # -*- coding: UTF-8 -*-
 
-from sympy import pi, cos, sin
-from sympy import S
-from sympy import Tuple
-from sympy import Matrix
-
-from sympde.core import Constant
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
-from sympde.calculus import laplace, hessian
-from sympde.topology import (dx, dy, dz)
-from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
-from sympde.topology import ProductSpace
+import os
+import numpy as np
+from scipy.linalg import eig as eig_solver
+
+from sympde.calculus import grad, dot
+from sympde.topology import ScalarFunctionSpace
 from sympde.topology import element_of
-from sympde.topology import Unknown
-from sympde.topology import InteriorDomain, Union
-from sympde.topology import Boundary, NormalVector, TangentVector
-from sympde.topology import Domain, Line, Square, Cube
-from sympde.topology import Trace, trace_0, trace_1
-from sympde.topology import Union
-from sympde.topology import Mapping
-from sympde.expr import BilinearForm, LinearForm, integral
-from sympde.expr import Norm
-from sympde.expr import find, EssentialBC
+from sympde.topology import Domain
+from sympde.expr     import BilinearForm, integral
 
 from gelato.expr import GltExpr
 
-from psydac.fem.basic   import FemField
-from psydac.fem.vector  import ProductFemSpace, VectorFemField
 from psydac.api.discretization import discretize
 
-from psydac.mapping.discrete import SplineMapping
-
-import numpy as np
-from scipy.linalg import eig as eig_solver
-from mpi4py import MPI
-import pytest
-
-import os
-
 # ... get the mesh directory
 try:
     mesh_dir = os.environ['PSYDAC_MESH_DIR']
@@ -91,7 +67,7 @@ def run_poisson_2d_dir(filename, comm=None):
     # ...
 
     # ... use eigenvalue solver
-    M = ah.assemble().tosparse().todense()
+    M = ah.assemble().toarray()
     w, v = eig_solver(M)
     eig = w.real
     eig.sort()
diff --git a/psydac/api/tests/test_api_glt_2d_vector.py b/psydac/api/tests/test_api_glt_2d_vector.py
index 648527797..a11991fbd 100644
--- a/psydac/api/tests/test_api_glt_2d_vector.py
+++ b/psydac/api/tests/test_api_glt_2d_vector.py
@@ -1,33 +1,18 @@
 # -*- coding: UTF-8 -*-
 
-from sympy import Tuple, Matrix
-from sympy import pi, cos, sin
-
-from sympde.core import Constant
-from sympde.calculus import grad, dot, inner, cross, rot, curl, div
-from sympde.calculus import laplace, hessian
-from sympde.topology import (dx, dy, dz)
-from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
-from sympde.topology import ProductSpace
+import numpy as np
+from scipy.linalg import eig as eig_solver
+
+from sympde.calculus import grad, inner
+from sympde.topology import VectorFunctionSpace
 from sympde.topology import element_of
-from sympde.topology import Boundary, NormalVector, TangentVector
-from sympde.topology import Domain, Line, Square, Cube
-from sympde.topology import Trace, trace_0, trace_1
-from sympde.topology import Union
-from sympde.expr import BilinearForm, LinearForm, integral
-from sympde.expr import Norm
-from sympde.expr import find, EssentialBC
+from sympde.topology import Square
+from sympde.expr     import BilinearForm, integral
 
 from gelato.expr import GltExpr
 
-from psydac.fem.vector  import VectorFemField
 from psydac.api.discretization import discretize
 
-import numpy as np
-from scipy.linalg import eig as eig_solver
-from mpi4py import MPI
-import pytest
-
 #==============================================================================
 def run_vector_poisson_2d_dir(ncells, degree):
 
@@ -73,7 +58,7 @@ def run_vector_poisson_2d_dir(ncells, degree):
     # ...
 
     # ... use eigenvalue solver
-    M = ah.assemble().tosparse().todense()
+    M = ah.assemble().toarray()
     w, v = eig_solver(M)
     eig = w.real
     eig.sort()
diff --git a/psydac/api/tests/utils.py b/psydac/api/tests/utils.py
index c6f7ece41..9c644b856 100644
--- a/psydac/api/tests/utils.py
+++ b/psydac/api/tests/utils.py
@@ -20,7 +20,7 @@ def assert_identical_coo(A, B):
         assert(allclose(A.col,  B.col))
         assert(allclose(A.data, B.data))
 
-#        assert(allclose(A.todense(), B.todense()))
+#        assert(allclose(A.toarray(), B.toarray()))
 
     else:
         raise TypeError('Wrong types for entries')
diff --git a/psydac/cad/geometry.py b/psydac/cad/geometry.py
index dbdf4f10b..bbccfab39 100644
--- a/psydac/cad/geometry.py
+++ b/psydac/cad/geometry.py
@@ -161,7 +161,7 @@ def read( self, filename, comm=MPI.COMM_WORLD ):
             kwargs = {}
 
         h5  = h5py.File( filename, mode='r', **kwargs )
-        yml = yaml.load( h5['geometry.yml'][()] )
+        yml = yaml.load( h5['geometry.yml'][()], Loader=yaml.SafeLoader )
 
         ldim = yml['ldim']
         pdim = yml['pdim']
diff --git a/psydac/cad/tests/test_geometry.py b/psydac/cad/tests/test_geometry.py
index 40690cd8e..1ac0b17c9 100644
--- a/psydac/cad/tests/test_geometry.py
+++ b/psydac/cad/tests/test_geometry.py
@@ -154,9 +154,24 @@ def test_geometry_1():
 #==============================================================================
 
 def teardown_module():
+    import os
     from sympy import cache
     cache.clear_cache()
 
+    # Remove HDF5 files generated by Geometry.export()
+    filenames = [
+        'geo.h5',
+        'geo_0.h5',
+        'geo_1.h5',
+        'quart_circle.h5',
+        'quart_circle_0.h5',
+        'quart_circle_1.h5',
+        'circle.h5'
+    ]
+    for fname in filenames:
+        if os.path.exists(fname):
+            os.remove(fname)
+
 def teardown_function():
     from sympy import cache
     cache.clear_cache()
diff --git a/psydac/feec/utilities.py b/psydac/feec/utilities.py
index b6c8706fd..d4916ee60 100644
--- a/psydac/feec/utilities.py
+++ b/psydac/feec/utilities.py
@@ -33,9 +33,9 @@ def solve(self, rhs, out=None):
        return out
        
     def toarray(self):
-        M1 = self.args[0].todense()
+        M1 = self.args[0].toarray()
         for M2 in self.args[1:]:
-            M1 = np.kron(M1, M2.todense())
+            M1 = np.kron(M1, M2.toarray())
         return M1
         
 
diff --git a/psydac/linalg/tests/test_stencil_matrix.py b/psydac/linalg/tests/test_stencil_matrix.py
index 384afa04b..86c9825bf 100644
--- a/psydac/linalg/tests/test_stencil_matrix.py
+++ b/psydac/linalg/tests/test_stencil_matrix.py
@@ -621,7 +621,7 @@ def test_stencil_matrix_2d_parallel_transpose( n1, n2, p1, p2, P1, P2, reorder )
     #...
     Ts_exact = Ts_exact.tocsr()
     for i, j in zip(*Ts_exact.nonzero()):
-        i1, i2 = np.unravel_index( i, dims=[n1, n2], order='C' )
+        i1, i2 = np.unravel_index( i, shape=[n1, n2], order='C' )
         if not (s1 <= i1 <= e1 and s2 <= i2 <= e2):
             Ts_exact[i, j] = 0.0
     Ts_exact = Ts_exact.tocoo()
diff --git a/psydac/version.py b/psydac/version.py
new file mode 100644
index 000000000..a4e2017f0
--- /dev/null
+++ b/psydac/version.py
@@ -0,0 +1 @@
+__version__ = "0.1"
diff --git a/requirements.txt b/requirements.txt
index 94b6fe43a..1543c677a 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,14 +1,3 @@
-# Packages from PyPi
-numpy>=1.13, <1.16
-scipy>=0.18
-matplotlib
-mpi4py
-h5py
-sympy>=1.3
-pytest
-pyyaml
-yamlloader
-
 # h5py must be built from source using MPI compiler
 # and linked to parallel HDF5 library. To do so set
 #
@@ -16,10 +5,10 @@ yamlloader
 # HDF5_MPI="ON"
 # HDF5_DIR=/usr/lib/x86_64-linux-gnu/hdf5/openmpi
 #
+numpy>=1.16
+mpi4py
+h5py
 --no-binary h5py
 
-# Igakit
-https://bitbucket.org/dalcinl/igakit/get/default.tar.gz
-
-# In addition, PSYDAC depends on "sympde" and "pyccel" libraries,
-# to be installed manually from their git repositories.
+# Igakit (TODO: remove library dependency)
+https://bitbucket.org/dalcinl/igakit/get/master.tar.gz
diff --git a/run_tests.sh b/run_tests.sh
deleted file mode 100755
index 2777d3c28..000000000
--- a/run_tests.sh
+++ /dev/null
@@ -1,4 +0,0 @@
-#!/bin/bash
-
-python3 -m pytest --pyargs spl/api/tests/ -m "not parallel"
-mpirun -np 2 python3 -m pytest --pyargs spl/api/tests/ -m "parallel"
diff --git a/setup.cfg b/setup.cfg
deleted file mode 100644
index fc404b260..000000000
--- a/setup.cfg
+++ /dev/null
@@ -1,17 +0,0 @@
-[config_fc]
-fcompiler=gnu95
-#f90flags=-g3 -O0 -fcheck=all -Wall -Wextra
-#noopt=True
-#noarch=True
-
-[build]
-#debug = 0
-
-[install]
-optimize = 1
-
-[sdist]
-force_manifest = 1
-
-[nosetests]
-where = test
diff --git a/setup.py b/setup.py
index b97b27f58..701b70881 100644
--- a/setup.py
+++ b/setup.py
@@ -1,15 +1,22 @@
 # -*- coding: UTF-8 -*-
 #! /usr/bin/python
 
-import sys
-import os
+from pathlib    import Path
+from importlib  import util
 from setuptools import find_packages
 from numpy.distutils.core import setup
 from numpy.distutils.core import Extension
-import psydac
+
+# ...
+# Load module 'psydac.version' without running file 'psydac.__init__'
+path = Path(__file__).parent / 'psydac' / 'version.py'
+spec = util.spec_from_file_location('version', str(path))
+mod  = util.module_from_spec(spec)
+spec.loader.exec_module(mod)
+# ...
 
 NAME    = 'psydac'
-VERSION = psydac.__version__
+VERSION = mod.__version__
 AUTHOR  = 'Ahmed Ratnani, Jalal Lakhlili, Yaman Güçlü'
 EMAIL   = 'ratnaniahmed@gmail.com'
 URL     = 'http://www.ahmed.ratnani.org'
@@ -36,12 +43,36 @@
 # ...
 
 # ...
-install_requires = ['numpy', 'scipy']
+install_requires = [
+
+    # Third-party packages from PyPi
+    'numpy>=1.16',
+    'scipy>=0.18',
+    'sympy>=1.2',
+    'matplotlib',
+    'pytest',
+    'pyyaml',
+    'yamlloader',
+
+    # Our packages from PyPi
+    'sympde',
+    'pyccel',
+    'gelato',
+
+    # In addition, we depend on mpi4py and h5py (MPI version).
+    # Since h5py must be built from source, we run the commands
+    #
+    # python3 -m pip install requirements.txt
+    # python3 -m pip install .
+    'mpi4py',
+    'h5py',
+]
+
 dependency_links = []
 # ...
 
 
-# ... bspline extension
+# ... bspline extension (TODO: remove Fortran files from library)
 bsp_ext = Extension(name    = 'psydac.core.bsp',
                         sources = ['psydac/core/external/bspline.F90',
                                    'psydac/core/external/pppack.F90',