function value = EvalClosureElement(xi1,xi2,xi3,u,N)

%
% value = EvalClosureElement(xi1,xi2,xi3,u,N)
%
% This function evaluates the value of the parameter u in a
% closure element with N perimeter nodes at the point
% (xi1,xi2,xi3).
%
% xi1 - the xi1 location
% xi2 - the xi2 location
% xi3 - the xi3 location
% u - the nodal parameters
% N - the number of perimeter nodes
%
% Mark Trew - July 2001
%

  Psi1 = CH(xi1);
  Psi2 = ClosCircBasis(xi2,N);
  Psi3 = CH(xi3);

  value = 0;
  for k=1:2,
    value = value + Psi1(1,1)*Psi2(1,1)*Psi3(k,1)*u(1+(k-1)*(N+1),1) + ... % u
                    Psi1(1,2)*Psi2(1,2)*Psi3(k,1)*u(1+(k-1)*(N+1),2) + ... % dudxi1
                    Psi1(1,2)*Psi2(1,3)*Psi3(k,1)*u(1+(k-1)*(N+1),3) + ... % dudxin
                    Psi1(1,2)*Psi2(1,3)*Psi3(k,1)*u(1+(k-1)*(N+1),5) + ... % d2udxi1dxin - always 0
                    Psi1(1,1)*Psi2(1,1)*Psi3(k,2)*u(1+(k-1)*(N+1),4) + ... % dudxi3
                    Psi1(1,2)*Psi2(1,2)*Psi3(k,2)*u(1+(k-1)*(N+1),6) + ... % d2udxi1dxi3
                    Psi1(1,2)*Psi2(1,3)*Psi3(k,2)*u(1+(k-1)*(N+1),7) + ... % d2udxindxi3
                    Psi1(1,2)*Psi2(1,3)*Psi3(k,2)*u(1+(k-1)*(N+1),8);      % d3udxi1dxindxi3 - 0??
  
    for p=2:N+1,
      value = value + Psi1(2,1)*Psi2(p,1)*Psi3(k,1)*u(p+(k-1)*(N+1),1) + ... % u
                      Psi1(2,2)*Psi2(p,1)*Psi3(k,1)*u(p+(k-1)*(N+1),2) + ... % dudxi1
                      Psi1(2,1)*Psi2(p,2)*Psi3(k,1)*u(p+(k-1)*(N+1),3) + ... % dudxi2
                      Psi1(2,2)*Psi2(p,2)*Psi3(k,1)*u(p+(k-1)*(N+1),5) + ... % d2udxi1dxi2
                      Psi1(2,1)*Psi2(p,1)*Psi3(k,2)*u(p+(k-1)*(N+1),4) + ... % dudxi3
                      Psi1(2,2)*Psi2(p,1)*Psi3(k,2)*u(p+(k-1)*(N+1),6) + ... % d2udxi1dxi3
                      Psi1(2,1)*Psi2(p,2)*Psi3(k,2)*u(p+(k-1)*(N+1),7) + ... % d2udxi2dxi3                  
                      Psi1(2,2)*Psi2(p,2)*Psi3(k,2)*u(p+(k-1)*(N+1),8);      % d3udxi1dxi2dxi3             
    end;
  end;

return;