function [ matrixvsfra,connectmf,fracnumber,fracstart,corevsfra,connect_infrac,N,T,zf,fcinff,vf,porf,fcff,T_convection,mat_frac,number_of_adding,cell_divided_by_fracture_flag,volume_ratio, fractureCell_matrixCell_addedMatrixCell] = connections_unstructured_EDFM_EX2_1to4_solutions( fracturemesh,pointvsregion,raopoint,f,nf,nx,ny,nz,kx,ky,kz,kf,wf,Porf,dxv,dyv,dzv,coord,nodes,NTG,addflag )
% 该函数用以给裂缝面网格编号，并存储网格间的连接关系及传导率，是前处理的关键步骤
% 不同于二维嵌入式离散裂缝模型，三维模型中，裂缝是二维平面，网格之间连接关系复杂，
%不能简单得从基质网格来判定，若所处基质网格不相邻，则裂缝网格必不相邻，若所处基质网格相邻，裂缝网格却不一定相邻
%上述原则无助于我们去确定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^nf)-1)/(2^nf));
        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);
mat_frac=[];
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;
        mat_frac=[mat_frac;i];
        %给这些裂缝点/裂缝单元按顺序编号
        A=matrixvsfra(i,:);
        A(A==0)=[];
        connectmf{cc,2} =A; 
end
end

%% 一条条裂缝计算裂缝网格面积、边界等等
% edgevsfra=cell(nf,1);
% Kffl=funtion_data.Kffl;
% Dffl=funtion_data.Dffl;
% kffl=zeros(q-1,1);
% dffl=zeros(q-1,1);
fractureCellFlag = (1:q-1)';
normalvec=zeros(q-1,3);%unit 法向量
areavsfra=zeros(q-1,1);%q-1即裂缝网格总数
lengthvsfra=zeros(q-1,10);
Kf=zeros(q-1,10);%裂缝网格渗透率
corevsfra=zeros(q-1,3);
disvsfra=zeros(q-1,10);
tran=zeros(q-1,10);
avtran=zeros(q-1,1);%裂缝网格平均传导，用于后续的简化计算
ave_disvsfra=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);
        Kf(i,1)=kf(k-1);
%         kffl(i,1)=Kffl(k-1);
%         dffl(i,1)=Dffl(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);
        ave_disvsfra(i)=sum(disvsfra(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); %记录相应的传导系数
Tff_convection=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(i)==0 hartran=0;hartran_convection=0;
            else
            tran1=kf(i)*wf(i)*lengthcat/disk;
            tran2=kf(i)*wf(i)*lengthcat/disj;
%             hartran=length/(disk+disj);%取调和平均
            hartran=tran1*tran2/(tran1+tran2);
            k_harmony= kf(i);
            hartran_convection=kf(i)/(disk+disj)^2;
            end
            Nff=[Nff; j,k]; 
            Tff=[Tff;hartran];
            Tff_convection=[Tff_convection;hartran_convection];
            %计算出每个裂缝网格重心到每条边的距离
        end
        end
        p=p+1;
        end
    connect_infrac{i,1}=raoconnect;
    end
end
Nff(1,:)=[]; %去除第一行的零行
Tff(1,:)=[];
Tff_convection(1,:)=[];
%% 在同一基质网格中的裂缝单元连接情况及传导率计算
nmc=size(matrixvsfra,1);
Nmff=zeros(1,2);
Tmff=zeros(1,1);
Tmff_convection=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];Tmff_convection=[Tmff_convection;0];
                else
                Tmff=[Tmff; fcff{1,7}(fc1)*fcff{1,7}(fc2)/sumtran];
                end
            end
        end
        end
    end
end
Nmff(1,:)=[]; %去掉开头的零行
Tmff(1,:)=[];
Tmff_convection(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是裂缝单元之间存在流体交换的总数
Kx=kx.*ones(nmc,1);
Ky=ky.*ones(nmc,1);
Kz=kz.*ones(nmc,1);
N = zeros(nmm, 2);%存储基质网格之间存在流体交换的网格编号
T = zeros(nmm, 1);%存储对应与N矩阵的传导系数
T_convection = zeros(nmm, 1);
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));
        T_convection(c) = 2/(1/Kx(index)+1/Kx(indexn))/(dxv(indexn)+dxv(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));
        T_convection(c) = 2/(1/Ky(index)+1/Ky(indexn))/(dyv(indexn)+dyv(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));
        T_convection(c) = 2/(1/Kz(index)+1/Kz(indexn))/(dzv(indexn)+dzv(index));
    end
end
end
%% 按照2014年的方法，将窜流处理为与上述相似的形式
syms x y z;
Ninterflow=[];
Tinterflow=[];
Tinterflow_convection=[];
matrix_cells_containing_fracture_cells = [];
volume_ratio = [];
fractureMatrixFlag = zeros(size(fractureCellFlag,1),2);
%% -------------------------------------------------------------------------此处是修改之重---------------------------------------------------------------
numberInGroup = 4;
%% -------------------------------------------------------------------------上方修改结束---------------------------------------------------------------
for i=1:nx*ny*nz  %基质网格编号
    for j=1:10
        if matrixvsfra(i,j)~=0  %this matrix cell contains a fracture cell
            fractureMatrixFlag(matrixvsfra(i,j),:) = [matrixvsfra(i,j),i];
            matrix_cells_containing_fracture_cells = [matrix_cells_containing_fracture_cells; i];
%% -------------------------------------------------------------------------此处是修改之重---------------------------------------------------------------
            if mod(matrixvsfra(i,j),numberInGroup) == 1 volume_ratio = [volume_ratio; 7/8];end % 右侧重新编号网格所占体积比例
            if mod(matrixvsfra(i,j),numberInGroup) == 2 volume_ratio = [volume_ratio; 5/8];end
            if mod(matrixvsfra(i,j),numberInGroup) == 3 volume_ratio = [volume_ratio; 3/8];end
            if mod(matrixvsfra(i,j),numberInGroup) == 0 volume_ratio = [volume_ratio; 1/8];end
%% -------------------------------------------------------------------------上方修改结束---------------------------------------------------------------
            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=2*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
    raoT=Knnc*Annc/Dn;
    raoT_convection=2*Knnc/Dn^2;
    Tinterflow=[Tinterflow;raoT];
    Tinterflow_convection=[Tinterflow_convection;raoT_convection];
        end
    end
end
%% 额外处理，此处并未通用化，有待进一步开展
nfc=size(fracnumber,1);
the_added_matrix_cells = [];
nmc=nx * ny*nz;
the_added_mf_connec = [];
the_added_mf_T = [];
number_of_adding = 0;
cell_divided_by_fracture_flag = matrix_cells_containing_fracture_cells;
% 记录基质网格被分成了哪两个网格,并删除原有的与这些基质网格相关的连接
matrixCellsDivided = zeros(size(matrix_cells_containing_fracture_cells,1),2);
deletedFlag = [];
for i = 1:size(matrix_cells_containing_fracture_cells,1)
    the_matrix_cell = matrix_cells_containing_fracture_cells(i,1);
    number_of_adding = number_of_adding+1;
    the_new_cell = number_of_adding+nmc;
    the_added_matrix_cells = [the_added_matrix_cells; the_new_cell];
    matrixCellsDivided(i,:) = [the_matrix_cell, the_new_cell];
    % m-m连接
    for j = 1:size(N,1)
        if N(j,1) == the_matrix_cell || N(j,2) == the_matrix_cell
            deletedFlag = [deletedFlag; j];
        end
% % %             if N(j,2) == the_matrix_cell+1 
% % %                 N(j,1) = the_new_cell; T(j) = ...
% % %                     ((Kx(the_matrix_cell)*dzv(the_matrix_cell)*dyv(the_matrix_cell)/(dxv(the_matrix_cell)/3))^(-1)+(Kx(the_matrix_cell+1)*dzv(the_matrix_cell+1)*dyv(the_matrix_cell+1)/(dxv(the_matrix_cell+1)/2))^(-1))^(-1);
% % %             end
            
% % %             if N(j,2) == the_matrix_cell-nx
% % %                 N(j,1) = the_new_cell; T(j) = ...
% % %                      ((Ky(the_matrix_cell)*dzv(the_matrix_cell)*dxv(the_matrix_cell)/(dyv(the_matrix_cell)/3))^(-1)+(Ky(the_matrix_cell-n)*dzv(the_matrix_cell-n)*dxv(the_matrix_cell-n)/(dxv(the_matrix_cell-n)/2))^(-1))^(-1);
% % %             end
        
% % %             if N(j,1) == the_matrix_cell+1 
% % %                 N(j,2) = the_new_cell; T(j) = ...
% % %                     ((Kx(the_matrix_cell)*dzv(the_matrix_cell)*dyv(the_matrix_cell)/(dxv(the_matrix_cell)/3))^(-1)+(Kx(the_matrix_cell+1)*dzv(the_matrix_cell+1)*dyv(the_matrix_cell+1)/(dxv(the_matrix_cell+1)/2))^(-1))^(-1);
% % %             end
     
% % %             if N(j,1) == the_matrix_cell-nx
% % %                 N(j,2) = the_new_cell; T(j) = ...
% % %                      ((Ky(the_matrix_cell)*dzv(the_matrix_cell)*dxv(the_matrix_cell)/(dyv(the_matrix_cell)/3))^(-1)+(Ky(the_matrix_cell-n)*dzv(the_matrix_cell-n)*dxv(the_matrix_cell-n)/(dxv(the_matrix_cell-n)/2))^(-1))^(-1);
% % %             end
       
    end
end
deletedFlag = unique(deletedFlag);
N(deletedFlag,:) = [];
T(deletedFlag,:) = [];
% 构建裂缝网格-基质网格-添加基质网格的对应矩阵，按照裂缝网格顺序排列，方便后续处理
fractureCell_matrixCell_addedMatrixCell = zeros(size(fractureCellFlag,1),3);
for xiang = 1:size(fractureCellFlag,1)
    theFractureCell = fractureCellFlag(xiang,1);
    aaa = find(fractureMatrixFlag(:,1) == theFractureCell);
    theMatrixCell = fractureMatrixFlag(aaa,2);
    bbb = find(matrixCellsDivided(:,1) == theMatrixCell);
    theAddedMatrixCells = matrixCellsDivided(bbb,2);
    fractureCell_matrixCell_addedMatrixCell(theFractureCell,:) = [theFractureCell, theMatrixCell, theAddedMatrixCells];
end
% 添加这些基质网格的真实连接
relatedN = [];
relatedT = [];
for xiang = 1:size(fractureCellFlag,1)/numberInGroup
%% -------------------------------------------------------------------------此处亦是修改之重-----------------------------------------------------------------------
        addedNSs =[
        % 左侧原有基质网格 x方向上的连接
        fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-3,2)-1, fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-3,2);
        fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-2,2)-1, fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-2,2);
        fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-1,2)-1, fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-1,2);
        fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-0,2)-1, fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-0,2);
        % 右侧添加基质网格 x方向上的连接
        fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-3,3), fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-3,2)+1;
        fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-2,3), fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-2,2)+1;
        fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-1,3), fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-1,2)+1;
        fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-0,3), fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-0,2)+1; 
        % 左侧原有基质网格 y方向上的连接
        fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-3,2), fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-2,2);
        fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-2,2), fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-1,2);
        fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-1,2), fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-0,2);
        fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-0,2), fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-0,2)+nx;
        % 右侧添加基质网格 y方向上的连接
        fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-3,2)-nx, fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-3,3);
        fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-3,3), fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-2,3);
        fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-2,3), fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-1,3);
        fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-1,3), fractureCell_matrixCell_addedMatrixCell(numberInGroup*xiang-0,3);
        ];
%% --------------------------------------------------------------此处修改结束-----------------------------------------------------------------
        relatedN = [relatedN;
            addedNSs
                 ];
        % 计算传导率
             addedTs = [];
             for dai = 1:size(addedNSs,1)
                 the_matrix_cell1 = addedNSs(dai,1); the_matrix_cell2 = addedNSs(dai,2);
                 if the_matrix_cell1> nmc
                     aaa = find(fractureCell_matrixCell_addedMatrixCell(:,3) == the_matrix_cell1);
                     the_matrix_cell1 = fractureCell_matrixCell_addedMatrixCell(aaa,2);
                 end
                 if the_matrix_cell2> nmc
                     aaa = find(fractureCell_matrixCell_addedMatrixCell(:,3) == the_matrix_cell2);
                     the_matrix_cell2 = fractureCell_matrixCell_addedMatrixCell(aaa,2);
                 end      
%% -------------------------------------------------------------------------此处亦是修改之重-----------------------------------------------------------------------
                 %
                 if dai == 1
                     theT = ((Kx(the_matrix_cell1)*dzv(the_matrix_cell1)*dyv(the_matrix_cell1)/(dxv(the_matrix_cell1)/2))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dyv(the_matrix_cell2)/(dxv(the_matrix_cell2)*1/16))^(-1))^(-1);
                 end
                 if dai == 2
                     theT = ((Kx(the_matrix_cell1)*dzv(the_matrix_cell1)*dyv(the_matrix_cell1)/(dxv(the_matrix_cell1)/2))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dyv(the_matrix_cell2)/(dxv(the_matrix_cell2)*3/16))^(-1))^(-1);
                 end
                 if dai == 3
                     theT = ((Kx(the_matrix_cell1)*dzv(the_matrix_cell1)*dyv(the_matrix_cell1)/(dxv(the_matrix_cell1)/2))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dyv(the_matrix_cell2)/(dxv(the_matrix_cell2)*5/16))^(-1))^(-1);
                 end
                 if dai == 4
                     theT = ((Kx(the_matrix_cell1)*dzv(the_matrix_cell1)*dyv(the_matrix_cell1)/(dxv(the_matrix_cell1)/2))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dyv(the_matrix_cell2)/(dxv(the_matrix_cell2)*7/16))^(-1))^(-1);
                 end
                 %
                 if dai == 5
                     theT = ((Kx(the_matrix_cell1)*dzv(the_matrix_cell1)*dyv(the_matrix_cell1)/(dxv(the_matrix_cell1)*7/16))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dyv(the_matrix_cell2)/(dxv(the_matrix_cell2)/2))^(-1))^(-1);
                 end
                 if dai == 6
                     theT = ((Kx(the_matrix_cell1)*dzv(the_matrix_cell1)*dyv(the_matrix_cell1)/(dxv(the_matrix_cell1)*5/16))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dyv(the_matrix_cell2)/(dxv(the_matrix_cell2)/2))^(-1))^(-1);
                 end
                 if dai == 7
                     theT = ((Kx(the_matrix_cell1)*dzv(the_matrix_cell1)*dyv(the_matrix_cell1)/(dxv(the_matrix_cell1)*3/16))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dyv(the_matrix_cell2)/(dxv(the_matrix_cell2)/2))^(-1))^(-1);
                 end
                 if dai == 8
                     theT = ((Kx(the_matrix_cell1)*dzv(the_matrix_cell1)*dyv(the_matrix_cell1)/(dxv(the_matrix_cell1)*1/16))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dyv(the_matrix_cell2)/(dxv(the_matrix_cell2)/2))^(-1))^(-1);
                 end
                 %
                 if dai == 9
                     theT = ((Ky(the_matrix_cell1)*dzv(the_matrix_cell1)*dxv(the_matrix_cell1)*1/4/(dyv(the_matrix_cell1)/2))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dxv(the_matrix_cell2)*1/4/(dyv(the_matrix_cell2)/2))^(-1))^(-1);
                 end
                 if dai == 10
                     theT = ((Ky(the_matrix_cell1)*dzv(the_matrix_cell1)*dxv(the_matrix_cell1)*2/4/(dyv(the_matrix_cell1)/2))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dxv(the_matrix_cell2)*2/4/(dyv(the_matrix_cell2)/2))^(-1))^(-1);
                 end
                 if dai == 11
                     theT = ((Ky(the_matrix_cell1)*dzv(the_matrix_cell1)*dxv(the_matrix_cell1)*3/4/(dyv(the_matrix_cell1)/2))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dxv(the_matrix_cell2)*3/4/(dyv(the_matrix_cell2)/2))^(-1))^(-1);
                 end
                 if dai == 12
                     theT = ((Ky(the_matrix_cell1)*dzv(the_matrix_cell1)*dxv(the_matrix_cell1)*4/4/(dyv(the_matrix_cell1)/2))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dxv(the_matrix_cell2)*4/4/(dyv(the_matrix_cell2)/2))^(-1))^(-1);
                 end
                 % 
                 if dai == 13
                     theT = ((Ky(the_matrix_cell1)*dzv(the_matrix_cell1)*dxv(the_matrix_cell1)*4/4/(dyv(the_matrix_cell1)/2))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dxv(the_matrix_cell2)*4/4/(dyv(the_matrix_cell2)/2))^(-1))^(-1);
                 end
                 if dai == 14
                     theT = ((Ky(the_matrix_cell1)*dzv(the_matrix_cell1)*dxv(the_matrix_cell1)*3/4/(dyv(the_matrix_cell1)/2))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dxv(the_matrix_cell2)*3/4/(dyv(the_matrix_cell2)/2))^(-1))^(-1);
                 end
                 if dai == 15
                     theT = ((Ky(the_matrix_cell1)*dzv(the_matrix_cell1)*dxv(the_matrix_cell1)*2/4/(dyv(the_matrix_cell1)/2))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dxv(the_matrix_cell2)*2/4/(dyv(the_matrix_cell2)/2))^(-1))^(-1);
                 end
                 if dai == 16
                     theT = ((Ky(the_matrix_cell1)*dzv(the_matrix_cell1)*dxv(the_matrix_cell1)*1/4/(dyv(the_matrix_cell1)/2))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dxv(the_matrix_cell2)*1/4/(dyv(the_matrix_cell2)/2))^(-1))^(-1);
                 end
%% -------------------------------------------------------------------------此处修改结束-----------------------------------------------------------------------
% % %                      ((Kx(the_matrix_cell1)*dzv(the_matrix_cell1)*dyv(the_matrix_cell1)/(dxv(the_matrix_cell1)*5/12))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dyv(the_matrix_cell2)/(dxv(the_matrix_cell2)/2))^(-1))^(-1);
% % %                      ((Kx(the_matrix_cell1)*dzv(the_matrix_cell1)*dyv(the_matrix_cell1)/(dxv(the_matrix_cell1)*1/4))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dyv(the_matrix_cell2)/(dxv(the_matrix_cell2)/2))^(-1))^(-1);
% % %                      ((Kx(the_matrix_cell1)*dzv(the_matrix_cell1)*dyv(the_matrix_cell1)/(dxv(the_matrix_cell1)*1/6))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dyv(the_matrix_cell2)/(dxv(the_matrix_cell2)/2))^(-1))^(-1);
% % %                      
% % %                      ((Ky(the_matrix_cell1)*dzv(the_matrix_cell1)*dxv(the_matrix_cell1)*1/3/(dyv(the_matrix_cell1)/2))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dxv(the_matrix_cell2)*1/3/(dyv(the_matrix_cell2)/2))^(-1))^(-1);
% % %                      ((Ky(the_matrix_cell1)*dzv(the_matrix_cell1)*dxv(the_matrix_cell1)*2/3/(dyv(the_matrix_cell1)/2))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dxv(the_matrix_cell2)*2/3/(dyv(the_matrix_cell2)/2))^(-1))^(-1);
% % %                      ((Ky(the_matrix_cell1)*dzv(the_matrix_cell1)*dxv(the_matrix_cell1)*3/3/(dyv(the_matrix_cell1)/2))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dxv(the_matrix_cell2)*3/3/(dyv(the_matrix_cell2)/2))^(-1))^(-1);
% % %                      
% % %                      ((Ky(the_matrix_cell1)*dzv(the_matrix_cell1)*dxv(the_matrix_cell1)*3/3/(dyv(the_matrix_cell1)/2))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dxv(the_matrix_cell2)*3/3/(dyv(the_matrix_cell2)/2))^(-1))^(-1);
% % %                      ((Ky(the_matrix_cell1)*dzv(the_matrix_cell1)*dxv(the_matrix_cell1)*2/3/(dyv(the_matrix_cell1)/2))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dxv(the_matrix_cell2)*2/3/(dyv(the_matrix_cell2)/2))^(-1))^(-1);
% % %                      ((Ky(the_matrix_cell1)*dzv(the_matrix_cell1)*dxv(the_matrix_cell1)*1/3/(dyv(the_matrix_cell1)/2))^(-1)+(Kx(the_matrix_cell2)*dzv(the_matrix_cell2)*dxv(the_matrix_cell2)*1/3/(dyv(the_matrix_cell2)/2))^(-1))^(-1);
                     
                 addedTs = [addedTs; theT];   
             end
        relatedT = [relatedT;addedTs];
end
N = [N; relatedN];
T = [T; relatedT];
%
Ninterflow(:,2) = Ninterflow(:,2)+ number_of_adding;
ori_nmc = nmc;
nmc = ori_nmc +number_of_adding;
% m-f连接
% number_of_adding = 0;
for i = 1:size(matrix_cells_containing_fracture_cells,1)
    the_matrix_cell = matrix_cells_containing_fracture_cells(i,1);
    aaa = find(fractureCell_matrixCell_addedMatrixCell(:,2) == the_matrix_cell);
    the_new_cell = fractureCell_matrixCell_addedMatrixCell(aaa,3);
    the_fracture_cell_flag = fractureCell_matrixCell_addedMatrixCell(aaa,1);
%     number_of_adding = number_of_adding+1;
%     the_new_cell = number_of_adding+ori_nmc;
for j = 1:size(Ninterflow,1)
      if Ninterflow(j,1) == the_matrix_cell
            the_added_mf_connec = [the_added_mf_connec;the_new_cell,Ninterflow(j,2)];
%% -------------------------------------------------------------------------此处亦是修改之重-----------------------------------------------------------------------
            if mod(the_fracture_cell_flag, numberInGroup) == 1
            the_added_mf_T = [the_added_mf_T; Tinterflow(j)/2*((1/4)/((1+3/4)/4))];
            Tinterflow(j) = Tinterflow(j)/2*((1/4)/((2-1-3/4)/4));
            end
            if mod(the_fracture_cell_flag, numberInGroup) == 2
            the_added_mf_T = [the_added_mf_T; Tinterflow(j)/2*((1/4)/((3/4+2/4)/4))];
            Tinterflow(j) = Tinterflow(j)/2*((1/4)/((2-3/4-2/4)/4));
            end
            if mod(the_fracture_cell_flag, numberInGroup) == 3
            the_added_mf_T = [the_added_mf_T; Tinterflow(j)/2*((1/4)/((2/4+1/4)/4))];
            Tinterflow(j) = Tinterflow(j)/2*((1/4)/((2-2/4-1/4)/4));
            end
            if mod(the_fracture_cell_flag, numberInGroup) == 0
            the_added_mf_T = [the_added_mf_T; Tinterflow(j)/2*((1/4)/((0+1/4)/4))];
            Tinterflow(j) = Tinterflow(j)/2*((1/4)/((2-0-1/4)/4));
            end
%% -------------------------------------------------------------------------此处亦是修改之重-----------------------------------------------------------------------
      end
end
end
Ninterflow = [Ninterflow; the_added_mf_connec];
Tinterflow = [Tinterflow; the_added_mf_T];
%% 按照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);
%             K1=Kf(matrixvsfra(i,j));d1=Wf(matrixvsfra(i,j));
%             K2=kffl(matrixvsfra(i,j));d2=dffl(matrixvsfra(i,j));
%             Annc=2*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
%     raoT=((Kcell*Annc/(Dn-d2))^(-1)+(K2*Annc/d2)^(-1)+(K1*Annc/d1)^(-1))^(-1);
%     Tinterflow=[Tinterflow;raoT];
%         end
%     end
% end
            
%% 把上述三种情况的矩阵分别叠加起来
N = [N; Nmff+nmc; Nff+nmc;Ninterflow;];
T = [T; Tmff; Tff;Tinterflow;];
T_convection = [T_convection;Tmff_convection; Tff_convection;Tinterflow_convection;];
nex=size(N,1);
end