function [h,ner,ratio,rner,rratio]=testMESDbc(f,s,x,c=0,M=5)
#testing convergence order of linear MES with Dirichlet bc
#-u''+cu=f on [a,b] 
#u(a)=al u(b)=be
#INPUT
#f- rhs function handle i.e. 
#s- function handle of the solution
#x- vector with meshpoints - a=x(1) b=x(length(N)) 
#c- coeff from the eq default =0
#M - no of refinement default M=5
#OUTPUT
#ner - Mx3 matrix 1st column with discrete max norm of I_hs- u_h
#                 2nd             L^2 norm of I_hs- u_h
#                 3rd             H^1 seminorm of I_hs- u_h
#rner - as ner but with real error on larger  mesh
#ratio  M-1 x 3 matrix of ratios ner(j,k)/ner(j+1,k)- err(h)/err(h/2)
#where err(h) given error on mesh x; err(h/2) error on refined mesh i.e.
#      by adding midpoints (x_k+x_{k+1})/2 to t he given mesh
#rratio -same as ratio but for rner

N=length(x)-1;
a=x(1);b=x(N+1);
al=s(a);be=s(b);

for k=1:M,
 [u,t,hh,A,B]=FEM1Dsolver(f,x(1),x(N+1),x,al,be,c);
 h(k)=hh;
 er=s(t)-u;
 der=er(1:N)-er(2:N+1);
 ner(k,1)=norm(er,'inf');#L^\infty norm error 
 ner(k,2)=sqrt(er'*B*er);#L^2 norm error 
 ner(k,3)=sqrt(er'*A*er);#H^1 seminorm error 
# ner(k,4)=norm(der,2)/sqrt(hh);
 [rner(k,1),rerm1,rner(k,2),rner(k,3)]=realerror(f,s,a,b,u,x,max(1000,20*N+1));
 y=0.5*(x(2:N+1)+x(1:N));
 xn=zeros(1,2*N+1);
 xn(1:2:2*N+1)=x;
 xn(2:2:2*N)=y;
 x=xn;N=2*N;
endfor

ratio= ner(1:M-1,:)./ner(2:M,:);
rratio= rner(1:M-1,:)./rner(2:M,:);
endfunction

function [rerm0,rerm1,rer0,rer1]=realerror(f,s,a,b,u,x,M=1000)
#real test
xs=linspace(a,b,M+1);
[pp]=linspinterp(x,u); #values of FEM solution on the mesh xs 
ys=ppleval(xs,pp);
#plot(xs,ys,";FEM;",xs,s(xs),";sol;");
err=ys-s(xs);
h=(b-a)/M;
derr=err(2:M+1)-err(1:M);derr=derr/h; #derivative
rerm0=norm(err,'inf');
rerm1=norm(derr,'inf');
rer0=norm(err,2)*sqrt(h);
rer1=norm(derr,2)*sqrt(h);
endfunction



function [pp]=linspinterp(w,wy)
#INPUT
#w - knots 
#wy - values at vector (at knots of w)
#OUTPUT
#pp - struct with coeff of linear spline interpolating wy at w


[M,N]=size(w);
if(M<N) w=w';;  endif;
[M,N]=size(wy);
if(M<N) wy=wy';  endif;
Nw=length(w)-1;#No of subintervals
hh=w(2:Nw+1)-w(1:Nw); #length of subintervals

if(min(hh)<=0)
 disp("wrong knots - not ordered; aborting");
 return;
endif
a=w(1);b=max(w);

c=zeros(Nw,2);
c(:,1)=wy(1:Nw);
c(:,2)=(wy(2:Nw+1)-wy(1:Nw))./hh;

pp=struct("Nw",Nw,"a",a,"b",b,"w",w,"hh",hh,"c",c);

endfunction



function ys=ppleval(xs,pp)
#computes value of linear spline from struct pp on the mesh xs
#xs must ne ordered i.e. be incresing vector; xs \subset [a,b]
#ys - valus of linear spline (from struct pp) at meshpoints of xs
 b=pp.b;a=pp.a;
 c=pp.c; 
 hh=pp.hh;
 w=pp.w;
 N=pp.Nw;
 M=length(xs);
 ys=zeros(size(xs));
 k=j=1;
 A=c(1,1); #value of pp at w(1)
 B=0;
 
 k=1;
 while(j<=M)
  A=c(k,1);B=c(k,2);
  while(xs(j)<=w(k+1))
    ys(j++)=A+B*(xs(j)-w(k));
    if(j>M) return; end
  endwhile 
  k++;
 endwhile

endfunction


function [u,gr,h,A,B,C,S,F]=FEM1Dsolver(f,a,b,x=[],al=0,be=0,c=0,ad=0,di=1.0)
#FEM solver for -di*u''+ad*u'+ c*u=f with
#Neumann/Robin bnd cond -u'(a)+cau(a)=al u'(b)+cb*u(b)=be
#or Dirichlet bnd : u(a)=al u(b)=be
#Input: f - function handle
#a,b- ends of interval
#x- no of mesh points OR a vector with mesh points : e.g. (xk=a+k*h) - default such with h=(b-a)/20 if x=[] then mesh equdistant
#al, beta - rhs of  bnd conditions
#c - either a value of c (Dirichelt bnd) or a vector [c,ca,cb] - then we have Robin bnd cond
#d - advection parameter 
#output:
#u - FEM1D linear  solution on these mesh - the rhs computed by trapezoidal rule
#gr =- grid
#h - max mesh parameter i.e. h=max h_k h_k=x_{k}-x_{k-1} k=1,..,N
#A,B,C,S - matrices - diffusion, mass, advection - final FEM matrix
#F - rhs of the FEM system 
#Ex:
#[u,t]=FEM1Dsolver(@sin,0,pi);plot(t,u)
#[u,t]=FEM1Dsolver(@sin,a=-1,b=2,[],sin(a),sin(b));plot(t,u) #nonhomogenous D bc
#[u,t]=FEM1Dsolver(@sin,a=-1,b=2,[],0.5*sin(a),0.5*sin(b),1.0);plot(t,u,t,sin(t)/2)
#[u,t]=FEM1Dsolver(@(x)2*sin(x)+cos(x),a=-1,b=2,[],sin(a),sin(b),1.0,1.0);plot(t,u,t,sin(t))
#[u,t]=FEM1Dsolver(@sin,0,pi,[],-cos(0),cos(pi),[0,0,0]);plot(t,u)#Neumanna bnd
#[u,t]=FEM1Dsolver(@sin,a=-2,b=1,[],-cos(a)+sin(a),cos(b)+sin(b),[0,1,1]);plot(t,u)#Robin bnd
#[u,t]=FEM1Dsolver(@(x) 1,0,1,100,0,-1,[0,0,0]);plot(t,u)#neumann  bnd 
#[u,t]=FEM1Dsolver(@(x) 1,0,1,100,0,0,[0,0,0]);# neumann bnd cond -  WRONG BND cond
#[u,t]=FEM1Dsolver(@(x) 100.0,0,2,100,0,1);plot(t,u) #very concave -u''=100
#[u,t]=FEM1Dsolver(@(x) -100.0,0,2,100,0,1);plot(t,u) #very convex -u''=-100
#[u,t]=FEM1Dsolver(@(x) (100+10*x)*(0.5-(x>0)),-1,2,100,-1,1);plot(t,u) #pu
#[u,t,h,A,B,S,F]=FEM1Dsolver(@sin,0,pi);plot(t,u)
#denser grid on ((b-a)/2,b] subinterval
#x1=linspace(0,pi/2,20);x2=linspace(pi/2,pi,1000);x=[x1,x2(2:999)];[u,t]=FEM1Dsolver(@sin,0,pi,x);plot(t,u-sin(t),";err;");
#Advection 
#eps=1e-1;f=@(x) (abs(x)>0.1)./eps;[u,t]=FEM1Dsolver(f,-1,1,100,0,0,0,ad=0);plot(t,u)
#eps=1e-1;f=@(x) (abs(x)>0.1)./eps;[u,t]=FEM1Dsolver(f,-1,1,100,0,0,0,ad=10);plot(t,u)
#eps=1e-1;f=@(x) (abs(x)>0.1)./eps;[u,t]=FEM1Dsolver(f,-1,1,100,0,0,0,ad=-10);plot(t,u)
#very small diffusion with respect to advection
#eps=1e-1;f=@(x) d*(abs(x)>0)./(eps);[u,t]=FEM1Dsolver(f,-1,1,100,0,0,0,-1,1e-3);plot(t,u)#ti does not work. Too small h... with respect to ad/di


#Code starts
type=0; #0 Dirichlet; 1 Robin; 2 neumann


if(isempty(x))
 N=round((b-a)*20);
 x=linspace(a,b,N+1);
else 
   if(length(x)==1)
     N=x; 
     x=linspace(a,b,1+N);
   else
     N=length(x)-1;
     a=x(1);
     b=x(N+1);
   endif  
endif
#N=length(x)-1; 
hh=zeros(1,N); #hh(k)=x(k+1)-x(k)
hh=x(2:N+1)-x(1:N);
h=max(hh);

mdA=1./hh(1:N-1)+1./hh(2:N);
mdA=[1/hh(1),mdA,1/hh(N)]; 
mdB=hh(1:N-1)+hh(2:N); 
mdB=[hh(1),mdB,hh(N)];
subd=supd=hh;
#stiffness
A=sparse(diag(-1.0./subd,-1)+ diag(mdA) + diag(-1.0./supd,1));
#mass
B=sparse(diag(subd./6,-1)+ diag(mdB/3) + diag(supd/6,1));
#\int u'v advection local (\phi_k,\phi_{k+1})=0.5 [-1,]
supd=ones(1,N);subd=-supd;
C=sparse(diag(0.5*subd,-1)+diag(0.5*supd,1));
C(1,1)=-0.5;C(N+1,N+1)=0.5;

F=zeros(N+1,1);
for k=0:N,
 F(k+1)=f(x(k+1));
endfor
#scaling
F=0.5*F.*mdB';
S=di*A+ad*C+c(1)*B;

if(length(c)==1)
 #Dirichlet bnd 
 F(2)-=S(2,1)*al; #-al*\int_a^{x_1} \phi_0'\phi_1'=al/h_0
 F(N)-=S(N,N+1)*be;
 S=S(2:N,2:N);
 F=F(2:N);
else 
 type=1;
#Robin bnd weak formul: 
#\int u'v' + cuvdt + u'(a)v(a)-u'(b)v(b)=\inf fv 
#\int..dt + cau(a)v(a)+cbu(b)v(b)=\intfv +alv(a)+ be v(b)
 F(1)+=al;F(N+1)+=be; #adding contribution of  alv(a)+ be v(b)
 S(1,1)+=c(2);
 S(N+1,N+1)+=c(3); # ca*u(a)v(a)+cb*u(b)v(b) to stiffness matrix
 if(norm([ad,c(2:3)],'inf')<1e-19) 
   type=2;
  ind=1:N;
   #case of pure Neumann condition we assume that u(N+1)=u(b)=0 to get a unique solution
   S=S(ind,ind);
   if(sum(F)>h)
     printf("Pure Neumann cond sum(F(k)<>0! Problem has no solution!\n");
   endif
   F=F(ind);
 endif
endif

u=S\F;
if(type==0)
 u=[al; u; be];
 else 
 if(type==2)
  #u coefficients of the solution in the basis (\phi_0,..\phi_{N-1}, 1)
  #here 
  u=[u;0];
 endif
endif
gr=x'; #returns grid 

endfunction

