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 mafilldm(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun, 20 & xboun,nboun, 21 & ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc, 22 & nforc,nelemload,sideload,xload,nload,xbody,ipobody,nbody,cgr, 23 & ad,au,nactdof,icol,jq,irow,neq,nzl,nmethod, 24 & ikmpc,ilmpc,ikboun,ilboun,elcon,nelcon,rhcon, 25 & nrhcon,alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_, 26 & t0,t1,ithermal,prestr, 27 & iprestr,vold,iperturb,sti,nzs,stx,adb,aub,iexpl,plicon, 28 & nplicon,plkcon,nplkcon,xstiff,npmat_,dtime, 29 & matname,mi,ncmat_,ttime,time,istep,iinc,ibody, 30 & clearini,mortar,springarea,pslavsurf,pmastsurf,reltime,nasym) 31! 32! filling the damping matrix in spare matrix format (sm) 33! 34 implicit none 35! 36 character*8 lakon(*) 37 character*20 sideload(*) 38 character*80 matname(*) 39! 40 integer kon(*),nodeboun(*),ndirboun(*),ipompc(*),nodempc(3,*), 41 & nodeforc(2,*),ndirforc(*),nelemload(2,*),icol(*),jq(*),ikmpc(*), 42 & ilmpc(*),ikboun(*),ilboun(*),mi(*),nactdof(0:mi(2),*),konl(20), 43 & irow(*),nelcon(2,*),nrhcon(*),nalcon(2,*),ielmat(mi(3),*), 44 & ielorien(mi(3),*),ipkon(*),ipobody(2,*),nbody,ibody(3,*), 45 & nk,ne,nboun,nmpc,nforc,nload,neq(2),nzl,nmethod,icolumn, 46 & ithermal(*),iprestr,iperturb(*),nzs(3),i,j,k,l,m,idist,jj, 47 & ll,id,id1,id2,ist,ist1,ist2,index,jdof1,jdof2,idof1,idof2, 48 & mpc1,mpc2,index1,index2,node1,node2,ne0,igauss,imat, 49 & ntmat_,indexe,nope,norien,iexpl,i0,ncmat_,istep,iinc, 50 & nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),npmat_,jfaces, 51 & kode,mortar,nasym 52! 53 real*8 co(3,*),xboun(*),coefmpc(*),xforc(*),xload(2,*),p1(3), 54 & p2(3),ad(*),au(*),bodyf(3),voldl(0:mi(2),26),clearini(3,9,*), 55 & t0(*),t1(*),prestr(6,mi(1),*),vold(0:mi(2),*),s(60,60), 56 & ff(60), 57 & sti(6,mi(1),*),sm(60,60),stx(6,mi(1),*),adb(*),aub(*), 58 & elcon(0:ncmat_,ntmat_,*),rhcon(0:1,ntmat_,*),elas(21), 59 & alcon(0:6,ntmat_,*),alzero(*),orab(7,*),xbody(7,*),cgr(4,*), 60 & plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*), 61 & xstiff(27,mi(1),*),om,value,dtime,ttime,time,elconloc(21), 62 & xl(3,26),springarea(2,*),pslavsurf(3,*),reltime,pmastsurf(6,*) 63! 64 i0=0 65! 66! determining nzl 67! 68 nzl=0 69 do i=neq(2),1,-1 70 if(icol(i).gt.0) then 71 nzl=i 72 exit 73 endif 74 enddo 75c! 76c! initializing the matrices 77c! 78c do i=1,neq(2) 79c ad(i)=0.d0 80c enddo 81c do i=1,nzs(2) 82c au(i)=0.d0 83c enddo 84! 85 if((ithermal(1).le.1).or.(ithermal(1).eq.3)) then 86! 87! mechanical analysis: loop over all elements 88! 89 ne0=0 90 do i=1,ne 91! 92 if(ipkon(i).lt.0) cycle 93 if(lakon(i)(1:2).eq.'ED') then 94 indexe=ipkon(i) 95 read(lakon(i)(8:8),'(i1)') nope 96 nope=nope+1 97! 98 do j=1,nope 99 konl(j)=kon(indexe+j) 100 enddo 101! 102 call e_damp(co,nk,konl,lakon(i),p1,p2,om,bodyf,nbody,s, 103 & sm,ff,i,elcon,nelcon,rhcon,nrhcon,alcon,nalcon, 104 & alzero,ielmat,ielorien,norien,orab,ntmat_, 105 & t0,t1,ithermal,vold,iperturb, 106 & nelemload,sideload,xload,nload,idist,sti,stx, 107 & iexpl,plicon,nplicon,plkcon,nplkcon,xstiff,npmat_, 108 & dtime,matname,mi(1),ncmat_,ttime,time,istep,iinc, 109 & nmethod) 110 elseif((lakon(i)(1:2).eq.'ES').and. 111 & (lakon(i)(7:7).eq.'C')) then 112 indexe=ipkon(i) 113 if(mortar.eq.0) then 114 read(lakon(i)(8:8),'(i1)') nope 115 nope=nope+1 116 konl(nope+1)=kon(indexe+nope+1) 117 elseif(mortar.eq.1) then 118 nope=kon(indexe) 119 endif 120 imat=ielmat(1,i) 121! 122! computation of the coordinates of the local nodes 123! 124 do k=1,nope 125 konl(k)=kon(indexe+k) 126 do j=1,3 127 xl(j,k)=co(j,konl(k)) 128 voldl(j,k)=vold(j,konl(k)) 129 enddo 130 enddo 131! 132! initialisation of s 133! 134 do k=1,3*nope 135 do j=1,3*nope 136 s(k,j)=0.d0 137 enddo 138 enddo 139! 140 kode=nelcon(1,imat) 141! 142! as soon as the first contact element is discovered ne0 is 143! determined and saved 144! 145 if(ne0.eq.0) ne0=i-1 146 if(mortar.eq.0) then 147 call springdamp_n2f(xl,elas,voldl,s,imat,elcon, 148 & ncmat_,ntmat_,nope,iperturb, 149 & springarea(1,konl(nope+1)),nmethod, 150 & mi,reltime,nasym) 151 elseif(mortar.eq.1) then 152 jfaces=kon(indexe+nope+2) 153 igauss=kon(indexe+nope+1) 154 call springdamp_f2f(xl,elas,voldl,s,imat,elcon, 155 & ncmat_,ntmat_,nope,lakon(i),iperturb, 156 & springarea(1,igauss), 157 & nmethod,mi,reltime,nasym,jfaces,igauss,pslavsurf, 158 & pmastsurf,clearini) 159 endif 160 else 161 cycle 162 endif 163! 164 if(nasym.eq.0) then 165! 166 do jj=1,3*nope 167! 168 j=(jj-1)/3+1 169 k=jj-3*(j-1) 170! 171 node1=kon(indexe+j) 172 jdof1=nactdof(k,node1) 173! 174 do ll=jj,3*nope 175! 176 l=(ll-1)/3+1 177 m=ll-3*(l-1) 178! 179 node2=kon(indexe+l) 180 jdof2=nactdof(m,node2) 181! 182! check whether one of the DOF belongs to a SPC or MPC 183! 184 if((jdof1.gt.0).and.(jdof2.gt.0)) then 185c write(*,*) 'mafilldm1',jdof1,jdof2,jj,ll,s(jj,ll) 186 call add_sm_st(au,ad,jq,irow,jdof1,jdof2, 187 & s(jj,ll),jj,ll) 188 elseif((jdof1.gt.0).or.(jdof2.gt.0)) then 189! 190! idof1: genuine DOF 191! idof2: nominal DOF of the SPC/MPC 192! 193 if(jdof1.le.0) then 194 idof1=jdof2 195c idof2=(node1-1)*8+k 196 idof2=jdof1 197 else 198 idof1=jdof1 199c idof2=(node2-1)*8+m 200 idof2=jdof2 201 endif 202 if(nmpc.gt.0) then 203c call nident(ikmpc,idof2,nmpc,id) 204c if((id.gt.0).and.(ikmpc(id).eq.idof2)) then 205 if(idof2.ne.2*(idof2/2)) then 206! 207! regular DOF / MPC 208! 209c id=ilmpc(id) 210 id=(-idof2+1)/2 211 ist=ipompc(id) 212 index=nodempc(3,ist) 213 if(index.eq.0) cycle 214 do 215 idof2=nactdof(nodempc(2,index), 216 & nodempc(1,index)) 217 value=-coefmpc(index)*s(jj,ll)/ 218 & coefmpc(ist) 219 if(idof1.eq.idof2) value=2.d0*value 220 if(idof2.gt.0) then 221c write(*,*) 'mafilldm2',idof1,idof2,i0,i0,value 222 call add_sm_st(au,ad,jq,irow,idof1, 223 & idof2,value,i0,i0) 224 endif 225 index=nodempc(3,index) 226 if(index.eq.0) exit 227 enddo 228 cycle 229 endif 230 endif 231 else 232c idof1=(node1-1)*8+k 233c idof2=(node2-1)*8+m 234 idof1=jdof1 235 idof2=jdof2 236 mpc1=0 237 mpc2=0 238 if(nmpc.gt.0) then 239c call nident(ikmpc,idof1,nmpc,id1) 240c if((id1.gt.0).and.(ikmpc(id1).eq.idof1)) mpc1=1 241c call nident(ikmpc,idof2,nmpc,id2) 242c if((id2.gt.0).and.(ikmpc(id2).eq.idof2)) mpc2=1 243 if(idof1.ne.2*(idof1/2)) mpc1=1 244 if(idof2.ne.2*(idof2/2)) mpc2=1 245 endif 246 if((mpc1.eq.1).and.(mpc2.eq.1)) then 247c id1=ilmpc(id1) 248c id2=ilmpc(id2) 249 id1=(-idof1+1)/2 250 id2=(-idof2+1)/2 251 if(id1.eq.id2) then 252! 253! MPC id1 / MPC id1 254! 255 ist=ipompc(id1) 256 index1=nodempc(3,ist) 257 if(index1.eq.0) cycle 258 do 259 idof1=nactdof(nodempc(2,index1), 260 & nodempc(1,index1)) 261 index2=index1 262 do 263 idof2=nactdof(nodempc(2,index2), 264 & nodempc(1,index2)) 265 value=coefmpc(index1)*coefmpc(index2)* 266 & s(jj,ll)/coefmpc(ist)/coefmpc(ist) 267 if((idof1.gt.0).and.(idof2.gt.0)) then 268c write(*,*) 'mafilldm3',idof1,idof2,i0,i0,value 269 call add_sm_st(au,ad,jq,irow, 270 & idof1,idof2,value,i0,i0) 271 endif 272! 273 index2=nodempc(3,index2) 274 if(index2.eq.0) exit 275 enddo 276 index1=nodempc(3,index1) 277 if(index1.eq.0) exit 278 enddo 279 else 280! 281! MPC id1 / MPC id2 282! 283 ist1=ipompc(id1) 284 index1=nodempc(3,ist1) 285 if(index1.eq.0) cycle 286 do 287 idof1=nactdof(nodempc(2,index1), 288 & nodempc(1,index1)) 289 ist2=ipompc(id2) 290 index2=nodempc(3,ist2) 291 if(index2.eq.0) then 292 index1=nodempc(3,index1) 293 if(index1.eq.0) then 294 exit 295 else 296 cycle 297 endif 298 endif 299 do 300 idof2=nactdof(nodempc(2,index2), 301 & nodempc(1,index2)) 302 value=coefmpc(index1)*coefmpc(index2)* 303 & s(jj,ll)/coefmpc(ist1)/ 304 & coefmpc(ist2) 305 if(idof1.eq.idof2) value=2.d0*value 306 if((idof1.gt.0).and.(idof2.gt.0)) then 307c write(*,*) 'mafilldm4',idof1,idof2,i0,i0,value 308 call add_sm_st(au,ad,jq,irow, 309 & idof1,idof2,value,i0,i0) 310 endif 311! 312 index2=nodempc(3,index2) 313 if(index2.eq.0) exit 314 enddo 315 index1=nodempc(3,index1) 316 if(index1.eq.0) exit 317 enddo 318 endif 319 endif 320 endif 321 enddo 322! 323 enddo 324! 325! asymmetrical calculations (due to friction or tangential damping...) 326! 327 else 328! 329 do jj=1,3*nope 330! 331 j=(jj-1)/3+1 332 k=jj-3*(j-1) 333! 334 node1=kon(indexe+j) 335 jdof1=nactdof(k,node1) 336! 337 do ll=1,3*nope 338! 339 l=(ll-1)/3+1 340 m=ll-3*(l-1) 341! 342 node2=kon(indexe+l) 343 jdof2=nactdof(m,node2) 344! 345! check whether one of the DOF belongs to a SPC or MPC 346! 347 if((jdof1.gt.0).and.(jdof2.gt.0)) then 348 call add_sm_st_as(au,ad,jq,irow,jdof1,jdof2, 349 & s(jj,ll),jj,ll,nzs) 350 elseif((jdof1.gt.0).or.(jdof2.gt.0)) then 351! 352! idof1: genuine DOF 353! idof2: nominal DOF of the SPC/MPC 354! 355 if(jdof1.le.0) then 356c idof1=(node1-1)*8+k 357 idof1=jdof1 358 idof2=jdof2 359 if(nmpc.gt.0) then 360c call nident(ikmpc,idof1,nmpc,id) 361c if((id.gt.0).and.(ikmpc(id).eq.idof1)) then 362 if(idof1.ne.2*(idof1/2)) then 363! 364! regular DOF / MPC 365! 366c id=ilmpc(id) 367 id=(-idof1+1)/2 368 ist=ipompc(id) 369 index=nodempc(3,ist) 370 if(index.eq.0) cycle 371 do 372 idof1=nactdof(nodempc(2,index), 373 & nodempc(1,index)) 374 value=-coefmpc(index)*s(jj,ll)/coefmpc(ist) 375 if(idof1.gt.0) then 376c write(*,*) 'mafilldm6',idof1,idof2,i0,i0,value 377 call add_sm_st_as(au,ad,jq,irow,idof1, 378 & idof2,value,i0,i0,nzs) 379 endif 380 index=nodempc(3,index) 381 if(index.eq.0) exit 382 enddo 383 cycle 384 endif 385 endif 386 else 387 idof1=jdof1 388c idof2=(node2-1)*8+m 389 idof2=jdof2 390 if(nmpc.gt.0) then 391c call nident(ikmpc,idof2,nmpc,id) 392c if((id.gt.0).and.(ikmpc(id).eq.idof2)) then 393 if(idof2.ne.2*(idof2/2)) then 394! 395! regular DOF / MPC 396! 397c id=ilmpc(id) 398 id=(-idof2+1)/2 399 ist=ipompc(id) 400 index=nodempc(3,ist) 401 if(index.eq.0) cycle 402 do 403 idof2=nactdof(nodempc(2,index), 404 & nodempc(1,index)) 405 value=-coefmpc(index)*s(jj,ll)/coefmpc(ist) 406 if(idof2.gt.0) then 407c write(*,*) 'mafilldm7',idof1,idof2,i0,i0,value 408 call add_sm_st_as(au,ad,jq,irow,idof1, 409 & idof2,value,i0,i0,nzs) 410 endif 411 index=nodempc(3,index) 412 if(index.eq.0) exit 413 enddo 414 cycle 415 endif 416 endif 417 endif 418 else 419c idof1=(node1-1)*8+k 420c idof2=(node2-1)*8+m 421 idof1=jdof1 422 idof2=jdof2 423 mpc1=0 424 mpc2=0 425 if(nmpc.gt.0) then 426c call nident(ikmpc,idof1,nmpc,id1) 427c if((id1.gt.0).and.(ikmpc(id1).eq.idof1)) mpc1=1 428c call nident(ikmpc,idof2,nmpc,id2) 429c if((id2.gt.0).and.(ikmpc(id2).eq.idof2)) mpc2=1 430 if(idof1.ne.2*(idof1/2)) mpc1=1 431 if(idof2.ne.2*(idof2/2)) mpc2=1 432 endif 433 if((mpc1.eq.1).and.(mpc2.eq.1)) then 434c id1=ilmpc(id1) 435c id2=ilmpc(id2) 436 id1=(-idof1+1)/2 437 id2=(-idof2+1)/2 438 if(id1.eq.id2) then 439! 440! MPC id1 / MPC id1 441! 442 ist1=ipompc(id1) 443 index1=nodempc(3,ist1) 444 if(index1.eq.0) cycle 445 do 446 idof1=nactdof(nodempc(2,index1), 447 & nodempc(1,index1)) 448c index2=index1 449 ist2=ipompc(id2) 450 index2=nodempc(3,ist2) 451 if(index2.eq.0) then 452 index1=nodempc(3,index1) 453 if(index1.eq.0) then 454 exit 455 else 456 cycle 457 endif 458 endif 459 do 460 idof2=nactdof(nodempc(2,index2), 461 & nodempc(1,index2)) 462c value=coefmpc(index1)*coefmpc(index2)* 463c & s(jj,ll)/coefmpc(ist)/coefmpc(ist) 464 value=coefmpc(index1)*coefmpc(index2)* 465 & s(jj,ll)/coefmpc(ist1)/coefmpc(ist2) 466 if((idof1.gt.0).and.(idof2.gt.0)) then 467c write(*,*) 'mafilldm8',idof1,idof2,i0,i0,value 468 call add_sm_st_as(au,ad,jq,irow, 469 & idof1,idof2,value,i0,i0,nzs) 470 endif 471! 472 index2=nodempc(3,index2) 473 if(index2.eq.0) exit 474 enddo 475 index1=nodempc(3,index1) 476 if(index1.eq.0) exit 477 enddo 478 else 479! 480! MPC id1 / MPC id2 481! 482 ist1=ipompc(id1) 483 index1=nodempc(3,ist1) 484 if(index1.eq.0) cycle 485 do 486 idof1=nactdof(nodempc(2,index1), 487 & nodempc(1,index1)) 488 ist2=ipompc(id2) 489 index2=nodempc(3,ist2) 490 if(index2.eq.0) then 491 index1=nodempc(3,index1) 492 if(index1.eq.0) then 493 exit 494 else 495 cycle 496 endif 497 endif 498 do 499 idof2=nactdof(nodempc(2,index2), 500 & nodempc(1,index2)) 501 value=coefmpc(index1)*coefmpc(index2)* 502 & s(jj,ll)/coefmpc(ist1)/coefmpc(ist2) 503 if((idof1.gt.0).and.(idof2.gt.0)) then 504c write(*,*) 'mafilldm9',idof1,idof2,i0,i0,value 505 call add_sm_st_as(au,ad,jq,irow, 506 & idof1,idof2,value,i0,i0,nzs) 507 endif 508! 509 index2=nodempc(3,index2) 510 if(index2.eq.0) exit 511 enddo 512 index1=nodempc(3,index1) 513 if(index1.eq.0) exit 514 enddo 515 endif 516 endif 517 endif 518 enddo 519 enddo 520 endif 521 enddo 522! 523 endif 524! 525c do i=1,neq(2) 526c write(*,*) i,ad(i) 527c enddo 528c do i=1,nzs(2) 529c write(*,*) i,au(i) 530c enddo 531 return 532 end 533 534