%--------------------------                            
% program Winkler.m
%--------------------------                            
 clear

%---------------- 
% INPUT DATA
%---------------- 
  load MESH   -ASCII
  load QUAD   -ASCII

% ---------------- 
% Define MESH Data
% ---------------- 
  NUMNP = MESH(1,1); 
  NUMEL = NUMNP-1;
  NUMEQ = 2*NUMNP;
  IB    = 4;
  for i=1:NUMNP;
     XORD(i)   = MESH(i+1,1);
     NPcode(i) = MESH(i+1,2);
  end

% ---------------- 
% Define QUAD Data
% ---------------- 
  NQPTS = QUAD(1,1);
  for i=1:NQPTS
     GPTS(i)=QUAD(i+1,1);
     GWTS(i)=QUAD(i+1,2);
  end

% --------------------------- 
% Declare sparse matrix
% --------------------------- 
  SK  = sparse(NUMEQ,NUMEQ);

% --------------------------- 
% General initialization
% --------------------------- 
  LHS(1:NUMEQ,1)=0.0;
  RHS(1:NUMEQ,1)=0.0;
  NPBC(1:NUMNP)=0.0;

% -------------------------------- 
% Include user's initialization
% -------------------------------- 
  INITIAL

% -------------------------------- 
% FORMATION OF STIFFNESS MATRIX   
% and RIGHT HAND SIDE
% -------------------------------- 

  for I=1:NUMEL
%    ----------------------------
%    Initialize Element Variables
%    ----------------------------
     S1(1:4,1:4)=0.0;
     S2(1:4,1:4)=0.0;
     Qe(1:4,1)=0.0;
%    ----------------------------
%    Define rows and columns for
%    this element
%    ----------------------------
     rc=[2*I-1,2*I,2*I+1,2*I+2];

% ----------------------------
% Calculate element coordinate
% information.
% ----------------------------
     Xa = XORD(I);
     Xb = XORD(I+1);
     RL=Xb-Xa; 
     DxDu=RL/2.0; 
     DuDx=2.0/RL;   

%    ------------------------- 
%    Begin Gaussian Quadrature 
%    ------------------------- 
     for J=1:NQPTS
       u  = GPTS(J);
       Wt = GWTS(J);
%      --------------------
%      Global coordinate of 
%      current Gauss point
%      --------------------
       Xg = (Xa+Xb)/2  + (RL/2.0)*u;

%      ----------------------------
%      INCLUDE COEF.m 
%      Defines: EIx, Kx, and Qx
%      ----------------------------
       COEF
%      ------------------------------
%      Calculate  shape functions wrt 
%      x at current Gauss point
%      ------------------------------
       No(1) = SF(0,1,u);
       No(2) = SF(0,2,u)*(DxDu);
       No(3) = SF(0,3,u);
       No(4) = SF(0,4,u)*(DxDu);

       Nxx(1) = SF(2,1,u)*(DuDx)^2;
       Nxx(2) = SF(2,2,u)*(DuDx)^2*(DxDu);
       Nxx(3) = SF(2,3,u)*(DuDx)^2;
       Nxx(4) = SF(2,4,u)*(DuDx)^2*(DxDu);

%      ----------------------------
%      Element stiffness matrices
%      ----------------------------
       Wt = GWTS(J);
       S1=S1+Wt*Nxx'*EIx*Nxx*DxDu;
       S2=S2+Wt*No'*Kx*No*DxDu;
       Qe=Qe+Wt*No'*Qx*DxDu;
     end   
%    ----- Quadrature now complete 

%    --------------------------
%    Assemble element matrices 
%    into global matrix 
%    -------------------------- 
     SK(rc,rc)=SK(rc,rc)+S1+S2;
     RHS(rc)=RHS(rc)+Qe;

   end
%  ---- Global Matrices are assembled

% -------------------------
% BOUNDARY CONDITIONS
% -------------------------
  B = 1.0E+06;
  for I=1:NUMNP
   if NPBC(I) == 1 | NPBC(I) == 3
     I1 = 2*I-1; 
     SK(I1,I1)=SK(I1,I1)*B;
     RHS(I1)=LHS(I1)*SK(I1,I1);
   end
   if NPBC(I) == 2 | NPBC(I) == 3
     I2=2*I;
     SK(I2,I2)=SK(I2,I2)*B;
     RHS(I2)=LHS(I2)*SK(I2,I2);
   end
  end

% -------------------------
% Solve equations:
% -------------------------
  LHS = SK\RHS;

% -------------------------------------- 
% Nodal values for w and dw/dx % are now 
% in LHS.  Use these values to calculate 
% shear and moment at center of each element.
%
% First calculate shape functions
% and their derivatives.  Ic is
% counter for number elements used.
% ------------------------------------ 
  Ic = 0;
  for I=1:NUMEL
    Xa = XORD(I);
    Xb = XORD(I+1);
    RL=Xb-Xa; 
    Xg=(Xa+Xb)/2;
    DuDx=2.0/RL;   
    DxDu=RL/2.0;   

    No(1) = SF(0,1,0);
    No(2) = SF(0,2,0)*(DxDu);
    No(3) = SF(0,3,0);
    No(4) = SF(0,4,0)*(DxDu);

    Nx(1) = SF(1,1,0)*DuDx;;
    Nx(2) = SF(1,2,0)*(DxDu)*DuDx;
    Nx(3) = SF(1,3,0)*DuDx;
    Nx(4) = SF(1,4,0)*(DxDu)*DuDx;

    Nxx(1)   = SF(2,1,0)*DuDx^2;
    Nxx(2)   = SF(2,2,0)*DxDu*DuDx^2;
    Nxx(3)   = SF(2,3,0)*DuDx^2;
    Nxx(4)   = SF(2,4,0)*DxDu*DuDx^2;

    Nxxx(1)   = SF(3,1,0)*DuDx^3;
    Nxxx(2)   = SF(3,2,0)*DxDu*DuDx^3;
    Nxxx(3)   = SF(3,3,0)*DuDx^3;
    Nxxx(4)   = SF(3,4,0)*DxDu*DuDx^3;

% --------------------------- 
%  Calculate deflection, slope
%  shear and moment at center
%  of element.
% --------------------------- 
    COEF

    W=0;
    S=0;
    V=0;
    M=0;
    for J=1:4
       W = W + No(J)*LHS(2*I-2+J);
       S = S + Nx(J)*LHS(2*I-2+J);
       M = M + EIx*Nxx(J)*LHS(2*I-2+J);
       V = V + EIx*Nxxx(J)*LHS(2*I-2+J)...
                +dEIx*Nxx(J)*LHS(2*I-2+J);
    end
   
% --------------------------- 
%  Save center values except 
%  for a pin-element.
% --------------------------- 
    if EIx  > 0
      Ic = Ic + 1;
      ANS(Ic,1)  =  Xg;
      ANS(Ic,2)  =  W;
      ANS(Ic,3)  =  S;
      ANS(Ic,4)  =  M;
      ANS(Ic,5)  =  V;
    end

  end
% --------------------------- 
% Solution is now complete.
% Save values.
% --------------------------- 

  save ANS ANS -ASCII

 NOTE = ...
 ['Solution is saved in file:   '; ...
  '        ANS                  '; ...
  'where the columns  are:      '; ...
  '                             '; ...
  '[ x | w | dw/dx | M | V ]    '; ...
  '                             '; ...
  'To plot, e.g.  w vs x,  type '; ...
  'plot(ANS(:,1),ANS(:,2))      ']

