Mixing fenics with fortran example not working

Question

The example found in http://hplgit.github.io/fenics-mixed/doc/web/index.html on the coupling of incompressible Navier-Stokes and Windkessel models (for a network representing a cerebral aneurysm) is not working.
Computed fluxes are always zero. Even trying to compute the area of outlets returns zero.
Any indications?
I am using version 2016.1.0.
Thanks in advanced,
Lucas.

Answer 1

Hi Lucas,

I experienced the same issue when testing this code for my project. Every boundary integral returned zero, which shouldn't be the case.

I see that it has been some time since you've asked this question, but you or others might still be helped with this solution. The problem seems to be that the boundary parts are not imported from the mesh after you execute "mesh = Mesh("aneurysm.xml.gz")".
You should add the following lines:

#import mesh 
mesh = Mesh("aneurysm.xml.gz")

#obtain boundary parts of the mesh
D = mesh.topology().dim()
boundaries = MeshFunction("size_t", mesh, D - 1, mesh.domains())

ds = ds(subdomain_data = boundaries)

You can now test the scipt by

#Define function space
Q = FunctionSpace(mesh, "Lagrange", 1)

#compute area of surface with identity 1 and 2
area_1 = assemble(interpolate(Constant(1), Q) * ds( 1 ))
area_2 = assemble(interpolate(Constant(1), Q) * ds( 2 ))
print(area_1)
print(area_1)

It should now return the area of the surface.

After the velocity solution is obtained, the fluxes can be calculated as follows, with u_ being the velocity solution.

n = FacetNormal(self.mesh)
flux1 = inner(u_, n)*ds(1)
flux2 = inner(u_, n)*ds(2)

Q1 = assemble(flux1)
Q2 = assemble(flux2)

Hope this helps,

Martijn