How do I make a UFL object acceptable to interpolate?

7 weeks ago by
I'm trying to make the following function work:
 ​def evenodd_functions(
    """Break a function into even and odd components

    Required parameters:
    omesh: the mesh on which the function is defined
    degree: the degree of the FunctionSpace 
    func: the Function. This has to be something that fe.interpolate
        can interpolate onto a FunctionSpace or that fe.project can
        project onto a FunctionSpace.
    width: the width of the domain on which func is defined. (If not
        provided, this will be determined from omesh.
    evenodd: the symmetries of the functions to be constructed
        evenodd_symmetries(dim) is used if this is not provided
    SS = FunctionSpace(omesh, 'CG', degree)
    dim = omesh.geometry().dim()
    if width is None:
        stats = mesh_stats(omesh)
        width = stats['xmax']
    if evenodd is None:
        evenodd = evenodd_symmetries(dim)
        f0 = fe.interpolate(func, SS)
    except TypeError:
        f0 = fe.project(func, SS)
    ffuncs = []
    flips = evenodd_symmetries(dim)
    for flip in (flips):
        fmesh = Mesh(omesh)
        SSf = FunctionSpace(fmesh, 'CG', degree)
        ffunc = fe.interpolate(f0, SSf)
        fmesh.coordinates()[:, :] = (width*flip
                                     + (1 - 2*flip)*fmesh.coordinates())
    E = evenodd_matrix(evenodd)
    components = matmul(2**(-dim)*E, ffuncs)
    cs = []
    for c in components:
            cs.append(fe.interpolate(c, SS))
        except TypeError:
            cs.append(fe.project(c, SS, solver_type='lu'))

matmul is just a general matrix multiplication function: it takes a matrix as its first argument and a list of whatevers as its second and produces a list of appropriate sums of products of numbers with whatevers. In this case the whatevers are FEniCS Functions. Thus components is a list of ufl.algebra.Sum's. For instance, here is repr(components[0], True):

Sum(Sum(Sum(Product(FloatValue(0.25), Coefficient(FunctionSpace(Mesh(VectorElem\
ent(FiniteElement('Lagrange', triangle, 1), dim=2), 74), FiniteElement('Lagrang\
e', triangle, 3)), 78)), Product(FloatValue(0.25), Coefficient(FunctionSpace(Me\
sh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 81), FiniteEle\
ment('Lagrange', triangle, 3)), 85))), Product(FloatValue(0.25), Coefficient(Fu\
nctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), \
88), FiniteElement('Lagrange', triangle, 3)), 92))), Product(FloatValue(0.25), \
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle\
, 1), dim=2), 95), FiniteElement('Lagrange', triangle, 3)), 99)))

str(components[0]) is more readable: it is '0.25 * f_78 + 0.25 * f_85 + 0.25 * f_92 + 0.25 * f_99'

OK, so that's all fine. The part I don't understand is the last thing in the function. When I try to interpolate this thing onto the CG FunctionSpace, I get a TypeError: in method 'Function_interpolate', argument 2 of type 'dolfin::GenericFunction const &'

OK, fine, so a ufl object is not acceptable to interpolate. This TypeError is caught and I try to project onto the FunctionSpace. That sometimes works, but sometimes fails with this error, which is not recoverable:


*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
*** -------------------------------------------------------------------------
*** Error:   Unable to successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 92 (See\
ation/linearsolvertable.html for possible LU and Cholesky solvers).
*** Where:   This error was encountered inside /feedstock_root/build_artefacts/\
*** Process: 0
*** DOLFIN version: 2017.1.0
*** Git changeset:  70e77054c0a013e3a47e7161fa01ca1b83930c26
*** -------------------------------------------------------------------------
So, here is my question. How do I make the ufl.algebra.Sum object into something I can interpolate? (If you care to answer by telling me of a different but far better approach to the problem, I will be equally happy with that.)

Muchas gracias.
Community: FEniCS Project

2 Answers

6 weeks ago by
I have actually solved this problem, but in a completely different way. I dug into the DofMap and figured out how to map the transformed DOFs to the original ones. Fortunately, for the simple meshes and transformations I'm using, there is a one-to-one correspondence. Finding it is a bit painful because of floating point error, but can be done. That done, I can do everything I need to efficiently in numpy, as long as I'm content to do it in a single process initially. Here's some code. The principle function is coord_remap.

def dtype2mpitype(dtype):
    """convert a numpy dtype to an mpi4py type

    This is based on an undocumented private MPI element, so it is
        mpitype = MPI.__TypeDict__[dtype.char]
    except AttributeError:
        mpitype = MPI._typedict[dtype.char]

def gather_array(array):
    """Concatenate ndarrays from different processes

    Required parameter:
    array -- the local segment of the array

    array is a one-dimensional ndarray, that differs from process to
    process. gather_array gathers all the arrays and returns an single
    one-dimensional array whose contents are the process arrays
    concatenated in rank order.
    comm = MPI.COMM_WORLD
    size = comm.size
    type = dtype2mpitype(array.dtype)
    outshape = (-1,) + array.shape[1:]
    farray = array.flatten()
    n = farray.size
    total = comm.allreduce(n)
    ns = np.empty(size, dtype=int)
    comm.Gather(np.array(n), ns)
    displs = np.zeros(size + 1, dtype=int);
    displs[1:] = ns
    displs = np.cumsum(displs)
    displs = tuple(displs[:-1])
    farrays = np.zeros(total, dtype=array.dtype)
    if comm.rank == 0:
        recvbuf = [farrays, tuple(ns), displs, type]
        recvbuf = [farrays, tuple(ns), displs, type]
    sendbuf = [farray, n]
    comm.Gatherv(sendbuf, recvbuf)
    arrays = farrays.reshape(outshape)

def gather_dof_coords(fs):
    """Gather the coordinates of the DOFs in the FunctionSpace."""
    dim = fs.mesh().geometry().dim()
    lcoords = np.reshape(fs.tabulate_dof_coordinates(), (-1, dim))
    gcoords = gather_array(lcoords)

def integerify_transform(array, tol=None):
    """Find spacing and base for transformation to integer."""
    if tol is None:
        tol = 1e-7
    xs = np.unique(array.flatten())
    dxs = np.unique(np.diff(xs))
    s1 = dxs.searchsorted(tol)
    s2 = dxs.searchsorted(dxs[s1] + tol, side='right')
    spacing = np.mean(dxs[s1:s2])
    bases = np.unique(tol + xs/spacing - np.floor(tol + xs/spacing)) - tol
    base = np.mean(bases)
    if (np.max(np.abs(bases - base)) > tol):
        raise KSDGException(
            "couldn't transform coordinates to integer"
    return((spacing, base))

def integerify(
    """Map elements of a numpy array to integers

    Required argument:
    array: the numpy array. This will usually have elements of some
    floating point type.

    Optional keyword argument:
    transform: This is a two element list (or tuple or array). The
    first element is the minimum spcaing, and the second the
    base. Elements of array are transformed by array/spacing - base,
    then rounded to integers
    If not provided, integerify_transform(array) is usewd.
    if transform is None:
        transform = integerify_transform(array, tol=tol)
    spacing = transform[0]
    base = transform[1]
    iarray = np.empty_like(array, dtyp=int)
    np.rint(array/spacing - base, iarray, casting='unsafe')

def local_dofs(fs):
    """List the dofs local to this process
    Required parameter:
    fs: the FunctionSpace to remap
    ksdg: the KSDGSolver object whose VS FunctionSpace to map.
    Returns an ndofs x dim+3 float array with one row for each local
    dof. The last element of the row is the global dof number. The
    first (i.e., 0th) is the global cell number in which this dof is
    located. The second is the number of the vector component
    associated with this dof. The remaining elements (2 through -2)
    are the coordinates of the dof.
    if not isinstance(fs, FunctionSpace):
        ksdg = fs
        fs = ksdg.sol.function_space()
    mesh = fs.mesh()
    dim = mesh.geometry().dim()
    fss = fs.split()
    if not fss:                 # scalar FunctionSpace
        fss = [fs]
    dofmap = fs.dofmap()
    dofmaps = [fsi.dofmap() for fsi in fss]
    # Now build a numpy array with the required info
    owneddofs = dofmap.ownership_range()
    nlocal = owneddofs[1] - owneddofs[0]
    doflist = np.zeros((nlocal, dim + 3), dtype=float) # float, for simplicity
    ltg = dofmap.tabulate_local_to_global_dofs() # local to global map
    doflist[:, -1] = ltg[:nlocal] # record global dof #s in last column
    gtl = np.argsort(ltg[:nlocal]) # global to local map
    for n,dofmapn in enumerate(dofmaps):
        doflist[gtl[dofmapn.dofs()-owneddofs[0]], 1] = float(n)
    cindex = mesh.topology().global_indices(dim)
    dofspercell = dofmap.cell_dofs(0).size
    logGATHER('dofmap.cell_dofs(0)', dofmap.cell_dofs(0))
    ncells = len(cindex)
    # Here we have a little problem. Some of the cells in our list may
    # contain dofs not owned by this process. I haven't found any list
    # in the dofmap or the mesh data structures that contains only
    # cells whose dofs are exclusive to this process. So we go with
    # EAFP and just try it and catch the IndexError if we get a dof we
    # don't own. This of course leaves open the possibility that some
    # dofs don't get a cell assigned to them, which would be a Very
    # Bad Thing. So we set everyone's cell identity to -1 first,
    # allowing us to catch such failures.
    doflist[:, 0] = -1          # to catch failure to find a cell
    for cell in range(ncells):
            doflist[dofmap.cell_dofs(cell), 0] = cindex[cell]
        except IndexError as ie:
            # would probably be safe to just break out of the loop here
            logGATHER('IndexError: cell,dofmap.cell_dofs(cell)',
    # encode coordinates as reverse cumulative sum, so sort order is
    # robust to floating point error
    coords = np.reshape(fs.tabulate_dof_coordinates(), (-1, dim))
    doflist[:, 2:-1] = coords

def remap_list_sort(ipdofs, igdofs, pdofs=None, gdofs=None, tol=1e-7):
    """Utility function to make remapping from lists"""
    psort = np.lexsort(np.rot90(ipdofs))
    gsort = np.lexsort(np.rot90(igdofs))
    if (pdofs is not None and gdofs is not None):
        error = np.abs(pdofs[psort] - gdofs[gsort])
        if (error > tol).any():
            raise KSDGException("discrepancy in DOf remapping")
    remap = psort[np.argsort(gsort)]

remap_list = remap_list_sort

def coord_remap(
    """Find a DOF remapping that effects a coordinate transformation.

    Required parameteters:
    fs: The FunctionSpace to be remapped
    new_coords: The coordinates of the remapped DOFs. This is a numpy
    matrix of shape (ndofs, dim), where ndofs is the number of degrees
    of freedom. 

    remap: an array of global dofs indexes such that dof i in the
    original FunctionSpace should be mapped to remap[i]. The remapping
    doesn't respect cell identities.
    old_coords = gather_dof_coords(fs)
    transform = integerify_transform(old_coords)
    spacing = transform[0]
    base = transform[1]
    ldofs = local_dofs(fs);
    pdofs = gather_array(ldofs)
    pdofs = pdofs[np.argsort(pdofs[:,-1])]
    pdofs = pdofs[:, 1:-1]      # delete cell no and global index cols
    ndofs = pdofs.copy()
    ndofs[:, 1:] = new_coords
    ipdofs = np.empty_like(pdofs,dtype=int)
    np.rint(pdofs[:, :1], out=ipdofs[:, :2], casting='unsafe')
    np.rint(pdofs[:, 1:]/spacing - base, out=ipdofs[:, 1:],
    indofs = np.empty_like(ndofs,dtype=int)
    np.rint(ndofs[:, :1], out=indofs[:, :1], casting='unsafe')
    np.rint(ndofs[:, 1:]/spacing - base, out=indofs[:, 1:],
    return(remap_list(ipdofs, indofs, pdofs, ndofs))
6 weeks ago by
I see that you've already identified the first workaround, but since I was just hunting around for exactly what you found, I wanted to highlight it as a potential solution to this problem.

Below is a minimal case that, as a total FEniCS newbie, I'd have thought would work. Adding in an extra `project` is helping my use cases:
from dolfin import *

mesh = UnitIntervalMesh(64)
vspace = FunctionSpace(mesh, 'P', 1)
two = Constant(1) + Constant(1)

# fails on FEniCS 2018.1.0 with:
# AttributeError: 'Sum' object has no attribute '_cpp_object'
interpolate(two, vspace)

# works for me, possibly silly?
interpolate(project(two, vspace), vspace)​
Yes. Unfortunately, although the project sometimes works, it sometimes fails with PETSc Error I transcribed.
written 6 weeks ago by Leon Avery  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »