function [A,b,addDof]=PortBCMatrixAssembly2(A,b,Nodes,Elements,EdgesOfElements,...
    Domains,mur,PortMesh,PortSolver,Dof,addDof)
%description: process Port boundary condition and assembly the
%corresponding matrix  
%beam wave

NbrFaces=length(PortMesh.OldFaces);

%gauss quadrature points
if PortSolver.type==1
    [xft,yft,pft,Integral2DOrder]=Get2DGuassPoints(2);
else
    [xft,yft,pft,Integral2DOrder]=Get2DGuassPoints(4);
end

%Initialization of matrix assembly variables
allAddDof=sum(addDof);%Offset for matrix assembly
newAddDof=PortSolver.modeNum;

NbrNonZeros=NbrFaces*3;
NumG=0;
Gj=zeros(NbrNonZeros,1);
Gv=complex(zeros(NbrNonZeros,1));
P=complex(zeros(1,1));
if PortSolver.type==1
    NumS=0;
    Si=zeros(NbrNonZeros,1);
    Sv=complex(zeros(NbrNonZeros,1));
    b1=complex(zeros(1,1));
else
    NumT=0;
    Ti=zeros(NbrNonZeros,1);
    Tv=complex(zeros(NbrNonZeros,1));
end


%loop over all faces
for n=1:NbrFaces
    %port mesh data 
    numEle=PortMesh.FacesConn(n,1);
    numFace=PortMesh.FacesConn(n,2);
    if numFace==1
        FaceNode=[1 2 3]';
        FaceBF=[1 2 4]';
    elseif numFace==2
        FaceNode=[1 2 4]';
        FaceBF=[1 3 5]';
    elseif numFace==3
        FaceNode=[1 3 4]';
        FaceBF=[2 3 6]';
    elseif numFace==4
        FaceNode=[2 3 4]';
        FaceBF=[4 5 6]';
    end

    %material
    domain=Domains(numEle);
    mu=1/mur(domain);

    %mesh data of object
    %coordinate of object nodes
    x=Nodes(Elements(numEle,:),1);
    y=Nodes(Elements(numEle,:),2);
    z=Nodes(Elements(numEle,:),3);
    %length of object edges
    l=zeros(6,1);
    l(1)=sqrt((x(1)-x(2))*(x(1)-x(2))+(y(1)-y(2))*(y(1)-y(2))+(z(1)-z(2))*(z(1)-z(2)));
    l(2)=sqrt((x(1)-x(3))*(x(1)-x(3))+(y(1)-y(3))*(y(1)-y(3))+(z(1)-z(3))*(z(1)-z(3)));
    l(3)=sqrt((x(1)-x(4))*(x(1)-x(4))+(y(1)-y(4))*(y(1)-y(4))+(z(1)-z(4))*(z(1)-z(4)));
    l(4)=sqrt((x(2)-x(3))*(x(2)-x(3))+(y(2)-y(3))*(y(2)-y(3))+(z(2)-z(3))*(z(2)-z(3)));
    l(5)=sqrt((x(2)-x(4))*(x(2)-x(4))+(y(2)-y(4))*(y(2)-y(4))+(z(2)-z(4))*(z(2)-z(4)));
    l(6)=sqrt((x(3)-x(4))*(x(3)-x(4))+(y(3)-y(4))*(y(3)-y(4))+(z(3)-z(4))*(z(3)-z(4)));
    %Jac matrix of object
    Jac=zeros(3,3);
    Jac(1,1)=x(1)-x(4);Jac(1,2)=y(1)-y(4);Jac(1,3)=z(1)-z(4);
    Jac(2,1)=x(2)-x(4);Jac(2,2)=y(2)-y(4);Jac(2,3)=z(2)-z(4);
    Jac(3,1)=x(3)-x(4);Jac(3,2)=y(3)-y(4);Jac(3,3)=z(3)-z(4);
    InvJac=inv(Jac);


    %compite basis function of faces
    %coordinate of faces nodes
    xx=x(FaceNode(:));
    yy=y(FaceNode(:));
    zz=z(FaceNode(:));
    %Jac matrix of faces
    fJac=zeros(3,3);
    fJac(1,1)=-xx(1)+xx(2);fJac(1,2)=-yy(1)+yy(2);
    fJac(2,1)=-xx(1)+xx(3);fJac(2,2)=-yy(1)+yy(3);fJac(3,3)=1;
    InvfJac=inv(fJac);
    fJacS=zeros(3,3);
    fJacS(1,1)=InvfJac(2,2);fJacS(1,2)=-InvfJac(2,1);
    fJacS(2,1)=-InvfJac(1,2);fJacS(2,2)=InvfJac(1,1);fJacS(3,3)=1;
    fTJac=fJac'/det(fJac);
    fDetJac=abs(det(fJac));
    %basis function for computing port function
    et=zeros(3,3,Integral2DOrder);ez=zeros(3,3,Integral2DOrder);
    curlEt=zeros(3,3,Integral2DOrder);curlEz=zeros(3,3,Integral2DOrder);
    temp=zeros(3,1);
    for i=1:3 %the 1-th basis function 
        for k=1:Integral2DOrder %the k-th point
            [temp(1),temp(2)]=BF_Et(i,xft(k),yft(k));temp(3)=0;
            et(:,i,k)=InvfJac*temp*l(FaceBF(i));
            temp(3)=BF_curlEt(i,xft(k),yft(k));temp(1)=0;temp(2)=0;
            curlEt(:,i,k)=fTJac*temp*l(FaceBF(i));
            temp(3)=BF_Ez(i,xft(k),yft(k));temp(1)=0;temp(2)=0;
            ez(:,i,k)=temp;
            [temp(1),temp(2)]=BF_curlEz(i,xft(k),yft(k));temp(3)=0;
            curlEz(:,i,k)=fJacS*temp;
        end
    end
    %port function of faces
    fE=complex(zeros(3,Integral2DOrder));fcurlE=complex(zeros(3,Integral2DOrder));
    if PortSolver.type==1 %if port input
        fN0=complex(zeros(3,Integral2DOrder));fcurlN0=complex(zeros(3,Integral2DOrder));
        fN1=complex(zeros(3,Integral2DOrder));fcurlN1=complex(zeros(3,Integral2DOrder));
    else
        fN2=complex(zeros(3,Integral2DOrder));fcurlN2=complex(zeros(3,Integral2DOrder));
    end
    fW=complex(zeros(3,Integral2DOrder));%test function
    for k=1:Integral2DOrder  %the k-th point
        for i=1:3  %the 1-th basis function 
            MappingEz=PortMesh.NewFaces(n,i);
            MappingEt=PortMesh.EdgesOfElements((n-1)*3+i);
            fE(:,k)=fE(:,k)+et(:,i,k)*PortSolver.Et(MappingEt)+ez(:,i,k)*PortSolver.Ez(MappingEz);
            fcurlE(:,k)=fcurlE(:,k)+curlEt(:,i,k)*PortSolver.Et(MappingEt)+curlEz(:,i,k)*PortSolver.Ez(MappingEz);
        end
        if PortSolver.type==1 %if port input
            fN0(:,k)=fE(:,k)*PortSolver.powerCoef;
            fcurlN0(:,k)=(fcurlE(:,k)+PortSolver.gamma*[fE(2,k);-fE(1,k);0])*PortSolver.powerCoef;
            fN1(:,k)=fE(:,k);
            fcurlN1(:,k)=fcurlE(:,k)-PortSolver.gamma*[fE(2,k);-fE(1,k);0];
        else
            fN2(:,k)=fE(:,k);
            fcurlN2(:,k)=fcurlE(:,k)+PortSolver.gamma*[fE(2,k);-fE(1,k);0];
        end
        fW(:,k)=(fE(:,k));
    end

    %compute basis function of object
    %Jac matrix 
    fJac2=zeros(2,3);
    fJac2(1,1)=-xx(1)+xx(2);fJac2(1,2)=-yy(1)+yy(2);
    fJac2(2,1)=-xx(1)+xx(3);fJac2(2,2)=-yy(1)+yy(3);
    vJac=inv(Jac);
    %compute basis function
    Alphaj=zeros(3,3,Integral2DOrder);
    for i=1:3
        for k=1:Integral2DOrder
            temp=[xft(k),yft(k)]*fJac2+[xx(1) yy(1) zz(1)];%2D reference -> 3D phyics
            temp=temp*vJac-[x(4) y(4) z(4)]*vJac;%3D phyics -> 3D reference
            temp=BF_Edge(1,FaceBF(i),temp(1),temp(2),temp(3));%basis function in 3D reference
            Alphaj(:,i,k)=InvJac*temp*l(FaceBF(i));%basis function in 3D phyics
        end
    end

    %compute submatrix fW fN0 fN1 fN2 Alphaj
    Ge=complex(zeros(3,1));
    Pe=complex(zeros(1,1));
    if PortSolver.type==1 %if port input
        Se=complex(zeros(3,1));
        be=complex(zeros(3,1));
        b1e=complex(zeros(1,1));
        for i=1:3
            for k=1:Integral2DOrder
                Ge(i,1)=Ge(i,1)+pft(k)*fDetJac*sum(fW(:,k).*[Alphaj(1,i,k);Alphaj(2,i,k);0]);

                temp=mu*cross(PortMesh.normal,fcurlN1(:,k));
                Se(i,1)=Se(i,1)+pft(k)*fDetJac*sum(Alphaj(:,i,k).*temp);

                temp=mu*cross(PortMesh.normal,fcurlN0(:,k));
                be(i,1)=be(i,1)+pft(k)*fDetJac*sum(Alphaj(:,i,k).*temp);
            end
        end
        for k=1:Integral2DOrder
            Pe=Pe+pft(k)*fDetJac*sum(fW(:,k).*[fN1(1,k);fN1(2,k);0]);

            b1e=b1e+pft(k)*fDetJac*sum(fW(:,k).*[fN0(1,k);fN0(2,k);0]);
        end
    else
        Te=complex(zeros(3,1));
        for i=1:3
            for k=1:Integral2DOrder
                temp=[Alphaj(1,i,k);Alphaj(2,i,k);0]*exp(-PortSolver.gamma*zz(1));
                Ge(i,1)=Ge(i,1)+pft(k)*fDetJac*sum(fW(:,k).*temp);

                temp=mu*cross(PortMesh.normal,fcurlN2(:,k));
                temp2=Alphaj(:,i,k)*exp(PortSolver.gamma*zz(1));
                Te(i,1)=Te(i,1)+pft(k)*fDetJac*sum(temp2.*temp);
            end
        end
        for k=1:Integral2DOrder
            Pe=Pe+pft(k)*fDetJac*sum(fW(:,k).*[fN2(1,k);fN2(2,k);0]);
        end
    end

    %Matrix assembly
    if PortSolver.type==1 %if port input
        for i=1:3
            MappingIndexi=EdgesOfElements((numEle-1)*6+FaceBF(i));
            NumG=NumG+1;
            Gj(NumG)=MappingIndexi;
            Gv(NumG)=Ge(i,1);
            NumS=NumS+1;
            Si(NumS)=MappingIndexi;
            Sv(NumS)=Se(i,1);
            b(MappingIndexi)=b(MappingIndexi)+be(i,1);
        end
        P=P+Pe;
        b1=b1+b1e;
    else
        for i=1:3
            MappingIndexi=EdgesOfElements((numEle-1)*6+FaceBF(i));
            NumG=NumG+1;
            Gj(NumG)=MappingIndexi;
            Gv(NumG)=Ge(i,1);
            NumT=NumT+1;
            Ti(NumT)=MappingIndexi;
            Tv(NumT)=Te(i,1);
        end
        P=P+Pe;
    end

end %loop over all faces

%Matrix assembly
if PortSolver.type==1 %if port input   allAddDof+Dof
    G=sparse(1,Gj,Gv,1,Dof);
    S=sparse(Si,1,Sv,Dof,1);
    P=sparse(1,1,P);
    A=[A S;G -P];
    b=[-b;b1];
%     A=A+S*G/P;
%     b=S*b1/P-b;
else
    G=sparse(1,Gj,Gv,1,Dof);
    T=sparse(Ti,1,Tv,Dof,1);
    P=sparse(1,1,P);
    A=[A [T;0];[G 0] -P];
    b=[b;0];
%     A=A+T*G/P;
%     b=b;
end

addDof=[addDof;newAddDof];

end