Extraction of Neumann boundary

[DanielBoigk]

Hello all,

My Problem is: I want to calculate both dirichlet-to-neumann maps and neumann-to-dirichlet maps. However I fail to reconstruct the values of the neumann boundary from the solution. If you guys know how the problem gets the neumann boundary in the first place, maybe you can give the hint on how to extract a neumann boundary properly. Basically for defining Neumann boundary conditions I do something like this:

Ω = Triangulation(mesh)
dΩ = Measure(Ω,2)
Γ = BoundaryTriangulation(mesh, tags= "boundary" )
dΓ = Measure(Γ,2) 
# Weak problem:
a(u, v) = ∫( γ * ∇(v) ⋅ ∇(u) )dΩ
l(v) = ∫( b * v )dΓ  #here it gets the Neumann condition

But there is no extract neumann boundary function. And sofar all my experiments on extracting it have failed. My input differs from my output. Like this:

I know how to calculate the gradient at a given point, but I don't know exactly how the "normal" vector for each boundary point is calculated.

[JordiManyer]

I would try something along the lines of

uh = solve(...)
n = get_normal_vector(Γ)

u_N = ∇(uh)⋅n

which returns a CellField on Γ that you can evaluate on your boundary.

[DanielBoigk]

My solution was to assemble the Neumann problem with:

op =AffineFEOperator(a,l,U,V)

and then get the Matrix with

K = op.op.matrix
uh = solve(op)
u = uh.free_values

Now i can extract the Neumann boundary was defined in the problem with:

f_out = K * u

if I select the entries corresponding to the boundary points I can get value that was put in out.

The Problem however was the following: Even though the $\gamma$ function integrated to zero in the analytical sense it didn't sum up to zero on the points where the boundary was evaluated. Thus the system became overdetermined. Is there a way to make enforce sum zero as an option in Gridap? While giving the boundary as a function is convenient it produces a faulty solution in this case, that need to be handled by explicitely editing:

f = op.op.vector

Many boundary value problems have the constraint that something sums up/integrates to zero.

[JordiManyer]

Is there a way to make enforce sum zero as an option in Gridap?

Yes, it's the :zerosum constraint (ZeroMeanFESpace). It does not work for spaces with dirichlet BCs, however. On the other hand, you can use a lagrange multiplier approach using a ConstantFESpace.