function Fe=FEMDbar(Fe)
% push forward the existing DLagr and add the initial stress component
% I x S.

    
    Fe.Dbar = zeros(Fe.ngauss,36);
    
    Findx = [1 4 7;     % addressing F
	     2 5 8;
	     3 6 9];
    
    % partial pushforward
    for kk=1:36,
      i=Fe.dirbar(kk,1);
      j=Fe.dirbar(kk,2);
      K=Fe.dirbar(kk,3);
      L=Fe.dirbar(kk,4);
      for I=1:3,
	for J=1:3;
	  fi = Findx(i,I);
	  fj = Findx(j,J);
	  ci = Fe.cindxbar(I,J,K,L);
	  Fe.Dbar(:,kk) = Fe.Dbar(:,kk) + Fe.F(:,fi).*Fe.F(:,fj).*Fe.DLagr(:,ci);
	end;
      end;
    end;

    Sindx = [1 2 3;
	     2 4 5;
	     3 5 6];
    
    for kk=1:36,
      I=Fe.dirbar(kk,1);
      J=Fe.dirbar(kk,2);
      K=Fe.dirbar(kk,3);
      L=Fe.dirbar(kk,4);
      Fe.Dbar(:,kk) = Fe.Dbar(:,kk) + (1/2)*(Fe.S(:,Sindx(K,L))*(I==J) + Fe.S(:,Sindx(I,J))*(K==L));
    end;
   
return;