function [edges, edge_elem] = constructEdgeElementConnectivity(self)
% constructEdgeElementConnectivity - Helper function to generate edge
% connectivity. 

% Copyright 2021 The MathWorks, Inc.

% Curl-element is not supported for quadratic mesh.
if isequal(self.Mesh.GeometricOrder,'quadratic')
    msgID = 'pde:ElectromagneticModel:meshNotLinearTetFor3DMS';
    if contains(self.AnalysisType,'harmonic')
        msgID = 'pde:ElectromagneticModel:meshNotLinearTetHarmonic';
    end
    error(message(msgID));
end

[p,~,t] = self.Mesh.meshToPet();
dim = size(p,1);
numElem = size(t,2);
DT = triangulation(t(1:dim+1,:)',p');
edges = DT.edges;
numEdges = size(edges,1);
tT = t';
edges=sortrows(sort(edges,2)); % Perhaps unnecessary
tz=sort(tT(:,1:dim+1),2);
eid = zeros(numEdges,1);
tid = (1:numElem)';
% rows of bm  form [v1, v2, 0] or [v1,v2,triangle index]
if self.IsTwoD
    bm = sortrows([edges,eid; tz(:,[1 2]),tid; tz(:,[1 3]),tid; tz(:,[2 3]),tid]);
else
    bm = sortrows([edges,eid; tz(:,[1 2]),tid; tz(:,[1 3]),tid; tz(:,[2 3]),tid; tz(:,[1 4]),tid; tz(:,[2 4]),tid; tz(:,[3 4]),tid]);
end
bm2 = [bm(:,3),cumsum(bm(:,3)==0)]; % counting 0s indicates edge index
bm2(bm2(:,1)==0,:)=[];
bm2 = sortrows(bm2);
edge_elem = (reshape(bm2(:,2),dim*(dim+1)/2,[]))';
end