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 to calculate all coupling matrix combinations and gap contributions of p,q for current slave face l and store it into contri. 21! calculate 22! \f$ \tilde{B}_d[p,q]=-<\tilde{\Phi}_q,\Psi_p> Id_d\f$ and \f$ \tilde{D}_d \f$. 23! see Phd-thesis Sitzmann Chapter 4 24! 25! Author:Saskia Sitzmann 26! 27! [in] ict current tie 28! [in] l current slave face 29! [in] gapmints (i) gap between slave surface and master surface in integration point i 30! [in] islavsurf islavsurf(1,i) slaveface i islavsurf(2,i) pointer into imastsurf and pmastsurf 31! [in] imastsurf index of masterface corresponding to integration point i 32! [in] pmastsurf field storing position xil, etal and weight for integration point on master side 33! [out] contr field containing B_d contributions for current face 34! [out] iscontr (i) slave node of contribution(i) 35! [out] imcontr (i) master node of contribution(i) 36! [out] dcontr field containing D_d contributions for current face 37! [out] idcontr1 (i) slave node of contribution(i) 38! [out] idcontr2 (i) master node of contribution(i) 39! [out] gcontr field containing gap contributions for current face 40! [out] igcontr (i) nodesf of contribution(i) 41! [in] pslavsurf field storing position xil, etal and weight for integration point on slave side 42! [in] pslavdual (:,i)coefficients \f$ \alpha_{ij}\f$, \f$ 1,j=1,..8\f$ for dual shape functions for face i 43! [in] nslavnode (i)pointer into field isalvnode for contact tie i 44! [in] islavnode field storing the nodes of the slave surface 45! [in] nmastnode (i)pointer into field imastnode for contact tie i 46! [in] imastnode field storing the nodes of the master surfaces 47! [out] icounter counter variable for contr 48! [out] icounter2 counter variable for dcontr 49! [in] islavact (i) indicates, if slave node i is active (=-3 no-slave-node, =-2 no-LM-node, =-1 no-gap-node, =0 inactive node, =1 sticky node, =2 slipping/active node) 50! [in] iflagdualquad flag indicating what mortar contact is used (=1 quad-lin, =2 quad-quad, =3 PG quad-lin, =4 PG quad-quad) 51! 52 subroutine createbd(ict,l,ipkon,kon,lakon,co, vold, gapmints, 53 & islavsurf,imastsurf,pmastsurf,contr,iscontr,imcontr, 54 & dcontr,idcontr1,idcontr2,gcontr,igcontr,mi, 55 & pslavsurf,pslavdual,nslavnode,islavnode,nmastnode,imastnode, 56 & icounter,icounter2,islavact,iflagdualquad) 57! 58 implicit none 59! 60 logical debug 61! 62 character*8 lakon(*) 63! 64 integer ipkon(*),kon(*),konl(20),iflag,m,l,j,jj, 65 & indexe,islavsurf(2,*), 66 & imastsurf(*),ifaces,nelemens,jfaces,ifacem, 67 & mint2d,indexf,nopes1,nopes2,nodem,nodesf, 68 & locs,locm,mi(*),ns,mint2dloc1,mint2dloc2, 69 & ifs,ifm,nope1,nope2, iscontr(*),imcontr(*),getlocno, 70 & jfacem,nelemenm,icounter,idummy,ifac,idcontr1(*),idcontr2(*), 71 & igcontr(*),icounter2,nslavnode(*),islavnode(*),ict,id, 72 & nmastnode(*),imastnode(*),islavact(*),iflagdualquad 73! 74 real*8 pmastsurf(2,*),co(3,*),gapmints(*), 75 & vold(0:mi(2),*),weight,dx,help, 76 & ets,xis,etm,xim,xl2s(3,8),xsj2s(3),xsj2s2(3), 77 & shp2s(7,8),xs2s(3,7),xl2m(3,8),xsj2m(3),shp2m(7,8),xs2m(3,7), 78 & contribution,pslavsurf(3,*),pslavdual(64,*),contr(*), 79 & dcontr(*),dcontribution,gcontr(*),gcontribution, 80 & shp2s2(7,8),xs2s2(3,7) 81! 82 debug=.false. 83 contribution = 0.d0 84 dcontribution = 0.d0 85 gcontribution = 0.d0 86 icounter=0 87 icounter2=0 88 ifaces=islavsurf(1,l) 89 nelemens = int(ifaces/10) 90 jfaces = ifaces - nelemens*10 91 indexe = ipkon(nelemens) 92 ict=ict+1 93 if(debug)write(*,*) 'createbd:l=',l,'tie',ict 94 call getnumberofnodes(nelemens,jfaces,lakon,nope1,nopes1,idummy) 95 mint2d=islavsurf(2,l+1)-islavsurf(2,l) 96 if(mint2d.eq.0) return 97 indexf=islavsurf(2,l) 98 if(debug)write(*,*) 'createbd:mint2d',mint2d 99! 100! loop over all nodesf of current slave face 101! 102 do j=1,nope1 103 konl(j)=kon(ipkon(nelemens)+j) 104 enddo 105 do m=1,nopes1 106 ifac=getlocno(m,jfaces,nope1) 107 do j=1,3 108 xl2s(j,m)=co(j,konl(ifac))+ 109 & vold(j,konl(ifac)) 110 enddo 111 enddo 112! 113 do j=1,nopes1 114 do jj=1,nopes1 115 dcontr(icounter2+nopes1*(j-1)+jj)=0.0 116 enddo 117 gcontr(icounter2+j)=0.0 118 enddo 119! 120 mint2dloc1=1 121 mint2dloc2=1 122 help=0.d0 123! 124! loop over all integration points created for current slave face 125! 126 do 127 if (mint2dloc2.ge.mint2d) exit 128! 129! find current master face 130! 131 ifacem=imastsurf(indexf+mint2dloc1) 132 nelemenm= int (ifacem/10); 133 jfacem=ifacem-nelemenm*10 134 call getnumberofnodes(nelemenm,jfacem,lakon, 135 & nope2,nopes2,idummy) 136! 137! find number of integration points belonging to master face 138! 139 do 140 if(ifacem.eq.imastsurf(indexf+mint2dloc2+1))then 141 mint2dloc2=mint2dloc2+1 142 if(mint2dloc2.eq.mint2d) exit 143 else 144 exit 145 endif 146 enddo 147 if(debug)write(*,*) 'createbd:MF,loc1, loc2',ifacem, 148 & mint2dloc1,mint2dloc2 149 help=0.0 150 do m=mint2dloc1,mint2dloc2 151 xis=pslavsurf(1,indexf+m) 152 ets=pslavsurf(2,indexf+m) 153 weight=pslavsurf(3,indexf+m) 154 ns=l 155 iflag = 2 156 if(nopes1.eq.8) then 157 if(iflagdualquad.eq.2 .or. iflagdualquad.eq.4)then 158 call dualshape8qtilde(xis,ets,xl2s,xsj2s,xs2s,shp2s, 159 & ns,pslavdual,iflag) 160 else 161 call dualshape8qtilde_lin(xis,ets,xl2s,xsj2s,xs2s, 162 & shp2s,ns,pslavdual,iflag) 163 endif 164 elseif(nopes1.eq.4) then 165 call dualshape4q(xis,ets,xl2s,xsj2s,xs2s,shp2s,ns, 166 & pslavdual,iflag) 167 elseif(nopes1.eq.6) then 168 if(iflagdualquad.eq.2 .or. iflagdualquad.eq.4)then 169 call dualshape6tritilde(xis,ets,xl2s,xsj2s,xs2s,shp2s, 170 & ns,pslavdual,iflag) 171 else 172 call dualshape6tritilde_lin(xis,ets,xl2s,xsj2s,xs2s, 173 & shp2s,ns,pslavdual,iflag) 174 endif 175 else 176 call dualshape3tri(xis,ets,xl2s,xsj2s,xs2s,shp2s,ns, 177 & pslavdual,iflag) 178 endif 179 if(nopes1.eq.8) then 180 if(iflagdualquad.eq.2 .or. iflagdualquad.eq.4)then 181 call shape8qtilde(xis,ets,xl2s,xsj2s2,xs2s2, 182 & shp2s2,iflag) 183 else 184 call shape8qtilde_lin(xis,ets,xl2s,xsj2s2,xs2s2, 185 & shp2s2,iflag) 186 endif 187 elseif(nopes1.eq.4) then 188 call shape4q(xis,ets,xl2s,xsj2s2,xs2s2, 189 & shp2s2,iflag) 190 elseif(nopes1.eq.6) then 191 if(iflagdualquad.eq.2 .or. iflagdualquad.eq.4)then 192 call shape6tritilde(xis,ets,xl2s,xsj2s2,xs2s2, 193 & shp2s2,iflag) 194 else 195 call shape6tritilde_lin(xis,ets,xl2s,xsj2s2,xs2s2, 196 & shp2s2,iflag) 197 endif 198 else 199 call shape3tri(xis,ets,xl2s,xsj2s2,xs2s2, 200 & shp2s2,iflag) 201 endif 202 xim = pmastsurf(1,indexf+m) 203 etm = pmastsurf(2,indexf+m) 204! 205 if(nopes2.eq.8) then 206 call shape8q(xim,etm,xl2m,xsj2m,xs2m,shp2m,iflag) 207 elseif(nopes2.eq.4) then 208 call shape4q(xim,etm,xl2m,xsj2m,xs2m,shp2m,iflag) 209 elseif(nopes2.eq.6) then 210 call shape6tri(xim,etm,xl2m,xsj2m,xs2m,shp2m,iflag) 211 else 212 call shape3tri(xim,etm,xl2m,xsj2m,xs2m,shp2m,iflag) 213 endif 214 if(m.eq.mint2dloc1)then 215 do j=1,nopes1 216 do jj=1,nopes2 217 contr(icounter+nopes2*(j-1)+jj)=0.0 218 enddo 219 enddo 220 endif 221 dx=dsqrt(xsj2s(1)**2+xsj2s(2)**2+xsj2s(3)**2) 222 do j=1,nopes1 223 ifs=getlocno(j,jfaces,nope1) 224 nodesf=kon(ipkon(nelemens)+ifs) 225 locs=j 226 gcontribution=shp2s(4,locs) 227 & *(gapmints(indexf+m)) 228 & *weight 229 & *dx 230 gcontr(icounter2+j)=gcontr(icounter2+j)+gcontribution 231 if(m.eq.1)then 232 call nident(islavnode(nslavnode(ict)+1), 233 & nodesf,(nslavnode(ict+1)-nslavnode(ict)),id) 234 if(islavnode(nslavnode(ict)+id).eq.nodesf) then 235 igcontr(icounter2+j)=nslavnode(ict)+id 236 else 237 write(*,*)'createbd: node',nodesf 238 write(*,*)'was not catalogued properly in', 239 & 'islavnode' 240 call exit(201) 241 endif 242 endif 243 do jj=1,nopes1 244 ifs=getlocno(jj,jfaces,nope1) 245 nodem=kon(ipkon(nelemens)+ifs) 246 if(m.eq.1)then 247 call nident(islavnode(nslavnode(ict)+1), 248 & nodesf,(nslavnode(ict+1)-nslavnode(ict)),id) 249 if(islavnode(nslavnode(ict)+id).eq.nodesf) then 250 idcontr1(icounter2+nopes1*(j-1)+jj)= 251 & nslavnode(ict)+id 252 else 253 write(*,*)'createbd: node',nodesf 254 write(*,*)'was not catalogued properly in', 255 & 'islavnode' 256 call exit(201) 257 endif 258! 259 call nident(islavnode(nslavnode(ict)+1), 260 & nodem,(nslavnode(ict+1)-nslavnode(ict)),id) 261 if(islavnode(nslavnode(ict)+id).eq.nodem) then 262 idcontr2(icounter2+nopes1*(j-1)+jj)= 263 & nslavnode(ict)+id 264 else 265 write(*,*)'createbd: node',nodem 266 write(*,*)'was not catalogued properly in', 267 & 'islavnode' 268 call exit(201) 269 endif 270! 271 endif 272 dcontribution=shp2s(4,locs)*shp2s2(4,jj) *weight*dx 273 dcontr(icounter2+nopes1*(j-1)+jj)= 274 & dcontr(icounter2+nopes1*(j-1)+jj)+dcontribution 275 enddo 276 do jj=1,nopes2 277 ifm=getlocno(jj,jfacem,nope2) 278 nodem=kon(ipkon(nelemenm)+ifm) 279 locm=jj 280 contribution=shp2s(4,locs)*shp2m(4,locm) 281 & *weight 282 & *dx 283 contr(icounter+nopes2*(j-1)+jj)= 284 & contr(icounter+nopes2*(j-1)+jj)+contribution 285 if(m.eq.mint2dloc1)then 286 call nident(islavnode(nslavnode(ict)+1), 287 & nodesf,(nslavnode(ict+1)-nslavnode(ict)),id) 288 if(islavnode(nslavnode(ict)+id).eq.nodesf) then 289 iscontr(icounter+nopes2*(j-1)+jj)= 290 & nslavnode(ict)+id 291 else 292 write(*,*)'createbd: node',nodesf 293 write(*,*)'was not catalogued properly in', 294 & ' islavnode' 295 call exit(201) 296 endif 297 call nident(imastnode(nmastnode(ict)+1), 298 & nodem,(nmastnode(ict+1)-nmastnode(ict)),id) 299 if(imastnode(nmastnode(ict)+id).eq.nodem) then 300 imcontr(icounter+nopes2*(j-1)+jj)= 301 & nmastnode(ict)+id 302 else 303 write(*,*)'createbd: node',nodem 304 write(*,*)'was not catalogued properly in', 305 & ' imastnode',nmastnode(ict)+id, 306 & imastnode(nmastnode(ict)+id), 307 & nmastnode(ict)+1, 308 & nmastnode(ict+1) 309 call exit(201) 310 endif 311! 312 endif 313 contribution=0.d0 314 enddo 315 enddo 316! 317 enddo 318 mint2dloc1=mint2dloc2+1 319 mint2dloc2=mint2dloc1 320 icounter=icounter+nopes1*nopes2 321 enddo 322 icounter2=icounter2+nopes1*nopes1 323! 324 debug=.false. 325 if(debug )then 326 write(*,*) 'createbd: gcontri,idcontr',l 327 do j=1, nopes1 328 write(*,*) gcontr(j),igcontr(j),islavact(igcontr(j)) 329 enddo 330 endif 331 ict=ict-1 332 return 333 end 334