[b_volumeforces, volumes]=rhs_CR_3D(elements3, coordinates, element2faces, faces)

function [b_volumeforces, volumes]=rhs_CR_3D(elements3, coordinates, element2faces, faces)

volumes=zeros(size(elements3,1),1);
for j = 1:size(elements3,1)
    I=elements3(j,:); 
    verticesTrans=coordinates(I,:)'; %new definition (transpose operation) for faster performance!!!!
    volume=det([1 1 1 1; verticesTrans])/6; 
    volumes(j)=abs(volume);
end

middle_points=face_middlepoint(coordinates,faces);

GaussPointsMatrix=[1/4 1/4 1/4 1/4];
GaussPoints1=zeros(size(elements3,1),3); 
%GaussPoints2=GaussPoints1;
%GaussPoints3=GaussPoints1; 

for j = 1:size(elements3,1)
    vertices=coordinates(elements3(j,:),:);
    GaussPoints=GaussPointsMatrix*vertices; 
%     vertices=middle_points(element2faces(j,:),:);
%     GaussPoints=GaussPointsMatrix*vertices
    GaussPoints1(j,:)=GaussPoints(1,:);
    %GaussPoints2(j,:)=GaussPoints(2,:);
    %GaussPoints3(j,:)=GaussPoints(3,:);
end
%f_GaussPoints=[f(GaussPoints1), f(GaussPoints2), f(GaussPoints3)];
f_GaussPoints=[f(GaussPoints1)];

% Assembly
%Z=[areas areas areas].*(f_GaussPoints*[-1 2 2; 2 -1 2; 2 2 -1]/9);
Z=[volumes volumes volumes volumes].*(f_GaussPoints*[1 1 1 1]/4);
b_volumeforces=sparse(element2faces,ones(size(element2faces,1),4),Z);