How to integrate an `Expression` over spatial domain?

1

At right-hand-side of my Poisson equation I have expression parametrized with two parameters: A and B. For a given A I need to select such B that my equation is compatible with Neumann BC.

I think the simplest way to do so is to integrate the expression over the spatial domain with A, B = 0, 1 and A, B = 1, 0 and calculate ratio of the results. How may I perform the integration with FEniCS?

A toy example of my code is below (in the real example I have a 3D mesh of complex geometry):

from fenics import *

# Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 1)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

f = Expression('A * exp(-0.5 * ((x[0] - X)*(x[0] - X) + (x[1] - Y)*(x[1] - Y) / S2)) - B', degree=1, X=0.5, Y=0.5, S2=0.04, A=1, B=0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx

# Compute solution
u = Function(V)
solver = KrylovSolver("cg", "ilu")
solver.parameters["maximum_iterations"] = 20000
solver.parameters["absolute_tolerance"] = 1E-8
solver.parameters["monitor_convergence"] = True
solver.solve(assemble(a), u.vector(), assemble(L))
2

Hey,

you should be able to do so with the following code

f.A = 1.0
f.B = 0.0
value_1 = assemble(f*dx)

f.A = 0.0
f.B = 1.0
value_2 = assemble(f*dx)

and that’s it. With this, you have both values and can work with them.

1 Like
3

@plugged I have (unsuccessfully) tried the expression f*dx before:

In [11]: value_1 = assemble(f*dx)
This integral is missing an integration domain.
---------------------------------------------------------------------------
UFLException                              Traceback (most recent call last)
<ipython-input-11-ca62c3a200cf> in <module>()
----> 1 value_1 = assemble(f*dx)

/usr/local/lib/python3.6/dist-packages/ufl/measure.py in __rmul__(self, integrand)
    439                 domain, = domains
    440             elif len(domains) == 0:
--> 441                 error("This integral is missing an integration domain.")
    442             else:
    443                 error("Multiple domains found, making the choice of integration domain ambiguous.")

/usr/local/lib/python3.6/dist-packages/ufl/log.py in error(self, *message)
    170         "Write error message and raise an exception."
    171         self._log.error(*message)
--> 172         raise self._exception_type(self._format_raw(*message))
    173 
    174     def begin(self, *message):

UFLException: This integral is missing an integration domain.
4

You should define the measure before, try using

dx = Measure('dx', mesh)

If you have subdomains (given by a meshfunction, for example) then you can specify them here via the subdomain_data keyword.

1 Like