function Res = solveNedelecMS3D(self,coefstruct)
% Internal Only

%   Copyright 2021-2022 The MathWorks, Inc.

[edges, edge_elem] = constructEdgeElementConnectivity(self);
numEdges = size(edges,1);
[M,K,F,H,R,facesWithBC] = discretizedModelNedelec(self,coefstruct,edges,edge_elem,numEdges,[]);
[B,Or] = pdenullorth(H);
ud = Or*((H*Or)\R);
Kc = B'*K*B;
Fc = B'*(F - K*ud);

% Form saddle point problem, with penalty DoFs that enforce divergence-free constraint on vector magnetic potential.
numNodalDOF = size(self.Mesh.Nodes,2); 
faceIDs = 1:self.Geometry.NumFaces;
nBC = findNodes(self.Mesh,'region','Face',faceIDs(facesWithBC));
iS = 1:numel(nBC);
jS = nBC;
v = ones(size(nBC));
S = sparse(iS,jS,v,numel(nBC), numNodalDOF);
[C,~] = pdenullorth(S);
Mc = B'*M*C; 
Zc = sparse(size(C,2),size(C,2));

% construct matrix
A = [Kc, Mc; ...
     Mc', Zc];
b = [Fc ; zeros(size(C,2),1)];
% solve system
sol = A\b;
% Select edge DoFs
sol = B*sol(1:size(B,2)) + ud;

% Postprocess to get nodal values from edge solution.
[p,~,t] = self.Mesh.meshToPet();
pI = pdeInterpolant(p,t,sol,edge_elem);
[Xvals,Yvals,Zvals] = pI.extrapolateEdgeToNodes();

Res = [Xvals;Yvals;Zvals];
end

