function self = generateBoundaryFacetData(self,interfaceConditionsPermitted, varargin)
%generateBoundaryFacetData Convert geometry-based BC definitions to mesh facets
% Mapping of geometry-based BC definitions to mesh faces is a one-time
% operation per solution call, as element facets are not permitted to change
% during solution.
% This undocumented function may be changed or removed in a future release.

%       Copyright 2020-2021 The MathWorks, Inc.

meshorder = 'linear';
msh = [];
if ~isempty(varargin)
    self.points = varargin{1};
    e = varargin{2};
    if isa(e,'pde.FEMesh') || isa(e,'pde.FEMeshAssociation')
        msh = e;
    end
else
    self.points = self.myPde.Mesh.Nodes;
    msh = self.myPde.Mesh;
    meshorder = msh.GeometricOrder;
end
self.dim = size(self.points,1);
self.numPoints = size(self.points,2);

numbcs =length(self.faceBoundaryConditions);
for i=numbcs:-1:1
bcsi = self.faceBoundaryConditions(i);
    if strcmp(bcsi.appRegionID,'default')
        continue
    end
    
    bcsi.elemFaces = [];
    
    if self.dim == 2
        if isa(msh,'pde.FEMesh') && interfaceConditionsPermitted
            % Interface condition are allowed for thermal-stress
            % from R2018b. Uses the new mesh query to get the
            % discrete edges.
            if strcmp(bcsi.appRegionType,'Edge')
                bcsi.elemFaces = msh.boundaryFacets('Edge',bcsi.appRegionID);
                bcsi.nodes = self.getNodesForFacets(bcsi.elemFaces);
            else
                bcsi.nodes = msh.findNodes('region','Vertex',bcsi.appRegionID);
            end
        else
            % Interface conditions are not allowed in pre-2018b models. Uses
            % the e-matrix to get the discrete edges.
            
            % Get rid of unwanted edges, that are part of internal
            % boundaries between subdomians.
            if isa(msh,'pde.FEMesh')
                [~,e,~] = msh.meshToPet;
            end
            ie=pdesde(e);
            e=e(:,ie);
            
            if strcmp(meshorder,'quadratic')
                bcsi.elemFaces = extractDiscreteEdgeQuadratic(e, bcsi.appRegionID);
            else
                bcsi.elemFaces = extractDiscreteEdgeLinear(e, bcsi.appRegionID);
            end
            bcsi.nodes = self.getNodesForFacets(bcsi.elemFaces);
        end
    else
        if isa(msh,'pde.FEMesh')
            if strcmp(bcsi.appRegionType,'Edge')
                bcsi.nodes = msh.findNodes('region','Edge',bcsi.appRegionID);
            elseif strcmp(bcsi.appRegionType,'Vertex')
                bcsi.nodes = msh.findNodes('region','Vertex',bcsi.appRegionID);
            else
                bcsi.elemFaces = msh.boundaryFacets('Face',bcsi.appRegionID);
                if isempty(bcsi.elemFaces)
                    error(message('pde:pde3DBCImpl:noMeshOnFace',bcsi.appRegionID))
                end
                bcsi.nodes = self.getNodesForFacets(bcsi.elemFaces);
            end
        else
            % msh is a FEMeshAssociation object.
            bcsi.elemFaces = msh.getElementFaces(bcsi.appRegionID);
            bcsi.nodes = self.getNodesForFacets(bcsi.elemFaces);
        end
    end

     self.faceBoundaryConditions(i) = bcsi;
end
end

function de = extractDiscreteEdgeQuadratic(e,eid)
% Given the e from a p, e, t mesh, extract the discrete edge
% corresponding to the geometry edge whos ID is eid.

% The e matrix in 2-D treats two segments of a quadratic
% element as two separate edges. The following code stitches
% adjacent edges to obtain a quadratic edge. Domain to the left
% or right of edge info, available in e matrix, is used to
% preserve orientation of quadratic edges.

idx = (e(5,:) == eid);
subDomainToLeft = e(6, idx);

facets = e([1, 2], idx);
numfacets = size(facets,2);
oddidx = 1:2:numfacets;
evenidx = 2:2:numfacets;
de = [facets(:,oddidx);facets(2,evenidx)];
sizeDe = size(de);
if (~all(subDomainToLeft))
    de = flipud(de(:));
end
de = reshape(de,sizeDe);
end

function de = extractDiscreteEdgeLinear(e,eid)
% Given the e from a p, e, t mesh, extract the discrete edge
% corresponding to the geometry edge whos ID is eid.

%Unlike quadratic elements, no further manipulation is needed
%in linear element case.
idx = (e(5,:) == eid);
de = e([1, 2], idx);
end