function S=stress_tensor_local(V,E,bd_box_sz)
%S - 3-dimensional array,
%    S(:,:,1) - (1,1) component of stress
%    S(:,:,2) - (2,2) component of stress
%    S(:,:,3) - (2,1) component of stress
%V - list of vertex indices
%E - List of connected edges associated with V
%bd_box_sz - size of box to average over to do local stress computation
%
%For bugs, comments, etc. please contact Kapil Krishan :
%kkrishan@uci.edu
%www.physics.uci.edu/~foams
%version: I

l_vec=V(E(:,1),:)-V(E(:,2),:);

volume=(bd_box_sz+1)^2;

for bxi=floor(min(V(E(:),1))):ceil(max(V(E(:),1)))
    for bxj=floor(min(V(E(:),2))):ceil(max(V(E(:),2)))
        coord=[bxi bxj];

        [V_idx_in,V_idx_out,edgs_bd_idx,edgs_in_idx]=sub_domain(coord,bd_box_sz,E,V);
        f_edge_in=[ones(size(edgs_in_idx,1),1) edgs_in_idx];
        
        f_edge_bd=f_edge_bdry(coord,bd_box_sz,E,V);
        f_edge=[f_edge_in;f_edge_bd];
        
        l_in=l_vec(f_edge(:,2),:);

        sigma_cap=(([f_edge(:,1) f_edge(:,1)]'.*l_in')*(l_in./[sqrt(sum(l_in.^2,2)) sqrt(sum(l_in.^2,2))]))/volume;
        sig_11(bxi,bxj)=sigma_cap(1,1);
        sig_22(bxi,bxj)=sigma_cap(2,2);
        sig_21(bxi,bxj)=sigma_cap(2,1);
    end
%     disp(bxi);
end
S(:,:,1)=sig_11;
S(:,:,2)=sig_22;
S(:,:,3)=sig_21;

function edgs=interior_edgs(C,bdry_cells)
ct=0;
for ic=setdiff([1:size(C,1)],bdry_cells')
    ct=ct+1;
    C_interior{ct}=C{ic};
end
C_interior=C_interior';

edgs=[0 0];
for ic=1:size(C_interior,1)
    vtxs=C_interior{ic};
    edgs=[edgs;setdiff(sort([vtxs;[vtxs(2:end) vtxs(1)]])',edgs,'rows')];
end
edgs(1,:)=[];

function [V_idx_in,V_idx_out,edgs_bd_idx,edgs_in_idx]=sub_domain(coord,bd_box_sz,edgs,V)

min_1=coord(1)-bd_box_sz;
max_1=coord(1)+bd_box_sz;
min_2=coord(2)-bd_box_sz;
max_2=coord(2)+bd_box_sz;

tvar=(V(:,1)<max_1).*(V(:,1)>min_1).*(V(:,2)<max_2).*(V(:,2)>min_2);
V_idx_in=find(tvar==1);

tvar=ismember(edgs,V_idx_in);
edgs_bd_idx=find(sum(tvar,2)==1);
edgs_in_idx=find(sum(tvar,2)==2);

V_all=unique(edgs([edgs_bd_idx;edgs_in_idx],:));
V_idx_out=setdiff(V_all(:),V_idx_in);

function f_edge=f_edge_bdry(coord,bd_box_sz,edgs,V)

edge_coord=[V(edgs(:,1),1) V(edgs(:,2),1) V(edgs(:,1),2) V(edgs(:,2),2)];

bx_edge=[coord(1)-bd_box_sz*[1 -1] coord(2)-bd_box_sz*[1 1]];
edge_coord_shft=edge_coord-ones(size(edge_coord,1),1)*[0 0 bx_edge(3)*[1 1]];
tvar=(edge_coord_shft(:,[3 4])>0)-(edge_coord_shft(:,[3 4])<0);
c_idx=find(sum(tvar,2)==0);
u=edge_coord(c_idx,:);
u_intercept=((u(:,1)-u(:,2))./(u(:,3)-u(:,4))).*(bx_edge(3)-u(:,3))+u(:,1);
u_intercept=[u_intercept ones(size(u_intercept,1),1)*bx_edge(3)];
d_idx=find((u_intercept(:,1)>bx_edge(1)).*(u_intercept(:,1)<bx_edge(2))==1);
edgs_no=c_idx(d_idx);
ttvar=(diff(tvar(c_idx(d_idx),:)')')/2;
e_idx=3-2.^((1-ttvar)/2);
u_intercept_edg=u_intercept(d_idx,:);
edgs_cord_inside=zeros(size(e_idx,1),2);
for ic=1:size(e_idx,1)
    edgs_cord_inside(ic,:)=V(edgs(c_idx(d_idx(ic)),e_idx(ic)),:);
end
edge_1=[u_intercept_edg edgs_no e_idx];

% line(bx_edge([1 2]),bx_edge([3 4]),'color','c');
% plot(edge_coord(c_idx,1),edge_coord(c_idx,3),'g.')
% plot(edge_coord(c_idx,2),edge_coord(c_idx,4),'g.')
% plot(u_intercept(:,1),u_intercept(:,2),'m.')
% plot(u_intercept_edg(:,1),u_intercept_edg(:,2),'k*')
% plot(edgs_cord_inside(:,1),edgs_cord_inside(:,2),'r*')
% drawnow;

bx_edge=[coord(1)-bd_box_sz*[1 -1] coord(2)+bd_box_sz*[1 1]];
edge_coord_shft=edge_coord-ones(size(edge_coord,1),1)*[0 0 bx_edge(3)*[1 1]];
tvar=(edge_coord_shft(:,[3 4])>0)-(edge_coord_shft(:,[3 4])<0);
c_idx=find(sum(tvar,2)==0);
u=edge_coord(c_idx,:);
u_intercept=((u(:,1)-u(:,2))./(u(:,3)-u(:,4))).*(bx_edge(3)-u(:,3))+u(:,1);
u_intercept=[u_intercept ones(size(u_intercept,1),1)*bx_edge(3)];
d_idx=find((u_intercept(:,1)>bx_edge(1)).*(u_intercept(:,1)<bx_edge(2))==1);
edgs_no=c_idx(d_idx);
ttvar=(diff(tvar(c_idx(d_idx),:)')')/2;
e_idx=2.^((1-ttvar)/2);
u_intercept_edg=u_intercept(d_idx,:);
edgs_cord_inside=zeros(size(e_idx,1),2);
for ic=1:size(e_idx,1)
    edgs_cord_inside(ic,:)=V(edgs(c_idx(d_idx(ic)),e_idx(ic)),:);
end
edge_2=[u_intercept_edg edgs_no e_idx];

% line(bx_edge([1 2]),bx_edge([3 4]),'color','y');
% plot(edge_coord(c_idx,1),edge_coord(c_idx,3),'g.')
% plot(edge_coord(c_idx,2),edge_coord(c_idx,4),'g.')
% plot(u_intercept(:,1),u_intercept(:,2),'m.')
% plot(edge_2(:,1),edge_2(:,2),'k*')
% plot(edgs_cord_inside(:,1),edgs_cord_inside(:,2),'r*')
% drawnow;

bx_edge=[coord(1)-bd_box_sz*[1 1] coord(2)+bd_box_sz*[1 -1]];
edge_coord_shft=edge_coord-ones(size(edge_coord,1),1)*[bx_edge(1)*[1 1] 0 0];
tvar=(edge_coord_shft(:,[1 2])>0)-(edge_coord_shft(:,[1 2])<0);
c_idx=find(sum(tvar,2)==0);
u=edge_coord(c_idx,:);
u_intercept=((u(:,3)-u(:,4))./(u(:,1)-u(:,2))).*(bx_edge(1)-u(:,1))+u(:,3);
u_intercept=[ones(size(u_intercept,1),1)*bx_edge(1) u_intercept];
d_idx=find((u_intercept(:,2)<bx_edge(3)).*(u_intercept(:,2)>bx_edge(4))==1);
edgs_no=c_idx(d_idx);
ttvar=(diff(tvar(c_idx(d_idx),:)')')/2;
e_idx=3-2.^((1-ttvar)/2);
u_intercept_edg=u_intercept(d_idx,:);
edgs_cord_inside=zeros(size(e_idx,1),2);
for ic=1:size(e_idx,1)
    edgs_cord_inside(ic,:)=V(edgs(c_idx(d_idx(ic)),e_idx(ic)),:);
end
edge_3=[u_intercept_edg edgs_no e_idx];

% line(bx_edge([1 2]),bx_edge([3 4]),'color','c');
% plot(edge_coord(c_idx,1),edge_coord(c_idx,3),'g.')
% plot(edge_coord(c_idx,2),edge_coord(c_idx,4),'g.')
% plot(u_intercept(:,1),u_intercept(:,2),'m.')
% plot(edge_3(:,1),edge_3(:,2),'k*')
% plot(edgs_cord_inside(:,1),edgs_cord_inside(:,2),'r*')
% drawnow;


bx_edge=[coord(1)+bd_box_sz*[1 1] coord(2)+bd_box_sz*[1 -1]];
edge_coord_shft=edge_coord-ones(size(edge_coord,1),1)*[bx_edge(1)*[1 1] 0 0];
tvar=(edge_coord_shft(:,[1 2])>0)-(edge_coord_shft(:,[1 2])<0);
c_idx=find(sum(tvar,2)==0);
u=edge_coord(c_idx,:);
u_intercept=((u(:,3)-u(:,4))./(u(:,1)-u(:,2))).*(bx_edge(1)-u(:,1))+u(:,3);
u_intercept=[ones(size(u_intercept,1),1)*bx_edge(1) u_intercept];
d_idx=find((u_intercept(:,2)<bx_edge(3)).*(u_intercept(:,2)>bx_edge(4))==1);
edgs_no=c_idx(d_idx);
ttvar=(diff(tvar(c_idx(d_idx),:)')')/2;
e_idx=2.^((1-ttvar)/2);
u_intercept_edg=u_intercept(d_idx,:);
edgs_cord_inside=zeros(size(e_idx,1),2);
for ic=1:size(e_idx,1)
    edgs_cord_inside(ic,:)=V(edgs(c_idx(d_idx(ic)),e_idx(ic)),:);
end
edge_4=[u_intercept_edg edgs_no e_idx];

% line(bx_edge([1 2]),bx_edge([3 4]),'color','y');
% plot(edge_coord(c_idx,1),edge_coord(c_idx,3),'g.')
% plot(edge_coord(c_idx,2),edge_coord(c_idx,4),'g.')
% plot(u_intercept(:,1),u_intercept(:,2),'m.')
% plot(edge_4(:,1),edge_4(:,2),'k*')
% plot(edgs_cord_inside(:,1),edgs_cord_inside(:,2),'r*')
% drawnow;

e=[edge_1;edge_2;edge_3;edge_4];
f_edge=[];
for ic=1:size(e,1)
    edg=e(ic,3);
    v_in=V(edgs(e(ic,3),e(ic,4)),:);
    v_out=V(edgs(e(ic,3),3-e(ic,4)),:);

    idx=find(e(:,3)==e(ic,3));
    v_bd_1=e(idx(1),[1 2]);
    if length(idx)==1
        f=norm(v_in-v_bd_1)/norm(v_in-v_out);
    end
    if length(idx)==2
        v_bd_1=e(idx(1),[1 2]);
        v_bd_2=e(idx(2),[1 2]);
        f=norm(v_bd_1-v_bd_2)/norm(v_in-v_out);
    end
    f_edge=[f_edge;f edg];
end
f_edge=unique(f_edge,'rows');