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!   Sets the value 1 to field islavact for nodes which have
20!   opposite master face.
21!
22      subroutine islavactive(tieset,ntie,itietri,
23     &  cg,straight,co,vold,xo,yo,zo,x,y,z,nx,ny,nz,
24     &  mi,imastop,nslavnode,islavnode,islavact)
25!
26      implicit none
27!
28      character*81 tieset(3,*)
29!
30      integer ntie,itietri(2,*),node,neigh(1),iflag,kneigh,i,j,k,l,
31     &  isol,itri,ll,kflag,n,nx(*),ny(*),mi(*),nz(*),nstart,
32     &  ifacew1(4,5),ifacew2(8,5),imastop(3,*),
33     &  itriangle(100),ntriangle,ntriangle_,itriold,itrinew,id,
34     &  nslavnode(*),islavnode(*),islavact(*),ifaceq(8,6),
35     &  ifacet(6,4)
36!
37      real*8 cg(3,*),straight(16,*),co(3,*),vold(0:mi(2),*),p(3),
38     &  dist,xo(*),yo(*),zo(*),x(*),y(*),z(*)
39!
40!     nodes per face for hex elements
41!
42      data ifaceq /4,3,2,1,11,10,9,12,
43     &            5,6,7,8,13,14,15,16,
44     &            1,2,6,5,9,18,13,17,
45     &            2,3,7,6,10,19,14,18,
46     &            3,4,8,7,11,20,15,19,
47     &            4,1,5,8,12,17,16,20/
48!
49!     nodes per face for tet elements
50!
51      data ifacet /1,3,2,7,6,5,
52     &             1,2,4,5,9,8,
53     &             2,3,4,6,10,9,
54     &             1,4,3,8,10,7/
55!
56!     nodes per face for linear wedge elements
57!
58      data ifacew1 /1,3,2,0,
59     &             4,5,6,0,
60     &             1,2,5,4,
61     &             2,3,6,5,
62     &             3,1,4,6/
63!
64!     nodes per face for quadratic wedge elements
65!
66      data ifacew2 /1,3,2,9,8,7,0,0,
67     &             4,5,6,10,11,12,0,0,
68     &             1,2,5,4,7,14,10,13,
69     &             2,3,6,5,8,15,11,14,
70     &             3,1,4,6,9,13,12,15/
71!
72      data iflag /2/
73!
74!
75!
76!     ***ISLAVACT***
77!
78      do i=1,ntie
79         if(tieset(1,i)(81:81).ne.'C') cycle
80         kneigh=1
81!
82!        search a master face for each slave node
83!
84         nstart=itietri(1,i)-1
85         n=itietri(2,i)-nstart
86         if(n.lt.kneigh) kneigh=n
87         do j=1,n
88            xo(j)=cg(1,nstart+j)
89            x(j)=xo(j)
90            nx(j)=j
91            yo(j)=cg(2,nstart+j)
92            y(j)=yo(j)
93            ny(j)=j
94            zo(j)=cg(3,nstart+j)
95            z(j)=zo(j)
96            nz(j)=j
97         enddo
98         kflag=2
99         call dsort(x,nx,n,kflag)
100         call dsort(y,ny,n,kflag)
101         call dsort(z,nz,n,kflag)
102!
103c         write(*,*) 'islavactive ntie= ',i
104         do j=nslavnode(i)+1,nslavnode(i+1)
105            node=islavnode(j)
106!
107            do k=1,3
108               p(k)=co(k,node)+vold(k,node)
109            enddo
110!
111!     determining the kneigh neighboring master contact
112!     triangle centers of gravity
113!
114c            write(*,*) 'islavactive j= ',j
115            call near3d(xo,yo,zo,x,y,z,nx,ny,nz,p(1),p(2),p(3),
116     &           n,neigh,kneigh)
117!
118            isol=0
119!
120            itriold=0
121            itri=neigh(1)+itietri(1,i)-1
122            ntriangle=0
123            ntriangle_=100
124!
125            loop1: do
126               do l=1,3
127                  ll=4*l-3
128                  dist=straight(ll,itri)*p(1)+
129     &                 straight(ll+1,itri)*p(2)+
130     &                 straight(ll+2,itri)*p(3)+
131     &                 straight(ll+3,itri)
132                  if(dist.gt.1.d-6) then
133                     itrinew=imastop(l,itri)
134                     if(itrinew.eq.0) then
135c     write(*,*) '**border reached'
136                        exit loop1
137                     elseif(itrinew.eq.itriold) then
138c     write(*,*) '**solution in between triangles'
139                        isol=itri
140                        exit loop1
141                     else
142                        call nident(itriangle,itrinew,ntriangle,id)
143                        if(id.gt.0) then
144                           if(itriangle(id).eq.itrinew) then
145c     write(*,*) '**circular path; no solution'
146                              exit loop1
147                           endif
148                        endif
149                        ntriangle=ntriangle+1
150                        if(ntriangle.gt.ntriangle_) then
151c     write(*,*) '**too many iterations'
152                           exit loop1
153                        endif
154                        do k=ntriangle,id+2,-1
155                           itriangle(k)=itriangle(k-1)
156                        enddo
157                        itriangle(id+1)=itrinew
158                        itriold=itri
159                        itri=itrinew
160                        cycle loop1
161                     endif
162                  elseif(l.eq.3) then
163c     write(*,*) '**regular solution'
164                     isol=itri
165                     exit loop1
166                  endif
167               enddo
168            enddo loop1
169!
170            if(isol.ne.0) then
171!     Active node
172               islavact(j)=1
173            else
174               islavact(j)=0
175            endif
176         enddo
177      enddo
178!
179      return
180      end
181
182