Dear users,

In order to obtain a partial reflective B.C. condition for a Helmholtz equation, I’d like to calculate the angle of incidence of a waterwave on a border.

I call u the solution of the equation, and the researched angle is $\gamma$.

Here is the formula (Steward - Panchang, 1999) :

$f(\gamma) = tan(\gamma) - \frac{\frac{1}{k} \frac{\partial arg(u)}{\partial x_s}}{ \frac{1}{k} \frac{\partial arg(u)}{\partial x_n} + \frac{2K_r(cos(k\beta)+K_r)}{1+2K_rcos(k\beta)+K_r^2}cos(\gamma)} $

$u$ is obtained by solving the equation a first time, with $\gamma= 0$.

I use a P1 triangulation.

Here I’d like to determine $\gamma$, knowing k, K_r, $\beta$, and $u$, buy solving $f(\gamma) = 0$ (with bracketing and disection).

x_n is the normal, and x_s is the tangent to the border, and arg(u) is the argument of u.

My problem is to calculate the derivatives of u.

I saw that I could obtain the tangent to the border by rotating the normal by pi/2.

Then on another topic, I read that I could collect the values of u on the border : getting data of specific points

So I tried to collect those values of u on the border, put them in a FE variable, take the argument of it, and apply the “dx” and “dy” operators to it. But I had an error.

Then my question is :

**Do you have any idea on how to obtain $darg(u)/dxs$ or $darg(u)/dxn$ on a specific border ?**

Here is a simple code I made to test the method :

//param

real k = 0.002; // wavenumber

real Kr = 0.4; //reflection coeff

real theta = 11*pi/8 ;
real k1 = k*cos(theta);

real k2 = k*sin(theta);

real internRadius = 2000;

real externRadius = 5000;

// interior disc

border Cint(t=0,2*pi){x= internRadius*cos(t) ; y = internRadius*sin(t) ;}

// exte disc

border CF(t=5/4.*pi, 3/2. pi){x = externRadiuscos(t) ; y = externRadius*sin(t);}

border Cext(t=-pi/2, 5/4.

*pi){x = externRadius*cos(t) ; y = externRadius*sin(t);}

mesh Th = buildmesh(Cint(-90) + CF(45) + Cext(90));

//Finite element space

fespace Vh(Th, P1);

Vh uh, vh, uinc = 0.5*exp(-1i*(k1*x+k2*y)); // uinc = incoming wave

Vh u,ur, uim; (for the real/imag solution)

Vh argu; // dNargu; // argu = argu of uh, dNargu = normal derivative ?

// Helmoltz

solve Helmoltz(uh,vh) = -int2d(Th)(dx(uh)*dx(vh) + dy(uh) dy(vh)*vh)

-k^2uh

-int1d(Th,Cint)(1i

*k*((1-Kr)/(1+Kr))

*uh*vh)

+on(CF, uh = -uinc);

//solution

// Obtain uh on the interior circle

// boundary dof

int ndof;

{ int1d(Th,Cint,qfe=qf1pE)((ndof++)*1.); }

// indexes on the boundary

int[int] i1(ndof);

{

int[int] a(2*ndof);
Vh index=0;
for(int i=0;i<Vh.ndof;i++) index[][i]=i;
int n=0;
int1d(Th,Cint,qfe=qf1pElump)((a(n++)=floor(index+0.5))*i);

*1.); Vh next=0; for(int i=0;i<ndof;i++) next[][a(2*i+1)]=a(2

i1(0)=a(1);

for(int i=1;i<ndof;i++) i1(i)=next[][i1(i-1)];

}

//display the value of uh on the boundary Cint

//for (int i=0;i<ndof;i++){

//cout<<i1(i)<<" uh("<<Th(i1(i)).x<<","<<Th(i1(i)).y<<")="<<uh[][i1(i)]<<endl;

//;

//}

//How to get darg(uh)/dn on Cint??

// plot the solution

u = real(uh);

plot(u);

Thank you all for reading,

Best