Replacing integral domains in form which has coefficients from QuadratureSpace

variational formulation
1

Dear Fenics Users,

I am trying to assemble a linear form on a submesh. Please consider the following code:

import dolfin as df
import numpy as np
import ufl
mesh = df.UnitIntervalMesh(8)
V = df.FunctionSpace(mesh, "CG", 1)
u = df.TrialFunction(V)
v = df.TestFunction(V)

quad_deg = 2
QE = df.FiniteElement("Quadrature", mesh.ufl_cell(), 
                       degree=quad_deg, quad_scheme="default")
Q = df.FunctionSpace(mesh, QE)

f = df.Expression("VALUE * x[0]", degree=1, VALUE=1.0)
c = df.Function(Q)
c.assign(df.Constant(1.))

L = f * v * df.ds + df.inner(df.grad(v)[0], c) * df.dx

dofmap = V.dofmap()
dof_indices = dofmap.dofs()
cell_indices = dof_indices[0::3]# cells for submesh definition

subdomain = df.MeshFunction("size_t", mesh, mesh.geometry().dim())
for ci in cell_indices:
    subdomain.set_value(ci, 78)
submesh = df.SubMesh(mesh, subdomain, 78)

V_sub = df.FunctionSpace(submesh, V.ufl_element())

# express L.arguments() as arguments of the new FunctionSpace V_sub
args = []
for arg in L.arguments():
    if arg.ufl_function_space() == V:
        arg = df.function.argument.Argument(V_sub, arg.number(), arg.part())
    args.append(arg)

# somewhat unsure about the next 3 lines:
M = ufl.replace_integral_domains(
    L(*args, coefficients={}),# just passing all arguments as *args
# do not want to change coefficients, so empty dict
    submesh.ufl_domain()
)

df.assemble(L)# works
df.assemble(M)# does not work

When I run this Code I get a RuntimeError:

RuntimeError                              Traceback (most recent call last)
~/python-projects/Fenics/snippets/replace_int_domains.py in <module>
----> 1 df.assemble(M)

~/miniconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/fem/assembling.py in assemble(form, tensor, form_compiler_parameters, add_values, finalize_tensor, keep_diagonal, backend)
    211 
    212     # Call C++ assemble
--> 213     assembler.assemble(tensor, dolfin_form)
    214 
    215     # Convert to float for scalars

RuntimeError: evaluate_basis: Function not supported/implemented for QuadratureElement.

  1. Since assembling L works, why should evaluate_basis be not implemented for QuadratureElements?

  2. Do I need to change/specify df.dx?

  3. Also when I use

c = df.Function(V)
c.assign(df.Constant(1.))

I can assemble M without error.

Can someone explain to me, what’s the problem here?

I’m using Linux Debian version 9.9 and have installed fenics via miniconda3.

In [1]: df.__version__                                                                                           
Out[1]: '2018.1.0'

Any help would be greatly appreciated.
Thanks,
Best,
Philipp