function femodel = imposeBCOnAxis(femodel)
% Explicitly enforce zero-dirichlet BC for magnetostatic case, as the
% conversion factor 1/r does not tend to zero at the axis of rotation.

%       Copyright 2022 The MathWorks, Inc.
nodesOnAxis = femodel.p(1,:)<100*eps;
numAxiSymmCons = nnz(nodesOnAxis);

if numAxiSymmCons == 0
    % No node on axis.
    return
end

if (size(femodel.Or,2)==0)
    % Model does not have any dirichlet BC and has an edge/vertex on
    % axis. Dirichlet BCs are implicitly enforced on nodes at axis.
    H = [];
    R = [];
else
    % Model has some dirichlet BC, these BCs may or may not be on the axis.
    % Append constraints to enforces zero dirichlet BC on the nodes at axis.
    H = femodel.H;
    R = femodel.R;
end
[numH, numColsH] = size(H);
if numColsH == 0
    % ensure size of H is correct if no BCs are present
    numColsH = size(femodel.K,2);
end
rowIdx = numH+1:numH+numAxiSymmCons;
H = [H; zeros(numAxiSymmCons, numColsH)];
H(rowIdx,nodesOnAxis) = eye(numAxiSymmCons);
R = [R;zeros(numAxiSymmCons,1)];
% if H has extra constraints for one dof (possible if BC specification
% intersects axisymmetric edge) remove BC and impose axisymmetric part
overconstrainedDofs = find(sum(H,1)>1);
for dof = 1:numel(overconstrainedDofs)
    row = find(H(:,overconstrainedDofs(dof)));
    H(row(1:end-1),:) = [];
    R(row(1:end-1)) = [];
end
femodel.H = H;
femodel.R = R;
[femodel.Nu,femodel.Or]=pdenullorth(H);
femodel.ud=full(femodel.Or*((H*femodel.Or)\R));
[femodel.totalNumEqns, femodel.numConstrainedEqns]=size(femodel.Nu);
femodel.B=femodel.Nu;
end