Paraview - incorporate solution computed for a certain mesh in a transformation of the same mesh

General Discussion
1

Dear Community Members,
I would like to ask for your help in merit to a recent problem I run up against.

I am trying to solve a boundary value problem defined over the square [0,pi] x [0,H], where H>0.
To define such a square I write:

mesh Square=square(50,50,[x*pi,y*H]);

The boundary condition is imposed along the edge of the square which is labelled as 1. When I run the code and export the solution u in ParaView, I obtain the following result

1632×1242 60.2 KB

The result is clearly consistent with the boundary condition I have imposed, since the solution is defined over a finite element space given over the mesh Square, and such a solution vanishes along the edge [0,pi] x {0}.

I would like to see the displacement along a cylindrical surface having radius R and height H, obtained by transforming the mesh Square as follows
meshS cylinder=movemesh23(Square, transfo=[R*cos(x), R*sin(x), y]);
and I then expect the solution to vanish along the cylinder arc corresponding to the side of the square having label 1.

However, this is not the case. When I incorporate the solution u in the deformed mesh cylinder via the command
savevtk("myfile.vtu", cylinder, u);
I obtain the following result

1632×1242 64.2 KB

which is not consistent with the boundary conditions I would like to impose.

I cannot see why the result is inconsistent. Could you please help me understand what I did wrong?

Many thanks and best regards,

Paolo

2

Could you please paste a full MWE? Thanks.

3

Hi @prj.
I am not sending you the original code because it would be very time consuming to read it all.
I managed to create a simpler example in the linear case where the same strange effect appears.

I do not know if it is mathematically meaningful, but for sure the final results do not coincide.

load "iovtk"
load "msh3"
real H=2.0;
real R=0.5;
mesh Square=square(50,50,[pi*x,H*y]);
fespace Vh(Square,[P1,P1,P1]);
Vh [f1,f2,f3]=[0.0,0.0,0.1];
Vh [u1,u2,u3], [v1,v2,v3];
solve PoissonVec([u1,u2,u3], [v1,v2,v3])=
int2d(Square)(dx(u1)*dx(v1)+dy(u1)*dy(v1)+dx(u2)*dx(v2)+dy(u2)*dy(v2)+dx(u3)*dx(v3)+dy(u3)*dy(v3))
-int2d(Square)(f1*v1 + f2* v2 +f3*v3)+on(1,u1=0,u2=0,u3=0);
savevtk("Original.vtu",Square,[u1,u2,u3]);
meshS cyl=movemesh23(Square,transfo=[R*cos(x),R*sin(x),y]);
savevtk("Cylinder.vtu",cyl,[u1,u2,u3]);

4

I’m not sure I fully understand why, but I think the following gives you what you are looking for.

load "iovtk"
load "msh3"
real H=2.0;
real R=0.5;
mesh Square=square(50,50,[pi*x-pi/2.0,H*y]);
fespace Vh(Square,[P1,P1,P1]);
Vh [f1,f2,f3]=[0.0,0.0,0.1];
Vh [u1,u2,u3], [v1,v2,v3];
solve PoissonVec([u1,u2,u3], [v1,v2,v3])=
int2d(Square)(dx(u1)*dx(v1)+dy(u1)*dy(v1)+dx(u2)*dx(v2)+dy(u2)*dy(v2)+dx(u3)*dx(v3)+dy(u3)*dy(v3))
-int2d(Square)(f1*v1 + f2* v2 +f3*v3)+on(1,u1=0,u2=0,u3=0);
savevtk("Original.vtu",Square,[u1,u2,u3]);
meshS cyl=movemesh23(Square,transfo=[R*cos(x),y,R*sin(x)]);
savevtk("Cylinder.vtu",cyl,[u1,u2,u3]);
5

Thank you for your help @prj. Your result is certainly correct and consistent, but I noticed that both the surface parametrisation and the square position were changed. Therefore I am wondering:

  1. What did I do wrong? Why does your solution work and mine does not?
  2. Would it be possible to obtain the same correct result using the parametrisation I originally proposed (i.e., [R*cos(x),R*sin(x),y])?

Thank you very much again for your help!

6
  1. I think it’s just a matter of defining what are the (x,y) and (x,y,z) axes.
  2. I’m not sure it will be possible. Your two first functions in your transformation are y-invariant, while your bidimensional function is not.

I’m not a movemesh23 expert, so take this with a grain of salt :slight_smile:

7

Hi, sorry if I insist again but I remember that few months ago I studied a different problem with the same cylindrical surface and the outcome was correct. At that time I was using version 4.4. of FreeFem.
On that occasion, I used the parametrisation [R*cos(x), R*sin(x), y] and I exported the result in the deprecated .vtk format.
Then I merged the mesh data and the solution to visualise the deformation.

I would like to try the same, but I noticed that the newer version of ParaView (version 5.8.0) does not process .vtk format any longer.
Do you know if there are any working scripts converting .vtk files into .vtu files?
I tried the script below, but it does not work

8

I know this script was working in the past. Be careful to use bin = false in savevtk when exporting using .vtk, I think last we checked it was bugged (but it is fine with .vtu).

9

Thanks for the script link but, unfortunately, I cannot run it because my PYTHONPATH variable is not properly configured.
I am quite rusty with Python so it will take me some time before I can fix the problem.

Anyhow, I tried using the old .vtk format with the option bin=false you suggested, and I implemented the same strategy as in my paper.
The result is now consistent with my original parametrisation

Do you think the problem was owing to the .vtu format itself?
Could you also please add few lines talking about the option bin=false in the FreeFem manual? I found no information in the current manual version.

Thanks again for your help. I will indicate your previous post as the solution to the issue I raised.

This is the working version of the example:

//beginning of code
load "iovtk"
load "msh3"
real H=2.0;
real R=0.5;
mesh Square=square(50,50,[pi*x,H*y]);
fespace Vh(Square,[P1,P1,P1]);
Vh [f1,f2,f3]=[0.0,0.0,0.1];
Vh [u1,u2,u3], [v1,v2,v3];
problem PoissonVec([u1,u2,u3], [v1,v2,v3])=
int2d(Square)(dx(u1)*dx(v1)+dy(u1)*dy(v1)+dx(u2)*dx(v2)+dy(u2)*dy(v2)+dx(u3)*dx(v3)+dy(u3)*dy(v3))
-int2d(Square)(f1*v1 + f2* v2 +f3*v3)+on(1,u1=0,u2=0,u3=0);
PoissonVec;
savevtk("SolutionXYZ_timestep0.vtk",Square,[u1,u2,u3],bin=false);
meshS cyl=movemesh23(Square,transfo=[R*cos(x),R*sin(x),y]);
savevtk("Solmesh3D.vtk",cyl,bin=false);
//end of code

10

No, I think the problem comes from .vtk being improperly exported. It’s deprecated for a good reason, you should stick to .vtu. If you add this at the bottom of your script, you’ll see that the result is still not what you expect.

[...]
//end of code
fespace VhS(cyl, P1);
VhS u3S = u3;
plot(u3S);

If you add it at the bottom of my modification to your script, you’ll get the desired mapping.

11

Yes, you are right. I noticed that the cylinder exported in .vtk is deformed in the wrong direction.
However, the square in .vtk and the square in .vtu also differ in the solution intensity.

I need to stick to the parametrisation [R*cos(x), R*sin(x), y] and I also need to define the finite element space over the square; I cannot define it on the deformed reference configuration (cylinder) because I have to write the code in curvilinear coordinates.

Anyhow, I tried your proposed solution in my research code and the problem persists in the deformed reference configuration (cylinder).

12

I think I have found a solution, which basically consists in transforming the mesh Square into a cylinder directly in ParaView. I tried this strategy and it works for nonzero volume meshes too.

This approach, however, makes no use of the movemesh-like functions since all the transformations are directly applied in ParaView now.

I would have very much preferred to avoid using ParaView filters because I am quite rusty with Python, but at least it works.