function [A_el,F_el] = cal_Melm_2D_b(iel,gradP,U);
//
//    calcul du systeme elementaire pour l'Ex 2
//
   global CONNEC COOR NUMER ADRESS UC Nelt Nno Npt NI;
   global XG PG Ng;

   zero = 1*10^(-5)

//
   x  = COOR(CONNEC(iel,:),1); y = COOR(CONNEC(iel,:),2);

   u_el = U(ADRESS(iel,:));
   
   A_el =zeros(Nno,Nno); F_el = zeros(Nno,1);

   for ig=1:Ng    // double boucle pour l'integration numerique
    xg = XG(ig);
    for jg=1:Ng
       yg = XG(jg);
//
// calcul de la matrice jacobienne
//
          [F,DxiFP,DetaFP] = Q1(xg,yg);
          j11 = x(1:4)'*DxiFP; 
          j12 = x(1:4)'*DetaFP;   
          j21 = y(1:4)'*DxiFP;  
          j22 = y(1:4)'*DetaFP;
//
          DT  = [j11 j12; j21 j22];
          jac = det(DT);
          B = inv(DT)';
//
       [FP,DxiFP,DetaFP] = Q2(xg,yg);
        
        DFP = B*[DxiFP'; DetaFP'];
        gradU = DFP * u_el;
        gp = norm(gradU);

        mu = cal_mu(gp);

        if  gp >= zero then
              mup = cal_Dmu(gp)/gp;  
        else   
              mup = 0 ;
        end;

       for i=1:Nno
         DFPi = DFP(:,i);
         F_el(i) = F_el(i) + PG(ig)*PG(jg)*(gradP*FP(i) - ...
                                            mu * gradU' * DFPi)*jac;
          for j=1:Nno
            DFPj = DFP(:,j);
            A_el(i,j) = A_el(i,j) + PG(ig)*PG(jg)*(mu* DFPj' * DFPi + ...
                                                   mup*(gradU'*DFPi)*(gradU'*DFPj) ) *jac;
          end;
       end;
    end;
  end;
endfunction

