### How do I make a UFL object acceptable to interpolate?

114
views
0
7 weeks ago by
I'm trying to make the following function work:
 ​def evenodd_functions(
omesh,
degree,
func,
width=None,
evenodd=None
):
"""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)
try:
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())
fmesh.bounding_box_tree().build(fmesh)
ffuncs.append(ffunc)
E = evenodd_matrix(evenodd)
components = matmul(2**(-dim)*E, ffuncs)
cs = []
for c in components:
try:
cs.append(fe.interpolate(c, SS))
except TypeError:
cs.append(fe.project(c, SS, solver_type='lu'))
return(cs)


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:

RuntimeError:

*** -------------------------------------------------------------------------
*** 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 http://www.mcs.anl.gov/petsc/document\
ation/linearsolvertable.html for possible LU and Cholesky solvers).
*** Where:   This error was encountered inside /feedstock_root/build_artefacts/\
fenics_1499695155604/work/dolfin-2017.1.0/dolfin/la/PETScKrylovSolver.cpp.
*** 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
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
fragile.
"""
try:
mpitype = MPI.__TypeDict__[dtype.char]
except AttributeError:
mpitype = MPI._typedict[dtype.char]
return(mpitype)

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]
else:
recvbuf = [farrays, tuple(ns), displs, type]
sendbuf = [farray, n]
comm.Gatherv(sendbuf, recvbuf)
arrays = farrays.reshape(outshape)
return(arrays)

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)
return(gcoords)

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(
array,
transform=None,
tol=None
):
"""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
OR
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):
try:
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)',
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
return(doflist)

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)]
return(remap)

remap_list = remap_list_sort

def coord_remap(
fs,
new_coords
):
"""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.

Returns:
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:],
casting='unsafe')
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:],
casting='unsafe')
return(remap_list(ipdofs, indofs, pdofs, ndofs))
​
0
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