function [M,K,F,H,R,faceWithBCs] = discretizedModelNedelec(self,coefstruct,edges,edge_elem,numEdges,A)
% Internal Only - should only be called within associated solve methods
% set up matrices 
% A is vector potential

% Copyright 2020-2022 The MathWorks, Inc.


[p,~,t] = self.Mesh.meshToPet();
numNodes = size(p,2);
for i = 1:numel(coefstruct.Coefficients.m)
    coefstruct.Coefficients.m{i} = 1;
end
% temporarily change system size to 3 to accomodate nonlinear c
self.PDESystemSize = 3;
femat = formGlobalMatrix(self,p,t,coefstruct,A,[],'kfm',edge_elem,numEdges);
self.PDESystemSize = 1;
K = femat.K;
F = femat.F;
% add magnetization to source if present
if startsWith(self.AnalysisType,'magnetostatic','IgnoreCase',true) && ...
        isfield(coefstruct.Coefficients,'fc')
    curlF = addCurlComponent(self,p,t,coefstruct,A,edge_elem,numEdges);
    F = F + curlF;
end

% Construct divergence constraint matrix.
[G,i1,i2] = formGlobalGradients3DCurl(edges,p,numEdges,numNodes);
G = sparse(i1,i2,G);
% 
M = femat.M*G;
[H,R,faceWithBCs] = assembleDirichletBCNedelec(self,edges);
end


