function f=pdehypf(t,u)
%PDEHYPF Right hand side for HYPERBOLIC.

%       Copyright 1994-2014 The MathWorks, Inc.

 global pdehyp
  
  if ~(pdehyp.vq || pdehyp.vg || pdehyp.vh || pdehyp.vr || ...
      pdehyp.vc || pdehyp.va || pdehyp.vf)
    f=-pdehyp.K*u+pdehyp.F;
    return
  end
  
  nu=pdehyp.numConstrainedEqns;
  if(nu == pdehyp.totalNumEqns)
    uFull = u(1:nu);
  else
    uFull = pdehyp.B*u(1:nu) + pdehyp.ud;
  end
  
  if pdehyp.vq || pdehyp.vg || pdehyp.vh || pdehyp.vr
    [pdehyp.Q,pdehyp.G,pdehyp.H,pdehyp.R]=assemb(pdehyp.b,pdehyp.p,pdehyp.e,uFull,t);
  end
  
  if pdehyp.vc || pdehyp.va || pdehyp.vf
    [pdehyp.K,pdehyp.M,pdehyp.F]=pdehyp.thePde.assema(pdehyp.p,pdehyp.t,pdehyp.c,pdehyp.a,pdehyp.f,uFull,t);
  end
 
  if ~(pdehyp.vh || pdehyp.vr)
    nu=pdehyp.numConstrainedEqns;
    KK = pdehyp.K + pdehyp.M + pdehyp.Q;
    K=[sparse(nu,nu) -speye(nu,nu); pdehyp.B'*KK*pdehyp.B sparse(nu,nu)];
    FF=pdehyp.B'*(pdehyp.F +pdehyp.G - KK*pdehyp.ud);
    F=[zeros(nu,1); FF];
  else
    dt=pdehyp.ts*eps^(1/3);
    t1=t+dt;
    dt=t1-t;
    tm1=t-dt;
    if pdehyp.vd
      [~,pdehyp.MM]=pdehyp.thePde.assema(pdehyp.p,pdehyp.t,0,pdehyp.d,zeros(pdehyp.N,1),uFull,t);
    end
    if(isempty(pdehyp.C))
      nrFull = size(uFull,1);
      dampForce = zeros(nrFull,1);
      dampMat = sparse(nrFull, nrFull);
    end
    KK = pdehyp.K + pdehyp.M + pdehyp.Q;
    FF = pdehyp.F+pdehyp.G;
    if pdehyp.vh
      [N,O]=pdenullorth(pdehyp.H);
      ud=O*((pdehyp.H*O)\pdehyp.R);
      [Q,G,H,R]=assemb(pdehyp.b,pdehyp.p,pdehyp.e,uFull,t1);
      [N1,O]=pdenullorth(H);
      ud1=O*((H*O)\R);
      [Q,G,H,R]=assemb(pdehyp.b,pdehyp.p,pdehyp.e,uFull,tm1);
      [Nm1,O]=pdenullorth(H);
      udm1=O*((H*O)\R);
      nu=size(N,2);
      if(~isempty(pdehyp.C))
        dampMat = pdehyp.C;
        udDot = (ud1 - udm1)/(2*dt);
        dampForce = dampMat*udDot;
      end
      BDot = (N1-Nm1)/(2*dt);
      BDotDot = (N1-2*N+Nm1)/(dt^2);
      udDotDot = (ud1-2*ud+udm1)/(dt^2);
      K=[sparse(nu,nu) -speye(nu,nu); ...
        N'*(KK*N+pdehyp.MM*BDotDot+dampMat*BDot) ...
        2*N'*pdehyp.MM*BDot];
      F=[zeros(nu,1); N'*(FF-...
                     KK*ud - ...
                      dampForce - pdehyp.MM*udDotDot)];
    else
      nu=size(pdehyp.Nu,2);
      HH=pdehyp.H*pdehyp.Or;
      ud=pdehyp.Or*(HH\pdehyp.R);
      [Q,G,H,R]=assemb(pdehyp.b,pdehyp.p,pdehyp.e,uFull,t1);
      ud1=pdehyp.Or*(HH\R);
      [Q,G,H,R]=assemb(pdehyp.b,pdehyp.p,pdehyp.e,uFull,tm1);
      udm1=pdehyp.Or*(HH\R);
      if(~isempty(pdehyp.C))
        udDot = (ud1 - udm1)/(2*dt);
        dampForce = pdehyp.C*udDot;
      end
      udDotDot = (ud1-2*ud+udm1)/(dt^2);
      K=[sparse(nu,nu) -speye(nu,nu); ...
        pdehyp.Nu'*KK*pdehyp.Nu sparse(nu,nu)];
      F=[zeros(nu,1); ...
        pdehyp.Nu'*(FF - KK*ud - ...
        dampForce - pdehyp.MM*udDotDot)];
    end
    pdehyp.ud = ud;
  end
  
  f=-K*u+F;
  
end
