function [M1, M2, F, b, x, e_conn] =  assemb_co(param)

% The FEM equation for the temp. dist at time t_{n+1} satisfies
%   (M_1 + M_2) z^{n+1} - M_1 z^n + F + b = 0
% The steady state solution is the solution of
%     M_2 z^{ss}  + F + b = 0


%% Initialization
%---------------------------------------------------------------------
  x_1    =    0.000;    y_1 =   0.000;
  x_2    =  param.L;    y_2 = param.W;

  n_nodesx   =  param.nodesx;    % number of nodes in the x direction
  n_nodesy   =  param.nodesy;    % number of nodes in the y direction

  n_gauss          = 7;  % number of point used in Gaussian integration


%%  Geometry Module
%-----------------------------------------------------------------------
   %  Generate a mesh of the unit square with quadratic elements
    n_nodes    = n_nodesx*n_nodesy;
    n_elements = 2*(n_nodesx-1)*(n_nodesy-1)/4;

    x_nodes = linspace(x_1, x_2, n_nodesx);
    y_nodes = linspace(y_1, y_2, n_nodesy);
%
%  Generate node coordinates.
%
    x = zeros(n_nodes, 2);
    for i=1:n_nodesx
      x_i = x_nodes(i); k = i - n_nodesx;
      for j=1:n_nodesy
        k = k + n_nodesx;
        x(k,:) = [ x_i  y_nodes(j)];
%         k = i+(j-1)*n_nodesx;
%         x(k,1) = x_nodes(i);
%         x(k,2) = y_nodes(j);
      end
    end

    %  Generate element connectivity
    e_conn = zeros(n_elements, 6);
    ie     = 0;
    for j=1:2:n_nodesy-1
      for i=1:2:n_nodesx-1
        ie = ie + 1;
        k = i + (j-1)*n_nodesx;
        e_conn(ie,:) = [k, k+2+2*n_nodesx, k+2*n_nodesx, ...
                        k + 1 + n_nodesx, k+1+2*n_nodesx, k+n_nodesx ];
        ie = ie + 1;
        e_conn(ie,:) = [k, k+2, k+2+2*n_nodesx, ...
                        k+1, k+2+n_nodesx, k+1+n_nodesx ];
      end
    end
%
%  Set up arrays for elements with faces on left/right boundary.
%
   Omega_left  =            1:(n_nodesx-1):n_elements;
   Omega_right = (n_nodesx-1):(n_nodesx-1):n_elements;
%
%  With Neumann/Robins bc all nodes are unknown.
%
   n_equations = n_nodes;
%
%  Comment this line for PMODE, uncomment if for SPMD.
%
   spmd
%
%  Set up codistributed structure
%
%  Column pointers and such for codistributed arrays.
%
  Vc  = codistributed.colon(1, n_equations);
  lP  = getLocalPart(Vc);
  lP_1 = lP(1);
  lP_end = lP(end);
  dPM = getCodistributor(Vc.Partition);
  col_shft = [0 cumsum(dPM(1:end-1))];
%
% local sparse arrays
%
  M1_part = sparse(n_equations, dPM(labindex));  M2_part = M1_part;
  b_part  = sparse(dPM(labindex), 1); F_part = b_part;

  fid  = fopen(['out' int2str(labindex)], 'w');
  fprintf(fid, 'Lab %d : column range %d  %d \n', labindex, lP_1, lP_end);

  [n_elements, nel_dof] = size(e_conn);
  [r, s, w] = twod_gauss(n_gauss);
%
%  Build the finite element matrices - Begin loop over elements
%
  n_el_lab = 0;
  for n_el=1:n_elements
    nodes_local     = e_conn(n_el,:);% which nodes are in this element
    % subset of nodes/columns on this lab
    lab_nodes_local = my_extract( nodes_local, lP_1, lP_end);
% fprintf(fid, '\n \n n_el = %d \n', n_el);
% fprintf(fid, 'lab_nodes_local has %d rows \n', size(lab_nodes_local, 1));

    if ~isempty(lab_nodes_local) % continue the calculation for this elmnt
    n_el_lab = n_el_lab + 1;
%
%  coordinates, shape functions & gradients evaluated at the Gauss pts.
%
    x_local                   = x(nodes_local,:);
    [x_g, w_g, phi, p_x, p_y] = twod_shape(x_local, r, s, w);


%-------------------------------------------------------------------
%% Element contributions to the weak form of the equations
%  First we evaluate the surface integral (2D) terms
%-------------------------------------------------------------------
    M1_loc =  twod_bilinear(      1  , phi, phi, w_g);
    M2_loc =  twod_bilinear(param.k_x, p_x, p_x, w_g) + ...
              twod_bilinear(param.k_y, p_y, p_y, w_g);
    src_g  =  source(x_g,  param); % distributed source term
    F_loc  =  twod_f_int( src_g , phi, w_g );
    % pre-allocate boundary arrays
    b_loc  =   zeros(nel_dof,1);
%---------------------------------------------------------------------
%%  Add boundary terms
%---------------------------------------------------------------------
%
% Add boundary-integral terms for left and right bounary elements
% Left  boundary n_el = [1:2*(n_nodesx-1):n_elements
% Right boundary n_el = [2*(n_nodesx-1):2*(n_nodesx-1):n_elements
%
      if any(Omega_left == n_el);
%
%  Here we impose Robin boundary term. The node points on the boundary
%  are nodes_local([3 6 1]) and  are at x_local([3 6 1],:)
%
        [phi_g1, ip, x_g1,  w_g1] = boundary2(3, x_local, 5); % n_gauss=5
        alpha_g1 = steps(x_g1(:,2), param.yh, param.ep, param.alpha );
        phi_a    = diag(alpha_g1)*phi_g1;  % mult. each row by alpha(y_g)
        beta_g1  = steps(x_g1(:,2), param.yh, param.ep, param.beta );
%
% 1D integrals are from y = w to y = 0
%
        for ii=1:3
            for jj=1:3
               M2_loc(ip(ii),ip(jj)) = M2_loc(ip(ii),ip(jj)) + ...
                            oned_f_int( phi_g1(:, jj), phi_a(:,ii), w_g1 );
            end
            b_loc(ip(ii),1) = b_loc(ip(ii),1) + ...
                                  oned_f_int( beta_g1, phi_a(:,ii), w_g1 );
        end

      elseif any(Omega_right == n_el);
%
%  Here we impose boundary term for a specified flux
%  The node points on the boundary are nodes_local([2 3 5]) and
%  are at x_local([2 3 5],:)
%  We need \int f(y) \phi(L, y) dy along the right boundary.
%
        [phi_g1, ip, x_g1,  w_g1] = boundary2(2, x_local, 5); % n_gauss=5;
        for ii=1:3
             b_loc(ip(ii),1) = b_loc(ip(ii),1) + ...
                             oned_f_int( param.flux, phi_g1(:,ii), w_g1 );
        end
      end

%-----------------------------------------------------------------
%% Assemble contributions into the global system matrices
%-----------------------------------------------------------------
%
      for n_t = 1:nel_dof               % local DOF  - test fcn
         t_glb  = nodes_local(n_t);     % global DOF - test fcn
         for n_u = 1:size(lab_nodes_local, 1)
            n_locj = lab_nodes_local(n_u, 1);  % local DOF in current n_el
            n_glbj = lab_nodes_local(n_u, 2) ...
                          -col_shft(labindex); % global DOF
            M1_part(t_glb,n_glbj) = M1_part(t_glb,n_glbj) ...
                                                    + M1_loc(n_t,n_locj);
            M2_part(t_glb,n_glbj) = M2_part(t_glb,n_glbj) ...
                                          +  param.dt*M2_loc(n_t,n_locj);
         end

%
         if t_glb >= lP_1 && t_glb <= lP_end % why is this needed ?
            t_loc = t_glb - col_shft(labindex);
            b_part(t_loc,1)  = b_part(t_loc,1)  - param.dt*b_loc(n_t,1);
            F_part(t_loc,1)  = F_part(t_loc,1)  - param.dt*F_loc(n_t,1);
         end
      end % for  n_t
%     fprintf(fid, 'b_loc(1:3) = %8.3f %8.3f %8.3f \n', b_loc(1:3));
%     fprintf(fid, 'b_loc(4:6) = %8.3f %8.3f %8.3f \n', b_loc(4:6));


    end % if not empty

  end  % n_el

  fprintf(fid,'Fraction of elements evaluated on this lab  is %8.6f \n',...
                                                      n_el_lab/n_elements);
%
% Assemble the local parts in a codistributed format.
%
  M1 = codistributed(M1_part, codistributor('1d', 2));
  M2 = codistributed(M2_part, codistributor('1d', 2));
  b  = codistributed(b_part , codistributor('1d', 1));
  F  = codistributed(F_part , codistributor('1d', 1));

   fclose(fid);
   end % spmd : comment this line for pmode, uncomment for spmd

  return
end