%算例1：混凝土立方块
% CSPD (Peridynamics based on Cosserat continuum) 、
%引入三段式本构
%引入应力比疲劳循环
clc;clear
%% 创建计算结果保存路径
directpath='G:\fatigue_work004\';
if (~exist(directpath,'dir') )
    mkdir(directpath)        % 若不存在，在当前目录中产生一个子目录
end
%% 参数输入
len=100.0d-3; %原模型横向长度
wid=100.0d-3; %原模型竖向长度
xdv=100;  %x向划分
ydv=100;  %y向划分
dl=len/xdv;  %节点之间水平或者竖直方向间距
dt=1.0d0;%虚拟时间步长
nt=60000;%总时间步数
allnum=xdv*ydv;  %节点总数
dot=zeros(allnum,2);  %将初始节点坐标（默认为零）存入dot矩阵
n=0; %节点编号初始化
%材料参数
thick =100.0d-3;%厚度
dens = 2400.0;%密度
emod = 24.6d9;   %弹性模量
emod2 = 200.0d0*emod;
pratio = 0.2d0;   %泊松比
G=emod/(2*(1+pratio));%剪切模量
delta=3.015*dl;    %近场范围尺寸
bc = 6.0d0 * emod / (pi * thick * (1-pratio)*(delta^3));
bc2 = 6.0d0 * emod2 / (pi * thick * (1-pratio)*(delta^3)); 
admg=zeros(nt,1);
nc=zeros(nt,1);
tao=4.0d0;
r=0.2353d0;
% scr0=sqrt((4.0*3488)/(bc*thick*(delta.^4)));%临界的变形,超过这个值就会破坏
einf=0;
Gc = 0.5d0*G;
value=1.0d0;
Gc2=Gc*200;
G2=G*200;
lc= 5.0d-5;
crackxl=zeros(nt,1);
fc=36.5d6;  %混凝土抗压强度
%迭代参数
a1=0.0184;
m1=19.53;
a2=0.0184;
m2=10.5;
scr1=0.0023;
scr2=0.006;
scr3=0.018718;
scr0=-0.00886;
s0=0.0045;
%% 下面将节点横纵坐标值按编号对应存入dot矩阵
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
allnum=n;
% 加入加载块 分左右部和上部
m=0;
%左右加载行
dwids=100.0d-3;
dlens=3.0d0*dl;
xbdvs=round(dlens/dl);
ybdvs=round(dwids/dl);
alldots=xbdvs*ybdvs; %每个加载行的节点数
dotlef=zeros(alldots,2);  %左块体节点坐标存入此矩阵
dotrit=zeros(alldots,2);  %右块体节点坐标存入此矩阵
%下面将节点横纵坐标值按编号对应存入dotlef和dotrig矩阵
for i=1:xbdvs
    for j=1:ybdvs
        m=m+1;
        dotlef(m,1)=(-1.0d0*len./2.0d0)-dlens-dl+(i-1).*dl;
        dotlef(m,2)=(-1.0d0*wid./2.0d0)+(j-1).*dl;
    end
end
dotlefn=allnum+m;
m=0;
for i=1:xbdvs
    for j=1:ybdvs
        m=m+1;
        dotrit(m,1)=(1.0d0*len./2.0d0)+dl+(i-1).*dl;
        dotrit(m,2)=(-1.0d0*wid./2.0d0)+(j-1).*dl;
    end
end
dotritn=dotlefn+m;
m=0;
%上加载行
dlenb=100.0d-3;
dwidb=3.0d0*dl;
xbdvb=round(dlenb/dl);
ybdvb=round(dwidb/dl);
alldotb=xbdvb*ybdvb; %每个加载行的节点数
dotup=zeros(alldotb,2);  %上块体节点坐标存入此矩阵
dotdown=zeros(alldotb,2);  %下块体节点坐标存入此矩阵
%下面将节点横纵坐标值按编号对应存入dotup和dotdown矩阵
for i=1:xbdvb
    for j=1:ybdvb
        m=m+1;
        dotup(m,1)=(-1.0d0*len./2.0d0)+(i-1).*dl;
        dotup(m,2)=max(dot(:,2))+2.0d0*dl+(j-1).*dl;
    end
end
dotupn=dotritn+m;
m=0;
for i=1:xbdvb
    for j=1:ybdvb
        m=m+1;
        dotdown(m,1)=(-1.0d0*len./2.0d0)+(i-1).*dl;
        dotdown(m,2)=min(dot(:,2))-dwidb+(j-1).*dl;
    end
end
dotdownn=dotupn+m;
dotall=[dot(:,:);dotlef;dotrit;dotup;dotdown];  %将所有所要的节点存入一个矩阵中
dotn=1:1:allnum; %内点编号
%% 观察整体构型 不用时可以注释掉
% k1=ones(dotdownn,1);
% scatter_color(dotall,k1)
%% 其他参数初始化
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,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); %每根键的应变
%损伤矩阵及键的寿命矩阵初始化
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;%加速度
bondforce=zeros(dotdownn,nhorizon); %每根键的力
%% 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) * bc / (dl) ;
    massvec(i,2) = 0.25d0 * dt * dt * (pi * (delta)^2 * thick) * bc / (dl) ; 
    massvec(i,3) = (1.0d3)*0.25d0 * dt * dt * (pi * (delta)^2 * thick) * bc / (dl) ;
end
flag2point=0;
%% 开始积分求解
switcht=0; %判定过程开关,0表示疲劳模拟过程，1表示裂纹扩展过程
appbf1 = 20.94d6/thick; %面力转为体力
appbf2 = 12.32d6/thick; %面力转为体力
fk=fc;
sbond=0.0d0; %只要有键的临界伸长率超过scr0，sbond就不为0
testbond=0.0d0;
timebond=zeros(nt,1); %记录每一次循环后的完好键数
kao=0;
for tt =  1:nt    
    fprintf('tt=%d\n',tt);    
    nc(tt,1)=tt/tao;
% 施加荷载  
    bforce(allnum+1:dotlefn,1)= appbf2;
    bforce(dotlefn+1:dotritn,1)= -appbf2;   
    bforce(dotritn+1:dotupn,2)= -appbf1;
    disp(dotupn+1:dotdownn,2)=0;
    disp(dotupn+1:dotdownn,1)=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);
            %体积修正
            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
            scx = (fncst(i,1) + fncst(cnode,1)) / 2.0d0;
            scy = (fncst(i,2) + fncst(cnode,2)) / 2.0d0;
            
            scr = 1.0d0/(((cos(theta)).^2/(scx).^2)+((sin(theta)).^2/(scy).^2));
            scr = sqrt(scr);
            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*emod)./(thick*pi*idist*(1-pratio)*(delta.^3)));
              k22=(((2*emod*(1-3*pratio))./(thick*pi*delta*(idist.^3)*(1-pratio.^2))));
              k23=(-(sqrt((2*Gc*emod*(1-3*pratio))./(1-pratio.^2))./(thick*pi*(idist.^3))));
              k24=(-(sqrt((2*Gc*emod*(1-3*pratio))./(1-pratio.^2))./(thick*pi*(idist.^3))));
              k32=(-(sqrt((2*Gc*emod*(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)));
              %加载块的物质点按刚体处理
              k110=((6*emod2)./(thick*pi*idist*(1-pratio)*(delta.^3)));
              k220=(((2*emod2*(1-3*pratio))./(thick*pi*delta*(idist.^3)*(1-pratio.^2))));
              k230=(-(sqrt((2*Gc*emod2*(1-3*pratio))./(1-pratio.^2))./(thick*pi*(idist.^3))));
              k240=(-(sqrt((2*Gc*emod2*(1-3*pratio))./(1-pratio.^2))./(thick*pi*(idist.^3))));
              k320=(-(sqrt((2*Gc*emod2*(1-3*pratio))./(1-pratio.^2))./(thick*pi*(idist.^3))));
              k330=(((Gc2*delta)./(thick*pi*(idist.^3))-(4*Gc2)./(thick*pi*(delta.^2))+(16*G2*lc.^2)./(thick*pi*delta.^4)));
              k340=(((Gc2*delta)./(thick*pi*(idist.^3))-(4*Gc2)./(thick*pi*(delta.^2))-(16*G2*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);
              F11=(k11*theat1.^2+k22*theat2.^2)*v1+(k11*theat1*theat2-k22*theat1*theat2)*v2-k23*theat2*v3-k24*theat2*v4;
              F21=(k11*theat1*theat2-k22*theat1*theat2)*v1+(k11*theat2.^2+k22*theat1.^2)*v2+k23*theat1*v3+k24*theat1*v4;
              F31=k32*(-theat2)*v1+k32*theat1*v2+k33*v3+k34*v4;  
              %加载块的物质点按刚体处理
              F110=(k110*theat1.^2+k220*theat2.^2)*v1+(k110*theat1*theat2-k220*theat1*theat2)*v2-k230*theat2*v3-k240*theat2*v4;
              F210=(k110*theat1*theat2-k220*theat1*theat2)*v1+(k110*theat2.^2+k220*theat1.^2)*v2+k230*theat1*v3+k240*theat1*v4;
              F310=k320*(-theat2)*v1+k320*theat1*v2+k330*v3+k340*v4;  
            if (fail(i,j)==1)%如果这两点间的键没有坏，则有近场力作用
                if i>allnum&&cnode>allnum
                   dforce1=F110* vol * scr * fac;
                   dforce2=F210* vol * scr * fac;
                   dforce3=F310* vol * fac;                      
                end
                if i<=allnum&&cnode<=allnum
                  if eps(i,j)<=-scr3
                       dforce1 = 0.0d0;
                       dforce2 = 0.0d0;
                       dforce3 = 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 * F11 * vol * fac;
                       dforce2 = k3 * F21 * vol * fac;
                       dforce3 = k3 * F31 * 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 * F11 * vol * fac;
                       dforce2 = k2 * F21 * vol * fac;
                       dforce3 = k2 * F31 * vol * fac;
                  elseif eps(i,j)<=scr1&&eps(i,j)>-scr1
                       dforce1 = F11 * vol * fac;
                       dforce2 = F21 * vol * fac;
                       dforce3 = F31 * 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 * F11 * vol * fac;
                       dforce2 = k2 * F21 * vol * fac;
                       dforce3 = k2 * F31 * 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 * F11 * vol * fac;
                       dforce2 = k3 * F21 * vol * fac;
                       dforce3 = k3 * F31 * vol * fac;
                  elseif eps(i,j)>scr3
                       dforce1 = 0.0d0;
                       dforce2 = 0.0d0;
                       dforce3 = 0.0d0;
                  end   
                end
            else%否则没有
                dforce1 = 0.0d0;
                dforce2 = 0.0d0;
                dforce3 = 0.0d0;
            end
            if i<=allnum&&cnode<=allnum
               bondforce(i,j)=sqrt(dforce1^2+dforce2^2+(lc*dforce3/(idist^2))^2);
            end            
            %求力矩阵，将每个点的作用力叠加起来
            pforce(i,1) = pforce(i,1) + dforce1;
            pforce(i,2) = pforce(i,2) + dforce2;
            pforce(i,3) = pforce(i,3) + dforce3;
           %进行键的伸长判断，键伸长达到一定值发生破坏,不去判断那些已经损伤的键
%            if (sqrt(((nlength - idist) / idist) .^2+(lc.*(disp(i,3)-disp(cnode,3))./idist).^2)>scr0)&&(fail(i,j)==1)
%              fail(i,j)=0;%0代表键已经坏了
%              flag(i,j)=0;     
%              fail(j,i)=0;
%              flag(j,i)=0;               
%              sbond=sqrt(((nlength - idist) / idist) .^2+(lc.*(disp(i,3)-disp(cnode,3))./idist));
%              numbond=numbond-1;
%              switcht=1.0d0;   %只要有键被拉坏，则跳过所有疲劳过程
%            end            
           %进行键的寿命判断，键的寿命在疲劳过程中进行迭代

        if (switcht<=0.5)&&(i<=allnum)&&(cnode<=allnum)  %判定是否还进行疲劳迭代过程
           if (tt>=2)&&(fail(i,j)==1)&&abs((nlength-idist)/idist)>einf
              if dmg(i,1)<=0.398
                    lamda(i,j)=lamda(i,j)-(nc(tt,1)-nc(tt-1,1))*a1*((bondforce(i,j)/fk).^(m1));          
              elseif dmg(i,1)>0.398
                  for k=1:nodefam(i,1)
                      if (fail(i,k)==1)&&(flag(i,k)==1) %只有1阶段,并且没有坏的键进入循环
                         lamda(i,k)=1.0d0;
                         flag(i,k)=2.0d0;   %这些键进入2阶段
                      end   
                  end                   
              end 
              if flag(i,j) > 1.5 %二阶段的键进行二阶段迭代
                 lamda(i,j)=lamda(i,j)-(nc(tt,1)-nc(tt-1,1))*a1*((bondforce(i,j)/fk).^(m1));
              end               
           end
        end
         if (switcht<=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
       
%判断键的寿命是否为0，并设置不坏区
%          if switcht<=0.5
             if (lamda(i,j)<=0)&&(fail(i,j)==1)
                 fail(i,j)=0;
                 flag(i,j)=0;
                 fail(j,i)=0;
                 flag(j,i)=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: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,100)==0)%只要时间步长是10的倍数，就进行输出，输出的文件里面有各点的位移和损伤值，若想输出其他的可自己定义。         
         str1=[directpath,'dispdmg',num2str(tt),'.mat'];
         dispdmg=[disp(dotn,:),dmg(dotn,1)];
         save(str1,'dispdmg')                
%% 简单计算应变分布
%   strx=zeros(allnum,3);
%   stry=zeros(allnum,3);
%   for i=1:allnum
%     strx(i,2)=i;
%     stry(i,2)=i;    
%     for j = 1:nodefam(i,1)
%         cnode = nodefam(i,j+1);
%         if (dotall(cnode,1)-dotall(i,1)<0)&&(dotall(cnode,1)-dotall(i,1)>=-1.1d0*dl)
%            strx(i,1)=cnode;
%         end
%         if (dotall(cnode,1)-dotall(i,1)>0)&&(dotall(cnode,1)-dotall(i,1)<=1.1d0*dl)
%            strx(i,3)=cnode;        
%         end
%         if (dotall(cnode,2)-dotall(i,2)<0)&&(dotall(cnode,2)-dotall(i,2)>=-1.1d0*dl)
%            stry(i,1)=cnode;
%         end
%         if (dotall(cnode,2)-dotall(i,2)>0)&&(dotall(cnode,2)-dotall(i,2)<=1.1d0*dl)
%            stry(i,3)=cnode;        
%         end        
%     end  
%   end
%   sttrx=zeros(allnum,1);
%   sttry=zeros(allnum,1);
%   for i=1:allnum 
%     if (strx(i,1)==0)&&(strx(i,3)~=0)
%        u1=disp(strx(i,2),1);
%        u3=disp(strx(i,3),1); 
%        sttrx(i,1)=(u3-u1)/dl;
%     elseif (strx(i,1)~=0)&&(strx(i,3)==0)
%        u2=disp(strx(i,1),1);
%        u1=disp(strx(i,2),1);       
%        sttrx(i,1)=(u1-u2)/dl;
%     elseif (strx(i,1)~=0)&&(strx(i,3)~=0)
%        u2=disp(strx(i,1),1);
%        u1=disp(strx(i,2),1);
%        u3=disp(strx(i,3),1);
%        sttrx(i,1)=0.5d0*((u1-u2)/dl+(u3-u1)/dl);       
%     end    
%   end
%   for i=1:allnum  
%     if (stry(i,1)==0)&&(stry(i,3)~=0)
%        v1=disp(stry(i,2),2);
%        v3=disp(stry(i,3),2);         
%        sttry(i,1)=(v3-v1)/dl;
%     elseif (stry(i,1)~=0)&&(stry(i,3)==0)
%        v2=disp(stry(i,1),2);
%        v1=disp(stry(i,2),2);        
%        sttry(i,1)=(v1-v2)/dl;
%     elseif (stry(i,1)~=0)&&(stry(i,3)~=0)
%        v2=disp(stry(i,1),2);
%        v1=disp(stry(i,2),2);
%        v3=disp(stry(i,3),2);  
%        sttry(i,1)=0.5d0*((v1-v2)/dl+(v3-v1)/dl);
%     end   
%   end
%     str2=[directpath,'sttrx',num2str(tt),'.mat'];
%     sttrx=[sttrx];
%     save(str2,'sttrx')  
%     str3=[directpath,'sttry',num2str(tt),'.mat'];
%     sttry=[sttry];
%     save(str3,'sttry')    
%   figure(2)
%   subplot(1,2,1)
%   scatter_color(dotall(1:allnum,:),sttrx(:,1))
%   subplot(1,2,2)
%   scatter_color(dotall(1:allnum,:),sttry(:,1))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    end
   timebond(tt,1)=numbond;
end
figure (1)
load('G:\fatigue_work004\dispdmg60000')
scatter_color(dotall(1:allnum,:),dispdmg(1:allnum,1))