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 adjustcontactnodes(tieset,ntie,itietri,cg,straight,
20     &  co,vold,xo,yo,zo,x,y,z,nx,ny,nz,istep,iinc,iit,
21     &  mi,imastop,nslavnode,islavnode,set,nset,istartset,
22     &  iendset,ialset,tietol,clearini,clearslavnode,itiefac,
23     &  ipkon,kon,lakon,islavsurf)
24!
25!     Adjusting contact nodes if requested by the user
26!
27      implicit none
28!
29      character*8 lakon(*)
30      character*81 tieset(3,*),slavset,set(*),noset
31!
32      integer ntie,
33     &  itietri(2,ntie),node,neigh(1),kneigh,i,j,k,l,m,isol,iset,
34     &  idummy,itri,ll,kflag,n,nx(*),ny(*),istep,iinc,mi(*),
35     &  nz(*),nstart,iit,imastop(3,*),itriangle(100),ntriangle,
36     &  ntriangle_,itriold,itrinew,id,nslavnode(*),islavnode(*),
37     &  ipos,nset,istartset(*),iendset(*),ialset(*),iclear,
38     &  istart,ilength,nope,nopes,nelems,jj,ifaces,itiefac(2,*),
39     &  jfaces,islavsurf(2,*),ifaceq(8,6),ifacet(6,4),ifacew1(4,5),
40     &  ifacew2(8,5),ipkon(*),kon(*)
41!
42      real*8 cg(3,*),straight(16,*),co(3,*),vold(0:mi(2),*),p(3),
43     &  dist,xo(*),yo(*),zo(*),x(*),y(*),z(*),tietol(3,*),adjust,
44     &  clearini(3,9,*),clearslavnode(3,*),clear
45!
46!     nodes per face for hex elements
47!
48      data ifaceq /4,3,2,1,11,10,9,12,
49     &            5,6,7,8,13,14,15,16,
50     &            1,2,6,5,9,18,13,17,
51     &            2,3,7,6,10,19,14,18,
52     &            3,4,8,7,11,20,15,19,
53     &            4,1,5,8,12,17,16,20/
54!
55!     nodes per face for tet elements
56!
57      data ifacet /1,3,2,7,6,5,
58     &             1,2,4,5,9,8,
59     &             2,3,4,6,10,9,
60     &             1,4,3,8,10,7/
61!
62!     nodes per face for linear wedge elements
63!
64      data ifacew1 /1,3,2,0,
65     &             4,5,6,0,
66     &             1,2,5,4,
67     &             2,3,6,5,
68     &             3,1,4,6/
69!
70!     nodes per face for quadratic wedge elements
71!
72      data ifacew2 /1,3,2,9,8,7,0,0,
73     &             4,5,6,10,11,12,0,0,
74     &             1,2,5,4,7,14,10,13,
75     &             2,3,6,5,8,15,11,14,
76     &             3,1,4,6,9,13,12,15/
77!
78      do i=1,ntie
79         if(tieset(1,i)(81:81).ne.'C') cycle
80         kneigh=1
81         slavset=tieset(2,i)
82!
83!        check whether a clearance has been specified by the user
84!
85c         tietol(3,i)=1.2357111317d0
86         if((tietol(3,i).gt.1.2357111316d0).and.
87     &      (tietol(3,i).lt.1.2357111318d0)) then
88            iclear=0
89         else
90            iclear=1
91            clear=tietol(3,i)
92         endif
93!
94!     check whether an adjust node set has been defined
95!     only checked in the first increment of the first step
96!
97         if(istep.eq.1) then
98            iset=0
99            if(tieset(1,i)(1:1).ne.' ') then
100               noset(1:80)=tieset(1,i)(1:80)
101               noset(81:81)=' '
102               ipos=index(noset,' ')
103               noset(ipos:ipos)='N'
104c               do iset=1,nset
105c                  if(set(iset).eq.noset) exit
106c               enddo
107               call cident81(set,noset,nset,id)
108               iset=nset+1
109               if(id.gt.0) then
110                 if(noset.eq.set(id)) then
111                   iset=id
112                 endif
113               endif
114               kflag=1
115               call isortii(ialset(istartset(iset)),idummy,
116     &              iendset(iset)-istartset(iset)+1,kflag)
117            endif
118         endif
119!
120!     search a master face for each slave node; determine the
121!     distance and adjust it if requested by the user
122!
123         nstart=itietri(1,i)-1
124         n=itietri(2,i)-nstart
125         if(n.lt.kneigh) kneigh=n
126         do j=1,n
127            xo(j)=cg(1,nstart+j)
128            x(j)=xo(j)
129            nx(j)=j
130            yo(j)=cg(2,nstart+j)
131            y(j)=yo(j)
132            ny(j)=j
133            zo(j)=cg(3,nstart+j)
134            z(j)=zo(j)
135            nz(j)=j
136         enddo
137         kflag=2
138         call dsort(x,nx,n,kflag)
139         call dsort(y,ny,n,kflag)
140         call dsort(z,nz,n,kflag)
141!
142         do j=nslavnode(i)+1,nslavnode(i+1)
143            node=islavnode(j)
144!
145            do k=1,3
146               p(k)=co(k,node)
147c               p(k)=co(k,node)+vold(k,node)
148            enddo
149!
150!     determining the kneigh neighboring master contact
151!     triangle centers of gravity
152!
153            call near3d(xo,yo,zo,x,y,z,nx,ny,nz,p(1),p(2),p(3),
154     &           n,neigh,kneigh)
155!
156            isol=0
157!
158            itriold=0
159            itri=neigh(1)+itietri(1,i)-1
160            ntriangle=0
161            ntriangle_=100
162!
163            loop1: do
164               do l=1,3
165                  ll=4*l-3
166                  dist=straight(ll,itri)*p(1)+
167     &                 straight(ll+1,itri)*p(2)+
168     &                 straight(ll+2,itri)*p(3)+
169     &                 straight(ll+3,itri)
170                  if(dist.gt.1.d-6) then
171                     itrinew=imastop(l,itri)
172                     if(itrinew.eq.0) then
173c     write(*,*) '**border reached'
174                        exit loop1
175                     elseif(itrinew.eq.itriold) then
176c     write(*,*) '**solution in between triangles'
177                        isol=itri
178                        exit loop1
179                     else
180                        call nident(itriangle,itrinew,ntriangle,id)
181                        if(id.gt.0) then
182                           if(itriangle(id).eq.itrinew) then
183c     write(*,*) '**circular path; no solution'
184                              exit loop1
185                           endif
186                        endif
187                        ntriangle=ntriangle+1
188                        if(ntriangle.gt.ntriangle_) then
189c     write(*,*) '**too many iterations'
190                           exit loop1
191                        endif
192                        do k=ntriangle,id+2,-1
193                           itriangle(k)=itriangle(k-1)
194                        enddo
195                        itriangle(id+1)=itrinew
196                        itriold=itri
197                        itri=itrinew
198                        cycle loop1
199                     endif
200                  elseif(l.eq.3) then
201c     write(*,*) '**regular solution'
202                     isol=itri
203                     exit loop1
204                  endif
205               enddo
206            enddo loop1
207!
208!     calculating the distance
209!
210            if(isol.ne.0) then
211               dist=straight(13,itri)*p(1)+
212     &              straight(14,itri)*p(2)+
213     &              straight(15,itri)*p(3)+
214     &              straight(16,itri)
215!
216!     check for an adjust parameter (only in the first
217!     increment of the first step)
218!
219               if(istep.eq.1) then
220                  if(iset.ne.0) then
221!
222!     check whether node belongs to the adjust node
223!     set
224!
225                     call nident(ialset(istartset(iset)),node,
226     &                    iendset(iset)-istartset(iset)+1,id)
227                     if(id.gt.0) then
228                        if(ialset(istartset(iset)+id-1).eq.node) then
229                           do k=1,3
230                              co(k,node)=co(k,node)-
231     &                             dist*straight(12+k,itri)
232                           enddo
233                           dist=0.d0
234                        endif
235                     endif
236                  elseif(dabs(tietol(1,i)).ge.2.d0) then
237!
238!     adjust parameter
239!
240                     adjust=dabs(tietol(1,i))-2.d0
241
242                     if(dist.le.adjust) then
243                        do k=1,3
244                           co(k,node)=co(k,node)-
245     &                          dist*straight(12+k,itri)
246                        enddo
247                        dist=0.d0
248                     endif
249                  endif
250               endif
251!
252!              clearance in each slave node
253!
254               if(iclear.eq.1) then
255                  do k=1,3
256                     clearslavnode(k,j)=(clear-dist)*straight(12+k,itri)
257                  enddo
258               endif
259!
260            endif
261         enddo
262!
263!        loop over all slave faces
264!
265         if(iclear.eq.1) then
266            do jj=itiefac(1,i), itiefac(2,i)
267               ifaces=islavsurf(1,jj)
268               nelems=int(ifaces/10)
269               jfaces=ifaces - nelems*10
270!
271               if(lakon(nelems)(4:5).eq.'8R') then
272                  nopes=4
273                  nope=8
274               elseif(lakon(nelems)(4:4).eq.'8') then
275                  nopes=4
276                  nope=8
277               elseif(lakon(nelems)(4:6).eq.'20R') then
278                  nopes=8
279                  nope=20
280               elseif(lakon(nelems)(4:5).eq.'20') then
281                  nopes=8
282                  nope=20
283               elseif(lakon(nelems)(4:5).eq.'10') then
284                  nopes=6
285                  nope=10
286               elseif(lakon(nelems)(4:4).eq.'4') then
287                  nopes=3
288                  nope=4
289!
290!     treatment of wedge faces
291!
292               elseif(lakon(nelems)(4:4).eq.'6') then
293                  nope=6
294                  if(jfaces.le.2) then
295                     nopes=3
296                  else
297                     nopes=4
298                  endif
299               elseif(lakon(nelems)(4:5).eq.'15') then
300                  nope=15
301                  if(jfaces.le.2) then
302                     nopes=6
303                  else
304                     nopes=8
305                  endif
306               endif
307!
308!     actual position of the nodes belonging to the
309!     slave surface
310!
311               istart=nslavnode(i)+1
312               ilength=nslavnode(i+1)-nslavnode(i)
313!
314               do m=1,nopes
315                  if((nope.eq.20).or.(nope.eq.8))then
316                     node=kon(ipkon(nelems)+ifaceq(m,jfaces))
317                  elseif((nope.eq.10).or.(nope.eq.4).or.(nope.eq.14))
318     &                    then
319                     node=kon(ipkon(nelems)+ifacet(m,jfaces))
320                  elseif(nope.eq.15) then
321                     node=kon(ipkon(nelems)+ifacew2(m,jfaces))
322                  else
323                     node=kon(ipkon(nelems)+ifacew1(m,jfaces))
324                  endif
325                  call nident(islavnode(istart),node,ilength,id)
326                  do k=1,3
327                     clearini(k,m,jj)=
328     &                    clearslavnode(k,nslavnode(i)+id)
329                  enddo
330               enddo
331            enddo
332         endif
333!
334!     transfer the clearance from the slave nodes in islavnode
335!     to the nodes per slave face in islavsurf
336!
337      enddo
338!
339      return
340      end
341
342