[M_CR Z]=stima_gradgradCR_3D(elements3,coordinates,element2faces,faces,coefficient_difussion,Z)

function [M_CR Z]=stima_gradgradCR_3D(elements3,coordinates,element2faces,faces,coefficient_difussion,Z)

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

if (nargin==4)
    coefficient_difussion=ones(size(elements3,1),1);
end

middle_points=face_middlepoint(coordinates,faces);

X=kron(ones(1,4),element2faces); Y=kron(element2faces,ones(1,4));
if (nargin~=5)
    %Crouzeix-Raviart stifness matrix computation
    Z=zeros(size(X));
    for j = 1:size(elements3,1)
        verticesTrans=middle_points(element2faces(j,:),:)';
        volume=volumes(j);
        PhiGrad = [1 1 1 1; verticesTrans]\[zeros(1,3);eye(3)];
        Alocal = volume * PhiGrad * PhiGrad';
        Z(j,:)= reshape(Alocal,1,numel(Alocal));
    end
end
M_CR=sparse(X,Y,spdiags(coefficient_difussion,0,size(coefficient_difussion,1),size(coefficient_difussion,1))*Z);