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!
20!     subroutine generating the first actice set
21!     and calculating the areas for all faces on slaveside
22!
23!     [in]  imastop        	(l,i) for edge l in triagle i neightbouring triangle
24!     [in]  islavsurf      	islavsurf(1,i) slaveface i islavsurf(2,i) pointer into imastsurf and pmastsurf
25!     [in]  itiefac     	pointer into field islavsurf: (1,i) beginning slave_i (2,i) end of slave_i
26!     [out]    areaslav    	(i)area of face i, ONLY HELP FIELD
27!     [in,out] ifree       	in=0 out=# active nodes
28!
29      subroutine genfirstactif(tieset,ntie,itietri,ipkon,kon,lakon,cg,
30     &     straight,co,vold,xo,yo,zo,x,y,z,nx,ny,nz,istep,iinc,iit,
31     &     mi,imastop,nslavnode,islavnode,islavsurf,itiefac,areaslav,
32     &     set,nset,istartset,iendset,ialset,islavact,ifree,tietol)
33!
34!     Initialization of the Active slave nodes set
35!
36!     Author: Samoela Rakotonanahary, Saskia Sitzmann
37!
38!     guido: areaslav not needed??
39!
40      implicit none
41!
42      character*8 lakon(*)
43      character*81 tieset(3,*),slavset,set(*),noset
44!
45      integer ntie,
46     &     itietri(2,ntie),ipkon(*),kon(*),node,neigh(1),iflag,kneigh,
47     &     i,j,k,l,isol,iset,idummy,itri,ll,kflag,n,nx(*),ny(*),istep,
48     &     iinc,mi(*),nz(*),nstart,iit,nope,iteller,ifaces,jfaces,
49     &     imastop(3,*), itriangle(100),ntriangle,ntriangle_,itriold,
50     &     itrinew,id,nslavnode(*),islavnode(*),islavsurf(2,*),
51     &     itiefac(2,*),konl(20),nelems,m,mint2d,nopes,
52     &     ipos,nset,istartset(*),iendset(*),
53     &     ialset(*),islavact(*),ifree,ifac,getlocno
54!
55      real*8 cg(3,*),straight(16,*),co(3,*),vold(0:mi(2),*),p(3),
56     &     dist,xo(*),yo(*),zo(*),x(*),y(*),z(*),c0,weight,
57     &     areaslav(*),xl2(3,8),area,xi,et,shp2(7,8),
58     &     xs2(3,2),xsj2(3),tietol(3,*),adjust
59!
60!     flag for shape functions
61!
62      data iflag /2/
63!
64      data iteller /0/
65      save iteller
66!
67      include "gauss.f"
68!
69      do i=1,ntie
70        if(tieset(1,i)(81:81).ne.'C') cycle
71        kneigh=1
72        slavset=tieset(2,i)
73!
74!     check whether an adjust node set has been defined
75!     only checked in the first increment of the first step
76!
77        if((istep.eq.1).and.(iinc.eq.1).and.(iit.le.1)) then
78          iset=0
79          if(tieset(1,i)(1:1).ne.' ') then
80            noset(1:80)=tieset(1,i)(1:80)
81            noset(81:81)=' '
82            ipos=index(noset,' ')
83            noset(ipos:ipos)='N'
84c            do iset=1,nset
85c              if(set(iset).eq.noset) exit
86c            enddo
87            call cident81(set,noset,nset,id)
88            iset=nset+1
89            if(id.gt.0) then
90              if(noset.eq.set(id)) then
91                iset=id
92              endif
93            endif
94            kflag=1
95            call isortii(ialset(istartset(iset)),idummy,
96     &           iendset(iset)-istartset(iset)+1,kflag)
97          endif
98        endif
99!
100!     determine the area of the slave surfaces
101!
102        do l=itiefac(1,i),itiefac(2,i)
103          ifaces=islavsurf(1,l)
104          nelems=int(ifaces/10)
105          jfaces=ifaces-nelems*10
106!
107!     Decide on the max integration points number, just consider 2D situation
108!
109          call getnumberofnodes(nelems,jfaces,lakon,nope,nopes,mint2d)
110!
111!     actual position of the nodes belonging to the
112!     slave surface
113!
114          do j=1,nope
115            konl(j)=kon(ipkon(nelems)+j)
116          enddo
117          do m=1,nopes
118            ifac=getlocno(m,jfaces,nope)
119            do j=1,3
120              xl2(j,m)=co(j,konl(ifac))+
121     &             vold(j,konl(ifac))
122            enddo
123          enddo
124!
125!     calculating the area of the slave face
126!
127          area=0.d0
128          do m=1,mint2d
129            if((lakon(nelems)(4:5).eq.'8R').or.
130     &           ((lakon(nelems)(4:4).eq.'6').and.(nopes.eq.4))) then
131              xi=gauss2d1(1,m)
132              et=gauss2d1(2,m)
133              weight=weight2d1(m)
134            elseif((lakon(nelems)(4:4).eq.'8').or.
135     &             (lakon(nelems)(4:6).eq.'20R').or.
136     &             ((lakon(nelems)(4:5).eq.'15').and.
137     &             (nopes.eq.8))) then
138              xi=gauss2d2(1,m)
139              et=gauss2d2(2,m)
140              weight=weight2d2(m)
141            elseif(lakon(nelems)(4:4).eq.'2') then
142              xi=gauss2d3(1,m)
143              et=gauss2d3(2,m)
144              weight=weight2d3(m)
145            elseif((lakon(nelems)(4:5).eq.'10').or.
146     &             ((lakon(nelems)(4:5).eq.'15').and.
147     &             (nopes.eq.6))) then
148              xi=gauss2d5(1,m)
149              et=gauss2d5(2,m)
150              weight=weight2d5(m)
151            elseif((lakon(nelems)(4:4).eq.'4').or.
152     &             ((lakon(nelems)(4:4).eq.'6').and.
153     &             (nopes.eq.3))) then
154              xi=gauss2d4(1,m)
155              et=gauss2d4(2,m)
156              weight=weight2d4(m)
157            endif
158!
159            if(nopes.eq.8) then
160              call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
161            elseif(nopes.eq.4) then
162              call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag)
163            elseif(nopes.eq.6) then
164              call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
165            else
166              call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
167            endif
168            area=area+weight*dsqrt(xsj2(1)**2+xsj2(2)**2+
169     &           xsj2(3)**2)
170          enddo
171          areaslav(l)=area
172        enddo
173!
174!     search a master face for each slave node and generate a contact
175!     spring element if successful
176!
177        nstart=itietri(1,i)-1
178        n=itietri(2,i)-nstart
179        if(n.lt.kneigh) kneigh=n
180        do j=1,n
181          xo(j)=cg(1,nstart+j)
182          x(j)=xo(j)
183          nx(j)=j
184          yo(j)=cg(2,nstart+j)
185          y(j)=yo(j)
186          ny(j)=j
187          zo(j)=cg(3,nstart+j)
188          z(j)=zo(j)
189          nz(j)=j
190        enddo
191        kflag=2
192        call dsort(x,nx,n,kflag)
193        call dsort(y,ny,n,kflag)
194        call dsort(z,nz,n,kflag)
195!
196        do j=nslavnode(i)+1,nslavnode(i+1)
197          node=islavnode(j)
198!
199          do k=1,3
200            p(k)=co(k,node)+vold(k,node)
201          enddo
202!
203!     determining the kneigh neighboring master contact
204!     triangle centers of gravity
205!
206          call near3d(xo,yo,zo,x,y,z,nx,ny,nz,p(1),p(2),p(3),
207     &         n,neigh,kneigh)
208!
209          isol=0
210!
211          itriold=0
212          itri=neigh(1)+itietri(1,i)-1
213          ntriangle=0
214          ntriangle_=100
215!
216          loop1: do
217            do l=1,3
218              ll=4*l-3
219              dist=straight(ll,itri)*p(1)+
220     &             straight(ll+1,itri)*p(2)+
221     &             straight(ll+2,itri)*p(3)+
222     &             straight(ll+3,itri)
223              if(dist.gt.1.d-6) then
224                itrinew=imastop(l,itri)
225                if(itrinew.eq.0) then
226c     write(*,*) '**border reached'
227                  exit loop1
228                elseif((itrinew.lt.itietri(1,i)).or.
229     &                 (itrinew.gt.itietri(2,i))) then
230c     write(*,*) '**border reached'
231                  exit loop1
232                elseif(itrinew.eq.itriold) then
233c     write(*,*) '**solution in between triangles'
234                  isol=itri
235                  exit loop1
236                else
237                  call nident(itriangle,itrinew,ntriangle,id)
238                  if(id.gt.0) then
239                    if(itriangle(id).eq.itrinew) then
240c     write(*,*) '**circular path; no solution'
241                      exit loop1
242                    endif
243                  endif
244                  ntriangle=ntriangle+1
245                  if(ntriangle.gt.ntriangle_) then
246c     write(*,*) '**too many iterations'
247                    exit loop1
248                  endif
249                  do k=ntriangle,id+2,-1
250                    itriangle(k)=itriangle(k-1)
251                  enddo
252                  itriangle(id+1)=itrinew
253                  itriold=itri
254                  itri=itrinew
255                  cycle loop1
256                endif
257              elseif(l.eq.3) then
258c     write(*,*) '**regular solution'
259                isol=itri
260                exit loop1
261              endif
262            enddo
263          enddo loop1
264!
265!     check whether distance is larger than c0:
266!     no element is generated
267!
268          if(isol.ne.0) then
269            dist=straight(13,itri)*p(1)+
270     &           straight(14,itri)*p(2)+
271     &           straight(15,itri)*p(3)+
272     &           straight(16,itri)
273!
274!     check for an adjust parameter (only in the first
275!     increment of the first step)
276!
277            if((istep.eq.1).and.(iinc.eq.1).and.(iit.le.1)) then
278              if(iset.ne.0) then
279!
280!     check whether node belongs to the adjust node
281!     set
282!
283                call nident(ialset(istartset(iset)),node,
284     &               iendset(iset)-istartset(iset)+1,id)
285                if(id.gt.0) then
286                  if(ialset(istartset(iset)+id-1).eq.node) then
287                    do k=1,3
288                      co(k,node)=co(k,node)-
289     &                     dist*straight(12+k,itri)
290                    enddo
291                    dist=0.d0
292                  endif
293                endif
294              elseif(dabs(tietol(1,i)).ge.2.d0) then
295!
296!     adjust parameter
297!
298                adjust=dabs(tietol(1,i))-2.d0
299
300                if(dist.le.adjust) then
301                  do k=1,3
302                    co(k,node)=co(k,node)-
303     &                   dist*straight(12+k,itri)
304                  enddo
305                  dist=0.d0
306                endif
307              endif
308            endif
309!
310            c0=1.d-10
311            if(dabs(tietol(1,i)).ge.2.d0) then
312              c0=dabs(tietol(1,i))-2.d0
313            endif
314            if(dist.gt.c0) then
315              isol=0
316!
317!     adjusting the bodies at the start of the
318!     calculation such that they touch
319!
320            endif
321          endif
322!
323          if(isol.ne.0) then
324!
325!     Active node
326            islavact(j)=2
327            ifree=ifree+1
328          else
329            islavact(j)=-1
330!
331          endif
332!
333        enddo
334      enddo
335!
336      return
337      end
338