function [ Qi, Gi ] = setNeumannBCOnFace( self, bci )
%setNeumannBCOnFace Calculate g and q entries for a Neumann BC
% This undocumented function may be changed or removed in a future release.
% This function handles the case when a 'g' and/or 'q' parameter
% has been included in a pdeBoundaryConditions entity.

%       Copyright 2014-2022 The MathWorks, Inc.


N2 = self.N*self.N;
neqn = self.N*self.numPoints;
appRegionID = bci.appRegionID;

isFuncQ = isa(bci.q, 'function_handle');
isFuncG = isa(bci.g, 'function_handle');
anyFunc = isFuncQ || isFuncG;
anyQ = isFuncQ || any(bci.q(:)); anyG = isFuncG || any(bci.g(:));

if(~ (anyQ || anyG))
  Qi = sparse(neqn,neqn);
  Gi = sparse(neqn,1);
  return;
end


if strcmp(bci.appRegionType,'Vertex')
    nodeID = bci.nodes;
    Qi = sparse(neqn,neqn);
    Gi = sparse(neqn,1);
    
    % Find the relevant DoF
    loadedDoF = nodeID + (0:self.dim-1)*self.numPoints;
    % Point load is aways packaged as function.
    uPoints = self.uN(:, nodeID);
    sPts = 0;
    xyz = self.points(:, nodeID);
    faceNormVec = ones(3, 1, 1);
    if isFuncG
        nodeG = self.callNeumannFuncOnFace(bci,xyz, sPts, bci.g, ...
            faceNormVec, uPoints, 1);
    else
        nodeG = bci.g(:);
    end
    Gi(loadedDoF) = nodeG;
    return
end

[numFaceNodes, numFaces] = size(bci.elemFaces);

% preallocate some space for sparse triplets
preAllocSize = 3*numFaces*self.N;
if(anyQ)
  qTriplets = zeros(3, preAllocSize);
  numQTriplets = 0;
end
if(anyG)
  gTriplets = zeros(2, preAllocSize);
  numGTriplets = 0;
end

% compute shape functions and derivatives at Gauss points
intWts = self.intRule.wts;
intPoints = self.intRule.points;
assert(self.dim>1)
dsPts = zeros(numFaceNodes, self.dim-1, self.numIntPoints);
for n=1:self.numIntPoints
  dsPts(:,:,n) = self.dsfCalc(intPoints(:,n));
end
sPts = self.sfCalc(intPoints);

if(anyG)
  Ge = zeros(numFaceNodes, self.N);
end

if(anyQ)
  Qe = zeros(numFaceNodes, numFaceNodes, N2);
end

if(~isempty(bci.q) && ~isFuncQ)
    qi = bci.q(:);
    if(any(qi) && length(qi) ~= N2)
        error(message('pde:pde3DBCImpl:invalidLenQ', appRegionID, N2));
    elseif(~any(qi) && length(qi) ~= N2 && length(qi) ~= 1)
        error(message('pde:pde3DBCImpl:invalidLenQ', appRegionID, N2));
    end
end
if(~isempty(bci.g) && ~isFuncG)
    gi = bci.g(:);
    if(any(gi) && length(gi) ~= self.N)
        error(message('pde:pde3DBCImpl:invalidLenG', numFaceNodes, self.N));
    elseif(~any(gi) && length(gi) ~= self.N && length(gi) ~= 1)
        error(message('pde:pde3DBCImpl:invalidLenG', numFaceNodes, self.N));
    end
end

%  calculate face normals
xyzAllFaceNodes = self.points(:, bci.elemFaces(:));
xyzAllFaceNodes = reshape(xyzAllFaceNodes, self.dim, numFaceNodes, numFaces);
elemLength = zeros(numFaces,1);
faceNorm = zeros(self.numIntPoints, numFaces);
if(anyFunc)
  uPoints = zeros(self.N, self.numIntPoints, numFaces);
  faceNormVec = zeros(3, self.numIntPoints, numFaces);
end
for i=1:numFaces
  xyz = xyzAllFaceNodes(:,:,i);
  if self.dim == 2
  elemLength(i) = getElemLength(xyz(:,1), xyz(:,end));
  end
  if(anyFunc)
    faceNodes = bci.elemFaces(:,i);
    uFaceNodes = self.uN(:, faceNodes);
    uPoints(:,:,i) = uFaceNodes*sPts;
  end
  for j=1:self.numIntPoints
    dxyzDrs = xyz*dsPts(:,:,j);
    if self.dim == 3
        v3 = cross3(dxyzDrs(:,1),dxyzDrs(:,2));
    else
        v3 = cross3([dxyzDrs;0],[0,0,1]);
    end
    faceNorm(j,i) = norm(v3);
    if(anyFunc)
      faceNormVec(:,j,i) = v3/faceNorm(j,i);
    end
  end
end
if(isFuncQ)
  faceQ = self.callNeumannFuncOnFace(bci,xyzAllFaceNodes, sPts, bci.q, ...,
    faceNormVec, uPoints, 2);
end
if(isFuncG)
  faceG = self.callNeumannFuncOnFace(bci,xyzAllFaceNodes, sPts, bci.g, ...
    faceNormVec, uPoints, 1);
end

fCol = 1;
% Iterate over all element faces in this BC
for i=1:numFaces
  coni = bci.elemFaces(:,i);
  % Iterate over the integration points
  for n=1:self.numIntPoints 
    shn = sPts(:,n);
    detJ = faceNorm(n, i); 
    if self.dim == 2
        detJ = elemLength(i);
    end
    wt = intWts(n);
    detWt = detJ*wt;
    if(anyG)
      %bci.g
      if(isFuncG)
        gi = faceG(:,fCol);
      end
      detShp = detWt*shn;
      for j=1:self.N
        Ge(:,j) = Ge(:,j) + gi(j)*detShp;
      end
    end
    if(anyQ)
      if(isFuncQ)
        qi = faceQ(:,fCol);
      end
      detShpShp = detWt*(shn*shn.');
      for j=1:N2
        Qe(:,:,j) = Qe(:,:,j) +qi(j)*detShpShp;
      end
    end
    fCol = fCol + 1;
  end
  if(anyG)
    % copy the face g-vector to triplets
    for j=1:self.N
      offset = (j-1)*self.numPoints;
      gj = Ge(:,j);
      nzInd = gj ~= 0;
      nnz = sum(nzInd);
      gTriplets(:, numGTriplets+1:numGTriplets+nnz) = [coni(nzInd).'+offset; gj(nzInd).'];
      numGTriplets = numGTriplets + nnz;
    end
    % set Ge back to zero for the next face
    Ge(:) = 0;
  end
  if(anyQ)
    % copy the face q-matrix to triplets
    jj = 1;
    for j1=1:self.N
      offset1 = (j1-1)*self.numPoints;
      for j2=1:self.N
        offset2 = (j2-1)*self.numPoints;
        qjk = reshape(Qe(:, :, jj).', [], 1);
        nzInd = qjk ~= 0;
        nnz = sum(nzInd);
        if(nnz)
          gRows = repmat(coni.' + offset2, numFaceNodes, 1); gRows = gRows(:);
          gCols = repmat(coni.' + offset1, 1, numFaceNodes);
          qTriplets(:,numQTriplets+1:numQTriplets + nnz) = [gRows(nzInd).'; gCols(nzInd); qjk(nzInd).'];
          numQTriplets = numQTriplets + nnz;
        end
        jj = jj + 1;
      end
    end
    % set Qe back to zero for the next face
    Qe(:) = 0;
  end
end

if(anyG && numGTriplets)
  Gi = sparse(gTriplets(1, 1:numGTriplets), 1, gTriplets(2, 1:numGTriplets), ...
    neqn,1);
else
  Gi = sparse(neqn,1);
end
if(anyQ && numQTriplets)
  Qi = sparse(qTriplets(1, 1:numQTriplets), qTriplets(2, 1:numQTriplets), ...
    qTriplets(3, 1:numQTriplets), neqn, neqn);
else
  Qi = sparse(neqn,neqn);
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

function elmLength=getElemLength(p1,p2)
% compute length of a line segment given two end point coordinates.
elmLength = sqrt((p2(1)-p1(1))^2 + (p2(2)-p1(2))^2);
end

