function L = edge_edge_laplacian_3D(V,E,F)

  n = length(V);
  f = size(F,1);
  k = length(E);
  L = sparse(k, k);
  
  assert(size(F, 2) == 4);
  
  reverse_edges = sparse(E(:,1), E(:,2), 1:k, k, k);
  reverse_edges = reverse_edges + reverse_edges';
  
%  reverse_edges_order = [
%    0, 3, 2, 4; 
%    3, 0, 1, 5; 
%    2, 1, 0, 6; 
%    4, 5, 6, 0];
  
  reverse_edges_order = [
    0, 6, 5, 1;
    6, 0, 4, 2;
    5, 4, 0, 3;
    1, 2, 3, 0];
    
%  reverse_edges_order = mod(reverse_edges_order + 2,6)+1;
    
    
  C = 4 * cotangent(V, F);
  
  for i=1:f
    
    % Vertex indices
    iV1 = F(i, 1);
    iV2 = F(i, 2);
    iV3 = F(i, 3);
    iV4 = F(i, 4);
    
    % Edges
    E1 = [iV2, iV3];
    E2 = [iV3, iV1];
    E3 = [iV1, iV2];
    E4 = [iV4, iV1];
    E5 = [iV4, iV2];
    E6 = [iV4, iV3];
    
    % Edge indices in F
    iE1 = reverse_edges(iV2,iV3);
    iE2 = reverse_edges(iV3,iV1);
    iE3 = reverse_edges(iV1,iV2);
    iE4 = reverse_edges(iV4,iV1);
    iE5 = reverse_edges(iV4,iV2);
    iE6 = reverse_edges(iV4,iV3);
    
    % Edge indices in C
    iE1c = reverse_edges_order(2,3);
    iE2c = reverse_edges_order(3,1);
    iE3c = reverse_edges_order(1,2);
    iE4c = reverse_edges_order(4,1);
    iE5c = reverse_edges_order(4,2);
    iE6c = reverse_edges_order(4,3);
    
    
    % Laplacian matrix filling
    L(iE1, iE1) = L(iE1, iE1) + C(i,iE2c) + C(i,iE3c) + C(i,iE5c) + C(i,iE6c);
    L(iE1, iE2) = L(iE1, iE2) - C(i,iE2c);
    L(iE1, iE3) = L(iE1, iE3) - C(i,iE3c);
%    L(iE1, iE4) = L(iE1, iE4);
    L(iE1, iE5) = L(iE1, iE5) - C(i,iE5c);
    L(iE1, iE6) = L(iE1, iE6) - C(i,iE6c);
    
    L(iE2, iE1) = L(iE2, iE1) - C(i,iE1c);
    L(iE2, iE2) = L(iE2, iE2) + C(i,iE1c) + C(i,iE3c) + C(i,iE4c) + C(i,iE6c);
    L(iE2, iE3) = L(iE2, iE3) - C(i,iE3c);
    L(iE2, iE4) = L(iE2, iE4) - C(i,iE4c);
%    L(iE2, iE5) = L(iE2, iE5);
    L(iE2, iE6) = L(iE2, iE6) - C(i,iE6c);
    
    L(iE3, iE1) = L(iE3, iE1) - C(i,iE1c);
    L(iE3, iE2) = L(iE3, iE2) - C(i,iE2c);
    L(iE3, iE3) = L(iE3, iE3) + C(i,iE1c) + C(i,iE2c) + C(i,iE4c) + C(i,iE5c);
    L(iE3, iE4) = L(iE3, iE4) - C(i,iE4c);
    L(iE3, iE5) = L(iE3, iE5) - C(i,iE5c);
%    L(iE3, iE6) = L(iE3, iE6);
    
%    L(iE4, iE1) = L(iE4, iE1);
    L(iE4, iE2) = L(iE4, iE2) - C(i,iE2c);
    L(iE4, iE3) = L(iE4, iE3) - C(i,iE3c);
    L(iE4, iE4) = L(iE4, iE4) + C(i,iE2c) + C(i,iE3c) + C(i,iE5c) + C(i,iE6c);
    L(iE4, iE5) = L(iE4, iE5) - C(i,iE5c);
    L(iE4, iE6) = L(iE4, iE6) - C(i,iE6c);
    
    L(iE5, iE1) = L(iE5, iE1) - C(i,iE1c);
%    L(iE5, iE2) = L(iE5, iE2);
    L(iE5, iE3) = L(iE5, iE3) - C(i,iE3c);
    L(iE5, iE4) = L(iE5, iE4) - C(i,iE4c);
    L(iE5, iE5) = L(iE5, iE5) + C(i,iE1c) + C(i,iE3c) + C(i,iE4c) + C(i,iE6c);
    L(iE5, iE6) = L(iE5, iE6) - C(i,iE6c);
    
    L(iE6, iE1) = L(iE6, iE1) - C(i,iE1c);
    L(iE6, iE2) = L(iE6, iE2) - C(i,iE2c);
%    L(iE6, iE3) = L(iE6, iE3);
    L(iE6, iE4) = L(iE6, iE4) - C(i,iE4c);
    L(iE6, iE5) = L(iE6, iE5) - C(i,iE5c);
    L(iE6, iE6) = L(iE6, iE6) + C(i,iE1c) + C(i,iE2c) + C(i,iE4c) + C(i,iE5c);
    
  endfor
  
  
endfunction
