function MyPLHeatEqSolver()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Solves the pde
%
%           du/dt + div q = f, q = - d grad u, 
%
%  with boundary condition
%
%             n q = k (u - g)
%
%  and initial condition
%
%             u = u0
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% defining geometry and meshing

  geom = geometry;

  [p,e,t] = initmesh(geom, 'HMax', .1);
%  [p,e,t] = refinemesh(geom, p, e, t);
  pdemesh(p, e, t); shg
  pause(1)


% defining data functions

  x1 = p(1,:);
  x2 = p(2,:);
  nn = size(p,2);
  d = ones(1,nn);
  m = ones(1,nn);
  f = zeros(1,nn);
  k = 100000 * ones(1,nn);
  g0 = x1.*x2.*(1-x2);              % preliminary, stationary part of  g
  u0 = zeros(nn,1);                 % initial values
  

% 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);
  U = u0;                            % storing


% 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;
    elLoad = f(elNodeInd);
    dF = elLoadVector(elCoords, elLoad);
    F(elNodeInd) = F(elNodeInd) + dF;
  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;
  u = u0;
  FinalTime = 3;
  dt = .03;

  A = M + dt/2 * D + dt/2 * K;
  L = M - dt/2 * D - dt/2 * K;


% time stepping

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

  while time < FinalTime

    G = zeros(nn,1);
    for bel = 1 : size(e, 2)
      belNodeInd = e(1: 2, bel);
      belCoords = p(:, belNodeInd);
      belPerm = k(belNodeInd);
      belLoad = g0(belNodeInd) * sin(2*pi*(time + dt/2));
      dG = elBdryVector(belCoords, belPerm, belLoad);
      G(belNodeInd) = G(belNodeInd) + dG;
    end


    b = L * U + dt * F + dt * G;
    U = A \ b;

    u = [u U];

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

  end

  delete(h)


% plotting

  set(gcf,'DoubleBuffer','on')
  for i = 1 : size(u, 2)
    pdeplot(p,e,t,'xydata',U,'zdata',u(:,i),'mesh','on');
    set(gca,'ZLim',[-.5 .5])
    xlabel('x_1')
    ylabel('x_2')
    zlabel('y')
    title(['solution  y = u(x_1, x_2, t),  t = ' num2str((i-1)*dt,'%1.1f')])
    drawnow
  end




% 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 = 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 geom = 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]';
  geom = [e1 e2 e3 e4];