Movemesh23 and move solution for visualisation

General Discussion
1

For a parameterised sphere,

if we solve a problem in the parameter space, for example (0, 2pi) x (0, pi)

we can use movemesh23() to generate the mesh of a sphere quickly,

my question is: how can we move the solution to the sphere for visualisation – ideally a node to node transfer such as us[] = u[] ?

To clarify my question: we have a flat 2D mesh Th

fespace Vh(Th, P1, periodic=[[2, y], [4, y]]);

Vh u.
...

meshS ThS = movemesh23(Th, transfo=[a*sin(y)cos(x),asin(y)sin(x),acos(y)],removeduplicate=1);

fespace VhS(ThS, P1);
VhS us;

How can we move u to us?

2

There is a method that works here

but it is related to the property that the mesh nodes are purely transported by th transformation.
In your case you have “removeduplicate=1” that breaks this property (the two meshes do not have the same number of nodes).
Maybe you can do it if you set “removeduplicate=0”

3

Thanks for your help. I wrote a piece of code to find the corresponding nodes, although it looks ugly:

Vh theta=x, phi=y; 

VhS us, xs=x, ys=y, zs=z;

real xx, yy, zz, dd;
for(int i=0; i<xs[].n; i++){
for(int j=0; j<theta[].n; j++){
xx = sin(phi[][j])cos(theta[][j]);
yy = sin(phi[][j])sin(theta[][j]);
zz = cos(phi[][j]);
dd = sqrt( (xs[][i]-xx)(xs[][i]-xx) + (ys[][i]-yy)(ys[][i]-yy) + (zs[][i]-zz)(zs[][i]-zz) );
if (dd < 1.e-5) us[][i]=u[][j];
}
}
4

You don’t need to find the nodes, they follow the map as long as you keep all the nodes (and do not remove some of them).
Thus you can do in several steps:

  1. define your square mesh and your initial periodic function u.
    Transport the mesh to a meshS by movemesh23.

  2. define a a slighly smaller square mesh, and interpolate your u on this mesh to ucut.
    Transpor the mesh to a meshS S1 which is not closed (so that you keep all the nodes).
    Transport ucut to us1

  3. finally interpolate us1 to us.

real a=1.;

mesh Th=square(10,10,[2.*pi*x,pi*y]);
cout << "Th.nt " << Th.nt << " Th.nv " << Th.nv << endl;;
fespace Vh(Th, P1, periodic=[[2, y], [4, y]]);
Vh u;
u=cos(x)*sin(y);

meshS ThS=movemesh23(Th, transfo=[a*sin(y)*cos(x),a*sin(y)*sin(x),a*cos(y)],removeduplicate=1);
cout << "ThS.nt " << ThS.nt << " ThS.nv " << ThS.nv << endl;
fespace VhS(ThS, P1);
VhS us;

mesh Thcut=square(10,10,[2.*pi*0.999*x,0.001+pi*0.998*y]);
cout << "Thcut.nt " << Thcut.nt << " Thcut.nv " << Thcut.nv << endl;;
fespace Vhcut(Thcut, P1);
Vhcut ucut;
ucut=u;

meshS ThS1=movemesh23(Thcut, transfo=[a*sin(y)*cos(x),a*sin(y)*sin(x),a*cos(y)],removeduplicate=0);
cout << "ThS1.nt " << ThS1.nt << " ThS1.nv " << ThS1.nv << endl;;
fespace VhS1(ThS1,P1);
VhS1 us1;
us1[]=ucut[];

us=us1;

plot(ThS);

However there are some warning messages “no manifold obj” for ThS, and “unfreed pointers, memory leak” that are a bit scary

1 Like