Exterior facet domains integrals only working when used in order

[Former user]

I experienced this in unicorn for the computation of the drag coefficient. See commit https://bitbucket.org/fenics-hpc/unicorn/commits/0f7d94e753adf4b68d9011b51f2dbd8bbd2ca3c3 for reference.

In there, I mark subdomains like this

MeshFunction<uint> thetadrag(mesh,mesh.topology().dim()-1); // create marker
thetadrag = 1;              // mark everything with 1
body.mark(thetadrag,0);     // mark integration domain with 0

where body is an instance of a subclass of SubDomain.

In the ufl file I then have

M = dot(dot(sigma(u,P),-n),ey)*ds(0)

which works as expected. However, now I have to explicitly assign the value 1 to the whole mesh (thetadrag = 1;), which was originally initialized to 0, in order to later mark the integration domain with 0 to match ds(0).

Using instead

M = dot(dot(sigma(u,P),-n),ey)*ds(1)

(and marking accordingly) segfaults, while using

M = dot(dot(sigma(u,P),-n),ey)*ds(0) + dot(dot(sigma(u,P),-n),ey)*ds(1)

works, so it's a matter of ds(1) being present when ds(0) is absent, rather than ds(1) being used at all.

[Niclas Jansson]

I assume this is in parallel? How do you deal with the "internal boundary" between PEs when defining the subdomain?

[Former user]

You assume correctly. I don't deal with internal boundaries in any particular way: in subclassing SubDomain I override bool inside(const real*, bool), just like we always do, and let dolfin take care of everything. Is this wrong?

[Niclas Jansson]

You have to make sure that inside doesn't give back false positives.

[Former user]

What is a false positive? And how is this related to swapping the numbers and the markers?

[Former user]

I think I found the bug. Consider cell integrals for example; we have

    // UFC.h
    ufc::cell_integral** cell_integrals;

    // UFC.cpp
    cell_integrals = new ufc::cell_integral*[form.num_cell_integrals()];
    for (uint i = 0; i < form.num_cell_integrals(); i++)
        cell_integrals[i] = form.create_cell_integral(i);

    // Assembler.cpp
    const uint domain = (*domains)(cell);
    if (domain < ufc.form.num_cell_integrals())
    integral = ufc.cell_integrals[domain];

so for a form like M = f*dx(2) what happens is that ufc.form.num_cell_integrals() will be 1 but const uint domain = (*domains)(cell) will be 2 and then integral = ufc.cell_integrals[domain]; segfaults.

A solution should be to change the double pointer to an std::map<uint,ufc::cell_integrals*> but now there is a problem of how to put the correct keys in place of for (uint i = 0; i < form.num_cell_integrals(); i++) which should be for i in list of used indexes. Any suggestion?