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 mafillsmcs(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,bb,nactdof,icol,jq,irow,neq,nzl,nmethod, 24 & ikmpc,ilmpc,ikboun,ilboun,elcon,nelcon,rhcon, 25 & nrhcon,alcon,nalcon,alzero,ielmat,ielorien,norien,orab, 26 & ntmat_,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,ics,cs,nm,ncmat_,labmpc,mass,stiffness,buckling, 30 & rhsi,intscheme,mcs,coriolis,ibody,xloadold,reltime,ielcs, 31 & veold,springarea,thicke,integerglob,doubleglob,tieset, 32 & istartset,iendset,ialset,ntie,nasym,pslavsurf,pmastsurf, 33 & mortar,clearini,ielprop,prop,ne0,kscale,xstateini,xstate, 34 & nstate_,set,nset) 35! 36! filling the stiffness matrix in spare matrix format (sm) 37! for cyclic symmetry calculations 38! 39 implicit none 40! 41 logical mass,stiffness,buckling,rhsi,coriolis 42! 43 character*8 lakon(*) 44 character*20 labmpc(*),sideload(*) 45 character*80 matname(*) 46 character*81 tieset(3,*),set(*) 47! 48 integer kon(*),nodeboun(*),ndirboun(*),ipompc(*),nodempc(3,*), 49 & nodeforc(2,*),ndirforc(*),nelemload(2,*),icol(*),jq(*), 50 & ilmpc(*),ikboun(*),ilboun(*),mi(*),nstate_,ne0,ielprop(*), 51 & nactdof(0:mi(2),*),irow(*),istartset(*),iendset(*),kscale, 52 & nelcon(2,*),nrhcon(*),nalcon(2,*),ielmat(mi(3),*),nset, 53 & ielorien(mi(3),*),integerglob(*),ialset(*),ntie,ikmpc(*), 54 & ipkon(*),ics(*),ij,ilength,lprev,ipobody(2,*),nbody, 55 & ibody(3,*),nk,ne,nboun,nmpc,nforc,nload,neq,nzl,nmethod, 56 & ithermal(*),iprestr,iperturb(*),nzs,i,j,k,l,m,idist,jj, 57 & ll,id,id1,id2,ist,ist1,ist2,index,jdof1,jdof2,idof1,idof2, 58 & mpc1,mpc2,index1,index2,node1,node2,kflag,nasym,mortar, 59 & ntmat_,indexe,nope,norien,iexpl,i0,nm,inode,icomplex, 60 & inode1,icomplex1,inode2,icomplex2,ner,ncmat_,intscheme,istep, 61 & iinc,mcs,ielcs(*),nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*), 62 & npmat_,islavelinv(1),irowtloc(1),jqtloc(1),mortartrafoflag, 63 & mscalmethod 64! 65 real*8 co(3,*),xboun(*),coefmpc(*),xforc(*),xload(2,*),p1(3), 66 & p2(3),ad(*),au(*),bodyf(3),bb(*),xbody(7,*),cgr(4,*),prop(*), 67 & t0(*),t1(*),prestr(6,mi(1),*),vold(0:mi(2),*),s(60,60), 68 & xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*),ff(60), 69 & sti(6,mi(1),*),sm(60,60),stx(6,mi(1),*),adb(*),aub(*), 70 & elcon(0:ncmat_,ntmat_,*),rhcon(0:1,ntmat_,*),xloadold(2,*), 71 & alcon(0:6,ntmat_,*),cs(17,*),alzero(*),orab(7,*),reltime, 72 & springarea(2,*),plicon(0:2*npmat_,ntmat_,*),smscale(1), 73 & plkcon(0:2*npmat_,ntmat_,*),thicke(mi(3),*),doubleglob(*), 74 & xstiff(27,mi(1),*),pi,theta,ti,tr,veold(0:mi(2),*),om,valu2, 75 & value,dtime,walue,walu2,time,ttime,clearini(3,9,*), 76 & pslavsurf(3,*),pmastsurf(6,*),autloc(1) 77! 78 mortartrafoflag=0 79! 80! calculating the scaling factors for the cyclic symmetry calculation 81! 82 pi=4.d0*datan(1.d0) 83! 84 do i=1,mcs 85 write(*,*) '*INFO in mafillsmcs: calculating nodal diameter',nm 86 write(*,*) ' for',cs(1,i),'sectors' 87 write(*,*) ' (cyclic symmetry definition number',i,')' 88 write(*,*) 89 theta=nm*2.d0*pi/cs(1,i) 90 cs(15,i)=dcos(theta) 91 cs(16,i)=dsin(theta) 92 enddo 93! 94 kflag=2 95 i0=0 96! 97! determining nzl 98! 99 nzl=0 100 do i=neq,1,-1 101 if(icol(i).gt.0) then 102 nzl=i 103 exit 104 endif 105 enddo 106! 107! initializing the matrices 108! 109 do i=1,neq 110 ad(i)=0.d0 111 enddo 112 do i=1,nzs 113 au(i)=0.d0 114 enddo 115! 116 do i=1,neq 117 adb(i)=0.d0 118 enddo 119 do i=1,nzs 120 aub(i)=0.d0 121 enddo 122! 123 ner=neq/2 124! 125! loop over all elements 126! 127! initialisation of the error parameter 128! 129c ne0=0 130 do i=1,ne 131! 132 if(ipkon(i).lt.0) cycle 133 indexe=ipkon(i) 134c Bernhardi start 135 if(lakon(i)(4:5).eq.'8I') then 136 nope=11 137c Bernhardi end 138 elseif(lakon(i)(4:5).eq.'20') then 139 nope=20 140 elseif(lakon(i)(4:4).eq.'2') then 141 nope=26 142 elseif(lakon(i)(4:4).eq.'8') then 143 nope=8 144 elseif(lakon(i)(4:5).eq.'10') then 145 nope=10 146 elseif(lakon(i)(4:4).eq.'4') then 147 nope=4 148 elseif(lakon(i)(4:5).eq.'15') then 149 nope=15 150 elseif(lakon(i)(4:4).eq.'6') then 151 nope=6 152 elseif(lakon(i)(1:2).eq.'ES') then 153c begin 16.07.2020 154 nope=ichar(lakon(i)(8:8))-47 155c read(lakon(i)(8:8),'(i1)') nope 156c nope=nope+1 157c end 16.07.2020 158! 159! local contact spring number 160! 161c write(*,*) 'nope before= ',nope 162 if(lakon(i)(7:7).eq.'C') then 163 if(nasym.eq.1) cycle 164c begin 16.07.2020 165 if(mortar.eq.1) nope=kon(indexe) 166c end 16.07.2020 167 endif 168c write(*,*) 'nope after= ',nope 169 elseif(lakon(i)(1:4).eq.'MASS') then 170 nope=1 171 else 172 cycle 173 endif 174! 175 om=0.d0 176! 177 if((nbody.gt.0).and.(lakon(i)(1:1).ne.'E')) then 178! 179! assigning centrifugal forces 180! 181 index=i 182 do 183 j=ipobody(1,index) 184 if(j.eq.0) exit 185 if(ibody(1,j).eq.1) then 186 om=xbody(1,j) 187 p1(1)=xbody(2,j) 188 p1(2)=xbody(3,j) 189 p1(3)=xbody(4,j) 190 p2(1)=xbody(5,j) 191 p2(2)=xbody(6,j) 192 p2(3)=xbody(7,j) 193 endif 194 index=ipobody(2,index) 195 if(index.eq.0) exit 196 enddo 197 endif 198! 199 call e_c3d(co,kon,lakon(i),p1,p2,om,bodyf,nbody,s,sm,ff,i, 200 & nmethod,elcon,nelcon,rhcon,nrhcon,alcon,nalcon, 201 & alzero,ielmat,ielorien,norien,orab,ntmat_, 202 & t0,t1,ithermal,vold,iperturb,nelemload,sideload,xload, 203 & nload,idist,sti,stx,iexpl,plicon, 204 & nplicon,plkcon,nplkcon,xstiff,npmat_, 205 & dtime,matname,mi(1),ncmat_,mass,stiffness,buckling,rhsi, 206 & intscheme,ttime,time,istep,iinc,coriolis,xloadold, 207 & reltime,ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,veold, 208 & springarea,nstate_,xstateini,xstate,ne0,ipkon,thicke, 209 & integerglob,doubleglob,tieset,istartset, 210 & iendset,ialset,ntie,nasym,pslavsurf,pmastsurf,mortar, 211 & clearini,ielprop,prop,kscale,smscale(1),mscalmethod, 212 & set,nset,islavelinv,autloc, 213 & irowtloc,jqtloc,mortartrafoflag) 214! 215 do jj=1,3*nope 216! 217 j=(jj-1)/3+1 218 k=jj-3*(j-1) 219! 220 node1=kon(indexe+j) 221 jdof1=nactdof(k,node1) 222! 223 do ll=jj,3*nope 224 if (mcs.gt.1)then 225 if(ielcs(i).gt.0) then 226 s(jj,ll)=(cs(1,(ielcs(i)+1))/cs(1,1))*s(jj,ll) 227 sm(jj,ll)=(cs(1,(ielcs(i)+1))/cs(1,1))*sm(jj,ll) 228 endif 229 endif 230! 231 l=(ll-1)/3+1 232 m=ll-3*(l-1) 233! 234 node2=kon(indexe+l) 235 jdof2=nactdof(m,node2) 236! 237! check whether one of the DOF belongs to a SPC or MPC 238! 239 if((jdof1.gt.0).and.(jdof2.gt.0)) then 240 call add_sm_ei(au,ad,aub,adb,jq,irow,jdof1,jdof2, 241 & s(jj,ll),sm(jj,ll),jj,ll) 242 call add_sm_ei(au,ad,aub,adb,jq,irow,jdof1+ner,jdof2+ner, 243 & s(jj,ll),sm(jj,ll),jj,ll) 244 elseif((jdof1.gt.0).or.(jdof2.gt.0)) then 245! 246! idof1: genuine DOF 247! idof2: nominal DOF of the SPC/MPC 248! 249 if(jdof1.le.0) then 250 idof1=jdof2 251 idof2=jdof1 252 else 253 idof1=jdof1 254 idof2=jdof2 255 endif 256! 257 if(nmpc.gt.0) then 258 if(idof2.ne.2*(idof2/2)) then 259! 260! regular DOF / MPC 261! 262 id1=(-idof2+1)/2 263 ist=ipompc(id1) 264 index=nodempc(3,ist) 265 if(index.eq.0) cycle 266 do 267 inode=nodempc(1,index) 268 icomplex=0 269 if(labmpc(id1)(1:6).eq.'CYCLIC') then 270 read(labmpc(id1)(7:20),'(i14)') icomplex 271 elseif(labmpc(id1)(1:9).eq.'SUBCYCLIC') then 272 do ij=1,mcs 273 ilength=int(cs(4,ij)) 274 lprev=int(cs(14,ij)) 275 call nident(ics(lprev+1),inode,ilength,id) 276 if(id.gt.0) then 277 if(ics(lprev+id).eq.inode) then 278 icomplex=ij 279 exit 280 endif 281 endif 282 enddo 283 endif 284 idof2=nactdof(nodempc(2,index),inode) 285 if(idof2.gt.0) then 286 value=-coefmpc(index)*s(jj,ll)/coefmpc(ist) 287 valu2=-coefmpc(index)*sm(jj,ll)/ 288 & coefmpc(ist) 289 if(idof1.eq.idof2) then 290 value=2.d0*value 291 valu2=2.d0*valu2 292 endif 293 if(icomplex.eq.0) then 294 call add_sm_ei(au,ad,aub,adb,jq,irow, 295 & idof1,idof2,value,valu2,i0,i0) 296 call add_sm_ei(au,ad,aub,adb,jq,irow, 297 & idof1+ner,idof2+ner,value,valu2,i0,i0) 298 else 299 walue=value*cs(15,icomplex) 300 walu2=valu2*cs(15,icomplex) 301 call add_sm_ei(au,ad,aub,adb,jq,irow, 302 & idof1,idof2,walue,walu2,i0,i0) 303 call add_sm_ei(au,ad,aub,adb,jq,irow, 304 & idof1+ner,idof2+ner,walue,walu2,i0,i0) 305 if(idof1.ne.idof2) then 306 walue=value*cs(16,icomplex) 307 walu2=valu2*cs(16,icomplex) 308 call add_sm_ei(au,ad,aub,adb,jq,irow, 309 & idof1,idof2+ner,walue,walu2,i0,i0) 310 walue=-walue 311 walu2=-walu2 312 call add_sm_ei(au,ad,aub,adb,jq,irow, 313 & idof1+ner,idof2,walue,walu2,i0,i0) 314 endif 315 endif 316 endif 317 index=nodempc(3,index) 318 if(index.eq.0) exit 319 enddo 320 cycle 321 endif 322 endif 323! 324 else 325 idof1=jdof1 326 idof2=jdof2 327! 328 mpc1=0 329 mpc2=0 330 if(nmpc.gt.0) then 331 if(idof1.ne.2*(idof1/2)) mpc1=1 332 if(idof2.ne.2*(idof2/2)) mpc2=1 333 endif 334 if((mpc1.eq.1).and.(mpc2.eq.1)) then 335 id1=(-idof1+1)/2 336 id2=(-idof2+1)/2 337 if(id1.eq.id2) then 338! 339! MPC id1 / MPC id1 340! 341 ist=ipompc(id1) 342 index1=nodempc(3,ist) 343 if(index1.eq.0) cycle 344 do 345 inode1=nodempc(1,index1) 346 icomplex1=0 347 if(labmpc(id1)(1:6).eq.'CYCLIC') then 348 read(labmpc(id1)(7:20),'(i14)') icomplex1 349 elseif(labmpc(id1)(1:9).eq.'SUBCYCLIC') then 350 do ij=1,mcs 351 ilength=int(cs(4,ij)) 352 lprev=int(cs(14,ij)) 353 call nident(ics(lprev+1),inode1, 354 & ilength,id) 355 if(id.gt.0) then 356 if(ics(lprev+id).eq.inode1) then 357 icomplex1=ij 358 exit 359 endif 360 endif 361 enddo 362 endif 363 idof1=nactdof(nodempc(2,index1),inode1) 364 index2=index1 365 do 366 inode2=nodempc(1,index2) 367 icomplex2=0 368 if(labmpc(id1)(1:6).eq.'CYCLIC') then 369 read(labmpc(id1)(7:20),'(i14)') icomplex2 370 elseif(labmpc(id1)(1:9).eq.'SUBCYCLIC') then 371 do ij=1,mcs 372 ilength=int(cs(4,ij)) 373 lprev=int(cs(14,ij)) 374 call nident(ics(lprev+1),inode2, 375 & ilength,id) 376 if(id.gt.0) then 377 if(ics(lprev+id).eq.inode2) then 378 icomplex2=ij 379 exit 380 endif 381 endif 382 enddo 383 endif 384 idof2=nactdof(nodempc(2,index2),inode2) 385 if((idof1.gt.0).and.(idof2.gt.0)) then 386 value=coefmpc(index1)*coefmpc(index2)* 387 & s(jj,ll)/coefmpc(ist)/coefmpc(ist) 388 valu2=coefmpc(index1)*coefmpc(index2)* 389 & sm(jj,ll)/coefmpc(ist)/coefmpc(ist) 390 if((icomplex1.eq.0).and. 391 & (icomplex2.eq.0)) then 392 call add_sm_ei(au,ad,aub,adb,jq, 393 & irow,idof1,idof2,value,valu2,i0,i0) 394 call add_sm_ei(au,ad,aub,adb,jq, 395 & irow,idof1+ner,idof2+ner,value, 396 & valu2,i0,i0) 397 elseif((icomplex1.ne.0).and. 398 & (icomplex2.ne.0)) then 399 if(icomplex1.eq.icomplex2) then 400 call add_sm_ei(au,ad,aub,adb,jq, 401 & irow,idof1,idof2,value,valu2,i0,i0) 402 call add_sm_ei(au,ad,aub,adb,jq, 403 & irow,idof1+ner,idof2+ner,value, 404 & valu2,i0,i0) 405 else 406 tr=cs(15,icomplex1)*cs(15,icomplex2) 407 & +cs(16,icomplex1)*cs(16,icomplex2) 408 walue=value*tr 409 walu2=valu2*tr 410 call add_sm_ei(au,ad,aub,adb,jq, 411 & irow,idof1,idof2,walue,walu2,i0,i0) 412 call add_sm_ei(au,ad,aub,adb,jq, 413 & irow,idof1+ner,idof2+ner,walue, 414 & walu2,i0,i0) 415 ti=cs(15,icomplex1)*cs(16,icomplex2) 416 & -cs(15,icomplex2)*cs(16,icomplex1) 417 walue=value*ti 418 walu2=valu2*ti 419 call add_sm_ei(au,ad,aub,adb,jq,irow 420 & ,idof1,idof2+ner,walue,walu2,i0,i0) 421 walue=-walue 422 walu2=-walu2 423 call add_sm_ei(au,ad,aub,adb,jq,irow 424 & ,idof1+ner,idof2,walue,walu2,i0,i0) 425 endif 426 elseif((icomplex1.eq.0).or. 427 & (icomplex2.eq.0)) then 428 if(icomplex2.ne.0) then 429 walue=value*cs(15,icomplex2) 430 walu2=valu2*cs(15,icomplex2) 431 else 432 walue=value*cs(15,icomplex1) 433 walu2=valu2*cs(15,icomplex1) 434 endif 435 call add_sm_ei(au,ad,aub,adb,jq,irow, 436 & idof1,idof2,walue,walu2,i0,i0) 437 call add_sm_ei(au,ad,aub,adb,jq,irow, 438 & idof1+ner,idof2+ner,walue,walu2,i0,i0) 439 if(icomplex2.ne.0) then 440 walue=value*cs(16,icomplex2) 441 walu2=valu2*cs(16,icomplex2) 442 else 443 walue=-value*cs(16,icomplex1) 444 walu2=-valu2*cs(16,icomplex1) 445 endif 446 call add_sm_ei(au,ad,aub,adb,jq,irow, 447 & idof1,idof2+ner,walue,walu2,i0,i0) 448 walue=-walue 449 walu2=-walu2 450 call add_sm_ei(au,ad,aub,adb,jq,irow, 451 & idof1+ner,idof2,walue,walu2,i0,i0) 452 endif 453 endif 454 index2=nodempc(3,index2) 455 if(index2.eq.0) exit 456 enddo 457 index1=nodempc(3,index1) 458 if(index1.eq.0) exit 459 enddo 460 else 461! 462! MPC id1 / MPC id2 463! 464 ist1=ipompc(id1) 465 index1=nodempc(3,ist1) 466 if(index1.eq.0) cycle 467 do 468 inode1=nodempc(1,index1) 469 icomplex1=0 470 if(labmpc(id1)(1:6).eq.'CYCLIC') then 471 read(labmpc(id1)(7:20),'(i14)') icomplex1 472 elseif(labmpc(id1)(1:9).eq.'SUBCYCLIC') then 473 do ij=1,mcs 474 ilength=int(cs(4,ij)) 475 lprev=int(cs(14,ij)) 476 call nident(ics(lprev+1),inode1, 477 & ilength,id) 478 if(id.gt.0) then 479 if(ics(lprev+id).eq.inode1) then 480 icomplex1=ij 481 exit 482 endif 483 endif 484 enddo 485 endif 486 idof1=nactdof(nodempc(2,index1),inode1) 487 ist2=ipompc(id2) 488 index2=nodempc(3,ist2) 489 if(index2.eq.0) then 490 index1=nodempc(3,index1) 491 if(index1.eq.0) then 492 exit 493 else 494 cycle 495 endif 496 endif 497 do 498 inode2=nodempc(1,index2) 499 icomplex2=0 500 if(labmpc(id2)(1:6).eq.'CYCLIC') then 501 read(labmpc(id2)(7:20),'(i14)') icomplex2 502 elseif(labmpc(id2)(1:9).eq.'SUBCYCLIC') then 503 do ij=1,mcs 504 ilength=int(cs(4,ij)) 505 lprev=int(cs(14,ij)) 506 call nident(ics(lprev+1),inode2, 507 & ilength,id) 508 if(id.gt.0) then 509 if(ics(lprev+id).eq.inode2) then 510 icomplex2=ij 511 exit 512 endif 513 endif 514 enddo 515 endif 516 idof2=nactdof(nodempc(2,index2),inode2) 517 if((idof1.gt.0).and.(idof2.gt.0)) then 518 value=coefmpc(index1)*coefmpc(index2)* 519 & s(jj,ll)/coefmpc(ist1)/coefmpc(ist2) 520 valu2=coefmpc(index1)*coefmpc(index2)* 521 & sm(jj,ll)/coefmpc(ist1)/coefmpc(ist2) 522 if(idof1.eq.idof2) then 523 value=2.d0*value 524 valu2=2.d0*valu2 525 endif 526 if((icomplex1.eq.0).and. 527 & (icomplex2.eq.0)) then 528 call add_sm_ei(au,ad,aub,adb,jq, 529 & irow,idof1,idof2,value,valu2,i0,i0) 530 call add_sm_ei(au,ad,aub,adb,jq, 531 & irow,idof1+ner,idof2+ner,value, 532 & valu2,i0,i0) 533 elseif((icomplex1.ne.0).and. 534 & (icomplex2.ne.0)) then 535 if(icomplex1.eq.icomplex2) then 536 call add_sm_ei(au,ad,aub,adb,jq, 537 & irow,idof1,idof2,value,valu2,i0,i0) 538 call add_sm_ei(au,ad,aub,adb,jq, 539 & irow,idof1+ner,idof2+ner,value, 540 & valu2,i0,i0) 541 else 542 tr=cs(15,icomplex1)*cs(15,icomplex2) 543 & +cs(16,icomplex1)*cs(16,icomplex2) 544c write(*,*) 'tr= ',tr 545 walue=value*tr 546 walu2=valu2*tr 547 call add_sm_ei(au,ad,aub,adb,jq, 548 & irow,idof1,idof2,walue,walu2,i0,i0) 549 call add_sm_ei(au,ad,aub,adb,jq, 550 & irow,idof1+ner,idof2+ner,walue, 551 & walu2,i0,i0) 552 ti=cs(15,icomplex1)*cs(16,icomplex2) 553 & -cs(15,icomplex2)*cs(16,icomplex1) 554 walue=value*ti 555 walu2=valu2*ti 556 call add_sm_ei(au,ad,aub,adb,jq,irow 557 & ,idof1,idof2+ner,walue,walu2,i0,i0) 558 walue=-walue 559 walu2=-walu2 560 call add_sm_ei(au,ad,aub,adb,jq,irow 561 & ,idof1+ner,idof2,walue,walu2,i0,i0) 562 endif 563 elseif((icomplex1.eq.0).or. 564 & (icomplex2.eq.0)) then 565 if(icomplex2.ne.0) then 566 walue=value*cs(15,icomplex2) 567 walu2=valu2*cs(15,icomplex2) 568 else 569 walue=value*cs(15,icomplex1) 570 walu2=valu2*cs(15,icomplex1) 571 endif 572 call add_sm_ei(au,ad,aub,adb,jq,irow, 573 & idof1,idof2,walue,walu2,i0,i0) 574 call add_sm_ei(au,ad,aub,adb,jq,irow, 575 & idof1+ner,idof2+ner,walue,walu2,i0,i0) 576 if(idof1.ne.idof2) then 577 if(icomplex2.ne.0) then 578 walue=value*cs(16,icomplex2) 579 walu2=valu2*cs(16,icomplex2) 580 else 581 walue=-value*cs(16,icomplex1) 582 walu2=-valu2*cs(16,icomplex1) 583 endif 584 call add_sm_ei(au,ad,aub,adb,jq, 585 & irow,idof1,idof2+ner,walue, 586 & walu2,i0,i0) 587 walue=-walue 588 walu2=-walu2 589 call add_sm_ei(au,ad,aub,adb,jq, 590 & irow,idof1+ner,idof2,walue, 591 & walu2,i0,i0) 592 endif 593 endif 594 endif 595 index2=nodempc(3,index2) 596 if(index2.eq.0) exit 597 enddo 598 index1=nodempc(3,index1) 599 if(index1.eq.0) exit 600 enddo 601 endif 602 endif 603 endif 604 enddo 605! 606 enddo 607 enddo 608! 609 return 610 end 611