function N = face_face_normal_derivative_3D(V,E,Tri,Tet)

  n = size(V,1);
  k = size(E,1);
  t = size(Tri,1);
  f = size(Tet,1);
  N = sparse(t, t);
%  N = sym('n', [t, t]);
%  N = 0*N;
  
  C = 9 * cotangent(V, Tet);
%  C = sym('c', [f, 6]);
%  C = 9*C;

  [~,on_b] = on_boundary(Tet);
  [on_b_f,on_b_c] = find(on_b);
  
  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);
  
  f_on_b=size(on_b_f,1);
  
  for i=1:f_on_b
    
    % Vertex subscript
    iiV1 = mod(on_b_c(i)+1,4)+1;
    iiV2 = mod(on_b_c(i)+2,4)+1;
    iiV3 = mod(on_b_c(i)+0,4)+1;
    iiV4 = mod(on_b_c(i)+3,4)+1;  % = on_b_c(i)
    
    % Vertex indices
    iV1 = Tet(on_b_f(i), iiV1);
    iV2 = Tet(on_b_f(i), iiV2);
    iV3 = Tet(on_b_f(i), iiV3);
    iV4 = Tet(on_b_f(i), iiV4);  % not on the boundary
    
    % Triangles
    T1 = sort([iV2, iV3, iV4]);
    T2 = sort([iV3, iV4, iV1]);
    T3 = sort([iV4, iV1, iV2]);
    T4 = sort([iV1, iV2, iV3]); % on the boundary
    
    % 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 = sort([iV2, iV3]);
%    E31 = sort([iV3, iV1]);
%    E12 = sort([iV1, iV2]);
%    E41 = sort([iV4, iV1]);
%    E42 = sort([iV4, iV2]);
%    E43 = sort([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);
%    
%    % Common edge between faces
%    common_edge = [
%         0, iE43, iE42, iE23;
%      iE43,    0, iE41, iE31;
%      iE42, iE41,    0, iE12;
%      iE23, iE31, iE12,    0 ];
    
    % Edge indices in C
    iE23c = reverse_edges_order(iiV2,iiV3);
    iE31c = reverse_edges_order(iiV3,iiV1);
    iE12c = reverse_edges_order(iiV1,iiV2);
    iE41c = reverse_edges_order(iiV4,iiV1);
    iE42c = reverse_edges_order(iiV4,iiV2);
    iE43c = reverse_edges_order(iiV4,iiV3);
    
    % 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;
    
    % Normal derivative matrix filling
    N(iT4, iT1) = N(iT4, iT1) + C(on_b_f(i), common_edge_indice(4,1));
    N(iT4, iT2) = N(iT4, iT2) + C(on_b_f(i), common_edge_indice(4,2));
    N(iT4, iT3) = N(iT4, iT3) + C(on_b_f(i), common_edge_indice(4,3));
    N(iT4, iT4) = N(iT4, iT4) - C(on_b_f(i), common_edge_indice(4,1)) - C(on_b_f(i), common_edge_indice(4,2)) - C(on_b_f(i), common_edge_indice(4,3));
    
  endfor
  
endfunction
