%算例1：混凝土立方块
% CSPD (Peridynamics based on Cosserat continuum) 
clc;clear
%% 创建计算结果保存路径
directpath='G:\fatigue_work001\';
if (~exist(directpath,'dir') )
    mkdir(directpath)        % 若不存在，在当前目录中产生一个子目录
end
%% 参数输入
len=100.0d-3; %原模型横向长度
wid=100.0d-3; %原模型竖向长度
hei=100.0d-3; %原模型高度
xdv=25;  %x向划分
ydv=25;  %y向划分
zdv=25; 
dl=len/xdv;  %节点之间水平或者竖直方向间距
dt=1.0d0;%虚拟时间步长
tao=1.0d0;
nt=20000*tao;%总时间步数
allnump=xdv*ydv;  %节点总数
dotp=zeros(allnump,2);  %将初始节点坐标（默认为零）存入dotp矩阵
n=0; %节点编号初始化
%材料参数
thick =100.0d-3;%厚度
dens = 2400.0;%密度
e1 = 24.6d9;   %弹性模量
emod2 = 200.0d0*e1;
pratio = 0.2d0;   %泊松比
G=e1/(2*(1+pratio));%剪切模量
delta=3.015*dl;    %近场范围尺寸
r=0.2353d0;
einf=0;
Gc = 0.5d0*G;
value=1.0d0;
Gc2=Gc*200;
G2=G*200;
lc= 5.0d-5;
crackxl=zeros(nt,1);

admg=zeros(nt,1);
nc=zeros(nt,1);
%迭代参数
a1=0.0611;
m1=19.53;
a2=1.0;
m2=3.1316;
%Cosserat参数
kn=6*e1/(pi*(delta.^4)*(1-pratio));
kt=2*e1*(1-3*pratio)/(pi*(delta^2)*(1-pratio^2));
ld=sqrt(Gc*(1-pratio^2)/(2*e1*(1-3*pratio)))*delta;
kr=4*Gc/(pi*delta^3)-16*G*(lc^2)/(pi*delta^5);

fc=24.6d6;  %混凝土抗压强度
ft=2.2d6;
fs=3.9d6;

%输入键的拉伸和压缩时的参数
GF=120; %能量释放率
Gf=GF/2.5; %参考文献中的定义
k2=ft^2/(2*Gf);
scr1t=sqrt((5.0*GF)/(6*e1*delta));%临界的变形,超过这个值就会破坏
fmaxt=scr1t*kn;
fmaxn=fmaxt*fc/ft;
f0=[fmaxt,fmaxn]; %第一个元素表示受拉，第二个元素表示受压
fmaxs=fmaxt*fs/ft;
scr2t=scr1t+0.75*fmaxt/k2;
bta=0.25;
gma=(3+2*bta)/(2*bta*(1-bta));
scr3t=scr1t+(10*gma*GF/(pi*delta^5*kn*scr1t*(1+gma*bta)));

scr1c=scr1t*(fc/ft);
s0c=scr1c*(-35)/8.94;
smc=scr1c*14/8.94;
scr2c=scr1c*23.23/8.94;
scr3c=scr1c*72/8.94;

%% 下面将节点横纵坐标值按编号对应存入dot矩阵
% 先建立平面
for i=1:xdv
    for j=1:ydv
        n=n+1;
        dotp(n,1)=(-1.0d0*len./2.0d0)+(j-1).*dl;
        dotp(n,2)=(-1.0d0*wid./2.0d0)+(i-1).*dl;
    end
end
allnum=allnump*zdv;
dot=zeros(allnum,3);
%按从左下至右上，从低到高的顺序存入坐标
n=0;
for j=1:zdv
  for i=1:allnump
      n=i+(j-1)*allnump;
      dot(n,1)=dotp(i,1);
      dot(n,2)=dotp(i,2);
      dot(n,3)=(j-1).*dl;
  end
end
% 加入加载块 分左右部和上部
m=0;
%左右加载行
dwids=100.0d-3;
dlens=3.0d0*dl;
dheis=100.0d-3;
xbdvs=round(dlens/dl);
ybdvs=round(dwids/dl);
zbdvs=round(dheis/dl);
dallnum=xbdvs*ybdvs;
allnump=xbdvs*ybdvs*zbdvs;
dotlef1=zeros(dallnum,2);
dotrig1=zeros(dallnum,2);
dotlef=zeros(allnump,3);
dotrig=zeros(allnump,3);
%下面将节点横纵坐标值按编号对应存入dotlef和dotrig矩阵
m=0;
for i=1:3
    for j=1:ybdvs
           m=m+1;
           dotlef1(m,1)=min(dot(:,1))-dlens+(i-1).*dl;
           dotlef1(m,2)=(-1.0d0*len./2.0d0)+(j-1).*dl;
    end
end
m=0;
for i=1:3
    for j=1:ybdvs
           m=m+1;
           dotrig1(m,1)=max(dot(:,1))+dl+(i-1).*dl;
           dotrig1(m,2)=(-1.0d0*len./2.0d0)+(j-1).*dl;
    end
end
%左
n=0;
for j=1:zdv
  for i=1:dallnum
      n=i+(j-1)*dallnum; 
      dotlef(n,1)=dotlef1(i,1);
      dotlef(n,2)=dotlef1(i,2);
      dotlef(n,3)=(j-1).*dl;
  end
end
dotlefn=allnum+n;
%右
n=0;
for j=1:zdv
  for i=1:dallnum
      n=i+(j-1)*dallnum; 
      dotrig(n,1)=dotrig1(i,1);
      dotrig(n,2)=dotrig1(i,2);
      dotrig(n,3)=(j-1).*dl;
  end
end
dotrign=dotlefn+n;
m=0;
%上下加载行
dzdv=3;   %加载板的z向划分
allnumd=xdv*ydv;
dotup=zeros(allnump,3);
dotdown=zeros(allnump,3);
%下面将节点横纵坐标值按编号对应存入dotup和dotdown矩阵
% 上加载板
n=0;
for j=1:dzdv
  for i=1:allnumd
      n=i+(j-1)*allnumd; 
      dotup(n,1)=dot(i,1);
      dotup(n,2)=dot(i,2);
      dotup(n,3)=max(dot(:,3))+j*dl;
  end
end
dotupn=dotrign+n;
% 下加载板
n=0;
for j=1:dzdv
  for i=1:allnumd
      n=i+(j-1)*allnumd;    
      dotdown(n,1)=dot(i,1);
      dotdown(n,2)=dot(i,2);
      dotdown(n,3)=min(dot(:,3))-j*dl;
  end
end
dotdownn=dotupn+n;
dotall=[dot;dotlef;dotrig;dotup;dotdown];  %将所有所要的节点存入一个矩阵中
dotn=1:1:allnum; %内点编号
%% 观察整体构型 不用时可以注释掉
% k1=ones(numel(1:dotdownn),1);
% scatter_color3d(dotall(1:dotdownn,:),k1)   
%% 其他参数初始化
area = dl * dl;%面积
vol = dl * dl * dl;%体积
% fk=vol*(bc1*scr1+bc2*(scr2-scr1));
fk=3.3d7/1.5;
fks=7.9d7/1.5;
idist = 0.0d0;%初始键长
fac = 0.0d0;%体积修正系数
radij = dl / 2.0d0;%点的半径
cnode = 0;%当前的点;
nlength  = 0.0d0;%两点间长度在变形后的长度
dforce1 = 0.0d0;%两点间的近场力分量x向
dforce2 = 0.0d0;%两点间的近场力分量y向
dforce3 = 0.0d0;%两点间的近场力分量z向
dforce4 = 0.0d0;%第一微转动角
dforce5 = 0.0d0;%第二微转动角
dforce6 = 0.0d0;%第三微转动角

nodefam=knn3D(dotall,delta,dl);%总的近场范围点的位置矩阵构造,用的是knn函数
%第一列是统计一共有多少个近场点的，所以循环都应该从第二个点开始
nhorizon=max(nodefam(:,1));%读取nodefam中的列数
fncst=ones(dotdownn,2);%表面修正系数
fail=zeros(dotdownn,nhorizon);%损伤  1代表没破坏，0代表破坏了
lamda=zeros(dotdownn,nhorizon);%键的剩余寿命 1代表完好，0代表破坏
flag=zeros(dotdownn,nhorizon); %判断每一根键的阶段 1代表1阶段，2代表2阶段
eps=zeros(dotdownn,nhorizon); %每根键的应变
bondforce=zeros(dotdownn,nhorizon); %每根键的法向作用
shearforce=zeros(dotdownn,nhorizon);%每根键的切向作用
switcht=zeros(dotdownn,nhorizon); %判定过程开关,0表示疲劳模拟过程，1表示裂纹扩展过程
%损伤矩阵及键的寿命矩阵初始化
numbond=0;
for i=1:size(fail,1)
    fail(i,1:nodefam(i,1))=1;
    lamda(i,1:nodefam(i,1))=1;
    flag(i,1:nodefam(i,1))=1;
    numbond=numbond+1;
end
basepoint=length(find(fail==1));
%下述分别为损伤矩阵，力，能量，位移，速度，加速度的初始化
dmg=zeros(dotdownn,1);%损伤
pforce=zeros(dotdownn,6);%PD力，包含fn,ft1,ft2,t,m1,m2
bforce=pforce;%外力
disp=pforce;%位移
fdisp=disp;
vel=pforce;%速度
acc=vel;%加速度
%% 选取若干键，使其一开始就断掉，形成初始缺陷
mpoint=0;
mdot=zeros(50,1);
for i = 1:allnum
%     if dotall(i,3)-0.5*(min(dotall(:,3))+max(dotall(:,3)))<15*dl
          mpoint=mpoint+1;
          mdot(mpoint)=i;
%     end
end
kb=ones(allnum,1);
anum=1:1:mpoint;
P=randperm(size(anum,2));
Pnum=mdot(P(1:800));
for i=1:800
    pm=Pnum(i);
    sbond=round(0.5*nodefam(pm,1));
    for j = 1:sbond%讨论每一点的进场范围内其他点的作用
        fail(pm,j)=0;
        fail(j,pm)=0;
    end
end
%% 设置不易坏区
%此算例中不设置不易坏区
nof=ones(allnum,1);
% nof=(1.0d-3)*ones(allnum,1);
%  for i=1:allnum
%      if abs(dotall(i,1)-0.5*(min(dotall(:,1))+max(dotall(:,1))))<=8*dl
%         if abs(dotall(i,2)-0.5*(min(dotall(:,2))+max(dotall(:,2))))<=8*dl
%            if abs(dotall(i,3)-0.5*(min(dotall(:,3))+max(dotall(:,3))))<=8*dl
%                nof(i)=1;
%            end
%         end
%      end
%  end
%% ADR算法变量初始化
velhalf=vel;%中间变量
velhalfold=vel;
pforceold=pforce;
massvec=zeros(dotdownn,6); %虚质量
for i = 1:dotdownn%虚拟质量   
    massvec(i,1) = 0.25d0 * dt * dt * (pi * (delta)^3) * bc1 / dl * 5.0d0;
    massvec(i,2) = 0.25d0 * dt * dt * (pi * (delta)^3) * bc1 / dl * 5.0d0;
    massvec(i,3) = 0.25d0 * dt * dt * (pi * (delta)^3) * bc1 / dl * 5.0d0;
    massvec(i,4) = 1.0d5 * 0.25d0 * dt * dt * (pi * (delta)^3) * bc1 / dl * 5.0d0;
    massvec(i,5) = 1.0d5 * 0.25d0 * dt * dt * (pi * (delta)^3) * bc1 / dl * 5.0d0;    
    massvec(i,6) = 1.0d5 * 0.25d0 * dt * dt * (pi * (delta)^3) * bc1 / dl * 5.0d0;
end
%% 开始积分求解
appbf1 = 20.94d6*len*hei/(vol*1200); %面力转为体力
appbf2 = 12.32d6*len*hei/(vol*1200); %面力转为体力
testbond=0.0d0;
timebond=zeros(nt,1); %记录每一次循环后的完好键数
aaamax=0;
for tt =  1:nt    
    fprintf('tt=%d\n',tt);    
    nc(tt,1)=tt/tao;
% 施加荷载  
    if tt<=5000
        bforce(allnum+1:dotlefn,1)= appbf2/5000;
        bforce(dotlefn+1:dotrign,1)= -appbf2/5000;   
        bforce(dotrign+1:dotupn,3)= -appbf1/5000;
    else
        bforce(allnum+1:dotlefn,1)= appbf2;
        bforce(dotlefn+1:dotrign,1)= -appbf2;   
        bforce(dotrign+1:dotupn,3)= -appbf1;
    end
    disp(dotupn+1:dotdownn,3)=0;
    disp(dotupn+1:dotdownn,2)=0;
    disp(dotupn+1:dotdownn,1)=0;  

    for i = 1:dotdownn  %对内部的每一个点都做循环     
        dmgpar1 = 0.0d0;%损伤变量初始化
        dmgpar2 = 0.0d0;
        pforce(i,1) = 0.0d0;%近场力初始化
        pforce(i,2) = 0.0d0;
        pforce(i,3) = 0.0d0;
        pforce(i,4) = 0.0d0;
        pforce(i,5) = 0.0d0; 
        pforce(i,6) = 0.0d0;
        for j = 1:nodefam(i,1)%讨论每一点的进场范围内其他点的作用           
            cnode = nodefam(i,j+1);
            idist = sqrt((dotall(cnode,1) - dotall(i,1))^2 + (dotall(cnode,2) - dotall(i,2))^2+(dotall(cnode,3) - dotall(i,3))^2);
            nlength = sqrt((dotall(cnode,1) + disp(cnode,1) - dotall(i,1) - disp(i,1))^2 + (dotall(cnode,2) + disp(cnode,2) - dotall(i,2) - disp(i,2))^2 + (dotall(cnode,3) + disp(cnode,3) - dotall(i,3) - disp(i,3))^2);
            eps(i,j)=(1-r)*((nlength - idist) / idist);
            nxylength=sqrt((dotall(cnode,1) + disp(cnode,1) - dotall(i,1) - disp(i,1))^2 + (dotall(cnode,2) + disp(cnode,2) - dotall(i,2) - disp(i,2))^2);%绕z轴转theta
            nzxlength=sqrt((dotall(cnode,1) + disp(cnode,1) - dotall(i,1) - disp(i,1))^2 + (dotall(cnode,3) + disp(cnode,3) - dotall(i,3) - disp(i,3))^2);%绕y轴转beta
            nyzlength=sqrt((dotall(cnode,2) + disp(cnode,2) - dotall(i,2) - disp(i,2))^2 + (dotall(cnode,3) + disp(cnode,3) - dotall(i,3) - disp(i,3))^2);%绕x轴转alpha
            %体积修正
            if (idist<=delta-radij)
                fac = 1.0d0;
            elseif (idist<=delta+radij)
                fac = (delta+radij-idist)/(2.0d0*radij);
            else
                fac = 0.0d0;
            end
            xl=(dotall(cnode,1) + disp(cnode,1) - dotall(i,1) - disp(i,1))/nlength;
            xm=(dotall(cnode,2) + disp(cnode,2) - dotall(i,2) - disp(i,2))/nlength;
            xn=(dotall(cnode,3) + disp(cnode,3) - dotall(i,3) - disp(i,3))/nlength;
            h=sqrt(xl*xl+xm*xm);
            if h~=0
               yl=xl*xn/h;
               ym=xm*xn/h;
               yn=-h;
            else
               yl=0;
               ym=0;
               yn=0;
            end
            zl=xm*yn-xn*ym;
            zm=xn*yl-xl*yn;
            zn=0;
            la11=xl;
            la12=xm;
            la13=xn;
            la21=yl;
            la22=ym;
            la23=yn;
            la31=zl;
            la32=zm;
            la33=zn;
            %位移分量
            v1=(disp(cnode,1)-disp(i,1));
            v2=(disp(cnode,2)-disp(i,2));
            v3=(disp(cnode,3)-disp(i,3));
            v4=disp(i,4);
            v5=disp(i,5);
            v6=disp(i,6);
            v7=disp(cnode,4);
            v8=disp(cnode,5);
            v9=disp(cnode,6);
            %局部刚度矩阵各分量
            k11=kn/idist;
            k22=kt/(idist.^3);
            k25=-ld*kt/(idist.^3);
            k28=k25;
            k33=k22;
            k36=k25;
            k39=k25;
            k44=kt/(12*(1+pratio)*idist);
            k47=-k44;
            k52=k25;
            k55=((ld^2)*kt/(idist^3))-kr;
            k58=((ld^2)*kt/(idist^3))+kr;
            k63=k25;
            k66=k55;
            k69=k58;
            %力的分量可以写为
            F1=v1*(k11*(la11^2)+k22*(la21^2)+k33*(la31^2))+...
               v2*(k11*la11*la12+k22*la21*la22+k33*la31*la32)+...
               v3*(k11*la11*la13+k22*la21*la23+k33*la31*la33)+...
               v4*(k25*la21^2+k36*la31^2)+v5*(k25*la21*la22+k36*la31*la32)+...
               v6*(k25*la21*la23+k36*la31*la33)+v7*(k28*la21^2+k39*la31^2)+...
               v8*(k28*la21*la22+k39*la31*la32)+v9*(k28*la21*la23+k39*la31*la33);
            F2=v1*(k11*la11*la12+k22*la21*la22+k33*la31*la32)+...
               v2*(k11*(la12^2)+k22*(la22^2)+k33*(la32^2))+...
               v3*(k11*la12*la13+k22*la22*la23+k33*la32*la33)+...
               v4*(k25*la21*la22+k36*la31*la32)+v5*(k25*la22^2+k36*la32^2)+...
               v6*(k25*la22*la23+k36*la32*la33)+v7*(k28*la21*la22+k39*la31*la32)+...
               v8*(k28*la22^2+k39*la32^2)+v9*(k28*la22*la23+k39*la32*la33);
            F3=v1*(k11*la11*la13+k22*la21*la23+k33*la31*la33)+...
               v2*(k11*la12*la13+k22*la22*la23+k33*la32*la33)+...
               v3*(k11*(la13^2)+k22*(la23^2)+k33*(la33^2))+...
               v4*(k25*la21*la23+k36*la31*la33)+v5*(k25*la22*la23+k36*la32*la33)+...
               v6*(k25*la32^2+k36*la33^2)+v7*(k28*la21*la23+k39*la31*la33)+...
               v8*(k28*la22*la23+k39*la32*la33)+v9*(k28*la23^2+k39*la33^2);
            F4=v1*(k52*la21^2+k63*la31^2)+v2*(k52*la21*la22+k63*la31*la32)+...
               v3*(k52*la21*la23+k63*la31*la33)+v4*(k44*la11^2+k55*la21^2+k66*la31^2)+...
               v5*(k44*la11*la12+k55*la21*la22+k66*la31*la32)+...
               v6*(k44*la11*la13+k55*la21*la23+k66*la31*la33)+...
               v7*(k47*la11^2+k58*la21^2+k69*la31^2)+...
               v8*(k47*la11*la12+k58*la21*la22+k69*la31*la32)+...
               v9*(k47*la11*la13+k58*la21*la23+k69*la31*la33);
            F5=v1*(k52*la21*la22+k63*la31*la32)+v2*(k52*la22^2+k63*la32^2)+...
               v3*(k52*la22*la23+k63*la32*la33)+v4*(k44*la11*la12+k55*la21*la22+k66*la31*la32)+...
               v5*(k44*la12^2+k55*la22^2+k66*la32^2)+v6*(k44*la12*la13+k55*la22*la23+k66*la32*la33)+...
               v7*(k47*la11*la12+k58*la21*la22+k69*la31*la32)+...
               v8*(k47*la12^2+k58*la22^2+k69*la32^2)+...
               v9*(k47*la12*la13+k58*la22*la23+k69*la32*la33);
            F6=v1*(k52*la21*la23+k63*la31*la33)+v2*(k52*la22*la23+k63*la32*la33)+...
               v3*(k52*la23^2+k63*la33^2)+v4*(k44*la11*la13+k55*la21*la23+k66*la31*la33)+...
               v5*(k44*la12*la13+k55*la22*la23+k66*la32*la33)+v6*(k44*la13^2+k55*la23^2+k66*la33^2)+...
               v7*(k47*la11*la13+k58*la21*la23+k69*la31*la33)+...
               v8*(k47*la12*la13+k58*la22*la23+k69*la32*la33)+...
               v9*(k47*la13^2+k58*la23^2+k69*la33^2);
            
            if (fail(i,j)==1)%如果这两点间的键没有坏，则有近场力作用
                if i>allnum||cnode>allnum
                      dforce1 = bc1 * eps(i,j) * vol * fac*(dotall(cnode,1)+disp(cnode,1)-dotall(i,1)-disp(i,1))/nlength;  
                      dforce2 = bc1 * eps(i,j) * vol * fac*(dotall(cnode,2)+disp(cnode,2)-dotall(i,2)-disp(i,2))/nlength;
                      dforce3 = bc1 * eps(i,j) * vol * fac*(dotall(cnode,3)+disp(cnode,3)-dotall(i,3)-disp(i,3))/nlength;
                   switcht(i,j)=1;                   
                end
                if i<=allnum&&cnode<=allnum
                  if eps(i,j)<=-scr3
                       dforce1 = 0.0d0;
                       dforce2 = 0.0d0;
                       dforce3 = 0.0d0;
                       dforce4 = 0.0d0;
                       dforce5 = 0.0d0;
                       dforce6 = 0.0d0;
                  elseif eps(i,j)<=-scr2&&eps(i,j)>=-scr3
                       k3=(s0/abs(eps(i,j)))*((scr3-abs(eps(i,j)))/(scr3-s0));          
                       dforce1 = k3 * F1 * vol * fac;
                       dforce2 = k3 * F2 * vol * fac;
                       dforce3 = k3 * F3 * vol * fac;
                       dforce4 = k3 * F4 * vol * fac;
                       dforce5 = k3 * F5 * vol * fac;
                       dforce6 = k3 * F6 * vol * fac;
                  elseif eps(i,j)<=-scr1&&eps(i,j)>-scr2
                       k2=(scr1*(abs(eps(i,j))-scr0))/((scr1-scr0)*abs(eps(i,j)));
                       dforce1 = k2 * F1 * vol * fac;
                       dforce2 = k2 * F2 * vol * fac;
                       dforce3 = k2 * F3 * vol * fac;
                       dforce4 = k2 * F4 * vol * fac;
                       dforce5 = k2 * F5 * vol * fac;
                       dforce6 = k2 * F6 * vol * fac;
                  elseif eps(i,j)<=scr1&&eps(i,j)>-scr1
                       dforce1 = F1 * vol * fac;
                       dforce2 = F2 * vol * fac;
                       dforce3 = F3 * vol * fac;
                       dforce4 = F4 * vol * fac;
                       dforce5 = F5 * vol * fac;
                       dforce6 = F6 * vol * fac;
                  elseif eps(i,j)>scr1&&eps(i,j)<=scr2
                       k2=(scr1*(abs(eps(i,j))-scr0))/((scr1-scr0)*abs(eps(i,j)));
                       dforce1 = k2 * F1 * vol * fac;
                       dforce2 = k2 * F2 * vol * fac;
                       dforce3 = k2 * F3 * vol * fac;
                       dforce4 = k2 * F4 * vol * fac;
                       dforce5 = k2 * F5 * vol * fac;
                       dforce6 = k2 * F6 * vol * fac;
                  elseif eps(i,j)>scr2&&eps(i,j)<=scr3
                       k3=(s0/abs(eps(i,j)))*((scr3-abs(eps(i,j)))/(scr3-s0));                        
                       dforce1 = k3 * F1 * vol * fac;
                       dforce2 = k3 * F2 * vol * fac;
                       dforce3 = k3 * F3 * vol * fac;
                       dforce4 = k3 * F4 * vol * fac;
                       dforce5 = k3 * F5 * vol * fac;     
                       dforce6 = k3 * F6 * vol * fac; 
                  elseif eps(i,j)>scr3
                       dforce1 = 0.0d0;
                       dforce2 = 0.0d0;
                       dforce3 = 0.0d0;
                       dforce4 = 0.0d0;
                       dforce5 = 0.0d0; 
                       dforce6 = 0.0d0;
                  end  
                end
            else%否则没有
                dforce1 = 0.0d0;
                dforce2 = 0.0d0;
                dforce3 = 0.0d0;
                dforce4 = 0.0d0;
                dforce5 = 0.0d0;
                dforce6 = 0.0d0;
            end
            dc1=dforce1*la11+dforce2*la21+dforce3*la31;
            dc2=dforce1*la12+dforce2*la22+dforce3*la32;
            dc3=dforce1*la13+dforce2*la23+dforce3*la33;
            dc4=dforce4*la11+dforce5*la21+dforce6*la31;
            dc5=dforce4*la12+dforce5*la22+dforce6*la32;
            dc6=dforce4*la13+dforce5*la23+dforce6*la33;            
            if i<=allnum&&cnode<=allnum
               bondforce(i,j)=sqrt(dc1^2)*vol/Ab;
               shearforce(i,j)=sqrt(dc2^2+dc3^2+(lc*(v7-v4)*dc4/(idist^2))^2+(lc*(v8-v5)*dc5/(idist^2))^2+(lc*(v9-v6)*dc6/(idist^2))^2)*vol/Ab;
               if eps(i,j)>=0
                  tell=1; %受拉
               else
                  tell=2; %受压
               end
            end
            %求力矩阵，将每个点的作用力叠加起来
            pforce(i,1) = pforce(i,1) + dforce1;
            pforce(i,2) = pforce(i,2) + dforce2;
            pforce(i,3) = pforce(i,3) + dforce3;
            pforce(i,4) = pforce(i,4) + dforce4;
            pforce(i,5) = pforce(i,5) + dforce5;
            pforce(i,6) = pforce(i,6) + dforce6;
            dforce1 = 0.0d0;
            dforce2 = 0.0d0;
            dforce3 = 0.0d0;
            dforce4 = 0.0d0;
            dforce5 = 0.0d0;
            dforce6 = 0.0d0;
           %进行键的伸长判断，键伸长达到一定值发生破坏,不去判断那些已经损伤的键     
%            if (i<=allnum)&&(cnode<=allnum)&&(abs(nlength - idist) / idist>scr3)&&(fail(i,j)==1)              
%                  fail(i,j)=0;%0代表键已经坏了
%                  flag(i,j)=0;
%                  switcht(i,j)=1.0d0;   %只要有键被拉坏，则跳过所有疲劳过程
%            end
           %进行键的寿命判断，键的寿命在疲劳过程中进行迭代
        if (i<=allnum)&&(cnode<=allnum)&&switcht(i,j)<=0.5   %判定是否还进行疲劳迭代过程
           if (tt>=2)&&(fail(i,j)==1)&&abs((nlength-idist)/idist)>einf&&nof(i)*nof(j)==1 
              if (dmg(i,1)<=0.38)&&(dmg(j,1)<=0.38)
                 aaa=sqrt(((bondforce(i,j)/ftc(tell)).^2+(shearforce(i,j)/fs).^2));
                 lamda(i,j)=lamda(i,j)-(nc(tt,1)-nc(tt-1,1))*a1*(aaa^m1); 
                 if aaamax<aaa
                     aaamax=aaa;
                 end 
              elseif ((dmg(i,1)>0.38)||(dmg(j,1)>0.38))&&(flag(i,j)<1.5)
                  %如果键的任意一端的点损伤值超过0.398，且这个键还未被判定为2阶段，则执行操作
                     lamda(i,j)=1.0d0;
                     flag(i,j)=2.0d0;   %这些键进入2阶段      
              end
              if flag(i,j) > 1.5 %二阶段的键进行二阶段迭代
                 aaa=sqrt(((bondforce(i,j)/ftc(tell)).^2+(shearforce(i,j)/fs).^2));
                 lamda(i,j)=lamda(i,j)-(nc(tt,1)-nc(tt-1,1))*a2*(aaa^m2);
              end               
           end
        end
        if (switcht(i,j)<=0.5)&&(i<=allnum)&&(cnode<=allnum)
            if (lamda(i,j)<=0)&&(fail(i,j)==1)
               fail(i,j)=0;
               flag(i,j)=0;
               lamda(i,j)=0;
            end
         end
         dmgpar1 = dmgpar1 + fail(i,j) *fac* vol ;%综合起来即为一点处的损伤值，通过损伤值的云图分布来反映具体的裂纹现象
         dmgpar2 = dmgpar2 + vol*fac  ;     
       end
        %计算损伤因子
        dmg(i,1) = 1.0d0 - dmgpar1 / dmgpar2;%损伤因子越大，损伤越大
    end
   %% ADR计算部分
    cn1 = 0.0d0;
    cn2 = 0.0d0;
    for i = 1:dotdownn
        if (velhalfold(i,1)~=0.0d0)
            cn1 = cn1 - disp(i,1) * disp(i,1) * (pforce(i,1) / massvec(i,1) - pforceold(i,1) / massvec(i,1)) / (dt * velhalfold(i,1));
        end
        if (velhalfold(i,2)~=0.0d0)
            cn1 = cn1 - disp(i,2) * disp(i,2) * (pforce(i,2) / massvec(i,2) - pforceold(i,2) / massvec(i,2)) / (dt * velhalfold(i,2));
        end
        if (velhalfold(i,3)~=0.0d0)
            cn1 = cn1 - disp(i,3) * disp(i,3) * (pforce(i,3) / massvec(i,3) - pforceold(i,3) / massvec(i,3)) / (dt * velhalfold(i,3));
        end
        if (velhalfold(i,4)~=0.0d0)
            cn1 = cn1 - disp(i,4) * disp(i,4) * (pforce(i,4) / massvec(i,4) - pforceold(i,4) / massvec(i,4)) / (dt * velhalfold(i,4));
        end
        if (velhalfold(i,5)~=0.0d0)
            cn1 = cn1 - disp(i,5) * disp(i,5) * (pforce(i,5) / massvec(i,5) - pforceold(i,5) / massvec(i,5)) / (dt * velhalfold(i,5));
        end
        if (velhalfold(i,6)~=0.0d0)
            cn1 = cn1 - disp(i,6) * disp(i,6) * (pforce(i,6) / massvec(i,6) - pforceold(i,6) / massvec(i,6)) / (dt * velhalfold(i,6));
        end        
        cn2 = cn2 + disp(i,1) * disp(i,1);
        cn2 = cn2 + disp(i,2) * disp(i,2);
        cn2 = cn2 + disp(i,3) * disp(i,3);
        cn2 = cn2 + disp(i,4) * disp(i,4);
        cn2 = cn2 + disp(i,5) * disp(i,5);
        cn2 = cn2 + disp(i,6) * disp(i,6);
    end
    
    if (cn2~=0.0d0)
        if ((cn1 / cn2) > 0.0d0)
            cn = 2.0d0 * sqrt(cn1 / cn2);
        else
            cn = 0.0d0;
        end
    else
        cn = 0.0d0;
    end
    
    if (cn > 2.0d0)
        cn = 1.9d0;
    end
    
    for i = 1:dotdownn
        if (tt==1)
            velhalf(i,1) = 1.0d0 * dt / massvec(i,1) * (pforce(i,1) + bforce(i,1)) / 2.0d0;
            velhalf(i,2) = 1.0d0 * dt / massvec(i,2) * (pforce(i,2) + bforce(i,2)) / 2.0d0;  
            velhalf(i,3) = 1.0d0 * dt / massvec(i,3) * (pforce(i,3) + bforce(i,3)) / 2.0d0;
            velhalf(i,4) = 1.0d0 * dt / massvec(i,4) * (pforce(i,4) + bforce(i,4)) / 2.0d0;
            velhalf(i,5) = 1.0d0 * dt / massvec(i,5) * (pforce(i,5) + bforce(i,5)) / 2.0d0;
            velhalf(i,6) = 1.0d0 * dt / massvec(i,6) * (pforce(i,6) + bforce(i,6)) / 2.0d0;
        else
            velhalf(i,1) = ((2.0d0 - cn * dt) * velhalfold(i,1) + 2.0d0 * dt / massvec(i,1) * (pforce(i,1) + bforce(i,1))) / (2.0d0 + cn * dt);
            velhalf(i,2) = ((2.0d0 - cn * dt) * velhalfold(i,2) + 2.0d0 * dt / massvec(i,2) * (pforce(i,2) + bforce(i,2))) / (2.0d0 + cn * dt);  
            velhalf(i,3) = ((2.0d0 - cn * dt) * velhalfold(i,3) + 2.0d0 * dt / massvec(i,3) * (pforce(i,3) + bforce(i,3))) / (2.0d0 + cn * dt);
            velhalf(i,4) = ((2.0d0 - cn * dt) * velhalfold(i,4) + 2.0d0 * dt / massvec(i,4) * (pforce(i,4) + bforce(i,4))) / (2.0d0 + cn * dt);
            velhalf(i,5) = ((2.0d0 - cn * dt) * velhalfold(i,5) + 2.0d0 * dt / massvec(i,5) * (pforce(i,5) + bforce(i,5))) / (2.0d0 + cn * dt);
            velhalf(i,6) = ((2.0d0 - cn * dt) * velhalfold(i,6) + 2.0d0 * dt / massvec(i,6) * (pforce(i,6) + bforce(i,6))) / (2.0d0 + cn * dt);
        end        
        vel(i,1) = 0.5d0 * (velhalfold(i,1) + velhalf(i,1))	;
        vel(i,2) = 0.5d0 * (velhalfold(i,2) + velhalf(i,2))	; 
        vel(i,3) = 0.5d0 * (velhalfold(i,3) + velhalf(i,3))	; 
        vel(i,4) = 0.5d0 * (velhalfold(i,4) + velhalf(i,4))	;
        vel(i,5) = 0.5d0 * (velhalfold(i,5) + velhalf(i,5))	;
        vel(i,6) = 0.5d0 * (velhalfold(i,6) + velhalf(i,6))	;
        disp(i,1) = disp(i,1) + velhalf(i,1) * dt	;
        disp(i,2) = disp(i,2) + velhalf(i,2) * dt	; 
        disp(i,3) = disp(i,3) + velhalf(i,3) * dt	;
        disp(i,4) = disp(i,4) + velhalf(i,4) * dt	;
        disp(i,5) = disp(i,5) + velhalf(i,5) * dt	;
        disp(i,6) = disp(i,6) + velhalf(i,6) * dt	;
        velhalfold(i,1) = velhalf(i,1);
        velhalfold(i,2) = velhalf(i,2); 
        velhalfold(i,3) = velhalf(i,3);
        velhalfold(i,4) = velhalf(i,4);
        velhalfold(i,5) = velhalf(i,5);
        velhalfold(i,6) = velhalf(i,6);
        pforceold(i,1) = pforce(i,1);
        pforceold(i,2) = pforce(i,2);  
        pforceold(i,3) = pforce(i,3); 
        pforceold(i,4) = pforce(i,4); 
        pforceold(i,5) = pforce(i,5); 
        pforceold(i,6) = pforce(i,6); 
    end
    % 基本数据输出，损伤分布以及位移分布
    if(mod(tt,10)==0)%只要时间步长是10的倍数，就进行输出，输出的文件里面有各点的位移和损伤值，若想输出其他的可自己定义。         
         str1=[directpath,'dispdmg',num2str(tt),'.mat'];
         dispdmg=[disp(dotn,:),dmg(dotn,1)];
         save(str1,'dispdmg') 
%          subplot(2,1,1)
%          scatter_color(dotall(1:allnum,:),dispdmg(1:allnum,1))
%          subplot(2,1,2)
%          scatter_color(dotall(1:allnum,:),dispdmg(1:allnum,2))
    end
   timebond(tt,1)=numbond;
end
% iii=zeros(1,800);
% ijk=0;
% for i=1:allnum
%     if dotall(i,2)==0
%         ijk=ijk+1;
%         iii(ijk)=i;
%     end
% end
% figure (1)
load('G:\fatigue_work001\dispdmg4300')
scatter_color3d(dotall(1:allnum,:),dispdmg(1:allnum,7))
scatter_color3d(dotall(1:allnum,:),pforce(1:allnum,3))
% scatter_color(dotall(iii(:),:),dispdmg(iii(:),7))