Solving PDE on sphere using example provided by FreeFem

General Discussion
1

Hi!

I am new to freefem and a bit confused. I try to solve a pde on the sphere. In order to generate the mesh, I use one of the examples that comes with FreeFem. I then write down my weak form as in the tutorials, but it does not work!
I solve “something” but the solution does not plot. It should be one of the spherical harmonics.

Also, if I have solved two PDEs on the same mesh, can I somehow sum the solutions?

See code below.

load “gsl”
load “msh3”
load “medit”
load “tetgen”
load “mmg”
load “iovtk”

func meshS Ellipsoide20 (real RX,real RY, real RZ,real hsize,int L,int orientation)
{
//===================================================================================
//Angles utiles
//===================================================================================
real tan3pi10 = sqrt(25.+10.sqrt(5.))/5.;//3pi/10 angle entre deux aretes du pentagone
real sin3pi10 = (sqrt(5)+1)/4;//3pi/10 angle entre deux aretes du pentagone
real cos3pi10 = sqrt(10-2sqrt(5))/4;//3pi/10 angle entre deux aretes du pentagone

real cosdiedre = sqrt(5)/3; //angle diedre de l’icosaedre -pi/2
real sindiedre = 2./3; //angle diedre de l’icosaedre -pi/2

real cosico = tan3pi10/sqrt(3); //angle entre une face de la pyramide pentagonale par rapport à l’horizontale
real sinico = sqrt(1-square(cosico)); //angle entre une face de la pyramide pentagonale par rapport à l’horizontale

real sin2pi5 = sqrt(10+2*sqrt(5))/4; //2pi/5 angle pour la rotation des aretes du pentagone
real cos2pi5 = (sqrt(5)-1)/4; //2pi/5 angle pour la rotation des aretes du pentagone

real cosicod = cosdiedrecosico+sindiedresinico;//angle diedre -pi/2 - ico
real sinicod = sindiedrecosico-cosdiedresinico;//angle diedre -pi/2 - ico

real sinpi3 = sqrt(3)/2; //angle du triangle equilateral

int n = max(RZ,max(RX,RY))2pi/6/hsize;

real sinpi5 = cos3pi10;//pi/5 angle de décalage entre deux demi icosaedre
real cospi5 = sin3pi10;//pi/5 angle de décalage entre deux demi icosaedre

real tanpi10 = sqrt(25.-10.sqrt(5.))/5.;//pi/10
real h = 0.5sqrt(3-square(tanpi10));//hauteur du prisme d’ordre 5;

//=================================================================================
//Construction du triangle equilateral en 2D
//=================================================================================
border a(t=0,1){x=t; y=0; label =1;};
border b(t=1,0.5){x=t; y=sqrt(3)(1-t); label =2;};
border c(t=0.5,0){x=t; y=sqrt(3)(t); label =3;};
mesh Triangle= buildmesh(a(n)+b(n)+c(n)); //traingle equilateral
func f = 1;
int[int] ref=[0,1];

meshS Triangle3 = movemesh23(Triangle,transfo=[x,0,y],label=ref);//trianglesup
meshS TriangleS = movemeshS(Triangle3,transfo=[x,sinicoy+cosicoz,-cosicoy+sinicoz],label=ref);//rotation de -(pi - diedre) par rapport à l’axe des x pour former une face de la pyramide pentagonale
meshS TriangleI = movemeshS(TriangleS,transfo=[x,-cosdiedrey+sindiedrez,-sindiedrey-cosdiedrez],orientation=-1,label=ref);//triangle inf rotation de l’angle diedre par rapport au triangle sup
meshS Triangles = TriangleI+TriangleS;
//medit(“face pyramide pentagonale + face antiprisme d’ordre 5”,Triangles);

meshS T1 = movemeshS(Triangles,transfo=[x-0.5,y-sinpi3cosico,z],label=ref,orientation=1);//translation pour que la figure soit sur le bord du pentagone
meshS T2 = movemeshS(T1,transfo=[cos2pi5x-sin2pi5y,sin2pi5x+cos2pi5y,z],label=ref,orientation=1);
meshS T3 = movemeshS(T2,transfo=[cos2pi5x-sin2pi5y,sin2pi5x+cos2pi5y,z],label=ref);
meshS T4 = movemeshS(T3,transfo=[cos2pi5x-sin2pi5y,sin2pi5x+cos2pi5y,z],label=ref);
meshS T5 = movemeshS(T4,transfo=[cos2pi5x-sin2pi5y,sin2pi5x+cos2pi5*y,z],label=ref);

meshS Tdemi= T1+T2+T3+T4+T5;//moitié de l’icosaedre

meshS Tdemi0 = movemeshS(Tdemi,transfo=[x,y,z+0.5h],label=ref);//moitié supérieure
meshS Tdemi1 = movemeshS(Tdemi0,transfo=[x,y,-z],label=ref);//moitié inférieure
meshS Tdemi1rot = movemeshS(Tdemi1,transfo=[cospi5x-sinpi5y,sinpi5x+cospi5*y,z],label=ref,orientation=-1);//rotation de la moitié inférieure pour les emboiter
meshS Ticosaedre = Tdemi0+Tdemi1rot;

//=================================================================================
//Construction de la sphere 3D
//=================================================================================
real aa =RX, bb=RY, cc =RZ;
func metric = sqrt(aaxx+bbyy+cczz);
meshS Th = movemeshS(Ticosaedre,transfo=[x/metric,y/metric,z/metric]);
meshS Thb=mmgs(Th,hmin=hsize,hmax=hsize,hgrad=1.1);
//mmgs(ThS,hausd=0.001);
return Thb;
}

//medit(“spheresurface”,Th);
//savemesh(Th,“T.mesh”);
//exec(“ffmedit T.mesh”);
real hsize=0.5;
meshS Thb=Ellipsoide20(1,1,1,hsize,1,1);
Thb=mmgs(Thb,hmin=hsize,hmax=hsize,hgrad=1.1);

//medit(“spheresurface”,Th);
//savemesh(Th,“T.mesh”);
//exec(“ffmedit T.mesh”);
real[int] domain = [0.,0.,0.,1,hsize^3/6];
mesh3 Th3sph=tetg(Thb,switch=“paAAQYY”,nbofregions=1,regionlist=domain);

plot(cmm=“sphere”,Th3sph,wait=1);

fespace Vh(Th3sph, P1);
Vh u, v,w;// Define u and v as piecewise-P1 continuous functions

macro Grad(u) [dx(u),dy(u),dz(u)] // EOM

solve Problem(u,v,solver=CG)
=int2d(Th3sph)(Grad(u)'Grad(v))-int2d(Th3sph)(1.0v*u);
Problem;
plot(Th3sph,u,wait=1);

Thank you in advance!