diff --git a/Makefile b/Makefile
index 296a6b19..c9c2805c 100644
--- a/Makefile
+++ b/Makefile
@@ -177,7 +177,7 @@ flake8:
 	rm flake8.txt
 
 cython-lint:
-	- cython-lint --max-line-length=160 drivers examples packageTools base metisCy fem multilevelSolver nl tests > cython-lint.txt 2>&1
+	- cython-lint --max-line-length=160 --ignore="E741" drivers examples packageTools base metisCy fem multilevelSolver nl tests > cython-lint.txt 2>&1
 	flake8_junit cython-lint.txt cython-lint.xml
 	rm cython-lint.txt
 	sed 's/name="flake8"/name="cython-lint"/g' cython-lint.xml > cython-lint2.xml
diff --git a/base/PyNucleus_base/factory.py b/base/PyNucleus_base/factory.py
index 342ad2dd..231e5c87 100644
--- a/base/PyNucleus_base/factory.py
+++ b/base/PyNucleus_base/factory.py
@@ -51,13 +51,6 @@ def numRegistered(self, countAliases=False):
             return len(self.classes) + len(self.aliases)
 
     def __str__(self):
-        s = ''
-        for canonical_name in self.classes:
-            a = [self.aliases[canonical_alias][0] for canonical_alias in self.aliases if self.aliases[canonical_alias][1] == canonical_name]
-            s += '{} {} {}\n'.format(canonical_name, a, self.classes[canonical_name])
-        return s
-
-    def __repr__(self):
         s = 'Available:\n'
         for canonical_name in self.classes:
             name = self.classes[canonical_name][0]
@@ -77,6 +70,3 @@ def __repr__(self):
             else:
                 s += '\'{}\', signature: \'{}\'\n'.format(name, sig)
         return s
-
-    def print(self):
-        print(repr(self))
diff --git a/base/PyNucleus_base/solver_factory.py b/base/PyNucleus_base/solver_factory.py
index 9ba3fffa..5b3f5724 100644
--- a/base/PyNucleus_base/solver_factory.py
+++ b/base/PyNucleus_base/solver_factory.py
@@ -95,12 +95,12 @@ def build(self, name, **kwargs):
                 solvers[k].setPreconditioner(solvers[k+1].asPreconditioner())
             return solvers[0]
 
-    def __repr__(self):
+    def __str__(self):
         s = ''
         if self.singleLevelSolverFactory.numRegistered() > 0:
             s += 'Single level solvers:\n'
-            s += repr(self.singleLevelSolverFactory)
+            s += str(self.singleLevelSolverFactory)
         if self.multiLevelSolverFactory.numRegistered() > 0:
             s += 'Multi level solvers:\n'
-            s += repr(self.multiLevelSolverFactory)
+            s += str(self.multiLevelSolverFactory)
         return s
diff --git a/base/PyNucleus_base/timestepping.py b/base/PyNucleus_base/timestepping.py
index 6cb2b3d2..8b88c578 100644
--- a/base/PyNucleus_base/timestepping.py
+++ b/base/PyNucleus_base/timestepping.py
@@ -15,11 +15,13 @@ class Stepper:
     """
     Solve
 
-    f(t, u, u_t) = A(t, u_t) + B(t, u) - g(t) = 0
+    f(t, u, u_t) = A(t, u_t) + B(t, u) - g(t) = 0.
 
-    mass:                  t, u  -> A(t, u)
-    explicit:              t, u  -> B(t, u)
-    forcing:               t     -> g(t)
+    We assume that A is linear in its second argument. Let
+
+    mass:                  t, u           -> A(t, u)
+    explicit:              t, u           -> B(t, u)
+    forcing:               t              -> g(t)
     implicitSolverBuilder: t, alpha, beta -> solver for alpha*A(t, u)+beta*B(t, u) = g
 
     """
@@ -51,9 +53,12 @@ def apply_jacobian(self, t, dt, ut, ut_plus_dt, residual, alpha=1., beta=1.):
 
 
 class CrankNicolson(Stepper):
+    """
+    1/dt*A(t+dt, u_{k+1}) + theta*B(t+dt, u_{k+1}) = (1-theta)*g(t_{k}) + theta*g(t_{k+1}) + 1/dt*A(t_{k}, u_{k}) - (1-theta)*B(t_{k}, u_{k})
+    """
     def __init__(self, dm, mass, solverBuilder, forcing, explicit=None, theta=0.5, dt=None,
                  solverIsTimeDependent=False, explicitIslinearAndTimeIndependent=False):
-        assert theta > 0 and theta <= 1.
+        assert theta >= 0 and theta <= 1.
         self.theta = theta
         super(CrankNicolson, self).__init__(dm, mass, solverBuilder, forcing, explicit, dt, solverIsTimeDependent, explicitIslinearAndTimeIndependent)
         if self.theta < 1.:
@@ -73,12 +78,16 @@ def getSolver(self, t, dt):
             return self.solverBuilder(t+dt, 1./dt, self.theta)
 
     def setRHS(self, t, dt, rhs):
-        self.forcing(t+dt, rhs)
-        if self.theta < 1.:
+        if self.theta > 0. and self.theta < 1.:
+            self.forcing(t+dt, rhs)
             rhs *= self.theta
             self.forcing(t, self.rhs2)
             self.rhs2 *= (1.-self.theta)
             rhs += self.rhs2
+        elif self.theta == 0.:
+            self.forcing(t, rhs)
+        elif self.theta == 1.:
+            self.forcing(t+dt, rhs)
 
     def step(self, t, dt, u, forcingVector=None):
         if dt is None:
@@ -105,6 +114,9 @@ def step(self, t, dt, u, forcingVector=None):
         return t+dt
 
     def residual(self, t, dt, ut, ut_plus_dt, residual, alpha=1., beta=1., forcingVector=None):
+        # alpha/dt*[A(t+dt, ut_plus_dt) - A(t, ut)]
+        # + beta*[(1-theta)*B(t, ut) + theta*B(t+dt, ut_plus_dt)]
+        # - (1-theta)*g(t) - theta*g(t+dt)
         if abs(alpha/dt) > 0:
             self.mass(t, ut, self.rhs)
             self.mass(t+dt, ut_plus_dt, self.rhs2)
@@ -134,6 +146,8 @@ def residual(self, t, dt, ut, ut_plus_dt, residual, alpha=1., beta=1., forcingVe
         residual -= self.rhs
 
     def apply_jacobian(self, t, dt, ut, ut_plus_dt, residual, alpha=1., beta=1.):
+        # alpha/dt*[A(t+dt, ut_plus_dt) - A(t, ut)]
+        # + beta*[(1-theta)*B(t, ut) + theta*B(t+dt, ut_plus_dt)]
         if abs(alpha/dt) > 0:
             self.mass(t, ut, self.rhs)
             self.mass(t+dt, ut_plus_dt, self.rhs2)
@@ -157,30 +171,21 @@ def apply_jacobian(self, t, dt, ut, ut_plus_dt, residual, alpha=1., beta=1.):
             residual += self.rhs2
 
 
-class ExplicitEuler(Stepper):
-    def __init__(self, dm, mass, solverBuilder, forcing, explicit, solverIsTimeDependent=False):
-        assert explicit is not None
-        super(ExplicitEuler, self).__init__(dm, mass, solverBuilder, forcing, explicit, 0., solverIsTimeDependent)
-        self.rhs = self.dm.zeros()
-        self.rhs2 = self.dm.zeros()
-
-    def step(self, t, dt, u):
-        self.forcing(t+dt, self.rhs)
-        self.rhs *= dt
-        self.explicit(t, u, self.rhs2)
-        self.rhs2 *= -dt
-        self.rhs += self.rhs2
-        if not self.solverIsTimeDependent:
-            solver = self.solver
-        else:
-            solver = self.solverBuilder(t+dt, 0.)
-        if isinstance(solver, iterative_solver):
-            solver.setInitialGuess(u)
-        solver(self.rhs, u)
-        return t+dt
+class ExplicitEuler(CrankNicolson):
+    """
+    1/dt*A(t+dt, u_{k+1}) = g(t_{k}) + 1/dt*A(t_{k}, u_{k}) - B(t_{k}, u_{k})
+    """
+    def __init__(self, dm, mass, solverBuilder, forcing, explicit=None, dt=None, solverIsTimeDependent=False):
+        super(ExplicitEuler, self).__init__(dm, mass, solverBuilder, forcing, explicit,
+                                            dt=dt,
+                                            theta=0.,
+                                            solverIsTimeDependent=solverIsTimeDependent)
 
 
 class ImplicitEuler(CrankNicolson):
+    """
+    1/dt*A(t+dt, u_{k+1}) + B(t+dt, u_{k+1}) = g(t_{k+1}) + 1/dt*A(t_{k}, u_{k})
+    """
     def __init__(self, dm, mass, solverBuilder, forcing, explicit=None, dt=None, solverIsTimeDependent=False):
         super(ImplicitEuler, self).__init__(dm, mass, solverBuilder, forcing, explicit,
                                             dt=dt,
diff --git a/base/PyNucleus_base/utilsFem.py b/base/PyNucleus_base/utilsFem.py
index 61c80b67..ce7a9250 100644
--- a/base/PyNucleus_base/utilsFem.py
+++ b/base/PyNucleus_base/utilsFem.py
@@ -310,7 +310,7 @@ def saveDictToHDF5(params, f, ignore=set()):
 def loadDictFromHDF5(f):
     import h5py
     from . linear_operators import LinearOperator
-    from PyNucleus.fem.DoFMaps import DoFMap
+    from PyNucleus_fem.DoFMaps import DoFMap
     params = {}
     for key in f.attrs:
         if isinstance(f.attrs[key], h5py.Empty):
@@ -352,7 +352,7 @@ def loadDictFromHDF5(f):
                 else:
                     params[key] = LinearOperator.HDF5read(f[key])
             elif 'vertices' in f[key] and 'cells' in f[key]:
-                from PyNucleus.fem.mesh import meshNd
+                from PyNucleus_fem.mesh import meshNd
                 params[key] = meshNd.HDF5read(f[key])
             else:
                 params[key] = loadDictFromHDF5(f[key])
@@ -369,7 +369,7 @@ def loadDictFromHDF5(f):
 
 
 def processDictForYaml(params):
-    from PyNucleus.fem.functions import function
+    from PyNucleus_fem.functions import function
     paramsNew = {}
     for key in params:
         if isinstance(params[key], dict):
@@ -411,7 +411,7 @@ def updateFromDefaults(params, defaults):
             updateFromDefaults(params[key], defaults[key])
 
 
-def getMPIinfo(grp):
+def getMPIinfo(grp, verbose=False):
     from sys import modules
     if 'mpi4py.MPI' in modules:
         import mpi4py
@@ -427,13 +427,15 @@ def getMPIinfo(grp):
         if MPI.COMM_WORLD.rank == 0:
             hosts = ','.join(set(hosts))
         grp.add('MPI library', '{}'.format(MPI.Get_library_version()[:-1]))
-        for label, value in [('MPI standard supported', MPI.Get_version()),
-                             ('Vendor', MPI.get_vendor()),
-                             ('Level of thread support', t[MPI.Query_thread()]),
-                             ('Is threaded', MPI.Is_thread_main()),
-                             ('Threads requested', mpi4py.rc.threads),
-                             ('Thread level requested', mpi4py.rc.thread_level),
-                             ('Hosts', hosts),
+        if verbose:
+            for label, value in [('MPI standard supported', MPI.Get_version()),
+                                 ('Vendor', MPI.get_vendor()),
+                                 ('Level of thread support', t[MPI.Query_thread()]),
+                                 ('Is threaded', MPI.Is_thread_main()),
+                                 ('Threads requested', mpi4py.rc.threads),
+                                 ('Thread level requested', mpi4py.rc.thread_level)]:
+                grp.add(label, value)
+        for label, value in [('Hosts', hosts),
                              ('Communicator size', MPI.COMM_WORLD.size)]:
             grp.add(label, value)
 
@@ -487,8 +489,8 @@ def getSystemInfo(grp, argv=None, envVars=[('OMP_NUM_THREADS', True)]):
             versions[(version, sha)].append(pkg)
         except KeyError:
             versions[(version, sha)] = [pkg]
-    for version in versions:
-        grp.add('PyNucleus_'+(','.join(versions[version])), version)
+    for version, sha in versions:
+        grp.add('PyNucleus_'+(','.join(versions[(version, sha)])), '{}, {}'.format(version, sha))
 
 
 class MPIFileHandler(logging.Handler):
@@ -1382,6 +1384,7 @@ def runDriver(path, py, python=None, timeout=600, ranks=None, cacheDir='',
             py += ['--overwriteCache']
     else:
         py += ['--test']
+    py += ['--disableFileLog']
     if extra is not None:
         plotDir.mkdir(exist_ok=True, parents=True)
         py += ['--plotFolder={}'.format(plotDir), '--plotFormat=png']
@@ -1444,7 +1447,7 @@ def __init__(self, name, params=[]):
                 fields.append('True|False')
             else:
                 raise NotImplementedError()
-        self.regexp = re.compile(name+'\('+','.join(['\s*(' + f + ')\s*' for f in fields])+'\)')
+        self.regexp = re.compile(name+'\(?'+','.join(['\s*(' + f + ')\s*' for f in fields])+'\)?')
 
     def match(self, s):
         return self.regexp.match(s) is not None
@@ -1841,7 +1844,13 @@ def processCmdline(self, params):
         pass
 
     def getIdentifier(self, params):
-        return ''
+        from sys import argv
+        identifier = '-'.join(argv)
+        pos = identifier.rfind('/')
+        if pos >= 0:
+            return identifier[pos+1:]
+        else:
+            return identifier
 
     def process(self, params):
         self.processCmdline(params)
diff --git a/docs/PyNucleus_base.rst b/docs/PyNucleus_base.rst
index 11ee8907..2c0d1831 100644
--- a/docs/PyNucleus_base.rst
+++ b/docs/PyNucleus_base.rst
@@ -1,111 +1,111 @@
 
-PyNucleus.base package
+PyNucleus_base package
 =======================
 
 .. automodule:: PyNucleus_base
 
 
-PyNucleus.base.blas module
+PyNucleus_base.blas module
 ---------------------------
 
 .. automodule:: PyNucleus_base.blas
 
 
-PyNucleus.base.convergence module
+PyNucleus_base.convergence module
 ----------------------------------
 
 .. automodule:: PyNucleus_base.convergence
 
 
-PyNucleus.base.factory module
+PyNucleus_base.factory module
 ------------------------------
 
 .. automodule:: PyNucleus_base.factory
 
 
-PyNucleus.base.intTuple module
+PyNucleus_base.intTuple module
 -------------------------------
 
 .. automodule:: PyNucleus_base.intTuple
 
 
-PyNucleus.base.ip.norm module
+PyNucleus_base.ip.norm module
 -------------------------------
 
 .. automodule:: PyNucleus_base.ip_norm
 
 
-PyNucleus.base.linalg module
+PyNucleus_base.linalg module
 -----------------------------
 
 .. automodule:: PyNucleus_base.linalg
 
 
-PyNucleus.base.linear.operators module
+PyNucleus_base.linear.operators module
 ----------------------------------------
 
 .. automodule:: PyNucleus_base.linear_operators
 
 
-PyNucleus.base.memProfile module
+PyNucleus_base.memProfile module
 ---------------------------------
 
 .. automodule:: PyNucleus_base.memProfile
 
 
-PyNucleus.base.myTypes module
+PyNucleus_base.myTypes module
 ------------------------------
 
 .. automodule:: PyNucleus_base.myTypes
 
 
-PyNucleus.base.performanceLogger module
+PyNucleus_base.performanceLogger module
 ----------------------------------------
 
 .. automodule:: PyNucleus_base.performanceLogger
 
 
-PyNucleus.base.plot.utils module
+PyNucleus_base.plot.utils module
 ----------------------------------
 
 .. automodule:: PyNucleus_base.plot_utils
 
 
-PyNucleus.base.solver.factory module
+PyNucleus_base.solver.factory module
 --------------------------------------
 
 .. automodule:: PyNucleus_base.solver_factory
 
 
-PyNucleus.base.solvers module
+PyNucleus_base.solvers module
 ------------------------------
 
 .. automodule:: PyNucleus_base.solvers
 
-PyNucleus.base.sparseGraph module
+PyNucleus_base.sparseGraph module
 ----------------------------------
 
 .. automodule:: PyNucleus_base.sparseGraph
 
-PyNucleus.base.sparsityPattern module
+PyNucleus_base.sparsityPattern module
 --------------------------------------
 
 .. automodule:: PyNucleus_base.sparsityPattern
 
 
-PyNucleus.base.tupleDict module
+PyNucleus_base.tupleDict module
 --------------------------------
 
 .. automodule:: PyNucleus_base.tupleDict
 
 
-PyNucleus.base.utilsCy module
+PyNucleus_base.utilsCy module
 ------------------------------
 
 .. automodule:: PyNucleus_base.utilsCy
 
 
-PyNucleus.base.utilsFem module
+PyNucleus_base.utilsFem module
 -------------------------------
 
 .. automodule:: PyNucleus_base.utilsFem
diff --git a/docs/PyNucleus_fem.rst b/docs/PyNucleus_fem.rst
index a6d86c7e..093ccc45 100644
--- a/docs/PyNucleus_fem.rst
+++ b/docs/PyNucleus_fem.rst
@@ -1,101 +1,101 @@
 
-PyNucleus.fem package
+PyNucleus_fem package
 ======================
 
 .. automodule:: PyNucleus_fem
 
 
-PyNucleus.fem.DoFMaps module
+PyNucleus_fem.DoFMaps module
 -----------------------------
 
 .. automodule:: PyNucleus_fem.DoFMaps
 
 
-PyNucleus.fem.algebraicOverlaps module
+PyNucleus_fem.algebraicOverlaps module
 ---------------------------------------
 
 .. automodule:: PyNucleus_fem.algebraicOverlaps
 
 
-PyNucleus.fem.boundaryLayerCy module
+PyNucleus_fem.boundaryLayerCy module
 -------------------------------------
 
 .. automodule:: PyNucleus_fem.boundaryLayerCy
 
 
-PyNucleus.fem.distributed.operators module
+PyNucleus_fem.distributed.operators module
 --------------------------------------------
 
 .. automodule:: PyNucleus_fem.distributed_operators
 
 
-PyNucleus.fem.femCy module
+PyNucleus_fem.femCy module
 ---------------------------
 
 .. automodule:: PyNucleus_fem.femCy
 
 
-PyNucleus.fem.functions module
+PyNucleus_fem.functions module
 -------------------------------
 
 .. automodule:: PyNucleus_fem.functions
 
 
-PyNucleus.fem.mesh module
+PyNucleus_fem.mesh module
 --------------------------
 
 .. automodule:: PyNucleus_fem.mesh
 
 
-PyNucleus.fem.meshConstruction module
+PyNucleus_fem.meshConstruction module
 --------------------------------------
 
 .. automodule:: PyNucleus_fem.meshConstruction
 
 
-PyNucleus.fem.meshCy module
+PyNucleus_fem.meshCy module
 ----------------------------
 
 .. automodule:: PyNucleus_fem.meshCy
 
 
-PyNucleus.fem.meshOverlaps module
+PyNucleus_fem.meshOverlaps module
 ----------------------------------
 
 .. automodule:: PyNucleus_fem.meshOverlaps
 
 
-PyNucleus.fem.meshPartitioning module
+PyNucleus_fem.meshPartitioning module
 --------------------------------------
 
 .. automodule:: PyNucleus_fem.meshPartitioning
 
 
-PyNucleus.fem.pdeProblems module
+PyNucleus_fem.pdeProblems module
 ---------------------------------
 
 .. automodule:: PyNucleus_fem.pdeProblems
 
 
-PyNucleus.fem.quadrature module
+PyNucleus_fem.quadrature module
 --------------------------------
 
 .. automodule:: PyNucleus_fem.quadrature
 
 
-PyNucleus.fem.repartitioner module
+PyNucleus_fem.repartitioner module
 -----------------------------------
 
 .. automodule:: PyNucleus_fem.repartitioner
 
 
-PyNucleus.fem.simplexMapper module
+PyNucleus_fem.simplexMapper module
 -----------------------------------
 
 .. automodule:: PyNucleus_fem.simplexMapper
 
 
-PyNucleus.fem.splitting module
+PyNucleus_fem.splitting module
 -------------------------------
 
 .. automodule:: PyNucleus_fem.splitting
diff --git a/docs/PyNucleus_metisCy.rst b/docs/PyNucleus_metisCy.rst
index 06a7afe2..ada11c90 100644
--- a/docs/PyNucleus_metisCy.rst
+++ b/docs/PyNucleus_metisCy.rst
@@ -1,17 +1,17 @@
 
-PyNucleus.metisCy package
+PyNucleus_metisCy package
 ==========================
 
 .. automodule:: PyNucleus_metisCy
 
 
-PyNucleus.metisCy.metisCy module
+PyNucleus_metisCy.metisCy module
 ---------------------------------
 
 .. automodule:: PyNucleus_metisCy.metisCy
 
 
-PyNucleus.metisCy.parmetisCy module
+PyNucleus_metisCy.parmetisCy module
 ------------------------------------
 
 .. automodule:: PyNucleus_metisCy.parmetisCy
diff --git a/docs/PyNucleus_multilevelSolver.rst b/docs/PyNucleus_multilevelSolver.rst
index c9beb982..64973cdc 100644
--- a/docs/PyNucleus_multilevelSolver.rst
+++ b/docs/PyNucleus_multilevelSolver.rst
@@ -1,53 +1,53 @@
 
-PyNucleus.multilevelSolver package
+PyNucleus_multilevelSolver package
 ===================================
 
 .. automodule:: PyNucleus_multilevelSolver
 
 
-PyNucleus.multilevelSolver.coarseSolvers module
+PyNucleus_multilevelSolver.coarseSolvers module
 ------------------------------------------------
 
 .. automodule:: PyNucleus_multilevelSolver.coarseSolvers
 
 
-PyNucleus.multilevelSolver.connectors module
+PyNucleus_multilevelSolver.connectors module
 ---------------------------------------------
 
 .. automodule:: PyNucleus_multilevelSolver.connectors
 
 
-PyNucleus.multilevelSolver.geometricMG module
+PyNucleus_multilevelSolver.geometricMG module
 ----------------------------------------------
 
 .. automodule:: PyNucleus_multilevelSolver.geometricMG
 
 
-PyNucleus.multilevelSolver.hierarchies module
+PyNucleus_multilevelSolver.hierarchies module
 ----------------------------------------------
 
 .. automodule:: PyNucleus_multilevelSolver.hierarchies
 
 
-PyNucleus.multilevelSolver.levels module
+PyNucleus_multilevelSolver.levels module
 -----------------------------------------
 
 .. automodule:: PyNucleus_multilevelSolver.levels
 
 
-PyNucleus.multilevelSolver.multigrid module
+PyNucleus_multilevelSolver.multigrid module
 --------------------------------------------
 
 .. automodule:: PyNucleus_multilevelSolver.multigrid
 
 
-PyNucleus.multilevelSolver.restrictionProlongation module
+PyNucleus_multilevelSolver.restrictionProlongation module
 ----------------------------------------------------------
 
 .. automodule:: PyNucleus_multilevelSolver.restrictionProlongation
 
 
-PyNucleus.multilevelSolver.smoothers module
+PyNucleus_multilevelSolver.smoothers module
 --------------------------------------------
 
 .. automodule:: PyNucleus_multilevelSolver.smoothers
diff --git a/docs/PyNucleus_nl.rst b/docs/PyNucleus_nl.rst
index 0140328b..13504b89 100644
--- a/docs/PyNucleus_nl.rst
+++ b/docs/PyNucleus_nl.rst
@@ -1,83 +1,83 @@
 
-PyNucleus.nl package
+PyNucleus_nl package
 =====================
 
 Submodules
 ----------
 
 
-PyNucleus.nl.clusterMethodCy module
+PyNucleus_nl.clusterMethodCy module
 ------------------------------------
 
 .. automodule:: PyNucleus_nl.clusterMethodCy
 
 
-PyNucleus.nl.config module
+PyNucleus_nl.config module
 ---------------------------
 
 .. automodule:: PyNucleus_nl.config
 
 
-PyNucleus.nl.fractionalLaplacian1D module
+PyNucleus_nl.fractionalLaplacian1D module
 ------------------------------------------
 
 .. automodule:: PyNucleus_nl.fractionalLaplacian1D
 
 
-PyNucleus.nl.fractionalLaplacian2D module
+PyNucleus_nl.fractionalLaplacian2D module
 ------------------------------------------
 
 .. automodule:: PyNucleus_nl.fractionalLaplacian2D
 
 
-PyNucleus.nl.twoPointFunctions module
+PyNucleus_nl.twoPointFunctions module
 --------------------------------------
 
 .. automodule:: PyNucleus_nl.twoPointFunctions
 
-PyNucleus.nl.fractionalOrders module
+PyNucleus_nl.fractionalOrders module
 -------------------------------------
 
 .. automodule:: PyNucleus_nl.fractionalOrders
 
 
-PyNucleus.nl.kernelNormalization module
+PyNucleus_nl.kernelNormalization module
 ----------------------------------------
 
 .. automodule:: PyNucleus_nl.kernelNormalization
 
 
-PyNucleus.nl.interactionDomains module
+PyNucleus_nl.interactionDomains module
 ---------------------------------------
 
 .. automodule:: PyNucleus_nl.interactionDomains
 
 
-PyNucleus.nl.kernels module
+PyNucleus_nl.kernels module
 ----------------------------
 
 .. automodule:: PyNucleus_nl.kernels
 
 
-PyNucleus.nl.kernelsCy module
+PyNucleus_nl.kernelsCy module
 -----------------------------
 
 .. automodule:: PyNucleus_nl.kernelsCy
 
 
-PyNucleus.nl.nonlocalOperator module
+PyNucleus_nl.nonlocalOperator module
 ------------------------------------------
 
 .. automodule:: PyNucleus_nl.nonlocalOperator
 
 
-PyNucleus.nl.nonlocalBuilder module
+PyNucleus_nl.nonlocalBuilder module
 --------------------------------------
 
 .. automodule:: PyNucleus_nl.nonlocalBuilder
 
 
-PyNucleus.nl.nonlocalProblems module
+PyNucleus_nl.nonlocalProblems module
 -------------------------------------
 
 .. automodule:: PyNucleus_nl.nonlocalProblems
diff --git a/docs/conf.py b/docs/conf.py
index dcb9aba6..b4de4320 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -98,4 +98,4 @@ def __call__(self, filename):
 # html_static_path = ['_static']
 
 
-# autodoc_mock_imports = ["PyNucleus.base"]
+# autodoc_mock_imports = ["PyNucleus_base"]
diff --git a/drivers/brusselator.py b/drivers/brusselator.py
index 6616eefa..751e0296 100755
--- a/drivers/brusselator.py
+++ b/drivers/brusselator.py
@@ -23,13 +23,13 @@
 from mpi4py import MPI
 import numpy as np
 from pathlib import Path
-from PyNucleus.base import driver, solverFactory
-from PyNucleus.base.linear_operators import TimeStepperLinearOperator
-from PyNucleus.fem.femCy import assembleNonlinearity
-from PyNucleus.multilevelSolver import hierarchyManager
-from PyNucleus.nl import paramsForFractionalHierarchy
-from PyNucleus.nl.nonlocalProblems import brusselatorProblem
-from PyNucleus.base.timestepping import EulerIMEX, ARS3, koto
+from PyNucleus_base import driver, solverFactory
+from PyNucleus_base.linear_operators import TimeStepperLinearOperator
+from PyNucleus_fem.femCy import assembleNonlinearity
+from PyNucleus_multilevelSolver import hierarchyManager
+from PyNucleus_nl import paramsForFractionalHierarchy
+from PyNucleus_nl.nonlocalProblems import brusselatorProblem
+from PyNucleus_base.timestepping import EulerIMEX, ARS3, koto
 import h5py
 
 ###############################################################################
diff --git a/drivers/brusselatorMovie.py b/drivers/brusselatorMovie.py
index d7244913..95afb22a 100755
--- a/drivers/brusselatorMovie.py
+++ b/drivers/brusselatorMovie.py
@@ -12,10 +12,10 @@
 from shutil import rmtree
 import h5py
 from subprocess import Popen
-from PyNucleus.base import driver
-from PyNucleus.fem import meshNd
-from PyNucleus.fem.DoFMaps import DoFMap
-from PyNucleus.nl.nonlocalProblems import brusselatorProblem
+from PyNucleus_base import driver
+from PyNucleus_fem import meshNd
+from PyNucleus_fem.DoFMaps import DoFMap
+from PyNucleus_nl.nonlocalProblems import brusselatorProblem
 
 d = driver()
 brusselatorProblem(d)
diff --git a/drivers/interfaceProblem.py b/drivers/interfaceProblem.py
index ad269332..08044cb5 100755
--- a/drivers/interfaceProblem.py
+++ b/drivers/interfaceProblem.py
@@ -7,15 +7,15 @@
 ###################################################################################
 
 import numpy as np
-from PyNucleus.base import REAL, driver, solverFactory
-from PyNucleus.base.ip_norm import norm_serial
-from PyNucleus.fem.mesh import simpleInterval, uniformSquare
-from PyNucleus.fem.DoFMaps import P1_DoFMap, fe_vector
-from PyNucleus.fem import functionFactory
-from PyNucleus.fem.femCy import getSurfaceDoFMap
-from PyNucleus.fem.functions import squareIndicator, constant, Lambda
-from PyNucleus.fem import NO_BOUNDARY, INTERIOR, PHYSICAL
-from PyNucleus.fem.splitting import meshSplitter, dofmapSplitter
+from PyNucleus_base import REAL, driver, solverFactory
+from PyNucleus_base.ip_norm import norm_serial
+from PyNucleus_fem.mesh import simpleInterval, uniformSquare
+from PyNucleus_fem.DoFMaps import P1_DoFMap, fe_vector
+from PyNucleus_fem import functionFactory
+from PyNucleus_fem.femCy import getSurfaceDoFMap
+from PyNucleus_fem.functions import squareIndicator, constant, Lambda
+from PyNucleus_fem import NO_BOUNDARY, INTERIOR, PHYSICAL
+from PyNucleus_fem.splitting import meshSplitter, dofmapSplitter
 
 
 d = driver()
diff --git a/drivers/runFractional.py b/drivers/runFractional.py
index 61dfeb67..c24df61e 100755
--- a/drivers/runFractional.py
+++ b/drivers/runFractional.py
@@ -8,7 +8,7 @@
 
 from mpi4py import MPI
 from PyNucleus import driver
-from PyNucleus.nl import (fractionalLaplacianProblem,
+from PyNucleus_nl import (fractionalLaplacianProblem,
                           discretizedNonlocalProblem)
 from PyNucleus_nl.fractionalOrders import singleVariableUnsymmetricFractionalOrder
 
@@ -16,8 +16,9 @@
 
 
 d = driver(MPI.COMM_WORLD)
-d.add('saveOperators', False)
-d.add('vtkOutput', "")
+io = d.addGroup('input/output').group
+d.add('saveOperators', False, help="add operators to data file", group=io)
+d.add('vtkOutput', "", help="write solution to VTK file", group=io)
 p = fractionalLaplacianProblem(d, False)
 discrProblem = discretizedNonlocalProblem(d, p)
 
diff --git a/drivers/runFractionalHeat.py b/drivers/runFractionalHeat.py
index e8616b15..a70cdb6d 100755
--- a/drivers/runFractionalHeat.py
+++ b/drivers/runFractionalHeat.py
@@ -8,7 +8,7 @@
 
 from mpi4py import MPI
 from PyNucleus import driver
-from PyNucleus.nl import (transientFractionalProblem,
+from PyNucleus_nl import (transientFractionalProblem,
                           discretizedTransientProblem)
 
 ##################################################
diff --git a/drivers/runHelmholtz.py b/drivers/runHelmholtz.py
index cec00d80..5c143e7b 100755
--- a/drivers/runHelmholtz.py
+++ b/drivers/runHelmholtz.py
@@ -8,16 +8,16 @@
 
 from mpi4py import MPI
 import numpy as np
-from PyNucleus.base import COMPLEX, driver, solverFactory
-from PyNucleus.base.linear_operators import wrapRealToComplexCSR
-from PyNucleus.fem.femCy import (assembleSurfaceMass,
+from PyNucleus_base import COMPLEX, driver, solverFactory
+from PyNucleus_base.linear_operators import wrapRealToComplexCSR
+from PyNucleus_fem.femCy import (assembleSurfaceMass,
                                  getSurfaceDoFMap)
-from PyNucleus.fem import (PHYSICAL,
+from PyNucleus_fem import (PHYSICAL,
                            NO_BOUNDARY,
                            helmholtzProblem)
-from PyNucleus.fem.mesh import plotManager
-from PyNucleus.fem.functions import real, imag
-from PyNucleus.multilevelSolver import (EmptyHierarchy,
+from PyNucleus_fem.mesh import plotManager
+from PyNucleus_fem.functions import real, imag
+from PyNucleus_multilevelSolver import (EmptyHierarchy,
                                         hierarchyManager,
                                         inputConnector,
                                         paramsForMG)
diff --git a/drivers/runNonlocal.py b/drivers/runNonlocal.py
index 97122936..b2c02b93 100755
--- a/drivers/runNonlocal.py
+++ b/drivers/runNonlocal.py
@@ -8,7 +8,7 @@
 
 from mpi4py import MPI
 from PyNucleus import driver
-from PyNucleus.nl import (nonlocalPoissonProblem,
+from PyNucleus_nl import (nonlocalPoissonProblem,
                           discretizedNonlocalProblem)
 
 ##################################################
diff --git a/drivers/runNonlocalInterface.py b/drivers/runNonlocalInterface.py
index 9397147d..fcb3bb5e 100755
--- a/drivers/runNonlocalInterface.py
+++ b/drivers/runNonlocalInterface.py
@@ -7,17 +7,17 @@
 ###################################################################################
 
 import numpy as np
-from PyNucleus.base import driver, solverFactory
-from PyNucleus.base.solvers import iterative_solver
-from PyNucleus.base.ip_norm import ip_serial, norm_serial
-from PyNucleus.fem import NO_BOUNDARY, HOMOGENEOUS_NEUMANN
-from PyNucleus.fem.functions import indicatorFunctor
-from PyNucleus.fem.quadrature import simplexXiaoGimbutas
-from PyNucleus.fem.DoFMaps import fe_vector, str2DoFMap
-from PyNucleus.fem.splitting import meshSplitter, dofmapSplitter
-from PyNucleus.nl.helpers import getFracLapl
-from PyNucleus.nl import FRACTIONAL
-from PyNucleus.nl.nonlocalProblems import nonlocalInterfaceProblem
+from PyNucleus_base import driver, solverFactory
+from PyNucleus_base.solvers import iterative_solver
+from PyNucleus_base.ip_norm import ip_serial, norm_serial
+from PyNucleus_fem import NO_BOUNDARY, HOMOGENEOUS_NEUMANN
+from PyNucleus_fem.functions import indicatorFunctor
+from PyNucleus_fem.quadrature import simplexXiaoGimbutas
+from PyNucleus_fem.DoFMaps import fe_vector, str2DoFMap
+from PyNucleus_fem.splitting import meshSplitter, dofmapSplitter
+from PyNucleus_nl.helpers import getFracLapl
+from PyNucleus_nl import FRACTIONAL
+from PyNucleus_nl.nonlocalProblems import nonlocalInterfaceProblem
 
 d = driver()
 nIP = nonlocalInterfaceProblem(d)
diff --git a/drivers/runParallelGMG.py b/drivers/runParallelGMG.py
index f0a7e06b..3df1b2aa 100755
--- a/drivers/runParallelGMG.py
+++ b/drivers/runParallelGMG.py
@@ -9,10 +9,10 @@
 from __future__ import division, print_function
 from mpi4py import MPI
 import numpy as np
-from PyNucleus.base import driver, solverFactory
-from PyNucleus.fem import diffusionProblem
-from PyNucleus.fem.mesh import accumulate2global
-from PyNucleus.multilevelSolver import (V, FMG_V,
+from PyNucleus_base import driver, solverFactory
+from PyNucleus_fem import diffusionProblem
+from PyNucleus_fem.mesh import accumulate2global
+from PyNucleus_multilevelSolver import (V, FMG_V,
                                         hierarchyManager,
                                         inputConnector,
                                         pCoarsenConnector,
diff --git a/drivers/runSerialGMG.py b/drivers/runSerialGMG.py
index 92c8ad3f..da63e47f 100755
--- a/drivers/runSerialGMG.py
+++ b/drivers/runSerialGMG.py
@@ -9,9 +9,9 @@
 from __future__ import division
 import mpi4py.MPI as MPI
 import numpy as np
-from PyNucleus.base import driver, solverFactory
-from PyNucleus.fem import diffusionProblem
-from PyNucleus.multilevelSolver import (V, FMG_V,
+from PyNucleus_base import driver, solverFactory
+from PyNucleus_fem import diffusionProblem
+from PyNucleus_multilevelSolver import (V, FMG_V,
                                         hierarchyManager,
                                         inputConnector,
                                         paramsForSerialMG)
diff --git a/drivers/testDistOp.py b/drivers/testDistOp.py
index 8cd239c3..293de6ba 100755
--- a/drivers/testDistOp.py
+++ b/drivers/testDistOp.py
@@ -7,7 +7,7 @@
 ###################################################################################
 
 from mpi4py import MPI
-from PyNucleus.base import driver
+from PyNucleus_base import driver
 from PyNucleus import (dofmapFactory,
                        solverFactory,
                        nonlocalMeshFactory,
@@ -15,9 +15,9 @@
                        DIRICHLET,
                        NEUMANN,
                        HOMOGENEOUS_NEUMANN)
-from PyNucleus.fem.mesh import plotManager
-from PyNucleus.base.utilsFem import TimerManager, timerOutputGroup
-from PyNucleus.nl import nonlocalPoissonProblem, FRACTIONAL
+from PyNucleus_fem.mesh import plotManager
+from PyNucleus_base.utilsFem import TimerManager, timerOutputGroup
+from PyNucleus_nl import nonlocalPoissonProblem, FRACTIONAL
 from PyNucleus_nl.clusterMethodCy import (DistributedH2Matrix_globalData,
                                           DistributedH2Matrix_localData)
 import numpy as np
diff --git a/drivers/variableOrder.py b/drivers/variableOrder.py
index f559db49..5c83a45a 100755
--- a/drivers/variableOrder.py
+++ b/drivers/variableOrder.py
@@ -8,19 +8,19 @@
 
 import numpy as np
 from numpy.linalg import norm as npnorm
-from PyNucleus.base.ip_norm import ip_serial, norm_serial
-from PyNucleus.base import driver, solverFactory, krylov_solver
-from PyNucleus.base.linear_operators import invDiagonal
-from PyNucleus.fem import (str2DoFMap,
+from PyNucleus_base.ip_norm import ip_serial, norm_serial
+from PyNucleus_base import driver, solverFactory, krylov_solver
+from PyNucleus_base.linear_operators import invDiagonal
+from PyNucleus_fem import (str2DoFMap,
                            meshFactory,
                            functionFactory)
-from PyNucleus.fem.mesh import plotManager
-from PyNucleus.nl.fractionalOrders import (constFractionalOrder,
+from PyNucleus_fem.mesh import plotManager
+from PyNucleus_nl.fractionalOrders import (constFractionalOrder,
                                            variableConstFractionalOrder,
                                            leftRightFractionalOrder,
                                            smoothedLeftRightFractionalOrder,
                                            innerOuterFractionalOrder)
-from PyNucleus.nl.kernels import getFractionalKernel
+from PyNucleus_nl.kernels import getFractionalKernel
 
 
 d = driver()
diff --git a/examples/example_nonlocal.py b/examples/example_nonlocal.py
index 48616867..a519a6e4 100755
--- a/examples/example_nonlocal.py
+++ b/examples/example_nonlocal.py
@@ -51,7 +51,7 @@
 import matplotlib.pyplot as plt
 from time import time
 from PyNucleus import kernelFactory
-kernelFactory.print()
+print(kernelFactory)
 
 # %%
 
diff --git a/examples/example_pde.py b/examples/example_pde.py
index c1870af5..90153911 100755
--- a/examples/example_pde.py
+++ b/examples/example_pde.py
@@ -19,7 +19,7 @@
 # ---------
 #
 # The creation of different groups of objects, such as finite element spaces or meshes, use factories.
-# The available classes that a factory provides can be displayed by calling the ``print()`` method of the factory.
+# The available classes that a factory provides can be displayed by calling the ``print()`` function on it.
 # An object is built by passing the name of the desired class and additional parameters to the factory.
 # If this sounds vague now, don't worry, the examples below will make it clear.
 #
@@ -31,7 +31,7 @@
 
 import matplotlib.pyplot as plt
 from PyNucleus import meshFactory
-meshFactory.print()
+print(meshFactory)
 
 # %%
 # We see what the different meshes that the ``meshFactory`` can construct and the default values for associated parameters.
@@ -60,7 +60,7 @@
 # We build a piecewise linear finite element space.
 
 from PyNucleus import dofmapFactory
-dofmapFactory.print()
+print(dofmapFactory)
 
 # %%
 # We use a piecewise linear continuous finite element space.
@@ -95,7 +95,7 @@
 #
 # Functions are created using the ``functionFactory``.
 from PyNucleus import functionFactory
-functionFactory.print()
+print(functionFactory)
 
 # %%
 # We will consider two different forcing functions :math:`f`.
@@ -154,7 +154,7 @@
 # Now that we have assembled our linear system, we want to solve it.
 
 from PyNucleus import solverFactory
-solverFactory.print()
+print(solverFactory)
 
 # %%
 # We choose to set up an LU direct solver.
diff --git a/fem/PyNucleus_fem/mesh.py b/fem/PyNucleus_fem/mesh.py
index 96416b74..e01aa608 100644
--- a/fem/PyNucleus_fem/mesh.py
+++ b/fem/PyNucleus_fem/mesh.py
@@ -11,7 +11,7 @@
 from mpi4py import MPI
 import numpy as np
 from PyNucleus_base.factory import factory
-from PyNucleus_base.myTypes import INDEX, REAL, TAG
+from PyNucleus_base.myTypes import INDEX, REAL, BOOL, TAG
 from PyNucleus_base.linear_operators import sparseGraph
 from PyNucleus_base import uninitialized, uninitialized_like
 from . meshCy import (meshBase,
@@ -581,7 +581,7 @@ def doubleSquareWithInteractionsCorners(ax=0., ay=0., bx=1., by=1., cx=2., cy=1.
                                         h=None,
                                         returnSketch=False,
                                         **kwargs):
-    from PyNucleus.fem.meshConstruction import (line,
+    from PyNucleus_fem.meshConstruction import (line,
                                                 polygon)
     assert horizon2 >= horizon1
     assert horizon1 >= 0
@@ -2993,7 +2993,7 @@ def stitchOverlappingMeshes(meshes, overlapManagers):
     globalCellLookup = []
     for mySubdomainNo in range(numPartitions):
         translate = -np.ones((meshes[mySubdomainNo].num_vertices), dtype=INDEX)
-        idx = np.ones((meshes[mySubdomainNo].cells.shape[0]), dtype=np.bool)
+        idx = np.ones((meshes[mySubdomainNo].cells.shape[0]), dtype=BOOL)
         lookup = -np.ones((meshes[mySubdomainNo].num_cells), dtype=INDEX)
         for otherSubdomainNo in range(mySubdomainNo):
             if otherSubdomainNo not in overlapManagers[mySubdomainNo].overlaps:
diff --git a/nl/PyNucleus_nl/discretizedProblems.py b/nl/PyNucleus_nl/discretizedProblems.py
index 855d87b2..04e28e8a 100644
--- a/nl/PyNucleus_nl/discretizedProblems.py
+++ b/nl/PyNucleus_nl/discretizedProblems.py
@@ -22,7 +22,6 @@
 from PyNucleus_multilevelSolver import hierarchyManager
 from copy import copy
 from . helpers import paramsForFractionalHierarchy
-from . kernelsCy import FRACTIONAL
 from . nonlocalProblems import (DIRICHLET,
                                 NEUMANN, HOMOGENEOUS_NEUMANN,
                                 transientFractionalProblem)
@@ -196,7 +195,7 @@ def plotSolution(self):
             plt.gca().set_aspect('equal')
 
     def plotSolutionComponents(self, plotDefaults={}):
-        from PyNucleus.fem.mesh import plotManager
+        from PyNucleus_fem.mesh import plotManager
         pm = plotManager(self.u.dm.scalarDM.mesh, self.u.dm.scalarDM,
                          defaults=plotDefaults)
         for c in range(self.u.dm.numComponents):
@@ -388,7 +387,7 @@ def buildMeshHierarchy(self, mesh, solverType, domainIndicator, fluxIndicator, n
             params['assemble'] = 'restrictionProlongation' if solverType.find('mg') >= 0 else 'dofmap only last'
             params['logging'] = True
             if self.debugAssemblyTimes:
-                from PyNucleus.base.utilsFem import TimerManager
+                from PyNucleus_base.utilsFem import TimerManager
                 tm = TimerManager(self.driver.logger, comm=self.driver.comm, memoryProfiling=False, loggingSubTimers=True)
                 params['PLogger'] = tm.PLogger
 
@@ -673,6 +672,8 @@ def adjointSolve(self, b, dm, dmInterior, P_interior, adjointSolver, tol, maxite
 
     def report(self, group):
         group.add('kernel', repr(self.continuumProblem.kernel))
+        group.add('problem', self.continuumProblem.problemDescription)
+        group.add('has analytic solution', self.continuumProblem.analyticSolution is not None)
         group.add('h', self.finalMesh.h)
         group.add('hmin', self.finalMesh.hmin)
         if self.continuumProblem.kernel is not None:
diff --git a/nl/PyNucleus_nl/nonlocalProblems.py b/nl/PyNucleus_nl/nonlocalProblems.py
index 6af77b50..ba65ce1e 100644
--- a/nl/PyNucleus_nl/nonlocalProblems.py
+++ b/nl/PyNucleus_nl/nonlocalProblems.py
@@ -525,15 +525,18 @@ def setDriverArgs(self):
         p = self.driver.addGroup('problem')
         self.setDriverFlag('domain', acceptedValues=['interval', 'disc', 'gradedInterval', 'gradedDisc', 'Lshape', 'square',
                                                      'cutoutCircle', 'disconnectedInterval', 'disconnectedDomain',
-                                                     'ball'], group=p)
+                                                     'ball'],
+                           help="spatial domain", group=p)
         self.setDriverFlag('problem', acceptedValues=['constant', 'notPeriodic', 'plateau',
-                                                      'sin', 'cos', 3, 'source', 'zeroFlux', 'Greens', 'knownSolution'], group=p)
-        self.setDriverFlag('element', acceptedValues=['P1', 'P2', 'P3', 'P0'], group=p)
+                                                      'sin', 'cos', 3, 'source', 'zeroFlux', 'Greens', 'knownSolution'],
+                           help="select a problem to solve", group=p)
+        self.setDriverFlag('element', acceptedValues=['P1', 'P2', 'P3', 'P0'], help="finite element space", group=p)
         self.setDriverFlag('adaptive', acceptedValues=['residualMelenk', 'residualNochetto',
                                                        'residual', 'hierarchical', 'knownSolution', None],
-                           argInterpreter=lambda v: None if v == 'None' else v, group=p)
-        self.setDriverFlag('noRef', -1, group=p)
-        self.setDriverFlag('targetDoFsAux', 0)
+                           argInterpreter=lambda v: None if v == 'None' else v,
+                           help="type of error indicator", group=p)
+        self.setDriverFlag('noRef', -1, help="number of uniform mesh refinements applied to initial mesh", group=p)
+        self.setDriverFlag('targetDoFsAux', 0, help="number of degrees of freedom for the fractional order (0 = do not discretize order)")
 
     def processCmdline(self, params):
         noRef = params['noRef']
@@ -617,7 +620,7 @@ def getDomainParams(self, domain):
                 'mesh_domain', 'mesh_params', 'tag', 'boundaryCondition',
                 'domainIndicator', 'interactionIndicator', 'fluxIndicator',
                 'zeroExterior',
-                'rhsData', 'dirichletData', 'fluxData'])
+                'rhsData', 'dirichletData', 'fluxData', 'problemDescription'])
     def processProblem(self, kernel, dim, domain, domainParams, problem, normalized):
         s = kernel.s
         self.analyticSolution = None
@@ -631,6 +634,7 @@ def processProblem(self, kernel, dim, domain, domainParams, problem, normalized)
             radius = 1.
 
             if problem == 'constant':
+                self.problemDescription = "Constant forcing, homogeneous Dirichlet volume condition"
                 self.rhs = constant(1.)
                 if (isinstance(s, (constFractionalOrder, variableConstFractionalOrder, constantNonSymFractionalOrder)) or
                         (isinstance(s, feFractionalOrder) and np.array(s.vec).min() == np.array(s.vec).max())):
@@ -643,10 +647,13 @@ def processProblem(self, kernel, dim, domain, domainParams, problem, normalized)
                     L2_ex = np.sqrt(C**2 * np.sqrt(np.pi) * Gamma(1+2*sValue)/Gamma(3/2+2*sValue) * radius**2)
                     self.analyticSolution = solFractional(sValue, dim, radius)
             elif problem == 'sin':
+                self.problemDescription = "Sin function forcing, homogeneous Dirichlet volume condition"
                 self.rhs = Lambda(lambda x: np.sin(np.pi*x[0]))
             elif problem == 'cos':
+                self.problemDescription = "Cos function forcing, homogeneous Dirichlet volume condition"
                 self.rhs = Lambda(lambda x: np.cos(np.pi*x[0]/2.))
             elif problem == 'plateau':
+                self.problemDescription = "Sign function forcing, homogeneous Dirichlet volume condition"
                 self.rhs = Lambda(np.sign)
 
                 # def e(n):
@@ -656,10 +663,12 @@ def processProblem(self, kernel, dim, domain, domainParams, problem, normalized)
                 # exactHsSquared = sum([e(n) for n in range(1000000)])
                 self.exactHsSquared = 2**(1-2*s) / (2*s+1) / Gamma(s+1)**2
             elif isinstance(problem, int):
+                self.problemDescription = "Family of forcings with known solutions, homogeneous Dirichlet volume condition"
                 self.rhs = rhsFractional1D(s, problem)
                 self.exactHsSquared = 2**(2*s)/(2*problem+s+0.5) * Gamma(1+s)**2 * binom(s+problem, problem)**2
                 self.analyticSolution = solFractional1D(s, problem)
             elif problem == 'zeroFlux':
+                self.problemDescription = "Linear solution, homogeneous Neumann volume condition"
                 boundaryCondition = HOMOGENEOUS_NEUMANN
 
                 if kernel.variable:
@@ -684,6 +693,7 @@ def fun(x):
                 self.analyticSolution = functionFactory('x0')
                 L2_ex = np.sqrt(2/3)
             elif problem == 'knownSolution':
+                self.problemDescription = "Known analytic solution for variable fractional order, homogeneous Dirichlet volume condition"
                 from scipy.special import hyp2f1
                 assert isinstance(s, (constFractionalOrder, variableConstFractionalOrder,
                                       constantNonSymFractionalOrder, singleVariableUnsymmetricFractionalOrder)), s
@@ -700,12 +710,14 @@ def fun(x):
                 self.analyticSolution = functionFactory('Lambda', lambda x: (1.-x[0]**2)**beta)
                 L2_ex = np.sqrt(np.sqrt(np.pi) * Gamma(1+2*beta)/Gamma(3/2+2*beta) * radius**2)
             elif problem == 'Greens':
+                self.problemDescription = "Narrow indicator function forcing, homogeneous Neumann volume condition"
                 boundaryCondition = HOMOGENEOUS_NEUMANN
                 self.rhs = functionFactory('squareIndicator', np.array([-0.1]), np.array([0.1]))
             else:
                 raise NotImplementedError(problem)
         elif domain == 'disconnectedInterval':
             if problem == 'constant':
+                self.problemDescription = "Constant forcing, homogeneous Dirichlet volume condition"
                 self.rhs = Lambda(lambda x: 1. if x[0] > 0.5 else 0.)
             else:
                 raise NotImplementedError(problem)
@@ -713,6 +725,7 @@ def fun(x):
             radius = 1.
 
             if problem == 'constant':
+                self.problemDescription = "Constant forcing, homogeneous Dirichlet volume condition"
                 self.rhs = constant(1.)
                 if isinstance(s, (constFractionalOrder, variableConstFractionalOrder, constantNonSymFractionalOrder)):
                     C = 2.**(-2.*s.value)*Gamma(dim/2.)/Gamma((dim+2.*s.value)/2.)/Gamma(1.+s.value)
@@ -720,6 +733,7 @@ def fun(x):
                     L2_ex = np.sqrt(C**2 * np.pi/(1+2*s.value)*radius**2)
                     self.analyticSolution = solFractional(s.value, dim, radius)
             elif problem == 'notPeriodic':
+                self.problemDescription = "Family of forcings with known solutions, homogeneous Dirichlet volume condition"
                 n = 2
                 freq = 2
                 self.exactHsSquared = 2**(2*s-1)/(2*n+s+freq+1) * Gamma(1+s+n)**2/Gamma(1+n)**2 * (np.pi+np.sin(4*np.pi*freq)/(4*freq))
@@ -729,6 +743,7 @@ def fun(x):
                 self.exactHsSquared += 2**(2*s-1)/(2*n+s+freq+1) * Gamma(1+s+n)**2/Gamma(1+n)**2 * (np.pi+np.sin(4*np.pi*freq)/(4*freq))
                 self.rhs = rhsFractional2D_nonPeriodic(s)
             elif problem == 'plateau':
+                self.problemDescription = "Sign function forcing, homogeneous Dirichlet volume condition"
                 self.rhs = Lambda(lambda x: x[0] > 0)
                 try:
                     from mpmath import meijerg
@@ -742,13 +757,16 @@ def fun(x):
                     for k in range(100000):
                         self.exactHsSquared += 2**(-2*s) / Gamma(s+3)**2 / (2*np.pi) * (2*k+s+2) * (k+1) / binom(k+s+1.5, s+2)**2
             elif isinstance(problem, tuple):
+                self.problemDescription = "Family of forcings with known solutions, homogeneous Dirichlet volume condition"
                 n, freq = problem
                 self.exactHsSquared = 2**(2*s-1)/(2*n+s+freq+1) * Gamma(1+s+n)**2/Gamma(1+n)**2 * (np.pi+np.sin(4*np.pi*freq)/(4*freq))
 
                 self.rhs = rhsFractional2D(s, n=n, l=freq)
             elif problem == 'sin':
+                self.problemDescription = "Radial sin function forcing, homogeneous Dirichlet volume condition"
                 self.rhs = Lambda(lambda x: np.sin(np.pi*(x[0]**2+x[1]**2)))
             elif problem == 'knownSolution':
+                self.problemDescription = "Known analytic solution for variable fractional order, homogeneous Dirichlet volume condition"
                 from scipy.special import hyp2f1
                 assert isinstance(s, (constFractionalOrder, variableConstFractionalOrder,
                                       constantNonSymFractionalOrder, singleVariableUnsymmetricFractionalOrder)), s
@@ -768,8 +786,10 @@ def fun(x):
                 raise NotImplementedError(problem)
         elif domain == 'square':
             if problem == 'constant':
+                self.problemDescription = "Constant forcing, homogeneous Dirichlet volume condition"
                 self.rhs = constant(1.)
             elif problem == 'sin':
+                self.problemDescription = "Tensor sin function forcing, homogeneous Dirichlet volume condition"
                 self.rhs = Lambda(lambda x: np.sin(np.pi*x[0])*np.sin(np.pi*x[1]))
             elif problem == 'source':
                 self.rhs = (functionFactory.build('radialIndicator', radius=0.3, center=np.array([0.2, 0.6], dtype=REAL)) -
@@ -778,21 +798,26 @@ def fun(x):
                 raise NotImplementedError(problem)
         elif domain == 'Lshape':
             if problem == 'constant':
+                self.problemDescription = "Constant forcing, homogeneous Dirichlet volume condition"
                 self.rhs = constant(1.)
             elif problem == 'sin':
+                self.problemDescription = "Tensor sin function forcing, homogeneous Dirichlet volume condition"
                 self.rhs = Lambda(lambda x: np.sin(np.pi*x[0])*np.sin(np.pi*x[1]))
             else:
                 raise NotImplementedError(problem)
         elif domain == 'cutoutCircle':
             if problem == 'constant':
+                self.problemDescription = "Constant forcing, homogeneous Dirichlet volume condition"
                 self.rhs = constant(1.)
             elif problem == 'sin':
+                self.problemDescription = "Radial sin function forcing, homogeneous Dirichlet volume condition"
                 self.rhs = Lambda(lambda x: np.sin(np.pi*(x[0]**2+x[1]**2)))
             else:
                 raise NotImplementedError(problem)
         elif domain == 'ball':
             radius = 1.
             if problem == 'constant':
+                self.problemDescription = "Constant forcing, homogeneous Dirichlet volume condition"
                 self.rhs = constant(1.)
                 if isinstance(s, (constFractionalOrder, variableConstFractionalOrder, constantNonSymFractionalOrder)):
                     C = 2.**(-2.*s.value)*Gamma(dim/2.)/Gamma((dim+2.*s.value)/2.)/Gamma(1.+s.value)
@@ -869,16 +894,6 @@ def constructAuxiliarySpace(self, dim, domain, domainParams, kernelType, horizon
             mesh = mesh.refine()
         self.dmAux = dofmapFactory('P1', mesh, NO_BOUNDARY)
 
-    def getIdentifier(self, params):
-        keys = ['domain', 'problem', 's', 'noRef', 'element', 'adaptive']
-        d = []
-        for k in keys:
-            try:
-                d.append((k, str(getattr(self, k))))
-            except KeyError:
-                d.append((k, str(params[k])))
-        return '-'.join(['fracLaplAdaptive'] + [key + '=' + v for key, v in d])
-
 
 class nonlocalPoissonProblem(nonlocalBaseProblem):
     def setDriverArgs(self):
@@ -919,7 +934,8 @@ def processCmdline(self, params):
                 'tag', 'zeroExterior', 'boundaryCondition',
                 'domainIndicator', 'fluxIndicator', 'interactionIndicator',
                 'rhs', 'rhsData', 'dirichletData', 'fluxData',
-                'analyticSolution', 'exactL2Squared', 'exactHsSquared'])
+                'analyticSolution', 'exactL2Squared', 'exactHsSquared',
+                'problemDescription'])
     def processProblem(self, kernel, domain, problem, normalized):
         if kernel is not None:
             kType = kernel.kernelType
@@ -938,6 +954,8 @@ def processProblem(self, kernel, domain, problem, normalized):
         self.exactL2Squared = None
         self.exactHsSquared = None
 
+        self.problemDescription = ""
+
         if problem in ('poly-Neumann', 'exact-sin-Neumann', 'zeroFlux'):
             self.boundaryCondition = NEUMANN
         elif self.parametrizedArg('indicator').match(problem):
@@ -1421,16 +1439,6 @@ def buildMesh(self, mesh_domain, mesh_params, hTarget):
         myMeshParams['hTarget'] = hTarget
         self.mesh, _ = nonlocalMeshFactory.build(mesh_domain, **myMeshParams)
 
-    def getIdentifier(self, params):
-        keys = ['domain', 'problem', 's', 'horizon', 'phi', 'noRef']
-        d = []
-        for k in keys:
-            try:
-                d.append((k, str(getattr(self, k))))
-            except KeyError:
-                d.append((k, str(params[k])))
-        return '-'.join(['nonlocal'] + [key + '=' + v for key, v in d])
-
 
 class transientFractionalProblem(fractionalLaplacianProblem):
     def __init__(self, driver, useMulti=False):
@@ -2254,16 +2262,6 @@ def flux_right_lam(x):
         self.local_H10ex_left = local_H10ex_left
         self.local_H10ex_right = local_H10ex_right
 
-    def getIdentifier(self, params):
-        keys = ['domain', 'problem', 'kernel1', 'kernel2']
-        d = []
-        for k in keys:
-            try:
-                d.append((k, str(getattr(self, k))))
-            except KeyError:
-                d.append((k, str(params[k])))
-        return '/'.join(['nonlocalInterface'] + [key + '=' + v for key, v in d])
-
 
 class brusselatorProblem(problem):
     """
@@ -2308,7 +2306,7 @@ def setDriverArgs(self):
                 'mesh',
                 'zeroExterior'])
     def processProblem(self, domain, bc, noRef, problem, T):
-        from PyNucleus.fem.femCy import brusselator
+        from PyNucleus_fem.femCy import brusselator
 
         if problem == 'spots':
             self.alpha = self.beta = 0.75
@@ -2408,14 +2406,4 @@ def initial_V(x):
                                                 boundaryCondition=self.boundaryCondition)
         self.zeroExterior = nI['zeroExterior']
 
-    def getIdentifier(self, params):
-        keys = ['domain', 'problem', 'alpha', 'beta', 'noRef', 'bc']
-        d = []
-        for k in keys:
-            try:
-                d.append((k, str(getattr(self, k))))
-            except KeyError:
-                d.append((k, str(params[k])))
-        return '-'.join(['brusselator'] + [key + '=' + v for key, v in d])
-
 
diff --git a/pytest.ini b/pytest.ini
new file mode 100644
index 00000000..c29cd39a
--- /dev/null
+++ b/pytest.ini
@@ -0,0 +1,8 @@
+[pytest]
+filterwarnings =
+    ignore:unit_vertices:DeprecationWarning
+markers =
+    slow: marks tests as slow (deselect with '-m "not slow"')
+testpaths =
+    tests
+generate_report_on_test = True
\ No newline at end of file
diff --git a/tests/bitArray.py b/tests/bitArray.py
index c302f337..59909b68 100644
--- a/tests/bitArray.py
+++ b/tests/bitArray.py
@@ -6,8 +6,8 @@
 ###################################################################################
 
 import numpy as np
-from PyNucleus.base import INDEX
-from PyNucleus.base.tupleDict import arrayIndexSet, bitArray
+from PyNucleus_base import INDEX
+from PyNucleus_base.tupleDict import arrayIndexSet, bitArray
 
 
 def test_arrayIndexSet():
diff --git a/tests/drivers_base.py b/tests/drivers_base.py
index b3d9abfd..f1fdbe02 100644
--- a/tests/drivers_base.py
+++ b/tests/drivers_base.py
@@ -6,7 +6,7 @@
 ###################################################################################
 
 from mpi4py import MPI
-from PyNucleus.base.utilsFem import runDriver
+from PyNucleus_base.utilsFem import runDriver
 import os
 import inspect
 import pytest
@@ -42,15 +42,15 @@ def symmetric(request):
 
 
 @pytest.mark.slow
-def testGMG(extra):
+def testGMG(extras):
     base = getPath()+'/../'
     py = 'runSerialGMG.py'
     path = base+'drivers'
     cacheDir = getPath()+'/'
-    runDriver(path, py, cacheDir=cacheDir, extra=extra)
+    runDriver(path, py, cacheDir=cacheDir, extra=extras)
 
 
-def testParallelGMG(ranks, domain, element, symmetric, extra):
+def testParallelGMG(ranks, domain, element, symmetric, extras):
     base = getPath()+'/../'
     py = ['runParallelGMG.py',
           '--domain', domain,
@@ -59,18 +59,18 @@ def testParallelGMG(ranks, domain, element, symmetric, extra):
         py.append('--symmetric')
     path = base+'drivers'
     cacheDir = getPath()+'/'
-    runDriver(path, py, ranks=ranks, cacheDir=cacheDir, relTol=3e-2, extra=extra)
+    runDriver(path, py, ranks=ranks, cacheDir=cacheDir, relTol=3e-2, extra=extras)
 
 
 ################################################################################
 # multigrid for Helmholtz
 
-def testHelmholtz(ranks, domain, extra):
+def testHelmholtz(ranks, domain, extras):
     base = getPath()+'/../'
     py = ['runHelmholtz.py', '--domain', domain]
     path = base+'drivers'
     cacheDir = getPath()+'/'
-    runDriver(path, py, ranks=ranks, cacheDir=cacheDir, extra=extra)
+    runDriver(path, py, ranks=ranks, cacheDir=cacheDir, extra=extras)
 
 
 ################################################################################
@@ -82,7 +82,7 @@ def domainNoRef(request):
     return request.param
 
 
-def testInterface(domainNoRef, extra):
+def testInterface(domainNoRef, extras):
     domain, noRef = domainNoRef
     base = getPath()+'/../'
     py = ['interfaceProblem.py',
@@ -90,4 +90,4 @@ def testInterface(domainNoRef, extra):
           '--noRef', str(noRef)]
     path = base+'drivers'
     cacheDir = getPath()+'/'
-    runDriver(path, py, ranks=1, cacheDir=cacheDir, relTol=5e-2, extra=extra)
+    runDriver(path, py, ranks=1, cacheDir=cacheDir, relTol=5e-2, extra=extras)
diff --git a/tests/test_drivers_intFracLapl.py b/tests/test_drivers_intFracLapl.py
index 4804bf7f..be42f887 100644
--- a/tests/test_drivers_intFracLapl.py
+++ b/tests/test_drivers_intFracLapl.py
@@ -5,7 +5,7 @@
 # If you want to use this code, please refer to the README.rst and LICENSE files. #
 ###################################################################################
 
-from PyNucleus.base.utilsFem import runDriver
+from PyNucleus_base.utilsFem import runDriver
 import os
 import inspect
 import pytest
@@ -52,8 +52,12 @@ def runNonlocal_params(request):
 
 
 @pytest.mark.slow
-def testNonlocal(runNonlocal_params, extra):
-    domain, kernel, problem, solver, matrixFormat = runNonlocal_params
+def testNonlocal(runNonlocal_params, extras):
+    if len(runNonlocal_params) == 5:
+        domain, kernel, problem, solver, matrixFormat = runNonlocal_params
+        interaction = None
+    else:
+        domain, kernel, problem, solver, matrixFormat, interaction = runNonlocal_params
     base = getPath()+'/../'
     py = ['runNonlocal.py',
           '--domain', domain,
@@ -61,12 +65,14 @@ def testNonlocal(runNonlocal_params, extra):
           '--problem', problem,
           '--solver', solver,
           '--matrixFormat', matrixFormat]
+    if interaction is not None:
+        py += ['--interaction', interaction]
     # if kernel != 'fractional':
     path = base+'drivers'
     cacheDir = getPath()+'/'
     if problem == 'poly-Neumann' and domain == 'square':
         return pytest.skip('not implemented')
-    runDriver(path, py, cacheDir=cacheDir, extra=extra)
+    runDriver(path, py, cacheDir=cacheDir, extra=extras)
 
 
 @pytest.fixture(scope='module', params=[
@@ -106,7 +112,7 @@ def runFractional_params(request):
 
 
 @pytest.mark.slow
-def testFractional(runFractional_params, extra):
+def testFractional(runFractional_params, extras):
     domain, s, problem, element, solver, matrixFormat, ranks = runFractional_params
     base = getPath()+'/../'
     py = ['runFractional.py',
@@ -120,11 +126,11 @@ def testFractional(runFractional_params, extra):
     cacheDir = getPath()+'/'
     if ranks == 1:
         ranks = None
-    runDriver(path, py, cacheDir=cacheDir, extra=extra, ranks=ranks)
+    runDriver(path, py, cacheDir=cacheDir, extra=extras, ranks=ranks)
 
 
 @pytest.mark.slow
-def testFractionalHeat(runFractional_params, extra):
+def testFractionalHeat(runFractional_params, extras):
     domain, s, problem, element, solver, matrixFormat, ranks = runFractional_params
     base = getPath()+'/../'
     py = ['runFractionalHeat.py',
@@ -138,16 +144,16 @@ def testFractionalHeat(runFractional_params, extra):
     cacheDir = getPath()+'/'
     if ranks == 1:
         ranks = None
-    runDriver(path, py, cacheDir=cacheDir, extra=extra, ranks=ranks)
+    runDriver(path, py, cacheDir=cacheDir, extra=extras, ranks=ranks)
 
 
 @pytest.mark.slow
-def testVariableOrder(extra):
+def testVariableOrder(extras):
     base = getPath()+'/../'
     py = 'variableOrder.py'
     path = base+'drivers'
     cacheDir = getPath()+'/'
-    runDriver(path, py, cacheDir=cacheDir, extra=extra)
+    runDriver(path, py, cacheDir=cacheDir, extra=extras)
 
 
 @pytest.fixture(scope='module', params=[
@@ -172,7 +178,7 @@ def runDistOp_params(request):
 
 
 @pytest.mark.slow
-def testMatvecs(runDistOp_params, extra):
+def testMatvecs(runDistOp_params, extras):
     base = getPath()+'/../'
     domain, fractionalOrder = runDistOp_params
     if domain == 'interval':
@@ -196,7 +202,7 @@ def testMatvecs(runDistOp_params, extra):
     py += ['--no-write']
     path = base+'drivers'
     cacheDir = getPath()+'/'
-    runDriver(path, py, ranks=4, cacheDir=cacheDir, extra=extra)
+    runDriver(path, py, ranks=4, cacheDir=cacheDir, extra=extras)
 
 
 @pytest.fixture(scope='module', params=[
@@ -219,7 +225,7 @@ def runNonlocalInterface_params(request):
 
 
 @pytest.mark.slow
-def testNonlocalInterface(runNonlocalInterface_params, extra):
+def testNonlocalInterface(runNonlocalInterface_params, extras):
     domain, kernel1, kernel2, s11, s22, horizon1, horizon2, problem = runNonlocalInterface_params
     s12 = s11
     s21 = s22
@@ -237,4 +243,4 @@ def testNonlocalInterface(runNonlocalInterface_params, extra):
           '--problem', problem]
     path = base+'drivers'
     cacheDir = getPath()+'/'
-    runDriver(path, py, cacheDir=cacheDir, extra=extra)
+    runDriver(path, py, cacheDir=cacheDir, extra=extras)
diff --git a/tests/test_fracLapl.py b/tests/test_fracLapl.py
index 0ccd49f2..0c9b638e 100644
--- a/tests/test_fracLapl.py
+++ b/tests/test_fracLapl.py
@@ -8,17 +8,17 @@
 from __future__ import division
 import numpy as np
 import numpy.linalg
-from PyNucleus.fem.mesh import simpleInterval, circle
-from PyNucleus.fem.DoFMaps import P1_DoFMap, P2_DoFMap
-from PyNucleus.fem.functions import constant
-from PyNucleus.nl.nonlocalAssembly import (assembleNonlocalOperator,
+from PyNucleus_fem.mesh import simpleInterval, circle
+from PyNucleus_fem.DoFMaps import P1_DoFMap, P2_DoFMap
+from PyNucleus_fem.functions import constant
+from PyNucleus_nl.nonlocalAssembly import (assembleNonlocalOperator,
                                            nonlocalBuilder)
-from PyNucleus.nl.clusterMethodCy import H2Matrix
-from PyNucleus.base.myTypes import REAL
+from PyNucleus_nl.clusterMethodCy import H2Matrix
+from PyNucleus_base.myTypes import REAL
 from scipy.special import gamma
-from PyNucleus.nl.kernels import getFractionalKernel
-from PyNucleus.nl.fractionalOrders import constFractionalOrder
-from PyNucleus.nl.kernelNormalization import variableFractionalLaplacianScaling
+from PyNucleus_nl.kernels import getFractionalKernel
+from PyNucleus_nl.fractionalOrders import constFractionalOrder
+from PyNucleus_nl.kernelNormalization import variableFractionalLaplacianScaling
 import pytest
 import logging
 LOGGER = logging.getLogger(__name__)
diff --git a/tests/test_h2finiteHorizon.py b/tests/test_h2finiteHorizon.py
index 19375952..fd076052 100644
--- a/tests/test_h2finiteHorizon.py
+++ b/tests/test_h2finiteHorizon.py
@@ -6,14 +6,14 @@
 ###################################################################################
 
 import numpy as np
-from PyNucleus.fem.mesh import intervalWithInteraction
-from PyNucleus.fem.DoFMaps import P1_DoFMap
-from PyNucleus.fem.functions import Lambda, constant
-from PyNucleus.nl.fractionalOrders import (constFractionalOrder,
+from PyNucleus_fem.mesh import intervalWithInteraction
+from PyNucleus_fem.DoFMaps import P1_DoFMap
+from PyNucleus_fem.functions import Lambda, constant
+from PyNucleus_nl.fractionalOrders import (constFractionalOrder,
                                            variableConstFractionalOrder)
-from PyNucleus.nl.nonlocalAssembly import nonlocalBuilder
-from PyNucleus.nl.kernelNormalization import variableFractionalLaplacianScaling
-from PyNucleus.nl.kernels import getFractionalKernel
+from PyNucleus_nl.nonlocalAssembly import nonlocalBuilder
+from PyNucleus_nl.kernelNormalization import variableFractionalLaplacianScaling
+from PyNucleus_nl.kernels import getFractionalKernel
 from scipy.linalg import solve
 import pytest
 
diff --git a/tests/test_nearField.py b/tests/test_nearField.py
index 73b32104..725b8a20 100644
--- a/tests/test_nearField.py
+++ b/tests/test_nearField.py
@@ -9,24 +9,24 @@
 
 from mpi4py import MPI
 import numpy as np
-from PyNucleus.base.myTypes import REAL, INDEX
-from PyNucleus.base import uninitialized
-from PyNucleus.base.tupleDict import arrayIndexSet
-from PyNucleus.fem.DoFMaps import P0_DoFMap, P1_DoFMap
-from PyNucleus.fem.functions import constant
-from PyNucleus.nl import nonlocalBuilder
-from PyNucleus.nl.kernels import getFractionalKernel
-from PyNucleus.nl.nonlocalAssembly import nearFieldClusterPair
-from PyNucleus.nl.clusterMethodCy import (H2Matrix,
+from PyNucleus_base.myTypes import REAL, INDEX
+from PyNucleus_base import uninitialized
+from PyNucleus_base.tupleDict import arrayIndexSet
+from PyNucleus_fem.DoFMaps import P0_DoFMap, P1_DoFMap
+from PyNucleus_fem.functions import constant
+from PyNucleus_nl import nonlocalBuilder
+from PyNucleus_nl.kernels import getFractionalKernel
+from PyNucleus_nl.nonlocalAssembly import nearFieldClusterPair
+from PyNucleus_nl.clusterMethodCy import (H2Matrix,
                                           getDoFBoxesAndCells,
                                           tree_node)
-from PyNucleus.nl.fractionalOrders import (constFractionalOrder,
+from PyNucleus_nl.fractionalOrders import (constFractionalOrder,
                                            variableConstFractionalOrder,
                                            leftRightFractionalOrder,
                                            layersFractionalOrder,
                                            singleVariableUnsymmetricFractionalOrder)
-from PyNucleus.base import driver
-from PyNucleus.nl.nonlocalProblems import nonlocalMeshFactory, HOMOGENEOUS_DIRICHLET
+from PyNucleus_base import driver
+from PyNucleus_nl.nonlocalProblems import nonlocalMeshFactory, HOMOGENEOUS_DIRICHLET
 import pytest
 
 epsRelDense = 3e-2
diff --git a/tests/tupleDict.py b/tests/tupleDict.py
index 4a2d5fe3..17a42344 100644
--- a/tests/tupleDict.py
+++ b/tests/tupleDict.py
@@ -6,8 +6,8 @@
 ###################################################################################
 
 import numpy as np
-from PyNucleus.base import INDEX
-from PyNucleus.base.tupleDict import tupleDictINDEX
+from PyNucleus_base import INDEX
+from PyNucleus_base.tupleDict import tupleDictINDEX
 from copy import deepcopy