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));
plot(Th,wait=1);
fespace Vh(Th,P1);
Vh u,v, gradu;
problem ver(u,v)=int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))
+on(2,u=sin(atan(y/(x+1e-10))*pi/beta))
+on(3,u=sin((atan(y/(x-1e-10))+pi)*pi/beta))
+on(4,u=sin((atan(y/(x+1e-10))+2*pi)*pi/beta))
+on(1,u=0)
;
for(int i=0;i<3;i++)
{
ver;
gradu=sqrt(dx(u)*dx(u)+dy(u)*dy(u));
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)
-1/(x+1e-10)*dx(u)*v
+1/(x+1e-10)^2*dy(u)*dy(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++)
{
tarea;
//
// SINGULARYTY APPEARS IN u_r (derivative with respect first argument)
//
//gradu=sqrt(dx(u)*dx(u)+dy(u)*dy(u));
//gradutrue=sqrt(dx(utrue)*dx(utrue)+dy(utrue)*dy(utrue));
//
dxu=dx(u);
dxutrue=dx(utrue);
gradu=abs(dx(u));
gradutrue=abs(dx(utrue));
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);
}