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 subroutine mafilldmss(co,nk,kon,ipkon,lakon,ne, 20 & ipompc,nodempc,coefmpc,nmpc, 21 & nelemload,sideload,xload,nload,xbody,ipobody,nbody,cgr, 22 & ad,au,nactdof,jq,irow,neq,nmethod, 23 & ikmpc,ilmpc,elcon,nelcon,rhcon, 24 & nrhcon,alcon,nalcon,alzero,ielmat,ielorien,norien,orab, 25 & ntmat_,t0,t1,ithermal,vold,iperturb,sti,stx,iexpl,plicon, 26 & nplicon,plkcon,nplkcon,xstiff,npmat_,dtime, 27 & matname,mi,ncmat_, 28 & physcon,ttime,time,istep,iinc, 29 & ibody,xloadold,reltime,veold,springarea,nstate_, 30 & xstateini,xstate,thicke,integerglob,doubleglob, 31 & tieset,istartset,iendset,ialset,ntie,nasym,pslavsurf, 32 & pmastsurf,mortar,clearini,ielprop,prop,ne0,nea,neb, 33 & freq,ndamp,dacon,set,nset) 34! 35! filling the stiffness matrix in spare matrix format (sm) 36! 37 implicit none 38! 39 integer mass(2),stiffness(2),buckling,rhsi,coriolis 40! 41 character*8 lakon(*) 42 character*20 sideload(*) 43 character*80 matname(*) 44 character*81 tieset(3,*),set(*) 45! 46 integer kon(*),ipompc(*),nodempc(3,*),nelemload(2,*),jq(*), 47 & ikmpc(*),ilmpc(*),mi(*),nstate_,ne0,nasym,konl(20), 48 & nactdof(0:mi(2),*),irow(*),icolumn,ialset(*),ielprop(*),one, 49 & nelcon(2,*),nrhcon(*),nalcon(2,*),ielmat(mi(3),*),ntie, 50 & jfaces,ielorien(mi(3),*),integerglob(*),istartset(*), 51 & iendset(*),ipkon(*),intscheme,ipobody(2,*),nbody,nset, 52 & ibody(3,*),nk,ne,nmpc,nload,neq(2),nmethod,mscalmethod, 53 & ithermal(*),iperturb(*),i,j,k,l,m,idist,jj, 54 & ll,id,id1,id2,ist,ist1,ist2,index,jdof1,jdof2,idof1,idof2, 55 & mpc1,mpc2,index1,index2,node1,node2,kflag,icalccg,ndamp, 56 & ntmat_,indexe,nope,norien,iexpl,i0,ncmat_,istep,iinc,imat, 57 & nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),npmat_,mortar,kode, 58 & nea,neb,kscale,ndof,ii,igauss,islavelinv(1),irowtloc(1), 59 & jqtloc(1),mortartrafoflag 60! 61 real*8 co(3,*),coefmpc(*),xload(2,*),p1(3),smscale(1), 62 & p2(3),ad(*),au(*),bodyf(3),xloadold(2,*),reltime, 63 & t0(*),t1(*),vold(0:mi(2),*),s(60,60), 64 & ff(60),xl(3,26),voldl(0:mi(2),26), 65 & sti(6,mi(1),*),sm(60,60),stx(6,mi(1),*), 66 & elcon(0:ncmat_,ntmat_,*),rhcon(0:1,ntmat_,*),springarea(2,*), 67 & alcon(0:6,ntmat_,*),physcon(*),prop(*), 68 & xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*), 69 & alzero(*),orab(7,*),xbody(7,*),cgr(4,*), 70 & plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*), 71 & xstiff(27,mi(1),*),veold(0:mi(2),*),om,value,dtime,ttime, 72 & time,thicke(mi(3),*),doubleglob(*),clearini(3,9,*),damping, 73 & pslavsurf(3,*),pmastsurf(6,*),freq,elas(21),dacon(*), 74 & autloc(1) 75! 76 mortartrafoflag=0 77 one=1 78! 79 kflag=2 80 i0=0 81 icalccg=0 82! 83 mass(1)=0 84 stiffness(1)=1 85 buckling=0 86 rhsi=0 87 intscheme=0 88 coriolis=0 89 kscale=1 90! 91 do i=nea,neb 92! 93 if(ipkon(i).lt.0) cycle 94! 95 if(ndamp.gt.0) then 96 damping=dacon(ielmat(1,i)) 97 if(mi(3).gt.1) then 98 do j=2,mi(3) 99 if(ielmat(j,i).gt.0) then 100 if(dacon(ielmat(j,i)).ne.damping) then 101 write(*,*) '*ERROR in mafilldmss: element',i 102 write(*,*) ' contains several layers' 103 write(*,*) ' with different damping ' 104 write(*,*) ' coefficients:' 105 write(*,*) ' for layer ',one,':',damping 106 write(*,*) ' for layer ',j,':', 107 & dacon(ielmat(j,i)) 108 call exit(201) 109 endif 110 endif 111 enddo 112 endif 113 else 114 damping=0.d0 115 endif 116 117 if(lakon(i)(1:2).eq.'ED') then 118 indexe=ipkon(i) 119 read(lakon(i)(8:8),'(i1)') nope 120 nope=nope+1 121! 122 do j=1,nope 123 konl(j)=kon(indexe+j) 124 enddo 125! 126 call e_damp(co,nk,konl,lakon(i),p1,p2,om,bodyf,nbody,s, 127 & sm,ff,i,elcon,nelcon,rhcon,nrhcon,alcon,nalcon, 128 & alzero,ielmat,ielorien,norien,orab,ntmat_, 129 & t0,t1,ithermal,vold,iperturb, 130 & nelemload,sideload,xload,nload,idist,sti,stx, 131 & iexpl,plicon,nplicon,plkcon,nplkcon,xstiff,npmat_, 132 & dtime,matname,mi(1),ncmat_,ttime,time,istep,iinc, 133 & nmethod) 134 ndof=3 135 do jj=1,ndof*nope 136 do ii=1,jj 137 s(ii,jj)=s(ii,jj)*freq 138 enddo 139 enddo 140 elseif((lakon(i)(1:2).eq.'ES').and. 141 & (lakon(i)(7:7).eq.'C')) then 142 indexe=ipkon(i) 143 if(mortar.eq.0) then 144 read(lakon(i)(8:8),'(i1)') nope 145 nope=nope+1 146 konl(nope+1)=kon(indexe+nope+1) 147 elseif(mortar.eq.1) then 148 nope=kon(indexe) 149 endif 150 imat=ielmat(1,i) 151! 152! computation of the coordinates of the local nodes 153! 154 do k=1,nope 155 konl(k)=kon(indexe+k) 156 do j=1,3 157 xl(j,k)=co(j,konl(k)) 158 voldl(j,k)=vold(j,konl(k)) 159 enddo 160 enddo 161! 162! initialisation of s 163! 164 do k=1,3*nope 165 do j=1,3*nope 166 s(k,j)=0.d0 167 enddo 168 enddo 169! 170 kode=nelcon(1,imat) 171! 172! as soon as the first contact element is discovered ne0 is 173! determined and saved 174! 175c if(ne0.eq.0) ne0=i-1 176 if(mortar.eq.0) then 177 call springdamp_n2f(xl,elas,voldl,s,imat,elcon, 178 & ncmat_,ntmat_,nope,iperturb, 179 & springarea(1,konl(nope+1)),nmethod, 180 & mi,reltime,nasym) 181 elseif(mortar.eq.1) then 182 jfaces=kon(indexe+nope+2) 183 igauss=kon(indexe+nope+1) 184 call springdamp_f2f(xl,elas,voldl,s,imat,elcon, 185 & ncmat_,ntmat_,nope,lakon(i),iperturb, 186 & springarea(1,igauss), 187 & nmethod,mi,reltime,nasym,jfaces,igauss,pslavsurf, 188 & pmastsurf,clearini) 189 endif 190 ndof=3 191 do jj=1,ndof*nope 192 do ii=1,jj 193 s(ii,jj)=s(ii,jj)*freq 194 enddo 195 enddo 196 elseif(damping.ne.0.d0) then 197 indexe=ipkon(i) 198c Bernhardi start 199 if(lakon(i)(1:5).eq.'C3D8I') then 200 nope=11 201 ndof=3 202 elseif(lakon(i)(4:5).eq.'20') then 203c Bernhardi end 204 nope=20 205 ndof=3 206 elseif(lakon(i)(4:4).eq.'8') then 207 nope=8 208 ndof=3 209 elseif(lakon(i)(4:5).eq.'10') then 210 nope=10 211 ndof=3 212 elseif(lakon(i)(4:4).eq.'4') then 213 nope=4 214 ndof=3 215 elseif(lakon(i)(4:5).eq.'15') then 216 nope=15 217 ndof=3 218 elseif(lakon(i)(4:4).eq.'6') then 219 nope=6 220 ndof=3 221 elseif((lakon(i)(1:2).eq.'ES').and. 222 & (lakon(i)(7:7).ne.'F')) then 223! 224! spring and contact spring elements (NO dashpot elements 225! = ED... elements) 226! 227 nope=ichar(lakon(i)(8:8))-47 228 ndof=3 229 elseif(lakon(i)(1:4).eq.'MASS') then 230 nope=1 231 ndof=3 232 elseif(lakon(i)(1:1).eq.'U') then 233 ndof=ichar(lakon(i)(7:7)) 234 nope=ichar(lakon(i)(8:8)) 235 else 236 cycle 237 endif 238! 239 om=0.d0 240! 241 if((nbody.gt.0).and.(lakon(i)(1:1).ne.'E')) then 242! 243! assigning centrifugal forces 244! 245 bodyf(1)=0.d0 246 bodyf(2)=0.d0 247 bodyf(3)=0.d0 248! 249 index=i 250 do 251 j=ipobody(1,index) 252 if(j.eq.0) exit 253 if(ibody(1,j).eq.1) then 254 om=xbody(1,j) 255 p1(1)=xbody(2,j) 256 p1(2)=xbody(3,j) 257 p1(3)=xbody(4,j) 258 p2(1)=xbody(5,j) 259 p2(2)=xbody(6,j) 260 p2(3)=xbody(7,j) 261! 262! assigning gravity forces 263! 264 elseif(ibody(1,j).eq.2) then 265 bodyf(1)=bodyf(1)+xbody(1,j)*xbody(2,j) 266 bodyf(2)=bodyf(2)+xbody(1,j)*xbody(3,j) 267 bodyf(3)=bodyf(3)+xbody(1,j)*xbody(4,j) 268! 269! assigning newton gravity forces 270! 271 elseif(ibody(1,j).eq.3) then 272 call newton(icalccg,ne,ipkon,lakon,kon,t0,co,rhcon, 273 & nrhcon,ntmat_,physcon,i,cgr,bodyf,ielmat, 274 & ithermal,vold,mi) 275 endif 276 index=ipobody(2,index) 277 if(index.eq.0) exit 278 enddo 279 endif 280! 281 if(lakon(i)(1:1).ne.'U') then 282 call e_c3d(co,kon,lakon(i),p1,p2,om,bodyf,nbody,s,sm,ff, 283 & i,nmethod,elcon,nelcon,rhcon,nrhcon,alcon,nalcon, 284 & alzero,ielmat,ielorien,norien,orab,ntmat_, 285 & t0,t1,ithermal,vold,iperturb,nelemload,sideload,xload, 286 & nload,idist,sti,stx,iexpl,plicon, 287 & nplicon,plkcon,nplkcon,xstiff,npmat_, 288 & dtime,matname,mi(1),ncmat_,mass(1),stiffness,buckling, 289 & rhsi,intscheme,ttime,time,istep,iinc,coriolis,xloadold, 290 & reltime,ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,veold, 291 & springarea,nstate_,xstateini,xstate,ne0,ipkon,thicke, 292 & integerglob,doubleglob,tieset,istartset, 293 & iendset,ialset,ntie,nasym,pslavsurf,pmastsurf,mortar, 294 & clearini,ielprop,prop,kscale,smscale(1),mscalmethod, 295 & set,nset,islavelinv, 296 & autloc,irowtloc,jqtloc,mortartrafoflag) 297 else 298 call e_c3d_u(co,kon,lakon(i),p1,p2,om,bodyf,nbody,s,sm, 299 & ff,i,nmethod,elcon,nelcon,rhcon,nrhcon,alcon,nalcon, 300 & alzero,ielmat,ielorien,norien,orab,ntmat_, 301 & t0,t1,ithermal,vold,iperturb,nelemload,sideload,xload, 302 & nload,idist,sti,stx,iexpl,plicon, 303 & nplicon,plkcon,nplkcon,xstiff,npmat_, 304 & dtime,matname,mi(1),ncmat_,mass(1),stiffness,buckling, 305 & rhsi,intscheme,ttime,time,istep,iinc,coriolis,xloadold, 306 & reltime,ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,veold, 307 & ne0,ipkon,thicke, 308 & integerglob,doubleglob,tieset,istartset, 309 & iendset,ialset,ntie,nasym, 310 & ielprop,prop) 311 endif 312 do jj=1,ndof*nope 313 do ii=1,jj 314 s(ii,jj)=s(ii,jj)*damping 315 enddo 316 enddo 317 else 318 cycle 319 endif 320! 321 if(nasym.eq.0) then 322 do jj=1,ndof*nope 323! 324 j=(jj-1)/ndof+1 325 k=jj-ndof*(j-1) 326! 327 node1=kon(indexe+j) 328 jdof1=nactdof(k,node1) 329! 330 do ll=jj,ndof*nope 331! 332 l=(ll-1)/ndof+1 333 m=ll-ndof*(l-1) 334! 335 node2=kon(indexe+l) 336 jdof2=nactdof(m,node2) 337! 338! check whether one of the DOF belongs to a SPC or MPC 339! 340 if((jdof1.gt.0).and.(jdof2.gt.0)) then 341 call add_sm_st(au,ad,jq,irow,jdof1,jdof2, 342 & s(jj,ll),jj,ll) 343 elseif((jdof1.gt.0).or.(jdof2.gt.0)) then 344! 345! idof1: genuine DOF 346! idof2: nominal DOF of the SPC/MPC 347! 348 if(jdof1.le.0) then 349 idof1=jdof2 350 idof2=jdof1 351 else 352 idof1=jdof1 353 idof2=jdof2 354 endif 355 if(nmpc.gt.0) then 356 if(idof2.ne.2*(idof2/2)) then 357! 358! regular DOF / MPC 359! 360 id=(-idof2+1)/2 361 ist=ipompc(id) 362 index=nodempc(3,ist) 363 if(index.eq.0) cycle 364 do 365 idof2=nactdof(nodempc(2,index), 366 & nodempc(1,index)) 367 value=-coefmpc(index)*s(jj,ll)/ 368 & coefmpc(ist) 369 if(idof1.eq.idof2) value=2.d0*value 370 if(idof2.gt.0) then 371 call add_sm_st(au,ad,jq,irow,idof1, 372 & idof2,value,i0,i0) 373 endif 374 index=nodempc(3,index) 375 if(index.eq.0) exit 376 enddo 377 cycle 378 endif 379 endif 380! 381! regular DOF / SPC 382! 383 if(rhsi.eq.1) then 384 elseif(nmethod.eq.2) then 385 value=s(jj,ll) 386 icolumn=neq(2)-idof2/2 387 call add_bo_st(au,jq,irow,idof1,icolumn,value) 388 endif 389 else 390 idof1=jdof1 391 idof2=jdof2 392 mpc1=0 393 mpc2=0 394 if(nmpc.gt.0) then 395 if(idof1.ne.2*(idof1/2)) mpc1=1 396 if(idof2.ne.2*(idof2/2)) mpc2=1 397 endif 398 if((mpc1.eq.1).and.(mpc2.eq.1)) then 399 id1=(-idof1+1)/2 400 id2=(-idof2+1)/2 401 if(id1.eq.id2) then 402! 403! MPC id1 / MPC id1 404! 405 ist=ipompc(id1) 406 index1=nodempc(3,ist) 407 if(index1.eq.0) cycle 408 do 409 idof1=nactdof(nodempc(2,index1), 410 & nodempc(1,index1)) 411 index2=index1 412 do 413 idof2=nactdof(nodempc(2,index2), 414 & nodempc(1,index2)) 415 value=coefmpc(index1)*coefmpc(index2)* 416 & s(jj,ll)/coefmpc(ist)/coefmpc(ist) 417 if((idof1.gt.0).and.(idof2.gt.0)) then 418 call add_sm_st(au,ad,jq,irow, 419 & idof1,idof2,value,i0,i0) 420 endif 421! 422 index2=nodempc(3,index2) 423 if(index2.eq.0) exit 424 enddo 425 index1=nodempc(3,index1) 426 if(index1.eq.0) exit 427 enddo 428 else 429! 430! MPC id1 / MPC id2 431! 432 ist1=ipompc(id1) 433 index1=nodempc(3,ist1) 434 if(index1.eq.0) cycle 435 do 436 idof1=nactdof(nodempc(2,index1), 437 & nodempc(1,index1)) 438 ist2=ipompc(id2) 439 index2=nodempc(3,ist2) 440 if(index2.eq.0) then 441 index1=nodempc(3,index1) 442 if(index1.eq.0) then 443 exit 444 else 445 cycle 446 endif 447 endif 448 do 449 idof2=nactdof(nodempc(2,index2), 450 & nodempc(1,index2)) 451 value=coefmpc(index1)*coefmpc(index2)* 452 & s(jj,ll)/coefmpc(ist1)/coefmpc(ist2) 453 if(idof1.eq.idof2) value=2.d0*value 454 if((idof1.gt.0).and.(idof2.gt.0)) then 455 call add_sm_st(au,ad,jq,irow, 456 & idof1,idof2,value,i0,i0) 457 endif 458! 459 index2=nodempc(3,index2) 460 if(index2.eq.0) exit 461 enddo 462 index1=nodempc(3,index1) 463 if(index1.eq.0) exit 464 enddo 465 endif 466 endif 467 endif 468 enddo 469! 470! 471 enddo 472 endif 473 enddo 474! 475 return 476 end 477 478