Hi,
I may be asking a stupid question, but how can I get the values of the fespace after the problem was solved. Let’s say, I implemented a diffusion pde. I can plot it with the plot function of freefem++. But how can I get the value of the function for each position in the space ? I would say I need to access every finite elements, its geometry and value. Or maybe there is a function to do that.
Mordokkai
HI,
I think, you question is a little be strange because, to get the valeur a point (x0,y0) in 2d of finite element function u
juste du u(x0,y0) (same in 3d with u(x0,y0,z0) )
Event for basis function like (sin(x)) you do not know the value at each point
this peace of code given a idea:
// value a each node of any Lagrange Finite element. here P2
mesh Th=square(10,10);
fespace Vh(Th,P2);
Vh xx=x, yy=x;
Vh u =sin(x)*cos(y);
for(int i=0; i<Vh.ndof; ++i)
cout << xx[][i] << " " << yy[][i] << " u = "<< u[][i] << endl;
Thanks for your answer.
I think, you question is a little be strange because, to get the valeur a point (x0,y0) in 2d of finite element function u
juste du u(x0,y0) (same in 3d with u(x0,y0,z0) )
That’s the first thing I tried but I got strange results.
Below you can find the result plot by FreeFem++: plot(u, wait=true, fill=true, dim=2);
Then I try to write everything in a text file:
ofstream out (“out.txt” );
out.fixed;
for(int i=0; i<dims[0]; i++){
for(int j=0; j<dims[1]; j++){
out << u(i,j) << " ";
}
out << endl;
}
And then I try to reopen it with python matplotlib:
import matplotlib.pyplot as plt
import numpy as np
mat = np.loadtxt(“out.txt”)
plt.matshow(mat)
plt.show()
And I get the following result:
If I take a basic example:
real[int] dims = [50,50];
border a(t=0, dims[0]){x=t; y=0; label=1;};
border b(t=0, dims[1]){x=dims[0]; y=t; label=2;};
border f(t=0, dims[0]){x=dims[0]-t; y=dims[1]; label=3;};
border d(t=0, dims[1]){x=0; y=dims[1]-t; label=4;};
border g(t=10, dims[0]-10){x=t; y=10; label=5;};
border h(t=10, dims[1]-10){x=dims[0]-10; y=t; label=6;};
border i(t=10, dims[0]-10){x=dims[0]-t; y=dims[1]-10; label=7;};
border j(t=10, dims[1]-10){x=10; y=dims[1]-t; label=8;};
int n=3;
mesh fm = buildmesh(a(n) + b(n) + f(n) + d(n) + g(-n) + h(-n) + i(-n) + j(-n));
plot(fm);
fespace Vh(fm, P1);
func real heatsource(){
return 0;
}
Vh u, v;
Vh myHeat = heatsource();
problem heat(u,v) = int2d(fm)(dx(u)*dx(v) + dy(u)dy(v)) - int2d(fm)(myHeatv)
+ on(1, u=0) + on(2, u=0) + on(3, u=0) + on(4, u=0) + on(5, u=0) + on(6, u=0)
+ on(7, u=10) + on(8, u=0);
heat;
plot(u, wait=true, fill=true, dim=2);
ofstream out (“map.txt” );
out.fixed;
for(int i=0; i<dims[0]; i++){
for(int j=0; j<dims[1]; j++){
out << u(i,j) << " ";
}
out << endl;
}
Then freefem display:
While matplotlib display:
It seems like some values are assigned inside the interior square, while not displayed by freefem
Hello
I often use something like (in the case of magnetic field created by a wire)
for (int i=0;i<Th.nt;i++)
{ for (int j=0; j <3; j++)
ff<<Th[i][j].x << " "<< Th[i][j].y<< " "<<Bx[][Vh(i,j)]<< " "<<By[][Vh(i,j)]<<endl;
ff<<Th[i][0].x << " "<< Th[i][0].y<< " “<<Bx[][Vh(i,0)]<< " “<<By[][Vh(i,0)]<<”\n\n\n”;
}
savemesh(Th,“wire.msh”);
HTH
Denis
Thanks for your help. I think my code in FreeFem++ was ok. It was just the way to display it with opencv or matplotlib that was wrong.
Sorry for that
Dear user
just a small remark,
freefem dowhen the point is outside the domain depend of the finite element.
2 case if the finite element is continuous then the we get a continuous prolongement
and the finite is discontinuous then the prolongement is zero.
In all your case you have a hole in the mesh so the finite element is not clearly difined
Remark: you can use the fonction chi(Th)(xx,yy) to know if the point (xx,yy) is in the mesh Th.
