Change _libwarpx.py functions to LibWarpX class methods (#2696)

* Initial attempt at moving function into LibWarpX

* Bugfix - undefined structs in get_particle_theta

This is untested, identified by IDE/visual inspection only

* Unify newlines

* Delete libwarpx.clight (duplicated picmi.py)

See picmi.constants.c

* Change function descriptive comment to docstring

Also remove extra newline in __init__ beginning.

* Replace pywarpx/_libwarpx calls appropriately

* Fix atexit register of finalize function

* Reorder WarpX.py imports for pre-commit

* Use picmi.Simulation.extension = libwarpx alias

This will hopefully preempt user script changes if we restructure things
to allow multiple WarpX runs to be done within a single Python script.

* Fix libwarpx ref in WarpX.getProbLo/getProbHi

* Update get_particle_boundary_buffer doc reference

6 files changed, 1879 insertions, 1961 deletions

diff --git a/Python/pywarpx/PGroup.py b/Python/pywarpx/PGroup.py
index f19c0d079..0dcb664e5 100644
--- a/Python/pywarpx/PGroup.py
+++ b/Python/pywarpx/PGroup.py
@@ -6,7 +6,8 @@
 
 import numpy as np
 
-from . import _libwarpx
+from ._libwarpx import libwarpx
+from ._picmi import constants
 
 
 class PGroup(object):
@@ -65,7 +66,7 @@ class PGroup(object):
         pass
 
     def getnpid(self):
-        return _libwarpx.get_nattr()
+        return libwarpx.get_nattr()
     npid = property(getnpid)
 
     def getnps(self):
@@ -85,45 +86,45 @@ class PGroup(object):
     npmax = property(getnpmax)
 
     def getxp(self):
-        return _libwarpx.get_particle_x(self.ispecie, self.level)[self.igroup]
+        return libwarpx.get_particle_x(self.ispecie, self.level)[self.igroup]
     xp = property(getxp)
 
     def getyp(self):
-        return _libwarpx.get_particle_y(self.ispecie, self.level)[self.igroup]
+        return libwarpx.get_particle_y(self.ispecie, self.level)[self.igroup]
     yp = property(getyp)
 
     def getrp(self):
-        return _libwarpx.get_particle_r(self.ispecie, self.level)[self.igroup]
+        return libwarpx.get_particle_r(self.ispecie, self.level)[self.igroup]
     rp = property(getrp)
 
     def getzp(self):
-        return _libwarpx.get_particle_z(self.ispecie, self.level)[self.igroup]
+        return libwarpx.get_particle_z(self.ispecie, self.level)[self.igroup]
     zp = property(getzp)
 
     def getuxp(self):
-        return _libwarpx.get_particle_ux(self.ispecie, self.level)[self.igroup]
+        return libwarpx.get_particle_ux(self.ispecie, self.level)[self.igroup]
     uxp = property(getuxp)
 
     def getuyp(self):
-        return _libwarpx.get_particle_uy(self.ispecie, self.level)[self.igroup]
+        return libwarpx.get_particle_uy(self.ispecie, self.level)[self.igroup]
     uyp = property(getuyp)
 
     def getuzp(self):
-        return _libwarpx.get_particle_uz(self.ispecie, self.level)[self.igroup]
+        return libwarpx.get_particle_uz(self.ispecie, self.level)[self.igroup]
     uzp = property(getuzp)
 
     def getw(self):
-        return _libwarpx.get_particle_weight(self.ispecie, self.level)[self.igroup]
+        return libwarpx.get_particle_weight(self.ispecie, self.level)[self.igroup]
 
     def getpid(self, id):
-        pid = _libwarpx.get_particle_arrays(self.ispecie, id, self.level)[self.igroup]
+        pid = libwarpx.get_particle_arrays(self.ispecie, id, self.level)[self.igroup]
         return np.array([pid]).T
 
     def getgaminv(self):
         uxp = self.getuxp()
         uyp = self.getuyp()
         uzp = self.getuzp()
-        return np.sqrt(1. - (uxp**2 + uyp**2 + uzp**2)/_libwarpx.clight**2)
+        return np.sqrt(1. - (uxp**2 + uyp**2 + uzp**2)/constants.c**2)
     gaminv = property(getgaminv)
 
     def getex(self):
@@ -151,7 +152,7 @@ class PGroup(object):
     bz = property(getbz)
 
     def gettheta(self):
-        return _libwarpx.get_particle_theta(self.ispecie, self.level)[self.igroup]
+        return libwarpx.get_particle_theta(self.ispecie, self.level)[self.igroup]
     theta = property(gettheta)
 
 class PGroups(object):
@@ -160,7 +161,7 @@ class PGroups(object):
         self.level = level
 
     def setuppgroups(self):
-        xall = _libwarpx.get_particle_x(self.ispecie, self.level)
+        xall = libwarpx.get_particle_x(self.ispecie, self.level)
         self.ngroups = len(xall)
 
         self._pgroups = []
diff --git a/Python/pywarpx/WarpX.py b/Python/pywarpx/WarpX.py
index 8c4900fab..ece837816 100644
--- a/Python/pywarpx/WarpX.py
+++ b/Python/pywarpx/WarpX.py
@@ -5,7 +5,7 @@
 #
 # License: BSD-3-Clause-LBNL
 
-from . import Particles, _libwarpx
+from . import Particles
 from .Algo import algo
 from .Amr import amr
 from .Boundary import boundary
@@ -19,6 +19,7 @@ from .Langmuirwave import langmuirwave
 from .Lasers import lasers, lasers_list
 from .PSATD import psatd
 from .Particles import particles, particles_list
+from ._libwarpx import libwarpx
 
 
 class WarpX(Bucket):
@@ -75,19 +76,19 @@ class WarpX(Bucket):
 
     def init(self, mpi_comm=None):
         argv = ['warpx'] + self.create_argv_list()
-        _libwarpx.initialize(argv, mpi_comm=mpi_comm)
+        libwarpx.initialize(argv, mpi_comm=mpi_comm)
 
     def evolve(self, nsteps=-1):
-        _libwarpx.evolve(nsteps)
+        libwarpx.evolve(nsteps)
 
     def finalize(self, finalize_mpi=1):
-        _libwarpx.finalize(finalize_mpi)
+        libwarpx.finalize(finalize_mpi)
 
     def getProbLo(self, direction):
-        return _libwarpx.libwarpx.libwarpx_so.warpx_getProbLo(direction)
+        return libwarpx.libwarpx_so.warpx_getProbLo(direction)
 
     def getProbHi(self, direction):
-        return _libwarpx.libwarpx.libwarpx_so.warpx_getProbHi(direction)
+        return libwarpx.libwarpx_so.warpx_getProbHi(direction)
 
     def write_inputs(self, filename='inputs', **kw):
         argv = self.create_argv_list()
diff --git a/Python/pywarpx/__init__.py b/Python/pywarpx/__init__.py
index 8751b07ff..a41b50dc9 100644
--- a/Python/pywarpx/__init__.py
+++ b/Python/pywarpx/__init__.py
@@ -17,4 +17,4 @@ from .Lasers import lasers
 from .PSATD import psatd
 from .Particles import electrons, newspecies, particles, positrons, protons
 from .WarpX import warpx
-from ._libwarpx import *  # noqa
+from ._libwarpx import libwarpx
diff --git a/Python/pywarpx/_libwarpx.py b/Python/pywarpx/_libwarpx.py
index 84e13d1af..52bebf46a 100755
--- a/Python/pywarpx/_libwarpx.py
+++ b/Python/pywarpx/_libwarpx.py
@@ -40,16 +40,19 @@ _libc = ctypes.CDLL(_path_libc)
 _LP_c_int = ctypes.POINTER(ctypes.c_int)
 _LP_c_char = ctypes.POINTER(ctypes.c_char)
 
+
 class LibWarpX():
+
     """This class manages the warpx shared object, the library from the compiled C++ code.
     It will only load the library when it is referenced, and this can only be done after
     the geometry is defined so that the version of the library that is needed can be determined.
     Once loaded, all of the setting of function call interfaces is setup.
     """
-    def __init__(self):
 
+    def __init__(self):
         # Track whether amrex and warpx have been initialized
         self.initialized = False
+        atexit.register(self.finalize)
 
     def __getattr__(self, attribute):
         if attribute == 'libwarpx_so':
@@ -338,2374 +341,2281 @@ class LibWarpX():
 
         return boundary_num
 
+    @staticmethod
+    def _array1d_from_pointer(pointer, dtype, size):
+        '''
+        Function for converting a ctypes pointer to a numpy array
+        '''
+        if not pointer:
+            raise Exception(f'_array1d_from_pointer: pointer is a nullptr')
+        if sys.version_info.major >= 3:
+            # from where do I import these? this might only work for CPython...
+            #PyBuf_READ  = 0x100
+            PyBUF_WRITE = 0x200
+            buffer_from_memory = ctypes.pythonapi.PyMemoryView_FromMemory
+            buffer_from_memory.argtypes = (ctypes.c_void_p, ctypes.c_int, ctypes.c_int)
+            buffer_from_memory.restype = ctypes.py_object
+            buf = buffer_from_memory(pointer, dtype.itemsize*size, PyBUF_WRITE)
+        else:
+            buffer_from_memory = ctypes.pythonapi.PyBuffer_FromReadWriteMemory
+            buffer_from_memory.restype = ctypes.py_object
+            buf = buffer_from_memory(pointer, dtype.itemsize*size)
+        return np.frombuffer(buf, dtype=dtype, count=size)
 
-libwarpx = LibWarpX()
-
-# --- Is there a better way of handling constants?
-clight = 2.99792458e+8 # m/s
-
-# this is a function for converting a ctypes pointer to a numpy array
-def _array1d_from_pointer(pointer, dtype, size):
-    if not pointer:
-        raise Exception(f'_array1d_from_pointer: pointer is a nullptr')
-    if sys.version_info.major >= 3:
-        # from where do I import these? this might only work for CPython...
-        #PyBuf_READ  = 0x100
-        PyBUF_WRITE = 0x200
-        buffer_from_memory = ctypes.pythonapi.PyMemoryView_FromMemory
-        buffer_from_memory.argtypes = (ctypes.c_void_p, ctypes.c_int, ctypes.c_int)
-        buffer_from_memory.restype = ctypes.py_object
-        buf = buffer_from_memory(pointer, dtype.itemsize*size, PyBUF_WRITE)
-    else:
-        buffer_from_memory = ctypes.pythonapi.PyBuffer_FromReadWriteMemory
-        buffer_from_memory.restype = ctypes.py_object
-        buf = buffer_from_memory(pointer, dtype.itemsize*size)
-    return np.frombuffer(buf, dtype=dtype, count=size)
-
-def getNProcs():
-    '''
-
-    Get the number of processors
-
-    '''
-    return libwarpx.libwarpx_so.warpx_getNProcs()
-
-def getMyProc():
-    '''
-
-    Get the number of the processor
-
-    '''
-    return libwarpx.libwarpx_so.warpx_getMyProc()
-
-def get_nattr():
-    '''
-
-    Get the number of extra attributes.
-
-    '''
-    # --- The -3 is because the comps include the velocites
-    return libwarpx.libwarpx_so.warpx_nComps() - 3
-
-def get_nattr_species(species_name):
-    '''
-
-    Get the number of real attributes for the given species.
-
-    '''
-    return libwarpx.libwarpx_so.warpx_nCompsSpecies(
-        ctypes.c_char_p(species_name.encode('utf-8')))
-
-def amrex_init(argv, mpi_comm=None):
-    # --- Construct the ctype list of strings to pass in
-    argc = len(argv)
-    argvC = (_LP_c_char * (argc+1))()
-    for i, arg in enumerate(argv):
-        enc_arg = arg.encode('utf-8')
-        argvC[i] = ctypes.create_string_buffer(enc_arg)
-
-    if mpi_comm is None or MPI is None:
-        libwarpx.libwarpx_so.amrex_init(argc, argvC)
-    else:
-        comm_ptr = MPI._addressof(mpi_comm)
-        comm_val = _MPI_Comm_type.from_address(comm_ptr)
-        libwarpx.libwarpx_so.amrex_init_with_inited_mpi(argc, argvC, comm_val)
-
-def initialize(argv=None, mpi_comm=None):
-    '''
-
-    Initialize WarpX and AMReX. Must be called before
-    doing anything else.
-
-    '''
-    if argv is None:
-        argv = sys.argv
-    amrex_init(argv, mpi_comm)
-    libwarpx.libwarpx_so.warpx_ConvertLabParamsToBoost()
-    libwarpx.libwarpx_so.warpx_ReadBCParams()
-    if libwarpx.geometry_dim == 'rz':
-        libwarpx.libwarpx_so.warpx_CheckGriddingForRZSpectral()
-    libwarpx.libwarpx_so.warpx_init()
-
-    libwarpx.initialized = True
-
-
-@atexit.register
-def finalize(finalize_mpi=1):
-    '''
-
-    Call finalize for WarpX and AMReX. Must be called at
-    the end of your script.
-
-    '''
-    if libwarpx.initialized:
-        libwarpx.libwarpx_so.warpx_finalize()
-        libwarpx.libwarpx_so.amrex_finalize(finalize_mpi)
-
-def getistep(level=0):
-    '''
-
-    Get the current time step number for the specified level
-
-    Parameter
-    ---------
-
-    level : the refinement level to reference
-    '''
-    return libwarpx.libwarpx_so.warpx_getistep(level)
-
-def gett_new(level=0):
-    '''
-
-    Get the next time for the specified level
-
-    '''
-    return libwarpx.libwarpx_so.warpx_gett_new(level)
-
-def evolve(num_steps=-1):
-    '''
-
-    Evolve the simulation for num_steps steps. If num_steps=-1,
-    the simulation will be run until the end as specified in the
-    inputs file.
+    def getNProcs(self):
+        '''
 
-    Parameters
-    ----------
+        Get the number of processors
 
-    num_steps: int, the number of steps to take
+        '''
+        return self.libwarpx_so.warpx_getNProcs()
 
-    '''
+    def getMyProc(self):
+        '''
 
-    libwarpx.libwarpx_so.warpx_evolve(num_steps);
+        Get the number of the processor
 
+        '''
+        return self.libwarpx_so.warpx_getMyProc()
 
-def getProbLo(direction):
-    assert 0 <= direction < libwarpx.dim, 'Inappropriate direction specified'
-    return libwarpx.libwarpx_so.warpx_getProbLo(direction)
+    def get_nattr(self):
+        '''
 
+        Get the number of extra attributes.
 
-def getProbHi(direction):
-    assert 0 <= direction < libwarpx.dim, 'Inappropriate direction specified'
-    return libwarpx.libwarpx_so.warpx_getProbHi(direction)
+        '''
+        # --- The -3 is because the comps include the velocites
+        return self.libwarpx_so.warpx_nComps() - 3
 
+    def get_nattr_species(self, species_name):
+        '''
 
-def getCellSize(direction, level=0):
-    assert 0 <= direction < 3, 'Inappropriate direction specified'
-    assert 0 <= level and level <= libwarpx.libwarpx_so.warpx_finestLevel(), 'Inappropriate level specified'
-    return libwarpx.libwarpx_so.warpx_getCellSize(direction, level)
+        Get the number of real attributes for the given species.
 
+        '''
+        return self.libwarpx_so.warpx_nCompsSpecies(
+            ctypes.c_char_p(species_name.encode('utf-8')))
+
+    def amrex_init(self, argv, mpi_comm=None):
+        # --- Construct the ctype list of strings to pass in
+        argc = len(argv)
+        argvC = (_LP_c_char * (argc+1))()
+        for i, arg in enumerate(argv):
+            enc_arg = arg.encode('utf-8')
+            argvC[i] = ctypes.create_string_buffer(enc_arg)
+
+        if mpi_comm is None or MPI is None:
+            self.libwarpx_so.amrex_init(argc, argvC)
+        else:
+            comm_ptr = MPI._addressof(mpi_comm)
+            comm_val = _MPI_Comm_type.from_address(comm_ptr)
+            self.libwarpx_so.amrex_init_with_inited_mpi(argc, argvC, comm_val)
 
-#def get_sigma(direction):
-#    '''
-#
-#    Return the 'sigma' PML coefficients for the electric field
-#    in a given direction.
-#
-#    '''
-#
-#    size = ctypes.c_int(0)
-#    data = libwarpx.libwarpx_so.warpx_getPMLSigma(direction, ctypes.byref(size))
-#    arr = np.ctypeslib.as_array(data, (size.value,))
-#    arr.setflags(write=1)
-#    return arr
-#
-#
-#def get_sigma_star(direction):
-#    '''
-#
-#    Return the 'sigma*' PML coefficients for the magnetic field
-#    in the given direction.
-#
-#    '''
-#
-#    size = ctypes.c_int(0)
-#    data = libwarpx.libwarpx_so.warpx_getPMLSigmaStar(direction, ctypes.byref(size))
-#    arr = np.ctypeslib.as_array(data, (size.value,))
-#    arr.setflags(write=1)
-#    return arr
-#
-#
-#def compute_pml_factors(lev, dt):
-#    '''
-#
-#    This recomputes the PML coefficients for a given level, using the
-#    time step dt. This needs to be called after modifying the coefficients
-#    from Python.
-#
-#    '''
-#
-#    libwarpx.libwarpx_so.warpx_ComputePMLFactors(lev, dt)
-
-def add_particles(species_name, x=None, y=None, z=None, ux=None, uy=None, uz=None, w=None,
-                  unique_particles=True, **kwargs):
-    '''
-
-    A function for adding particles to the WarpX simulation.
-
-    Parameters
-    ----------
-
-    species_name     : the species to add the particle to
-    x, y, z          : arrays or scalars of the particle positions (default = 0.)
-    ux, uy, uz       : arrays or scalars of the particle momenta (default = 0.)
-    w                : array or scalar of particle weights (default = 0.)
-    unique_particles : whether the particles are unique or duplicated on
-                       several processes. (default = True)
-    kwargs           : dictionary containing an entry for all the extra particle
-                       attribute arrays. If an attribute is not given it will be
-                       set to 0.
-
-    '''
-
-    # --- Get length of arrays, set to one for scalars
-    lenx = np.size(x)
-    leny = np.size(y)
-    lenz = np.size(z)
-    lenux = np.size(ux)
-    lenuy = np.size(uy)
-    lenuz = np.size(uz)
-    lenw = np.size(w)
-
-    # --- Find the max length of the parameters supplied
-    maxlen = 0
-    if x is not None:
-        maxlen = max(maxlen, lenx)
-    if y is not None:
-        maxlen = max(maxlen, leny)
-    if z is not None:
-        maxlen = max(maxlen, lenz)
-    if ux is not None:
-        maxlen = max(maxlen, lenux)
-    if uy is not None:
-        maxlen = max(maxlen, lenuy)
-    if uz is not None:
-        maxlen = max(maxlen, lenuz)
-    if w is not None:
-        maxlen = max(maxlen, lenw)
-
-    # --- Make sure that the lengths of the input parameters are consistent
-    assert x is None or lenx==maxlen or lenx==1, "Length of x doesn't match len of others"
-    assert y is None or leny==maxlen or leny==1, "Length of y doesn't match len of others"
-    assert z is None or lenz==maxlen or lenz==1, "Length of z doesn't match len of others"
-    assert ux is None or lenux==maxlen or lenux==1, "Length of ux doesn't match len of others"
-    assert uy is None or lenuy==maxlen or lenuy==1, "Length of uy doesn't match len of others"
-    assert uz is None or lenuz==maxlen or lenuz==1, "Length of uz doesn't match len of others"
-    assert w is None or lenw==maxlen or lenw==1, "Length of w doesn't match len of others"
-    for key, val in kwargs.items():
-        assert np.size(val)==1 or len(val)==maxlen, f"Length of {key} doesn't match len of others"
-
-    # --- If the length of the input is zero, then quietly return
-    # --- This is not an error - it just means that no particles are to be injected.
-    if maxlen == 0:
-        return
-
-    # --- Broadcast scalars into appropriate length arrays
-    # --- If the parameter was not supplied, use the default value
-    if lenx == 1:
-        x = np.full(maxlen, (x or 0.), float)
-    if leny == 1:
-        y = np.full(maxlen, (y or 0.), float)
-    if lenz == 1:
-        z = np.full(maxlen, (z or 0.), float)
-    if lenux == 1:
-        ux = np.full(maxlen, (ux or 0.), float)
-    if lenuy == 1:
-        uy = np.full(maxlen, (uy or 0.), float)
-    if lenuz == 1:
-        uz = np.full(maxlen, (uz or 0.), float)
-    if lenw == 1:
-        w = np.full(maxlen, (w or 0.), float)
-    for key, val in kwargs.items():
-        if np.size(val) == 1:
-            kwargs[key] = np.full(maxlen, val, float)
-
-    # --- The -3 is because the comps include the velocites
-    nattr = get_nattr_species(species_name) - 3
-    attr = np.zeros((maxlen, nattr))
-    attr[:,0] = w
-
-    for key, vals in kwargs.items():
-        # --- The -3 is because components 1 to 3 are velocities
-        attr[:,get_particle_comp_index(species_name, key)-3] = vals
-
-    libwarpx.libwarpx_so.warpx_addNParticles(
-        ctypes.c_char_p(species_name.encode('utf-8')), x.size,
-        x, y, z, ux, uy, uz, nattr, attr, unique_particles
-    )
-
-
-def get_particle_count(species_name):
-    '''
-
-    This returns the number of particles of the specified species in the
-    simulation.
-
-    Parameters
-    ----------
-
-        species_name : the species name that the number will be returned for
-
-    Returns
-    -------
-
-        An integer count of the number of particles
-
-    '''
-    return libwarpx.libwarpx_so.warpx_getNumParticles(
-        ctypes.c_char_p(species_name.encode('utf-8'))
-    )
-
-
-def get_particle_structs(species_name, level):
-    '''
-
-    This returns a list of numpy arrays containing the particle struct data
-    on each tile for this process. The particle data is represented as a structured
-    numpy array and contains the particle 'x', 'y', 'z', 'id', and 'cpu'.
-
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
-
-    Parameters
-    ----------
-
-        species_name : the species name that the data will be returned for
-
-    Returns
-    -------
-
-        A List of numpy arrays.
-
-    '''
-
-    particles_per_tile = _LP_c_int()
-    num_tiles = ctypes.c_int(0)
-    data = libwarpx.libwarpx_so.warpx_getParticleStructs(
-        ctypes.c_char_p(species_name.encode('utf-8')), level,
-        ctypes.byref(num_tiles), ctypes.byref(particles_per_tile)
-    )
-
-    particle_data = []
-    for i in range(num_tiles.value):
-        if particles_per_tile[i] == 0:
-            continue
-        arr = _array1d_from_pointer(data[i], libwarpx._p_dtype, particles_per_tile[i])
-        particle_data.append(arr)
+    def initialize(self, argv=None, mpi_comm=None):
+        '''
 
-    _libc.free(particles_per_tile)
-    _libc.free(data)
-    return particle_data
+        Initialize WarpX and AMReX. Must be called before
+        doing anything else.
 
+        '''
+        if argv is None:
+            argv = sys.argv
+        self.amrex_init(argv, mpi_comm)
+        self.libwarpx_so.warpx_ConvertLabParamsToBoost()
+        self.libwarpx_so.warpx_ReadBCParams()
+        if self.geometry_dim == 'rz':
+            self.libwarpx_so.warpx_CheckGriddingForRZSpectral()
+        self.libwarpx_so.warpx_init()
+
+        self.initialized = True
+
+    def finalize(self, finalize_mpi=1):
+        '''
 
-def get_particle_arrays(species_name, comp_name, level):
-    '''
+        Call finalize for WarpX and AMReX. Registered to run at program exit.
 
-    This returns a list of numpy arrays containing the particle array data
-    on each tile for this process.
+        '''
+        if self.initialized:
+            self.libwarpx_so.warpx_finalize()
+            self.libwarpx_so.amrex_finalize(finalize_mpi)
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+    def getistep(self, level=0):
+        '''
 
-    Parameters
-    ----------
+        Get the current time step number for the specified level
 
-        species_name   : the species name that the data will be returned for
-        comp_name      : the component of the array data that will be returned.
+        Parameter
+        ---------
 
-    Returns
-    -------
+        level : the refinement level to reference
+        '''
+        return self.libwarpx_so.warpx_getistep(level)
 
-        A List of numpy arrays.
+    def gett_new(self, level=0):
+        '''
 
-    '''
+        Get the next time for the specified level
 
-    particles_per_tile = _LP_c_int()
-    num_tiles = ctypes.c_int(0)
-    data = libwarpx.libwarpx_so.warpx_getParticleArrays(
-        ctypes.c_char_p(species_name.encode('utf-8')),
-        ctypes.c_char_p(comp_name.encode('utf-8')),
-        level, ctypes.byref(num_tiles), ctypes.byref(particles_per_tile)
-    )
+        '''
+        return self.libwarpx_so.warpx_gett_new(level)
 
-    particle_data = []
-    for i in range(num_tiles.value):
-        if particles_per_tile[i] == 0:
-            continue
-        if not data[i]:
-            raise Exception(f'get_particle_arrays: data[i] for i={i} was not initialized')
-        arr = np.ctypeslib.as_array(data[i], (particles_per_tile[i],))
-        try:
-            # This fails on some versions of numpy
-            arr.setflags(write=1)
-        except ValueError:
-            pass
-        particle_data.append(arr)
+    def evolve(self, num_steps=-1):
+        '''
 
-    _libc.free(particles_per_tile)
-    _libc.free(data)
-    return particle_data
+        Evolve the simulation for num_steps steps. If num_steps=-1,
+        the simulation will be run until the end as specified in the
+        inputs file.
 
+        Parameters
+        ----------
 
-def get_particle_x(species_name, level=0):
-    '''
+        num_steps: int, the number of steps to take
 
-    Return a list of numpy arrays containing the particle 'x'
-    positions on each tile.
+        '''
 
-    '''
-    structs = get_particle_structs(species_name, level)
-    if libwarpx.geometry_dim == '3d' or libwarpx.geometry_dim == '2d':
-        return [struct['x'] for struct in structs]
-    elif libwarpx.geometry_dim == 'rz':
-        return [struct['x']*np.cos(theta) for struct, theta in zip(structs, get_particle_theta(species_name))]
+        self.libwarpx_so.warpx_evolve(num_steps);
+
+    def getProbLo(self, direction):
+        assert 0 <= direction < self.dim, 'Inappropriate direction specified'
+        return self.libwarpx_so.warpx_getProbLo(direction)
+
+    def getProbHi(self, direction):
+        assert 0 <= direction < self.dim, 'Inappropriate direction specified'
+        return self.libwarpx_so.warpx_getProbHi(direction)
+
+    def getCellSize(self, direction, level=0):
+        assert 0 <= direction < 3, 'Inappropriate direction specified'
+        assert 0 <= level and level <= self.libwarpx_so.warpx_finestLevel(), 'Inappropriate level specified'
+        return self.libwarpx_so.warpx_getCellSize(direction, level)
+
+    #def get_sigma(self, direction):
+    #    '''
+    #
+    #    Return the 'sigma' PML coefficients for the electric field
+    #    in a given direction.
+    #
+    #    '''
+    #
+    #    size = ctypes.c_int(0)
+    #    data = self.libwarpx_so.warpx_getPMLSigma(direction, ctypes.byref(size))
+    #    arr = np.ctypeslib.as_array(data, (size.value,))
+    #    arr.setflags(write=1)
+    #    return arr
+    #
+    #
+    #def get_sigma_star(self, direction):
+    #    '''
+    #
+    #    Return the 'sigma*' PML coefficients for the magnetic field
+    #    in the given direction.
+    #
+    #    '''
+    #
+    #    size = ctypes.c_int(0)
+    #    data = self.libwarpx_so.warpx_getPMLSigmaStar(direction, ctypes.byref(size))
+    #    arr = np.ctypeslib.as_array(data, (size.value,))
+    #    arr.setflags(write=1)
+    #    return arr
+    #
+    #
+    #def compute_pml_factors(self, lev, dt):
+    #    '''
+    #
+    #    This recomputes the PML coefficients for a given level, using the
+    #    time step dt. This needs to be called after modifying the coefficients
+    #    from Python.
+    #
+    #    '''
+    #
+    #    self.libwarpx_so.warpx_ComputePMLFactors(lev, dt)
+
+    def add_particles(self, species_name, x=None, y=None, z=None, ux=None, uy=None, uz=None, w=None,
+                    unique_particles=True, **kwargs):
+        '''
 
+        A function for adding particles to the WarpX simulation.
 
-def get_particle_y(species_name, level=0):
-    '''
+        Parameters
+        ----------
 
-    Return a list of numpy arrays containing the particle 'y'
-    positions on each tile.
+        species_name     : the species to add the particle to
+        x, y, z          : arrays or scalars of the particle positions (default = 0.)
+        ux, uy, uz       : arrays or scalars of the particle momenta (default = 0.)
+        w                : array or scalar of particle weights (default = 0.)
+        unique_particles : whether the particles are unique or duplicated on
+                        several processes. (default = True)
+        kwargs           : dictionary containing an entry for all the extra particle
+                        attribute arrays. If an attribute is not given it will be
+                        set to 0.
 
-    '''
-    structs = get_particle_structs(species_name, level)
-    if libwarpx.geometry_dim == '3d' or libwarpx.geometry_dim == '2d':
-        return [struct['y'] for struct in structs]
-    elif libwarpx.geometry_dim == 'rz':
-        return [struct['x']*np.sin(theta) for struct, theta in zip(structs, get_particle_theta(species_name))]
+        '''
 
+        # --- Get length of arrays, set to one for scalars
+        lenx = np.size(x)
+        leny = np.size(y)
+        lenz = np.size(z)
+        lenux = np.size(ux)
+        lenuy = np.size(uy)
+        lenuz = np.size(uz)
+        lenw = np.size(w)
+
+        # --- Find the max length of the parameters supplied
+        maxlen = 0
+        if x is not None:
+            maxlen = max(maxlen, lenx)
+        if y is not None:
+            maxlen = max(maxlen, leny)
+        if z is not None:
+            maxlen = max(maxlen, lenz)
+        if ux is not None:
+            maxlen = max(maxlen, lenux)
+        if uy is not None:
+            maxlen = max(maxlen, lenuy)
+        if uz is not None:
+            maxlen = max(maxlen, lenuz)
+        if w is not None:
+            maxlen = max(maxlen, lenw)
+
+        # --- Make sure that the lengths of the input parameters are consistent
+        assert x is None or lenx==maxlen or lenx==1, "Length of x doesn't match len of others"
+        assert y is None or leny==maxlen or leny==1, "Length of y doesn't match len of others"
+        assert z is None or lenz==maxlen or lenz==1, "Length of z doesn't match len of others"
+        assert ux is None or lenux==maxlen or lenux==1, "Length of ux doesn't match len of others"
+        assert uy is None or lenuy==maxlen or lenuy==1, "Length of uy doesn't match len of others"
+        assert uz is None or lenuz==maxlen or lenuz==1, "Length of uz doesn't match len of others"
+        assert w is None or lenw==maxlen or lenw==1, "Length of w doesn't match len of others"
+        for key, val in kwargs.items():
+            assert np.size(val)==1 or len(val)==maxlen, f"Length of {key} doesn't match len of others"
+
+        # --- If the length of the input is zero, then quietly return
+        # --- This is not an error - it just means that no particles are to be injected.
+        if maxlen == 0:
+            return
+
+        # --- Broadcast scalars into appropriate length arrays
+        # --- If the parameter was not supplied, use the default value
+        if lenx == 1:
+            x = np.full(maxlen, (x or 0.), float)
+        if leny == 1:
+            y = np.full(maxlen, (y or 0.), float)
+        if lenz == 1:
+            z = np.full(maxlen, (z or 0.), float)
+        if lenux == 1:
+            ux = np.full(maxlen, (ux or 0.), float)
+        if lenuy == 1:
+            uy = np.full(maxlen, (uy or 0.), float)
+        if lenuz == 1:
+            uz = np.full(maxlen, (uz or 0.), float)
+        if lenw == 1:
+            w = np.full(maxlen, (w or 0.), float)
+        for key, val in kwargs.items():
+            if np.size(val) == 1:
+                kwargs[key] = np.full(maxlen, val, float)
+
+        # --- The -3 is because the comps include the velocites
+        nattr = self.get_nattr_species(species_name) - 3
+        attr = np.zeros((maxlen, nattr))
+        attr[:,0] = w
+
+        for key, vals in kwargs.items():
+            # --- The -3 is because components 1 to 3 are velocities
+            attr[:,self.get_particle_comp_index(species_name, key)-3] = vals
+
+        self.libwarpx_so.warpx_addNParticles(
+            ctypes.c_char_p(species_name.encode('utf-8')), x.size,
+            x, y, z, ux, uy, uz, nattr, attr, unique_particles
+        )
 
-def get_particle_r(species_name, level=0):
-    '''
+    def get_particle_count(self, species_name):
+        '''
 
-    Return a list of numpy arrays containing the particle 'r'
-    positions on each tile.
+        This returns the number of particles of the specified species in the
+        simulation.
 
-    '''
-    structs = get_particle_structs(species_name, level)
-    if libwarpx.geometry_dim == 'rz':
-        return [struct['x'] for struct in structs]
-    elif libwarpx.geometry_dim == '3d':
-        return [np.sqrt(struct['x']**2 + struct['y']**2) for struct in structs]
-    elif libwarpx.geometry_dim == '2d':
-        raise Exception('get_particle_r: There is no r coordinate with 2D Cartesian')
+        Parameters
+        ----------
 
+            species_name : the species name that the number will be returned for
 
-def get_particle_z(species_name, level=0):
-    '''
+        Returns
+        -------
 
-    Return a list of numpy arrays containing the particle 'z'
-    positions on each tile.
+            An integer count of the number of particles
 
-    '''
-    structs = get_particle_structs(species_name, level)
-    if libwarpx.geometry_dim == '3d':
-        return [struct['z'] for struct in structs]
-    elif libwarpx.geometry_dim == 'rz' or libwarpx.geometry_dim == '2d':
-        return [struct['y'] for struct in structs]
+        '''
+        return self.libwarpx_so.warpx_getNumParticles(
+            ctypes.c_char_p(species_name.encode('utf-8'))
+        )
 
+    def get_particle_structs(self, species_name, level):
+        '''
 
-def get_particle_id(species_name, level=0):
-    '''
+        This returns a list of numpy arrays containing the particle struct data
+        on each tile for this process. The particle data is represented as a structured
+        numpy array and contains the particle 'x', 'y', 'z', 'id', and 'cpu'.
 
-    Return a list of numpy arrays containing the particle 'id'
-    positions on each tile.
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-    '''
-    structs = get_particle_structs(species_name, level)
-    return [struct['id'] for struct in structs]
+        Parameters
+        ----------
 
+            species_name : the species name that the data will be returned for
 
-def get_particle_cpu(species_name, level=0):
-    '''
+        Returns
+        -------
 
-    Return a list of numpy arrays containing the particle 'cpu'
-    positions on each tile.
+            A List of numpy arrays.
 
-    '''
-    structs = get_particle_structs(species_name, level)
-    return [struct['cpu'] for struct in structs]
+        '''
 
+        particles_per_tile = _LP_c_int()
+        num_tiles = ctypes.c_int(0)
+        data = self.libwarpx_so.warpx_getParticleStructs(
+            ctypes.c_char_p(species_name.encode('utf-8')), level,
+            ctypes.byref(num_tiles), ctypes.byref(particles_per_tile)
+        )
 
-def get_particle_weight(species_name, level=0):
-    '''
+        particle_data = []
+        for i in range(num_tiles.value):
+            if particles_per_tile[i] == 0:
+                continue
+            arr = self._array1d_from_pointer(data[i], self._p_dtype, particles_per_tile[i])
+            particle_data.append(arr)
 
-    Return a list of numpy arrays containing the particle
-    weight on each tile.
+        _libc.free(particles_per_tile)
+        _libc.free(data)
+        return particle_data
 
-    '''
+    def get_particle_arrays(self, species_name, comp_name, level):
+        '''
 
-    return get_particle_arrays(species_name, 'w', level)
+        This returns a list of numpy arrays containing the particle array data
+        on each tile for this process.
 
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-def get_particle_ux(species_name, level=0):
-    '''
+        Parameters
+        ----------
 
-    Return a list of numpy arrays containing the particle
-    x momentum on each tile.
+            species_name   : the species name that the data will be returned for
+            comp_name      : the component of the array data that will be returned.
 
-    '''
+        Returns
+        -------
 
-    return get_particle_arrays(species_name, 'ux', level)
+            A List of numpy arrays.
 
+        '''
 
-def get_particle_uy(species_name, level=0):
-    '''
+        particles_per_tile = _LP_c_int()
+        num_tiles = ctypes.c_int(0)
+        data = self.libwarpx_so.warpx_getParticleArrays(
+            ctypes.c_char_p(species_name.encode('utf-8')),
+            ctypes.c_char_p(comp_name.encode('utf-8')),
+            level, ctypes.byref(num_tiles), ctypes.byref(particles_per_tile)
+        )
 
-    Return a list of numpy arrays containing the particle
-    y momentum on each tile.
+        particle_data = []
+        for i in range(num_tiles.value):
+            if particles_per_tile[i] == 0:
+                continue
+            if not data[i]:
+                raise Exception(f'get_particle_arrays: data[i] for i={i} was not initialized')
+            arr = np.ctypeslib.as_array(data[i], (particles_per_tile[i],))
+            try:
+                # This fails on some versions of numpy
+                arr.setflags(write=1)
+            except ValueError:
+                pass
+            particle_data.append(arr)
+
+        _libc.free(particles_per_tile)
+        _libc.free(data)
+        return particle_data
+
+    def get_particle_x(self, species_name, level=0):
+        '''
 
-    '''
+        Return a list of numpy arrays containing the particle 'x'
+        positions on each tile.
 
-    return get_particle_arrays(species_name, 'uy', level)
+        '''
+        structs = self.get_particle_structs(species_name, level)
+        if self.geometry_dim == '3d' or self.geometry_dim == '2d':
+            return [struct['x'] for struct in structs]
+        elif self.geometry_dim == 'rz':
+            return [struct['x']*np.cos(theta) for struct, theta in zip(structs, self.get_particle_theta(species_name))]
 
+    def get_particle_y(self, species_name, level=0):
+        '''
 
-def get_particle_uz(species_name, level=0):
-    '''
+        Return a list of numpy arrays containing the particle 'y'
+        positions on each tile.
 
-    Return a list of numpy arrays containing the particle
-    z momentum on each tile.
+        '''
+        structs = self.get_particle_structs(species_name, level)
+        if self.geometry_dim == '3d' or self.geometry_dim == '2d':
+            return [struct['y'] for struct in structs]
+        elif self.geometry_dim == 'rz':
+            return [struct['x']*np.sin(theta) for struct, theta in zip(structs, self.get_particle_theta(species_name))]
 
-    '''
+    def get_particle_r(self, species_name, level=0):
+        '''
 
-    return get_particle_arrays(species_name, 'uz', level)
+        Return a list of numpy arrays containing the particle 'r'
+        positions on each tile.
 
+        '''
+        structs = self.get_particle_structs(species_name, level)
+        if self.geometry_dim == 'rz':
+            return [struct['x'] for struct in structs]
+        elif self.geometry_dim == '3d':
+            return [np.sqrt(struct['x']**2 + struct['y']**2) for struct in structs]
+        elif self.geometry_dim == '2d':
+            raise Exception('get_particle_r: There is no r coordinate with 2D Cartesian')
 
-def get_particle_theta(species_name, level=0):
-    '''
+    def get_particle_z(self, species_name, level=0):
+        '''
 
-    Return a list of numpy arrays containing the particle
-    theta on each tile.
+        Return a list of numpy arrays containing the particle 'z'
+        positions on each tile.
 
-    '''
+        '''
+        structs = self.get_particle_structs(species_name, level)
+        if self.geometry_dim == '3d':
+            return [struct['z'] for struct in structs]
+        elif self.geometry_dim == 'rz' or self.geometry_dim == '2d':
+            return [struct['y'] for struct in structs]
 
-    if libwarpx.geometry_dim == 'rz':
-        return get_particle_arrays(species_name, 'theta', level)
-    elif libwarpx.geometry_dim == '3d':
-        return [np.arctan2(struct['y'], struct['x']) for struct in structs]
-    elif libwarpx.geometry_dim == '2d':
-        raise Exception('get_particle_r: There is no theta coordinate with 2D Cartesian')
+    def get_particle_id(self, species_name, level=0):
+        '''
 
+        Return a list of numpy arrays containing the particle 'id'
+        positions on each tile.
 
-def get_particle_comp_index(species_name, pid_name):
-    '''
+        '''
+        structs = self.get_particle_structs(species_name, level)
+        return [struct['id'] for struct in structs]
 
-    Get the component index for a given particle attribute. This is useful
-    to get the corrent ordering of attributes when adding new particles using
-    `add_particles()`.
+    def get_particle_cpu(self, species_name, level=0):
+        '''
 
-    Parameters
-    ----------
+        Return a list of numpy arrays containing the particle 'cpu'
+        positions on each tile.
 
-        species_name   : the species name that the data will be returned for
-        pid_name       : string that is used to identify the new component
+        '''
+        structs = self.get_particle_structs(species_name, level)
+        return [struct['cpu'] for struct in structs]
 
-    Returns
-    -------
+    def get_particle_weight(self, species_name, level=0):
+        '''
 
-        Integer corresponding to the index of the requested attribute
+        Return a list of numpy arrays containing the particle
+        weight on each tile.
 
-    '''
-    return libwarpx.libwarpx_so.warpx_getParticleCompIndex(
-        ctypes.c_char_p(species_name.encode('utf-8')),
-        ctypes.c_char_p(pid_name.encode('utf-8'))
-    )
+        '''
 
+        return self.get_particle_arrays(species_name, 'w', level)
 
-def add_real_comp(species_name, pid_name, comm=True):
-    '''
+    def get_particle_ux(self, species_name, level=0):
+        '''
 
-    Add a real component to the particle data array.
+        Return a list of numpy arrays containing the particle
+        x momentum on each tile.
 
-    Parameters
-    ----------
+        '''
 
-        species_name   : the species name for which the new component will be added
-        pid_name       : string that is used to identify the new component
-        comm           : should the component be communicated
+        return self.get_particle_arrays(species_name, 'ux', level)
 
-    '''
-    libwarpx.libwarpx_so.warpx_addRealComp(
-        ctypes.c_char_p(species_name.encode('utf-8')),
-        ctypes.c_char_p(pid_name.encode('utf-8')), comm
-    )
+    def get_particle_uy(self, species_name, level=0):
+        '''
 
+        Return a list of numpy arrays containing the particle
+        y momentum on each tile.
 
-def get_particle_boundary_buffer_size(species_name, boundary):
-    '''
+        '''
 
-    This returns the number of particles that have been scraped so far in the simulation
-    from the specified boundary and of the specified species.
+        return self.get_particle_arrays(species_name, 'uy', level)
 
-    Parameters
-    ----------
+    def get_particle_uz(self, species_name, level=0):
+        '''
 
-        species_name   : return the number of scraped particles of this species
-        boundary       : the boundary from which to get the scraped particle data.
-                         In the form x/y/z_hi/lo
+        Return a list of numpy arrays containing the particle
+        z momentum on each tile.
 
-    Returns
-    -------
+        '''
 
-        The number of particles scraped so far from a boundary and of a species.
+        return self.get_particle_arrays(species_name, 'uz', level)
 
-    '''
-    return libwarpx.libwarpx_so.warpx_getParticleBoundaryBufferSize(
-        ctypes.c_char_p(species_name.encode('utf-8')),
-        libwarpx.get_boundary_number(boundary)
-    )
+    def get_particle_theta(self, species_name, level=0):
+        '''
 
+        Return a list of numpy arrays containing the particle
+        theta on each tile.
 
-def get_particle_boundary_buffer_structs(species_name, boundary, level):
-    '''
+        '''
 
-    This returns a list of numpy arrays containing the particle struct data
-    for a species that has been scraped by a specific simulation boundary. The
-    particle data is represented as a structured numpy array and contains the
-    particle 'x', 'y', 'z', 'id', and 'cpu'.
+        if self.geometry_dim == 'rz':
+            return self.get_particle_arrays(species_name, 'theta', level)
+        elif self.geometry_dim == '3d':
+            structs = self.get_particle_structs(species_name, level)
+            return [np.arctan2(struct['y'], struct['x']) for struct in structs]
+        elif self.geometry_dim == '2d':
+            raise Exception('get_particle_r: There is no theta coordinate with 2D Cartesian')
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+    def get_particle_comp_index(self, species_name, pid_name):
+        '''
 
-    Parameters
-    ----------
+        Get the component index for a given particle attribute. This is useful
+        to get the corrent ordering of attributes when adding new particles using
+        `add_particles()`.
 
-        species_name : the species name that the data will be returned for
-        boundary     : the boundary from which to get the scraped particle data.
-                       In the form x/y/z_hi/lo or eb.
-        level        : Which AMR level to retrieve scraped particle data from.
+        Parameters
+        ----------
 
-    Returns
-    -------
+            species_name   : the species name that the data will be returned for
+            pid_name       : string that is used to identify the new component
 
-        A List of numpy arrays.
+        Returns
+        -------
 
-    '''
+            Integer corresponding to the index of the requested attribute
 
-    particles_per_tile = _LP_c_int()
-    num_tiles = ctypes.c_int(0)
-    data = libwarpx.libwarpx_so.warpx_getParticleBoundaryBufferStructs(
+        '''
+        return self.libwarpx_so.warpx_getParticleCompIndex(
             ctypes.c_char_p(species_name.encode('utf-8')),
-            libwarpx.get_boundary_number(boundary), level,
-            ctypes.byref(num_tiles), ctypes.byref(particles_per_tile)
-    )
-
-    particle_data = []
-    for i in range(num_tiles.value):
-        if particles_per_tile[i] == 0:
-            continue
-        arr = _array1d_from_pointer(data[i], libwarpx._p_dtype, particles_per_tile[i])
-        particle_data.append(arr)
-
-    _libc.free(particles_per_tile)
-    _libc.free(data)
-    return particle_data
-
-
-def get_particle_boundary_buffer(species_name, boundary, comp_name, level):
-    '''
-
-    This returns a list of numpy arrays containing the particle array data
-    for a species that has been scraped by a specific simulation boundary.
-
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+            ctypes.c_char_p(pid_name.encode('utf-8'))
+        )
 
-    Parameters
-    ----------
+    def add_real_comp(self, species_name, pid_name, comm=True):
+        '''
 
-        species_name   : the species name that the data will be returned for.
-        boundary       : the boundary from which to get the scraped particle data.
-                         In the form x/y/z_hi/lo or eb.
-        comp_name      : the component of the array data that will be returned.
-                         If "step_scraped" the special attribute holding the
-                         timestep at which a particle was scraped will be
-                         returned.
-        level          : Which AMR level to retrieve scraped particle data from.
+        Add a real component to the particle data array.
 
-    Returns
-    -------
+        Parameters
+        ----------
 
-        A List of numpy arrays.
+            species_name   : the species name for which the new component will be added
+            pid_name       : string that is used to identify the new component
+            comm           : should the component be communicated
 
-    '''
-    particles_per_tile = _LP_c_int()
-    num_tiles = ctypes.c_int(0)
-    if comp_name == 'step_scraped':
-        data = libwarpx.libwarpx_so.warpx_getParticleBoundaryBufferScrapedSteps(
-            ctypes.c_char_p(species_name.encode('utf-8')),
-            libwarpx.get_boundary_number(boundary), level,
-            ctypes.byref(num_tiles), ctypes.byref(particles_per_tile)
-        )
-    else:
-        data = libwarpx.libwarpx_so.warpx_getParticleBoundaryBuffer(
+        '''
+        self.libwarpx_so.warpx_addRealComp(
             ctypes.c_char_p(species_name.encode('utf-8')),
-            libwarpx.get_boundary_number(boundary), level,
-            ctypes.byref(num_tiles), ctypes.byref(particles_per_tile),
-            ctypes.c_char_p(comp_name.encode('utf-8'))
+            ctypes.c_char_p(pid_name.encode('utf-8')), comm
         )
 
-    particle_data = []
-    for i in range(num_tiles.value):
-        if particles_per_tile[i] == 0:
-            continue
-        if not data[i]:
-            raise Exception(f'get_particle_arrays: data[i] for i={i} was not initialized')
-        arr = np.ctypeslib.as_array(data[i], (particles_per_tile[i],))
-        try:
-            # This fails on some versions of numpy
-            arr.setflags(write=1)
-        except ValueError:
-            pass
-        particle_data.append(arr)
-
-    _libc.free(particles_per_tile)
-    _libc.free(data)
-    return particle_data
-
-
-def clearParticleBoundaryBuffer():
-    '''
-
-    Clear the buffer that holds the particles lost at the boundaries.
-
-    '''
-    libwarpx.libwarpx_so.warpx_clearParticleBoundaryBuffer()
-
-
-def _get_mesh_field_list(warpx_func, level, direction, include_ghosts):
-    """
-     Generic routine to fetch the list of field data arrays.
-    """
-    shapes = _LP_c_int()
-    size = ctypes.c_int(0)
-    ncomps = ctypes.c_int(0)
-    ngrowvect = _LP_c_int()
-    if direction is None:
-        data = warpx_func(level,
-                          ctypes.byref(size), ctypes.byref(ncomps),
-                          ctypes.byref(ngrowvect), ctypes.byref(shapes))
-    else:
-        data = warpx_func(level, direction,
-                          ctypes.byref(size), ctypes.byref(ncomps),
-                          ctypes.byref(ngrowvect), ctypes.byref(shapes))
-    if not data:
-        raise Exception('object was not initialized')
-
-    ngvect = [ngrowvect[i] for i in range(libwarpx.dim)]
-    grid_data = []
-    shapesize = libwarpx.dim
-    if ncomps.value > 1:
-        shapesize += 1
-    for i in range(size.value):
-        shape = tuple([shapes[shapesize*i + d] for d in range(shapesize)])
-        # --- The data is stored in Fortran order, hence shape is reversed and a transpose is taken.
-        if shape[::-1] == 0:
-            continue
-        if not data[i]:
-            raise Exception(f'get_particle_arrays: data[i] for i={i} was not initialized')
-        arr = np.ctypeslib.as_array(data[i], shape[::-1]).T
-        try:
-            # This fails on some versions of numpy
-            arr.setflags(write=1)
-        except ValueError:
-            pass
-        if include_ghosts:
-            grid_data.append(arr)
-        else:
-            grid_data.append(arr[tuple([slice(ngvect[d], -ngvect[d]) for d in range(libwarpx.dim)])])
-
-    _libc.free(shapes)
-    _libc.free(data)
-    return grid_data
-
-
-def get_mesh_electric_field(level, direction, include_ghosts=True):
-    '''
-
-    This returns a list of numpy arrays containing the mesh electric field
-    data on each grid for this process.
-
-    This version is for the full "auxiliary" solution on the given level.
-
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
-
-    Parameters
-    ----------
-
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
-
-    Returns
-    -------
-
-        A List of numpy arrays.
-
-    '''
-
-    return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getEfield, level, direction, include_ghosts)
-
-
-def get_mesh_electric_field_cp(level, direction, include_ghosts=True):
-    '''
-
-    This returns a list of numpy arrays containing the mesh electric field
-    data on each grid for this process. This version returns the field on
-    the coarse patch for the given level.
-
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
-
-    Parameters
-    ----------
-
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
-
-    Returns
-    -------
-
-        A List of numpy arrays.
-
-    '''
-
-    return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getEfieldCP, level, direction, include_ghosts)
-
-
-def get_mesh_electric_field_fp(level, direction, include_ghosts=True):
-    '''
-
-    This returns a list of numpy arrays containing the mesh electric field
-    data on each grid for this process. This version returns the field on
-    the fine patch for the given level.
-
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+    def get_particle_boundary_buffer_size(self, species_name, boundary):
+        '''
 
-    Parameters
-    ----------
+        This returns the number of particles that have been scraped so far in the simulation
+        from the specified boundary and of the specified species.
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+        Parameters
+        ----------
 
-    Returns
-    -------
+            species_name   : return the number of scraped particles of this species
+            boundary       : the boundary from which to get the scraped particle data.
+                            In the form x/y/z_hi/lo
 
-        A List of numpy arrays.
+        Returns
+        -------
 
-    '''
+            The number of particles scraped so far from a boundary and of a species.
 
-    return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getEfieldFP, level, direction, include_ghosts)
+        '''
+        return self.libwarpx_so.warpx_getParticleBoundaryBufferSize(
+            ctypes.c_char_p(species_name.encode('utf-8')),
+            self.get_boundary_number(boundary)
+        )
 
+    def get_particle_boundary_buffer_structs(self, species_name, boundary, level):
+        '''
 
-def get_mesh_electric_field_cp_pml(level, direction, include_ghosts=True):
-    '''
+        This returns a list of numpy arrays containing the particle struct data
+        for a species that has been scraped by a specific simulation boundary. The
+        particle data is represented as a structured numpy array and contains the
+        particle 'x', 'y', 'z', 'id', and 'cpu'.
 
-    This returns a list of numpy arrays containing the mesh electric field
-    data on each grid for this process. This version returns the field on
-    the coarse patch for the PML for the given level.
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+        Parameters
+        ----------
 
-    Parameters
-    ----------
+            species_name : the species name that the data will be returned for
+            boundary     : the boundary from which to get the scraped particle data.
+                        In the form x/y/z_hi/lo or eb.
+            level        : Which AMR level to retrieve scraped particle data from.
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+        Returns
+        -------
 
-    Returns
-    -------
+            A List of numpy arrays.
 
-        A List of numpy arrays.
+        '''
 
-    '''
+        particles_per_tile = _LP_c_int()
+        num_tiles = ctypes.c_int(0)
+        data = self.libwarpx_so.warpx_getParticleBoundaryBufferStructs(
+                ctypes.c_char_p(species_name.encode('utf-8')),
+                self.get_boundary_number(boundary), level,
+                ctypes.byref(num_tiles), ctypes.byref(particles_per_tile)
+        )
 
-    try:
-        return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getEfieldCP_PML, level, direction, include_ghosts)
-    except ValueError:
-        raise Exception('PML not initialized')
+        particle_data = []
+        for i in range(num_tiles.value):
+            if particles_per_tile[i] == 0:
+                continue
+            arr = self._array1d_from_pointer(data[i], self._p_dtype, particles_per_tile[i])
+            particle_data.append(arr)
 
+        _libc.free(particles_per_tile)
+        _libc.free(data)
+        return particle_data
 
-def get_mesh_electric_field_fp_pml(level, direction, include_ghosts=True):
-    '''
+    def get_particle_boundary_buffer(self, species_name, boundary, comp_name, level):
+        '''
 
-    This returns a list of numpy arrays containing the mesh electric field
-    data on each grid for this process. This version returns the field on
-    the fine patch for the PML for the given level.
+        This returns a list of numpy arrays containing the particle array data
+        for a species that has been scraped by a specific simulation boundary.
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-    Parameters
-    ----------
+        Parameters
+        ----------
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+            species_name   : the species name that the data will be returned for.
+            boundary       : the boundary from which to get the scraped particle data.
+                            In the form x/y/z_hi/lo or eb.
+            comp_name      : the component of the array data that will be returned.
+                            If "step_scraped" the special attribute holding the
+                            timestep at which a particle was scraped will be
+                            returned.
+            level          : Which AMR level to retrieve scraped particle data from.
 
-    Returns
-    -------
+        Returns
+        -------
 
-        A List of numpy arrays.
+            A List of numpy arrays.
 
-    '''
+        '''
+        particles_per_tile = _LP_c_int()
+        num_tiles = ctypes.c_int(0)
+        if comp_name == 'step_scraped':
+            data = self.libwarpx_so.warpx_getParticleBoundaryBufferScrapedSteps(
+                ctypes.c_char_p(species_name.encode('utf-8')),
+                self.get_boundary_number(boundary), level,
+                ctypes.byref(num_tiles), ctypes.byref(particles_per_tile)
+            )
+        else:
+            data = self.libwarpx_so.warpx_getParticleBoundaryBuffer(
+                ctypes.c_char_p(species_name.encode('utf-8')),
+                self.get_boundary_number(boundary), level,
+                ctypes.byref(num_tiles), ctypes.byref(particles_per_tile),
+                ctypes.c_char_p(comp_name.encode('utf-8'))
+            )
 
-    try:
-        return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getEfieldFP_PML, level, direction, include_ghosts)
-    except ValueError:
-        raise Exception('PML not initialized')
+        particle_data = []
+        for i in range(num_tiles.value):
+            if particles_per_tile[i] == 0:
+                continue
+            if not data[i]:
+                raise Exception(f'get_particle_arrays: data[i] for i={i} was not initialized')
+            arr = np.ctypeslib.as_array(data[i], (particles_per_tile[i],))
+            try:
+                # This fails on some versions of numpy
+                arr.setflags(write=1)
+            except ValueError:
+                pass
+            particle_data.append(arr)
+
+        _libc.free(particles_per_tile)
+        _libc.free(data)
+        return particle_data
+
+    def clearParticleBoundaryBuffer(self):
+        '''
 
+        Clear the buffer that holds the particles lost at the boundaries.
 
-def get_mesh_magnetic_field(level, direction, include_ghosts=True):
-    '''
+        '''
+        self.libwarpx_so.warpx_clearParticleBoundaryBuffer()
+
+    def _get_mesh_field_list(self, warpx_func, level, direction, include_ghosts):
+        """
+        Generic routine to fetch the list of field data arrays.
+        """
+        shapes = _LP_c_int()
+        size = ctypes.c_int(0)
+        ncomps = ctypes.c_int(0)
+        ngrowvect = _LP_c_int()
+        if direction is None:
+            data = warpx_func(level,
+                            ctypes.byref(size), ctypes.byref(ncomps),
+                            ctypes.byref(ngrowvect), ctypes.byref(shapes))
+        else:
+            data = warpx_func(level, direction,
+                            ctypes.byref(size), ctypes.byref(ncomps),
+                            ctypes.byref(ngrowvect), ctypes.byref(shapes))
+        if not data:
+            raise Exception('object was not initialized')
+
+        ngvect = [ngrowvect[i] for i in range(self.dim)]
+        grid_data = []
+        shapesize = self.dim
+        if ncomps.value > 1:
+            shapesize += 1
+        for i in range(size.value):
+            shape = tuple([shapes[shapesize*i + d] for d in range(shapesize)])
+            # --- The data is stored in Fortran order, hence shape is reversed and a transpose is taken.
+            if shape[::-1] == 0:
+                continue
+            if not data[i]:
+                raise Exception(f'get_particle_arrays: data[i] for i={i} was not initialized')
+            arr = np.ctypeslib.as_array(data[i], shape[::-1]).T
+            try:
+                # This fails on some versions of numpy
+                arr.setflags(write=1)
+            except ValueError:
+                pass
+            if include_ghosts:
+                grid_data.append(arr)
+            else:
+                grid_data.append(arr[tuple([slice(ngvect[d], -ngvect[d]) for d in range(self.dim)])])
 
-    This returns a list of numpy arrays containing the mesh magnetic field
-    data on each grid for this process.
+        _libc.free(shapes)
+        _libc.free(data)
+        return grid_data
 
-    This version is for the full "auxiliary" solution on the given level.
+    def get_mesh_electric_field(self, level, direction, include_ghosts=True):
+        '''
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+        This returns a list of numpy arrays containing the mesh electric field
+        data on each grid for this process.
 
-    Parameters
-    ----------
+        This version is for the full "auxiliary" solution on the given level.
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-    Returns
-    -------
+        Parameters
+        ----------
 
-        A List of numpy arrays.
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-    '''
+        Returns
+        -------
 
-    return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getBfield, level, direction, include_ghosts)
+            A List of numpy arrays.
 
+        '''
 
-def get_mesh_magnetic_field_cp(level, direction, include_ghosts=True):
-    '''
+        return self._get_mesh_field_list(self.libwarpx_so.warpx_getEfield, level, direction, include_ghosts)
 
-    This returns a list of numpy arrays containing the mesh magnetic field
-    data on each grid for this process. This version returns the field on
-    the coarse patch for the given level.
+    def get_mesh_electric_field_cp(self, level, direction, include_ghosts=True):
+        '''
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+        This returns a list of numpy arrays containing the mesh electric field
+        data on each grid for this process. This version returns the field on
+        the coarse patch for the given level.
 
-    Parameters
-    ----------
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+        Parameters
+        ----------
 
-    Returns
-    -------
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-        A List of numpy arrays.
+        Returns
+        -------
 
-    '''
+            A List of numpy arrays.
 
-    return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getBfieldCP, level, direction, include_ghosts)
+        '''
 
+        return self._get_mesh_field_list(self.libwarpx_so.warpx_getEfieldCP, level, direction, include_ghosts)
 
-def get_mesh_magnetic_field_fp(level, direction, include_ghosts=True):
-    '''
+    def get_mesh_electric_field_fp(self, level, direction, include_ghosts=True):
+        '''
 
-    This returns a list of numpy arrays containing the mesh magnetic field
-    data on each grid for this process. This version returns the field on
-    the fine patch for the given level.
+        This returns a list of numpy arrays containing the mesh electric field
+        data on each grid for this process. This version returns the field on
+        the fine patch for the given level.
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-    Parameters
-    ----------
+        Parameters
+        ----------
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-    Returns
-    -------
+        Returns
+        -------
 
-        A List of numpy arrays.
+            A List of numpy arrays.
 
-    '''
+        '''
 
-    return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getBfieldFP, level, direction, include_ghosts)
+        return self._get_mesh_field_list(self.libwarpx_so.warpx_getEfieldFP, level, direction, include_ghosts)
 
+    def get_mesh_electric_field_cp_pml(self, level, direction, include_ghosts=True):
+        '''
 
-def get_mesh_magnetic_field_cp_pml(level, direction, include_ghosts=True):
-    '''
+        This returns a list of numpy arrays containing the mesh electric field
+        data on each grid for this process. This version returns the field on
+        the coarse patch for the PML for the given level.
 
-    This returns a list of numpy arrays containing the mesh magnetic field
-    data on each grid for this process. This version returns the field on
-    the coarse patch for the PML for the given level.
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+        Parameters
+        ----------
 
-    Parameters
-    ----------
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+        Returns
+        -------
 
-    Returns
-    -------
+            A List of numpy arrays.
 
-        A List of numpy arrays.
+        '''
 
-    '''
+        try:
+            return self._get_mesh_field_list(self.libwarpx_so.warpx_getEfieldCP_PML, level, direction, include_ghosts)
+        except ValueError:
+            raise Exception('PML not initialized')
 
-    try:
-        return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getBfieldCP_PML, level, direction, include_ghosts)
-    except ValueError:
-        raise Exception('PML not initialized')
+    def get_mesh_electric_field_fp_pml(self, level, direction, include_ghosts=True):
+        '''
 
+        This returns a list of numpy arrays containing the mesh electric field
+        data on each grid for this process. This version returns the field on
+        the fine patch for the PML for the given level.
 
-def get_mesh_magnetic_field_fp_pml(level, direction, include_ghosts=True):
-    '''
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-    This returns a list of numpy arrays containing the mesh magnetic field
-    data on each grid for this process. This version returns the field on
-    the fine patch for the PML for the given level.
+        Parameters
+        ----------
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-    Parameters
-    ----------
+        Returns
+        -------
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+            A List of numpy arrays.
 
-    Returns
-    -------
+        '''
 
-        A List of numpy arrays.
+        try:
+            return self._get_mesh_field_list(self.libwarpx_so.warpx_getEfieldFP_PML, level, direction, include_ghosts)
+        except ValueError:
+            raise Exception('PML not initialized')
 
-    '''
+    def get_mesh_magnetic_field(self, level, direction, include_ghosts=True):
+        '''
 
-    try:
-        return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getBfieldFP_PML, level, direction, include_ghosts)
-    except ValueError:
-        raise Exception('PML not initialized')
+        This returns a list of numpy arrays containing the mesh magnetic field
+        data on each grid for this process.
 
+        This version is for the full "auxiliary" solution on the given level.
 
-def get_mesh_current_density(level, direction, include_ghosts=True):
-    '''
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-    This returns a list of numpy arrays containing the mesh current density
-    data on each grid for this process.
+        Parameters
+        ----------
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-    Parameters
-    ----------
+        Returns
+        -------
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+            A List of numpy arrays.
 
-    Returns
-    -------
+        '''
 
-        A List of numpy arrays.
+        return self._get_mesh_field_list(self.libwarpx_so.warpx_getBfield, level, direction, include_ghosts)
 
-    '''
+    def get_mesh_magnetic_field_cp(self, level, direction, include_ghosts=True):
+        '''
 
-    return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getCurrentDensity, level, direction, include_ghosts)
+        This returns a list of numpy arrays containing the mesh magnetic field
+        data on each grid for this process. This version returns the field on
+        the coarse patch for the given level.
 
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-def get_mesh_current_density_cp(level, direction, include_ghosts=True):
-    '''
+        Parameters
+        ----------
 
-    This returns a list of numpy arrays containing the mesh current density
-    data on each grid for this process. This version returns the density for
-    the coarse patch on the given level.
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+        Returns
+        -------
 
-    Parameters
-    ----------
+            A List of numpy arrays.
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+        '''
 
-    Returns
-    -------
+        return self._get_mesh_field_list(self.libwarpx_so.warpx_getBfieldCP, level, direction, include_ghosts)
 
-        A List of numpy arrays.
+    def get_mesh_magnetic_field_fp(self, level, direction, include_ghosts=True):
+        '''
 
-    '''
+        This returns a list of numpy arrays containing the mesh magnetic field
+        data on each grid for this process. This version returns the field on
+        the fine patch for the given level.
 
-    return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getCurrentDensityCP, level, direction, include_ghosts)
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
+        Parameters
+        ----------
 
-def get_mesh_current_density_fp(level, direction, include_ghosts=True):
-    '''
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-    This returns a list of numpy arrays containing the mesh current density
-    data on each grid for this process. This version returns the density on
-    the fine patch for the given level.
+        Returns
+        -------
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+            A List of numpy arrays.
 
-    Parameters
-    ----------
+        '''
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+        return self._get_mesh_field_list(self.libwarpx_so.warpx_getBfieldFP, level, direction, include_ghosts)
 
-    Returns
-    -------
+    def get_mesh_magnetic_field_cp_pml(self, level, direction, include_ghosts=True):
+        '''
 
-        A List of numpy arrays.
+        This returns a list of numpy arrays containing the mesh magnetic field
+        data on each grid for this process. This version returns the field on
+        the coarse patch for the PML for the given level.
 
-    '''
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-    return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getCurrentDensityFP, level, direction, include_ghosts)
+        Parameters
+        ----------
 
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-def get_mesh_current_density_cp_pml(level, direction, include_ghosts=True):
-    '''
+        Returns
+        -------
 
-    This returns a list of numpy arrays containing the mesh current density
-    data on each grid for this process. This version returns the density for
-    the coarse patch for the PML for the given level.
+            A List of numpy arrays.
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+        '''
 
-    Parameters
-    ----------
+        try:
+            return self._get_mesh_field_list(self.libwarpx_so.warpx_getBfieldCP_PML, level, direction, include_ghosts)
+        except ValueError:
+            raise Exception('PML not initialized')
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+    def get_mesh_magnetic_field_fp_pml(self, level, direction, include_ghosts=True):
+        '''
 
-    Returns
-    -------
+        This returns a list of numpy arrays containing the mesh magnetic field
+        data on each grid for this process. This version returns the field on
+        the fine patch for the PML for the given level.
 
-        A List of numpy arrays.
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-    '''
+        Parameters
+        ----------
 
-    try:
-        return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getCurrentDensityCP_PML, level, direction, include_ghosts)
-    except ValueError:
-        raise Exception('PML not initialized')
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
+        Returns
+        -------
 
-def get_mesh_current_density_fp_pml(level, direction, include_ghosts=True):
-    '''
+            A List of numpy arrays.
 
-    This returns a list of numpy arrays containing the mesh current density
-    data on each grid for this process. This version returns the density on
-    the fine patch for the PML for the given level.
+        '''
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+        try:
+            return self._get_mesh_field_list(self.libwarpx_so.warpx_getBfieldFP_PML, level, direction, include_ghosts)
+        except ValueError:
+            raise Exception('PML not initialized')
 
-    Parameters
-    ----------
+    def get_mesh_current_density(self, level, direction, include_ghosts=True):
+        '''
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+        This returns a list of numpy arrays containing the mesh current density
+        data on each grid for this process.
 
-    Returns
-    -------
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-        A List of numpy arrays.
+        Parameters
+        ----------
 
-    '''
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-    try:
-        return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getCurrentDensityFP_PML, level, direction, include_ghosts)
-    except ValueError:
-        raise Exception('PML not initialized')
+        Returns
+        -------
 
+            A List of numpy arrays.
 
-def get_mesh_charge_density_cp(level, include_ghosts=True):
-    '''
+        '''
 
-    This returns a list of numpy arrays containing the mesh charge density
-    data on each grid for this process. This version returns the density for
-    the coarse patch on the given level.
+        return self._get_mesh_field_list(self.libwarpx_so.warpx_getCurrentDensity, level, direction, include_ghosts)
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+    def get_mesh_current_density_cp(self, level, direction, include_ghosts=True):
+        '''
 
-    Parameters
-    ----------
+        This returns a list of numpy arrays containing the mesh current density
+        data on each grid for this process. This version returns the density for
+        the coarse patch on the given level.
 
-        level          : the AMR level to get the data for
-        include_ghosts : whether to include ghost zones or not
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-    Returns
-    -------
+        Parameters
+        ----------
 
-        A List of numpy arrays.
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-    '''
+        Returns
+        -------
 
-    return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getChargeDensityCP, level, None, include_ghosts)
+            A List of numpy arrays.
 
+        '''
 
-def get_mesh_charge_density_fp(level, include_ghosts=True):
-    '''
+        return self._get_mesh_field_list(self.libwarpx_so.warpx_getCurrentDensityCP, level, direction, include_ghosts)
 
-    This returns a list of numpy arrays containing the mesh charge density
-    data on each grid for this process. This version returns the density on
-    the fine patch for the given level.
+    def get_mesh_current_density_fp(self, level, direction, include_ghosts=True):
+        '''
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+        This returns a list of numpy arrays containing the mesh current density
+        data on each grid for this process. This version returns the density on
+        the fine patch for the given level.
 
-    Parameters
-    ----------
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-        level          : the AMR level to get the data for
-        include_ghosts : whether to include ghost zones or not
+        Parameters
+        ----------
 
-    Returns
-    -------
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-        A List of numpy arrays.
+        Returns
+        -------
 
-    '''
+            A List of numpy arrays.
 
-    return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getChargeDensityFP, level, None, include_ghosts)
+        '''
 
+        return self._get_mesh_field_list(self.libwarpx_so.warpx_getCurrentDensityFP, level, direction, include_ghosts)
 
-def get_mesh_phi_fp(level, include_ghosts=True):
-    '''
+    def get_mesh_current_density_cp_pml(self, level, direction, include_ghosts=True):
+        '''
 
-    This returns a list of numpy arrays containing the mesh electrostatic
-    potential data on each grid for this process. This version returns the
-    potential on the fine patch for the given level.
+        This returns a list of numpy arrays containing the mesh current density
+        data on each grid for this process. This version returns the density for
+        the coarse patch for the PML for the given level.
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-    Parameters
-    ----------
+        Parameters
+        ----------
 
-        level          : the AMR level to get the data for
-        include_ghosts : whether to include ghost zones or not
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-    Returns
-    -------
+        Returns
+        -------
 
-        A List of numpy arrays.
+            A List of numpy arrays.
 
-    '''
-    return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getPhiFP, level, None, include_ghosts)
+        '''
 
+        try:
+            return self._get_mesh_field_list(self.libwarpx_so.warpx_getCurrentDensityCP_PML, level, direction, include_ghosts)
+        except ValueError:
+            raise Exception('PML not initialized')
 
-def get_mesh_F_cp(level, include_ghosts=True):
-    '''
+    def get_mesh_current_density_fp_pml(self, level, direction, include_ghosts=True):
+        '''
 
-    This returns a list of numpy arrays containing the mesh F field
-    data on each grid for this process. This version returns the F field for
-    the coarse patch on the given level.
+        This returns a list of numpy arrays containing the mesh current density
+        data on each grid for this process. This version returns the density on
+        the fine patch for the PML for the given level.
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-    Parameters
-    ----------
+        Parameters
+        ----------
 
-        level          : the AMR level to get the data for
-        include_ghosts : whether to include ghost zones or not
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-    Returns
-    -------
+        Returns
+        -------
 
-        A List of numpy arrays.
+            A List of numpy arrays.
 
-    '''
+        '''
 
-    return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getFfieldCP, level, None, include_ghosts)
+        try:
+            return self._get_mesh_field_list(self.libwarpx_so.warpx_getCurrentDensityFP_PML, level, direction, include_ghosts)
+        except ValueError:
+            raise Exception('PML not initialized')
 
+    def get_mesh_charge_density_cp(self, level, include_ghosts=True):
+        '''
 
-def get_mesh_F_fp(level, include_ghosts=True):
-    '''
+        This returns a list of numpy arrays containing the mesh charge density
+        data on each grid for this process. This version returns the density for
+        the coarse patch on the given level.
 
-    This returns a list of numpy arrays containing the mesh F field
-    data on each grid for this process. This version returns the F field on
-    the fine patch for the given level.
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+        Parameters
+        ----------
 
-    Parameters
-    ----------
+            level          : the AMR level to get the data for
+            include_ghosts : whether to include ghost zones or not
 
-        level          : the AMR level to get the data for
-        include_ghosts : whether to include ghost zones or not
+        Returns
+        -------
 
-    Returns
-    -------
+            A List of numpy arrays.
 
-        A List of numpy arrays.
+        '''
 
-    '''
+        return self._get_mesh_field_list(self.libwarpx_so.warpx_getChargeDensityCP, level, None, include_ghosts)
 
-    return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getFfieldFP, level, None, include_ghosts)
+    def get_mesh_charge_density_fp(self, level, include_ghosts=True):
+        '''
 
+        This returns a list of numpy arrays containing the mesh charge density
+        data on each grid for this process. This version returns the density on
+        the fine patch for the given level.
 
-def get_mesh_F_fp_pml(level, include_ghosts=True):
-    '''
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-    This returns a list of numpy arrays containing the mesh F field
-    data on each grid for this process. This version returns the F field on
-    the fine patch for the PML on the given level.
+        Parameters
+        ----------
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+            level          : the AMR level to get the data for
+            include_ghosts : whether to include ghost zones or not
 
-    Parameters
-    ----------
+        Returns
+        -------
 
-        level          : the AMR level to get the data for
-        include_ghosts : whether to include ghost zones or not
+            A List of numpy arrays.
 
-    Returns
-    -------
+        '''
 
-        A List of numpy arrays.
+        return self._get_mesh_field_list(self.libwarpx_so.warpx_getChargeDensityFP, level, None, include_ghosts)
 
-    '''
-    try:
-        return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getFfieldFP_PML, level, None, include_ghosts)
-    except ValueError:
-        raise Exception('PML not initialized')
+    def get_mesh_phi_fp(self, level, include_ghosts=True):
+        '''
 
+        This returns a list of numpy arrays containing the mesh electrostatic
+        potential data on each grid for this process. This version returns the
+        potential on the fine patch for the given level.
 
-def get_mesh_F_cp_pml(level, include_ghosts=True):
-    '''
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-    This returns a list of numpy arrays containing the mesh F field
-    data on each grid for this process. This version returns the F field on
-    the coarse patch for the PML on the given level.
+        Parameters
+        ----------
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+            level          : the AMR level to get the data for
+            include_ghosts : whether to include ghost zones or not
 
-    Parameters
-    ----------
+        Returns
+        -------
 
-        level          : the AMR level to get the data for
-        include_ghosts : whether to include ghost zones or not
+            A List of numpy arrays.
 
-    Returns
-    -------
+        '''
+        return self._get_mesh_field_list(self.libwarpx_so.warpx_getPhiFP, level, None, include_ghosts)
 
-        A List of numpy arrays.
+    def get_mesh_F_cp(self, level, include_ghosts=True):
+        '''
 
-    '''
-    try:
-        return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getFfieldCP_PML, level, None, include_ghosts)
-    except ValueError:
-        raise Exception('PML not initialized')
+        This returns a list of numpy arrays containing the mesh F field
+        data on each grid for this process. This version returns the F field for
+        the coarse patch on the given level.
 
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-def get_mesh_G_cp(level, include_ghosts=True):
-    '''
+        Parameters
+        ----------
 
-    This returns a list of numpy arrays containing the mesh G field
-    data on each grid for this process. This version returns the G field for
-    the coarse patch on the given level.
+            level          : the AMR level to get the data for
+            include_ghosts : whether to include ghost zones or not
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+        Returns
+        -------
 
-    Parameters
-    ----------
+            A List of numpy arrays.
 
-        level          : the AMR level to get the data for
-        include_ghosts : whether to include ghost zones or not
+        '''
 
-    Returns
-    -------
+        return self._get_mesh_field_list(self.libwarpx_so.warpx_getFfieldCP, level, None, include_ghosts)
 
-        A List of numpy arrays.
+    def get_mesh_F_fp(self, level, include_ghosts=True):
+        '''
 
-    '''
+        This returns a list of numpy arrays containing the mesh F field
+        data on each grid for this process. This version returns the F field on
+        the fine patch for the given level.
 
-    return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getGfieldCP, level, None, include_ghosts)
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
+        Parameters
+        ----------
 
-def get_mesh_G_fp(level, include_ghosts=True):
-    '''
+            level          : the AMR level to get the data for
+            include_ghosts : whether to include ghost zones or not
 
-    This returns a list of numpy arrays containing the mesh G field
-    data on each grid for this process. This version returns the G field on
-    the fine patch for the given level.
+        Returns
+        -------
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+            A List of numpy arrays.
 
-    Parameters
-    ----------
+        '''
 
-        level          : the AMR level to get the data for
-        include_ghosts : whether to include ghost zones or not
+        return self._get_mesh_field_list(self.libwarpx_so.warpx_getFfieldFP, level, None, include_ghosts)
 
-    Returns
-    -------
+    def get_mesh_F_fp_pml(self, level, include_ghosts=True):
+        '''
 
-        A List of numpy arrays.
+        This returns a list of numpy arrays containing the mesh F field
+        data on each grid for this process. This version returns the F field on
+        the fine patch for the PML on the given level.
 
-    '''
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-    return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getGfieldFP, level, None, include_ghosts)
+        Parameters
+        ----------
 
+            level          : the AMR level to get the data for
+            include_ghosts : whether to include ghost zones or not
 
-def get_mesh_G_cp_pml(level, include_ghosts=True):
-    '''
+        Returns
+        -------
 
-    This returns a list of numpy arrays containing the mesh G field
-    data on each grid for this process. This version returns the G field on
-    the coarse patch for the PML on the given level.
+            A List of numpy arrays.
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+        '''
+        try:
+            return self._get_mesh_field_list(self.libwarpx_so.warpx_getFfieldFP_PML, level, None, include_ghosts)
+        except ValueError:
+            raise Exception('PML not initialized')
 
-    Parameters
-    ----------
+    def get_mesh_F_cp_pml(self, level, include_ghosts=True):
+        '''
 
-        level          : the AMR level to get the data for
-        include_ghosts : whether to include ghost zones or not
+        This returns a list of numpy arrays containing the mesh F field
+        data on each grid for this process. This version returns the F field on
+        the coarse patch for the PML on the given level.
 
-    Returns
-    -------
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-        A List of numpy arrays.
+        Parameters
+        ----------
 
-    '''
-    try:
-        return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getGfieldCP_PML, level, None, include_ghosts)
-    except ValueError:
-        raise Exception('PML not initialized')
+            level          : the AMR level to get the data for
+            include_ghosts : whether to include ghost zones or not
 
+        Returns
+        -------
 
-def get_mesh_G_fp_pml(level, include_ghosts=True):
-    '''
+            A List of numpy arrays.
 
-    This returns a list of numpy arrays containing the mesh G field
-    data on each grid for this process. This version returns the G field on
-    the fine patch for the PML on the given level.
+        '''
+        try:
+            return self._get_mesh_field_list(self.libwarpx_so.warpx_getFfieldCP_PML, level, None, include_ghosts)
+        except ValueError:
+            raise Exception('PML not initialized')
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+    def get_mesh_G_cp(self, level, include_ghosts=True):
+        '''
 
-    Parameters
-    ----------
+        This returns a list of numpy arrays containing the mesh G field
+        data on each grid for this process. This version returns the G field for
+        the coarse patch on the given level.
 
-        level          : the AMR level to get the data for
-        include_ghosts : whether to include ghost zones or not
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-    Returns
-    -------
+        Parameters
+        ----------
 
-        A List of numpy arrays.
+            level          : the AMR level to get the data for
+            include_ghosts : whether to include ghost zones or not
 
-    '''
-    try:
-        return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getGfieldFP_PML, level, None, include_ghosts)
-    except ValueError:
-        raise Exception('PML not initialized')
+        Returns
+        -------
 
+            A List of numpy arrays.
 
-def get_mesh_edge_lengths(level, direction, include_ghosts=True):
-    '''
+        '''
 
-    This returns a list of numpy arrays containing the mesh edge lengths
-    data on each grid for this process. This version returns the density on
-    the fine patch for the given level.
+        return self._get_mesh_field_list(self.libwarpx_so.warpx_getGfieldCP, level, None, include_ghosts)
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+    def get_mesh_G_fp(self, level, include_ghosts=True):
+        '''
 
-    Parameters
-    ----------
+        This returns a list of numpy arrays containing the mesh G field
+        data on each grid for this process. This version returns the G field on
+        the fine patch for the given level.
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-    Returns
-    -------
+        Parameters
+        ----------
 
-        A List of numpy arrays.
+            level          : the AMR level to get the data for
+            include_ghosts : whether to include ghost zones or not
 
-    '''
+        Returns
+        -------
 
-    return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getEdgeLengths, level, direction, include_ghosts)
+            A List of numpy arrays.
 
+        '''
 
-def get_mesh_face_areas(level, direction, include_ghosts=True):
-    '''
+        return self._get_mesh_field_list(self.libwarpx_so.warpx_getGfieldFP, level, None, include_ghosts)
 
-    This returns a list of numpy arrays containing the mesh face areas
-    data on each grid for this process. This version returns the density on
-    the fine patch for the given level.
+    def get_mesh_G_cp_pml(self, level, include_ghosts=True):
+        '''
 
-    The data for the numpy arrays are not copied, but share the underlying
-    memory buffer with WarpX. The numpy arrays are fully writeable.
+        This returns a list of numpy arrays containing the mesh G field
+        data on each grid for this process. This version returns the G field on
+        the coarse patch for the PML on the given level.
 
-    Parameters
-    ----------
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+        Parameters
+        ----------
 
-    Returns
-    -------
+            level          : the AMR level to get the data for
+            include_ghosts : whether to include ghost zones or not
 
-        A List of numpy arrays.
+        Returns
+        -------
 
-    '''
+            A List of numpy arrays.
 
-    return _get_mesh_field_list(libwarpx.libwarpx_so.warpx_getFaceAreas, level, direction, include_ghosts)
+        '''
+        try:
+            return self._get_mesh_field_list(self.libwarpx_so.warpx_getGfieldCP_PML, level, None, include_ghosts)
+        except ValueError:
+            raise Exception('PML not initialized')
 
+    def get_mesh_G_fp_pml(self, level, include_ghosts=True):
+        '''
 
-def _get_mesh_array_lovects(level, direction, include_ghosts=True, getlovectsfunc=None):
-    assert(0 <= level and level <= libwarpx.libwarpx_so.warpx_finestLevel())
+        This returns a list of numpy arrays containing the mesh G field
+        data on each grid for this process. This version returns the G field on
+        the fine patch for the PML on the given level.
 
-    size = ctypes.c_int(0)
-    ngrowvect = _LP_c_int()
-    if direction is None:
-        data = getlovectsfunc(level, ctypes.byref(size), ctypes.byref(ngrowvect))
-    else:
-        data = getlovectsfunc(level, direction, ctypes.byref(size), ctypes.byref(ngrowvect))
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-    if not data:
-        raise Exception('object was not initialized')
+        Parameters
+        ----------
 
-    lovects_ref = np.ctypeslib.as_array(data, (size.value, libwarpx.dim))
+            level          : the AMR level to get the data for
+            include_ghosts : whether to include ghost zones or not
 
-    # --- Make a copy of the data to avoid memory problems
-    # --- Also, take the transpose to give shape (dims, number of grids)
-    lovects = lovects_ref.copy().T
+        Returns
+        -------
 
-    ng = []
-    if include_ghosts:
-        for d in range(libwarpx.dim):
-            ng.append(ngrowvect[d])
-    else:
-        for d in range(libwarpx.dim):
-            ng.append(0)
-            lovects[d,:] += ngrowvect[d]
+            A List of numpy arrays.
 
-    del lovects_ref
-    _libc.free(data)
-    return lovects, ng
+        '''
+        try:
+            return self._get_mesh_field_list(self.libwarpx_so.warpx_getGfieldFP_PML, level, None, include_ghosts)
+        except ValueError:
+            raise Exception('PML not initialized')
 
+    def get_mesh_edge_lengths(self, level, direction, include_ghosts=True):
+        '''
 
-def get_mesh_electric_field_lovects(level, direction, include_ghosts=True):
-    '''
+        This returns a list of numpy arrays containing the mesh edge lengths
+        data on each grid for this process. This version returns the density on
+        the fine patch for the given level.
 
-    This returns a list of the lo vectors of the arrays containing the mesh electric field
-    data on each grid for this process.
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-    This version is for the full "auxiliary" solution on the given level.
+        Parameters
+        ----------
 
-    Parameters
-    ----------
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+        Returns
+        -------
 
-    Returns
-    -------
+            A List of numpy arrays.
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+        '''
 
-    '''
-    return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.libwarpx_so.warpx_getEfieldLoVects)
+        return self._get_mesh_field_list(self.libwarpx_so.warpx_getEdgeLengths, level, direction, include_ghosts)
 
+    def get_mesh_face_areas(self, level, direction, include_ghosts=True):
+        '''
 
-def get_mesh_electric_field_cp_lovects(level, direction, include_ghosts=True):
-    '''
+        This returns a list of numpy arrays containing the mesh face areas
+        data on each grid for this process. This version returns the density on
+        the fine patch for the given level.
 
-    This returns a list of the lo vectors of the arrays containing the mesh electric field
-    data on each grid for this process.
+        The data for the numpy arrays are not copied, but share the underlying
+        memory buffer with WarpX. The numpy arrays are fully writeable.
 
-    Parameters
-    ----------
+        Parameters
+        ----------
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-    Returns
-    -------
+        Returns
+        -------
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+            A List of numpy arrays.
 
-    '''
-    return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.libwarpx_so.warpx_getEfieldCPLoVects)
+        '''
 
+        return self._get_mesh_field_list(self.libwarpx_so.warpx_getFaceAreas, level, direction, include_ghosts)
 
-def get_mesh_electric_field_fp_lovects(level, direction, include_ghosts=True):
-    '''
+    def _get_mesh_array_lovects(self, level, direction, include_ghosts=True, getlovectsfunc=None):
+        assert(0 <= level and level <= self.libwarpx_so.warpx_finestLevel())
 
-    This returns a list of the lo vectors of the arrays containing the mesh electric field
-    data on each grid for this process.
+        size = ctypes.c_int(0)
+        ngrowvect = _LP_c_int()
+        if direction is None:
+            data = getlovectsfunc(level, ctypes.byref(size), ctypes.byref(ngrowvect))
+        else:
+            data = getlovectsfunc(level, direction, ctypes.byref(size), ctypes.byref(ngrowvect))
 
-    Parameters
-    ----------
+        if not data:
+            raise Exception('object was not initialized')
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+        lovects_ref = np.ctypeslib.as_array(data, (size.value, self.dim))
 
-    Returns
-    -------
+        # --- Make a copy of the data to avoid memory problems
+        # --- Also, take the transpose to give shape (dims, number of grids)
+        lovects = lovects_ref.copy().T
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+        ng = []
+        if include_ghosts:
+            for d in range(self.dim):
+                ng.append(ngrowvect[d])
+        else:
+            for d in range(self.dim):
+                ng.append(0)
+                lovects[d,:] += ngrowvect[d]
 
-    '''
-    return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.libwarpx_so.warpx_getEfieldFPLoVects)
+        del lovects_ref
+        _libc.free(data)
+        return lovects, ng
 
+    def get_mesh_electric_field_lovects(self, level, direction, include_ghosts=True):
+        '''
 
-def get_mesh_electric_field_cp_lovects_pml(level, direction, include_ghosts=True):
-    '''
+        This returns a list of the lo vectors of the arrays containing the mesh electric field
+        data on each grid for this process.
 
-    This returns a list of the lo vectors of the arrays containing the mesh electric field
-    data on each PML grid for this process.
+        This version is for the full "auxiliary" solution on the given level.
 
-    Parameters
-    ----------
+        Parameters
+        ----------
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-    Returns
-    -------
+        Returns
+        -------
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-    '''
-    try:
-        return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.libwarpx_so.warpx_getEfieldCPLoVects_PML)
-    except ValueError:
-        raise Exception('PML not initialized')
+        '''
+        return self._get_mesh_array_lovects(level, direction, include_ghosts, self.libwarpx_so.warpx_getEfieldLoVects)
 
+    def get_mesh_electric_field_cp_lovects(self, level, direction, include_ghosts=True):
+        '''
 
-def get_mesh_electric_field_fp_lovects_pml(level, direction, include_ghosts=True):
-    '''
+        This returns a list of the lo vectors of the arrays containing the mesh electric field
+        data on each grid for this process.
 
-    This returns a list of the lo vectors of the arrays containing the mesh electric field
-    data on each PML grid for this process.
+        Parameters
+        ----------
 
-    Parameters
-    ----------
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+        Returns
+        -------
 
-    Returns
-    -------
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+        '''
+        return self._get_mesh_array_lovects(level, direction, include_ghosts, self.libwarpx_so.warpx_getEfieldCPLoVects)
 
-    '''
-    try:
-        return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.libwarpx_so.warpx_getEfieldFPLoVects_PML)
-    except ValueError:
-        raise Exception('PML not initialized')
+    def get_mesh_electric_field_fp_lovects(self, level, direction, include_ghosts=True):
+        '''
 
+        This returns a list of the lo vectors of the arrays containing the mesh electric field
+        data on each grid for this process.
 
-def get_mesh_magnetic_field_lovects(level, direction, include_ghosts=True):
-    '''
+        Parameters
+        ----------
 
-    This returns a list of the lo vectors of the arrays containing the mesh electric field
-    data on each grid for this process.
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-    This version is for the full "auxiliary" solution on the given level.
+        Returns
+        -------
 
-    Parameters
-    ----------
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+        '''
+        return self._get_mesh_array_lovects(level, direction, include_ghosts, self.libwarpx_so.warpx_getEfieldFPLoVects)
 
-    Returns
-    -------
+    def get_mesh_electric_field_cp_lovects_pml(self, level, direction, include_ghosts=True):
+        '''
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+        This returns a list of the lo vectors of the arrays containing the mesh electric field
+        data on each PML grid for this process.
 
-    '''
-    return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.libwarpx_so.warpx_getBfieldLoVects)
+        Parameters
+        ----------
 
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-def get_mesh_magnetic_field_cp_lovects(level, direction, include_ghosts=True):
-    '''
+        Returns
+        -------
 
-    This returns a list of the lo vectors of the arrays containing the mesh electric field
-    data on each grid for this process.
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-    Parameters
-    ----------
+        '''
+        try:
+            return self._get_mesh_array_lovects(level, direction, include_ghosts, self.libwarpx_so.warpx_getEfieldCPLoVects_PML)
+        except ValueError:
+            raise Exception('PML not initialized')
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+    def get_mesh_electric_field_fp_lovects_pml(self, level, direction, include_ghosts=True):
+        '''
 
-    Returns
-    -------
+        This returns a list of the lo vectors of the arrays containing the mesh electric field
+        data on each PML grid for this process.
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+        Parameters
+        ----------
 
-    '''
-    return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.libwarpx_so.warpx_getBfieldCPLoVects)
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
+        Returns
+        -------
 
-def get_mesh_magnetic_field_fp_lovects(level, direction, include_ghosts=True):
-    '''
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-    This returns a list of the lo vectors of the arrays containing the mesh electric field
-    data on each grid for this process.
+        '''
+        try:
+            return self._get_mesh_array_lovects(level, direction, include_ghosts, self.libwarpx_so.warpx_getEfieldFPLoVects_PML)
+        except ValueError:
+            raise Exception('PML not initialized')
 
-    Parameters
-    ----------
+    def get_mesh_magnetic_field_lovects(self, level, direction, include_ghosts=True):
+        '''
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+        This returns a list of the lo vectors of the arrays containing the mesh electric field
+        data on each grid for this process.
 
-    Returns
-    -------
+        This version is for the full "auxiliary" solution on the given level.
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+        Parameters
+        ----------
 
-    '''
-    return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.libwarpx_so.warpx_getBfieldFPLoVects)
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
+        Returns
+        -------
 
-def get_mesh_magnetic_field_cp_lovects_pml(level, direction, include_ghosts=True):
-    '''
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-    This returns a list of the lo vectors of the arrays containing the mesh electric field
-    data on each PML grid for this process.
+        '''
+        return self._get_mesh_array_lovects(level, direction, include_ghosts, self.libwarpx_so.warpx_getBfieldLoVects)
 
-    Parameters
-    ----------
+    def get_mesh_magnetic_field_cp_lovects(self, level, direction, include_ghosts=True):
+        '''
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+        This returns a list of the lo vectors of the arrays containing the mesh electric field
+        data on each grid for this process.
 
-    Returns
-    -------
+        Parameters
+        ----------
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-    '''
-    try:
-        return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.libwarpx_so.warpx_getBfieldCPLoVects_PML)
-    except ValueError:
-        raise Exception('PML not initialized')
+        Returns
+        -------
 
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-def get_mesh_magnetic_field_fp_lovects_pml(level, direction, include_ghosts=True):
-    '''
+        '''
+        return self._get_mesh_array_lovects(level, direction, include_ghosts, self.libwarpx_so.warpx_getBfieldCPLoVects)
 
-    This returns a list of the lo vectors of the arrays containing the mesh electric field
-    data on each PML grid for this process.
+    def get_mesh_magnetic_field_fp_lovects(self, level, direction, include_ghosts=True):
+        '''
 
-    Parameters
-    ----------
+        This returns a list of the lo vectors of the arrays containing the mesh electric field
+        data on each grid for this process.
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+        Parameters
+        ----------
 
-    Returns
-    -------
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+        Returns
+        -------
 
-    '''
-    try:
-        return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.libwarpx_so.warpx_getBfieldFPLoVects_PML)
-    except ValueError:
-        raise Exception('PML not initialized')
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
+        '''
+        return self._get_mesh_array_lovects(level, direction, include_ghosts, self.libwarpx_so.warpx_getBfieldFPLoVects)
 
-def get_mesh_current_density_lovects(level, direction, include_ghosts=True):
-    '''
+    def get_mesh_magnetic_field_cp_lovects_pml(self, level, direction, include_ghosts=True):
+        '''
 
-    This returns a list of the lo vectors of the arrays containing the mesh electric field
-    data on each grid for this process.
+        This returns a list of the lo vectors of the arrays containing the mesh electric field
+        data on each PML grid for this process.
 
-    Parameters
-    ----------
+        Parameters
+        ----------
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-    Returns
-    -------
+        Returns
+        -------
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-    '''
-    return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.libwarpx_so.warpx_getCurrentDensityLoVects)
+        '''
+        try:
+            return self._get_mesh_array_lovects(level, direction, include_ghosts, self.libwarpx_so.warpx_getBfieldCPLoVects_PML)
+        except ValueError:
+            raise Exception('PML not initialized')
 
+    def get_mesh_magnetic_field_fp_lovects_pml(self, level, direction, include_ghosts=True):
+        '''
 
-def get_mesh_current_density_cp_lovects(level, direction, include_ghosts=True):
-    '''
+        This returns a list of the lo vectors of the arrays containing the mesh electric field
+        data on each PML grid for this process.
 
-    This returns a list of the lo vectors of the arrays containing the mesh electric field
-    data on each grid for this process.
+        Parameters
+        ----------
 
-    Parameters
-    ----------
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+        Returns
+        -------
 
-    Returns
-    -------
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+        '''
+        try:
+            return self._get_mesh_array_lovects(level, direction, include_ghosts, self.libwarpx_so.warpx_getBfieldFPLoVects_PML)
+        except ValueError:
+            raise Exception('PML not initialized')
 
-    '''
-    return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.libwarpx_so.warpx_getCurrentDensityCPLoVects)
+    def get_mesh_current_density_lovects(self, level, direction, include_ghosts=True):
+        '''
 
-def get_mesh_current_density_fp_lovects(level, direction, include_ghosts=True):
-    '''
+        This returns a list of the lo vectors of the arrays containing the mesh electric field
+        data on each grid for this process.
 
-    This returns a list of the lo vectors of the arrays containing the mesh electric field
-    data on each grid for this process.
+        Parameters
+        ----------
 
-    Parameters
-    ----------
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+        Returns
+        -------
 
-    Returns
-    -------
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+        '''
+        return self._get_mesh_array_lovects(level, direction, include_ghosts, self.libwarpx_so.warpx_getCurrentDensityLoVects)
 
-    '''
-    return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.libwarpx_so.warpx_getCurrentDensityFPLoVects)
+    def get_mesh_current_density_cp_lovects(self, level, direction, include_ghosts=True):
+        '''
 
+        This returns a list of the lo vectors of the arrays containing the mesh electric field
+        data on each grid for this process.
 
-def get_mesh_current_density_cp_lovects_pml(level, direction, include_ghosts=True):
-    '''
+        Parameters
+        ----------
 
-    This returns a list of the lo vectors of the arrays containing the mesh electric field
-    data on each PML grid for this process.
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-    Parameters
-    ----------
+        Returns
+        -------
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-    Returns
-    -------
+        '''
+        return self._get_mesh_array_lovects(level, direction, include_ghosts, self.libwarpx_so.warpx_getCurrentDensityCPLoVects)
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+    def get_mesh_current_density_fp_lovects(self, level, direction, include_ghosts=True):
+        '''
 
-    '''
-    try:
-        return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.libwarpx_so.warpx_getCurrentDensityCPLoVects_PML)
-    except ValueError:
-        raise Exception('PML not initialized')
+        This returns a list of the lo vectors of the arrays containing the mesh electric field
+        data on each grid for this process.
 
-def get_mesh_current_density_fp_lovects_pml(level, direction, include_ghosts=True):
-    '''
+        Parameters
+        ----------
 
-    This returns a list of the lo vectors of the arrays containing the mesh electric field
-    data on each PML grid for this process.
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-    Parameters
-    ----------
+        Returns
+        -------
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-    Returns
-    -------
+        '''
+        return self._get_mesh_array_lovects(level, direction, include_ghosts, self.libwarpx_so.warpx_getCurrentDensityFPLoVects)
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+    def get_mesh_current_density_cp_lovects_pml(self, level, direction, include_ghosts=True):
+        '''
 
-    '''
-    try:
-        return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.libwarpx_so.warpx_getCurrentDensityFPLoVects_PML)
-    except ValueError:
-        raise Exception('PML not initialized')
+        This returns a list of the lo vectors of the arrays containing the mesh electric field
+        data on each PML grid for this process.
 
+        Parameters
+        ----------
 
-def get_mesh_charge_density_cp_lovects(level, include_ghosts=True):
-    '''
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-    This returns a list of the lo vectors of the arrays containing the mesh electric field
-    data on each grid for this process.
+        Returns
+        -------
 
-    Parameters
-    ----------
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-        level          : the AMR level to get the data for
-        include_ghosts : whether to include ghost zones or not
+        '''
+        try:
+            return self._get_mesh_array_lovects(level, direction, include_ghosts, self.libwarpx_so.warpx_getCurrentDensityCPLoVects_PML)
+        except ValueError:
+            raise Exception('PML not initialized')
 
-    Returns
-    -------
+    def get_mesh_current_density_fp_lovects_pml(self, level, direction, include_ghosts=True):
+        '''
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+        This returns a list of the lo vectors of the arrays containing the mesh electric field
+        data on each PML grid for this process.
 
-    '''
-    return _get_mesh_array_lovects(level, None, include_ghosts, libwarpx.libwarpx_so.warpx_getChargeDensityCPLoVects)
+        Parameters
+        ----------
 
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-def get_mesh_charge_density_fp_lovects(level, include_ghosts=True):
-    '''
+        Returns
+        -------
 
-    This returns a list of the lo vectors of the arrays containing the mesh
-    charge density data on each grid for this process.
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-    Parameters
-    ----------
+        '''
+        try:
+            return self._get_mesh_array_lovects(level, direction, include_ghosts, self.libwarpx_so.warpx_getCurrentDensityFPLoVects_PML)
+        except ValueError:
+            raise Exception('PML not initialized')
 
-        level          : the AMR level to get the data for
-        include_ghosts : whether to include ghost zones or not
+    def get_mesh_charge_density_cp_lovects(self, level, include_ghosts=True):
+        '''
 
-    Returns
-    -------
+        This returns a list of the lo vectors of the arrays containing the mesh electric field
+        data on each grid for this process.
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+        Parameters
+        ----------
 
-    '''
-    return _get_mesh_array_lovects(level, None, include_ghosts, libwarpx.libwarpx_so.warpx_getChargeDensityFPLoVects)
+            level          : the AMR level to get the data for
+            include_ghosts : whether to include ghost zones or not
 
+        Returns
+        -------
 
-def get_mesh_phi_fp_lovects(level, include_ghosts=True):
-    '''
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-    This returns a list of the lo vectors of the arrays containing the mesh
-    electrostatic potential data on each grid for this process.
+        '''
+        return self._get_mesh_array_lovects(level, None, include_ghosts, self.libwarpx_so.warpx_getChargeDensityCPLoVects)
 
-    Parameters
-    ----------
+    def get_mesh_charge_density_fp_lovects(self, level, include_ghosts=True):
+        '''
 
-        level          : the AMR level to get the data for
-        include_ghosts : whether to include ghost zones or not
+        This returns a list of the lo vectors of the arrays containing the mesh
+        charge density data on each grid for this process.
 
-    Returns
-    -------
+        Parameters
+        ----------
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+            level          : the AMR level to get the data for
+            include_ghosts : whether to include ghost zones or not
 
-    '''
-    return _get_mesh_array_lovects(level, None, include_ghosts, libwarpx.libwarpx_so.warpx_getPhiFPLoVects)
+        Returns
+        -------
 
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-def get_mesh_F_cp_lovects(level, include_ghosts=True):
-    '''
+        '''
+        return self._get_mesh_array_lovects(level, None, include_ghosts, self.libwarpx_so.warpx_getChargeDensityFPLoVects)
 
-    This returns a list of the lo vectors of the arrays containing the mesh F field
-    data on each grid for this process.
+    def get_mesh_phi_fp_lovects(self, level, include_ghosts=True):
+        '''
 
-    Parameters
-    ----------
+        This returns a list of the lo vectors of the arrays containing the mesh
+        electrostatic potential data on each grid for this process.
 
-        level          : the AMR level to get the data for
-        include_ghosts : whether to include ghost zones or not
+        Parameters
+        ----------
 
-    Returns
-    -------
+            level          : the AMR level to get the data for
+            include_ghosts : whether to include ghost zones or not
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+        Returns
+        -------
 
-    '''
-    return _get_mesh_array_lovects(level, None, include_ghosts, libwarpx.libwarpx_so.warpx_getFfieldCPLoVects)
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
+        '''
+        return self._get_mesh_array_lovects(level, None, include_ghosts, self.libwarpx_so.warpx_getPhiFPLoVects)
 
-def get_mesh_F_fp_lovects(level, include_ghosts=True):
-    '''
+    def get_mesh_F_cp_lovects(self, level, include_ghosts=True):
+        '''
 
-    This returns a list of the lo vectors of the arrays containing the mesh F field
-    data on each grid for this process.
+        This returns a list of the lo vectors of the arrays containing the mesh F field
+        data on each grid for this process.
 
-    Parameters
-    ----------
+        Parameters
+        ----------
 
-        level          : the AMR level to get the data for
-        include_ghosts : whether to include ghost zones or not
+            level          : the AMR level to get the data for
+            include_ghosts : whether to include ghost zones or not
 
-    Returns
-    -------
+        Returns
+        -------
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-    '''
-    return _get_mesh_array_lovects(level, None, include_ghosts, libwarpx.libwarpx_so.warpx_getFfieldFPLoVects)
+        '''
+        return self._get_mesh_array_lovects(level, None, include_ghosts, self.libwarpx_so.warpx_getFfieldCPLoVects)
 
+    def get_mesh_F_fp_lovects(self, level, include_ghosts=True):
+        '''
 
-def get_mesh_F_cp_lovects_pml(level, include_ghosts=True):
-    '''
+        This returns a list of the lo vectors of the arrays containing the mesh F field
+        data on each grid for this process.
 
-    This returns a list of the lo vectors of the arrays containing the mesh F field
-    data on each PML grid for this process.
+        Parameters
+        ----------
 
-    Parameters
-    ----------
+            level          : the AMR level to get the data for
+            include_ghosts : whether to include ghost zones or not
 
-        level          : the AMR level to get the data for
-        include_ghosts : whether to include ghost zones or not
+        Returns
+        -------
 
-    Returns
-    -------
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+        '''
+        return self._get_mesh_array_lovects(level, None, include_ghosts, self.libwarpx_so.warpx_getFfieldFPLoVects)
 
-    '''
-    try:
-        return _get_mesh_array_lovects(level, None, include_ghosts, libwarpx.libwarpx_so.warpx_getFfieldCPLoVects_PML)
-    except ValueError:
-        raise Exception('PML not initialized')
+    def get_mesh_F_cp_lovects_pml(self, level, include_ghosts=True):
+        '''
 
+        This returns a list of the lo vectors of the arrays containing the mesh F field
+        data on each PML grid for this process.
 
-def get_mesh_F_fp_lovects_pml(level, include_ghosts=True):
-    '''
+        Parameters
+        ----------
 
-    This returns a list of the lo vectors of the arrays containing the mesh F field
-    data on each PML grid for this process.
+            level          : the AMR level to get the data for
+            include_ghosts : whether to include ghost zones or not
 
-    Parameters
-    ----------
+        Returns
+        -------
 
-        level          : the AMR level to get the data for
-        include_ghosts : whether to include ghost zones or not
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-    Returns
-    -------
+        '''
+        try:
+            return self._get_mesh_array_lovects(level, None, include_ghosts, self.libwarpx_so.warpx_getFfieldCPLoVects_PML)
+        except ValueError:
+            raise Exception('PML not initialized')
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+    def get_mesh_F_fp_lovects_pml(self, level, include_ghosts=True):
+        '''
 
-    '''
-    try:
-        return _get_mesh_array_lovects(level, None, include_ghosts, libwarpx.libwarpx_so.warpx_getFfieldFPLoVects_PML)
-    except ValueError:
-        raise Exception('PML not initialized')
+        This returns a list of the lo vectors of the arrays containing the mesh F field
+        data on each PML grid for this process.
 
+        Parameters
+        ----------
 
-def get_mesh_G_cp_lovects(level, include_ghosts=True):
-    '''
+            level          : the AMR level to get the data for
+            include_ghosts : whether to include ghost zones or not
 
-    This returns a list of the lo vectors of the arrays containing the mesh G field
-    data on each grid for this process.
+        Returns
+        -------
 
-    Parameters
-    ----------
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-        level          : the AMR level to get the data for
-        include_ghosts : whether to include ghost zones or not
+        '''
+        try:
+            return self._get_mesh_array_lovects(level, None, include_ghosts, self.libwarpx_so.warpx_getFfieldFPLoVects_PML)
+        except ValueError:
+            raise Exception('PML not initialized')
 
-    Returns
-    -------
+    def get_mesh_G_cp_lovects(self, level, include_ghosts=True):
+        '''
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+        This returns a list of the lo vectors of the arrays containing the mesh G field
+        data on each grid for this process.
 
-    '''
-    return _get_mesh_array_lovects(level, None, include_ghosts, libwarpx.libwarpx_so.warpx_getGfieldCPLoVects)
+        Parameters
+        ----------
 
+            level          : the AMR level to get the data for
+            include_ghosts : whether to include ghost zones or not
 
-def get_mesh_G_fp_lovects(level, include_ghosts=True):
-    '''
+        Returns
+        -------
 
-    This returns a list of the lo vectors of the arrays containing the mesh G field
-    data on each grid for this process.
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-    Parameters
-    ----------
+        '''
+        return self._get_mesh_array_lovects(level, None, include_ghosts, self.libwarpx_so.warpx_getGfieldCPLoVects)
 
-        level          : the AMR level to get the data for
-        include_ghosts : whether to include ghost zones or not
+    def get_mesh_G_fp_lovects(self, level, include_ghosts=True):
+        '''
 
-    Returns
-    -------
+        This returns a list of the lo vectors of the arrays containing the mesh G field
+        data on each grid for this process.
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+        Parameters
+        ----------
 
-    '''
-    return _get_mesh_array_lovects(level, None, include_ghosts, libwarpx.libwarpx_so.warpx_getGfieldFPLoVects)
+            level          : the AMR level to get the data for
+            include_ghosts : whether to include ghost zones or not
 
+        Returns
+        -------
 
-def get_mesh_G_cp_lovects_pml(level, include_ghosts=True):
-    '''
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-    This returns a list of the lo vectors of the arrays containing the mesh G field
-    data on each PML grid for this process.
+        '''
+        return self._get_mesh_array_lovects(level, None, include_ghosts, self.libwarpx_so.warpx_getGfieldFPLoVects)
 
-    Parameters
-    ----------
+    def get_mesh_G_cp_lovects_pml(self, level, include_ghosts=True):
+        '''
 
-        level          : the AMR level to get the data for
-        include_ghosts : whether to include ghost zones or not
+        This returns a list of the lo vectors of the arrays containing the mesh G field
+        data on each PML grid for this process.
 
-    Returns
-    -------
+        Parameters
+        ----------
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+            level          : the AMR level to get the data for
+            include_ghosts : whether to include ghost zones or not
 
-    '''
-    try:
-        return _get_mesh_array_lovects(level, None, include_ghosts, libwarpx.libwarpx_so.warpx_getGfieldCPLoVects_PML)
-    except ValueError:
-        raise Exception('PML not initialized')
+        Returns
+        -------
 
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-def get_mesh_G_fp_lovects_pml(level, include_ghosts=True):
-    '''
+        '''
+        try:
+            return self._get_mesh_array_lovects(level, None, include_ghosts, self.libwarpx_so.warpx_getGfieldCPLoVects_PML)
+        except ValueError:
+            raise Exception('PML not initialized')
 
-    This returns a list of the lo vectors of the arrays containing the mesh G field
-    data on each PML grid for this process.
+    def get_mesh_G_fp_lovects_pml(self, level, include_ghosts=True):
+        '''
 
-    Parameters
-    ----------
+        This returns a list of the lo vectors of the arrays containing the mesh G field
+        data on each PML grid for this process.
 
-        level          : the AMR level to get the data for
-        include_ghosts : whether to include ghost zones or not
+        Parameters
+        ----------
 
-    Returns
-    -------
+            level          : the AMR level to get the data for
+            include_ghosts : whether to include ghost zones or not
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+        Returns
+        -------
 
-    '''
-    try:
-        return _get_mesh_array_lovects(level, None, include_ghosts, libwarpx.libwarpx_so.warpx_getGfieldFPLoVects_PML)
-    except ValueError:
-        raise Exception('PML not initialized')
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
+        '''
+        try:
+            return self._get_mesh_array_lovects(level, None, include_ghosts, self.libwarpx_so.warpx_getGfieldFPLoVects_PML)
+        except ValueError:
+            raise Exception('PML not initialized')
 
-def get_mesh_edge_lengths_lovects(level, direction, include_ghosts=True):
-    '''
+    def get_mesh_edge_lengths_lovects(self, level, direction, include_ghosts=True):
+        '''
 
-    This returns a list of the lo vectors of the arrays containing the mesh edge lengths
-    data on each grid for this process.
+        This returns a list of the lo vectors of the arrays containing the mesh edge lengths
+        data on each grid for this process.
 
-    Parameters
-    ----------
+        Parameters
+        ----------
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-    Returns
-    -------
+        Returns
+        -------
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-    '''
-    return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.libwarpx_so.warpx_getEdgeLengthsLoVects)
+        '''
+        return self._get_mesh_array_lovects(level, direction, include_ghosts, self.libwarpx_so.warpx_getEdgeLengthsLoVects)
 
+    def get_mesh_face_areas_lovects(self, level, direction, include_ghosts=True):
+        '''
 
-def get_mesh_face_areas_lovects(level, direction, include_ghosts=True):
-    '''
+        This returns a list of the lo vectors of the arrays containing the mesh face areas
+        data on each grid for this process.
 
-    This returns a list of the lo vectors of the arrays containing the mesh face areas
-    data on each grid for this process.
+        Parameters
+        ----------
 
-    Parameters
-    ----------
+            level          : the AMR level to get the data for
+            direction      : the component of the data you want
+            include_ghosts : whether to include ghost zones or not
 
-        level          : the AMR level to get the data for
-        direction      : the component of the data you want
-        include_ghosts : whether to include ghost zones or not
+        Returns
+        -------
 
-    Returns
-    -------
+            A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
 
-        A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids)
+        '''
+        return self._get_mesh_array_lovects(level, direction, include_ghosts, self.libwarpx_so.warpx_getFaceAreasLoVects)
 
-    '''
-    return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.libwarpx_so.warpx_getFaceAreasLoVects)
+    def _get_nodal_flag(self, getdatafunc):
+        data = getdatafunc()
+        if not data:
+            raise Exception('object was not initialized')
 
+        nodal_flag_ref = np.ctypeslib.as_array(data, (self.dim,))
 
-def _get_nodal_flag(getdatafunc):
-    data = getdatafunc()
-    if not data:
-        raise Exception('object was not initialized')
+        # --- Make a copy of the data to avoid memory problems
+        nodal_flag = nodal_flag_ref.copy()
 
-    nodal_flag_ref = np.ctypeslib.as_array(data, (libwarpx.dim,))
+        del nodal_flag_ref
+        _libc.free(data)
+        return nodal_flag
 
-    # --- Make a copy of the data to avoid memory problems
-    nodal_flag = nodal_flag_ref.copy()
+    def get_Ex_nodal_flag(self):
+        '''
+        This returns a 1d array of the nodal flags for Ex along each direction. A 1 means node centered, and 0 cell centered.
+        '''
+        return self._get_nodal_flag(self.libwarpx_so.warpx_getEx_nodal_flag)
 
-    del nodal_flag_ref
-    _libc.free(data)
-    return nodal_flag
+    def get_Ey_nodal_flag(self):
+        '''
+        This returns a 1d array of the nodal flags for Ey along each direction. A 1 means node centered, and 0 cell centered.
+        '''
+        return self._get_nodal_flag(self.libwarpx_so.warpx_getEy_nodal_flag)
 
+    def get_Ez_nodal_flag(self):
+        '''
+        This returns a 1d array of the nodal flags for Ez along each direction. A 1 means node centered, and 0 cell centered.
+        '''
+        return self._get_nodal_flag(self.libwarpx_so.warpx_getEz_nodal_flag)
 
-def get_Ex_nodal_flag():
-    '''
-    This returns a 1d array of the nodal flags for Ex along each direction. A 1 means node centered, and 0 cell centered.
-    '''
-    return _get_nodal_flag(libwarpx.libwarpx_so.warpx_getEx_nodal_flag)
+    def get_Bx_nodal_flag(self):
+        '''
+        This returns a 1d array of the nodal flags for Bx along each direction. A 1 means node centered, and 0 cell centered.
+        '''
+        return self._get_nodal_flag(self.libwarpx_so.warpx_getBx_nodal_flag)
 
+    def get_By_nodal_flag(self):
+        '''
+        This returns a 1d array of the nodal flags for By along each direction. A 1 means node centered, and 0 cell centered.
+        '''
+        return self._get_nodal_flag(self.libwarpx_so.warpx_getBy_nodal_flag)
 
-def get_Ey_nodal_flag():
-    '''
-    This returns a 1d array of the nodal flags for Ey along each direction. A 1 means node centered, and 0 cell centered.
-    '''
-    return _get_nodal_flag(libwarpx.libwarpx_so.warpx_getEy_nodal_flag)
+    def get_Bz_nodal_flag(self):
+        '''
+        This returns a 1d array of the nodal flags for Bz along each direction. A 1 means node centered, and 0 cell centered.
+        '''
+        return self._get_nodal_flag(self.libwarpx_so.warpx_getBz_nodal_flag)
 
+    def get_Jx_nodal_flag(self):
+        '''
+        This returns a 1d array of the nodal flags for Jx along each direction. A 1 means node centered, and 0 cell centered.
+        '''
+        return self._get_nodal_flag(self.libwarpx_so.warpx_getJx_nodal_flag)
 
-def get_Ez_nodal_flag():
-    '''
-    This returns a 1d array of the nodal flags for Ez along each direction. A 1 means node centered, and 0 cell centered.
-    '''
-    return _get_nodal_flag(libwarpx.libwarpx_so.warpx_getEz_nodal_flag)
+    def get_Jy_nodal_flag(self):
+        '''
+        This returns a 1d array of the nodal flags for Jy along each direction. A 1 means node centered, and 0 cell centered.
+        '''
+        return self._get_nodal_flag(self.libwarpx_so.warpx_getJy_nodal_flag)
 
+    def get_Jz_nodal_flag(self):
+        '''
+        This returns a 1d array of the nodal flags for Jz along each direction. A 1 means node centered, and 0 cell centered.
+        '''
+        return self._get_nodal_flag(self.libwarpx_so.warpx_getJz_nodal_flag)
 
-def get_Bx_nodal_flag():
-    '''
-    This returns a 1d array of the nodal flags for Bx along each direction. A 1 means node centered, and 0 cell centered.
-    '''
-    return _get_nodal_flag(libwarpx.libwarpx_so.warpx_getBx_nodal_flag)
+    def get_Rho_nodal_flag(self):
+        '''
+        This returns a 1d array of the nodal flags for Rho along each direction. A 1 means node centered, and 0 cell centered.
+        '''
+        return self._get_nodal_flag(self.libwarpx_so.warpx_getRho_nodal_flag)
 
+    def get_Phi_nodal_flag(self):
+        '''
+        This returns a 1d array of the nodal flags for Phi along each direction. A 1 means node centered, and 0 cell centered.
+        '''
+        return self._get_nodal_flag(self.libwarpx_so.warpx_getPhi_nodal_flag)
 
-def get_By_nodal_flag():
-    '''
-    This returns a 1d array of the nodal flags for By along each direction. A 1 means node centered, and 0 cell centered.
-    '''
-    return _get_nodal_flag(libwarpx.libwarpx_so.warpx_getBy_nodal_flag)
+    def get_F_nodal_flag(self):
+        '''
+        This returns a 1d array of the nodal flags for F along each direction. A 1 means node centered, and 0 cell centered.
+        '''
+        return self._get_nodal_flag(self.libwarpx_so.warpx_getF_nodal_flag)
 
+    def get_G_nodal_flag(self):
+        '''
+        This returns a 1d array of the nodal flags for G along each direction. A 1 means node centered, and 0 cell centered.
+        '''
+        return self._get_nodal_flag(self.libwarpx_so.warpx_getG_nodal_flag)
 
-def get_Bz_nodal_flag():
-    '''
-    This returns a 1d array of the nodal flags for Bz along each direction. A 1 means node centered, and 0 cell centered.
-    '''
-    return _get_nodal_flag(libwarpx.libwarpx_so.warpx_getBz_nodal_flag)
+    def get_edge_lengths_x_nodal_flag(self):
+        '''
+        This returns a 1d array of the nodal flags for the x edge lengths along each direction. A 1 means node centered, and 0 cell centered.
+        '''
+        return self._get_nodal_flag(self.libwarpx_so.warpx_get_edge_lengths_x_nodal_flag)
 
+    def get_edge_lengths_y_nodal_flag(self):
+        '''
+        This returns a 1d array of the nodal flags for the y edge lengths along each direction. A 1 means node centered, and 0 cell centered.
+        '''
+        return self._get_nodal_flag(self.libwarpx_so.warpx_get_edge_lengths_y_nodal_flag)
 
-def get_Jx_nodal_flag():
-    '''
-    This returns a 1d array of the nodal flags for Jx along each direction. A 1 means node centered, and 0 cell centered.
-    '''
-    return _get_nodal_flag(libwarpx.libwarpx_so.warpx_getJx_nodal_flag)
+    def get_edge_lengths_z_nodal_flag(self):
+        '''
+        This returns a 1d array of the nodal flags for the z edge lengths along each direction. A 1 means node centered, and 0 cell centered.
+        '''
+        return self._get_nodal_flag(self.libwarpx_so.warpx_get_edge_lengths_z_nodal_flag)
 
+    def get_face_areas_x_nodal_flag(self):
+        '''
+        This returns a 1d array of the nodal flags for the x face areas along each direction. A 1 means node centered, and 0 cell centered.
+        '''
+        return self._get_nodal_flag(self.libwarpx_so.warpx_get_face_areas_x_nodal_flag)
 
-def get_Jy_nodal_flag():
-    '''
-    This returns a 1d array of the nodal flags for Jy along each direction. A 1 means node centered, and 0 cell centered.
-    '''
-    return _get_nodal_flag(libwarpx.libwarpx_so.warpx_getJy_nodal_flag)
+    def get_face_areas_y_nodal_flag(self):
+        '''
+        This returns a 1d array of the nodal flags for the y face areas along each direction. A 1 means node centered, and 0 cell centered.
+        '''
+        return self._get_nodal_flag(self.libwarpx_so.warpx_get_face_areas_y_nodal_flag)
 
+    def get_face_areas_z_nodal_flag(self):
+        '''
+        This returns a 1d array of the nodal flags for the z face areas along each direction. A 1 means node centered, and 0 cell centered.
+        '''
+        return self._get_nodal_flag(self.libwarpx_so.warpx_get_face_areas_z_nodal_flag)
 
-def get_Jz_nodal_flag():
-    '''
-    This returns a 1d array of the nodal flags for Jz along each direction. A 1 means node centered, and 0 cell centered.
-    '''
-    return _get_nodal_flag(libwarpx.libwarpx_so.warpx_getJz_nodal_flag)
+    def get_F_pml_nodal_flag(self):
+        '''
+        This returns a 1d array of the nodal flags for F in the PML along each direction. A 1 means node centered, and 0 cell centered.
+        '''
+        return self._get_nodal_flag(self.libwarpx_so.warpx_getF_pml_nodal_flag)
 
+    def get_G_pml_nodal_flag(self):
+        '''
+        This returns a 1d array of the nodal flags for G in the PML along each direction. A 1 means node centered, and 0 cell centered.
+        '''
+        return self._get_nodal_flag(self.libwarpx_so.warpx_getG_pml_nodal_flag)
 
-def get_Rho_nodal_flag():
-    '''
-    This returns a 1d array of the nodal flags for Rho along each direction. A 1 means node centered, and 0 cell centered.
-    '''
-    return _get_nodal_flag(libwarpx.libwarpx_so.warpx_getRho_nodal_flag)
 
-def get_Phi_nodal_flag():
-    '''
-    This returns a 1d array of the nodal flags for Phi along each direction. A 1 means node centered, and 0 cell centered.
-    '''
-    return _get_nodal_flag(libwarpx.libwarpx_so.warpx_getPhi_nodal_flag)
-
-def get_F_nodal_flag():
-    '''
-    This returns a 1d array of the nodal flags for F along each direction. A 1 means node centered, and 0 cell centered.
-    '''
-    return _get_nodal_flag(libwarpx.libwarpx_so.warpx_getF_nodal_flag)
-
-def get_G_nodal_flag():
-    '''
-    This returns a 1d array of the nodal flags for G along each direction. A 1 means node centered, and 0 cell centered.
-    '''
-    return _get_nodal_flag(libwarpx.libwarpx_so.warpx_getG_nodal_flag)
-
-def get_edge_lengths_x_nodal_flag():
-    '''
-    This returns a 1d array of the nodal flags for the x edge lengths along each direction. A 1 means node centered, and 0 cell centered.
-    '''
-    return _get_nodal_flag(libwarpx.libwarpx_so.warpx_get_edge_lengths_x_nodal_flag)
-
-def get_edge_lengths_y_nodal_flag():
-    '''
-    This returns a 1d array of the nodal flags for the y edge lengths along each direction. A 1 means node centered, and 0 cell centered.
-    '''
-    return _get_nodal_flag(libwarpx.libwarpx_so.warpx_get_edge_lengths_y_nodal_flag)
-
-def get_edge_lengths_z_nodal_flag():
-    '''
-    This returns a 1d array of the nodal flags for the z edge lengths along each direction. A 1 means node centered, and 0 cell centered.
-    '''
-    return _get_nodal_flag(libwarpx.libwarpx_so.warpx_get_edge_lengths_z_nodal_flag)
-
-def get_face_areas_x_nodal_flag():
-    '''
-    This returns a 1d array of the nodal flags for the x face areas along each direction. A 1 means node centered, and 0 cell centered.
-    '''
-    return _get_nodal_flag(libwarpx.libwarpx_so.warpx_get_face_areas_x_nodal_flag)
-
-def get_face_areas_y_nodal_flag():
-    '''
-    This returns a 1d array of the nodal flags for the y face areas along each direction. A 1 means node centered, and 0 cell centered.
-    '''
-    return _get_nodal_flag(libwarpx.libwarpx_so.warpx_get_face_areas_y_nodal_flag)
-
-def get_face_areas_z_nodal_flag():
-    '''
-    This returns a 1d array of the nodal flags for the z face areas along each direction. A 1 means node centered, and 0 cell centered.
-    '''
-    return _get_nodal_flag(libwarpx.libwarpx_so.warpx_get_face_areas_z_nodal_flag)
-
-def get_F_pml_nodal_flag():
-    '''
-    This returns a 1d array of the nodal flags for F in the PML along each direction. A 1 means node centered, and 0 cell centered.
-    '''
-    return _get_nodal_flag(libwarpx.libwarpx_so.warpx_getF_pml_nodal_flag)
-
-def get_G_pml_nodal_flag():
-    '''
-    This returns a 1d array of the nodal flags for G in the PML along each direction. A 1 means node centered, and 0 cell centered.
-    '''
-    return _get_nodal_flag(libwarpx.libwarpx_so.warpx_getG_pml_nodal_flag)
+libwarpx = LibWarpX()
diff --git a/Python/pywarpx/fields.py b/Python/pywarpx/fields.py
index 85b0f2f23..56ab6fd2e 100644
--- a/Python/pywarpx/fields.py
+++ b/Python/pywarpx/fields.py
@@ -22,7 +22,7 @@ try:
 except ImportError:
     npes = 1
 
-from . import _libwarpx
+from ._libwarpx import libwarpx
 
 
 class _MultiFABWrapper(object):
@@ -43,7 +43,7 @@ class _MultiFABWrapper(object):
         self.level = level
         self.include_ghosts = include_ghosts
 
-        self.dim = _libwarpx.libwarpx.dim
+        self.dim = libwarpx.dim
 
         # overlaps is one along the axes where the grid boundaries overlap the neighboring grid,
         # which is the case with node centering.
@@ -86,16 +86,16 @@ class _MultiFABWrapper(object):
         """
 
         try:
-            if _libwarpx.libwarpx.geometry_dim == '3d':
+            if libwarpx.geometry_dim == '3d':
                 idir = ['x', 'y', 'z'].index(direction)
                 celldir = idir
-            elif _libwarpx.libwarpx.geometry_dim == '2d':
+            elif libwarpx.geometry_dim == '2d':
                 idir = ['x', 'z'].index(direction)
                 celldir = 2*idir
-            elif _libwarpx.libwarpx.geometry_dim == 'rz':
+            elif libwarpx.geometry_dim == 'rz':
                 idir = ['r', 'z'].index(direction)
                 celldir = 2*idir
-            elif _libwarpx.geometry_dim == '1d':
+            elif libwarpx.geometry_dim == '1d':
                 idir = ['z'].index(direction)
                 celldir = idir
         except ValueError:
@@ -108,7 +108,7 @@ class _MultiFABWrapper(object):
             nn = comm_world.allreduce(nn, op=mpi.MAX)
 
         # --- Cell size in the direction
-        dd = _libwarpx.getCellSize(celldir, self.level)
+        dd = libwarpx.getCellSize(celldir, self.level)
 
         # --- Get the nodal flag along direction
         nodal_flag = self.get_nodal_flag()[idir]
@@ -300,7 +300,7 @@ class _MultiFABWrapper(object):
         if ncomps > 1 and ic is None:
             sss = tuple(list(sss) + [ncomps])
         # --- Create the array to be returned.
-        resultglobal = np.zeros(sss, dtype=_libwarpx.libwarpx._numpy_real_dtype)
+        resultglobal = np.zeros(sss, dtype=libwarpx._numpy_real_dtype)
 
         datalist = []
         for i in range(len(fields)):
@@ -435,457 +435,457 @@ class _MultiFABWrapper(object):
 
 def ExWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=0,
-                            get_lovects=_libwarpx.get_mesh_electric_field_lovects,
-                            get_fabs=_libwarpx.get_mesh_electric_field,
-                            get_nodal_flag=_libwarpx.get_Ex_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_electric_field_lovects,
+                            get_fabs=libwarpx.get_mesh_electric_field,
+                            get_nodal_flag=libwarpx.get_Ex_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def EyWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=1,
-                            get_lovects=_libwarpx.get_mesh_electric_field_lovects,
-                            get_fabs=_libwarpx.get_mesh_electric_field,
-                            get_nodal_flag=_libwarpx.get_Ey_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_electric_field_lovects,
+                            get_fabs=libwarpx.get_mesh_electric_field,
+                            get_nodal_flag=libwarpx.get_Ey_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def EzWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=2,
-                            get_lovects=_libwarpx.get_mesh_electric_field_lovects,
-                            get_fabs=_libwarpx.get_mesh_electric_field,
-                            get_nodal_flag=_libwarpx.get_Ez_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_electric_field_lovects,
+                            get_fabs=libwarpx.get_mesh_electric_field,
+                            get_nodal_flag=libwarpx.get_Ez_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def BxWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=0,
-                            get_lovects=_libwarpx.get_mesh_magnetic_field_lovects,
-                            get_fabs=_libwarpx.get_mesh_magnetic_field,
-                            get_nodal_flag=_libwarpx.get_Bx_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_magnetic_field_lovects,
+                            get_fabs=libwarpx.get_mesh_magnetic_field,
+                            get_nodal_flag=libwarpx.get_Bx_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def ByWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=1,
-                            get_lovects=_libwarpx.get_mesh_magnetic_field_lovects,
-                            get_fabs=_libwarpx.get_mesh_magnetic_field,
-                            get_nodal_flag=_libwarpx.get_By_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_magnetic_field_lovects,
+                            get_fabs=libwarpx.get_mesh_magnetic_field,
+                            get_nodal_flag=libwarpx.get_By_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def BzWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=2,
-                            get_lovects=_libwarpx.get_mesh_magnetic_field_lovects,
-                            get_fabs=_libwarpx.get_mesh_magnetic_field,
-                            get_nodal_flag=_libwarpx.get_Bz_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_magnetic_field_lovects,
+                            get_fabs=libwarpx.get_mesh_magnetic_field,
+                            get_nodal_flag=libwarpx.get_Bz_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def JxWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=0,
-                            get_lovects=_libwarpx.get_mesh_current_density_lovects,
-                            get_fabs=_libwarpx.get_mesh_current_density,
-                            get_nodal_flag=_libwarpx.get_Jx_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_current_density_lovects,
+                            get_fabs=libwarpx.get_mesh_current_density,
+                            get_nodal_flag=libwarpx.get_Jx_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def JyWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=1,
-                            get_lovects=_libwarpx.get_mesh_current_density_lovects,
-                            get_fabs=_libwarpx.get_mesh_current_density,
-                            get_nodal_flag=_libwarpx.get_Jy_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_current_density_lovects,
+                            get_fabs=libwarpx.get_mesh_current_density,
+                            get_nodal_flag=libwarpx.get_Jy_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def JzWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=2,
-                            get_lovects=_libwarpx.get_mesh_current_density_lovects,
-                            get_fabs=_libwarpx.get_mesh_current_density,
-                            get_nodal_flag=_libwarpx.get_Jz_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_current_density_lovects,
+                            get_fabs=libwarpx.get_mesh_current_density,
+                            get_nodal_flag=libwarpx.get_Jz_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def ExCPWrapper(level=1, include_ghosts=False):
     assert level>0, Exception('Coarse patch only available on levels > 0')
     return _MultiFABWrapper(direction=0,
-                            get_lovects=_libwarpx.get_mesh_electric_field_cp_lovects,
-                            get_fabs=_libwarpx.get_mesh_electric_field_cp,
-                            get_nodal_flag=_libwarpx.get_Ex_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_electric_field_cp_lovects,
+                            get_fabs=libwarpx.get_mesh_electric_field_cp,
+                            get_nodal_flag=libwarpx.get_Ex_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def EyCPWrapper(level=1, include_ghosts=False):
     assert level>0, Exception('Coarse patch only available on levels > 0')
     return _MultiFABWrapper(direction=1,
-                            get_lovects=_libwarpx.get_mesh_electric_field_cp_lovects,
-                            get_fabs=_libwarpx.get_mesh_electric_field_cp,
-                            get_nodal_flag=_libwarpx.get_Ey_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_electric_field_cp_lovects,
+                            get_fabs=libwarpx.get_mesh_electric_field_cp,
+                            get_nodal_flag=libwarpx.get_Ey_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def EzCPWrapper(level=1, include_ghosts=False):
     assert level>0, Exception('Coarse patch only available on levels > 0')
     return _MultiFABWrapper(direction=2,
-                            get_lovects=_libwarpx.get_mesh_electric_field_cp_lovects,
-                            get_fabs=_libwarpx.get_mesh_electric_field_cp,
-                            get_nodal_flag=_libwarpx.get_Ez_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_electric_field_cp_lovects,
+                            get_fabs=libwarpx.get_mesh_electric_field_cp,
+                            get_nodal_flag=libwarpx.get_Ez_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def BxCPWrapper(level=1, include_ghosts=False):
     assert level>0, Exception('Coarse patch only available on levels > 0')
     return _MultiFABWrapper(direction=0,
-                            get_lovects=_libwarpx.get_mesh_magnetic_field_cp_lovects,
-                            get_fabs=_libwarpx.get_mesh_magnetic_field_cp,
-                            get_nodal_flag=_libwarpx.get_Bx_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_magnetic_field_cp_lovects,
+                            get_fabs=libwarpx.get_mesh_magnetic_field_cp,
+                            get_nodal_flag=libwarpx.get_Bx_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def ByCPWrapper(level=1, include_ghosts=False):
     assert level>0, Exception('Coarse patch only available on levels > 0')
     return _MultiFABWrapper(direction=1,
-                            get_lovects=_libwarpx.get_mesh_magnetic_field_cp_lovects,
-                            get_fabs=_libwarpx.get_mesh_magnetic_field_cp,
-                            get_nodal_flag=_libwarpx.get_By_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_magnetic_field_cp_lovects,
+                            get_fabs=libwarpx.get_mesh_magnetic_field_cp,
+                            get_nodal_flag=libwarpx.get_By_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def BzCPWrapper(level=1, include_ghosts=False):
     assert level>0, Exception('Coarse patch only available on levels > 0')
     return _MultiFABWrapper(direction=2,
-                            get_lovects=_libwarpx.get_mesh_magnetic_field_cp_lovects,
-                            get_fabs=_libwarpx.get_mesh_magnetic_field_cp,
-                            get_nodal_flag=_libwarpx.get_Bz_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_magnetic_field_cp_lovects,
+                            get_fabs=libwarpx.get_mesh_magnetic_field_cp,
+                            get_nodal_flag=libwarpx.get_Bz_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def JxCPWrapper(level=1, include_ghosts=False):
     assert level>0, Exception('Coarse patch only available on levels > 0')
     return _MultiFABWrapper(direction=0,
-                            get_lovects=_libwarpx.get_mesh_current_density_cp_lovects,
-                            get_fabs=_libwarpx.get_mesh_current_density_cp,
-                            get_nodal_flag=_libwarpx.get_Jx_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_current_density_cp_lovects,
+                            get_fabs=libwarpx.get_mesh_current_density_cp,
+                            get_nodal_flag=libwarpx.get_Jx_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def JyCPWrapper(level=1, include_ghosts=False):
     assert level>0, Exception('Coarse patch only available on levels > 0')
     return _MultiFABWrapper(direction=1,
-                            get_lovects=_libwarpx.get_mesh_current_density_cp_lovects,
-                            get_fabs=_libwarpx.get_mesh_current_density_cp,
-                            get_nodal_flag=_libwarpx.get_Jy_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_current_density_cp_lovects,
+                            get_fabs=libwarpx.get_mesh_current_density_cp,
+                            get_nodal_flag=libwarpx.get_Jy_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def JzCPWrapper(level=1, include_ghosts=False):
     assert level>0, Exception('Coarse patch only available on levels > 0')
     return _MultiFABWrapper(direction=2,
-                            get_lovects=_libwarpx.get_mesh_current_density_cp_lovects,
-                            get_fabs=_libwarpx.get_mesh_current_density_cp,
-                            get_nodal_flag=_libwarpx.get_Jz_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_current_density_cp_lovects,
+                            get_fabs=libwarpx.get_mesh_current_density_cp,
+                            get_nodal_flag=libwarpx.get_Jz_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def RhoCPWrapper(level=1, include_ghosts=False):
     assert level>0, Exception('Coarse patch only available on levels > 0')
     return _MultiFABWrapper(direction=None,
-                            get_lovects=_libwarpx.get_mesh_charge_density_cp_lovects,
-                            get_fabs=_libwarpx.get_mesh_charge_density_cp,
-                            get_nodal_flag=_libwarpx.get_Rho_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_charge_density_cp_lovects,
+                            get_fabs=libwarpx.get_mesh_charge_density_cp,
+                            get_nodal_flag=libwarpx.get_Rho_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def FCPWrapper(level=1, include_ghosts=False):
     assert level>0, Exception('Coarse patch only available on levels > 0')
     return _MultiFABWrapper(direction=None,
-                            get_lovects=_libwarpx.get_mesh_F_cp_lovects,
-                            get_fabs=_libwarpx.get_mesh_F_cp,
-                            get_nodal_flag=_libwarpx.get_F_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_F_cp_lovects,
+                            get_fabs=libwarpx.get_mesh_F_cp,
+                            get_nodal_flag=libwarpx.get_F_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def GCPWrapper(level=1, include_ghosts=False):
     assert level>0, Exception('Coarse patch only available on levels > 0')
     return _MultiFABWrapper(direction=None,
-                            get_lovects=_libwarpx.get_mesh_G_cp_lovects,
-                            get_fabs=_libwarpx.get_mesh_G_cp,
-                            get_nodal_flag=_libwarpx.get_G_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_G_cp_lovects,
+                            get_fabs=libwarpx.get_mesh_G_cp,
+                            get_nodal_flag=libwarpx.get_G_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def ExFPWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=0,
-                            get_lovects=_libwarpx.get_mesh_electric_field_fp_lovects,
-                            get_fabs=_libwarpx.get_mesh_electric_field_fp,
-                            get_nodal_flag=_libwarpx.get_Ex_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_electric_field_fp_lovects,
+                            get_fabs=libwarpx.get_mesh_electric_field_fp,
+                            get_nodal_flag=libwarpx.get_Ex_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def EyFPWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=1,
-                            get_lovects=_libwarpx.get_mesh_electric_field_fp_lovects,
-                            get_fabs=_libwarpx.get_mesh_electric_field_fp,
-                            get_nodal_flag=_libwarpx.get_Ey_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_electric_field_fp_lovects,
+                            get_fabs=libwarpx.get_mesh_electric_field_fp,
+                            get_nodal_flag=libwarpx.get_Ey_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def EzFPWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=2,
-                            get_lovects=_libwarpx.get_mesh_electric_field_fp_lovects,
-                            get_fabs=_libwarpx.get_mesh_electric_field_fp,
-                            get_nodal_flag=_libwarpx.get_Ez_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_electric_field_fp_lovects,
+                            get_fabs=libwarpx.get_mesh_electric_field_fp,
+                            get_nodal_flag=libwarpx.get_Ez_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def BxFPWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=0,
-                            get_lovects=_libwarpx.get_mesh_magnetic_field_fp_lovects,
-                            get_fabs=_libwarpx.get_mesh_magnetic_field_fp,
-                            get_nodal_flag=_libwarpx.get_Bx_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_magnetic_field_fp_lovects,
+                            get_fabs=libwarpx.get_mesh_magnetic_field_fp,
+                            get_nodal_flag=libwarpx.get_Bx_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def ByFPWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=1,
-                            get_lovects=_libwarpx.get_mesh_magnetic_field_fp_lovects,
-                            get_fabs=_libwarpx.get_mesh_magnetic_field_fp,
-                            get_nodal_flag=_libwarpx.get_By_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_magnetic_field_fp_lovects,
+                            get_fabs=libwarpx.get_mesh_magnetic_field_fp,
+                            get_nodal_flag=libwarpx.get_By_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def BzFPWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=2,
-                            get_lovects=_libwarpx.get_mesh_magnetic_field_fp_lovects,
-                            get_fabs=_libwarpx.get_mesh_magnetic_field_fp,
-                            get_nodal_flag=_libwarpx.get_Bz_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_magnetic_field_fp_lovects,
+                            get_fabs=libwarpx.get_mesh_magnetic_field_fp,
+                            get_nodal_flag=libwarpx.get_Bz_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def JxFPWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=0,
-                            get_lovects=_libwarpx.get_mesh_current_density_fp_lovects,
-                            get_fabs=_libwarpx.get_mesh_current_density_fp,
-                            get_nodal_flag=_libwarpx.get_Jx_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_current_density_fp_lovects,
+                            get_fabs=libwarpx.get_mesh_current_density_fp,
+                            get_nodal_flag=libwarpx.get_Jx_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def JyFPWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=1,
-                            get_lovects=_libwarpx.get_mesh_current_density_fp_lovects,
-                            get_fabs=_libwarpx.get_mesh_current_density_fp,
-                            get_nodal_flag=_libwarpx.get_Jy_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_current_density_fp_lovects,
+                            get_fabs=libwarpx.get_mesh_current_density_fp,
+                            get_nodal_flag=libwarpx.get_Jy_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def JzFPWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=2,
-                            get_lovects=_libwarpx.get_mesh_current_density_fp_lovects,
-                            get_fabs=_libwarpx.get_mesh_current_density_fp,
-                            get_nodal_flag=_libwarpx.get_Jz_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_current_density_fp_lovects,
+                            get_fabs=libwarpx.get_mesh_current_density_fp,
+                            get_nodal_flag=libwarpx.get_Jz_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def RhoFPWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=None,
-                            get_lovects=_libwarpx.get_mesh_charge_density_fp_lovects,
-                            get_fabs=_libwarpx.get_mesh_charge_density_fp,
-                            get_nodal_flag=_libwarpx.get_Rho_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_charge_density_fp_lovects,
+                            get_fabs=libwarpx.get_mesh_charge_density_fp,
+                            get_nodal_flag=libwarpx.get_Rho_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def PhiFPWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=None,
-                            get_lovects=_libwarpx.get_mesh_phi_fp_lovects,
-                            get_fabs=_libwarpx.get_mesh_phi_fp,
-                            get_nodal_flag=_libwarpx.get_Phi_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_phi_fp_lovects,
+                            get_fabs=libwarpx.get_mesh_phi_fp,
+                            get_nodal_flag=libwarpx.get_Phi_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def FFPWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=None,
-                            get_lovects=_libwarpx.get_mesh_F_fp_lovects,
-                            get_fabs=_libwarpx.get_mesh_F_fp,
-                            get_nodal_flag=_libwarpx.get_F_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_F_fp_lovects,
+                            get_fabs=libwarpx.get_mesh_F_fp,
+                            get_nodal_flag=libwarpx.get_F_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def GFPWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=None,
-                            get_lovects=_libwarpx.get_mesh_G_fp_lovects,
-                            get_fabs=_libwarpx.get_mesh_G_fp,
-                            get_nodal_flag=_libwarpx.get_G_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_G_fp_lovects,
+                            get_fabs=libwarpx.get_mesh_G_fp,
+                            get_nodal_flag=libwarpx.get_G_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def EdgeLengthsxWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=0,
-                            get_lovects=_libwarpx.get_mesh_edge_lengths_lovects,
-                            get_fabs=_libwarpx.get_mesh_edge_lengths,
-                            get_nodal_flag=_libwarpx.get_edge_lengths_x_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_edge_lengths_lovects,
+                            get_fabs=libwarpx.get_mesh_edge_lengths,
+                            get_nodal_flag=libwarpx.get_edge_lengths_x_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def EdgeLengthsyWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=1,
-                            get_lovects=_libwarpx.get_mesh_edge_lengths_lovects,
-                            get_fabs=_libwarpx.get_mesh_edge_lengths,
-                            get_nodal_flag=_libwarpx.get_edge_lengths_y_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_edge_lengths_lovects,
+                            get_fabs=libwarpx.get_mesh_edge_lengths,
+                            get_nodal_flag=libwarpx.get_edge_lengths_y_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def EdgeLengthszWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=2,
-                            get_lovects=_libwarpx.get_mesh_edge_lengths_lovects,
-                            get_fabs=_libwarpx.get_mesh_edge_lengths,
-                            get_nodal_flag=_libwarpx.get_edge_lengths_z_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_edge_lengths_lovects,
+                            get_fabs=libwarpx.get_mesh_edge_lengths,
+                            get_nodal_flag=libwarpx.get_edge_lengths_z_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def FaceAreasxWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=0,
-                            get_lovects=_libwarpx.get_mesh_face_areas_lovects,
-                            get_fabs=_libwarpx.get_mesh_face_areas,
-                            get_nodal_flag=_libwarpx.get_face_areas_x_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_face_areas_lovects,
+                            get_fabs=libwarpx.get_mesh_face_areas,
+                            get_nodal_flag=libwarpx.get_face_areas_x_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def FaceAreasyWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=1,
-                            get_lovects=_libwarpx.get_mesh_face_areas_lovects,
-                            get_fabs=_libwarpx.get_mesh_face_areas,
-                            get_nodal_flag=_libwarpx.get_face_areas_y_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_face_areas_lovects,
+                            get_fabs=libwarpx.get_mesh_face_areas,
+                            get_nodal_flag=libwarpx.get_face_areas_y_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def FaceAreaszWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=2,
-                            get_lovects=_libwarpx.get_mesh_face_areas_lovects,
-                            get_fabs=_libwarpx.get_mesh_face_areas,
-                            get_nodal_flag=_libwarpx.get_face_areas_z_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_face_areas_lovects,
+                            get_fabs=libwarpx.get_mesh_face_areas,
+                            get_nodal_flag=libwarpx.get_face_areas_z_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def ExCPPMLWrapper(level=1, include_ghosts=False):
     assert level>0, Exception('Coarse patch only available on levels > 0')
     return _MultiFABWrapper(direction=0,
-                            get_lovects=_libwarpx.get_mesh_electric_field_cp_lovects_pml,
-                            get_fabs=_libwarpx.get_mesh_electric_field_cp_pml,
-                            get_nodal_flag=_libwarpx.get_Ex_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_electric_field_cp_lovects_pml,
+                            get_fabs=libwarpx.get_mesh_electric_field_cp_pml,
+                            get_nodal_flag=libwarpx.get_Ex_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def EyCPPMLWrapper(level=1, include_ghosts=False):
     assert level>0, Exception('Coarse patch only available on levels > 0')
     return _MultiFABWrapper(direction=1,
-                            get_lovects=_libwarpx.get_mesh_electric_field_cp_lovects_pml,
-                            get_fabs=_libwarpx.get_mesh_electric_field_cp_pml,
-                            get_nodal_flag=_libwarpx.get_Ey_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_electric_field_cp_lovects_pml,
+                            get_fabs=libwarpx.get_mesh_electric_field_cp_pml,
+                            get_nodal_flag=libwarpx.get_Ey_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def EzCPPMLWrapper(level=1, include_ghosts=False):
     assert level>0, Exception('Coarse patch only available on levels > 0')
     return _MultiFABWrapper(direction=2,
-                            get_lovects=_libwarpx.get_mesh_electric_field_cp_lovects_pml,
-                            get_fabs=_libwarpx.get_mesh_electric_field_cp_pml,
-                            get_nodal_flag=_libwarpx.get_Ez_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_electric_field_cp_lovects_pml,
+                            get_fabs=libwarpx.get_mesh_electric_field_cp_pml,
+                            get_nodal_flag=libwarpx.get_Ez_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def BxCPPMLWrapper(level=1, include_ghosts=False):
     assert level>0, Exception('Coarse patch only available on levels > 0')
     return _MultiFABWrapper(direction=0,
-                            get_lovects=_libwarpx.get_mesh_magnetic_field_cp_lovects_pml,
-                            get_fabs=_libwarpx.get_mesh_magnetic_field_cp_pml,
-                            get_nodal_flag=_libwarpx.get_Bx_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_magnetic_field_cp_lovects_pml,
+                            get_fabs=libwarpx.get_mesh_magnetic_field_cp_pml,
+                            get_nodal_flag=libwarpx.get_Bx_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def ByCPPMLWrapper(level=1, include_ghosts=False):
     assert level>0, Exception('Coarse patch only available on levels > 0')
     return _MultiFABWrapper(direction=1,
-                            get_lovects=_libwarpx.get_mesh_magnetic_field_cp_lovects_pml,
-                            get_fabs=_libwarpx.get_mesh_magnetic_field_cp_pml,
-                            get_nodal_flag=_libwarpx.get_By_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_magnetic_field_cp_lovects_pml,
+                            get_fabs=libwarpx.get_mesh_magnetic_field_cp_pml,
+                            get_nodal_flag=libwarpx.get_By_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def BzCPPMLWrapper(level=1, include_ghosts=False):
     assert level>0, Exception('Coarse patch only available on levels > 0')
     return _MultiFABWrapper(direction=2,
-                            get_lovects=_libwarpx.get_mesh_magnetic_field_cp_lovects_pml,
-                            get_fabs=_libwarpx.get_mesh_magnetic_field_cp_pml,
-                            get_nodal_flag=_libwarpx.get_Bz_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_magnetic_field_cp_lovects_pml,
+                            get_fabs=libwarpx.get_mesh_magnetic_field_cp_pml,
+                            get_nodal_flag=libwarpx.get_Bz_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def JxCPPMLWrapper(level=1, include_ghosts=False):
     assert level>0, Exception('Coarse patch only available on levels > 0')
     return _MultiFABWrapper(direction=0,
-                            get_lovects=_libwarpx.get_mesh_current_density_cp_lovects_pml,
-                            get_fabs=_libwarpx.get_mesh_current_density_cp_pml,
-                            get_nodal_flag=_libwarpx.get_Jx_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_current_density_cp_lovects_pml,
+                            get_fabs=libwarpx.get_mesh_current_density_cp_pml,
+                            get_nodal_flag=libwarpx.get_Jx_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def JyCPPMLWrapper(level=1, include_ghosts=False):
     assert level>0, Exception('Coarse patch only available on levels > 0')
     return _MultiFABWrapper(direction=1,
-                            get_lovects=_libwarpx.get_mesh_current_density_cp_lovects_pml,
-                            get_fabs=_libwarpx.get_mesh_current_density_cp_pml,
-                            get_nodal_flag=_libwarpx.get_Jy_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_current_density_cp_lovects_pml,
+                            get_fabs=libwarpx.get_mesh_current_density_cp_pml,
+                            get_nodal_flag=libwarpx.get_Jy_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def JzCPPMLWrapper(level=1, include_ghosts=False):
     assert level>0, Exception('Coarse patch only available on levels > 0')
     return _MultiFABWrapper(direction=2,
-                            get_lovects=_libwarpx.get_mesh_current_density_cp_lovects_pml,
-                            get_fabs=_libwarpx.get_mesh_current_density_cp_pml,
-                            get_nodal_flag=_libwarpx.get_Jz_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_current_density_cp_lovects_pml,
+                            get_fabs=libwarpx.get_mesh_current_density_cp_pml,
+                            get_nodal_flag=libwarpx.get_Jz_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def FCPPMLWrapper(level=1, include_ghosts=False):
     assert level>0, Exception('Coarse patch only available on levels > 0')
     return _MultiFABWrapper(direction=None,
-                            get_lovects=_libwarpx.get_mesh_F_cp_lovects_pml,
-                            get_fabs=_libwarpx.get_mesh_F_cp_pml,
-                            get_nodal_flag=_libwarpx.get_F_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_F_cp_lovects_pml,
+                            get_fabs=libwarpx.get_mesh_F_cp_pml,
+                            get_nodal_flag=libwarpx.get_F_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def GCPPMLWrapper(level=1, include_ghosts=False):
     assert level>0, Exception('Coarse patch only available on levels > 0')
     return _MultiFABWrapper(direction=None,
-                            get_lovects=_libwarpx.get_mesh_G_cp_lovects_pml,
-                            get_fabs=_libwarpx.get_mesh_G_cp_pml,
-                            get_nodal_flag=_libwarpx.get_G_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_G_cp_lovects_pml,
+                            get_fabs=libwarpx.get_mesh_G_cp_pml,
+                            get_nodal_flag=libwarpx.get_G_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def ExFPPMLWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=0,
-                            get_lovects=_libwarpx.get_mesh_electric_field_fp_lovects_pml,
-                            get_fabs=_libwarpx.get_mesh_electric_field_fp_pml,
-                            get_nodal_flag=_libwarpx.get_Ex_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_electric_field_fp_lovects_pml,
+                            get_fabs=libwarpx.get_mesh_electric_field_fp_pml,
+                            get_nodal_flag=libwarpx.get_Ex_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def EyFPPMLWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=1,
-                            get_lovects=_libwarpx.get_mesh_electric_field_fp_lovects_pml,
-                            get_fabs=_libwarpx.get_mesh_electric_field_fp_pml,
-                            get_nodal_flag=_libwarpx.get_Ey_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_electric_field_fp_lovects_pml,
+                            get_fabs=libwarpx.get_mesh_electric_field_fp_pml,
+                            get_nodal_flag=libwarpx.get_Ey_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def EzFPPMLWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=2,
-                            get_lovects=_libwarpx.get_mesh_electric_field_fp_lovects_pml,
-                            get_fabs=_libwarpx.get_mesh_electric_field_fp_pml,
-                            get_nodal_flag=_libwarpx.get_Ez_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_electric_field_fp_lovects_pml,
+                            get_fabs=libwarpx.get_mesh_electric_field_fp_pml,
+                            get_nodal_flag=libwarpx.get_Ez_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def BxFPPMLWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=0,
-                            get_lovects=_libwarpx.get_mesh_magnetic_field_fp_lovects_pml,
-                            get_fabs=_libwarpx.get_mesh_magnetic_field_fp_pml,
-                            get_nodal_flag=_libwarpx.get_Bx_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_magnetic_field_fp_lovects_pml,
+                            get_fabs=libwarpx.get_mesh_magnetic_field_fp_pml,
+                            get_nodal_flag=libwarpx.get_Bx_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def ByFPPMLWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=1,
-                            get_lovects=_libwarpx.get_mesh_magnetic_field_fp_lovects_pml,
-                            get_fabs=_libwarpx.get_mesh_magnetic_field_fp_pml,
-                            get_nodal_flag=_libwarpx.get_By_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_magnetic_field_fp_lovects_pml,
+                            get_fabs=libwarpx.get_mesh_magnetic_field_fp_pml,
+                            get_nodal_flag=libwarpx.get_By_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def BzFPPMLWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=2,
-                            get_lovects=_libwarpx.get_mesh_magnetic_field_fp_lovects_pml,
-                            get_fabs=_libwarpx.get_mesh_magnetic_field_fp_pml,
-                            get_nodal_flag=_libwarpx.get_Bz_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_magnetic_field_fp_lovects_pml,
+                            get_fabs=libwarpx.get_mesh_magnetic_field_fp_pml,
+                            get_nodal_flag=libwarpx.get_Bz_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def JxFPPMLWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=0,
-                            get_lovects=_libwarpx.get_mesh_current_density_fp_lovects_pml,
-                            get_fabs=_libwarpx.get_mesh_current_density_fp_pml,
-                            get_nodal_flag=_libwarpx.get_Jx_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_current_density_fp_lovects_pml,
+                            get_fabs=libwarpx.get_mesh_current_density_fp_pml,
+                            get_nodal_flag=libwarpx.get_Jx_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def JyFPPMLWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=1,
-                            get_lovects=_libwarpx.get_mesh_current_density_fp_lovects_pml,
-                            get_fabs=_libwarpx.get_mesh_current_density_fp_pml,
-                            get_nodal_flag=_libwarpx.get_Jy_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_current_density_fp_lovects_pml,
+                            get_fabs=libwarpx.get_mesh_current_density_fp_pml,
+                            get_nodal_flag=libwarpx.get_Jy_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def JzFPPMLWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=2,
-                            get_lovects=_libwarpx.get_mesh_current_density_fp_lovects_pml,
-                            get_fabs=_libwarpx.get_mesh_current_density_fp_pml,
-                            get_nodal_flag=_libwarpx.get_Jz_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_current_density_fp_lovects_pml,
+                            get_fabs=libwarpx.get_mesh_current_density_fp_pml,
+                            get_nodal_flag=libwarpx.get_Jz_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def FFPPMLWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=None,
-                            get_lovects=_libwarpx.get_mesh_F_fp_lovects_pml,
-                            get_fabs=_libwarpx.get_mesh_F_fp_pml,
-                            get_nodal_flag=_libwarpx.get_F_pml_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_F_fp_lovects_pml,
+                            get_fabs=libwarpx.get_mesh_F_fp_pml,
+                            get_nodal_flag=libwarpx.get_F_pml_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
 
 def GFPPMLWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(direction=None,
-                            get_lovects=_libwarpx.get_mesh_G_fp_lovects_pml,
-                            get_fabs=_libwarpx.get_mesh_G_fp_pml,
-                            get_nodal_flag=_libwarpx.get_G_pml_nodal_flag,
+                            get_lovects=libwarpx.get_mesh_G_fp_lovects_pml,
+                            get_fabs=libwarpx.get_mesh_G_fp_pml,
+                            get_nodal_flag=libwarpx.get_G_pml_nodal_flag,
                             level=level, include_ghosts=include_ghosts)
diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py
index cbe01baa6..49aac44dc 100644
--- a/Python/pywarpx/picmi.py
+++ b/Python/pywarpx/picmi.py
@@ -924,6 +924,12 @@ class EmbeddedBoundary(picmistandard.base._ClassWithInit):
 
 
 class Simulation(picmistandard.PICMI_Simulation):
+
+    # Set the C++ WarpX interface (see _libwarpx.LibWarpX) as an extension to
+    # Simulation objects. In the future, LibWarpX objects may actually be owned
+    # by Simulation objects to permit multiple WarpX runs simultaneously.
+    extension = pywarpx.libwarpx
+
     def init(self, kw):
 
         self.current_deposition_algo = kw.pop('warpx_current_deposition_algo', None)