% beam1.m: probleme  (au'')'' + (bu')' +cu = f      0 <= x  <= 1
% d'apres  beam.m de      femlab-1d.tar.Z         a(x)>0  (ou a=0 et b<0)
%  with boundary conditions
%
%          u = u0,  u' = u1   (clamped)
%                                  seul le cas encastre est au point, hum
%  or
%          u'' = 0,  (au'')' = 0  (free)
%
%  at  x = 0  and  x = 1,  respectively.
%
%  using  C1  piecewice qubic elements and mesh size  h.
%
%
%echo off
clear;
fnom=input('nom du fichier donnees (entre apostrophes)  ');
%  formules de a, b,c et f  dans un fichier
% exemple:   beamd1.d
%1
%2.10/(4*14)
%14*(x^3-x)
%14*(3*x^2-1)
%c
%-1
%0
%c
%1
%0
fid=fopen(fnom);
acopy=fscanf(fid,'%s',1);bcopy=fscanf(fid,'%s',1);
ccopy=fscanf(fid,'%s',1);fcopy=fscanf(fid,'%s',1);
bctype=[1 1];
bounda=fscanf(fid,'%s',1);      % c=clamped , f=free
if bounda=='c'; bctype(1)=0; u0(1)=eval(fscanf(fid,'%s',1));
                             u1(1)=eval(fscanf(fid,'%s',1));
elseif bounda=='f'; bctype(1)=1;
end
bounda=fscanf(fid,'%s',1);
if bounda=='c'; bctype(2)=0; u0(2)=eval(fscanf(fid,'%s',1));
                             u1(2)=eval(fscanf(fid,'%s',1));
elseif bounda=='f'; bctype(2)=1;
end
refine=1;       %Boolean used for refinement loop
first=1;        %Boolean used to indicate if turn in first loop
while refine
 %Clear old data structures (otherwise, selecting next h > previous h
 % will cause an error)
 clear h hh ni nn ndf rhs f fi a b c u ucomp Ducomp DDucomp  xplot;
 disp('Give desired mesh size h')
 h=input('h = ');
 ni=round(1/h);  h = 1/ni; %Modify h according to ni
 disp(['Number of sub-intervals = ',int2str(ni)])
 nn=(ni+1);
 disp(['Number of nodes = ',int2str(nn)])
 ndf=2*nn;
 disp(['Number of degrees of freedom = ',int2str(ndf-2*sum(bctype))])
 %
 %disp('Before you use this progam, enter formulas for the four basis functions and their second derivatives on the unit interval [0,1]');
     b1copy='1-3*x^2+2*x^3';
     b2copy='x-2*x^2+x^3';
     b3copy='1-3*(1-x)^2+2*(1-x)^3';
     b4copy='-(1-x)+2*(1-x)^2-(1-x)^3';
     Db1copy='-6*x+6*x^2';
     Db2copy='1-4*x+3*x^2';
     Db3copy='6*(1-x)-6*(1-x)^2';
     Db4copy='1-4*(1-x)+3*(1-x)^2';
     DDb1copy='-6+12*x';
     DDb2copy='-4+6*x';
     DDb3copy='-6+12*(1-x)';
     DDb4copy='4-6*(1-x)';
 %Set up six-diagonal stiffnes matrix
 %
 gauss=sqrt(3/5);
 g(1)=1/2-gauss/2;
 g(2)=1/2;
 g(3)=1/2+gauss/2;
 w(1)=5/18;
 w(2)=8/18;
 w(3)=5/18;
 %
 %
  diag=zeros(size(1:ndf));
  subdiag=zeros(size(1:ndf-1));
  sub2diag=zeros(size(1:ndf-2));
  sub3diag=zeros(size(1:ndf-3));
  rhs=zeros(size(1:ndf));
  integJ=0;
   %
   %
   %
  for gp=1:3;
  %
  x=g(gp);
  %
  b1=eval(b1copy);
  b2=eval(b2copy);
  b3=eval(b3copy);
  b4=eval(b4copy);
  %
  Db1=eval(Db1copy)/h;
  Db2=eval(Db2copy)/h;
  Db3=eval(Db3copy)/h;
  Db4=eval(Db4copy)/h;
  %
  DDb1=eval(DDb1copy)/h^2;
  DDb2=eval(DDb2copy)/h^2;
  DDb3=eval(DDb3copy)/h^2;
  DDb4=eval(DDb4copy)/h^2;
  %
  for i=1:ni;
   x=(i-1)*h+g(gp)*h;
   a(i)=eval(acopy);
   b(i)=eval(bcopy);
   c(i)=eval(ccopy);
   f(i)=eval(fcopy);
  end;
  %
  hh=h*w(gp);
  %
  diag(3:2:ndf-1)=diag(3:2:ndf-1)+(a.*DDb3.*DDb3-b.*Db3.*Db3+c.*b3.*b3).*hh;
  diag(1:2:ndf-3)=diag(1:2:ndf-3)+(a.*DDb1.*DDb1-b.*Db1.*Db1+c.*b1.*b1).*hh;
  diag(4:2:ndf)=diag(4:2:ndf)+(a.*DDb4.*DDb4-b.*Db4.*Db4+c.*b4.*b4).*hh;
  diag(2:2:ndf-2)=diag(2:2:ndf-2)+(a.*DDb2.*DDb2-b.*Db2.*Db2+c.*b2.*b2).*hh;
  %
  subdiag(3:2:ndf-1)=subdiag(3:2:ndf-1)+(a.*DDb3.*DDb4-b.*Db3.*Db4+c.*b3.*b4).*hh;
  subdiag(1:2:ndf-3)=subdiag(1:2:ndf-3)+(a.*DDb1.*DDb2-b.*Db1.*Db2+c.*b3.*b4).*hh;
  subdiag(2:2:ndf-2)=subdiag(2:2:ndf-2)+(a.*DDb2.*DDb3-b.*Db2.*Db3+c.*b3.*b4).*hh;
  %
  sub2diag(1:2:ndf-3)=sub2diag(1:2:ndf-3)+(a.*DDb1.*DDb3-b.*Db1.*Db3+c.*b1.*b3).*hh;
  sub2diag(2:2:ndf-2)=sub2diag(2:2:ndf-2)+(a.*DDb2.*DDb4-b.*Db2.*Db4+c.*b2.*b4).*hh;
  %
  sub3diag(1:2:ndf-3)=sub3diag(1:2:ndf-3)+(a.*DDb1.*DDb4-b.*Db1.*Db4+c.*b1.*b4).*hh;
  %
   rhs(3:2:ndf-1)=rhs(3:2:ndf-1)+f.*b3.*hh;
   rhs(1:2:ndf-3)=rhs(1:2:ndf-3)+f.*b1.*hh;
   rhs(4:2:ndf)=rhs(4:2:ndf)+f.*b4.*hh;
   rhs(2:2:ndf-2)=rhs(2:2:ndf-2)+f.*b2.*hh;
  %
  end %gp
  %
 disp('Start of stiffness matrix:');
 [ diag(1)    subdiag(1) sub2diag(1) sub3diag(1) 0 ; ...
   subdiag(1)   diag(2)  subdiag(2)  sub2diag(2) sub3diag(2) ; ...
  sub2diag(1) subdiag(2)   diag(3)   subdiag(3)  sub2diag(3) ; ...
  sub3diag(1) sub2diag(2) subdiag(3)  diag(4)     subdiag(4) ; ...
       0      sub3diag(2) sub2diag(3) subdiag(4)   diag(5)
 ],
  %
  ndff=ndf;
  if bctype(2)==0;
% rhs=rhs-stiff(:,ndf-1)*u0(2)-stiff(:,ndf)*u1(2)*h;
   rhs(ndf-2)=rhs(ndf-2)-subdiag(ndf-2)*u0(2)-sub2diag(ndf-2)*u1(2)*h;
   rhs(ndf-3)=rhs(ndf-3)-sub2diag(ndf-3)*u0(2)-sub3diag(ndf-3)*u1(2)*h;
   rhs(ndf-4)=rhs(ndf-4)-sub3diag(ndf-4)*u0(2);
   diag=diag(1:ndf-2);
   subdiag=subdiag(1:ndf-3);sub2diag=sub2diag(1:ndf-4);sub3diag=sub3diag(1:ndf-5);
   ndff=ndf-2;
  end
  if bctype(1)==0;
% rhs=rhs-stiff(:,1)*u0(1)-stiff(:,2)*u1(1)*h;
   rhs(3)=rhs(3)-sub2diag(1)*u0(1)-subdiag(2)*u1(1)*h;
   rhs(4)=rhs(4)-sub3diag(1)*u0(1)-sub2diag(2)*u1(1)*h;
   rhs(5)=rhs(5)-sub3diag(2)*u1(1)*h;
  rhs=rhs(3:ndff);
  diag=diag(3:ndff);
  subdiag=subdiag(3:ndff-1);sub2diag=sub2diag(3:ndff-2);sub3diag=sub3diag(3:ndff-3);
  ndff=ndff-2;
  end
  %
 disp('Stiffness matrix after boundary cond.:');
 [ diag(1)    subdiag(1) sub2diag(1) sub3diag(1) 0 ; ...
   subdiag(1)   diag(2)  subdiag(2)  sub2diag(2) sub3diag(2) ; ...
  sub2diag(1) subdiag(2)   diag(3)   subdiag(3)  sub2diag(3) ; ...
  sub3diag(1) sub2diag(2) subdiag(3)  diag(4)     subdiag(4) ; ...
       0      sub3diag(2) sub2diag(3) subdiag(4)   diag(5)
 ],
  %
  %
  %
% Cholesky
%
%  [ d1  s1  ss1  sss1       ] [c1               ][c1 cs2 css3 csss4     ]
%  [ s1  d2   s2   ss2 sss2  ] [cs2  c2          ][   c2  cs3  css4 csss5]
%  [ss1  s2   d3    s3  ss3  ]=[css3 cs3 c3      ][       c3   cs4  css5 ]
%  ...
  endchol=0;
  rhs=[0 0 0 rhs];
  subdiag=[subdiag 0];sub2diag=[sub2diag 0 0];sub3diag=[sub3diag 0 0 0];
  sub1chol=zeros(size(1:ndff+1));  sub2chol=zeros(size(1:ndff+2));
  sub3chol=zeros(size(1:ndff+3));
  for j=1:ndff,
   chol(j)=diag(j)-sub1chol(j)^2-sub2chol(j)^2-sub3chol(j)^2;
   if chol(j)<0, endchol=1;break; end;
   chol(j)=sqrt(chol(j));
   sub1chol(j+1)=(subdiag(j)-sub1chol(j)*sub2chol(j+1)-sub2chol(j)*sub3chol(j+1) )/chol(j);
   sub2chol(j+2)=(sub2diag(j)-sub1chol(j)*sub3chol(j+2) )/chol(j);
   sub3chol(j+3)= sub3diag(j) /chol(j);
   rhs(j+3)=(rhs(j+3)-sub1chol(j)*rhs(j+2)-sub2chol(j)*rhs(j+1)  ...
                                 -sub3chol(j)*rhs(j) )/chol(j);
  end
  %
 disp('Start of Cholesky upper factor:');
 [ chol(1)    sub1chol(2) sub2chol(3) sub3chol(4) 0 ; ...
     0          chol(2)  sub1chol(3)  sub2chol(4) sub3chol(5) ; ...
     0           0         chol(3)   sub1chol(4)  sub2chol(5) ; ...
     0           0           0        chol(4)     sub1chol(5) ; ...
     0           0           0           0         chol(5)
 ],
  %
% u=stiff\rhs;
  rhs=rhs(4:ndff+3);
  u(ndff)=rhs(ndff)/chol(ndff);
  u(ndff-1)=(rhs(ndff-1)-sub1chol(ndff)*u(ndff))/chol(ndff-1);
  u(ndff-2)=(rhs(ndff-2)-sub1chol(ndff-1)*u(ndff-1)-sub2chol(ndff)*u(ndff))/chol(ndff-2);
  for j=ndff-3:-1:1,
   u(j)=(rhs(j)-sub1chol(j+1)*u(j+1)-sub2chol(j+2)*u(j+2)-sub3chol(j+3)*u(j+3))/chol(j);
  end
if bctype(2)==0;   u=[u u0(2) u1(2)*h];  end;
if bctype(1)==0;   u=[u0(1) u1(1)*h u ];  end;
  %
  %  plot
  %
  % Graphics preparations
  %
  for gp=1:3;
  x=g(gp);
  b1=eval(b1copy);
  b2=eval(b2copy);
  b3=eval(b3copy);
  b4=eval(b4copy);
  Db1=eval(Db1copy)/h;
  Db2=eval(Db2copy)/h;
  Db3=eval(Db3copy)/h;
  Db4=eval(Db4copy)/h;
  DDb1=eval(DDb1copy)/h^2;
  DDb2=eval(DDb2copy)/h^2;
  DDb3=eval(DDb3copy)/h^2;
  DDb4=eval(DDb4copy)/h^2;
  ucomp(gp:3:3*ni-3+gp)=u(1:2:2*ni-1)*b1+u(2:2:2*ni)*b2+u(3:2:2*ni+1)*b3+u(4:2:2*ni+2)*b4;
  Ducomp(gp:3:3*ni-3+gp)=u(1:2:2*ni-1)*Db1+u(2:2:2*ni)*Db2+u(3:2:2*ni+1)*Db3+u(4:2:2*ni+2)*Db4;
  DDucomp(gp:3:3*ni-3+gp)=u(1:2:2*ni-1)*DDb1+u(2:2:2*ni)*DDb2+u(3:2:2*ni+1)*DDb3+u(4:2:2*ni+2)*DDb4;
%
   integJ=integJ+sum(a.*DDucomp(gp:3:3*ni-3+gp).^2 ...
             -b.*Ducomp(gp:3:3*ni-3+gp).^2+c.*ucomp(gp:3:3*ni-3+gp))*hh;
  end
ucomp=[u(1) ucomp u(2*ni+1)];
 disp('Valeurs de u en 0.25, 0.5 et 0.75:');
 for j=1:3,
  xj=0.25*j;ixj=floor(xj/h);x=xj/h-ixj;
  b1=eval(b1copy);b2=eval(b2copy);b3=eval(b3copy);b4=eval(b4copy);
  uc(j)=u(2*ixj+1)*b1+u(2*ixj+2)*b2+u(2*ixj+3)*b3+u(2*ixj+4)*b4;
 end;
 format long;uc, format short;
  rn=1:ni;
  rn=rn*h;
  xplot(3:3:3*ni)=rn-ones(size(rn))*g(1)*h;
  xplot(2:3:3*ni-1)=rn-ones(size(rn))*g(2)*h;
  xplot(1:3:3*ni-2)=rn-ones(size(rn))*g(3)*h;
  xplot=[0 xplot 1];
    %Erase old plots?
  figure(gcf)
    if first
      clgwin=1;     %If first time clear window and plot uexact, uh
      first=~first;
    else
      disp('Do you want to clear the graph window?');
      clgwin = (input('answer y or n ','s')=='y');
    end
    if clgwin
  %   if uexactflag
  %     hold off;  plot(xxplot,uplot,':',xplot,ucomp,'-');  hold on;
  %     title(['u(x) = ' uexact ' dotted;     uh solid']);
  %   else
        hold off;  plot(xplot,ucomp,'-');  hold on;
        title('uh');
  %   end
    else
      plot(xplot,ucomp,'-');
    end
 disp(' Cholesky  min max :');
  format long;  [ min(chol(1:ndff)), max(chol(1:ndff))], format short ;
 disp(' integrale:');
  format long;   integJ, format short ;
 disp('Do you want to refine the mesh?');
 refine=(input('Answer y or n ','s')=='y');
 end; %of refinement while loop
hold  off
disp('If  you want a hardcopy of the plot, enter the command print');
disp('to close all: fclose(''all'')');
