function L = face_face_laplacian_3D(V, E,Tri,Tet)

  n = size(V,1);
  k = size(E,1);
  t = size(Tri,1);
  f = size(Tet,1);
  L = sparse(t, t);
%  L = sym('l', [t, t]);
%  L = 0*L;
  
  assert(size(Tet, 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];
  
  hash = Tri * [1;n+1;0];
  reverse_triangles = sparse(Tri(:,3), hash, (1:t)',n+1,(n+1)^2);
  
  
  C = 9 * cotangent(V, Tet);
%  C = sym('c', [f, 6]);
%  C = 9*C;
  
  
  for i=1:f
    
    % Vertex indices
    iV1 = Tet(i, 1);
    iV2 = Tet(i, 2);
    iV3 = Tet(i, 3);
    iV4 = Tet(i, 4);
    
    % Triangles
    T1 = sort([iV2, iV3, iV4]);
    T2 = sort([iV3, iV4, iV1]);
    T3 = sort([iV4, iV1, iV2]);
    T4 = sort([iV1, iV2, iV3]);
    
    % Triangle indices in Tri
    iT1 = reverse_triangles(T1(3),T1 * [1;n+1;0]);
    iT2 = reverse_triangles(T2(3),T2 * [1;n+1;0]);
    iT3 = reverse_triangles(T3(3),T3 * [1;n+1;0]);
    iT4 = reverse_triangles(T4(3),T4 * [1;n+1;0]);
    
%    % Edges
%    E23 = [iV2, iV3];
%    E31 = [iV3, iV1];
%    E12 = [iV1, iV2];
%    E41 = [iV4, iV1];
%    E42 = [iV4, iV2];
%    E43 = [iV4, iV3];
%    
%    % Edge indices in F
%    iE23 = reverse_edges(iV2,iV3);
%    iE31 = reverse_edges(iV3,iV1);
%    iE12 = reverse_edges(iV1,iV2);
%    iE41 = reverse_edges(iV4,iV1);
%    iE42 = reverse_edges(iV4,iV2);
%    iE43 = reverse_edges(iV4,iV3);
    
    % Edge indices in C
    iE23c = reverse_edges_order(2,3);
    iE31c = reverse_edges_order(3,1);
    iE12c = reverse_edges_order(1,2);
    iE41c = reverse_edges_order(4,1);
    iE42c = reverse_edges_order(4,2);
    iE43c = reverse_edges_order(4,3);
    
    % Common edge between faces
    common_edge_indice = [
          0, iE43c, iE42c, iE23c;
      iE43c,     0, iE41c, iE31c;
      iE42c, iE41c,     0, iE12c;
      iE23c, iE31c, iE12c,     0 ];
    
    % Getting opposite edge (edges are stored like that in C)
    common_edge_indice = mod(common_edge_indice+2,6)+1;
    
    
    % Laplacian matrix filling
    L(iT1, iT1) = L(iT1, iT1) + C(i,common_edge_indice(1,2)) + C(i,common_edge_indice(1,3)) + C(i,common_edge_indice(1,4));
    L(iT1, iT2) = L(iT1, iT2) - C(i,common_edge_indice(1,2));
    L(iT1, iT3) = L(iT1, iT3) - C(i,common_edge_indice(1,3));
    L(iT1, iT4) = L(iT1, iT4) - C(i,common_edge_indice(1,4));
    
    L(iT2, iT1) = L(iT2, iT1) - C(i,common_edge_indice(2,1));
    L(iT2, iT2) = L(iT2, iT2) + C(i,common_edge_indice(2,1)) + C(i,common_edge_indice(2,3)) + C(i,common_edge_indice(2,4));
    L(iT2, iT3) = L(iT2, iT3) - C(i,common_edge_indice(2,3));
    L(iT2, iT4) = L(iT2, iT4) - C(i,common_edge_indice(2,4));
    
    L(iT3, iT1) = L(iT3, iT1) - C(i,common_edge_indice(3,1));
    L(iT3, iT2) = L(iT3, iT2) - C(i,common_edge_indice(3,2));
    L(iT3, iT3) = L(iT3, iT3) + C(i,common_edge_indice(3,1)) + C(i,common_edge_indice(3,2)) + C(i,common_edge_indice(3,4));
    L(iT3, iT4) = L(iT3, iT4) - C(i,common_edge_indice(3,4));
    
    L(iT4, iT1) = L(iT4, iT1) - C(i,common_edge_indice(4,1));
    L(iT4, iT2) = L(iT4, iT2) - C(i,common_edge_indice(4,2));
    L(iT4, iT3) = L(iT4, iT3) - C(i,common_edge_indice(4,3));
    L(iT4, iT4) = L(iT4, iT4) + C(i,common_edge_indice(4,1)) + C(i,common_edge_indice(4,2)) + C(i,common_edge_indice(4,3));
    
  endfor
  
  
endfunction
