function A = Matrix_Assemble(meshs,Matrix)

%% Assemble the matrix corresponding the mesh, Base1,Base2 and the
% coifficient matrix coeff : A = Matrix_Assemble(mesh,Matrix,R_mesh)
%
% INPUT:  Base1, Base2, the left and right hand side base functions, coord3
% and coord4 is the quadrature matrix on the triangle and quadrilateral,
% mesh is the corresponding the triangulation, coeff is the coefficient
% matrix
% mesh: this is the working mesh for the left and right fem (or basis function). 
% R_mesh: when the mesh for the right fem (base) is different (coarser) from 
% the working mesh (always for the left fem)
%
% OUTPUT: A: the assembled matrix
% Some assumptions: 
%           1. There are most two meshes which are nested. 
%           2. Can be used for the multilevel correction method and
%           nonlinear problems.
%           3. If Aux_FEM has a mesh which should be the fine mesh and will
%           be used as the working mesh and then the right and left meshes
%           are both coarse mesh. 
%           4. If Aux_FEM is a cell and there is a mesh in it, the mesh
%           should be sace in the first element of cell.
% HEHU XIE: 02-10-2010
fprintf('Assemble the matrix\n');
tic
%--------get the element matrix size---------------------------------------
if size(Matrix.Coord3,1)
    Nr_int = size(Matrix.Coord3,1);
else
    fprintf('There is no quadrature information\n');
    return
end
%--------compute the base value at the integral points---------------------
%--------check the meshes in the system------------------------------------
if size(meshs,2) == 0
    fprintf('There is no mesh inout!\n');
    return
end
N_mesh = size(meshs,2);
Nr_elem = size(meshs(1).elem,1);
% %----get the number of triangles-------------------------------------------
% II3 = find(meshs(1).elem(:,4)<0);
% Nr_elem3 = length(II3);
% II4 = setdiff([1:Nr_elem]',II3);
% Nr_elem4 = length(II4);
%--------compute the base value at the integral points---------------------
Nr_elem = size(meshs(1).elem,1);
%----------compute the freedom degree vector-------------------------------
%Matrix.Left_Base
if isstruct(Matrix.L_Base)
    [Free_deg3C,Dim_C] = DOF_Mangement(meshs(Matrix.L_Mesh),Matrix.L_Base.DOF_3);
    m_C = size(Matrix.L_Base.gb_3([0,0]),2);
    C_B = 1;  % indicate there is basis function for column
else 
    Dim_C = 1;
    m_C = 1;
    Free_deg3C = ones(Nr_elem,1);
    C_B = 0;   % indicate there is no bases function for row
end
if isstruct(Matrix.R_Base)
    [Free_deg3R,Dim_R] = DOF_Mangement(meshs(Matrix.R_Mesh),Matrix.R_Base.DOF_3);
    m_R = size(Matrix.R_Base.gb_3([0,0]),2);
    R_B = 1; % indicate there is basis function for row 
else
    Dim_R = 1;
    m_R = 1;
    Free_deg3R = ones(Nr_elem,1);
    R_B = 0;  % indicate there is no basis function for row
end
%-----check the auxiliary fem funtion for nonlinear problem----------------
Aux_FEM_Ind = 0;
if isstruct(Matrix.Aux_FEM) 
    Aux_FEM_Ind = 1;   % 1: means there exist Aux_FEMs
    Num_Aux_FEM = size(Matrix.Aux_FEM,2);
    Free_deg3A = {};
    for ind = 1 : Num_Aux_FEM
        [Free_deg3A{ind},Dim_A(ind)] = ...
            DOF_Mangement(meshs(Matrix.Aux_FEM(ind).Mesh_Pos),Matrix.Aux_FEM(ind).Base.DOF_3);
        m_A(ind) = size(Matrix.Aux_FEM(ind).Base.gb_3([0,0]),2); 
    end      
end
%--------------------------------------------------------------------------------
A = sparse(Dim_R,Dim_C);
coord3 = Matrix.Coord3;
%-------------produce the stiffness matrix---------------------------------
if Nr_elem > 0    
    [area,inv_transf] = Transform_information(meshs(1),1,1);
    if N_mesh > 1
        [area2,inv_transf2] = Transform_information(meshs(2),1,1);
        Node1 = meshs(1).elemlocalnode(:,2)*[1 0] + meshs(1).elemlocalnode(:,3)*[0 1];
        Node2 = meshs(1).elemlocalnode(:,5)*[1 0] + meshs(1).elemlocalnode(:,6)*[0 1];
        Node3 = meshs(1).elemlocalnode(:,8)*[1 0] + meshs(1).elemlocalnode(:,9)*[0 1];
    end  

    for kk = 1 : Nr_int
        %----------Compute the coefficient matrix at the integral point----
        X = (1-coord3(kk,1)-coord3(kk,2))*meshs(1).node(meshs(1).elem(:,1),1:2)+...
                           coord3(kk,1)*meshs(1).node(meshs(1).elem(:,2),1:2)+...
                           coord3(kk,2)*meshs(1).node(meshs(1).elem(:,3),1:2);
        if C_B == 1
            if Matrix.L_Mesh == 1
                if Matrix.L_Diff(1) == 1
                    L0 =  Matrix.L_Base.b_3([0 0]*(1-coord3(kk,1)-coord3(kk,2))+[1 0]*coord3(kk,1)+[0 1]*coord3(kk,2));
                end
                if sum(Matrix.L_Diff(2:3)) > 1
                    LB_D = Matrix.L_Base.gb_3([0 0]*(1-coord3(kk,1)-coord3(kk,2))+[1 0]*coord3(kk,1)+[0 1]*coord3(kk,2));
                end
            else
                 if Matrix.L_Diff(1) == 1
                    L0 = Matrix.L_Base.b_3(Node1*(1-coord3(kk,1)-coord3(kk,2))+Node2*coord3(kk,1)+Node3*coord3(kk,2));
                end
                if sum(Matrix.L_Diff(2:3)) > 1
                    [LB_1x,LB_1y] = Matrix.L_Base.gb_multinode_3(Node1*(1-coord3(kk,1)-coord3(kk,2))+Node2*coord3(kk,1)+Node3*coord3(kk,2));
                end                
            end
        end %end for C_B == 1
        if R_B == 1
            if Matrix.R_Mesh == 1
                if Matrix.R_Diff(1) == 1
                    R0 = Matrix.R_Base.b_3([0 0]*(1-coord3(kk,1)-coord3(kk,2))+[1 0]*coord3(kk,1)+[0 1]*coord3(kk,2)); 
                end
                if sum(Matrix.R_Diff(2:3)) > 1
                    RB_D = Matrix.R_Base.gb_3([0 0]*(1-coord3(kk,1)-coord3(kk,2))+[1 0]*coord3(kk,1)+[0 1]*coord3(kk,2)); 
                end
            else
                if Matrix.R_Diff(1) == 1
                    R0 = Matrix.R_Base.b_3(Node1*(1-coord3(kk,1)-coord3(kk,2))+Node2*coord3(kk,1)+Node3*coord3(kk,2));
                end
                if sum(Matrix.R_Diff(2:3)) > 1
                    %-----produce two vectors----------------
                    [RB_1x,RB_1y] = Matrix.R_Base.gb_multinode_3((1-coord3(kk,1)-coord3(kk,2))*Node1+coord3(kk,1)*Node2+coord3(kk,2)*Node3);
                end
            end % end if Matrix.R_Mesh == 1
        end  % end if R_B == 1                  
        Aux_FEM = [];
        %------for multi aux fem function        
        if Aux_FEM_Ind == 1    
                 %   Num_Aux_FEM
        
            for ind = 1 : Num_Aux_FEM                
                dim_aux = size(Matrix.Aux_FEM(ind).Vec,2);    
                %dim_aux
                if Matrix.Aux_FEM(ind).Diff(1) == 1
                    Aux0 =  Matrix.Aux_FEM(ind).Base.b_3([0 0]*(1-coord3(kk,1)-coord3(kk,2))+[1 0]*coord3(kk,1)+[0 1]*coord3(kk,2));               
                end
                if sum(Matrix.Aux_FEM(ind).Diff(2:3)) > 1
                    Aux_D = Matrix.Aux_FEM(ind).Base.gb_3([0 0]*(1-coord3(kk,1)-coord3(kk,2))+[1 0]*coord3(kk,1)+[0 1]*coord3(kk,2));
                end
                for aux_ind = 1 : dim_aux
                    vec_tmp = Matrix.Aux_FEM(ind).Vec(:,aux_ind);
                    if Matrix.Aux_FEM(ind).Diff(1) == 1
                        % Aux_L0
                        tmp = vec_tmp(Free_deg3A{ind})*Aux0';
                        if Matrix.Aux_FEM(ind).Mesh_Pos == 1
                            Aux_FEM = [Aux_FEM,tmp];
                        else
                            Aux_FEM = [Aux_FEM,tmp(meshs(1).owner,1)];
                        end
                    end % end if
                    if Matrix.Aux_FEM(ind).Diff(2) == 1
                        %VLX
                        tmp = vec_tmp(Free_deg3A{cell_id})*(inv_transf(:,1,1)*Aux_D(1,:)+inv_transf(:,2,1)*Aux_D(2,:))';
                        if Matrix.Aux_FEM(ind).Mesh_Pos == 1
                            Aux_FEM = [Aux_FEM,tmp];
                        else
                            Aux_FEM = [Aux_FEM,tmp(meshs(1).owner,1)];
                        end                        
                    end
                    if Matrix.Aux_FEM(ind).Diff(3) == 1
                        % VLY
                        tmp = vec_tmp(Free_deg3A{cell_id})*(inv_transf(:,1,2)*Aux_D(1,:)+inv_transf(:,2,2)*Aux_D(2,:))'; 
                         if Matrix.Aux_FEM(ind).Mesh_Pos == 1
                            Aux_FEM = [Aux_FEM,tmp];
                        else
                            Aux_FEM = [Aux_FEM,tmp(meshs(1).owner,1)];
                         end  % if Matrix.Aux_FEM(ind).Mesh_Pos == 1
                    end %end if Matrix.Aux_FEM(ind).Diff(3) == 1
                end % end for aux_ind = 1 : dim_aux
            end %end for ind = 1: Num_Aux_FEM
        end % end if Aux_FEM_Ind == 1        
        %: Now come to produce the matrix values    
        %: iteration for column indices 
        for i = 1 : m_C            
            VL = [];
            if C_B == 1
                if Matrix.L_Diff(1) == 1
                    if Matrix.L_Mesh == 1
                        % VL0
                        VL = [VL,L0(i)*ones(Nr_elem,1)];
                    else
                        VL = [VL,L0(:,i)];
                    end % end if Matrix.L_Mesh == 1
                end  % end  if Matrix.L_Diff(1) == 1
                if Matrix.L_Diff(2) == 1
                    if Matrix.L_Mesh == 1
                        %VLX
                        VL = [VL,inv_transf(:,1,1)*LB_D(1,i)+inv_transf(:,2,1)*LB_D(2,i)];
                    else
                        VL = [VL,LB_1x(:,i).*inv_transf2(meshs(1).owner,1,1)+LB_1y(:,i).*inv_transf2(meshs(1).owner,2,1)];
                    end  % end if Matrix.L_Mesh == 1
                end % end if Matrix.L_Diff(2) == 1
                if Matrix.L_Diff(3) == 1
                    if Matrix.L_Mesh == 1
                        % VLY
                        VL = [VL,inv_transf(:,1,2)*LB_D(1,i)+inv_transf(:,2,2)*LB_D(2,i)];
                    else
                        VL = [VL,LB_1x(:,i).*inv_transf2(meshs(1).owner,1,2)+LB_1y(:,i).*inv_transf2(meshs(1).owner,2,2)];                       
                    end  % end if Matrix.L_Mesh == 1
                end % end if Matrix.L_Diff(3) == 1
            end  %end for if C_B == 1
            %----iteration for the row indices------
            for j = 1 : m_R
                VR = [];
                if R_B == 1
                    %-------The right hand side term ------------------------
                    if Matrix.R_Diff(1) == 1
                        if Matrix.R_Mesh == 1
                            %VR0
                            VR = [VR,R0(j)*ones(Nr_elem,1)];                            
                        else
                            %VR0
                            VR = [VR,R0(:,j)];
                        end %end if Matrix.R_Mesh == 1
                    end  % end  if Matrix.R_Diff(1) == 1
                    if Matrix.R_Diff(2) == 1
                        if Matrix.R_Mesh == 1
                            %VRX
                            VR = [VR,inv_transf(:,1,1)*RB_D(1,j)+inv_transf(:,2,1)*RB_D(2,j)];                            
                        else
                            %VRX
                            VR = [VR,RB_1x(:,j).*inv_transf2(meshs(1).owner,1,1)+RB_1y(:,j).*inv_transf2(meshs(1).owner,2,1)];
                        end % end if Matrix.R_Mesh == 1
                    end % end if Matrix.R_Diff(2) == 1
                    if Matrix.R_Diff(3) == 1
                        if Matrix.R_Mesh == 1
                            %VRY
                            VR = [VR,inv_transf(:,1,2)*RB_D(1,j)+inv_transf(:,2,2)*RB_D(2,j)];                            
                        else
                            %VRY
                            VR = [VR,RB_1x(:,j).*inv_transf2(meshs(1).owner,1,2)+RB_1y(:,j).*inv_transf2(meshs(1).owner,2,2)];
                        end %end if Matrix.R_Mesh == 1
                    end  % end if Matrix.R_Diff(3) == 1
                end  %end for R_B==1
                value = Matrix.Discrete_Form(VL,VR,X,Aux_FEM,Matrix);
                if Matrix.L_Mesh == 2
                    if Matrix.R_Mesh == 2
                        A = A + sparse(Free_deg3R(meshs(1).owner,j),Free_deg3C(meshs(1).owner,i),coord3(kk,3)*value.*area, Dim_R,Dim_C);                    
                    else
                       A = A + sparse(Free_deg3R(II3,j),Free_deg3C(meshs(1).owner,i),coord3(kk,3)*value.*area, Dim_R,Dim_C);                           
                    end % end if Matrix.R_Mesh == 2
                else
                    if Matrix.R_Mesh == 2
                        A = A + sparse(Free_deg3R(meshs(1).owner,j),Free_deg3C(:,i),coord3(kk,3)*value.*area, Dim_R,Dim_C);                    
                    else
                       A = A + sparse(Free_deg3R(:,j),Free_deg3C(:,i),coord3(kk,3)*value.*area, Dim_R,Dim_C);                           
                    end % end if Matrix.R_Mesh == 2
                end %end if Matrix.L_Mesh == 2
            end %end for j = 1:M_R
        end %end for i = 1:M_C     
    end %end for kk
end %end for Nr_elem
toc