function [ matrixvsfra,connectmf,fracnumber,fracstart,corevsfra,connect_infrac,N,T,zf,fcinff,vf,porf,fcff] = connections_PEDFM_EX1_solutions( fracturemesh,pointvsregion,raopoint,f,nf,nx,ny,nz,kx,ky,kz,kf,wf,Porf,dxv,dyv,dzv,coord,nodes,NTG,addflag,frac_information,coordinates)
% 该函数用以给裂缝面网格编号，并存储网格间的连接关系及传导率，是前处理的关键步骤
% 不同于二维嵌入式离散裂缝模型，三维模型中，裂缝是二维平面，网格之间连接关系复杂，
%不能简单得从基质网格来判定，若所处基质网格不相邻，则裂缝网格必不相邻，若所处基质网格相邻，裂缝网格却不一定相邻
%上述原则无助于我们去确定connections，故拟采取以下方案：
%第一步:按照xink矩阵去掉零行后，按照行数依次给某裂缝面网格编号，同时存储裂缝网格所在的基质网格序号
%第二步:存储每个裂缝网格的面积，边界，边界长度，中心点坐标，中心点到各边界的距离，为后续操作做准备
%第三步:按照裂缝网格序号顺序，按照边界进行搜索，如果有相同的边界，则两个裂缝网格是相邻的
%matrixvsfra矩阵代表基质网格中包含裂缝网格编号的情况
%fracnumber矩阵，行序号表示该裂缝网格编号，行内容表示该网格包含的点序号
%areavsfra矩阵，行序号表示该裂缝网格编号，行内容表示该裂缝网格的面积
%每相邻两列是一条边，然后将顺序反过来，又是两列一条边
%lengthvsfra矩阵，行序号表示该裂缝网格编号，行内容表示该裂缝网格的各边长度
%corevsfra矩阵，行序号表示该裂缝网格编号，行内容表示该裂缝网格的重心坐标
%disvsfra矩阵，行序号表示该裂缝网格编号，行内容表示该裂缝网格的重心到各边的距离
matrixvsfra=zeros(nx*ny*nz,10);%为了减少该矩阵的大小，考虑实际情况，一个基质网格中，一般不会有超过10个裂缝网格
fracnumber=zeros(1000,20);%根据实际情况，裂缝网格一般不会超过1000个，也可根据实际情况改变
%m=size(xink,1);
rownum=zeros(nf,1);
for i=1:nf
    rownum(i)=size(fracturemesh{i,1},1);
end
%% 将所有裂缝面网格连起来一起编号，并将编号扔给相应的基质网格
q=1;
fracstart=zeros(nf+1,1);%这个表示在总编号中，每条裂缝起始的号码
for j=1:nf
    fracstart(j)=q;
    numofmesh=pointvsregion{j,1};
    for i=1:1:rownum(j)
        if(norm(fracturemesh{j,1}(i,:))~=0)
            raoindice=size(fracturemesh{j,1}(i,:),2);
            fracnumber(q,1:raoindice)=fracturemesh{j,1}(i,:);%q即是该裂缝编号
            m=floor((i+(2^addflag(j)-1))/(2^addflag(j)));
            indice=find(matrixvsfra(numofmesh(m),:)~=0);
            %         indice=find(matrixvsfra(numofmesh(i),:)~=0);
            if (size(indice,2)==0) xiang=1;
            else
                rao=max(indice');
                xiang=rao(1)+1;
            end
            matrixvsfra(numofmesh(m),xiang)=q;%将裂缝编号附给相应的基质网格
            q=q+1;
        end
    end
end
fracstart(nf+1)=q;%裂缝网格总数量+1
fracnumber(q:1000,:)=[];


%% 将包含有裂缝单元的基质网格筛选出来,并进行处理
m=size(matrixvsfra,1);
connectmf=cell(1,4);
% for i=1:m
%     if (norm(matrixvsfra(i,:))~=0)
% ConnecS = cell(1,4);
cc = 0;
cl = 1;
for i = 1 : m
    %筛选出包含裂缝的基质网格
    %     if ~isempty(matrixvsfra(i,:))
    if (norm(matrixvsfra(i,:)))~=0
        cc = cc + 1;
        %存储该基质网格编号
        connectmf{cc,1} = i;
        %给这些裂缝点/裂缝单元按顺序编号
        A=matrixvsfra(i,:);
        A(A==0)=[];
        connectmf{cc,2} =A;
    end
end

%% 一条条裂缝计算裂缝网格面积、边界等等
% edgevsfra=cell(nf,1);
%
fractureCellType = zeros(q-1,1);
theNumber = 3;
fractureCellType = repmat((1:theNumber)',(q-1)/theNumber,1);
%
normalvec=zeros(q-1,3);%unit 法向量
areavsfra=zeros(q-1,1);%q-1即裂缝网格总数
lengthvsfra=zeros(q-1,10);
Kf=zeros(q-1,1);%裂缝网格渗透率
Wf=zeros(q-1,1);%裂缝单元缝宽
corevsfra=zeros(q-1,3);
disvsfra=zeros(q-1,10);
tran=zeros(q-1,10);
avtran=zeros(q-1,1);%裂缝网格平均传导，用于后续的简化计算
fcinff=zeros(q-1,1);%表示裂缝单元所在的裂缝面序号
zf=zeros(q-1,1);%裂缝单元高度（以下表面为基准）
vf=zeros(q-1,1);%裂缝单元的体积
porf=zeros(q-1,1);%裂缝单元的孔隙度
for i=1:1:q-1 %对裂缝网格进行操作
    if(norm(fracnumber(i,:))==0) break;
    else
        for k=2:(nf+1)
            if(i<fracstart(k)) fcinff(i)=k-1; break;%判断该点处在哪个裂缝平面上,k-1号裂缝面
            end
        end
        d3intersection=raopoint{k-1,1};
        anothersection=raopoint{k-1,2};
        fraindice=find(fracnumber(i,:)~=0);
        rx=size(fraindice,2);
        for j=1:1:(rx-1)
            %edgevsfra(i,(4*j-3):(4*j))=[fracnumber(i,j),fracnumber(i,j+1),fracnumber(i,j+1),fracnumber(i,j)];%完成对edgevsfra矩阵赋值
            %%因为是凸多边形，因此只要代表裂缝网格的两个行向量的共有元素有2个，即有公共边；
            %如果多于2个，则是相同的裂缝网格；如果小于两个，则没有共有的边线
            normalvec(i,:)=cross(f(5*(k-1)-3,:),f(5*(k-1)-2,:))/norm(cross(f(5*(k-1)-3,:),f(5*(k-1)-2,:)));
            areavsfra(i)=area(anothersection(fracnumber(i,1:(rx-1))',:) ,f(5*(k-1)-3,:),f(5*(k-1)-2,:) );%计算出每个裂缝网格的面积
            vf(i)=areavsfra(i)*wf(k-1);
            porf(i)=Porf(k-1);
            Wf(i,1)=wf(k-1);
            Kf(i,1)=kf(k-1);
            lengthvsfra(i,j)=norm(d3intersection(fracnumber(i,j+1),:)-d3intersection(fracnumber(i,j),:));%计算裂缝网格每条边的长度
            corevsfra(i,:)=1/(rx-1)*d3intersection(fracnumber(i,j),:)+corevsfra(i,:);%计算每个裂缝网格重心的坐标
            zf(i)=corevsfra(i,3);%裂缝网格计算深度（以下平面为基准计算得到的高度，值为正数）
            rao1=corevsfra(i,:)-d3intersection(fracnumber(i,j),:);
            rao2=d3intersection(fracnumber(i,j+1),:)-d3intersection(fracnumber(i,j),:);
            disvsfra(i,j)=norm(cross(rao1',rao2'))/norm(rao2);%计算出每个裂缝网格重心到每条边的距离
            tran(i,j)=kf(k-1)*wf(k-1)*lengthvsfra(i,j)./disvsfra(i,j);
        end
        avtran(i)=sum(tran(i,:))/(rx-1);
    end
end

% %% 生成裂缝单元及所在裂缝面的cell数组及其它形状信息
% fcff=cell(1,8);
% fcff{1,1}=fracnumber; fcff{1,2}=fcinff;fcff{1,3}=lengthvsfra;fcff{1,4}=corevsfra;fcff{1,5}=disvsfra;fcff{1,6}=tran;fcff{1,7}=avtran;fcff{1,8}=areavsfra;
%% 为了防止出现某裂缝网格面积相对过小，引起计算出错，故将面积很小的裂缝网格给去掉
maxarea=max(areavsfra,2);
% minlength=min(lengthvsfra,[],2);
% maxlength=max(lengthvsfra,[],2);
newnum=zeros(q-1,1);
[deindex1,~]=find(areavsfra<0*maxarea);
deindex2=[];
for i=1:(q-1)
    zeroindex=find(lengthvsfra(i,:)~=0);
    minlength=min(lengthvsfra(i,zeroindex));
    maxlength=max(lengthvsfra(i,:));
    if minlength<0*maxlength
        deindex2=[deindex2;i];
    end
end
deindex=[deindex1;deindex2];
deindex=unique(deindex);
deindex=sort(deindex);
hui=length(deindex);
if hui>0
    for i=1:(q-1)
        [dogindex,~]=find(deindex==i);
        if length(dogindex)~=0
            newnum(i)=0;
        else
            [catindex,~]=find(deindex<i);
            rao=length(catindex);
            newnum(i)=i-rao;
        end
    end
    fracnumber(deindex,:)=[];
    fcinff(deindex,:)=[];
    lengthvsfra(deindex,:)=[];
    corevsfra(deindex,:)=[];
    disvsfra(deindex,:)=[];
    tran(deindex,:)=[];
    avtran(deindex,:)=[];
    areavsfra(deindex,:)=[];
    zf(deindex,:)=[];
    vf(deindex,:)=[];
    porf(deindex,:)=[];
    % change matrixvsfra
    for i=1:m%matrixvsfra矩阵的行数
        for j=1:10
            if matrixvsfra(i,j)~=0
                matrixvsfra(i,j)=newnum(matrixvsfra(i,j));
            end
        end
    end
    % change fracstart
    fracstart(nf+1,1)=length(fracnumber)+1;
    for i=1:nf
        [daiindex,~]=find(deindex<fracstart(i,1));
        dai=length(daiindex);
        fracstart(i)=fracstart(i)-dai;
    end
    
end
%% 生成裂缝单元及所在裂缝面的cell数组及其它形状信息
fcff=cell(1,8);
fcff{1,1}=fracnumber; fcff{1,2}=fcinff;fcff{1,3}=lengthvsfra;fcff{1,4}=corevsfra;fcff{1,5}=disvsfra;fcff{1,6}=tran;fcff{1,7}=avtran;fcff{1,8}=areavsfra;
%% 裂缝单元的连接情况

%% 在同一裂缝面上，裂缝单元的连接情况及传导率计算(不包含在同一基质网格中的相邻裂缝单元）
connect_infrac=cell(nf,1);
Nff=zeros(1,2); %记录同一裂缝面上的网格相邻情况，第一列和第二列分别是相邻网格的编号
Tff=zeros(1,1); %记录相应的传导系数
for i=1:nf
    p=1;
    raoconnect=zeros(fracstart(i+1)-fracstart(i),10);%一般来说，裂缝网格不超过10条边
    d3intersection=raopoint{i,1};
    anothersection=raopoint{i,2};
    for j=fracstart(i):((fracstart(i+1)-1)-1)%第i条裂缝的编号范围
        q=1;
        for k=j:(fracstart(i+1)-1)
            A=fracnumber(j,:);
            B=fracnumber(k,:);
            %将两个向量很可能都有的0去掉，这个0没有什么意义
            A(A==0)=[];
            B(B==0)=[];
            [a,~]=find(matrixvsfra==j); [b,~]=find(matrixvsfra==k);
            if length(a)~=1
                heihei=1;
            end
            raoflag=intersect(A,B);
            if((size(raoflag,2)>2)||((size(raoflag,2)<2))) continue;%因为是凸多边形，因此只要代表裂缝网格的行向量的共有元素有2个，即有公共边；
                %如果多于2个，则是相同的裂缝网格，因为如果不是，则必有一个是凹多边形，矛盾，证毕；如果小于两个，则没有共有的边线
            else if (a==b) continue; %说明这两个裂缝单元在同一个基质网格中，因此不应该算在此类中
                else
                    raoconnect(p,q)=k;q=q+1;
                    lengthcat=norm(d3intersection(raoflag(1),:)-d3intersection(raoflag(2),:));%计算裂缝单元公共边长度
                    rao1k=corevsfra(k,:)-d3intersection(raoflag(1),:);
                    rao1j=corevsfra(j,:)-d3intersection(raoflag(1),:);
                    rao2=d3intersection(raoflag(1),:)-d3intersection(raoflag(2),:);
                    %分别计算两裂缝单元中心到公共边的距离
                    disk=norm(cross(rao1k',rao2'))/norm(rao2);
                    disj=norm(cross(rao1j',rao2'))/norm(rao2);
                    if Kf(j)==0 || Kf(k)==0
                        hartran=0;
                    else
                        tran1=Kf(j)*wf(i)*lengthcat/disk;
                        tran2=Kf(k)*wf(i)*lengthcat/disj;
                        %             hartran=length/(disk+disj);%取调和平均的一半
                        hartran=tran1*tran2/(tran1+tran2);
                    end
                    Nff=[Nff; j,k];
                    Tff=[Tff;hartran];
                    %计算出每个裂缝网格重心到每条边的距离
                end
            end
            p=p+1;
        end
        connect_infrac{i,1}=raoconnect;
    end
end
Nff(1,:)=[]; %去除第一行的零行
Tff(1,:)=[];
%% 在同一基质网格中的裂缝单元连接情况及传导率计算
nmc=size(matrixvsfra,1);
Nmff=zeros(1,2);
Tmff=zeros(1,1);
for i=1:nmc
    if(norm(matrixvsfra(i,:))~=0)
        indice=find(matrixvsfra(i,:)~=0);
        n=size(indice,2);
        if (n>1) %此时表明要采取下述的简化算法
            %         lengthvsfra(matrixvsfra(i,indice))./disvsfra(matrixvsfra(i,indice))
            for j=1:n-1
                for k=(j+1):n
                    Nmff=[Nmff;matrixvsfra(i,indice(j)), matrixvsfra(i,indice(k))];
                    % 此时的计算方式是先求每个裂缝网格的算术平均，再算包含在该基质网格中的所有裂缝单元的算术平均
                    fc1=matrixvsfra(i,indice(j)); fc2=matrixvsfra(i,indice(k));
                    allfc=matrixvsfra(i,1:n); sumtran=sum(fcff{1,7}(allfc));
                    if sumtran==0 Tmff=[Tmff; 0];
                    else
                        Tmff=[Tmff; fcff{1,7}(fc1)*fcff{1,7}(fc2)/sumtran];
                    end
                end
            end
        end
    end
end
Nmff(1,:)=[]; %去掉开头的零行
Tmff(1,:)=[];
% N = [N; Nmff+nmc];
% T = [T; Tmff];

%% 基质网格的连接情况及传导系数计算
nmc=nx * ny*nz;%基质网格数目
nfc=size(fracnumber,1);%裂缝单元数目
rpt = [ones(nmc, 1); zeros(nfc, 1)];
% r.rpt = rpt;
nmm = (nx-1)*ny*nz+(ny-1)*nx*nz+(nz-1)*nx*ny;%nmm是基质网格之间存在流体交换的总数
% r.nf = nmm + nff;%nff是裂缝单元之间存在流体交换的总数
S_ex=zeros(nmm,1);% 流体交换面面积

Kx=kx.*ones(nmc,1);
Ky=ky.*ones(nmc,1);
Kz=kz.*ones(nmc,1);
N = zeros(nmm, 2);%存储基质网格之间存在流体交换的网格编号
T = zeros(nmm, 1);%存储对应与N矩阵的传导系数
c = 0;
for k= 1 : nz
    for j = 1 : ny
        for i = 1 : nx - 1
            c = c + 1;
            index = i + (j - 1) * nx +  (k - 1) * nx * ny;
            indexn = index + 1;
            N(c, :) = [index, indexn];
            T(c) = 2 * dzv(index)*dyv(index)*Kx(index)*Kx(indexn)*NTG(index)*NTG(indexn)/(Kx(index)*dxv(indexn)*NTG(index) + Kx(indexn)*dxv(index)*NTG(indexn));
            %         2 * dzv(index)*dyv(index)*Kx(index)*Kx(indexn)/(Kx(index)*dxv(indexn) + Kx(indexn)*dxv(index));
            S_ex(c)=dzv(index)*dyv(index);
        end
    end
end

for k= 1 : nz
    for i = 1 : nx
        for j = 1 : ny - 1
            c = c + 1;
            index = i + (j - 1) * nx +  (k - 1) * nx * ny;
            indexn = index + nx;
            N(c, :) = [index, indexn];
            T(c) = 2 * dzv(index)*dxv(index)*Ky(index)*Ky(indexn)/(Ky(index)*dyv(indexn) + Ky(indexn)*dyv(index));
            S_ex(c)=dzv(index)*dxv(index);
        end
    end
end

for i = 1 : nx
    for j = 1 : ny
        for k= 1: nz-1
            c = c + 1;
            index = i + (j - 1) * nx + (k - 1) * nx * ny;
            indexn = index + nx * ny;
            N(c, :) = [index, index + nx*ny];
            T(c) = 2 * dxv(index)*dyv(index)*Kz(index)*Kz(indexn)/(Kz(index)*dzv(indexn) + Kz(indexn)*dzv(index));
            S_ex(c)=dxv(index)*dyv(index);
        end
    end
end




%% 按照2014年的方法，将窜流处理为与上述相似的形式
syms x y z;
Ninterflow=[];
Tinterflow=[];
for i=1:nx*ny*nz  %基质网格编号
    for j=1:10
        if matrixvsfra(i,j)~=0  %this matrix cell contains a fracture cell
            Ninterflow=[Ninterflow;i,matrixvsfra(i,j)+nmc];
            kcell=(kx(i)*ky(i)*kz(i))^(1/3);
            Knnc=kcell*Kf(matrixvsfra(i,j))/(kcell+Kf(matrixvsfra(i,j)));
            Annc=areavsfra(matrixvsfra(i,j));% 此处与原始嵌入式不一样，原始是两倍
            fracore=corevsfra(matrixvsfra(i,j),:);
            matnodes=nodes(i,:);
            verco=coord(matnodes,:);norvec=normalvec(matrixvsfra(i,j),:);
            matcore=mean(verco);d0=matcore-fracore;
            dn=(x+d0(1))*norvec(1)+(y+d0(2))*norvec(2)+(z+d0(3))*norvec(3);
            %     dn=abs(dn);%影响数值积分效率
            dn=sqrt(dn^2);
            dn=matlabFunction(dn);
            if norvec(1)~=0 && norvec(2)~=0 && norvec(3)~=0
                Dn=integral3(dn,-dxv(i)/2,dxv(i)/2,-dyv(i)/2,dyv(i)/2,-dzv(i)*NTG(i)/2,dzv(i)*NTG(i)/2)/(dxv(i)*dyv(i)*dzv(i)*NTG(i));
            else if norvec(1)~=0 && norvec(2)~=0
                    Dn=dzv(i)*NTG(i)*integral2(dn,-dxv(i)/2,dxv(i)/2,-dyv(i)/2,dyv(i)/2)/(dxv(i)*dyv(i)*dzv(i)*NTG(i));
                else if norvec(1)~=0 && norvec(3)~=0
                        Dn=dyv(i)*integral2(dn,-dxv(i)/2,dxv(i)/2,-dzv(i)*NTG(i)/2,dzv(i)*NTG(i)/2)/(dxv(i)*dyv(i)*dzv(i)*NTG(i));
                    else if norvec(2)~=0 && norvec(3)~=0
                            Dn=dxv(i)*integral2(dn,-dyv(i)/2,dyv(i)/2,-dzv(i)*NTG(i)/2,dzv(i)*NTG(i)/2)/(dxv(i)*dyv(i)*dzv(i)*NTG(i));
                        else if norvec(1)~=0
                                Dn=dyv(i)*dzv(i)*NTG(i)*integral(dn,-dxv(i)/2,dxv(i)/2)/(dxv(i)*dyv(i)*dzv(i)*NTG(i));
                            else if norvec(2)~=0
                                    Dn=dxv(i)*dzv(i)*NTG(i)*integral(dn,-dyv(i)/2,dyv(i)/2)/(dxv(i)*dyv(i)*dzv(i)*NTG(i));
                                else if norvec(3)~=0
                                        Dn=dxv(i)*dyv(i)*integral(dn,-dzv(i)*NTG(i)/2,dzv(i)*NTG(i)/2)/(dxv(i)*dyv(i)*dzv(i)*NTG(i));
                                    end
                                end
                            end
                        end
                    end
                end
            end
            % 为论文做修改
            %     Dn=0.9;
            %     Dn=2.5;
            raoT=Knnc*Annc/Dn;
            Tinterflow=[Tinterflow;raoT];
        end
    end
end

%% 投影嵌入式处理
nmc=size(matrixvsfra,1);
%基质网格流体交换面上的裂缝单元投影面之和
S_p=zeros(nmm,1);
% 因为投影处理新加的F-M connections
Nmf_projection=[];
Tmf_projection=[];
for i=1:nmc
    matnodes=nodes(i,:);
    dx_grid=dxv(i);
    dy_grid=dyv(i);
    dz_grid=dzv(i);
    grid_core=mean(coord(matnodes,:));%该基质网格中心
    if(norm(matrixvsfra(i,:))~=0)
        indice=find(matrixvsfra(i,:)~=0);
        n=size(indice,2);
        if (n>=1) %先仅考虑n=1时的情况
            selected_face=zeros(1,3);
            the_frac_cell=matrixvsfra(i,indice);%该裂缝单元编号
            %             if Kf(the_frac_cell) <= (kx(i)*ky(i)*kz(i))^(1/3);%
            %             else
            which_frac=fcff{1,2}(the_frac_cell,1);%该裂缝单元所在的裂缝面
            core_frac_cell=fcff{1,4}(the_frac_cell,1:3);%该裂缝单元的中心坐标
            area_frac_cell=fcff{1,8}(the_frac_cell,1);%该裂缝单元面积
            rao_struct=[i,the_frac_cell+nmc];
            pos_rao=find_row(rao_struct,Ninterflow);
            %该裂缝面法向量
            D_value=core_frac_cell-grid_core;% 向量OF
            vector_1=f(5*which_frac-3,:);
            vector_2=f(5*which_frac-2,:);
            n_vector_1=cross(vector_1,vector_2);
            n_vector_2=-cross(vector_1,vector_2);
            if (dot(n_vector_1,D_value))>=0
                n_vector=n_vector_1;
            else 
                n_vector=n_vector_2;
            end
            % 方向投影面积
            area_frac_cell_x=abs(area_frac_cell*dot([1,0,0],n_vector)/norm(n_vector));
            area_frac_cell_y=abs(area_frac_cell*dot([0,1,0],n_vector)/norm(n_vector));
            area_frac_cell_z=abs(area_frac_cell*dot([0,0,1],n_vector)/norm(n_vector));
            %% ----------------------% 根据fractureCellType选取投影面--------------------------------------
            % selected_face(1)=4; nei_grid_x=i+1;
            % selected_face(2)=5; nei_grid_y=i+nx;
            % selected_face(3)=6;nei_grid_z=i+nx*ny;
            %% 选择哪种方法
            selectedMethod = 1;
            %% method 1
            if selectedMethod == 1
            if fractureCellType(the_frac_cell) == 1
            nei_grid_x=i-1;nei_grid_y=i-nx;nei_grid_z=i+nx*ny;
            end
            if fractureCellType(the_frac_cell) == 2
            nei_grid_x=i-1;nei_grid_y=i-nx;nei_grid_z=i+nx*ny;           
            end
            if fractureCellType(the_frac_cell) == 3
            nei_grid_x=i+1;nei_grid_y=i+nx;nei_grid_z=i+nx*ny;
            end
            struc_x=[i,nei_grid_x];
            struc_y=[i,nei_grid_y];
            struc_z=[i,nei_grid_z];
            end
            %% method 2
            if selectedMethod == 2
            if fractureCellType(the_frac_cell) == 1
            nei_grid_x=i-1;nei_grid_y=i+nx;nei_grid_z=i+nx*ny;
            end
            if fractureCellType(the_frac_cell) == 2
            nei_grid_x=i-1;nei_grid_y=i+nx;nei_grid_z=i+nx*ny;
            end
            if fractureCellType(the_frac_cell) == 3
            nei_grid_x=i+1;nei_grid_y=i-nx;nei_grid_z=i+nx*ny;
            end
            struc_x=[i,nei_grid_x];
            struc_y=[i,nei_grid_y];
            struc_z=[i,nei_grid_z];
            end
            %% method 3
            if selectedMethod == 3
            if fractureCellType(the_frac_cell) == 1
            nei_grid_x=i-1;nei_grid_y=i+nx;nei_grid_z=i+nx*ny;
            struc_x=[i,nei_grid_x];
            struc_y=[i,i];% 相当于没有
            struc_z=[i,nei_grid_z];
            end
            if fractureCellType(the_frac_cell) == 2
            nei_grid_x=i-1;nei_grid_y=i+nx;nei_grid_z=i+nx*ny;
            struc_x=[i,nei_grid_x];
            struc_y=[i,nei_grid_y];area_frac_cell_y = area_frac_cell_y*2;% 加上了fractureCellType(the_frac_cell) == 1时的area_frac_cell_y
            struc_z=[i,nei_grid_z];
            end
            if fractureCellType(the_frac_cell) == 3
            nei_grid_x=i+1;nei_grid_y=i-nx;nei_grid_z=i+nx*ny;
            struc_x=[i,nei_grid_x];
            struc_y=[i,nei_grid_y];
            struc_z=[i,nei_grid_z];
            end
            end            
            %% ----------------------% 各方向投影面积叠加及传导率计算---------------------------------------------------------------
            %% x direction
%             struc_x=[i,nei_grid_x];
            pos=find_row(struc_x,N);
            if size(pos,1)==1 %说明需要削减这两个基质网格间的传导系数
                S_p(pos)=S_p(pos)+area_frac_cell_x;
                kcell=(kx(i)*ky(i)*kz(i))^(1/3);%基质网格渗透率
                knei=(kx(nei_grid_x)*ky(nei_grid_x)*kz(nei_grid_x))^(1/3);%相邻那个基质网格渗透率
                added_matnodes=nodes(nei_grid_x,:);
                nei_grid_core=mean(coord(added_matnodes,:));%该基质网格中心
                %                 Knnc=2*kcell*Kf(the_frac_cell)/(kcell+Kf(the_frac_cell));%调和平均后
                %                 dis_added=norm(nei_grid_core-core_frac_cell);
                %                 T_added=Knnc*area_frac_cell_x/dis_added;
                % 采用另一种方式计算几何因子
                dis_added=norm(nei_grid_core-core_frac_cell);
                if Kf(the_frac_cell)==0
                    T_added=0;
                else
                    % 单独裂缝单元几何因子
                    Geof=Kf(the_frac_cell)*area_frac_cell/(Wf(the_frac_cell)/2);
                    % 单独相邻基质网格向裂缝单元几何因子
                    Geofnei=knei*area_frac_cell_x/dis_added;
                    % 单独该基质网格向裂缝单元几何因子
                    Geofcell=Tinterflow(pos_rao,1);
                    
                    T_added=1/(1/Geofnei+1/Geofcell+1/Geof);
                    %                 T_added=0;%特殊情况
                    %                 T_added=0;%特殊情况
                    %                 T_added=0.0039604;%特殊情况
                end
                Nmf_projection=[Nmf_projection;nei_grid_x,the_frac_cell+nmc];
                Tmf_projection=[Tmf_projection;T_added];
            end
            %% y direction
%             struc_y=[i,nei_grid_y];
            pos=find_row(struc_y,N);
            if size(pos,1)==1 %说明需要削减这两个基质网格间的传导系数
                S_p(pos)=S_p(pos)+area_frac_cell_y;
                kcell=(kx(i)*ky(i)*kz(i))^(1/3);%基质网格渗透率
                knei=(kx(nei_grid_y)*ky(nei_grid_y)*kz(nei_grid_y))^(1/3);%相邻那个基质网格渗透率
                added_matnodes=nodes(nei_grid_y,:);
                nei_grid_core=mean(coord(added_matnodes,:));%该基质网格中心
                %                 Knnc=2*kcell*Kf(the_frac_cell)/(kcell+Kf(the_frac_cell));%调和平均后
                %                 dis_added=norm(nei_grid_core-core_frac_cell);
                %                 T_added=Knnc*area_frac_cell_y/dis_added;
                % 采用另一种方式计算几何因子
                dis_added=norm(nei_grid_core-core_frac_cell);
                if Kf(the_frac_cell)==0
                    T_added=0;
                else
                    % 单独裂缝单元几何因子
                    Geof=Kf(the_frac_cell)*area_frac_cell/(Wf(the_frac_cell)/2);
                    % 单独相邻基质网格向裂缝单元几何因子
                    Geofnei=knei*area_frac_cell_y/dis_added;
                    % 单独该基质网格向裂缝单元几何因子
                    Geofcell=Tinterflow(pos_rao,1);
                    
                    T_added=1/(1/Geofnei+1/Geofcell+1/Geof);
                    %                 T_added=0;%特殊情况
                    %                 T_added=0.002;%特殊情况
                    %                 T_added=0.0039604;%特殊情况
                end
                Nmf_projection=[Nmf_projection;nei_grid_y,the_frac_cell+nmc];
                Tmf_projection=[Tmf_projection;T_added];
            end
            
            %% z direction
%             struc_z=[i,nei_grid_z];
            pos=find_row(struc_z,N);
            if size(pos,1)==1 %说明需要削减这两个基质网格间的传导系数
                S_p(pos)=S_p(pos)+area_frac_cell_z;
                kcell=(kx(i)*ky(i)*kz(i))^(1/3);%基质网格渗透率
                knei=(kx(nei_grid_z)*ky(nei_grid_z)*kz(nei_grid_z))^(1/3);%相邻那个基质网格渗透率
                added_matnodes=nodes(nei_grid_z,:);
                nei_grid_core=mean(coord(added_matnodes,:));%该基质网格中心
                %                 Knnc=2*kcell*Kf(the_frac_cell)/(kcell+Kf(the_frac_cell));%调和平均后
                %                 dis_added=norm(nei_grid_core-core_frac_cell);
                %                 T_added=Knnc*area_frac_cell_z/dis_added;
                % 采用另一种方式计算几何因子
                dis_added=norm(nei_grid_core-core_frac_cell);
                if Kf(the_frac_cell)==0
                    T_added=0;
                else
                    % 单独裂缝单元几何因子
                    Geof=Kf(the_frac_cell)*area_frac_cell/(Wf(the_frac_cell)/2);
                    % 单独相邻基质网格向裂缝单元几何因子
                    Geofnei=knei*area_frac_cell_z/dis_added;
                    % 单独该基质网格向裂缝单元几何因子
                    Geofcell=Tinterflow(pos_rao,1);
                    
                    T_added=1/(1/Geofnei+1/Geofcell+1/Geof);
                    %                 T_added=0;%特殊情况
                    %                 T_added=0.0;%特殊情况
                    %                 T_added=0.0039604;%特殊情况
                end
                Nmf_projection=[Nmf_projection;nei_grid_z,the_frac_cell+nmc];
                Tmf_projection=[Tmf_projection;T_added];
            end
            
        end
    end
end
%% 额外添加的改进
im_N=[];
im_T=[];
% % % q=length(Kf)+1;
% % % for i=1:(q-1)/2%q-1是裂缝单元数量
% % %     for j=(q-1)/2+1:(q-1)
% % %         if abs(corevsfra(j,2)-corevsfra(i,2))<=1e-2
% % %             im_N=[im_N;i+nmc,j+nmc];
% % % %             im_T=[im_T;0];
% % % %             im_T=[im_T;0.0039604];
% % % %               im_T=[im_T;(((Kf(i))^(-1)+(1e-3)^(-1)+(1e-3)^(-1)+(Kf(j))^(-1))/4)^(-1)*2*2/0.4];
% % %               im_T=[im_T;(((Kf(i))^(-1)+(1e-3)^(-1)+(1e-3)^(-1)+(Kf(j))^(-1))/4)^(-1)*0.5*2/0.4];
% % %             %
% % %         end
% % %     end
% % % end
%% 用投影面积修正基质网格间传导
for i=1:nmm
    % 为避免数值误差，导致没有真正封堵住，故
    if (S_ex(i)-S_p(i))/S_ex(i)<=1e-2
        T(i)=0;
    else
        T(i)= T(i)*(S_ex(i)-S_p(i))/S_ex(i);
    end
end
%% 实用型pEDFM额外处理
% % % num_frac = size(frac_information,1)/3;
% % % for i = 1:num_frac
% % %     if frac_information(3*i,1) <= (kx(1)*ky(1)*kz(1))^(1/3)
% % %         % 裂缝的两端点
% % %         one_end = frac_information(3*i-2,:);
% % %         the_other_end = frac_information(3*i-1,:);
% % %         % 两网格之间的连接
% % %         for j = 1:size(N,1)
% % %         added_matnodes_1=nodes(N(j,1),:);
% % %         one_node_connect=mean(coord(added_matnodes_1,:));%该基质网格中心
% % %         added_matnodes_2=nodes(N(j,2),:);
% % %         the_other_node_connect=mean(coord(added_matnodes_2,:));
% % %         coef_matrix = [
% % %             the_other_end(1)-one_end(1),-(the_other_node_connect(1)-one_node_connect(1));
% % %             the_other_end(2)-one_end(2),-(the_other_node_connect(2)-one_node_connect(2));];
% % %         coef_vector = [one_node_connect(1)-one_end(1);one_node_connect(2)-one_end(2);];
% % %         if abs(det(coef_matrix))<=1e-12
% % %         else
% % %             parameters = coef_matrix\coef_vector;
% % %             if parameters(1)>=0 && parameters(1)<=1 && parameters(2)>=0 && parameters(2)<=1
% % %                 T(j) = (T(j)^-1 + frac_information(3*i,1)^-1)^-1;
% % %             end
% % %         end
% % %         end
% % %     end
% % % end
%% 把上述三种情况的矩阵分别叠加起来
N = [N; Nmff+nmc; Nff+nmc;Ninterflow;Nmf_projection;im_N];
T = [T; Tmff; Tff;Tinterflow;Tmf_projection;im_T];
% N = [N; Nmff+nmc; Nff+nmc;Ninterflow;Nmf_projection];
% T = [T; Tmff; Tff;Tinterflow;Tmf_projection];
nex=size(N,1);
end