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