1! 2! CalculiX - A 3-dimensional finite element program 3! Copyright (C) 1998-2021 Guido Dhondt 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU General Public License as 7! published by the Free Software Foundation(version 2); 8! 9! 10! This program is distributed in the hope that it will be useful, 11! but WITHOUT ANY WARRANTY; without even the implied warranty of 12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13! GNU General Public License for more details. 14! 15! You should have received a copy of the GNU General Public License 16! along with this program; if not, write to the Free Software 17! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 18! 19 subroutine cracklength_smoothing(nnfront,isubsurffront, 20 & istartfront,iendfront,ifrontrel,costruc,dist,a, 21 & istartcrackfro,iendcrackfro,acrack,amin,nfront,ncrack, 22 & idist,nstep) 23! 24! smoothes the crack length 25! calculates the minimum crack length for each crack 26! 27 implicit none 28! 29 integer i,j,nnfront,isubsurffront(*),istart,iend,istartfront(*), 30 & iendfront(*),ifrontrel(*),jforward,jbackward,n,nmax,node, 31 & istartcrackfro(*),iendcrackfro(*),nfront,ncrack,k,idist(*), 32 & nstep,m 33! 34 real*8 xp,yp,zp,costruc(3,*),distforward,distbackward,dist(*), 35 & radius,sum,temp,acrack(nstep,*),amin(*),alpha,a(nstep,*) 36! 37! copy the crack length in temporary field a 38! 39 do i=1,nfront 40 do m=1,nstep 41 a(m,i)=acrack(m,i) 42 enddo 43 enddo 44! 45 do i=1,nnfront 46! 47! target number of nodes in the circle of influence 48! 49 nmax=max(1,int((iendfront(i)-istartfront(i)+1)/1.1d0)) 50 if(nmax.le.2) cycle 51! 52 if(isubsurffront(i).eq.1) then 53 istart=istartfront(i) 54 iend=iendfront(i) 55 else 56 istart=istartfront(i)+1 57 iend=iendfront(i)-1 58 endif 59! 60 do j=istart,iend 61 xp=costruc(1,ifrontrel(j)) 62 yp=costruc(2,ifrontrel(j)) 63 zp=costruc(3,ifrontrel(j)) 64! 65! actual forward node and its distance from node 66! actual backward node and its distance from node 67! actual number of nodes n in the circle of influence 68! 69 jforward=j 70 jbackward=j 71 n=1 72 dist(j)=0.d0 73 idist(1)=j 74! 75! calculating one position forward 76! 77 jforward=jforward+1 78 if(isubsurffront(i).eq.1) then 79 if(jforward.gt.iendfront(i)) jforward=istartfront(i) 80 endif 81! 82 if((isubsurffront(i).eq.0).and. 83 & (jforward.gt.iendfront(i))) then 84 distforward=1.d30 85 else 86 distforward= 87 & dsqrt((costruc(1,ifrontrel(jforward))-xp)**2+ 88 & (costruc(2,ifrontrel(jforward))-yp)**2+ 89 & (costruc(3,ifrontrel(jforward))-zp)**2) 90 endif 91! 92! calculating one position backward 93! 94 jbackward=jbackward-1 95 if(isubsurffront(i).eq.1) then 96 if(jbackward.lt.istartfront(i)) jbackward=iendfront(i) 97 endif 98! 99 if((isubsurffront(i).eq.0).and. 100 & (jbackward.lt.istartfront(i))) then 101 distbackward=1.d30 102 else 103 distbackward= 104 & dsqrt((costruc(1,ifrontrel(jbackward))-xp)**2+ 105 & (costruc(2,ifrontrel(jbackward))-yp)**2+ 106 & (costruc(3,ifrontrel(jbackward))-zp)**2) 107 endif 108! 109! start infinite loop 110! 111 do 112 if(distforward.lt.distbackward) then 113 n=n+1 114 idist(n)=jforward 115 dist(jforward)=distforward 116 radius=distforward 117 if(n.eq.nmax) exit 118! 119 jforward=jforward+1 120! 121! calculate the next forward position 122! 123 if(isubsurffront(i).eq.1) then 124 if(jforward.gt.iendfront(i)) jforward=istartfront(i) 125 endif 126! 127 if((isubsurffront(i).eq.0).and. 128 & (jforward.gt.iendfront(i))) then 129 distforward=1.d30 130 else 131 distforward= 132 & dsqrt((costruc(1,ifrontrel(jforward))-xp)**2+ 133 & (costruc(2,ifrontrel(jforward))-yp)**2+ 134 & (costruc(3,ifrontrel(jforward))-zp)**2) 135 endif 136 else 137 n=n+1 138 idist(n)=jbackward 139 dist(jbackward)=distbackward 140 radius=distbackward 141 if(n.eq.nmax) exit 142! 143 jbackward=jbackward-1 144! 145! calculate the next backward position 146! 147 if(isubsurffront(i).eq.1) then 148 if(jbackward.lt.istartfront(i)) jbackward=iendfront(i) 149 endif 150! 151 if((isubsurffront(i).eq.0).and. 152 & (jbackward.lt.istartfront(i))) then 153 distbackward=1.d30 154 else 155 distbackward= 156 & dsqrt((costruc(1,ifrontrel(jbackward))-xp)**2+ 157 & (costruc(2,ifrontrel(jbackward))-yp)**2+ 158 & (costruc(3,ifrontrel(jbackward))-zp)**2) 159 endif 160 endif 161 enddo 162! 163! calculate the weighted mean 164! 165 do m=1,nstep 166 sum=0.d0 167 temp=0.d0 168 do k=1,n 169 node=idist(k) 170 alpha=1.d0-dist(node)/radius 171 sum=sum+alpha 172 temp=temp+alpha*a(m,node) 173 enddo 174 acrack(m,j)=temp/sum 175 enddo 176 enddo 177! 178 if(isubsurffront(i).eq.0) then 179 do m=1,nstep 180 acrack(m,istartfront(i))=acrack(m,istartfront(i)+1) 181 acrack(m,iendfront(i))=acrack(m,iendfront(i)-1) 182 enddo 183 endif 184 enddo 185! 186! determine the minimum crack length for each crack 187! 188c write(*,*) 189c write(*,*) 'cracklength after smoothing' 190c write(*,*) 191c do i=1,ncrack 192c amin(i)=1.d30 193c do j=istartcrackfro(i),iendcrackfro(i) 194c do m=1,nstep 195c amin(i)=min(amin(i),acrack(m,j)) 196c write(*,*) 'cracklength ',j,m,acrack(m,j) 197c enddo 198c enddo 199c enddo 200! 201 return 202 end 203