# Integral over internal boundary using facet normal

**URL:** <https://fenicsproject.discourse.group/t/integral-over-internal-boundary-using-facet-normal/363>\
**Category:** Uncategorized\
**Created:** [March 13, 2019, 4:24pm UTC](https://fenicsproject.discourse.group/t/integral-over-internal-boundary-using-facet-normal/363 "2019-03-13T16:24:33Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![alastair](https://avatars.discourse-cdn.com/v4/letter/a/74df32/32.png) [@alastair](https://fenicsproject.discourse.group/u/alastair)\
**Post date:** [March 13, 2019, 4:24pm UTC](https://fenicsproject.discourse.group/t/integral-over-internal-boundary-using-facet-normal/363/1 "2019-03-13T16:24:33Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [March 13, 2019, 6:25pm UTC](https://fenicsproject.discourse.group/t/integral-over-internal-boundary-using-facet-normal/363/2 "2019-03-13T18:25:18Z")

</div>

See @MiroK’s answer here:

> [@Integrating over an interior surface](https://fenicsproject.discourse.group/t/integrating-over-an-interior-surface/247/4):
>
> Hi, consider using restriction using + and -. Here foo('+') restricts foo to the positive cell and FacetNormal('+') is the outer normal of a positive cell. Following the discussion [here](https://bitbucket.org/fenics-project/dolfin/pull-requests/199/remove-facet_orientation-and-make-ds-pick/diff) the +, - side can be changed by providing a form which has a dx measure term which has specified subdomain\_data by some cell function g. Then a positive cell for a facet is the one for which value of g is higher. This allows for the following from dolfin import \* from mshr import \* omega1 = Rectangle(Point(0, 0…