sandbox/easystab/multi_domain_2D_tri.m

    Just like for multi_domain_1D.m and multi_domain_2D.m but with an interesting domain (that we could not do otherwise than using multi-domains). Please check these other codes for the detailed explanations.

    clear all; clf; 

% parameters 
nx=20; % number of grid nodes in x
ny=20; % number of grid nodes in y
method='cheb'; % which discretization to use
pts=9;  % how mnay points in the finite difference stencil

% differentiation 
[d.x,d.xx,d.wx,x]=dif1D(method,0,1,nx,pts);
[d.y,d.yy,d.wy,y]=dif1D(method,0,1,ny,pts);
[D1,l1,X1,Y1,Z1,I1,NN1]=dif2D(d,x,y);

[d.x,d.xx,d.wx,x]=dif1D(method,1,2,nx,pts);
[d.y,d.yy,d.wy,y]=dif1D(method,0,1,ny,pts);
[D2,l2,X2,Y2,Z2,I2,NN2]=dif2D(d,x,y);

[d.x,d.xx,d.wx,x]=dif1D(method,0,1,nx,pts);
[d.y,d.yy,d.wy,y]=dif1D(method,1,2,ny,pts);
[D3,l3,X3,Y3,Z3,I3,NN3]=dif2D(d,x,y);

D.x=blkdiag(D1.x,D2.x,D3.x);
D.y=blkdiag(D1.y,D2.y,D3.y);
D.xx=blkdiag(D1.xx,D2.xx,D3.xx);
D.yy=blkdiag(D1.yy,D2.yy,D3.yy);
I=blkdiag(I1,I2,I3);
Z=blkdiag(Z1,Z2,Z3);

NN=NN1+NN2+NN3;
f1=(1:NN1)'; 
f2=(f1(end)+(1:NN2))';
f3=(f2(end)+(1:NN3))';

X=[X1(:);X2(:);X3(:)];
Y=[Y1(:);Y2(:);Y3(:)];

    We solve simply \displaystyle f_{xx}+f_{yy}=0 

    % System matrices
A=D.xx+D.yy;
b=0*X;

    For validation, we use an analytical solution of the Poisson problem with Dirichlet boundary conditions. \displaystyle f(x,y)=3x^2y-y^3 

    % the analytical solution
ftheo=3*X.^2.*Y-Y.^3;

    To get the connections, it is good to draw a sketch of how the domains are located with respect to each others.

    % Boundary conditions
dir=[f1([l1.left; l1.bot]); ...
    f2([l2.right; l2.top; l2.bot]); ...
    f3([l3.left; l3.top; l3.right])];
loc=[dir; f1(l1.right);f2(l2.left);f1(l1.top);f3(l3.bot)];
C=[I(dir,:); ...
    I(f1(l1.right),:)-I(f2(l2.left),:); ...
    D.x(f1(l1.right),:)-D.x(f2(l2.left),:); ...
    I(f1(l1.top),:)-I(f3(l3.bot),:); ...
    D.y(f1(l1.top),:)-D.y(f3(l3.bot),:)];

A(loc,:)=C;
b(loc)=C*ftheo;

% solve system
f=A\b;

% validation
err=norm(f-ftheo,2)

mesh(X1,Y1,reshape(f(f1),ny,nx),'edgecolor','b');
hold on
mesh(X2,Y2,reshape(f(f2),ny,nx),'edgecolor','r');
mesh(X3,Y3,reshape(f(f3),ny,nx),'edgecolor','m');
view(72,50)

xlabel('x');ylabel('y');
legend('domain 1','domain 2','domain 3');

set(gcf,'paperpositionmode','auto');
print('-dpng','-r80','multi_domain_2D_tri2.png');

    Here is the figure produced by the code as it is:

    And here is an other test with \displaystyle f_{xx}+f{yy}=-1 and homogeneous boundary conditions.

    %}