# Solving PDE on sphere using example provided by FreeFem

**URL:** <https://community.freefem.org/t/solving-pde-on-sphere-using-example-provided-by-freefem/1049>\
**Category:** General Discussion\
**Created:** [June 18, 2021, 1:16pm UTC](https://community.freefem.org/t/solving-pde-on-sphere-using-example-provided-by-freefem/1049 "2021-06-18T13:16:17Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![ejkc](https://avatars.discourse-cdn.com/v4/letter/e/ecccb3/32.png) [@ejkc](https://community.freefem.org/u/ejkc)\
**Post date:** [June 18, 2021, 1:16pm UTC](https://community.freefem.org/t/solving-pde-on-sphere-using-example-provided-by-freefem/1049/1 "2021-06-18T13:16:17Z")

</div>

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-2_sqrt(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 = cosdiedre_cosico+sindiedre_sinico;//angle diedre -pi/2 - ico  
real sinicod = sindiedre_cosico-cosdiedre_sinico;//angle diedre -pi/2 - ico

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

int n = max(RZ,max(RX,RY))_2_pi/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.5_sqrt(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,sinico_y+cosico_z,-cosico_y+sinico_z],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,-cosdiedre_y+sindiedre_z,-sindiedre_y-cosdiedre_z],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-sinpi3_cosico,z],label=ref,orientation=1);//translation pour que la figure soit sur le bord du pentagone  
meshS T2 = movemeshS(T1,transfo=[cos2pi5_x-sin2pi5_y,sin2pi5_x+cos2pi5_y,z],label=ref,orientation=1);  
meshS T3 = movemeshS(T2,transfo=[cos2pi5_x-sin2pi5_y,sin2pi5_x+cos2pi5_y,z],label=ref);  
meshS T4 = movemeshS(T3,transfo=[cos2pi5_x-sin2pi5_y,sin2pi5_x+cos2pi5_y,z],label=ref);  
meshS T5 = movemeshS(T4,transfo=[cos2pi5_x-sin2pi5_y,sin2pi5_x+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.5_h],label=ref);//moitié supérieure  
meshS Tdemi1 = movemeshS(Tdemi0,transfo=[x,y,-z],label=ref);//moitié inférieure  
meshS Tdemi1rot = movemeshS(Tdemi1,transfo=[cospi5_x-sinpi5_y,sinpi5_x+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(aa_x_x+bb_y_y+cc_z_z);  
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.0_v\*u);  
Problem;  
plot(Th3sph,u,wait=1);

Thank you in advance!