Hello to every one,
I would like to ask please,
Suppose I created a square using the “square” command, and I want to define on a certain edge of this square that on half-edge Dirichlet boundary condition, u=0, and in the other half-edge Neumann boundary condition, u_n=0.
Is it possible to define it in FreeFem++ using the “on” command?
Best Regards,
Mordechai.
Yes, but you need to change the label of the respective edges first.
Thank you prj,
Could you please send me an example how to write it in FreeFem++, if you have?
Best Regards,
Mordechai.
int[int] lab = [1,2,10,5];
mesh ThBot = square(20, 10, [x, 0.5*y], label = lab);
lab = [10,2,3,4];
mesh ThTop = square(20, 10, [x, 0.5+0.5*y], label = lab);
mesh Th = ThBot + ThTop;
Th = change(Th, rmInternalEdges = true);
fespace Vh(Th, P1);
varf Poisson(u, v) = int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v)) + int2d(Th)(v) + on(5, u = 0.0);
matrix A = Poisson(Vh, Vh);
real[int] rhs = Poisson(0, Vh);
Vh sol;
sol[] = A^-1 * rhs;
plot(sol);
Thank you prj!
I appreciate it very much
Hi prj,
Just to be sure, you have taken in your example,
-L (u)=1, where L means the Laplacian operator
then you premultiply it with a test function v and use the Divergence theorem thus, in the variational formulation we obtain:
varf Poisson(u, v) = int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v)) + int2d(Th)(v) + on(5, u = 0.0);
, right?
Thank you so much it really helped me.
Mordechai.
Right, that is correct.
Hi Prj ,
I would like to ask you please two questions,
I have written a FF++ code on the 2D wave equation in a square.
On the boundary I imposed the Neumann boundary condition, u_n=0, on three edges, and on the left bottom edge I imposed the Dirichlet boundary condition, u=0, and in the top half-edge, I imposed Neumann boundary condition, u_n=0.
I used your code example to impose these boundary conditions for the wave equation however, I don’t know why FF++ doesn’t show me these separate boundary conditions in the left edge of the square when I ask it to plot the mesh?
I have a subdomain in my original domain which is also a square, in this subdomain, I impose the Dirichlet boundary condition, u=0, on all the boundaries. (you can see this restriction of the B.C of the subdomain in the weak formulation), however, when I ask FF++ to plot the solution, u , I get a warning comment
“Warning manifold obj nb:48896 adj 97792 of dim =2”, in each time of the simulation. Could you please tell me how can I fix this warning? and where is the problem in my code?
Thank you,
Mordechai.
Forward problem EG.edp.edp (4.0 KB)
Your meshes ThBot and ThTop are exactly the same, this doesn’t make much sense, does it?
Thank so much prj,
Yes you are right. I changed these lines code to
int[int] lab = [1,2,10,5];
mesh ThBot = square(N,N,[x0+(xF-x0)x,y0+0.5(yF-y0)*y] , label = lab);
lab = [10,2,3,4];
mesh ThTop = square(N,N,[x0+(xF-x0)x,y0+0.5(yF-y0)y+0.5pi], label = lab);
mesh Th = ThBot + ThTop;
Th = change(Th, rmInternalEdges = true);
Now it works and doesn’t give me the warning.
you made my day thanks again!
Hi, prj
I have written a code in FF++ to generate a mesh with a hole using the “trunc” command in FF++:
real x0=0,y0=0;
real xF=pi,yF=pi;
int N =2^7;
real xS=48pi/N,yS=48pi/N;
real xSF=80pi/N,ySF=80pi/N;
// domain of problem
int[int] lab = [1,2,10,5];
mesh ThBot = square(N,N,[x0+(xF-x0)x,y0+0.5(yF-y0)*y] , label = lab);
lab = [10,2,3,4];
mesh ThTop = square(N,N,[x0+(xF-x0)x,y0+0.5(yF-y0)y+0.5pi], label = lab);
mesh Th = ThBot + ThTop;
Th = change(Th, rmInternalEdges = true);
mesh Thwithole = trunc(Th, (x<xS | x>xSF|y<yS) | (x<xS | x>xSF| y>ySF ));
I would like to ask please, for the mesh with the hole I have an inner boundary and outer boundary:
How can I ask FF++ to assign the same label number for the four edges in the inner boundary?
Is “lab” is the label number assigned to the edges in the outer boundary?
Thank you,
Mordechai.
label = what you want in the trunc functionHi,
prj thank you for your reply,
for the mesh with the hole I have an inner boundary and an outer boundary:
real x0=0,y0=0;
real xF=pi,yF=pi;
int N =2^7;
real xS=48 pi/N,yS=48 pi/N;
real xSF=80 pi/N,ySF=80 pi/N;
// domain of problem
int[int] lab = [1,2,10,5];
mesh ThBot = square(N,N,[x0+(xF-x0) x,y0+0.5 (yF-y0)*y] , label = lab);
lab = [10,2,3,4];
mesh ThTop = square(N,N,[x0+(xF-x0) x,y0+0.5 (yF-y0) y+0.5 pi], label = lab);
mesh Th = ThBot + ThTop;
mesh Thwithole = trunc(Th, (x<xS | x>xSF|y<yS) | (x<xS | x>xSF| y>ySF ));
Because I want to keep the label of the outer boundary the same as it was before using the “trunc” command.
// Mesh
mesh Th1 = square(10, 10);
mesh Th2 = square(20, 10, [x+1, y]);
int[int] r1=[2,0];
plot(Th1, wait=true);
Th1 = change(Th1, label=r1); //change the label of Edges 2 in 0.
plot(Th1, wait=true);
// boundary label: 1 -> 1 bottom, 2 -> 1 right, 3->1 top, 4->1 left boundary label is 1
int[int] re=[1,1, 2,1, 3,1, 4,1]
Th2=change(Th2,refe=re);
plot(Th2,wait=1) ;
But I am not sure if this is what you meant for the labeling order?
Best Regards,
Mordechai.
Just call trunc twice, with a different label each time.
Hi,
Thank you so much prj!
I add the changes to my code (lines 40-44),
If I put in comment lines (48-127) and plot the mesh it works.
I hope these lines that I add do what you meant.
I am worried that line 41 gives the same label to all the edges on the outer boundary, wherein I have different B.C, I have at the bottom of the left half-edge Dirichlet B.C, u=0 and in other edges Neumann condition, u_n=0.
Thus, in the weak formulation if I use the command "on(outerboundary ,u=0) ", I want it to be only at the bottom of the left half-edge and not for all the edges in the outer boundary.
I have another problem,
I wrote two codes, the first is for the forward problem of the 2D wave equation defined on the square wherein I have a subdomain which is also a square. The forward problem gives the solution to the 2D wave equation at each point of the mesh.
I save the values of the solution u on each point of the mesh as a finite element function in each time of the simulation using “ofstream” command, these values of u are given to the second code using the “ifstream” command. (there is no problem in the forward code)
My second code is the backward problem for the 2D wave equation wherein I have a mesh with a hole, so when I use the “ifstream” command to read the solution u from the forward problem, FF++ sends en error that it is impossible to read the data from the forward problem to the backward problem since that the finite element functions that are defined in the backward problem have a different amount of degrees of freedom than the finite element functions that are defined in the forward problem.
This true because in the forward problem the finite element functions are defined on a finite element space whose mesh is without a hole, whereas in the backward problem the element space is defined on a mesh with a hole so there are less DOF in the backward code.
can I fix it in some way? in order to be able to read the data from the forward problem the backward problem.
Forward problem EG.edp.edp (2.3 KB) Backward EG.edp (4.0 KB)
Thank you so much for your help!
Best Regards,
mordechai.
Unless I’m mistaken, the data you save after the forward problem are for the full two-domain square, while you expect to read on the square with the hole. The meshes are different in both case right? Thus the fespaces and degrees of freedom are also different…
Hi, Julien thanks for your reply!
yes, the meshes are different and thus the fespaces and the DOF of the finite element functions are different.
So you are suggesting that for the forward problem, I should also define a mesh with a hole, and thus the fespace and finite element functions would have the same DOF?
Best Regards,
Mordechai.
Yes. depending on the efficiency of the code your looking for…
Alternatively, in your backward code you can load the whole full mesh and the corresponding solution. And then interpolate on the mesh with hole. Like in this pseudo-code
mesh ThWithHole …;
mesh ThNoHole = readmesh(“full.msh”);
fespace VhFull(ThNoHole P1);
fespace Vh(Thwithole,P1);
VhFull r1READ;
{ifstream s1249(“Results wave at iteration=1249 for N=128.txt”);
s1249>>r1READ;
r1=r1READ;
}
the last line interpolates from full mesh to mesh with hole