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