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:
Here are the relevant snippets of the source code:
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:
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:
grad_u= ∇(u)
Then I construct points:
grid_points = [Point((b_coordinates[i][1], b_coordinates[i][2])) for i in 1:n]
and do:
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.