dolfin team mailing list archive

[David Maxwell <question133939@xxxxxxxxxxxxxxxxxxxxx>]

New question #133939 on DOLFIN:
https://answers.launchpad.net/dolfin/+question/133939

I'd like to extract Neumann data from a solution to an elliptic PDE.

For example, if u is a weak solution of -Laplacian(u) = f, then its Neumann data (\partial_n u) on the boundary is well defined as follows:

Given sufficiently regular v defined on the boundary, extend v to a function in H^1 in the domain.  Then

<\partial_n u, v> = \int_\Omega \nabla u \nabla v - f v

Since v belongs to L^2 of the boundary, I can interpret \partial_n as an element of L^2 of the boundary.

In truth, I really only want the Neumann data on part of the boundary.

I thought I might proceed as follows.  In the following, 'u' is the solution of the PDE in function space V, and mesh_function is a mesh function on the edges that equals 1 on the 'top' boundary, 2 on the remainder of the boundary, and 0 on all other edges.

My first attempt went something like: 

uu = TestFunction(V)
vv = TrialFunction(V)

inside=0; top=1; bottom=2

a = uu*vv*ds(bottom)
rhs = inner(grad(u),grad(vv))*dx - vv*f*dx

top_bc = DirichletBC(V,Constant(0),mesh_function,top)
interior_bc = DirichletBC(V,Constant(0),mesh_function,inside) # I thought I was being clever here!

u_n = VariationalProblem(a,rhs,bcs=[top_bc,interior_bc],exterior_facet_domains=mesh_function).solve()

This crashes on applying the boundary conditions using 'ident', so I tried setting the 'use_ident' parameter to False on the boundary conditions.  

This still doesn't work because the interior_bc is (rightfully) too enthusiastic -- every vertex is joined to some interior edge and so I pick up the zero solution.

I tried some other tacts, but they are perhaps too foolish to describe here.

And I have the sense that I really ought to be solving a problem on a BoundaryMesh (rather than a highly constrained problem on the original mesh) but I don't know how to extend Functions on a boundary mesh to Functions on the original mesh.

Help!

David Maxwell


You received this question notification because you are a member of
DOLFIN Team, which is an answer contact for DOLFIN.