2d integration on cross section of 3d mesh

1

Hi everyone, I’m new from FEniCS.
I’ve already run a 3D a stationary Navier-Stokes simulation. Now, calling x the direction of propagation of the flow and y and z the transversal one, I would like to obtain the flow Q integrating the velocity on a section of my channel defined as:

Q=\int^{h}_{0}\int^{L}_{0} v(x,y,z)dydz

I’m trying to define a 2d surface and integrating with the following code:

    class Plane(SubDomain):
       def inside(self, x, on_boundary):
          return near(x[1], 1162)
    facets  = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
    Plane().mark(facets, 1)
    dS = dS(subdomain_data = facets)

    n = Constant((0.0,1.0,0.0))
    print(dot(u, n))
    flux = assemble(dot(u, n)*dS(1))`

where u is the velocity previously obtained. Unfortunately I obtain unreasonable results, can someone help me with the code?

2

I don’t know how to do it exactly but I think you should be able to do it with Paraview.
With your approach based on Subdomain you will not be able to do it except if there is a plane of element facets exactly lying at y=1162.

3

Another alternative would be to create a manual integration loop, by manually defining quadrature points and weights for the rectangular slice, and evaluate the function at these points.

4

thank you so much! I think I’ll try with the manual integration!