# Extracting neumann boundary in Gridap - discrepancy

**URL:** <https://discourse.julialang.org/t/extracting-neumann-boundary-in-gridap-discrepancy/114179>\
**Category:** Modelling & Simulations\
**Tags:** package, finite-element, gridap\
**Created:** [May 12, 2024, 10:12pm UTC](https://discourse.julialang.org/t/extracting-neumann-boundary-in-gridap-discrepancy/114179 "2024-05-12T22:12:54Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![DnBgk](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnbgk/32/208819_2.png) [@DnBgk](https://discourse.julialang.org/u/DnBgk)\
**Post date:** [May 12, 2024, 10:12pm UTC](https://discourse.julialang.org/t/extracting-neumann-boundary-in-gridap-discrepancy/114179/1 "2024-05-12T22:12:54Z")

</div>

So I want to model the Neumann-to-Dirichlet map in a Poisson like equation.

If I can give it my neumann boundary right away is there an implemented method that also extracts the neumann boundary values?

If no. I’m wondering about the differences between the neumann boundary I put in and the boundary I read out. Why is there a discrepancy?  
I tried different solvers but still got the same result.  
I tried a lot, but I cannot reproduce the values I put in:  
 ![neumann](https://global.discourse-cdn.com/julialang/original/3X/a/9/a9962fd1462c95505f62b80306aa99e625af02e6.png)

Here are the relevant snippets of the source code:

```julia
using Gridap
using LineSearches: BackTracking
using Plots

domain = (0.0, 1.0, 0.0, 1.0) 
n = (99, 99) 

order = 1
model = CartesianDiscreteModel(domain, n)
reffe = ReferenceFE(lagrangian, Float64, order)
V = TestFESpace(model, reffe, conformity=:H1)
U = V 

Ω = Triangulation(model)
dΩ = Measure(Ω,2)
Γ = BoundaryTriangulation(model)
dΓ = Measure(Γ,2) 

# Two very simple functions that are still giving the error:
function γ(x) 
    return 1.001 + sin(x[1]) + cos(x[2])
end
function g(x)
    return x[1] + x[2] - 1
end

# Linear Problem:
a(u, v) = ∫( γ * ∇(v) ⋅ ∇(u) )dΩ
l(v) = ∫( g * v )dΓ # Neumann condition

lop = AffineFEOperator(a,l,U,V)
lsolver = LUSolver()
ls= LinearFESolver(lsolver)
luh= solve(ls, lop)

# Same problem non linear solution:
res(u,v) = ∫( γ * ∇(v) ⋅ ∇(u) )dΩ - ∫( g * v )dΓ
jac(u,du,v) = ∫( γ * ∇(v) ⋅ ∇(du) )dΩ
op = FEOperator(res,jac,U,V)
nsolver = NLSolver(show_trace=true, method=:newton, linesearch=BackTracking())
nls = FESolver(nsolver)
luh= solve(ls, lop)

```

The result is:  
 ![neumann](https://global.discourse-cdn.com/julialang/original/3X/e/a/ea94541187e14a535b404f4beba1dc18711bfba2.png)  
The two jumps can be explained by my imperfect readout, but else it should be a perfect triangle.

Now I have various methods to readout the gradient but it doesn’t work to extract the same values. What do I need to consider?

To readout the gradient I use this:

```julia
grad_u= ∇(u)

```

Then I construct points:

```julia
 grid_points = [Point((b_coordinates[i][1], b_coordinates[i][2])) for i in 1:n]

```

and do:

```julia
Grid_Gradient = grad_u.(grid_points)

```

I then I take the points in either x or y direction from the generated variable since I’m working on a square mesh.