function [H, R, faceWithBCs] = assembleDirichletBCNedelec(self,edges)

% Copyright 2020-2022 The MathWorks, Inc.

if self.ParentModel.IsTwoD
    if isempty(self.ParentModel.EdgeBC)
        error(message('pde:ElectromagneticModel:noBCs'))
    end
else
    if isempty(self.ParentModel.FaceBC)
        error(message('pde:ElectromagneticModel:noBCs'))
    end
end

numEdge = size(edges,1);
Hmap = containers.Map('KeyType','uint32','ValueType','any');
if contains(self.ParentModel.AnalysisType,'magnetostatic')
    BCType = 'MagneticPotential';
elseif contains(self.ParentModel.AnalysisType,'electricHarmonic') 
    BCType = 'ElectricField';
else % magneticHarmonic
    BCType = 'MagneticField';
end

isAxisymm = (self.ParentModel.Geometry.Dim  ==  2)...
    && (self.ParentModel.PlanarType == "axisymmetric");

if self.ParentModel.IsTwoD 
    numFacets = self.ParentModel.Geometry.NumEdges;
    faceWithBCs = false(numFacets,1);
    R = zeros(numEdge,1);
    for boundaryFacet = 1:numFacets
        bc = self.ParentModel.EdgeBC(boundaryFacet);
        if ~isempty(bc) && ~isempty(bc.FarField)
            % set far field BC to zeros here
            bc.ElectricField = [0;0];
            bc.MagneticField = [0;0];
        end
        if isempty(bc) || isempty(bc.(BCType))
            continue;
        end
        if isAxisymm
            bc.(BCType) = axisymmetricModification(bc.(BCType));
        end
        faceWithBCs(boundaryFacet) = true;
        tt = self.Mesh.boundaryFacets('Edge',boundaryFacet);
        for iEdge = 1:size(tt,2)
            dofInd = getBoundaryEdgeDof(iEdge,tt,edges);
            nodeOneInd = self.Mesh.Nodes(:,edges(dofInd,1)); 
            nodeTwoInd = self.Mesh.Nodes(:,edges(dofInd,2));
            if isa(bc.(BCType),'function_handle')
                bcVal =  computeNedDOF2D(bc.(BCType),nodeOneInd,nodeTwoInd);
            else
                dr = nodeTwoInd - nodeOneInd;
                bcVal = dot(bc.(BCType),dr);
            end
            if ~Hmap.isKey(dofInd)
                Hmap(dofInd) = 1;
                R(dofInd) = bcVal;
            end
        end
    end
    mH = size(Hmap,1);
    nH = size(edges,1);
    dofInd = double(cell2mat(Hmap.keys));
    H = sparse(1:mH,dofInd,ones(mH,1),mH,nH);
    R = R(dofInd);
    emptyDirichletBCError(self.ParentModel.AnalysisType,R)
    return
end

numFacets = self.ParentModel.Geometry.NumFaces;
faceWithBCs = false(numFacets,1);
R = zeros(numEdge,1);
% compute h and r 
for boundaryFacet = 1:numFacets
    bc = self.ParentModel.FaceBC(boundaryFacet);
    if ~isempty(bc) && ~isempty(bc.FarField)
        % set far field BC to zeros here
        bc.ElectricField = [0;0;0];
        bc.MagneticField = [0;0;0];
    end
    if isempty(bc) || isempty(bc.(BCType))
        continue;
    end
    faceWithBCs(boundaryFacet) = true;
    tt = self.Mesh.boundaryFacets('Face',boundaryFacet);
    boundaryEdges = [tt(1,:),tt(2,:),tt(3,:) ; tt(2,:),tt(3,:),tt(1,:)];
    boundaryEdges = unique(sort(boundaryEdges)','rows');
    boundaryEdges = boundaryEdges';
    for iEdge = 1:size(boundaryEdges,2)
        dofInd = getBoundaryEdgeDof(iEdge,boundaryEdges,edges);
        %Compute geometric edge vector.
        nodeOneInd = self.Mesh.Nodes(:,edges(dofInd,1)); 
        nodeTwoInd = self.Mesh.Nodes(:,edges(dofInd,2));
        
        if isa(bc.(BCType),'function_handle')
            bcVal =  computeNedDOF(bc.(BCType),nodeOneInd,nodeTwoInd);
        else
            dr = nodeTwoInd - nodeOneInd;
            bcVal = dot(bc.(BCType),dr);
        end
        
        if ~Hmap.isKey(dofInd)
            Hmap(dofInd) = 1;
            R(dofInd) = bcVal;
        end
    end
end
mH = size(Hmap,1);
nH = size(edges,1);
dofInd = double(cell2mat(Hmap.keys));
H = sparse(1:mH,dofInd,ones(mH,1),mH,nH);
R = R(dofInd);
emptyDirichletBCError(self.ParentModel.AnalysisType,R)
end

function dof = getBoundaryEdgeDof(edgeInd,e,edges)
% auxiliary function needed to compute boundary conditions
inds = edges == sort([e(1,edgeInd),e(2,edgeInd)]);
%inds = edges == [e(1,edgeInd),e(2,edgeInd)];
inds = inds(:,1)&inds(:,2);
dof = find(inds);
end

function alpha = computeNedDOF2D(boundaryFunc,x1,x2)
    dr = x2 - x1;
    xm = (x1+x2)./2;
    region.x = xm(1);
    region.y = xm(2);
    state.u = []';
    state.uy = [];
    state.ux = [];

    alpha = dot(boundaryFunc(region,state),dr)/norm(dr);
end

function alpha = computeNedDOF(boundaryFunc,x1,x2)
% auxiliary function needed to compute boundary conditions
% map a function to a Nedelec dof alpha
% x1 and x2 are the endpoints of the edge in global index order to
% maintain correct dof orientation
    dr = x2 - x1;
    xm = (x1+x2)./2;    
    region.x = xm(1);
    region.y = xm(2);
    region.z = xm(3);
    state.u = []';
    state.uy = [];
    state.ux = [];
    
    alpha = dot(boundaryFunc(region,state),dr)/norm(dr);
end

function modcoef = axisymmetricModification(bcIn)
    if isa(bcIn,'function_handle')
        modcoef = @axisymmetricFcn;
    else
        modcoef = @(region,state) 2*pi*region.x.*bcIn;
    end
    function q = axisymmetricFcn(region,state)
        [regionAxisymm, stateAxisymm] = pde.EquationModel.getAxisymmStructs(region,state);
        q = 2*pi*region.x.*bcIn(regionAxisymm,stateAxisymm);
    end
end

function emptyDirichletBCError(AnalysisType,R)
    if ~isempty(R)
        return
    end
    
    if contains(AnalysisType,'magnetostatic')
        error(message('pde:ElectromagneticModel:bcMissingFromMS'));
    end
    
    if contains(AnalysisType,'electricHarmonic')
        error(message('pde:ElectromagneticModel:bcMissingElectricHarmonic'));
    end
    
    if contains(AnalysisType,'magneticHarmonic') 
        error(message('pde:ElectromagneticModel:bcMissingMagneticHarmonic'));
    end
end
