Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Extraction of Neumann boundary #1008

Open
DanielBoigk opened this issue May 18, 2024 · 3 comments
Open

Extraction of Neumann boundary #1008

DanielBoigk opened this issue May 18, 2024 · 3 comments

Comments

@DanielBoigk
Copy link

DanielBoigk commented May 18, 2024

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:
image

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
Copy link
Member

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
Copy link
Author

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
Copy link
Member

JordiManyer commented Aug 13, 2024

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants