function u = solveNedelecHarmonic(self,coefstruct,frequency)
% Internal only.
% function called in solve when the harmonic analysis type is selected.

%   Copyright 2022 The MathWorks, Inc.

[edges, edge_elem] = constructEdgeElementConnectivity(self);
numEdges = size(edges,1);

% construct the discretized model and apply boundary conditions
% if using magnetic field, update coefstruct to take curl of source
if isequal(self.FieldType,"magnetic")
    coefstruct.ComputeCurlSource = true;
end
[M,K,F,H,R,~] = discretizedModelHarmonicNedelec(self,coefstruct,edges,edge_elem,numEdges,frequency);
A = M+K;
[B,Or] = pdenullorth(H);
ud = Or*((H*Or)\R);
Ac = B'*A*B;
Fc = B'*(F - A*ud);

% solve the system of equations
sol = Ac\Fc;
sol = B*sol(1:size(B,2)) + ud;

% use the PDE interpolant object to extrapolate to nodal degrees of freedom
[p,~,t] = self.Mesh.meshToPet();
pI = pdeInterpolant(p,t,sol,edge_elem);
if self.IsTwoD
    [Xvals,Yvals] = pI.extrapolateEdgeToNodes();
    u = [Xvals;Yvals];
else
    [Xvals,Yvals,Zvals] = pI.extrapolateEdgeToNodes();
    u = [Xvals;Yvals;Zvals];
end

end