Use integral coefficient in ufl formulation

Question

Hello everybody,

I am trying to solve a diffusion problem where the (constant) diffusivity is the result of the integration of another function over the domain.
That is, something like:
$$\int_{\Omega} k \nabla u \nabla v = \int_{\Omega} f v$$
where
$$k=\int_{\Omega} g$$
for a given $g$.

I know I can do this by writing a ufl file to compute the integral of $g$ and then passing this value along to the ufl file that solves the problem, but I would like to do it in just one shot. For example:

element = FiniteElement("Lagrange", triangle, 1)

u = TrialFunction(element)
v = TestFunction(element)
f = Coefficient(element)
g = Coefficient(element)
gg = g*dx

a = inner(gg * grad(u), grad(v))*dx
L = f*v*dx

This of course does not work.

How would I go about integrating a coefficient in a ufl code?

Thanks a lot!

Mattia

Comments

Have you managed to implement this?

---

No luck, unfortunately. I ended up using two different ufl files, one for the integration of the coefficient and one for the weak problem.

---

I'd love to see how you handled the coefficient in the weak problem, I have troubles choosing an appropriate class for the coefficient (probably Coefficient, but I couldn't use the constructor correctly).

---

I'm not sure I understand what you mean, but I declared the coefficient as a Constant, since the integration eliminates any spatial dependency of the integrand function.

Here is a MWE I put together to show you what I did:
main.cpp

#include <dolfin.h>
#include "poisson.h"
#include "integral.h"

class G : public dolfin::Expression
{
    void eval(dolfin::Array<double>& values, const dolfin::Array<double>& x) const
    {
        values[0] = cos(x[0]);
    }
};

class DirichletBoundary : public dolfin::SubDomain
{
    bool inside(const dolfin::Array<double>& x, bool on_boundary) const
    {
        return on_boundary;
    }
};

int main()
{
    dolfin::UnitSquareMesh mesh (100, 100);
    poisson::FunctionSpace V (mesh);

    dolfin::Constant zero (0.0);
    DirichletBoundary boundary;
    dolfin::DirichletBC bc (V, zero, boundary);

    poisson::BilinearForm a (V, V);
    poisson::LinearForm L (V);

    dolfin::Constant one (1.0);
    L.f = one;

    G g;
    integral::Form_integral integrator (mesh);
    integrator.g = g;
    dolfin::Constant k (dolfin::assemble(integrator));
    a.k = k;

    dolfin::Function u(V);
    dolfin::solve(a == L, u, bc);

    dolfin::plot(u);
    dolfin::interactive ();

    return 0;
}

poisson.ufl

element = FiniteElement("Lagrange", triangle, 1)

u = TrialFunction(element)
v = TestFunction(element)
f = Coefficient(element)
k = Constant(triangle)

a = inner(k * grad(u), grad(v))*dx
L = f * v * dx

integral.ufl

element = FiniteElement("Lagrange", triangle, 1)

g = Coefficient(element)
integral = g * dx

forms = [integral]

---

Thank you, that's exactly what I was looking for.

---

No problem! Glad I could help!

Answer 1

Did you try:

gg = assemble(g * dx)

Comments

Hello!

Thanks for the reply!

Unfortunately, though, that command does not seem to work:

mattia@alien:test_fenics$ ffc -l dolfin a.ufl 
This is FFC, the FEniCS Form Compiler, version 1.5.0.
For further information, visit http://www.fenics.org/ffc/.

An exception occured during evaluation of form file.
To help you find the location of the error, a temporary script
'a_debug.py'
has been created and will now be executed with debug output enabled:
Traceback (most recent call last):
  File "/usr/bin/ffc", line 213, in <module>
    sys.exit(main(sys.argv[1:]))
  File "/usr/bin/ffc", line 170, in main
    ufd = load_ufl_file(filename)
  File "/usr/lib/python2.7/dist-packages/ufl/algorithms/formfiles.py", line 192, in load_ufl_file
    namespace = execute_ufl_code(uflcode, filename)
  File "/usr/lib/python2.7/dist-packages/ufl/algorithms/formfiles.py", line 96, in execute_ufl_code
    m = __import__(basename)
  File "/home/mattia/Dropbox/polimi/test_fenics/a_debug.py", line 10, in <module>
    gg = assemble(g * dx)
NameError: name 'assemble' is not defined
mattia@alien:test_fenics$

---

Oh, sorry, I am used to using dolfin as an API for everything.
In the documentation for ufl (http://fenicsproject.org/documentation/ufl/1.0-beta2/ufl.html) they mention

class ufl.integral.Integral(integrand, measure)

so maybe try Integral(g, dx)? If that doesn't work, I am afraid you will have to wait until someone with more experience comes along.

---

Your advice seems very promising, thanks!

The only problem is (and the blame is on me) I did not specify I am using the C++ API to dolfin, so I cannot use the variables defined in my main file in the ufl formulation. Specifically, the constructor of Integral takes 7 arguments (http://fenicsproject.org/documentation/ufl/1.5.0/ufl.html#ufl.classes.Integral), among which is the integration domain, which I defined in main.cpp but cannot define in a pure ufl file.

I willl try to see what I can do with this Integral object and maybe I can come up with something.

Or did I completely misunderstood your advice?