clear *
r=1.0;
R=2.5;
global h=0.01;
H=0.5;

N=201;

area=0.0001;

xy=[0,0,1; R,0,2; R,H,2; r,H,1; r,h,1; 0,h,2];

CreateMeshTriangle("cap",xy,area)
[nodes,elem,edges]=ReadMeshTriangle("./cap.1");
# ShowMesh(nodes,elem)
# disp("hit return");pause();

clear volt

function res = volt(xy)
global h;
[n,m]=size(xy);
res=zeros(n,1);
for k=1:n
  if (xy(k,2)<h/2) res(k)=0;
  else res(k)=1;endif
endfor
endfunction

[A,b,n2d]=FEMEquation(nodes,elem,edges,nodes(:,1),0,0,'volt',0);
u=FEMSolveSym(nodes,A,b,n2d,'volt');
# ShowSolutionMTV(nodes,elem,u,"u.mtv");
# system("plotmtv -nodate u.mtv -o volt.eps")
ShowLevelCurves(nodes,elem,u,linspace(0,1,11));
#print("volt.eps","-depsc");
disp("hit return");pause();

x=linspace(0,R,N); 
y=0.0*ones(1,N);
[uVal,grad]=FEMValue([x;y]',nodes,elem,u);

plot(x,grad(:,2))

# trapezoidal integration
flux = 2*(x*grad(:,2)-x(1)*grad(1,2)/2-x(N)*grad(N,2)/2)*R/(N-1)*(h/r**2)
