function circle_test_dirichlet(h, curveit)
% circle_test_dirichlet(h,curveit)
% Solve a Dirichlet problem on a disk with a known analytic solution of
%  u = exp(x)*sin(y).
% Use quadratic 6-node triangle elements to solve this problem
% also using the finite element method.
% Print out the difference between the analytic solution and the
% finite element solution using the L^infty norm.
% Input arguments:
%   h = mesh size (passed to mesh generator)
%   curveit = 0 not to curve boundary elements
%   curveit = 1 to curve boundary elements
%

[xy,trilist,outeredge] = circle_t6_meshgen(h,curveit);
[p,scrap] = size(outeredge);
if scrap ~= 2
  error ('outeredge argument returned by mesh generator should have 2 columns');
end
  
[nn,scrap] = size(xy);
if scrap ~= 2
  error('xy argument returned by mesh generator should have 2 columns')
end

[ntri,scrap] = size(trilist);
if scrap ~= 6
  error('trilist argument returned by mesh generator should have 6 columns')
end

if any(outeredge(:,1)<1) | any(outeredge(:,1)>ntri)
  error('outeredge argument returned by mesh generator should have all entries in column 1 between 1 and ntri')
end

if any(outeredge(:,2)<1) | any(outeredge(:,2)>3)
  error('outeredge argument returned by mesh generator should have all entries in column 2 between 1 and 3')
end

if any(any(trilist<1)) | any(any(trilist>nn))
  error('trilist entries should be between 1 and nn')
end


% Set up dirichlet bc's
% Make a list of boundary nodes.
isbdry = logical(zeros(nn,1));
for j = 1 : p
  if outeredge(j,2) == 1
    vlist = [1,2,4];
  elseif outeredge(j,2) == 2
    vlist = [4,5,6];
  else
    vlist = [1,3,6];
  end
  isbdry(trilist(outeredge(j,1), vlist)) = [1;1;1];
end
bdrynodes = find(isbdry);
normxy = sqrt(xy(:,1).^2 + xy(:,2).^2);

xy2 = zeros(nn,2);
xy2(bdrynodes,:) = xy(bdrynodes,:) ./ (normxy(bdrynodes) * [1,1]);
val = exp(xy2(isbdry,1)) .* sin(xy2(isbdry,2));

dirichlet_bc = [find(isbdry),val];
neumann_bc = [];


numnode = nn;
numtri = ntri;




% Here are the 6 shape functions of xi,eta, one shape function per row
% of this matrix.

if nargin < 3
   sg = 0;
end


% Data for computing the element stiffness matrices.
% First are the six coefficients of each of the six
% shape functions.

shapecoef = [ 
%  1   xi   eta   xi^2 xi*eta eta^2
   0   0   -1      0     0     2
   0   0    4      0    -4    -4
   0   0    0      0     4     0
   1  -3   -3      2     4     2
   0   4    0     -4    -4     0
   0  -1    0      2     0     0
];

% Derivatives of shape functions
shapecoef_xi = [shapecoef(:,2), 2*shapecoef(:,4), shapecoef(:,5)];
shapecoef_eta = [shapecoef(:,3), shapecoef(:,5), 2*shapecoef(:,6)];

% Gauss quadrature points and weights for a triangle.

nq = 4;
Gqpoints = [
   1/3, 1/3
   0.6, 0.2
   0.2, 0.6
   0.2, 0.2
];

Gqweights = [
   -27/48
   25/48
   25/48
   25/48
] / 2;


% Make the global assembled stiffness matrix.
K = sparse([], [], [], numnode, numnode);

for j = 1 : numtri
   % set thistri to the vertex indices of triangle j (6-element vector)
   thistri = trilist(j,:);
   % Set rowlist to the list of rows spanned by this triangle 
   %  (6-element vector);
   vertices = xy(thistri, :);
   % compute the t6 laplacian element stiffness matrix (6-by-6)


   % Initialize element stiffness matrix to zeros
   Ke = zeros(6,6);

   % Loop on quadrature points
   for qn = 1 : nq
   
      % Get coordinates and weights of current point in ref triangle.
      xi = Gqpoints(qn,1);
      eta = Gqpoints(qn,2);
      qw = Gqweights(qn);
      v = [1;xi;eta];
   
      % Compute the Jacobian of the geometric mapping
      Jc = [vertices(:,1)' * shapecoef_xi * v, ...
            vertices(:,1)' * shapecoef_eta * v; ...
            vertices(:,2)' * shapecoef_xi * v, ...
            vertices(:,2)' * shapecoef_eta * v];
      detJc =  Jc(1,1) * Jc(2,2) - Jc(1,2) * Jc(2,1);
      if  sign(detJc) < 0
        warning('negative Jacobian sign -- inverted element')
     end
     detJc = abs(detJc);
     invJc = inv(Jc);
   
     % I and J loop over shape functions.
     for I = 1 : 6     
        dNI = [shapecoef_xi(I,:) * v, shapecoef_eta(I,:) * v] * invJc;
        for J = 1 : 6    
           dNJ = [shapecoef_xi(J,:) * v, shapecoef_eta(J,:) * v] * invJc;
	   Ke(I,J) = Ke(I,J) + (dNI(1) * dNJ(1) + dNI(2) * dNJ(2)) ...
		   * qw * detJc;
        end
     end
  end
  K = K + sparse(kron(thistri', ones(6,1)), ...
      kron(ones(6,1), thistri'), Ke(:), numnode, numnode);
end

% force vector, initially all zeros.
fdir = zeros(numnode,1);
numdisp = length(bdrynodes);


ud = zeros(numnode,1);

% Compute load vector due to Dirichlet BC's.

for k = 1 : numdisp
   vn = dirichlet_bc(k,1);
   val = dirichlet_bc(k,2);
   ud(vn) = val;
   fdir = fdir - full(K(:,vn)) * val;
end



% Form the reduced stiffness matrix (delete Dir. rows & cols) and solve.

f = fdir;
Kr = K(~isbdry,~isbdry);
fr = f(~isbdry);
ur = Kr \ fr;

u = zeros(numnode, 1);
u(~isbdry) = ur;
u(isbdry) = ud(isbdry);



uanalytic = exp(xy(:,1)) .* sin(xy(:,2));

disp(sprintf('max difference: u - uanalytic = %e', max(abs(u-uanalytic))))

% The xi-eta coordinates of the six nodes
nodal_xi_eta = [
  0 1
  0 .5
  .5 .5
  0 0
  .5 0
  1 0];

gradu = zeros(nn,2);
for j = 1 : ntri
  utri = u(trilist(j,:));
  xytri = xy(trilist(j,:), :);
  for vi = 1 : 6
    xi = nodal_xi_eta(vi,1);
    eta = nodal_xi_eta(vi,2);
    v = [1;xi;eta];
    vnum = trilist(j,vi);
    dudxi = utri' * shapecoef_xi * v;
    dudeta = utri' * shapecoef_eta * v;
    J = [xytri' * shapecoef_xi * v, xytri' * shapecoef_eta * v];
    gradu(vnum,:) = [dudxi,dudeta] * inv(J);
  end
end

graduanalytic = [uanalytic,exp(xy(:,1)) .* cos(xy(:,2))];

disp(sprintf('max difference: gradu - graduanalytic = %e', ...
  max(max(abs(gradu-graduanalytic)))))