sandbox/easystab/integration_2D.m

    This code is a demonstration of the spatial integration in 2D and with a mapped mesh. To learn more about mapped meshes, please see README#differentiation-with-a-non-rectangular-mesh. We test the integration in 1D in integration.m.

    Dependency

    clear all; clf
format compact 
format long

% parameters
Lx=2*pi;
Ly=2*pi;
Nx=20;
Ny=20;

% differentiation 
[d.x,d.xx,d.wx,x]=dif1D('cheb',0,Lx,Nx);
[d.y,d.yy,d.wy,y]=dif1D('fd',0,Ly,Ny,5);
[D,l]=dif2D(d);

    Integration in 1D

    The variable d is a structure. The integration weights in x are in d.wx, and the integration weights in y are in d.wy. The grid vectors are x and y, and are column vectors. The integration weights are row-vectors such that the matrix multiplication with column vectors is the scalar product: the product of each elements and the sum of them \displaystyle \left(\begin{array}{cccc} w_1 & w_2 & \dots & w_N \end{array}\right) \left(\begin{array}{c} f_1 \\ f_2 \\ \vdots \\ f_N \end{array}\right) = w_1f_1+w_2f_2+\dots+w_Nf_N

    Stretching in y

    Here we test the 2D integration when there is a stretching of the mesh. We do two cases of simple stetching, the first one we stretch in y and the second we stretch in x. For the stretching, we simply scale down the grid by building the vector eta as long as the x grid.

    % stretching in y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' ');disp('stretching in y');
[X,Y]=meshgrid(x,y);
eta=0*x+1+0.5*x.*(x-Lx)*4/Lx^2; 
Y=Y.*repmat(eta',Ny,1);  
Dm=map2D(X,Y,D);

    We chose to integrate something very simple: a constant function f of value 1, such that the surface integral of f is equal to the total area of the mesh.

    % the function to integrate
f=0*X+1;
subplot(2,2,1);
mesh(X,Y,f); view(2); axis tight
xlabel('x');ylabel('y');title('mapped domain')

    Then to perform the integration is simple, we multiply element by element the integration weights in D.w by the elements of f and we sum. The theoretical surface is he surface of the square L_xL_y minus the surface of he parabola defined by eta.

    % surface integral
itheo=2*Lx^2/3
inum=sum(sum(Dm.w.*(0*X+1)))
err=abs(itheo-inum)

    For the integration in y, we multiply element by element the integration weights D.wx with the function f, then we sum on the direction on which we integrate. Since here we integrate along x, the sum is done along the dimension 2 of the arrays. This value should be equal of Ly*eta.

    % integration in y
disp('Integration in y');
iy=sum(Dm.wy.*f,1);
subplot(2,2,2);
plot(x,iy,'b',x,Ly*eta,'r.')
xlabel('x');ylabel('integral in y');
axis([0,Lx,0,Ly]); grid on

    Stretching in x

    We chose something simpler for the stertching in x, we just displace the mesh to the right with the same function eta. The total surface is thus unchanged from L_xL_y.

    % stretching in x
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' ');disp('stretching in x');
[X,Y]=meshgrid(x,y);
eta=1+y.*(y-Ly)*4/Ly^2; 
X=X-repmat(eta,1,Nx);  
Dm=map2D(X,Y,D);

% the function to integrate
f=0*X+1;
subplot(2,2,3);
mesh(X,Y,f); view(2); axis tight
xlabel('x');ylabel('y');title('mapped domain')

% surface integral
itheo=d.wy*(eta*Lx)
inum=sum(sum(Dm.w.*(0*X+1)))
err=itheo-inum

% integration in x
ix=sum(Dm.wx.*f,2);
subplot(2,2,4);
plot(y,ix,'b',y,Lx,'r.')
xlabel('y');ylabel('integral in x');
xlim([0,Ly]); grid on

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

    The figure

    The shape of the domains and the integrals

    The shape of the domains and the integrals

    Exercices/Contributions

    • Please test other types of stretchings
    • Please integrate without stretching
    • Please integrate a non-constant function