I am currently working on a simulation using FreeFem, and I am facing difficulties in implementing a nonlinear boundary condition for a PDE I am solving. Specifically, I am trying to impose a boundary condition that is dependent on the solution to the equation at the boundary, which makes the problem nonlinear.
Here’s the outline of the problem:
I am solving a second-order elliptic PDE for an unknown function u in a 2D domain using the finite element method.
The boundary condition I need to implement is of the form: ∂u/∂n = f(u) on the boundary, where f(u) is a known nonlinear function of the solution u.
In my current approach, I am trying to define the boundary condition using the on keyword for the boundary and then applying the nonlinear term. However, I am unsure about the correct syntax and handling of the nonlinear boundary condition, especially in terms of how to update the boundary condition during each iteration of the solver.
I have tried a few variations of the boundary condition and reformulations of the problem, but I am not getting the expected results. My guess is that I might be misinterpreting how to handle nonlinear tableau boundary conditions within FreeFem.
Has anyone worked with nonlinear boundary conditions in FreeFem before? Any suggestions or examples would be greatly appreciated. Additionally, if there’s a better approach for implementing this kind of condition in FreeFem, I would love to hear about it.
Hello,
Here is a basic example for a problem with nonhomogeneous Neuman BC: u-\Delta u=f in \Omega, \partial u/\partial n=g on \partial\Omega,
with \Omega=(-1,1)\times(-1,1), f=xy, g=|x+y|-1 and the exact solution u_{ex}=xy.
The variational formulation is \int_\Omega uv+\int_\Omega\nabla u\cdot\nabla v-\int_\Omega fv-\int_{\partial\Omega}gv=0, for all test functions v.
int nn=30;
mesh Th=square(nn,nn,[2.*x-1.,2.*y-1.]);
fespace Vh(Th,P1);
Vh u,v;// u: unknown, v: test function
Vh f,g;// f: rhs, g:nonhommogeneous Neumann BC
f=x*y;
g=abs(x+y)-1.;
Vh uex;
uex=x*y;
solve laplacenhN(u,v)=
int2d(Th)(u*v)
+int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))
-int2d(Th)(f*v)
-int1d(Th)(g*v)
;
real error=sqrt(int2d(Th)((u-uex)^2));
cout << "error " << error << endl;
plot(u,wait=1,value =1);
Then for your problem you could consider a sequence of approximations u_k,
and solve the problem for u_{k+1} with boundary condition \partial u_{k+1}/\partial n=F(u_k). I hope that F is nonincreasing in order to have a well-posed problem.
If you want to use the Newton method you have to solve \partial u_{k+1}/\partial n=F(u_k)+F'(u_k)(u_{k+1}-u_k)