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,
Mordechai.
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/(2dt);
// 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
uvold=f;
uold=uvold;
u=uold;
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)(idt2uv)-int2d(Th)(2idt2uoldv)+int2d
(Th)(idt2uvoldv)+int2d(Th)(c^2*grad(uold)’*grad(v))
+on(5,6,7,8,9, u = 0);
//initial conditions
verbosity = getARGV("-verbosity", 0);
wave;
uvold=uold;
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;
}