function MyPLWaveEqSolver()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Solves the pde system
%
%     du1/dt = u2,  du2/dt + div ( - grad u1) = 0
%
%  with boundary conditions
%
%            n qi = ki (ui - gi),
%
%  and initial conditions
%
%            ui = ui0,  i = 1,2.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% defining geometry and meshing

  g = geometry;

  [p,e,t] = initmesh(g,'HMax',.1);


% defining data functions

  x1 = p(1,:);
  x2 = p(2,:);
  nn = size(p,2);
  d = .1*ones(1,nn);
  m = ones(1,nn);
  k = zeros(1,nn);
  g = zeros(1,nn);        % has no effect because k1 = 0


% initialising matrises
  
  D = sparse(zeros(nn,nn));
  M = sparse(zeros(nn,nn));
  K = sparse(zeros(nn,nn));
  F = zeros(nn,1);
  G = zeros(nn,1);
  U1 = .2*sin(2*pi*x1.*x2)';
  U2 = zeros(nn,1);
  O = zeros(nn,1);


% assembling element contributions

  for el = 1 : size(t,2)
    elNodeInd = t(1: 3, el);
    elCoords = p(:,elNodeInd);
    elDiff = d(elNodeInd);
    dD = elDiffMatrix(elCoords, elDiff);
    D(elNodeInd,elNodeInd) = D(elNodeInd,elNodeInd) + dD;
    elMass = m(elNodeInd);
    dM = elMassMatrix(elCoords, elMass);
    M(elNodeInd,elNodeInd) = M(elNodeInd,elNodeInd) + dM;
  end


% assembling boundary element conributions

  for bel = 1 : size(e, 2)
    belNodeInd = e(1: 2, bel);
    belCoords = p(:, belNodeInd);
    belPerm = k(belNodeInd);
    dK = elBdryMatrix(belCoords, belPerm);
    K(belNodeInd, belNodeInd) = K(belNodeInd, belNodeInd) + dK;
    belLoad = g(belNodeInd);
    dG = elBdryVector(belCoords, belPerm, belLoad);
    G(belNodeInd) = G(belNodeInd) + dG;
  end


% initializations

  time = 0;
  u1 = U1;                             % storing
  u2 = U2;                             % storing
  FinalTime = 2;
  dt = .01;
  ekin = [];
  epot = [];
  etot = [];

  A = [M, -dt/2*M;  dt/2*D+dt/2*K, M+0*dt/2*M];
  L = [M,  dt/2*M; -dt/2*D-dt/2*K, M-0*dt/2*M];
  W = [U1; U2];


% time stepping

  h = waitbar(0,'Just a minute..');

  while time < FinalTime

     F = zeros(nn, 1);
%    F1 = zeros(nn,1);
%    F2 = zeros(nn,1);
    for el = 1 : size(t, 2)
      elNodeInd = t(1: 3, el);
      elCoords = p(:, elNodeInd);
%      f1 = U2(elNodeInd);
      f = -U2(elNodeInd);
%      f2 = 0;
%      dF1 = elLoadVector(elCoords, f1);
%      dF2 = elLoadVector(elCoords, f2);
%      F1(elNodeInd) = F1(elNodeInd) + dF1;
%      F2(elNodeInd) = F2(elNodeInd) + dF2;
    end

    W = [U1; U2];
    b = L * W + dt * [O; F] + dt * [O; G];
    W = A \ b;
    U1 = W(1:nn);;
    U2 = W(nn+1:2*nn);

% checking energy
    ekin = [ekin U2'*M*U2];
    epot = [epot U1'*D*U1];
    etot = [etot U2'*M*U2+U1'*D*U1];
    u1 = [u1 U1];
    u2 = [u2 U2];

    time = time + dt;
    waitbar(time/FinalTime)

  end

  delete(h)


% plotting

  set(gcf,'DoubleBuffer','on')
  for i = 1 : size(u1, 2)
    subplot(121)
    pdeplot(p,e,t,'xydata',U1,'zdata',u1(:,i),'mesh','on');
    set(gca,'ZLim',[-1 1])
    grid on
    zlabel('displacement')
    title(['time = ' num2str((i-1)*dt)])
    subplot(122)
    pdeplot(p,e,t,'xydata',U1,'zdata',u2(:,i),'mesh','on');
    set(gca,'ZLim',[-1 1])
    grid on
    zlabel('velocity')
    title(['time = ' num2str((i-1)*dt)])
    drawnow
  end
  
  plot(dt*(1:length(etot)),[ekin;epot;etot])
  xlabel('time')
  title('energies, total red, pot. green, kin. blue')

% subroutines


function dD = elDiffMatrix(coords, diffusion)

  dx = area(coords);
  x = sum(coords, 2) / 3;
  Dphi = Dphi(coords);
  d = sum(diffusion)/length(diffusion);
  dD = Dphi' * (d * Dphi) * dx;


function dM = elMassMatrix(coords, mass)

  dx = area(coords);
  x = sum(coords, 2) / 3;
  phi = [1/3 1/3 1/3];
  m = sum(mass)/length(mass);
  dM = phi' * (m * phi) * dx;


function dF = elLoadVector(coords, domainLoad)

  dx = area(coords);
  x = sum(coords, 2) / 3;
  phi = [1/3 1/3 1/3];
%  f = feval(Myf, u1, u2);
  f = sum(domainLoad)/length(domainLoad);
  dF = phi' * f * dx;


function dK = elBdryMatrix(coords, bdryPerm)

  v = coords(:, 2) - coords(:, 1);
  ds = sqrt(v' * v);
  p = 0.57735026918963;
  points = [.5-p/2; .5+p/2];
  weights = [.5 .5];
  x = sum(coords, 2) / 2;
  phi = [1/2 1/2];
  phi1 = [1-points(1) points(1)];
  phi2 = [1-points(2) points(2)];
  k = sum(bdryPerm)/length(bdryPerm);
  dK = phi1' * (k * phi1) * ds / 2 ...
     + phi2' * (k * phi2) * ds / 2;


function dG = elBdryVector(coords, bdryPerm, bdryLoad)

  v = coords(:, 2) - coords(:, 1);
  ds = sqrt(v' * v);
  x = sum(coords, 2) / 2;
  phi = [1/2 1/2];
  k = sum(bdryPerm)/length(bdryPerm);
  g = sum(bdryLoad)/length(bdryLoad);
  dG = phi' * k * g * ds;  


function dx = area(coords)

  v = coords(:, 2) - coords(:, 1);
  w = coords(:, 3) - coords(:, 1);
  dx = (v(1) * w(2) - v(2) * w(1)) / 2;


function Dphi = Dphi(coords)

  v = coords(:, 2) - coords(:, 1);
  w = coords(:, 3) - coords(:, 1);

  n2v = [-v(2) v(1)]';                    % normal to v
  n2v = n2v / sqrt(n2v' * n2v);           % normalizing
  pw = w' * n2v;                          % height of triangle
  Dphi(:, 3) = n2v / pw;

  n2w = [w(2) -w(1)]';
  n2w = n2w / sqrt(n2w' * n2w);
  pv = v' * n2w;
  Dphi(:, 2) = n2w / pv;

  Dphi(:, 1) = - Dphi(:, 2) - Dphi(:, 3);



% data functions


function g = geometry()
  e1 = [2 0 1 0 0 1 0]';
  e2 = [2 1 1 0 1 1 0]';
  e3 = [2 1 0 1 1 1 0]';
  e4 = [2 0 0 1 0 1 0]';
  g = [e1 e2 e3 e4];