function MyPoissonSolver()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Solves the pde
%
%           div q = f, q = - d grad u, 
%
%  with boundary conditions
%
%             n q = k (u - g)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% 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,:);
  f = 2*sin(pi*x1).*sin(pi*x2)*pi^2;
  nn = size(p,2);                        % number of nodes/degrees of freedom
  d = ones(1,nn);
  k = 100000 * ones(1,nn);               % large k to force u = g
  g = zeros(1,nn);
  

% initialising matrises


  
  D = sparse(zeros(nn,nn));
  K = sparse(zeros(nn,nn));
  F = zeros(nn,1);
  G = zeros(nn,1);


% assembling element contributions

  for el = 1 : size(t,2)
    elNodeInd = t(1: 3, el);
    elCoords = p(:,elNodeInd);             % element vertex coords (2 by 3 matrix)
    elDiff = d(elNodeInd);                 % element "diffusivity" (1 by 3 matrix)
    dD = elDiffMatrix(elCoords, elDiff);
    D(elNodeInd, elNodeInd) = D(elNodeInd, elNodeInd) + dD;
    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);           % bdry elem coords (2 by 2)
    belPerm = k(belNodeInd);                % bdry elem "permeability"
    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


% solving

  U = (D + K)\(F + G);


% plotting

  h = pdeplot(p, e, t, 'xydata', U, 'zdata', U, 'mesh', 'on');
  xlabel('x_1')
  ylabel('x_2')
  zlabel('y')
  grid on
  title('solution  y = u ( x_1, x_2)')


% test

%  u = sin(pi * x1) .* sin(pi * x2);  % exact sol. if domain changed to the unit square
%  maxerr = max( abs(U - u') )





% subroutines


function dD = elDiffMatrix(coords, diffusion)

% computes the local stiffness matrix on a triangle with
% vertex coordinates coords (2 by 3 matrix) and 
% vertex diffusivity diffusion (1 by 3 matrix)

  dx = area(coords);
  x = sum(coords, 2) / 3;                 % element midpoint (center of gravity)
  Dph = Dphi(coords);
  d = sum(diffusion)/length(diffusion);   % midpoint diffusivity
  dD = Dph' * (d * Dph) * dx;


function dF = elLoadVector(coords, domainLoad)

  dx = area(coords);
  x = sum(coords, 2) / 3;
  phi = [1/3 1/3 1/3];      % element midpoint values of the three basis functions
  f = sum(domainLoad)/length(domainLoad);
  dF = phi' * f * dx;


function dK = elBdryMatrix(coords, bdryPerm)

% computes the local "mass" matrix of a boundary element
% (line segment) with coordinates coords (2 by 2 matrix)
% and "permeability" bdryPerm (1 by 2 matrix)

  v = coords(:, 2) - coords(:, 1);
  ds = sqrt(v' * v);           % the (arc) length of the bdry element
  x = sum(coords, 2) / 2;      % element midpoint
  phi = [1/2 1/2];             % midpoint value of the two local basis functions
  k = sum(bdryPerm)/length(bdryPerm);   % midpoint k value
  dK = phi' * (k * phi) * ds;


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)

% computes the area of a triangle with vertex coordinates coords
% (stored in a 2 by 3 matrix)

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


function Dph = Dphi(coords)

% computes the gradients of the three local basis functions on a
% triangle with coordinates 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
  Dph(:, 3) = n2v / pw;

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

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



% data functions


function geom = geometry()

% define bdry consisting of 4 line segments
% (indicated by the initial numbers 2, see $>>$ help decsg),
% the first of which extends from 0 to 1 in the x direction
% and from 0 to -1 in the y direction, while the second
% goes from 1 to 2 in the x direction, and from -1 to 1 in y, etc,
% and giving number 1 to the domain to "the left" of the lines
% while the last 0's means that the domain to "the right" of
% the lines is the "exterior" (not part of the computational
% domain).

  e1 = [2 0 1 0 -1 1 0]';
  e2 = [2 1 2 -1 1 1 0]';
  e3 = [2 2 0 1 1 1 0]';
  e4 = [2 0 0 1 0 1 0]';
  geom = [e1 e2 e3 e4];       % the bdry edge data stored as
                              % columns in a geometry matrix geom