function faceNormVecAtNodes = calcFaceNormalsAtNodes(self, faces, nodes)
%calcFaceNormalsAtNodes Calculate normal to a face at each of its nodes
% This undocumented function may be changed or removed in a future release.
% Given a set of faces and face nodes, this function calculates the face 
% normal at each node. 

%       Copyright 2014-2020 The MathWorks, Inc.

import pde.internal.elements.*;

[numFaceNodes, numFaces] = size(faces);
numNodes = length(nodes);
numFacesAtNode = zeros(numNodes, 1);
faceNormVecAtNodes = zeros(3, numNodes);

if isempty(faces) % facets empty indicates vertex bc. Don't compute normals..
    return;
end

% Collect the shape function derivatives for all face nodes
dsNodes = zeros(numFaceNodes, self.dim-1, numFaceNodes);
 for n=1:numFaceNodes
  dsNodes(:,:,n) = self.dsfCalc(self.rs(:,n));
end

nodeIndexVec = zeros(self.numPoints, 1);
nodeIndexVec(nodes) = 1:numNodes;

xyzAllFaceNodes = self.points(:, faces(:));
xyzAllFaceNodes = reshape(xyzAllFaceNodes, self.dim, numFaceNodes, numFaces);
faceNormVec = zeros(3, numFaceNodes);
for i=1:numFaces
  xyz = xyzAllFaceNodes(:,:,i);
  for j=1:self.numNormalsPerFacet
    % Compute the two vectors dX/dr and dX/ds
    dxyzDrs = xyz*dsNodes(:,:,j);
    % Their cross product is the normal at the point
    if self.dim == 3
        v3 = cross3(dxyzDrs(:,1),dxyzDrs(:,2));
    else
        v3 = cross3([dxyzDrs;0],[0,0,1]);
    end
    faceNormVec(:,j) = v3/norm(v3);
  end
  if(self.numNormalsPerFacet==1)
    faceNormVec = repmat(faceNormVec(:,1), 1, numFaceNodes);
  end
  for j=1:numFaceNodes
    nj = nodeIndexVec(faces(j,i));
    numFacesAtNode(nj) = numFacesAtNode(nj) + 1;
    faceNormVecAtNodes(:,nj) = faceNormVecAtNodes(:,nj) + faceNormVec(:,j);
  end
end

% The normal at a node is the average of the normals of all the facets
% connected to that node.
for i=1:numNodes
  faceNormVecAtNodes(:,i) = faceNormVecAtNodes(:,i)/numFacesAtNode(i);
end

end

% This highly specialized version of the cross function with no error
% checking is used to improve performance.
function c=cross3(a,b) 
c=[a(2)*b(3)-a(3)*b(2); 
   a(3)*b(1)-a(1)*b(3); 
   a(1)*b(2)-a(2)*b(1)];
end

