Wave equation in a square

Hi to every one,
I would like to ask please,
I written a code in FreeFem++ on the 2D wave equation in a square wherein at a certain edge I imposed on half-edge of the square Dirichlet boundary condition u=0 and in other half-edge Neumann boundary condition u_n=0, on the other three edges I imposed Neumann boundary condition u_n=0, at the end of the code I ask FF++ to plot the solution but FF++ sends me this warning,
" Warning manifold obj nb:48896 adj 97792 of dim =2" ,
how can I fix it?

I attached my code, if anyone runs it please put it in a folder because I ask FF++ to extract in each time the solution.

Best Regards,

include “getARGV.idp”
real x0=0,y0=0;
real xF=pi,yF=pi;
int n=7;
int N =2^n;// partitions of the space
real h=pi/N;// the space step size
real xI=0.15,yI=0.15;
real xIF=0.75,yIF=0.75;
real xS=1.05,yS=1.05;
real xSF=1.95,ySF=1.95;
func c1=1e-5*((x>xS)(x<xSF)(y>yS)(y<ySF));// velocity of the wave in a square region where there is an obstacle
func c2=1
((x<xS | x>xSF|y<yS) | (x<xS | x>xSF| y>ySF ));// velocity of the wave in a region without the obstacle
func c=c1+c2;
macro grad(u) [dx(u), dy(u)]//EOM
real t,dt,td,T;
dt=(0.163h);// step in time
td=floor(3/dt);// time duration
real idt2=1/(dt^2);
real idt1=1/(2

// domain of problem
int[int] lab = [1,2,10,5];
mesh ThBot = square(N,N,[x0+(xF-x0)*x,y0+(yF-y0)*y] , label = lab);
lab = [10,2,3,4];
mesh ThTop = square(N,N,[x0+(xF-x0)*x,y0+(yF-y0)*y], label = lab);
mesh Th = ThBot + ThTop;
Th = change(Th, rmInternalEdges = true);
mesh ThI=square(N,N,[xI+(xIF-xI)*x,yI+(yIF-yI)*y] );// subdomain for the I.C
int [int] lab1=[6,7,8,9];
mesh ThS=square(N,N,[xS+(xSF-xS)*x,yS+(ySF-yS)*y],label=lab1 );// subdomain for the obstacle
Th=Th+ThS;// to include the subdomain
// Finite element space
fespace Vh(Th,P1);
Vh u,v,uold,uvold,v0,g,f,k;
fespace Ph(Th,P0);
Ph cvelocity = c;
plot(cvelocity,cmm=" velocity in the domain", wait=1, fill=1,value=1);

// functions to define Hermite cubic bell for the I.C
g=(1-15x+66.6667x^2-74.0741x^3)(x>=0.15 && x<=0.45)+(-12.5+75x-133.333x^2+74.0741x^3)(x>=0.45 && x<=0.75);
k=(1-15y+66.6667y^2-74.0741y^3)(y>=0.15 && y<=0.45)+(-12.5+75y-133.333y^2+74.0741y^3)(y>=0.45 && y<=0.75);

// initialisation
f=gk((x>xI)(x<xIF)(y>yI)*(y<yIF));// I.C: Hermite cubic bell on ThI
plot(Th,uvold,cmm= “initial time u0”,wait=true,fill=1,dim=2,value=true);

// weak formulation and a loop on the time for a given mesh
for(int s=0;s<=td;s++)
T =sdt;
problem wave(u,v)=int2d(Th)(idt2
+on(5,6,7,8,9, u = 0);
//initial conditions
verbosity = getARGV("-verbosity", 0);
uold =u;
plot(Th,u,cmm= “Wave t=”+T+", N = “+N,fill=1,dim=2,value=true);
ofstream s(“Results wave at iteration=”+s+” for N="+N+".txt");
s<<u[]<< endl;

your code don’t work. I think formatting is messing things up (see e.g func c). try to upload the file instead of copy paste…

p.s. I think simply adding

Th=Th+ThS;// to include the subdomain

is strange it superimposes the two grids…
have a look to the documentation
fig 59 seems to do correctly what you’re after

Hi julien,
Thank you for your comments , I attached the code
please put it in a folder before you run it.

I used the “+” operator in the following line code
Since I want to simulate the square on a uniform mesh. If I use the “border” and the “buildmesh” commands, then the mesh will not be uniform as it using the square command.

I have changed a little bit the code, I defined the obstacle region using the points of the mesh, and similarly, I defined the region of the compact support for the initial condition to be defined on points of the mesh, and change the I.C to a simpler one, I taught that doing so it would avoid the warning line that FF++ sends me at each time of the simulation, unfortunately it keeps sending me this warning.

Best Regards,

Forward problem EG.edp.edp (2.1 KB)