function [Ke,Fe] = GaussintKF(h,D)
Ke = zeros(24,24);
Fe = zeros(24,6);
x0 = 0; y0 = 0; z0=0;
numgp=2;
gpx = getgp(x0,x0+h,numgp);
gpy = getgp(y0,y0+h,numgp);
gpz = getgp(z0,z0+h,numgp);
for ix = 1:numgp*2
    x=gpx(ix);
    for iy = 1:numgp*2 
        y=gpy(iy);
        for iz = 1:numgp*2
            z=gpz(iz);
            B = LocalB(x,y,z,x0,y0,z0,h);
            Ke = Ke + B.'*D*B;
            Fe = Fe + B.'*D;
        end
    end
end
Ke = Ke.*h^3/(numgp*2)^3;
Fe = Fe.*h^3/(numgp*2)^3;
end