Having calculated the pressure p and velocity field v for a fluid structure interaction problem, defined over a mesh that encompasses both the fluid and the solid, I would like to find the drag acting on the solid. Previously I have done this by creating a mesh for the fluid part of the domain only (e.g. using (mesh_fluid=SubMesh(mesh, mesh_fcn, 1)) and integrating the pressure times n[0] (n is the FacetNormal) over the relevant part of the boundary of this mesh, e.g.
assemble(p*n[0]*ds(1))
However, this requires the use of SubMesh, and also requires me to project the pressure solution from the original mesh onto the (non-matching) fluid mesh. Both of these processes do not work in parallel, and so I am trying to perform the drag calculation using the original mesh.
I can define a boundary function that marks the facets corresponding the boundary between the fluid and the solid, but if I try to use the FacetNormal (n=FacetNormal(mesh)) to integrate over this boundary, I need to decide whether to use the positive or negative side of the facet, i.e. n('+') or n('-'), since the boundary is now internal to the mesh. My issue is that which side of the facet is positive and which is negative is chosen arbitrarily, and hence is not consistent over the whole of the internal boundary I am interested in, and hence if I use,
assemble(p*n('+')[0]*dS(1))
then I do not get the right answer. Is there a way of defining the FacetNormal on the internal boundary such that it will always point into or out of the fluid domain?