treat_Robin_2d.m

function [w,R] = treat_Robin_2d(cfun,rfun,qfun,Gauss_coefficient_reference_1D,Gauss_point_reference_1D,boundary_edges,local_basis_trial,local_basis_test,P,T,Tb,vector_size,matrix_size,basis_type_trial,basis_type_test,der_order_trial_x,der_order_trial_y,der_order_test_x,der_order_test_y)

nbe = size(boundary_edges,2);

w = sparse(vector_size,1);

R = sparse(matrix_size,matrix_size);

for k = 1:nbe
    
    if boundary_edges(1,k) == -3
              
        nk = boundary_edges(2,k);
        
        end_point_1 = P(:,boundary_edges(3,k));
        
        end_point_2 = P(:,boundary_edges(4,k));
        
        [Gauss_nodes_2d_line,Gauss_weights_2d_line] = generate_2d_line_Gauss(Gauss_coefficient_reference_1D,Gauss_point_reference_1D,end_point_1,end_point_2);
        
        vertices = P(:,T(:,nk));
        
        for beta = 1:local_basis_test
            
            temp = Gauss_quadrature_2D_line_test(cfun,qfun,Gauss_weights_2d_line,Gauss_nodes_2d_line,vertices,basis_type_test,der_order_test_x,der_order_test_y,beta);
                
            w(Tb(beta,nk),1) = w(Tb(beta,nk),1) + temp; 
            
        end
        
        for alpha = 1:local_basis_trial
            
            for beta = 1:local_basis_test
            
                temp = Gauss_quadrature_2D_line_trial_test(cfun,rfun,Gauss_weights_2d_line,Gauss_nodes_2d_line,vertices,basis_type_trial,der_order_trial_x,der_order_trial_y,alpha,vertices,basis_type_test,der_order_test_x,der_order_test_y,beta);
                
                R(Tb(beta,nk),Tb(alpha,nk)) = R(Tb(beta,nk),Tb(alpha,nk)) + temp; 
            
            end
            
        end
        
    end
    
end