function [bccode, nangle, boundn, boundt] = assignbcs(numnp,x,y,np,e,code,normbc,tangbc)
% function to take info from mesh, input file and assign the proper
% boundary conditions and angles to each nodal point.

% make sure the surface edges are all oriented clockwise around the
% boundary

bccode = zeros(1,numnp); % type of boundary condition
nangle = zeros(1,numnp); % angle btwn x-y and norm-tang systems
boundt = zeros(1,numnp); % tangential boundary conditions
boundn = zeros(1,numnp); % normal boundary conditions
e=FixE(e,np');

%
% Number of nodes on the boundary is same as the number of
% edges on the boundary, nbnodes = number of boundary nodes
nbnodes = size(e,2);

% bnodelist is the list of boundary nodes,
% and the global edge numbers to which the node belongs
% first row contains the node number
% second row contains the global edge number containing that node
% third row contains the other global edge which has the current node in it
bnodelist = zeros(3,nbnodes);
% n_1 is the nodal array of first edge normals
% n_2 is the nodal array of second edge normals
n_1 = zeros(2,nbnodes); n_2 = zeros(2,nbnodes);
e_1 = zeros(2,nbnodes); e_2 = zeros(2,nbnodes);
e1_len=zeros(1,nbnodes);% lengths of first edges
e2_len=zeros(1,nbnodes);% lengths of second edges
for i = 1:nbnodes
    bnodelist(1,i) = e(1,i); % global node number
    n_1(:,i) = findnormal(e(1,i),e(2,i),x,y); % store the normal to first edge
    e_1(:,i) = findtangent(e(1,i),e(2,i),x,y); % store the tangent to first edge
    e1_len(i)= findlength(e(1,i),e(2,i),x,y);
    bnodelist(2,i) = e(5,i); % global boundary number on one side
    otheredge=find(e(2,:)==e(1,i));
    bnodelist(3,i) = e(5,otheredge); % global boundary number on the other side
    n_2(:,i) = ...% store the normal to second edge
        findnormal(e(1,otheredge),e(2,otheredge),x,y);
    e_2(:,i) = ...% store the normal to second edge
        findtangent(e(1,otheredge),e(2,otheredge),x,y);
    e2_len(i)= findlength(e(1,otheredge),e(2,otheredge),x,y);
end
% Now we have the normals on both the edges and the boundary numbers
% Here we take care of corner bc's and
% find the normals which should be assigned to the node.
% the boundary conditions assigned to the node then have to be rotated to
% reflect the angle with which the local stiffness matrix will be changed
node_normal = zeros(2,nbnodes); % the normal to be assigned to the node
for i = 1:nbnodes % loop over the boundary nodes
    node =bnodelist(1,i); % global node number
    edge1=bnodelist(2,i); % boundary on one side of the node
    edge2=bnodelist(3,i); % boundary on the other side of the node
    % node in the middle of a boundary? choose those BC's
    if (edge1==edge2) % it's in the middle of a boundary, most frequent case
        bccode(node) = code(edge1); % set the bc code to that of one edge
        node_normal(:,i) = (n_1(:,i)+n_2(:,i))/2; % Choose mean of normals
        switch bccode(node)
            case 1
                boundn(node)=normbc(edge1)*e1_len(i)/2+normbc(edge2)*e2_len(i)/2;
                boundt(node)=tangbc(edge1)*e1_len(i)/2+tangbc(edge2)*e2_len(i)/2;
            case 2
                boundn(node)=normbc(edge1);
                boundt(node)=tangbc(edge1)*e1_len(i)/2+tangbc(edge2)*e2_len(i)/2;
            case 3
                boundn(node)=normbc(edge1)*e1_len(i)/2+normbc(edge2)*e2_len(i)/2;
                boundt(node)=tangbc(edge1);
            case 4
                boundn(node)=normbc(edge1);
                boundt(node)=tangbc(edge1);
        end % done with mid-boundary nodes
    else % edges on either side of node are on different global boundaries
        if (code(edge1)==4) % If edge-1 is fully constrained
            bccode(node)=4;
            node_normal(:,i) = n_1(:,i);
            boundn(node)=normbc(edge1);
            boundt(node)=tangbc(edge1);
        elseif (code(edge2)==4) % If edge-2 is fully constrained
            bccode(node)=4;
            node_normal(:,i) = n_2(:,i);
            boundn(node)=normbc(edge2);
            boundt(node)=tangbc(edge2);
        else % neither edge is fully constrained
            curang=dot(n_1(:,i),n_2(:,i));
            if abs(curang)<=1e-8
                curang=0;
            elseif (abs(curang)<(1+1e-8))&(abs(curang)>(1-1e-8))
                curang=1;
            end
            switch curang % make 3 cases based on how the edges are aligned
                case 1 % The two edges are parallel
                    node_normal(:,i) = n_1(:,i); % choose the normal of first edge
                    switch code(edge1) % always choose essential bcs over natural bcs
                        case 1
                            switch code(edge2)
                                case 1
                                    bccode(node) = 1;
                                    boundn(node) = normbc(edge1)*e1_len(i)/2+normbc(edge2)*e2_len(i)/2;
                                    boundt(node) = tangbc(edge1)*e1_len(i)/2+tangbc(edge2)*e2_len(i)/2;
                                case 2
                                    bccode(node) = 2;
                                    boundn(node) = normbc(edge2); % essential of edge2 wins
                                    boundt(node) = tangbc(edge1)*e1_len(i)/2+tangbc(edge2)*e2_len(i)/2;
                                case 3
                                    bccode(node) = 3;
                                    boundn(node) = normbc(edge1)*e1_len(i)/2+normbc(edge2)*e2_len(i)/2;
                                    boundt(node) = tangbc(edge2);
                            end %end of case e1=1
                        case 2
                            switch code(edge2)
                                case 1
                                    bccode(node) = 2;
                                    boundn(node) = normbc(edge1); % essential of edge-1 wins
                                    boundt(node) = tangbc(edge1)*e1_len(i)/2+tangbc(edge2)*e2_len(i)/2;
                                case 2
                                    bccode(node) = 2;
                                    boundn(node) = normbc(edge1);
                                    boundt(node) = tangbc(edge1)*e1_len(i)/2+tangbc(edge2)*e2_len(i)/2;
                                case 3
                                    bccode(node) = 4;
                                    boundn(node) = normbc(edge1); % essential of edge1
                                    boundt(node) = tangbc(edge2); % essential of edge2
                            end %end of case e1=2
                        case 3
                            switch code(edge2)
                                case 1
                                    bccode(node) = 3;
                                    boundn(node) = normbc(edge1)*e1_len(i)/2+normbc(edge2)*e2_len(i)/2;
                                    boundt(node) = tangbc(edge1);
                                case 2
                                    bccode(node) = 4;
                                    boundn(node) = normbc(edge2); % essential of edge2 wins
                                    boundt(node) = tangbc(edge1); % essential of edge1 wins
                                case 3
                                    bccode(node) = 3;
                                    boundn(node) = normbc(edge1)*e1_len(i)/2+normbc(edge2)*e2_len(i)/2;
                                    boundt(node) = tangbc(edge1);
                            end %end of case e1=3
                    end % end of all cases with parallel edges

                case 0 %two edges perpendicular
                    node_normal(:,i) = n_1(:,i); % choose the normal of first edge
                    switch code(edge1)
                        case 1
                            switch code(edge2)
                                case 1
                                    bccode(node) = 1;
                                    boundn(node) = normbc(edge1)*e1_len(i)/2+tangbc(edge2)*e2_len(i)/2;
                                    boundt(node) = tangbc(edge1)*e1_len(i)/2-normbc(edge2)*e2_len(i)/2;
                                case 2
                                    bccode(node) = 3;
                                    boundn(node) = normbc(edge1)*e1_len(i)/2+tangbc(edge2)*e2_len(i)/2;
                                    boundt(node) = -normbc(edge2);
                                case 3
                                    bccode(node) = 2;
                                    boundn(node) = tangbc(edge2);
                                    boundt(node) = tangbc(edge1)*e1_len(i)/2-normbc(edge2)*e2_len(i)/2;
                            end
                        case 2
                            switch code(edge2)
                                case 1
                                    bccode(node) = 2;
                                    boundn(node) = normbc(edge1);
                                    boundt(node) = tangbc(edge1)*e1_len(i)/2-normbc(edge2)*e2_len(i)/2;
                                case 2
                                    bccode(node) = 4;
                                    boundn(node) = normbc(edge1);
                                    boundt(node) = -normbc(edge2);
                                case 3
                                    bccode(node) = 2;
                                    boundn(node) = normbc(edge1);
                                    boundt(node) = tangbc(edge1)*e1_len(i)/2-normbc(edge2)*e2_len(i)/2;
                            end
                        case 3
                            switch code(edge2)
                                case 1
                                    bccode(node) = 3;
                                    boundn(node) = normbc(edge1)*e1_len(i)/2+tangbc(edge2)*e2_len(i)/2;
                                    boundt(node) = tangbc(edge1);
                                case 2
                                    bccode(node) = 3;
                                    boundn(node) = normbc(edge1)*e1_len(i)/2+tangbc(edge2)*e2_len(i)/2;
                                    boundt(node) = tangbc(edge1);% essential bc from other
                                    % edge same for continuity
                                case 3
                                    bccode(node) = 4;
                                    boundn(node) = tangbc(edge2); % essential bc from other
                                    % edge is also same
                                    boundt(node) = tangbc(edge1);% essential bc from other
                                    % edge same for continuity
                            end
                    end % end of all cases with perpendicular edges

                otherwise %two edges are at some arbitrary angle
                    switch code(edge1)
                        case 1
                            switch code(edge2)
                                case 1
                                    node_normal(:,i) = (n_1(:,i)+n_2(:,i))/2; % Choose mean of normals
                                    mynorm=node_normal(:,i);
                                    mytang=(e_1(:,i)+e_2(:,i))/2;
                                    bccode(node) = 1;
                                    % components of edge2 along edge1,
                                    % both normal and tangential are resolved along that direction
                                    boundn(node) = ...
                                        dot(e_2(:,i),mynorm)*tangbc(edge2)*e2_len(i)/2+...
                                        dot(n_2(:,i),mynorm)*normbc(edge2)*e2_len(i)/2+...
                                        dot(e_1(:,i),mynorm)*tangbc(edge1)*e1_len(i)/2+...
                                        dot(n_1(:,i),mynorm)*normbc(edge1)*e1_len(i)/2;
                                    boundt(node) = ...
                                        dot(e_2(:,i),mytang)*tangbc(edge2)*e2_len(i)/2+...
                                        dot(n_2(:,i),mytang)*normbc(edge2)*e2_len(i)/2+...
                                        dot(e_1(:,i),mytang)*tangbc(edge1)*e1_len(i)/2+...
                                        dot(n_1(:,i),mytang)*normbc(edge1)*e1_len(i)/2;
                                case 2
                                    node_normal(:,i)=n_2(:,i);
                                    mynorm=node_normal(:,i);
                                    mytang=e_2(:,i);
                                    bccode(node) = 2;
                                    boundn(node) = normbc(edge2);
                                    boundt(node) =  ...
                                        tangbc(edge2)*e2_len(i)/2+...
                                        dot(e_1(:,i),mytang)*tangbc(edge1)*e1_len(i)/2+...
                                        dot(n_1(:,i),mytang)*normbc(edge1)*e1_len(i)/2;
                                case 3
                                    node_normal(:,i)=n_2(:,i);
                                    mynorm=node_normal(:,i);
                                    mytang=e_2(:,i);
                                    bccode(node) = 3;
                                    boundn(node) =  ...
                                        normbc(edge2)*e2_len(i)/2+...
                                        dot(e_1(:,i),mynorm)*tangbc(edge1)*e1_len(i)/2+...
                                        dot(n_1(:,i),mynorm)*normbc(edge1)*e1_len(i)/2;
                                    boundt(node) = tangbc(edge2);
                            end % end of case e1=1
                        case 2
                            switch code(edge2)
                                case 1
                                    node_normal(:,i) = n_1(:,i);
                                    mynorm=node_normal(:,i);
                                    mytang=e_1(:,i);
                                    bccode(node) = 2;

                                    boundn(node) = normbc(edge1);
                                    boundt(node) = ...
                                        dot(e_2(:,i),mytang)*tangbc(edge2)*e2_len(i)/2+...
                                        dot(n_2(:,i),mytang)*normbc(edge2)*e2_len(i)/2+...
                                        tangbc(edge1);
                                case 2
                                    node_normal(:,i) = n_1(:,i);
                                    mynorm=node_normal(:,i);
                                    mytang=e_1(:,i);
                                    bccode(node) = 4;
                                    boundn(node) = normbc(edge1);
                                    boundt(node) = dot(n_2(:,i),mytang)*normbc(edge2);
                                case 3
                                    bccode(node) = 4;
                                    node_normal(:,i) = n_1(:,i);
                                    mynorm=node_normal(:,i);
                                    mytang=e_1(:,i);
                                    boundn(node) = normbc(edge1);
                                    boundt(node) = dot(e_2(:,i),mytang)*tangbc(edge2);
                            end  % end of case e1=2
                        case 3
                            switch code(edge2)
                                case 1
                                    node_normal(:,i) = n_1(:,i);
                                    mynorm=node_normal(:,i);
                                    mytang=e_1(:,i);
                                    bccode(node) = 3;
                                    boundn(node) = ...
                                        dot(e_2(:,i),mynorm)*tangbc(edge2)*e2_len(i)/2+...
                                        dot(n_2(:,i),mynorm)*normbc(edge2)*e2_len(i)/2+...
                                        normbc(edge1)*e1_len(i)/2;
                                    boundt(node) =tangbc(edge1);
                                case 2
                                    node_normal(:,i) = n_1(:,i); % Choose mean of normals
                                    mynorm=node_normal(:,i);
                                    mytang=e_1(:,i);
                                    bccode(node) = 4;
                                    boundn(node) =dot(n_2(:,i),mynorm)*normbc(edge2);
                                    boundt(node) = tangbc(edge1);
                                case 3
                                    node_normal(:,i) = n_1(:,i);
                                    mynorm=node_normal(:,i);
                                    mytang=e_1(:,i);
                                    bccode(node) = 4;
                                    % components of edge2 along edge1,
                                    % both normal and tangential are resolved along that direction
                                    boundn(node) = dot(e_2(:,i),mynorm)*tangbc(edge2);
                                    boundt(node) =tangbc(edge1);
                            end % end of case e1=3
                    end % end of all cases with edges at arbitrary angle
            end %end the main switch statement on dot product of normals
        end % ends the if (code(edge1)==4) % If edge-1 is fully constrained
    end % end of all the cases- ends the loop if (edge1==edge2)
    v=[node_normal(1,i), node_normal(2,i)]*[ 0 -1; 1 0];

    nangle(node)=atan2(v(2),v(1));

end % done with all boundary nodes

function length = findlength(i,j,x,y)
%global x y
% i, j are the nodes on the edge whose normal we should find
% call the first node as 1 and the second node as 2
vx=x(j)-x(i); vy=y(j)-y(i);
length=sqrt(vx^2+vy^2);


function normal = findnormal(i,j,x,y)
% i, j are the nodes on the edge whose normal we should find
% call the first node as 1 and the second node as 2
x1 = x(i); y1 = y(i); x2 = x(j); y2 = y(j);
% To find the normal, shift the origin to (x1,y1) and then rotate the vector
% connecting point 1 to point 2 by 90 degree, rotation matrix
normalp = [x2-x1 y2-y1]*[ 0 1; -1 0]; % the matrix here is the 90 degree rotation matrix
% need to normalize the above
% and write it as a row vector
% the length of the above vector is vect_len
vect_len = sqrt(normalp(1)^2+normalp(2)^2);
normal = normalp'/vect_len;


function tangentv = findtangent(i,j,x,y)
% i, j are the nodes on the edge whose tangent we should find
% call the first node as 1 and the second node as 2
x1 = x(i); y1 = y(i); x2 = x(j); y2 = y(j);
tangentp = [x2-x1 y2-y1];
% need to normalize the length
divlen = sqrt((x2-x1)^2+(y2-y1)^2);
tangentv = tangentp'/divlen;


function earray=FixE(e,con);
% takes the edge array and rearranges it so that it's all
% counter-clockwise
numedges=size(e,2);
[f, m]=MyMeshFaces(con);
eordered=f(:, m==1);
a_o_e=[eordered, e(1:2,:)];
bnums=[e(5,:)];
a_o_e_srt=sort(a_o_e);
[a_o_e_unq, i, j] = unique(a_o_e_srt', 'rows');
j1=j(1:numedges);
j2=j(numedges+1:2*numedges);
[trash, j2inv] = sort(j2);
earray=[eordered(2,:);
    eordered(1,:);
    e(3:7,j2inv)];

function [faces, mult] = MyMeshFaces(con)
%   INPUT:
%   con is d x n,
%       the mesh connectivity (tetrahedra or triangles)
%   OUTPUT:
%   faces is (d-1) x m, (integer)
%         the connectivity of the surface
%   mult is 1 x m,  (integer, optional)
%        the multiplicity of each face; this should
%        be either 1 (for a surface face) or 2 (for
%        an interior face)
%
d = size(con, 1); % dimension plus one
n = size(con, 2); % number of elements
dface = d - 1;    % dimension of face
nallf = d*n;      % total number of faces, including interior ones
%
%  efaces is the connectivity for faces within an element
%
efaces = [1 2 3 ; 2 3 1 ; 3 1 2];
efaces = efaces(1:dface, :);
%
f = [];
for i=1:d;
    f = [f    con(efaces(:, i), :)];
end
%
fsrt = sort(f); % fsrt = f(:, iorig)
[allfaces, ifsrt] = sortrows(fsrt'); allfaces = allfaces'; %
% allfaces=
% fsrt(:, ifsrt)
[u, f_of, a_of] = unique(allfaces', 'rows'); % u = allfaces(:, f_of)
utest = f(:, ifsrt(f_of));

faces = utest;
%
%  Multiplicity
nfaces = size(faces, 2);
wunz = ones(1, nallf);
%
mult = sparse(wunz, a_of, wunz, 1, nfaces);
mult = full(mult);   % note sparse accumulates entries
%
