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 adjacentbounodes(ifront,ifrontrel,nfront,iedno,iedg,
20     &     nnfront,ifront2,ifrontrel2,ibounnod,nbounnod,istartfront,
21     &     iendfront,ibounedg,istartcrackfro,iendcrackfro,ncrack,
22     &     ibounnod2,istartcrackbou,iendcrackbou,isubsurffront,iedno2,
23     &     stress,stress2,iresort,ieled,kontri,costruc,costruc2,temp,
24     &     temp2,nstep,ier)
25!
26!     sorting the crack boundary nodes according to adjacency
27!
28      implicit none
29!
30      integer ifront(*),ifrontrel(*),nfront,iedno(2,*),iedg(3,*),
31     &     ichange,nnfront,nfront2,node,noderel,iedge,ibounnod(*),
32     &     nbounnod,i,j,k,istartfront(*),iendfront(*),ifront2(*),
33     &     ifrontrel2(*),id,ibounedg(*),nodestart,istartcrackfro(*),
34     &     iendcrackfro(*),ncrack,iact,nodeprev,noderelprev,kflag,
35     &     ibounnod2(*),istartcrackbou(*),iendcrackbou(*),nbounnod2,
36     &     isubsurffront(*),isubsurf,iedno2(2,*),iresort(*),itri,
37     &     nodenext,ieled(2,*),kontri(3,*),nstep,ier
38!
39      real*8 stress(6,nstep,*),stress2(6,nstep,*),costruc(3,*),
40     &     costruc2(3,*),temp(nstep,*),temp2(nstep,*)
41!
42!     nnfront is the number of distinct fronts
43!
44      nnfront=0
45      nfront2=0
46      nbounnod2=0
47      ncrack=0
48      kflag=2
49c      do i=1,nfront
50c        write(*,*) 'adjacentbounodes ',i,ifront(i)
51c      enddo
52!
53      do
54!
55!     for all recorded front nodes i the entry in ifrontrel(i) is
56!     set to 0. If no zero is left, the job is done.
57!
58        ichange=0
59!
60        loop: do i=1,nfront
61        if(ifrontrel(i).gt.0) then
62          ichange=1
63!
64!     find the end of a front: start with an arbitrary node
65!
66c          write(*,*) 'look for new front '
67          node=ifront(i)
68c          write(*,*) node
69          nodestart=node
70          noderel=ifrontrel(i)
71!
72!     take one of the edges to which the node belongs
73!     check whether the nodes belonging to the edge are in the
74!     same order as in the element topology of the crack
75!
76          iedge=ibounedg(iedno(1,noderel))
77!
78!     crack triangle belonging to the edge
79!
80          itri=ieled(1,iedge)
81!
82!     finding the other node of the edge
83!
84          if(iedg(1,iedge).ne.node) then
85            nodenext=iedg(1,iedge)
86          else
87            nodenext=iedg(2,iedge)
88          endif
89!
90!     node and nodenext also belong to itri; check whether the
91!     order in the element topology is nodenext -> node
92!     (because the order is reverted as soon as an external node
93!     is found (subsurface crack) or the start node is retreated
94!     (surface crack)
95!
96          do j=1,3
97            if(kontri(j,itri).eq.node) exit
98          enddo
99          k=j+1
100          if(k.gt.3) k=1
101c          k=mod(j+1,3)
102!
103!     take the other edge is the order is not the same
104!
105          if(kontri(k,itri).eq.nodenext) then
106            iedge=ibounedg(iedno(2,noderel))
107          endif
108!
109          do
110!
111!     find the other node on the edge
112!
113            if(iedg(1,iedge).ne.node) then
114              node=iedg(1,iedge)
115            else
116              node=iedg(2,iedge)
117            endif
118c            write(*,*) node
119            if(node.eq.nodestart) then
120!
121!     subsurface crack
122!
123              call nident(ifront,node,nfront,id)
124              ifrontrel(id)=0
125              isubsurf=1
126              exit
127            endif
128!
129!     check whether this is a front node
130!
131            call nident(ifront,node,nfront,id)
132            if(id.gt.0) then
133              if(ifront(id).eq.node) then
134!
135!     front node! search for other edge adjacent to the node
136!
137                noderel=ifrontrel(id)
138                if(ibounedg(iedno(1,noderel)).eq.iedge) then
139                  iedge=ibounedg(iedno(2,noderel))
140                else
141                  iedge=ibounedg(iedno(1,noderel))
142                endif
143                cycle
144              endif
145            endif
146!
147            isubsurf=0
148            exit
149          enddo
150!
151!     if surface crack:
152!     node is a boundary node which is external to the structure
153!     but adjacent to a front node (i.e. a boundary node internal
154!     to the structure)
155!
156!     if subsurface crack:
157!     node is an arbitary node belonging to the crack front
158!
159!     finding the relative node number
160!
161          call nident(ibounnod,node,nbounnod,noderel)
162c          write(*,*) 'adjacentbounodes start of new front'
163!
164!     defining a new crack
165!
166          ncrack=ncrack+1
167          nodestart=node
168!
169!     defining a new front segment
170!
171          nnfront=nnfront+1
172          iact=1
173!
174!     new node belonging to a segment (a segment contains the internal
175!     nodes plus the immediately adjacent external nodes
176!
177          nfront2=nfront2+1
178          ifront2(nfront2)=node
179!
180          nbounnod2=nbounnod2+1
181          ibounnod2(nbounnod2)=node
182!
183c          write(*,*) node
184          ifrontrel2(nfront2)=noderel
185          istartfront(nnfront)=nfront2
186          isubsurffront(nnfront)=isubsurf
187          istartcrackfro(ncrack)=nfront2
188          istartcrackbou(ncrack)=nbounnod2
189!
190!     looking for other nodes belonging to the segment
191!
192          do
193!
194!     find the other node on the edge
195!
196            if(iedg(1,iedge).ne.node) then
197              nodeprev=node
198              noderelprev=noderel
199              node=iedg(1,iedge)
200            else
201              nodeprev=node
202              noderelprev=noderel
203              node=iedg(2,iedge)
204            endif
205!
206!     reached starting node
207!
208            if(node.eq.nodestart) then
209              iendcrackbou(ncrack)=nbounnod2
210              iendfront(nnfront)=nfront2
211              exit
212            endif
213!
214!     check whether this is a front node
215!
216            call nident(ifront,node,nfront,id)
217            if(id.gt.0) then
218              if(ifront(id).eq.node) then
219!
220!     front node !
221!
222!     check whether new front is to be activated
223!
224                if(iact.eq.0) then
225!
226!     new front
227!
228                  iact=1
229                  nnfront=nnfront+1
230                  isubsurffront(nnfront)=isubsurf
231!
232!     first node of new front
233!
234                  nfront2=nfront2+1
235                  ifront2(nfront2)=nodeprev
236                  ifrontrel2(nfront2)=noderelprev
237                  istartfront(nnfront)=nfront2
238                endif
239!
240                noderel=ifrontrel(id)
241                if(ibounedg(iedno(1,noderel)).eq.iedge) then
242                  iedge=ibounedg(iedno(2,noderel))
243                else
244                  iedge=ibounedg(iedno(1,noderel))
245                endif
246!
247                nfront2=nfront2+1
248                ifront2(nfront2)=node
249!
250                nbounnod2=nbounnod2+1
251                ibounnod2(nbounnod2)=node
252!
253c                write(*,*) node
254                ifrontrel2(nfront2)=noderel
255!
256!     setting the entry in frontrel for a node which was treated
257!     to zero
258!
259                ifrontrel(id)=0
260                cycle
261              endif
262            endif
263!
264!     no front node!
265!
266            call nident(ibounnod,node,nbounnod,noderel)
267            if(ibounedg(iedno(1,noderel)).eq.iedge) then
268              iedge=ibounedg(iedno(2,noderel))
269            else
270              iedge=ibounedg(iedno(1,noderel))
271            endif
272!
273            nbounnod2=nbounnod2+1
274            ibounnod2(nbounnod2)=node
275!
276            if(iact.eq.1) then
277!
278!     last node of the segment (external boundary node immediately
279!     adjacent to an internal node)
280!
281              nfront2=nfront2+1
282              ifront2(nfront2)=node
283c              write(*,*) node
284              ifrontrel2(nfront2)=noderel
285!
286!     finishing the actual segment
287!
288              iendfront(nnfront)=nfront2
289              iact=0
290            endif
291          enddo
292!
293!     end of crack
294!
295          iendcrackfro(ncrack)=iendfront(nnfront)
296        endif
297      enddo loop                ! end loop over the front nodes
298!
299!     if no front node is left untreated, leave
300!
301      if(ichange.eq.0) exit
302      enddo
303!
304!     adjust the dependent fields of ibounnod
305!
306      do i=1,nbounnod2
307        node=ibounnod2(i)
308        call nident(ibounnod,node,nbounnod,id)
309        do j=1,2
310          iedno2(j,i)=iedno(j,id)
311        enddo
312        do k=1,nstep
313          temp2(k,i)=temp(k,id)
314        enddo
315        do k=1,nstep
316          do j=1,6
317            stress2(j,k,i)=stress(j,k,id)
318          enddo
319        enddo
320        do j=1,3
321          costruc2(j,i)=costruc(j,id)
322        enddo
323      enddo
324!
325!     copy the ibounnod2 field and its dependent fields
326!
327      do i=1,nbounnod2
328        ibounnod(i)=ibounnod2(i)
329        do j=1,2
330          iedno(j,i)=iedno2(j,i)
331        enddo
332        do k=1,nstep
333          temp(k,i)=temp2(k,i)
334        enddo
335        do k=1,nstep
336          do j=1,6
337            stress(j,k,i)=stress2(j,k,i)
338          enddo
339        enddo
340        do j=1,3
341          costruc(j,i)=costruc2(j,i)
342        enddo
343      enddo
344!
345!     calculating in which order ibounnod was resorted
346!
347      do i=1,nbounnod2
348        iresort(i)=i
349      enddo
350!
351      nfront=nfront
352      call isortii(ibounnod2,iresort,nbounnod2,kflag)
353!
354!     copy the temporary fields (ifront2 and ifrontrel2)
355!
356      do i=1,nfront2
357        ifront(i)=ifront2(i)
358        ifrontrel(i)=iresort(ifrontrel2(i))
359      enddo
360      nfront=nfront2
361!
362c      write(*,*)
363c      write(*,*) 'adjacentbounodes: front nodes'
364c      write(*,*)
365c      do i=1,nfront
366c        write(*,*) i,ifront(i)
367c      enddo
368!
369      if(nfront.eq.0) then
370        write(*,*) '*ERROR in adjacentbounodes: no front node found'
371        ier=1
372c         call exit(201)
373       endif
374!
375      return
376      end
377
378