n=4; coor=[1 0;2 0;3 2;0 2];
coor=[coor;[coor(1,:)]];

N1=3; xi1=[-sqrt(3/5),0,sqrt(3/5)];A1=[5/9,8/9,5/9];
N2=20; xi2=-1+1/20:2/20:1-1/20;A2=2/20*ones(1,N2);
##N2=50; xi2=-1+1/50:2/50:1-1/50;A2=2/50*ones(1,N2);
##N2=150; xi2=-1+1/150:2/150:1-1/150;A2=2/150*ones(1,N2);

pint1=[(1-xi1)/2;(1+xi1)/2];
pint2=[(1-xi2)/2;(1+xi2)/2];

d=[coor(2:n+1,:)-coor(1:n,:)];ll=norm(d,"rows");jac=ll/2;

normal=[d(:,2)./ll,-d(:,1)./ll]

N=N1; xi=xi1; A=A1; pint=pint1;
perim=0;
for i=1:n
 perim=perim+A*ones(N,1)*jac(i);
endfor 

perim

res=0;
for i=1:n
 for k=1:N
  x=pint(1,k)*coor(i,:)'+pint(2,k)*coor(i+1,:)';
  res=res+A(k)*1/2*normal(i,1:2)*x*jac(i);
 endfor
endfor 

volume=res

res=[0;0];
for i=1:n
 for k=1:N
  res=res+A(k)*normal(i,:)'*jac(i);
 endfor
endfor 

intnormal=res

res=[0;0];
for i=1:n
 for k=1:N
  x=pint(1,k)*coor(i,:)'+pint(2,k)*coor(i+1,:)';
  xpen=normal(i,:)*x;
  res=res+A(k)*xpen*x*jac(i);
 endfor
endfor 

cg=res/(3*volume)

dt=0.05;
empux=[]; tempo=[];
for t=-1:dt:3

Z=@(t,x1) t;
##Z=@(t,x1) 1+exp(-(x1-0.1*t+3)^2/0.04);
##Z=@(t,x1) 1+0.3*sin(kk*(x1-0.1*t));
p=@(t,x) max(0,Z(t,x(1))-x(2));
##p=@(t,x) max(0,Z(t,x(1))-x(2));
res=[0;0];
for i=1:n
 for k=1:N
  x=pint(1,k)*coor(i,:)'+pint(2,k)*coor(i+1,:)';
  res=res+A(k)*normal(i,1:2)'*(-p(t,x))*jac(i);
 endfor
endfor 
empux=[empux,res];
tempo=[tempo,t];

endfor
empux1=empux;

##################################
dt=0.5;
empux=[];tempo2=[];
for t=0:dt:70

##Z=@(t,x1) t;
Z=@(t,x1) 1+exp(-(x1-0.1*t+3)^2/0.01);
##Z=@(t,x1) 1+0.3*sin(kk*(x1-0.1*t));
pz=@(t,x) Z(t,x(1))-x(2);
##p=@(t,x) max(0,Z(t,x(1))-x(2));
res=[0;0];
for i=1:n
 N=N1; xi=xi1; A=A1; pint=pint1;
 if (sign(pz(t,coor(i,:)'))!=sign(pz(t,coor(i+1,:)'))) 
  N=N2; xi=xi2; A=A2; pint=pint2;
 endif
 for k=1:N
  x=pint(1,k)*coor(i,:)'+pint(2,k)*coor(i+1,:)';
  res=res+A(k)*normal(i,1:2)'*(-max(0,pz(t,x)))*jac(i);
 endfor
endfor 
empux=[empux,res];
tempo2=[tempo2,t];

endfor
empux2=empux;

##################################
dt=0.5;
tor=[];tempo3=[];
for t=0:dt:70

##Z=@(t,x1) t;
Z=@(t,x1) 1.+exp(-(x1-0.1*t+3)^2/0.01);
##Z=@(t,x1) 1+0.3*sin(kk*(x1-0.1*t));
pz=@(t,x) Z(t,x(1))-x(2);
##p=@(t,x) max(0,Z(t,x(1))-x(2));
res=0;
for i=1:n
 N=N1; xi=xi1; A=A1; pint=pint1;
 if (sign(pz(t,coor(i,:)'))!=sign(pz(t,coor(i+1,:)'))) 
  N=N1; xi=xi1; A=A1; pint=pint1;
 endif
 for k=1:N
  x=pint(1,k)*coor(i,:)'+pint(2,k)*coor(i+1,:)';
  ff=(x(1)-cg(1))*normal(i,2)-(x(2)-cg(2))*normal(i,1);
  res=res+A(k)*ff*(-max(0,pz(t,x)))*jac(i);
 endfor
endfor 
tor=[tor,res];
tempo3=[tempo3,t];

endfor
tor1=tor;

##################################
dt=0.5;
tor=[];tempo3=[];
for t=0:dt:70

##Z=@(t,x1) t;
Z=@(t,x1) 1.+exp(-(x1-0.1*t+3)^2/0.01);
##Z=@(t,x1) 1+0.3*sin(kk*(x1-0.1*t));
pz=@(t,x) Z(t,x(1))-x(2);
##p=@(t,x) max(0,Z(t,x(1))-x(2));
res=0;
for i=1:n
 N=N2; xi=xi2; A=A2; pint=pint2;
 if (sign(pz(t,coor(i,:)'))!=sign(pz(t,coor(i+1,:)'))) 
  N=N2; xi=xi2; A=A2; pint=pint2;
 endif
 for k=1:N
  x=pint(1,k)*coor(i,:)'+pint(2,k)*coor(i+1,:)';
  ff=(x(1)-cg(1))*normal(i,2)-(x(2)-cg(2))*normal(i,1);
  res=res+A(k)*ff*(-max(0,pz(t,x)))*jac(i);
 endfor
endfor 
tor=[tor,res];
tempo3=[tempo3,t];

endfor
tor2=tor;

##################################
dt=0.5;
tor=[];tempo3=[];
for t=0:dt:70

##Z=@(t,x1) t;
Z=@(t,x1) 1.+exp(-(x1-0.1*t+3)^2/0.01);
##Z=@(t,x1) 1+0.3*sin(kk*(x1-0.1*t));
pz=@(t,x) Z(t,x(1))-x(2);
##p=@(t,x) max(0,Z(t,x(1))-x(2));
res=0;
for i=1:n
 N=N1; xi=xi1; A=A1; pint=pint1;
 if (sign(pz(t,coor(i,:)'))!=sign(pz(t,coor(i+1,:)'))) 
  N=N2; xi=xi2; A=A2; pint=pint2;
 endif
 for k=1:N
  x=pint(1,k)*coor(i,:)'+pint(2,k)*coor(i+1,:)';
  ff=(x(1)-cg(1))*normal(i,2)-(x(2)-cg(2))*normal(i,1);
  res=res+A(k)*ff*(-max(0,pz(t,x)))*jac(i);
 endfor
endfor 
tor=[tor,res];
tempo3=[tempo3,t];

endfor
tor3=tor;

##################################