%混凝土三点弯曲梁算例
%CSPD
%tao=1.0
%设置约束处和不坏区
clc;clear
%% 创建计算结果保存路径
directpath='F:\fatigue_01_01\';
if (~exist(directpath,'dir') )
    mkdir(directpath)        % 若不存在，在当前目录中产生一个子目录
end
%% 参数输入
len=406.4d-3; %原模型横向长度
wid=152.4d-3; %原模型竖向长度
thick=38.1d-3; %原模型高度
xdv=200;  %x向划分
ydv=75;  %y向划分
dl=len/xdv;  %节点之间水平或者竖直方向间距
allnum=xdv*ydv;  %平面节点总数 p表示plane
dot=zeros(allnum,2);  %将初始节点坐标（默认为零）存入dotp矩阵
n=0; %节点编号初始化
% 材料参数
tao=4.0d0;
nt=3000+1000*tao;
dt=1.0d0;
e1 = 27.12d9;   %弹性模量
pratio = 0.2d0;   %泊松比
delta=3.015*dl;    %近场范围尺寸
bc1 = 6.0d0 * e1 / (pi *  (1-pratio)*(delta^4));
admg=zeros(nt,1);
nc=zeros(nt,1);
r=0.1d0;
einf=0;

fc=28.9624d6;  %混凝土抗压强度
ft=fc/10.1267;
%迭代参数
% a1=0.6812;
% m1=37.037;
a1=0.1492;
m1=21.7391;
% a2=2.7;
a2=1.0;
m2=9.27;
%Cosserat参数
G=e1/(2*(1+pratio));%剪切模量
Gc = 0.5d0.*G;
lc= 1.0d-4;
kn=6*e1/(thick*pi*(delta.^3)*(1-pratio));
kt=2*e1*(1-3*pratio)/(thick*pi*delta*(1-pratio^2));
ld=sqrt(Gc*(1-pratio^2)/(2*e1*(1-3*pratio)))*delta;
kr=4*Gc/(thick*pi*delta^2)-16*G*(lc^2)/(thick*pi*delta^4);

%输入键的拉伸和压缩时的参数
GF=0.0422*1.6289; %能量释放率
scr1t=sqrt(((1-2*pratio)*pi*GF)/(e1*delta));%临界的变形,超过这个值就会破坏
fmaxt=scr1t*kn;%%%
fmaxn=fmaxt*fc/ft;
%输入混凝土强度包络线的三个系数
mf=0.011094;
nf=0.09974;
pf=0.10907;
%% 建模
%主体
for i=1:xdv
    for j=1:ydv
        n=n+1;
        dot(n,1)=(-1.0d0*len./2.0d0)+(i-1).*dl;
        dot(n,2)=(-1.0d0*wid./2.0d0)+(j-1).*dl;
        %从左下角往右上角
    end
end
mu = zeros(1,n);
%设置预制裂纹
for i = 1:numel(mu)
    mu(1,i)=0;
%     排除裂纹处的物质点
    if abs(dot(i,1)-0.5*(max(dot(:,1))+min(dot(:,1))))>=0.6*dl
       mu(1,i)=i;
    else 
       if dot(i,2)-min(dot(1:allnum,2))>=wid/6
          mu(1,i)=i;
       end
    end
end
mu(find(mu==0))=[];
dot2 = zeros(numel(mu),2);
for i=1:numel(mu)
    dot2(i,1)=dot(mu(1,i),1);
    dot2(i,2)=dot(mu(1,i),2);
end
allnum = numel(mu);
% allnum = n;
dot0=zeros(xdv*3,2);
n=0;
% 加入垫片
for i=1:xdv
    for j=1:3
        n=n+1;
        dot0(n,1)=(-1.0d0*len./2.0d0)+(i-1).*dl;
        dot0(n,2)=max(dot(:,2))+dl+(j-1).*dl;
        %从左下角往右上角
    end
end
allnum0=allnum+n;
n=0;
%加入加载块
dlens=5.0d0*dl;
dwids=5.0d0*dl;
xbdvs=round(dlens/dl);
ybdvs=round(dwids/dl);
alldots=xbdvs*ybdvs; %每个加载行的节点数
dotload=zeros(11*3,2);
%将节点横纵坐标按编号存入dotload矩阵
m=0;
for i=1:11
    for j=1:3
        m=m+1;
        dotload(m,1)=0.5*(min(dot(:,1))+max(dot(:,1)))-5*dl+(i-1).*dl;
        dotload(m,2)=max(dot0(:,2))+dl+(j-1).*dl;
    end
end
dotloadn=allnum0+m;
%加入约束块,分左下及右下块
m=0;
dotlef=zeros(alldots,2);  %左块体节点坐标存入此矩阵
dotrit=zeros(alldots,2);  %右块体节点坐标存入此矩阵
%下面将节点横纵坐标值按编号对应存入dotlef和dotrig矩阵
for i=1:xbdvs
    for j=1:ybdvs
        m=m+1;
        dotlef(m,1)=min(dot(:,1))+(i-1).*dl;
        dotlef(m,2)=min(dot(:,2))-5*dl+(j-1).*dl;
    end
end
dotlefn=dotloadn+m;
m=0;
for i=1:xbdvs
    for j=1:ybdvs
        m=m+1;
        dotrit(m,1)=max(dot(:,1))-4*dl+(i-1).*dl;
        dotrit(m,2)=min(dot(:,2))-5*dl+(j-1).*dl;
    end
end
dotritn=dotlefn+m;
dotdownn=dotritn;
dotall=[dot2;dot0;dotload;dotlef;dotrit];  %将所有所要的节点存入一个矩阵中
dotn=1:1:allnum; %内点编号
%% 观察整体构型 不用时可以注释掉
% k1=ones(dotritn,1);
% scatter_color(dotall,k1)
% scatter_color(dotall(1:allnum),nof)
%% 其他参数初始化
area = dl * dl;%面积
vol = area * thick;%体积
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;

nodefam=knn(dotall,delta,dl);%总的近场范围点的位置矩阵构造,用的是knn函数
%第一列是统计一共有多少个近场点的，所以循环都应该从第二个点开始
nhorizon=max(nodefam(:,1));%读取nodefam中的列数
fncst=ones(dotdownn,nhorizon);%表面修正系数
fail=zeros(dotdownn,nhorizon);%损伤  1代表没破坏，0代表破坏了
lamda=zeros(dotdownn,nhorizon);%键的剩余寿命 1代表完好，0代表破坏
flag=zeros(dotdownn,nhorizon); %判断每一根键的阶段 1代表1阶段，2代表2阶段
eps=zeros(dotdownn,nhorizon); %每根键的应变
sn=zeros(dotdownn,nhorizon); %每根键的法向作用
st=zeros(dotdownn,nhorizon);%每根键的切向作用
aaa=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,3);%PD力
bforce=pforce;%外力
disp=pforce;%位移
vel=pforce;%速度
acc=vel;%加速度
aaamax=0;
aaamaxi=0;
aaamaxj=0;
% %% 预置裂纹
% %写出裂纹中心点的x，y坐标
% crackx=0.5*(max(dot(:,1))+min(dot(:,1)));
% crlength=wid/6;
% cracky=min(dot(:,2))+wid/12-0.5*dl;
% for i = 1:allnum  
%     cxx=dotall(i,1)-crackx;
%     cyy=dotall(i,2)-cracky;  
%     for j = 1:nodefam(i,1)
%         cnode =nodefam(i,j+1);
%         cnxx=dotall(cnode,1)-crackx;
%         cnyy=dotall(cnode,2)-cracky;
%         if (cxx*cnxx<0)
%             if(abs(cnyy)<=crlength / 2.0d0)
%                 fail(i,j)=0;
%             elseif(abs(cyy)<=crlength / 2.0d0)
%                 fail(i,j)=0;
%             end
%         end
%     end
% end
%% ADR算法变量初始化
velhalf=vel;%中间变量
velhalfold=vel;
pforceold=pforce;
massvec=zeros(dotdownn,3); %虚质量
for i = 1:dotdownn%虚拟质量   
    massvec(i,1) = 0.25d0 * dt * dt * (pi * (delta)^2 * thick) * bc1 / (dl) ;
    massvec(i,2) = 0.25d0 * dt * dt * (pi * (delta)^2 * thick) * bc1 / (dl) ; 
    massvec(i,3) = (1.0d3)*0.25d0 * dt * dt * (pi * (delta)^2 * thick) * bc1 / (dl) ;
end
flag2point=0;
%% 设置不坏区
% 此算例中不设置不坏区
nof=ones(allnum,1);
 for i=1:allnum
   if (dotall(i,2)>=max(dotall(1:allnum,2))-24*dl)
      nof(i)=0; 
   end
   if (dotall(i,2)<min(dotall(1:allnum,2))+12*dl)&&abs(dotall(i,1)-0.5*(max(dotall(1:allnum,1))+min(dotall(1:allnum,1))))>=6*dl
      nof(i)=0;
   end   
   if abs(dotall(i,1)-0.5*(max(dotall(1:allnum,1))+min(dotall(1:allnum,1))))>60*dl
      nof(i)=0;
   end        
 end
%  scatter_color(dotall(1:allnum,:),nof)
%% 开始积分求解
Abond=0;
alldmg=zeros(nt,1);
appbf1 = 5184.2*0.8/(vol*(dotloadn-allnum0));
% displacement=0.004*25.4d-3;
sbond=0.0d0; %只要有键的临界伸长率超过scr0，sbond就不为0
testbond=0.0d0;
timebond=zeros(nt,1); %记录每一次循环后的完好键数
kao=0;
tell=0;
crack=zeros(nt,1);
beta = 0.3;
%加入底座支反力-时间曲线。
lefforce=zeros(nt,1);
ritforce=zeros(nt,1);
for tt =  1:nt    
    fprintf('tt=%d\n',tt);    
    nc(tt,1)=tt/tao;
% 施加荷载与约束
    if tt<=3000
%        disp(allnum0+1:dotloadn,2)=-displacement*tt/10000;
       bforce(allnum0+1:dotloadn,2)=-appbf1*tt/3000;  
    else
%        disp(allnum0+1:dotloadn,2)=-displacement;
       bforce(allnum0+1:dotloadn,2)=-appbf1; 
    end
       disp(dotloadn+1:dotritn,2)=0;

    if sbond > 0.0d0
        time=tt;    %记录哪一步时进入静载破坏阶段
    end
    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;
        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);
            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);
            eps(i,j)=(1-r)*((nlength - idist) / idist);
            Abond=pi*delta^2*thick/(464*idist);
            %体积修正
            if (idist<=delta-radij)
                fac = 1.0d0;
            elseif (idist<=delta+radij)
                fac = (delta+radij-idist)/(2.0d0*radij);
            else
                fac = 0.0d0;
            end
            if (abs(dotall(cnode,2) - dotall(i,2))<=1.0d-10)
                theta = 0.0d0;
            elseif (abs(dotall(cnode,1) - dotall(i,1))<=1.0d-10)
                theta = 90.0d0 * pi / 180.0d0;
            else
                theta = atan(abs(dotall(cnode,2) - dotall(i,2)) / abs(dotall(cnode,1) - dotall(i,1)));
            end
            
            fncst(i,j)=2*nhorizon/(nodefam(i,1)+nodefam(cnode,1));
            scr=fncst(i,j);
            
            theat1=(dotall(cnode,1)+disp(cnode,1)-dotall(i,1)-disp(i,1))/nlength;%键的角度值，cos
            theat2=(dotall(cnode,2)+disp(cnode,2)-dotall(i,2)-disp(i,2))/nlength;%键的角度值，sin
              k11=((6*e1)./(thick*pi*idist*(1-pratio)*(delta.^3)));
              k22=(((2*e1*(1-3*pratio))./(thick*pi*delta*(idist.^3)*(1-pratio.^2))));
              k23=(-(sqrt((2*Gc*e1*(1-3*pratio))./(1-pratio.^2))./(thick*pi*(idist.^3))));
              k24=(-(sqrt((2*Gc*e1*(1-3*pratio))./(1-pratio.^2))./(thick*pi*(idist.^3))));
              k32=(-(sqrt((2*Gc*e1*(1-3*pratio))./(1-pratio.^2))./(thick*pi*(idist.^3))));
              k33=(((Gc*delta)./(thick*pi*(idist.^3))-(4*Gc)./(thick*pi*(delta.^2))+(16*G*lc.^2)./(thick*pi*delta.^4)));
              k34=(((Gc*delta)./(thick*pi*(idist.^3))-(4*Gc)./(thick*pi*(delta.^2))-(16*G*lc.^2)./(thick*pi*delta.^4)));           
              v1=(disp(cnode,1)-disp(i,1));
              v2=(disp(cnode,2)-disp(i,2));
              v3=disp(i,3);
              v4=disp(cnode,3);
              F1=(k11*theat1.^2+k22*theat2.^2)*v1+(k11*theat1*theat2-k22*theat1*theat2)*v2-k23*theat2*v3-k24*theat2*v4;
              F2=(k11*theat1*theat2-k22*theat1*theat2)*v1+(k11*theat2.^2+k22*theat1.^2)*v2+k23*theat1*v3+k24*theat1*v4;
              F3=k32*(-theat2)*v1+k32*theat1*v2+k33*v3+k34*v4;   
              F1c=k11*v1*theat1+k11*v2*theat2;
              F2c=-k22*v1*theat2+k22*theat1*v2+k23*v3+k24*v4;
              F3c=-k32*v1*theat2+k32*theat1*v2+k33*v3+k34*v4;
            if (fail(i,j)==1)%如果这两点间的键没有坏，则有近场力作用
               dforce1 = F1 * vol * fac * scr;
               dforce2 = F2 * vol * fac * scr;
               dforce3 = F3 * vol * fac;
               dc1 = F1c * fac * scr;
               dc2 = F2c * fac * scr;
               dc3 = F3c * fac;                      
            else%否则没有
                dforce1 = 0.0d0;
                dforce2 = 0.0d0;
                dforce3 = 0.0d0;
                dc1 = 0.0d0;
                dc2 = 0.0d0;
                dc3 = 0.0d0;             
            end

            % 预计算迭代因子aaa
            if i<=allnum&&cnode<=allnum
               sn(i,j)=sqrt(dc1^2)/fmaxn;
               st(i,j)=beta*sqrt(dc2^2+(dc3/(idist))^2)/fmaxn; 
            end

            if sn(i,j)==0&&st(i,j)==0%表示处于坐标原点
               aaa(i,j)=0;
            elseif sn(i,j)==0&&st(i,j)~=0% sigma n =0,按纯剪切考虑
               aaa(i,j)=st(i,j)/sqrt(mf);
            elseif eps(i,j)>0&&st(i,j)==0% sigma t =0, 按纯拉考虑
               aaa(i,j)=-sn(i,j)/((nf/2/pf)-sqrt((nf/2/pf)^2+mf/pf));
            elseif eps(i,j)<=0&&st(i,j)==0% sigma t =0, 按纯压考虑
               aaa(i,j)=sn(i,j)/((nf/2/pf)+sqrt((nf/2/pf)^2+mf/pf));
            elseif eps(i,j)>0&&st(i,j)~=0&&sn(i,j)~=0%受拉为负
               k=-st(i,j)/sn(i,j);
               stmax=0.5*k*nf/(k^2+pf)+sqrt((0.5*k*nf/(k^2+pf)).^2+mf*k^2/(k^2+pf));
               snmax=stmax/k;
               aaa(i,j)=-sn(i,j)/snmax;
            elseif eps(i,j)<=0&&st(i,j)~=0&&sn(i,j)~=0%受压为正
               k=st(i,j)/sn(i,j);
               stmax=0.5*k*nf/(k^2+pf)+sqrt((0.5*k*nf/(k^2+pf)).^2+mf*k^2/(k^2+pf));
               snmax=stmax/k;
               aaa(i,j)=sn(i,j)/snmax;                            
            end

            if aaamax<aaa(i,j)
               aaamax=aaa(i,j);
               aaamaxi=i;
               aaamaxj=j;
            end 
            %求力矩阵，将每个点的作用力叠加起来
            pforce(i,1) = pforce(i,1) + dforce1;
            pforce(i,2) = pforce(i,2) + dforce2;
            pforce(i,3) = pforce(i,3) + dforce3;
                 
           %进行键的伸长判断，键伸长达到一定值发生破坏,不去判断那些已经损伤的键     
%            if tt>100&&(fail(i,j)==1)&&i<=allnum&&eps(i,j)>=scr1t&&nof(i)>0            
%                  fail(i,j)=0;%0代表键已经坏了
% %                  flag(i,j)=0;
%                  switcht(i,j)=1.0d0;   %只要有键被拉坏，则跳过所有疲劳过程
%            end
%            进行键的寿命判断，键的寿命在疲劳过程中进行迭代
        if tt>=3000&&(fail(i,j)==1)&&(i<=allnum)&&(cnode<=allnum)&&switcht(i,j)<=0.5   %判定是否还进行疲劳迭代过程
           if dmg(i)>0.398||dmg(cnode)>0.398  %有问题
              for k=1:nodefam(i,1)
                 if (flag(i,k)>0)&&(flag(i,k)<1.5)%1阶段的键进入循环
                     lamda(i,k)=1.0;
                     flag(i,k)=2.0;%这些键进入2阶段
                 end   
              end   
           end
        end
       end
    end
    for i = 1:dotdownn  %对内部的每一个点都做循环
        dmgpar1 = 0.0d0;%损伤变量初始化
        dmgpar2 = 0.0d0;
      for j = 1:nodefam(i,1)%讨论每一点的进场范围内其他点的作用
        if tt>=3000&&(flag(i,j)>0)&&(flag(i,j)<1.5)&&i<=allnum
            lamda(i,j)=lamda(i,j)-nof(i)*(nc(tt,1)-nc(tt-1,1))*a1*((aaa(i,j))^m1);           
        elseif flag(i,j)>1.5 %二阶段的键进行二阶段迭代
            lamda(i,j)=lamda(i,j)-nof(i)*(nc(tt,1)-nc(tt-1,1))*a2*((aaa(i,j))^m2);                
        end
              
        if (lamda(i,j)<=0)&&(fail(i,j)==1)
            fail(i,j)=0;
            flag(i,j)=0;
        end

         dmgpar1 = dmgpar1 + fail(i,j) *fac* vol ;%综合起来即为一点处的损伤值，通过损伤值的云图分布来反映具体的裂纹现象
         dmgpar2 = dmgpar2 + vol*fac  ;     
       end
        %计算损伤因子
        dmg(i,1) = 1.0d0 - dmgpar1 / dmgpar2;%损伤因子越大，损伤越大
    end    
%% 找裂纹长度的代码
    snumber=0;
    jud_cl=zeros(10000,1);
    for ii=1:allnum
        if dotall(ii,2) >= min(dotall(:,2))+wid/6
           if (dmg(ii,1)>0.0001)
              snumber=snumber+1;
              jud_cl(snumber,1)=ii;
           end
        end
    end
    jud_cl=jud_cl(jud_cl~=0);
     if (tt==1)
         crack(tt,1)=min(dotall(1:allnum,2));
     end
     if (tt~=1)
         if(max(jud_cl)~=0)
             crack(tt,1)=max(dotall(jud_cl,2))-crack(1,1);
         else
             crack(tt,1)=0.0d0;
         end
     end
    %% ADR计算部分
    cn1 = 0.0d0;
    cn2 = 0.0d0;
    
    for i = 1:allnum
        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      
        cn2 = cn2 + disp(i,1) * disp(i,1);
        cn2 = cn2 + disp(i,2) * disp(i,2);
        cn2 = cn2 + disp(i,3) * disp(i,3);
    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;
        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);
        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))	; 
        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	;
        velhalfold(i,1) = velhalf(i,1);
        velhalfold(i,2) = velhalf(i,2);   
        velhalfold(i,3) = velhalf(i,3);
        pforceold(i,1) = pforce(i,1);
        pforceold(i,2) = pforce(i,2);  
        pforceold(i,3) = pforce(i,3); 
    end
    % 基本数据输出，损伤分布以及位移分布
    if(mod(tt,10)==0)%只要时间步长是10的倍数，就进行输出，输出的文件里面有各点的位移和损伤值，若想输出其他的可自己定义。         
         str1=[directpath,'dispdmg',num2str(tt),'.mat'];
         dispdmg=[disp(dotn,:),dmg(dotn,1)];
         save(str1,'dispdmg')  
%          if tt>=3000
%              figure(10)
%              subplot(2,1,1)
%              scatter_color(dotall(1:allnum,:),dispdmg(:,2))  
%              subplot(2,1,2)
%              scatter_color(dotall(1:allnum,:),dispdmg(:,4))    
%          end
    end
   timebond(tt,1)=numbond;
   for i=dotloadn+1:dotlefn
       lefforce(tt)=lefforce(tt)+pforce(i,2)*vol;
   end
   for i=dotlefn+1:dotritn
       ritforce(tt)=ritforce(tt)+pforce(i,2)*vol;
   end
end
crack(1,1)=crack(2,1);
save([directpath,'crack_y_length01',],'crack');
save([directpath,'lefforce01',],'lefforce');
save([directpath,'ritforce01',],'ritforce');
% figure (100)
% load('F:\fatigue_01_02\dispdmg3500')         
% scatter_color(dotall(1:allnum,:),dispdmg(1:allnum,4))
% set(gca,'FontSize',16,'FontName','Times New Roman');
% scatter_color(dotall(1:allnum,:),log(abs(pforce(1:allnum,2))))
% scatter_color(dotall(1:allnum,:),log(abs(pforce(1:allnum,2)))) 
load('F:\fatigue_01_01\crack_y_length01')
t=1:1:4000;
% crack(1:2000,1)=wid/6;
for i=1:7000
    if crack(i,1)<=wid/5
       crack(i,1)=wid/5;
    end
end
plot(t/4,crack(3001:7000,1)/wid)
hold on
% xe=[100,200,300,400,500,600,700,800,900,1000];
% ye=[0.2026,0.2053,0.2078,0.2146,0.2200,0.2256,0.2378,0.2597,0.3001,0.3301];
% plot(xe,ye)
xlabel('Cycle')
ylabel(' Crack length(m)')
set(gca,'FontSize',16,'FontName','Times New Roman');
% plot(1:2000,-ritforce(1:2000))