# 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

# Exercices/Contributions

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