- edited description
Exterior facet domains integrals only working when used in order
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.
Comments (6)
-
Account Deleted reporter -
I assume this is in parallel? How do you deal with the "internal boundary" between PEs when defining the subdomain?
-
Account Deleted reporter You assume correctly. I don't deal with internal boundaries in any particular way: in subclassing
SubDomain
I overridebool inside(const real*, bool)
, just like we always do, and letdolfin
take care of everything. Is this wrong? -
You have to make sure that inside doesn't give back false positives.
-
Account Deleted reporter What is a false positive? And how is this related to swapping the numbers and the markers?
-
Account Deleted reporter 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 thatufc.form.num_cell_integrals()
will be1
butconst uint domain = (*domains)(cell)
will be2
and thenintegral = 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 offor (uint i = 0; i < form.num_cell_integrals(); i++)
which should befor i in list of used indexes
. Any suggestion? - Log in to comment