Corner singularity on polar coordinates

(Eliseo Chacon) #1

Dear all,

the first code below shows the standard and well-known singularity on the corner domain.

My surprise is when solving the same problem using polar coordinates and the singularity
on the origin seems to unfold to the straight line r=0 times theta in (0,beta) and is not well visualized.

The second code is the problem solved straight using polar coordinates.

Any idea on this? Thanks!

// Code using x-y coordinates
real beta=4.8;
border a(t=0,1) { x=t; y=0; label=1;}
border b2(t=0,pi/2) { x=cos(t); y=sin(t); label=2;}
border b3(t=pi/2,3*pi/2) { x=cos(t); y=sin(t); label=3;}
border b4(t=3*pi/2,beta) { x=cos(t); y=sin(t); label=4;}
border c(t=cos(beta),0) { x=t; y=sin(beta)/cos(beta)*t; label=1;}

int n=20;
mesh Th = buildmesh(a(n)+b2(n)+b3(n)+b4(n)+c(n));

fespace Vh(Th,P1);
Vh u,v, gradu;
problem ver(u,v)=int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))
for(int i=0;i<3;i++)
plot(u,dim=3,fill=1,value=1,cmm="u--->Mesh adaptation n= "+i,wait=1);
plot(Th,gradu,fill=1,dim=3,value=1,cmm="gradu-->Mesh adaptation n= "+i,wait=1);
Th = adaptmesh(Th,gradu);

// Code using r-theta coordinates
real beta=4;  //primero resolvemos con beta<pi
border a1(t=0,1) {x=t;y=0;label=1;}
border a2(t=0,beta) {x=1;y=t;label=2;}
border a3(t=0,1) {x=1-t;y=beta;label=3;}
border a4(t=0,beta) {x=0;y=beta-t;label=4;}

int n=10;
mesh th0=buildmesh(a1(n)+a2(n)+a3(n)+a4(n));
fespace Vh(th0,P1);
func ff=pow(x,pi/beta)*sin(pi/beta*y);
Vh u,v,gradu,gradutrue,utrue=ff,dxu,dxutrue;

problem tarea(u,v) = int2d(th0)(dx(u)*dx(v)
					+ on(2,u=sin(pi*y/beta))
					+ on(1,u=0)
					+ on(3,u=0)
					+ on(4,u=0);

for(int i=0;i<6;i++)
// SINGULARYTY APPEARS IN  u_r (derivative with respect first argument)

plot(u,dim=3,fill=1,value=1,cmm="u--->Mesh adaptation n= "+i,wait=1);
plot(utrue,dim=3,fill=1,value=1,cmm="utrue--->Mesh adaptation n= "+i,wait=1);
plot(th0,dxu,fill=1,dim=3,value=1,cmm="dxu-->Mesh adaptation n= "+i,wait=1);
plot(th0,dxu,fill=1,dim=3,value=1,cmm="dxutrue-->Mesh adaptation n= "+i,wait=1);
th0 = adaptmesh(th0,dxu);