Say $\Omega$ is the domain and V a FEM space over $\Omega$. Say $\Omega_i$ is a subdomain and f_i is a function with support in $\Omega_i$, i.e. f is zero outside $\Omega_i$.
I want to assemble the form that is related to the functional $f(v) := \int_\Omega f\cdot v dx$ on V. So, let v be a TrialFunction(V).
If I simply compute inner(f, v)*dx, the computed integrals may be a poor approximation since the quadrature undersamples f. This happens in particular at the boundary of $\Omega_i$ if it does not align with the cell boundaries of V.
I was thinking of defining $\Omega_i$ as a subdomain with the measure dx(1) and compute inner(f, v)*dx(1).
Thus, my question is: How does Fenics handle these cases, where the boundary of the domain of integration does not aligh with the actual cell boundaries?
Or is there a better way to avoid the undersampling described above?
EDIT
As Jan Blechta has pointed out, integration in FEniCS only considers complete cells. The following script gives an example on the interval $\Omega : = [0,1]$ divided into 5 segments, with the FEM space Y, and a function ffunc that is $1$ on $(0.1, 0.3)$ and zero elsewhere. Say '1' is the constant mesh function in Y that is always one.
If one defines the subdomain funcdom as the interval $(0.1, 0.3)$ with the measure dx(1), one obtains that inner('1', ffunc)*dx(1) is $0$. Probably because, the subdomain does not share a cell with Y.
The integral over the whole domain gives the correct value 0.2, probably, because FEniCS (mis)interpretes the function as the triangle marked in red.
import dolfin
class LocFun(dolfin.Expression):
""" a locally defined constant function"""
def __init__(self, xmin, xmax, fvalue=1):
self.xmin = xmin
self.xmax = xmax
self.fvalue = fvalue
def eval(self, value, x):
value[:] = self.fvalue if self.xmin < x[0] < self.xmax else 0
class LocFuncDom(dolfin.SubDomain):
""" domain of definition as subdomain """
def __init__(self, xmin, xmax):
dolfin.SubDomain.__init__(self)
self.xmin, self.xmax = xmin, xmax
def inside(self, x, on_boundary):
return (dolfin.between(x[0], (self.xmin, self.xmax)))
# extension of the domain where ffunc is not 0
smin, smax = 0.1, 0.3
mesh = dolfin.IntervalMesh(5, 0, 1)
# dolfin.plot(mesh)
Y = dolfin.FunctionSpace(mesh, 'CG', 1)
y = dolfin.Expression('1')
y = dolfin.interpolate(y, Y)
# dolfin.plot(y)
domains = dolfin.CellFunction('uint', Y.mesh())
domains.set_all(0)
# the support of ffunc as subdomain
ffunc = LocFun(smin, smax)
funcdom = LocFuncDom(smin, smax)
funcdom.mark(domains, 1)
dx = dolfin.Measure('dx')[domains]
# dolfin.plot(ffunc, mesh)
# dolfin.interactive(True)
# ffunc integrated over the whole domain
y1 = dolfin.inner(y, ffunc) * dolfin.dx
# ffunc integrated on over its support
y2 = dolfin.inner(y, ffunc) * dx(1)
print dolfin.assemble(y1), dolfin.assemble(y2)
