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 mafillem(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,fext,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_,mass,stiffness,buckling,rhsi,intscheme, 30 & physcon,shcon,nshcon,cocon,ncocon,ttime,time,istep,iinc, 31 & coriolis,ibody,xloadold,reltime,veold,springarea,nstate_, 32 & xstateini,xstate,thicke,integerglob,doubleglob, 33 & tieset,istartset,iendset,ialset,ntie,nasym,iactive,h0, 34 & pslavsurf,pmastsurf,mortar,clearini,ielprop,prop, 35 & iponoel,inoel,network) 36! 37! filling the stiffness matrix in spare matrix format (sm) 38! 39! domain 1: phi-domain (air) 40! domain 2: A,V-domain (body) 41! domain 3: A-domain (air, the union of domain 2 and 3 should 42! be simple connected) 43! 44 implicit none 45! 46 logical mass(2),stiffness,buckling,rhsi,stiffonly(2),coriolis 47! 48 character*8 lakon(*) 49 character*20 sideload(*) 50 character*80 matname(*) 51 character*81 tieset(3,*) 52! 53 integer kon(*),nodeboun(*),ndirboun(*),ipompc(*),nodempc(3,*), 54 & nodeforc(2,*),ndirforc(*),nelemload(2,*),icol(*),jq(*),ikmpc(*), 55 & ilmpc(*),ikboun(*),ilboun(*),mi(*),nstate_,ne0,nasym, 56 & nactdof(0:mi(2),*),konl(26),irow(*),icolumn,ialset(*), 57 & nelcon(2,*),nrhcon(*),nalcon(2,*),ielmat(mi(3),*),ntie, 58 & ielorien(mi(3),*),integerglob(*),istartset(*),iendset(*), 59 & ipkon(*),intscheme,ncocon(2,*),nshcon(*),ipobody(2,*),nbody, 60 & ibody(3,*),nk,ne,nboun,nmpc,nforc,nload,neq(2),nzl,nmethod, 61 & ithermal(*),iprestr,iperturb(*),nzs(3),i,j,k,l,m,idist,jj, 62 & ll,id,id1,id2,ist,ist1,ist2,index,jdof1,jdof2,idof1,idof2, 63 & mpc1,mpc2,index1,index2,jdof,node1,node2,kflag,icalccg, 64 & ntmat_,indexe,nope,norien,iexpl,i0,ncmat_,istep,iinc,mortar, 65 & nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),npmat_,iactive(3), 66 & ielprop(*),iponoel(*),inoel(2,*),network 67! 68 real*8 co(3,*),xboun(*),coefmpc(*),xforc(*),xload(2,*),p1(3), 69 & p2(3),ad(*),au(*),bodyf(3),fext(*),xloadold(2,*),reltime, 70 & t0(*),t1(*),prestr(6,mi(1),*),vold(0:mi(2),*),s(100,100), 71 & sti(6,mi(1),*),sm(100,100),stx(6,mi(1),*),adb(*),aub(*), 72 & elcon(0:ncmat_,ntmat_,*),rhcon(0:1,ntmat_,*),springarea(2,*), 73 & alcon(0:6,ntmat_,*),physcon(*),cocon(0:6,ntmat_,*),ff(100), 74 & xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*), 75 & shcon(0:3,ntmat_,*),alzero(*),orab(7,*),xbody(7,*),cgr(4,*), 76 & plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*), 77 & xstiff(27,mi(1),*),veold(0:mi(2),*),om,valu2,value,dtime,ttime, 78 & time,thicke(mi(3),*),doubleglob(*),h0(3,*),st(60,60),smt(60,60), 79 & pslavsurf(3,*),pmastsurf(6,*),clearini(3,9,*),prop(*) 80! 81 kflag=2 82 i0=0 83 icalccg=0 84! 85 if(stiffness.and.(.not.mass(1))) then 86 stiffonly(1)=.true. 87 else 88 stiffonly(1)=.false. 89 endif 90 if(stiffness.and.(.not.mass(2))) then 91 stiffonly(2)=.true. 92 else 93 stiffonly(2)=.false. 94 endif 95! 96! determining nzl 97! 98 nzl=0 99 do i=neq(2),1,-1 100 if(icol(i).gt.0) then 101 nzl=i 102 exit 103 endif 104 enddo 105! 106! initializing the matrices 107! 108 do i=1,neq(2) 109 ad(i)=0.d0 110 enddo 111 do i=1,nzs(3) 112 au(i)=0.d0 113 enddo 114! 115 if(rhsi) then 116 do i=1,neq(2) 117 fext(i)=0.d0 118 enddo 119 endif 120! 121 if(mass(1)) then 122 do i=1,neq(1) 123 adb(i)=0.d0 124 enddo 125 do i=1,nzs(1) 126 aub(i)=0.d0 127 enddo 128 endif 129 if(mass(2)) then 130 do i=neq(1)+1,neq(2) 131 adb(i)=0.d0 132 enddo 133 do i=nzs(1)+1,nzs(2) 134 aub(i)=0.d0 135 enddo 136 endif 137! 138! electromagnetic force should always be taken into account 139! 140 if(rhsi) idist=1 141! 142 if((ithermal(1).le.1).or.(ithermal(1).eq.3)) then 143! 144! electromagnetic analysis: loop over all elements 145! 146 ne0=0 147 do i=1,ne 148! 149 if(ipkon(i).lt.0) cycle 150 indexe=ipkon(i) 151 if(lakon(i)(4:5).eq.'20') then 152 nope=20 153 elseif(lakon(i)(4:4).eq.'8') then 154 nope=8 155 elseif(lakon(i)(4:5).eq.'10') then 156 nope=10 157 elseif(lakon(i)(4:4).eq.'4') then 158 nope=4 159 elseif(lakon(i)(4:5).eq.'15') then 160 nope=15 161 elseif(lakon(i)(4:4).eq.'6') then 162 nope=6 163 else 164 cycle 165 endif 166! 167 do j=1,nope 168 konl(j)=kon(indexe+j) 169 enddo 170! 171 call e_c3d_em(co,konl,lakon(i),s,sm,ff,i,nmethod, 172 & ielmat,ntmat_,t1,ithermal,vold, 173 & idist,matname,mi,mass(1),rhsi, 174 & ncmat_,elcon,nelcon,h0,iactive, 175 & alcon,nalcon,istartset,iendset,ialset) 176! 177 do jj=1,5*nope 178! 179 j=(jj-1)/5+1 180 k=jj-5*(j-1) 181! 182 node1=kon(indexe+j) 183 jdof1=nactdof(k,node1) 184! 185 do ll=jj,5*nope 186! 187 l=(ll-1)/5+1 188 m=ll-5*(l-1) 189! 190 node2=kon(indexe+l) 191 jdof2=nactdof(m,node2) 192! 193! check whether one of the DOF belongs to a SPC or MPC 194! 195 if((jdof1.gt.0).and.(jdof2.gt.0)) then 196 if(stiffonly(1)) then 197 call add_sm_st(au,ad,jq,irow,jdof1,jdof2, 198 & s(jj,ll),jj,ll) 199 else 200 call add_sm_ei(au,ad,aub,adb,jq,irow,jdof1,jdof2, 201 & s(jj,ll),sm(jj,ll),jj,ll) 202 endif 203 elseif((jdof1.gt.0).or.(jdof2.gt.0)) then 204! 205! idof1: genuine DOF 206! idof2: nominal DOF of the SPC/MPC 207! 208 if(jdof1.le.0) then 209 idof1=jdof2 210c idof2=(node1-1)*8+k 211 idof2=jdof1 212 else 213 idof1=jdof1 214c idof2=(node2-1)*8+m 215 idof2=jdof2 216 endif 217 if(nmpc.gt.0) then 218c call nident(ikmpc,idof2,nmpc,id) 219c if((id.gt.0).and.(ikmpc(id).eq.idof2)) then 220 if(idof2.ne.2*(idof2/2)) then 221! 222! regular DOF / MPC 223! 224c id=ilmpc(id) 225 id=(-idof2+1)/2 226 ist=ipompc(id) 227 index=nodempc(3,ist) 228 if(index.eq.0) cycle 229 do 230 idof2=nactdof(nodempc(2,index),nodempc(1,index)) 231 value=-coefmpc(index)*s(jj,ll)/coefmpc(ist) 232 if(idof1.eq.idof2) value=2.d0*value 233 if(idof2.gt.0) then 234 if(stiffonly(1)) then 235 call add_sm_st(au,ad,jq,irow,idof1, 236 & idof2,value,i0,i0) 237 else 238 valu2=-coefmpc(index)*sm(jj,ll)/ 239 & coefmpc(ist) 240c 241 if(idof1.eq.idof2) valu2=2.d0*valu2 242c 243 call add_sm_ei(au,ad,aub,adb,jq,irow, 244 & idof1,idof2,value,valu2,i0,i0) 245 endif 246 endif 247 index=nodempc(3,index) 248 if(index.eq.0) exit 249 enddo 250 cycle 251 endif 252 endif 253! 254! regular DOF / SPC 255! 256 if(rhsi) then 257 elseif(nmethod.eq.2) then 258 value=s(jj,ll) 259c call nident(ikboun,idof2,nboun,id) 260c icolumn=neq(2)+ilboun(id) 261 icolumn=neq(2)-idof2/2 262 call add_bo_st(au,jq,irow,idof1,icolumn,value) 263 endif 264 else 265c idof1=(node1-1)*8+k 266c idof2=(node2-1)*8+m 267 idof1=jdof1 268 idof2=jdof2 269 mpc1=0 270 mpc2=0 271 if(nmpc.gt.0) then 272c call nident(ikmpc,idof1,nmpc,id1) 273c if((id1.gt.0).and.(ikmpc(id1).eq.idof1)) mpc1=1 274c call nident(ikmpc,idof2,nmpc,id2) 275c if((id2.gt.0).and.(ikmpc(id2).eq.idof2)) mpc2=1 276 if(idof1.ne.2*(idof1/2)) mpc1=1 277 if(idof2.ne.2*(idof2/2)) mpc2=1 278 endif 279 if((mpc1.eq.1).and.(mpc2.eq.1)) then 280c id1=ilmpc(id1) 281c id2=ilmpc(id2) 282 id1=(-idof1+1)/2 283 id2=(-idof2+1)/2 284 if(id1.eq.id2) then 285! 286! MPC id1 / MPC id1 287! 288 ist=ipompc(id1) 289 index1=nodempc(3,ist) 290 if(index1.eq.0) cycle 291 do 292 idof1=nactdof(nodempc(2,index1), 293 & nodempc(1,index1)) 294 index2=index1 295 do 296 idof2=nactdof(nodempc(2,index2), 297 & nodempc(1,index2)) 298 value=coefmpc(index1)*coefmpc(index2)* 299 & s(jj,ll)/coefmpc(ist)/coefmpc(ist) 300 if((idof1.gt.0).and.(idof2.gt.0)) then 301 if(stiffonly(1)) then 302 call add_sm_st(au,ad,jq,irow, 303 & idof1,idof2,value,i0,i0) 304 else 305 valu2=coefmpc(index1)*coefmpc(index2)* 306 & sm(jj,ll)/coefmpc(ist)/coefmpc(ist) 307 call add_sm_ei(au,ad,aub,adb,jq, 308 & irow,idof1,idof2,value,valu2,i0,i0) 309 endif 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 else 319! 320! MPC id1 / MPC id2 321! 322 ist1=ipompc(id1) 323 index1=nodempc(3,ist1) 324 if(index1.eq.0) cycle 325 do 326 idof1=nactdof(nodempc(2,index1), 327 & nodempc(1,index1)) 328 ist2=ipompc(id2) 329 index2=nodempc(3,ist2) 330 if(index2.eq.0) then 331 index1=nodempc(3,index1) 332 if(index1.eq.0) then 333 exit 334 else 335 cycle 336 endif 337 endif 338 do 339 idof2=nactdof(nodempc(2,index2), 340 & nodempc(1,index2)) 341 value=coefmpc(index1)*coefmpc(index2)* 342 & s(jj,ll)/coefmpc(ist1)/coefmpc(ist2) 343 if(idof1.eq.idof2) value=2.d0*value 344 if((idof1.gt.0).and.(idof2.gt.0)) then 345 if(stiffonly(1)) then 346 call add_sm_st(au,ad,jq,irow, 347 & idof1,idof2,value,i0,i0) 348 else 349 valu2=coefmpc(index1)*coefmpc(index2)* 350 & sm(jj,ll)/coefmpc(ist1)/coefmpc(ist2) 351c 352 if(idof1.eq.idof2) valu2=2.d0*valu2 353c 354 call add_sm_ei(au,ad,aub,adb,jq, 355 & irow,idof1,idof2,value,valu2,i0,i0) 356 endif 357 endif 358! 359 index2=nodempc(3,index2) 360 if(index2.eq.0) exit 361 enddo 362 index1=nodempc(3,index1) 363 if(index1.eq.0) exit 364 enddo 365 endif 366 endif 367 endif 368 enddo 369! 370 if(rhsi) then 371! 372! distributed forces 373! 374 if(idist.ne.0) then 375 if(jdof1.le.0) then 376 if(nmpc.ne.0) then 377c idof1=(node1-1)*8+k 378 idof1=jdof1 379c call nident(ikmpc,idof1,nmpc,id) 380c if((id.gt.0).and.(ikmpc(id).eq.idof1)) then 381 if(idof1.ne.2*(idof1/2)) then 382c id=ilmpc(id) 383 id=(-idof1+1)/2 384 ist=ipompc(id) 385 index=nodempc(3,ist) 386 if(index.eq.0) cycle 387 do 388 jdof1=nactdof(nodempc(2,index), 389 & nodempc(1,index)) 390 if(jdof1.gt.0) then 391 fext(jdof1)=fext(jdof1) 392 & -coefmpc(index)*ff(jj) 393 & /coefmpc(ist) 394 endif 395 index=nodempc(3,index) 396 if(index.eq.0) exit 397 enddo 398 endif 399 endif 400 cycle 401 endif 402 fext(jdof1)=fext(jdof1)+ff(jj) 403 endif 404 endif 405! 406 enddo 407 enddo 408! 409 endif 410 if(ithermal(1).gt.1) then 411! 412! thermal analysis: loop over all elements 413! 414 do i=1,ne 415! 416 if(ipkon(i).lt.0) cycle 417! 418! only elements belonging to the A-V-domain should be 419! included in the thermal analysis 420! 421 if(int(elcon(2,1,ielmat(1,i))).ne.2) cycle 422 indexe=ipkon(i) 423 if(lakon(i)(4:5).eq.'20') then 424 nope=20 425 elseif(lakon(i)(4:4).eq.'8') then 426 nope=8 427 elseif(lakon(i)(4:5).eq.'10') then 428 nope=10 429 elseif(lakon(i)(4:4).eq.'4') then 430 nope=4 431 elseif(lakon(i)(4:5).eq.'15') then 432 nope=15 433 elseif(lakon(i)(4:4).eq.'6') then 434 nope=6 435 elseif((lakon(i)(1:1).eq.'E').and.(lakon(i)(7:7).ne.'A')) then 436! 437! advection elements 438! 439 read(lakon(i)(8:8),'(i1)') nope 440 nope=nope+1 441 elseif((lakon(i)(1:2).eq.'D ').or. 442 & ((lakon(i)(1:1).eq.'D').and.(network.eq.1))) then 443! 444! asymmetrical contribution -> mafillsmas.f 445! 446 cycle 447 else 448 cycle 449 endif 450! 451 call e_c3d_th(co,nk,kon,lakon(i),st,smt, 452 & ff,i,nmethod,rhcon,nrhcon,ielmat,ielorien,norien,orab, 453 & ntmat_,t0,t1,ithermal,vold,iperturb,nelemload, 454 & sideload,xload,nload,idist,iexpl,dtime, 455 & matname,mi(1),mass(2),stiffness,buckling,rhsi,intscheme, 456 & physcon,shcon,nshcon,cocon,ncocon,ttime,time,istep,iinc, 457 & xstiff,xloadold,reltime,ipompc,nodempc,coefmpc,nmpc,ikmpc, 458 & ilmpc,springarea,plkcon,nplkcon,npmat_,ncmat_,elcon,nelcon, 459 & lakon,pslavsurf,pmastsurf,mortar,clearini,plicon,nplicon, 460 & ipkon,ielprop,prop,iponoel,inoel,sti,xstateini,xstate, 461 & nstate_,network,ipobody,xbody,ibody) 462! 463 do jj=1,nope 464! 465 j=jj 466! 467 node1=kon(indexe+j) 468 jdof1=nactdof(0,node1) 469! 470 do ll=jj,nope 471! 472 l=ll 473! 474 node2=kon(indexe+l) 475 jdof2=nactdof(0,node2) 476! 477! check whether one of the DOF belongs to a SPC or MPC 478! 479 if((jdof1.gt.0).and.(jdof2.gt.0)) then 480 if(stiffonly(2)) then 481 call add_sm_st(au,ad,jq,irow,jdof1,jdof2, 482 & st(jj,ll),jj,ll) 483 else 484 call add_sm_ei(au,ad,aub,adb,jq,irow,jdof1,jdof2, 485 & st(jj,ll),smt(jj,ll),jj,ll) 486 endif 487 elseif((jdof1.gt.0).or.(jdof2.gt.0)) then 488! 489! idof1: genuine DOF 490! idof2: nominal DOF of the SPC/MPC 491! 492 if(jdof1.le.0) then 493 idof1=jdof2 494c idof2=(node1-1)*8 495 idof2=jdof1 496 else 497 idof1=jdof1 498c idof2=(node2-1)*8 499 idof2=jdof2 500 endif 501 if(nmpc.gt.0) then 502c call nident(ikmpc,idof2,nmpc,id) 503c if((id.gt.0).and.(ikmpc(id).eq.idof2)) then 504 if(idof2.ne.2*(idof2/2)) then 505! 506! regular DOF / MPC 507! 508c id=ilmpc(id) 509 id=(-idof2+1)/2 510 ist=ipompc(id) 511 index=nodempc(3,ist) 512 if(index.eq.0) cycle 513 do 514 idof2=nactdof(nodempc(2,index),nodempc(1,index)) 515 value=-coefmpc(index)*st(jj,ll)/coefmpc(ist) 516 if(idof1.eq.idof2) value=2.d0*value 517 if(idof2.gt.0) then 518 if(stiffonly(2)) then 519 call add_sm_st(au,ad,jq,irow,idof1, 520 & idof2,value,i0,i0) 521 else 522 valu2=-coefmpc(index)*smt(jj,ll)/ 523 & coefmpc(ist) 524! 525 if(idof1.eq.idof2) valu2=2.d0*valu2 526! 527 call add_sm_ei(au,ad,aub,adb,jq,irow, 528 & idof1,idof2,value,valu2,i0,i0) 529 endif 530 endif 531 index=nodempc(3,index) 532 if(index.eq.0) exit 533 enddo 534 cycle 535 endif 536 endif 537! 538! regular DOF / SPC 539! 540 if(rhsi) then 541 elseif(nmethod.eq.2) then 542 value=st(jj,ll) 543c call nident(ikboun,idof2,nboun,id) 544c icolumn=neq(2)+ilboun(id) 545 icolumn=neq(2)-idof2/2 546 call add_bo_st(au,jq,irow,idof1,icolumn,value) 547 endif 548 else 549c idof1=(node1-1)*8 550c idof2=(node2-1)*8 551 idof1=jdof1 552 idof2=jdof2 553 mpc1=0 554 mpc2=0 555 if(nmpc.gt.0) then 556c call nident(ikmpc,idof1,nmpc,id1) 557c if((id1.gt.0).and.(ikmpc(id1).eq.idof1)) mpc1=1 558c call nident(ikmpc,idof2,nmpc,id2) 559c if((id2.gt.0).and.(ikmpc(id2).eq.idof2)) mpc2=1 560 if(idof1.ne.2*(idof1/2)) mpc1=1 561 if(idof2.ne.2*(idof2/2)) mpc2=1 562 endif 563 if((mpc1.eq.1).and.(mpc2.eq.1)) then 564c id1=ilmpc(id1) 565c id2=ilmpc(id2) 566 id1=(-idof1+1)/2 567 id2=(-idof2+1)/2 568 if(id1.eq.id2) then 569! 570! MPC id1 / MPC id1 571! 572 ist=ipompc(id1) 573 index1=nodempc(3,ist) 574 if(index1.eq.0) cycle 575 do 576 idof1=nactdof(nodempc(2,index1), 577 & nodempc(1,index1)) 578 index2=index1 579 do 580 idof2=nactdof(nodempc(2,index2), 581 & nodempc(1,index2)) 582 value=coefmpc(index1)*coefmpc(index2)* 583 & st(jj,ll)/coefmpc(ist)/coefmpc(ist) 584 if((idof1.gt.0).and.(idof2.gt.0)) then 585 if(stiffonly(2)) then 586 call add_sm_st(au,ad,jq,irow, 587 & idof1,idof2,value,i0,i0) 588 else 589 valu2=coefmpc(index1)*coefmpc(index2)* 590 & smt(jj,ll)/coefmpc(ist)/coefmpc(ist) 591 call add_sm_ei(au,ad,aub,adb,jq, 592 & irow,idof1,idof2,value,valu2,i0,i0) 593 endif 594 endif 595! 596 index2=nodempc(3,index2) 597 if(index2.eq.0) exit 598 enddo 599 index1=nodempc(3,index1) 600 if(index1.eq.0) exit 601 enddo 602 else 603! 604! MPC id1 / MPC id2 605! 606 ist1=ipompc(id1) 607 index1=nodempc(3,ist1) 608 if(index1.eq.0) cycle 609 do 610 idof1=nactdof(nodempc(2,index1), 611 & nodempc(1,index1)) 612 ist2=ipompc(id2) 613 index2=nodempc(3,ist2) 614 if(index2.eq.0) then 615 index1=nodempc(3,index1) 616 if(index1.eq.0) then 617 exit 618 else 619 cycle 620 endif 621 endif 622 do 623 idof2=nactdof(nodempc(2,index2), 624 & nodempc(1,index2)) 625 value=coefmpc(index1)*coefmpc(index2)* 626 & st(jj,ll)/coefmpc(ist1)/coefmpc(ist2) 627 if(idof1.eq.idof2) value=2.d0*value 628 if((idof1.gt.0).and.(idof2.gt.0)) then 629 if(stiffonly(2)) then 630 call add_sm_st(au,ad,jq,irow, 631 & idof1,idof2,value,i0,i0) 632 else 633 valu2=coefmpc(index1)*coefmpc(index2)* 634 & smt(jj,ll)/coefmpc(ist1)/coefmpc(ist2) 635! 636 if(idof1.eq.idof2) valu2=2.d0*valu2 637! 638 call add_sm_ei(au,ad,aub,adb,jq, 639 & irow,idof1,idof2,value,valu2,i0,i0) 640 endif 641 endif 642! 643 index2=nodempc(3,index2) 644 if(index2.eq.0) exit 645 enddo 646 index1=nodempc(3,index1) 647 if(index1.eq.0) exit 648 enddo 649 endif 650 endif 651 endif 652 enddo 653! 654 if(rhsi) then 655! 656! distributed forces 657! 658 if(idist.ne.0) then 659 if(jdof1.le.0) then 660 if(nmpc.ne.0) then 661c idof1=(node1-1)*8 662 idof1=jdof1 663c call nident(ikmpc,idof1,nmpc,id) 664c if((id.gt.0).and.(ikmpc(id).eq.idof1)) then 665 if(idof1.ne.2*(idof1/2)) then 666c id=ilmpc(id) 667 id=(-idof1+1)/2 668 ist=ipompc(id) 669 index=nodempc(3,ist) 670 if(index.eq.0) cycle 671 do 672 jdof1=nactdof(nodempc(2,index), 673 & nodempc(1,index)) 674 if(jdof1.gt.0) then 675 fext(jdof1)=fext(jdof1) 676 & -coefmpc(index)*ff(jj) 677 & /coefmpc(ist) 678 endif 679 index=nodempc(3,index) 680 if(index.eq.0) exit 681 enddo 682 endif 683 endif 684 cycle 685 endif 686 fext(jdof1)=fext(jdof1)+ff(jj) 687 endif 688 endif 689! 690 enddo 691 enddo 692! 693 endif 694! 695 return 696 end 697