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