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 couplings(inpc,textpart,set,istartset,iendset, 20 & ialset,nset,nboun,nk,ipompc,nodempc,coefmpc,nmpc,nmpc_, 21 & mpcfree,ikboun,ikmpc,ilmpc,co,labmpc,istat,n,iline,ipol, 22 & inl,ipoinp,inp,ipoinpc,norien,orname,orab,irstrt,ipkon, 23 & kon,lakon,istep,ics,dcs,nk_,nboun_,nodeboun,ndirboun, 24 & typeboun,ilboun,xboun,ier) 25! 26! reading the input deck: *COUPLING in combination with 27! *KINEMATIC or *DISTRIBUTING 28! 29 implicit none 30! 31 character*1 inpc(*),surfkind,typeboun(*) 32 character*8 lakon(*) 33 character*20 labmpc(*),label 34 character*80 orname(*),orientation 35 character*81 set(*),surfset 36 character*132 textpart(16),name 37! 38 integer istartset(*),iendset(*),ialset(*),norien,irstrt(*),nface, 39 & iorientation,iface,jface,nset,nboun,istat,n,i,j,k,ibounstart, 40 & ibounend,key,nk,ipompc(*),nodempc(3,*),nmpc,nmpc_,mpcfree, 41 & ikboun(*),ikmpc(*),ilmpc(*),ipos,m,node,iline,ipol,inl,nope, 42 & ipoinp(2,*),inp(3,*),ipoinpc(0:*),jsurf,irefnode,indexe,nopes, 43 & ifaceq(8,6),ifacet(6,4),ifacew1(4,5),ifacew2(8,5),ipkon(*), 44 & kon(*),istep,nelem,ics(2,*),nodef(8),nboun_,nodeboun(*), 45 & npt,mint,index1,mpcfreeold,id,idof,iflag,inhomnode,nk_,kk, 46 & irotnode(3),l,index2,indexold(3),indexnew(3),idirold,idirmax, 47 & idir,indexmax,nodeold,node1,nodemax,ndirboun(*),ilboun(*), 48 & irotnode_kin,idupnode,ier 49! 50 real*8 coefmpc(*),co(3,*),orab(7,*),dcs(*),areanodal(8),xl2(3,8), 51 & shp2(7,8),xsj2(3),xsj,xi,et,weight,xs2(3,2),area,a(3,3),cgx(3), 52 & aa(3),pi(3),c1,c4,c9,c10,amax,xcoef,coef(3),xboun(*),stdev 53! 54 include "gauss.f" 55! 56! nodes per face for hex elements 57! 58 data ifaceq /4,3,2,1,11,10,9,12, 59 & 5,6,7,8,13,14,15,16, 60 & 1,2,6,5,9,18,13,17, 61 & 2,3,7,6,10,19,14,18, 62 & 3,4,8,7,11,20,15,19, 63 & 4,1,5,8,12,17,16,20/ 64! 65! nodes per face for tet elements 66! 67 data ifacet /1,3,2,7,6,5, 68 & 1,2,4,5,9,8, 69 & 2,3,4,6,10,9, 70 & 1,4,3,8,10,7/ 71! 72! nodes per face for linear wedge elements 73! 74 data ifacew1 /1,3,2,0, 75 & 4,5,6,0, 76 & 1,2,5,4, 77 & 2,3,6,5, 78 & 3,1,4,6/ 79! 80! nodes per face for quadratic wedge elements 81! 82 data ifacew2 /1,3,2,9,8,7,0,0, 83 & 4,5,6,10,11,12,0,0, 84 & 1,2,5,4,7,14,10,13, 85 & 2,3,6,5,8,15,11,14, 86 & 3,1,4,6,9,13,12,15/ 87! 88! flag for shape functions 89! 90 data iflag /2/ 91! 92 if((istep.gt.0).and.(irstrt(1).ge.0)) then 93 write(*,*) '*ERROR reading *COUPLING: *COUPLING' 94 write(*,*)' should be placed before all step definitions' 95 ier=1 96 return 97 endif 98! 99 label=' ' 100 orientation=' 101 & ' 102 do i=1,81 103 surfset(i:i)=' ' 104 enddo 105! 106 name(1:1)=' ' 107 do i=2,n 108 if(textpart(i)(1:8).eq.'REFNODE=') then 109 read(textpart(i)(9:18),'(i10)',iostat=istat) irefnode 110 if(istat.gt.0) then 111 call inputerror(inpc,ipoinpc,iline, 112 & "*COUPLING%",ier) 113 return 114 endif 115 if(irefnode.gt.nk) then 116 write(*,*) '*ERROR reading *COUPLING: ref node',irefnode 117 write(*,*) ' has not been defined' 118 ier=1 119 return 120 endif 121 else if(textpart(i)(1:8).eq.'SURFACE=') then 122 surfset(1:80)=textpart(i)(9:88) 123! 124 ipos=index(surfset,' ') 125 surfkind='S' 126 surfset(ipos:ipos)=surfkind 127c do j=1,nset 128c if(set(j).eq.surfset) exit 129c enddo 130 call cident81(set,surfset,nset,id) 131 j=nset+1 132 if(id.gt.0) then 133 if(surfset.eq.set(id)) then 134 j=id 135 endif 136 endif 137 if(j.gt.nset) then 138 surfkind='T' 139 surfset(ipos:ipos)=surfkind 140c do j=1,nset 141c if(set(j).eq.surfset) exit 142c enddo 143 call cident81(set,surfset,nset,id) 144 j=nset+1 145 if(id.gt.0) then 146 if(surfset.eq.set(id)) then 147 j=id 148 endif 149 endif 150 if(j.gt.nset) then 151 write(*,*) '*ERROR reading *COUPLING:' 152 write(*,*) ' surface ',surfset 153 write(*,*) ' has not yet been defined.' 154 ier=1 155 return 156 endif 157 endif 158 jsurf=j 159 elseif(textpart(i)(1:12).eq.'ORIENTATION=') then 160 orientation=textpart(i)(13:92) 161 elseif(textpart(i)(1:15).eq.'CONSTRAINTNAME=') then 162 name(1:117)=textpart(i)(16:132) 163 else 164 write(*,*) 165 & '*WARNING reading *COUPLING: parameter not recognized:' 166 write(*,*) ' ', 167 & textpart(i)(1:index(textpart(i),' ')-1) 168 call inputwarning(inpc,ipoinpc,iline, 169 &"*COUPLING%") 170 endif 171 enddo 172! 173 if(name(1:1).eq.' ') then 174 write(*,*) 175 & '*ERROR reading *COUPLING: no CONTRAINT NAME given' 176 write(*,*) ' ' 177 call inputerror(inpc,ipoinpc,iline, 178 & "*COUPLING%",ier) 179 return 180 endif 181! 182 if(orientation.eq.' ') then 183 iorientation=0 184 else 185 do i=1,norien 186 if(orname(i).eq.orientation) exit 187 enddo 188 if(i.gt.norien) then 189 write(*,*) 190 & '*ERROR reading *COUPLING: nonexistent orientation' 191 write(*,*) ' ' 192 call inputerror(inpc,ipoinpc,iline, 193 & "*COUPLING%",ier) 194 return 195 endif 196 iorientation=i 197 endif 198! 199! next keyword should be *KINEMATIC or *DISTRIBUTING 200! 201 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 202 & ipoinp,inp,ipoinpc) 203 if(textpart(1)(2:10).eq.'KINEMATIC') then 204! 205! catalogueing the nodes 206! 207 npt=0 208! 209 do j=istartset(jsurf),iendset(jsurf) 210 if(ialset(j).gt.0) then 211 if(surfkind.eq.'T') then 212! 213! facial surface 214! 215 iface=ialset(j) 216 nelem=int(iface/10) 217 jface=iface-nelem*10 218 indexe=ipkon(nelem) 219! 220 if(lakon(nelem)(4:5).eq.'20') then 221 nopes=8 222 elseif(lakon(nelem)(4:4).eq.'2') then 223 nopes=9 224 elseif(lakon(nelem)(4:4).eq.'8') then 225 nopes=4 226 elseif(lakon(nelem)(4:5).eq.'10') then 227 nopes=6 228 elseif(lakon(nelem)(4:4).eq.'4') then 229 nopes=3 230 endif 231! 232 if(lakon(nelem)(4:4).eq.'6') then 233 if(jface.le.2) then 234 nopes=3 235 else 236 nopes=4 237 endif 238 endif 239 if(lakon(nelem)(4:5).eq.'15') then 240 if(jface.le.2) then 241 nopes=6 242 else 243 nopes=8 244 endif 245 endif 246 else 247! 248! nodal surface 249! 250 nopes=1 251 endif 252! 253 do m=1,nopes 254 if(surfkind.eq.'T') then 255 if((lakon(nelem)(4:4).eq.'2').or. 256 & (lakon(nelem)(4:4).eq.'8')) then 257 node=kon(indexe+ifaceq(m,jface)) 258 elseif((lakon(nelem)(4:4).eq.'4').or. 259 & (lakon(nelem)(4:5).eq.'10')) then 260 node=kon(indexe+ifacet(m,jface)) 261 elseif(lakon(nelem)(4:4).eq.'6') then 262 node=kon(indexe+ifacew1(m,jface)) 263 elseif(lakon(nelem)(4:5).eq.'15') then 264 node=kon(indexe+ifacew2(m,jface)) 265 endif 266 else 267 node =ialset(j) 268 endif 269! 270 call nident2(ics,node,npt,id) 271 if(id.gt.0) then 272 if(ics(1,id).eq.node) then 273 cycle 274 endif 275 endif 276! 277! updating ics 278! 279 npt=npt+1 280 do l=npt,id+2,-1 281 ics(1,l)=ics(1,l-1) 282 enddo 283 ics(1,id+1)=node 284 enddo 285 else 286! 287! if a negative value occurs the surface has to be 288! nodal 289! 290 k=ialset(j-2) 291 do 292 k=k-ialset(j) 293 if(k.ge.ialset(j-1)) exit 294 node=k 295 call nident2(ics,node,npt,id) 296 if(id.gt.0) then 297 if(ics(1,id).eq.node) then 298 cycle 299 endif 300 endif 301! 302! updating ics 303! 304 npt=npt+1 305 do l=npt,id+2,-1 306 ics(1,l)=ics(1,l-1) 307 enddo 308 ics(1,id+1)=node 309 enddo 310 endif 311 enddo 312! 313! generating a rotational node and connecting the 314! rotational dofs of the reference node with the 315! translational dofs of the rotational node 316! 317! generating a rotational reference node 318! 319 nk=nk+1 320 if(nk.gt.nk_) then 321 write(*,*) 322 & '*ERROR reading *KINEMATIC: increase nk_' 323 ier=1 324 return 325 endif 326 irotnode_kin=nk 327 do l=1,3 328 co(l,nk)=co(l,irefnode) 329 enddo 330! 331! generating connecting MPCs between the rotational 332! dofs of irefnode and the translational dofs of 333! irotnode_kin 334! 335 do k=1,3 336! 337 nmpc=nmpc+1 338 if(nmpc.gt.nmpc_) then 339 write(*,*) 340 & '*ERROR reading *KINEMATIC: increase nmpc_' 341 ier=1 342 return 343 endif 344! 345! the internal dofs for rotation are 4, 5 and 6 346! 347 ipompc(nmpc)=mpcfree 348 labmpc(nmpc)='ROTTRACOUPLING ' 349 idof=8*(irefnode-1)+k+3 350 call nident(ikmpc,idof,nmpc-1,id) 351 do l=nmpc,id+2,-1 352 ikmpc(l)=ikmpc(l-1) 353 ilmpc(l)=ilmpc(l-1) 354 enddo 355 ikmpc(id+1)=idof 356 ilmpc(id+1)=nmpc 357! 358 nodempc(1,mpcfree)=irefnode 359 nodempc(2,mpcfree)=k+4 360 coefmpc(mpcfree)=1.d0 361 mpcfree=nodempc(3,mpcfree) 362! 363 nodempc(1,mpcfree)=irotnode_kin 364 nodempc(2,mpcfree)=k 365 coefmpc(mpcfree)=-1.d0 366 mpcfreeold=mpcfree 367 mpcfree=nodempc(3,mpcfree) 368 nodempc(3,mpcfreeold)=0 369 enddo 370! 371 if(iorientation.gt.0) then 372! 373! duplicating the nodes 374! generating rigid body MPC's for all dofs in the new nodes 375! 376 do m=1,npt 377 nk=nk+1 378 if(nk.gt.nk_) then 379 write(*,*) 380 & '*ERROR reading *KINEMATIC: increase nk_' 381 ier=1 382 return 383 endif 384 ics(2,m)=nk 385 do k=1,3 386 co(k,nk)=co(k,ics(1,m)) 387 enddo 388 node=nk 389 ibounstart=1 390 ibounend=3 391 call rigidmpc(ipompc,nodempc,coefmpc,irefnode, 392 & irotnode_kin, 393 & labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,nk,nk_, 394 & nodeboun,ndirboun,ikboun,ilboun,nboun,nboun_,node, 395 & typeboun,co,ibounstart,ibounend) 396 enddo 397 endif 398! 399! reading the degrees of freedom 400! 401 do 402 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 403 & ipoinp,inp,ipoinpc) 404 if((istat.lt.0).or.(key.eq.1)) return 405! 406 read(textpart(1)(1:10),'(i10)',iostat=istat) ibounstart 407 if(istat.gt.0) then 408 call inputerror(inpc,ipoinpc,iline, 409 & "*KINEMATIC%",ier) 410 return 411 endif 412 if(ibounstart.lt.1) then 413 write(*,*) '*ERROR reading *KINEMATIC' 414 write(*,*) ' starting degree of freedom cannot' 415 write(*,*) ' be less than 1' 416 write(*,*) ' ' 417 call inputerror(inpc,ipoinpc,iline, 418 & "*KINEMATIC%",ier) 419 return 420 endif 421! 422 if(textpart(2)(1:1).eq.' ') then 423 ibounend=ibounstart 424 else 425 read(textpart(2)(1:10),'(i10)',iostat=istat) ibounend 426 if(istat.gt.0) then 427 call inputerror(inpc,ipoinpc,iline, 428 & "*BOUNDARY%",ier) 429 return 430 endif 431 endif 432 if(ibounend.gt.3) then 433 write(*,*) '*ERROR reading *KINEMATIC' 434 write(*,*) ' final degree of freedom cannot' 435 write(*,*) ' exceed 3' 436 write(*,*) ' ' 437 call inputerror(inpc,ipoinpc,iline, 438 & "*KINEMATIC%",ier) 439 return 440 elseif(ibounend.lt.ibounstart) then 441 write(*,*) '*ERROR reading *KINEMATIC' 442 write(*,*) ' initial degree of freedom cannot' 443 write(*,*) ' exceed final degree of freedom' 444 write(*,*) ' ' 445 call inputerror(inpc,ipoinpc,iline, 446 & "*KINEMATIC%",ier) 447 return 448 endif 449! 450! generating the MPCs 451! 452 if(iorientation.eq.0) then 453! 454! generating rigid body MPC's for the appropriate dofs 455! 456 do j=1,npt 457 node=ics(1,j) 458 call rigidmpc(ipompc,nodempc,coefmpc,irefnode, 459 & irotnode_kin, 460 & labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,nk,nk_, 461 & nodeboun,ndirboun,ikboun,ilboun,nboun,nboun_, 462 & node,typeboun,co,ibounstart,ibounend) 463 enddo 464 else 465! 466! connecting the original nodes with the duplicated nodes 467! for the appropriate dofs 468! 469 do j=1,npt 470 node=ics(1,j) 471 idupnode=ics(2,j) 472 call mpcadd(node,ibounstart,ibounend,nboun,ipompc, 473 & nodempc,coefmpc,nmpc,nmpc_,mpcfree,orab,ikboun, 474 & ikmpc,ilmpc,co,labmpc,label,idupnode, 475 & iorientation) 476 enddo 477 endif 478 enddo 479 elseif(textpart(1)(2:13).eq.'DISTRIBUTING') then 480 if(surfkind.eq.'S') then 481 write(*,*) '*ERROR reading *DISTRIBUTING' 482 write(*,*) ' a nodal surface is not allowed' 483 write(*,*) ' please use a facial surface on' 484 write(*,*) ' the *COUPLING card' 485 ier=1 486 return 487 endif 488! 489 npt=0 490 area=0.d0 491! 492! catalogueing the nodes belonging to the surface (ics(1,*)) 493! catalogueing the area (dcs(*)) 494! 495 do k=istartset(jsurf),iendset(jsurf) 496! 497! facial surface 498! 499 iface=ialset(k) 500 nelem=int(iface/10) 501 jface=iface-nelem*10 502 indexe=ipkon(nelem) 503! 504! nodes: #nodes in the face 505! the nodes are stored in nodef(*) 506! 507 if(lakon(nelem)(4:4).eq.'2') then 508 nopes=8 509 nface=6 510 elseif(lakon(nelem)(3:4).eq.'D8') then 511 nopes=4 512 nface=6 513 elseif(lakon(nelem)(4:5).eq.'10') then 514 nopes=6 515 nface=4 516 nope=10 517 elseif(lakon(nelem)(4:4).eq.'4') then 518 nopes=3 519 nface=4 520 nope=4 521 elseif(lakon(nelem)(4:5).eq.'15') then 522 if(jface.le.2) then 523 nopes=6 524 else 525 nopes=8 526 endif 527 nface=5 528 nope=15 529 elseif(lakon(nelem)(3:4).eq.'D6') then 530 if(jface.le.2) then 531 nopes=3 532 else 533 nopes=4 534 endif 535 nface=5 536 nope=6 537 else 538 cycle 539 endif 540! 541! determining the nodes of the face 542! 543 if(nface.eq.4) then 544 do i=1,nopes 545 nodef(i)=kon(indexe+ifacet(i,jface)) 546 enddo 547 elseif(nface.eq.5) then 548 if(nope.eq.6) then 549 do i=1,nopes 550 nodef(i)=kon(indexe+ifacew1(i,jface)) 551 enddo 552 elseif(nope.eq.15) then 553 do i=1,nopes 554 nodef(i)=kon(indexe+ifacew2(i,jface)) 555 enddo 556 endif 557 elseif(nface.eq.6) then 558 do i=1,nopes 559 nodef(i)=kon(indexe+ifaceq(i,jface)) 560 enddo 561 endif 562! 563! loop over the nodes belonging to the face 564! ics(1,*): pretension node 565! dcs(*): area corresponding to pretension node 566! 567 do i=1,nopes 568 node=nodef(i) 569 call nident2(ics,node,npt,id) 570 if(id.gt.0) then 571 if(ics(1,id).eq.node) then 572 cycle 573 endif 574 endif 575! 576! updating ics 577! 578 npt=npt+1 579 do j=npt,id+2,-1 580 ics(1,j)=ics(1,j-1) 581 dcs(j)=dcs(j-1) 582 enddo 583 ics(1,id+1)=node 584 dcs(id+1)=0.d0 585 enddo 586! 587! calculating the area of the face and its contributions 588! to the facial nodes 589! 590! number of integration points 591! 592 if(lakon(nelem)(3:5).eq.'D8R') then 593 mint=1 594 elseif(lakon(nelem)(3:4).eq.'D8') then 595 mint=4 596 elseif(lakon(nelem)(4:6).eq.'20R') then 597 mint=4 598 elseif(lakon(nelem)(4:4).eq.'2') then 599 mint=9 600 elseif(lakon(nelem)(4:5).eq.'10') then 601 mint=3 602 elseif(lakon(nelem)(4:4).eq.'4') then 603 mint=1 604 elseif(lakon(nelem)(3:4).eq.'D6') then 605 mint=1 606 elseif(lakon(nelem)(4:5).eq.'15') then 607 if(jface.le.2) then 608 mint=3 609 else 610 mint=4 611 endif 612 endif 613! 614 do i=1,nopes 615 areanodal(i)=0.d0 616 do j=1,3 617 xl2(j,i)=co(j,nodef(i)) 618 enddo 619 enddo 620! 621 do m=1,mint 622 if((lakon(nelem)(3:5).eq.'D8R').or. 623 & ((lakon(nelem)(3:4).eq.'D6').and.(nopes.eq.4))) then 624 xi=gauss2d1(1,m) 625 et=gauss2d1(2,m) 626 weight=weight2d1(m) 627 elseif((lakon(nelem)(3:4).eq.'D8').or. 628 & (lakon(nelem)(4:6).eq.'20R').or. 629 & ((lakon(nelem)(4:5).eq.'15').and. 630 & (nopes.eq.8))) then 631 xi=gauss2d2(1,m) 632 et=gauss2d2(2,m) 633 weight=weight2d2(m) 634 elseif(lakon(nelem)(4:4).eq.'2') then 635 xi=gauss2d3(1,m) 636 et=gauss2d3(2,m) 637 weight=weight2d3(m) 638 elseif((lakon(nelem)(4:5).eq.'10').or. 639 & ((lakon(nelem)(4:5).eq.'15').and. 640 & (nopes.eq.6))) then 641 xi=gauss2d5(1,m) 642 et=gauss2d5(2,m) 643 weight=weight2d5(m) 644 elseif((lakon(nelem)(4:4).eq.'4').or. 645 & ((lakon(nelem)(3:4).eq.'D6').and. 646 & (nopes.eq.3))) then 647 xi=gauss2d4(1,m) 648 et=gauss2d4(2,m) 649 weight=weight2d4(m) 650 endif 651! 652 if(nopes.eq.8) then 653 call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag) 654 elseif(nopes.eq.4) then 655 call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag) 656 elseif(nopes.eq.6) then 657 call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 658 elseif(nopes.eq.3) then 659 call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 660 endif 661! 662! calculating the total area and nodal area 663! 664 xsj=weight*dsqrt(xsj2(1)**2+xsj2(2)**2+xsj2(3)**2) 665 area=area+xsj 666 do i=1,nopes 667 areanodal(i)=areanodal(i)+xsj*shp2(4,i) 668 enddo 669! 670 enddo 671! 672! inserting the nodal area into field dcs 673! 674 do i=1,nopes 675 node=nodef(i) 676 call nident2(ics,node,npt,id) 677 dcs(id)=dcs(id)+areanodal(i) 678 enddo 679! 680 enddo 681! 682! create for each node in the surface a new node. 683! the ratio of the displacements between both nodes 684! is governed by the area for which each node is 685! representative 686! 687! initializing the location of the center of gravity 688! 689 do k=1,3 690 cgx(k)=0.d0 691 enddo 692! 693 do m=1,npt 694 nk=nk+1 695 if(nk.gt.nk_) then 696 write(*,*) 697 & '*ERROR reading *DISTRIBUTING: increase nk_' 698 ier=1 699 return 700 endif 701 node=ics(1,m) 702 ics(1,m)=nk 703 do k=1,3 704 co(k,nk)=co(k,node) 705 enddo 706! 707! generate the connecting MPC's 708! 709 do k=1,3 710! 711! contribution to the location of the center of gravity 712! 713 cgx(k)=cgx(k)+dcs(m)*co(k,node) 714c cgx(k)=cgx(k)+co(k,node) 715! 716 nmpc=nmpc+1 717 if(nmpc.gt.nmpc_) then 718 write(*,*) 719 & '*ERROR reading *DISTRIBUTING: increase nmpc_' 720 ier=1 721 return 722 endif 723! 724! MPC: u(old node)= 725! (mean area)/(area_m) * u(new node) (in all directions) 726! 727! check whether MPC was already used 728! 729 ipompc(nmpc)=mpcfree 730 labmpc(nmpc)=' ' 731 idof=8*(node-1)+k 732 call nident(ikmpc,idof,nmpc-1,id) 733 if(id.gt.0) then 734 if(ikmpc(id).eq.idof) then 735 write(*,*) '*ERROR reading *COUPLING: dof',k 736 write(*,*) ' in node ',node 737 write(*,*) ' was already used' 738 ier=1 739 return 740 endif 741 endif 742 do l=nmpc,id+2,-1 743 ikmpc(l)=ikmpc(l-1) 744 ilmpc(l)=ilmpc(l-1) 745 enddo 746 ikmpc(id+1)=idof 747 ilmpc(id+1)=nmpc 748! 749! generating the terms of the MPC 750! 751 nodempc(1,mpcfree)=node 752 nodempc(2,mpcfree)=k 753 coefmpc(mpcfree)=1.d0 754 mpcfree=nodempc(3,mpcfree) 755! 756 nodempc(1,mpcfree)=nk 757 nodempc(2,mpcfree)=k 758 coefmpc(mpcfree)=-area/(npt*dcs(m)) 759 mpcfreeold=mpcfree 760 mpcfree=nodempc(3,mpcfree) 761 nodempc(3,mpcfreeold)=0 762 enddo 763 enddo 764! 765 do k=1,3 766 cgx(k)=cgx(k)/area 767c cgx(k)=cgx(k)/npt 768 enddo 769! 770! calculating a standard deviation; this quantity will 771! serve as a limit for checking the closeness of individual 772! nodes to the center of gravity 773! 774 stdev=0.d0 775 do m=1,npt 776 node=ics(1,m) 777 do k=1,3 778c stdev=stdev+(dcs(m)*co(k,node)-cgx(k))**2 779 stdev=stdev+(co(k,node)-cgx(k))**2 780 enddo 781 enddo 782 stdev=stdev/npt 783! 784! generating the translational MPC's (default) 785! 786 do m=1,3 787 node=ics(1,1) 788 idof=8*(node-1)+m 789 call nident(ikmpc,idof,nmpc,id) 790! 791 nmpc=nmpc+1 792 if(nmpc.gt.nmpc_) then 793 write(*,*) '*ERROR reading *COUPLING: increase nmpc_' 794 ier=1 795 return 796 endif 797 labmpc(nmpc)=' ' 798 ipompc(nmpc)=mpcfree 799! 800! updating ikmpc and ilmpc 801! 802 do j=nmpc,id+2,-1 803 ikmpc(j)=ikmpc(j-1) 804 ilmpc(j)=ilmpc(j-1) 805 enddo 806 ikmpc(id+1)=idof 807 ilmpc(id+1)=nmpc 808! 809 do j=1,npt 810 node=ics(1,j) 811 nodempc(1,mpcfree)=node 812 nodempc(2,mpcfree)=m 813 coefmpc(mpcfree)=1.d0 814 mpcfree=nodempc(3,mpcfree) 815 enddo 816 nodempc(1,mpcfree)=irefnode 817 nodempc(2,mpcfree)=m 818 coefmpc(mpcfree)=-npt 819 mpcfreeold=mpcfree 820 mpcfree=nodempc(3,mpcfree) 821 nodempc(3,mpcfreeold)=0 822 enddo 823! 824! generating the rotational MPC's 825! 826 irotnode(1)=0 827! 828! reading the degrees of freedom 829! 830 do 831 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 832 & ipoinp,inp,ipoinpc) 833 if((istat.lt.0).or.(key.eq.1)) return 834! 835 read(textpart(1)(1:10),'(i10)',iostat=istat) ibounstart 836 if(istat.gt.0) then 837 call inputerror(inpc,ipoinpc,iline, 838 & "*KINEMATIC%",ier) 839 return 840 endif 841 if(ibounstart.lt.1) then 842 write(*,*) '*ERROR reading *KINEMATIC' 843 write(*,*) ' starting degree of freedom cannot' 844 write(*,*) ' be less than 1' 845 write(*,*) ' ' 846 call inputerror(inpc,ipoinpc,iline, 847 & "*KINEMATIC%",ier) 848 return 849 endif 850! 851 if(textpart(2)(1:1).eq.' ') then 852 ibounend=ibounstart 853 else 854 read(textpart(2)(1:10),'(i10)',iostat=istat) ibounend 855 if(istat.gt.0) then 856 call inputerror(inpc,ipoinpc,iline, 857 & "*BOUNDARY%",ier) 858 return 859 endif 860 endif 861! 862 ibounstart=max(4,ibounstart) 863 if(ibounend.gt.6) then 864 write(*,*) '*ERROR reading *DISTRIBUTING' 865 write(*,*) ' final degree of freedom cannot' 866 write(*,*) ' exceed 6' 867 write(*,*) ' ' 868 call inputerror(inpc,ipoinpc,iline, 869 & "*DISTRIBUTING%",ier) 870 return 871 elseif(ibounend.lt.ibounstart) then 872 cycle 873 endif 874! 875 if((ibounend.gt.3).and.(irotnode(1).eq.0)) then 876! 877! check the orientation definition: the local reference 878! system is not allowed to be cylindrical 879! 880 if(iorientation.ne.0) then 881 if(orab(7,iorientation).lt.0.d0) then 882 write(*,*) '*ERROR reading *DISTRIBUTING' 883 write(*,*) ' a cylindrical local coordinate' 884 write(*,*) ' system is not allowed' 885 ier=1 886 return 887 endif 888! 889 call transformatrix(orab(1,iorientation), 890 & co(1,irefnode),a) 891 endif 892! 893! generating a rotational reference node 894! 895 do k=1,3 896 nk=nk+1 897 if(nk.gt.nk_) then 898 write(*,*) 899 & '*ERROR reading *DISTRIBUTING: increase nk_' 900 ier=1 901 return 902 endif 903 irotnode(k)=nk 904 do l=1,3 905 co(l,nk)=co(l,irefnode) 906 enddo 907 enddo 908! 909! generating connecting MPCs between the rotational 910! dofs of irefnode and the translational dofs of 911! irotnode 912! 913 do k=1,3 914! 915 nmpc=nmpc+1 916 if(nmpc.gt.nmpc_) then 917 write(*,*) 918 & '*ERROR reading *DISTRIBUTING: increase nmpc_' 919 ier=1 920 return 921 endif 922! 923! the internal dofs for rotation are 4, 5 and 7 924! 925 ipompc(nmpc)=mpcfree 926 labmpc(nmpc)='ROTTRACOUPLING ' 927 idof=8*(irefnode-1)+k+3 928 call nident(ikmpc,idof,nmpc-1,id) 929 do l=nmpc,id+2,-1 930 ikmpc(l)=ikmpc(l-1) 931 ilmpc(l)=ilmpc(l-1) 932 enddo 933 ikmpc(id+1)=idof 934 ilmpc(id+1)=nmpc 935! 936 nodempc(1,mpcfree)=irefnode 937 nodempc(2,mpcfree)=k+4 938 coefmpc(mpcfree)=1.d0 939 mpcfree=nodempc(3,mpcfree) 940! 941 nodempc(1,mpcfree)=irotnode(k) 942 nodempc(2,mpcfree)=1 943 coefmpc(mpcfree)=-1.d0 944 mpcfreeold=mpcfree 945 mpcfree=nodempc(3,mpcfree) 946 nodempc(3,mpcfreeold)=0 947 enddo 948! 949! generating a node for the inhomogeneous term 950! 951 nk=nk+1 952 if(nk.gt.nk_) then 953 write(*,*) 954 & '*ERROR reading *DISTRIBUTING: increase nk_' 955 ier=1 956 return 957 endif 958 inhomnode=nk 959! 960! fictituous center of gravity (for compatibility 961! reasons with usermpc.f) 962! 963c do k=1,3 964c cgx(k)=co(k,irefnode) 965c enddo 966 endif 967! 968! generating rotational MPC's 969! 970 do kk=ibounstart,ibounend 971 k=kk-3 972! 973! axis of rotation 974! 975 if(iorientation.eq.0) then 976 do j=1,3 977 aa(j)=0.d0 978 enddo 979 aa(k)=1.d0 980 else 981 do j=1,3 982 aa(j)=a(j,k) 983 enddo 984 endif 985c write(*,*) 'couplings aa ',(aa(j),j=1,3) 986! 987! storing the axis of rotation as coordinates of the 988! rotation nodes 989! 990 do j=1,3 991 co(j,irotnode(k))=aa(j) 992 enddo 993! 994 nmpc=nmpc+1 995 if(nmpc.gt.nmpc_) then 996 write(*,*) 997 & '*ERROR reading *DISTRIBUTING: increase nmpc_' 998 ier=1 999 return 1000 endif 1001! 1002 ipompc(nmpc)=mpcfree 1003c labmpc(nmpc)='MEANROTDISTRIB ' 1004 labmpc(nmpc)=' ' 1005! 1006! defining the terms of the MPC 1007! 1008c write(*,*) 'couplings, npt ',npt 1009 do m=1,npt 1010 do j=1,3 1011 nodempc(1,mpcfree)=ics(1,m) 1012 nodempc(2,mpcfree)=0 1013 coefmpc(mpcfree)=0.d0 1014 mpcfree=nodempc(3,mpcfree) 1015 enddo 1016c write(*,*) 'couplings, co',m,(co(j,ics(1,m)),j=1,3) 1017 enddo 1018 nodempc(1,mpcfree)=irotnode(k) 1019 nodempc(2,mpcfree)=1 1020 coefmpc(mpcfree)=0.d0 1021 mpcfree=nodempc(3,mpcfree) 1022 nodempc(1,mpcfree)=inhomnode 1023 nodempc(2,mpcfree)=k 1024c changed on 2. november 2018 1025 coefmpc(mpcfree)=1.d0 1026c coefmpc(mpcfree)=0.d0 1027 mpcfreeold=mpcfree 1028 mpcfree=nodempc(3,mpcfree) 1029 nodempc(3,mpcfreeold)=0 1030! 1031! storing the inhomogeneous value as a boundary 1032! condition in the inhomogeneous node 1033! 1034 idof=8*(inhomnode-1)+k 1035 call nident(ikboun,idof,nboun,id) 1036 nboun=nboun+1 1037 if(nboun.gt.nboun_) then 1038 write(*,*) '*ERROR reading *COUPLING: increase nboun_' 1039 ier=1 1040 return 1041 endif 1042 nodeboun(nboun)=inhomnode 1043 ndirboun(nboun)=k 1044 xboun(nboun)=0.d0 1045 typeboun(nboun)='U' 1046 do l=nboun,id+2,-1 1047 ikboun(l)=ikboun(l-1) 1048 ilboun(l)=ilboun(l-1) 1049 enddo 1050 ikboun(id+1)=idof 1051 ilboun(id+1)=nboun 1052c write(*,*) 'couplings cgx',(cgx(l),l=1,3) 1053c write(*,*) 'couplings stdev',stdev 1054! 1055! calculating the coefficients of the rotational MPC 1056! 1057! loop over all nodes belonging to the mean rotation MPC 1058! 1059 index2=ipompc(nmpc) 1060 do 1061 node=nodempc(1,index2) 1062 if(node.eq.irotnode(k)) exit 1063! 1064! relative positions 1065! 1066 do j=1,3 1067 pi(j)=co(j,node)-cgx(j) 1068 enddo 1069c write(*,*) 'couplings pi',(pi(j),j=1,3) 1070! 1071! projection on a plane orthogonal to the rotation vector 1072! 1073 c1=pi(1)*aa(1)+pi(2)*aa(2)+pi(3)*aa(3) 1074 do j=1,3 1075 pi(j)=pi(j)-c1*aa(j) 1076 enddo 1077c write(*,*) 'couplings c1 pi',c1,(pi(j),j=1,3) 1078! 1079 c1=pi(1)*pi(1)+pi(2)*pi(2)+pi(3)*pi(3) 1080c if(c1.lt.1.d-20) then 1081 if(c1.lt.stdev*1.d-10) then 1082 write(*,*) 1083 & '*WARNING reading *DISTRIBUTING: node ',node 1084 write(*,*) ' is very close to the ' 1085 write(*,*) ' rotation axis through the' 1086 write(*,*) ' center of gravity of' 1087 write(*,*) ' the nodal cloud in a' 1088 write(*,*) ' mean rotation MPC.' 1089 write(*,*) ' This node is not taken' 1090 write(*,*) ' into account in the MPC' 1091 index2=nodempc(3,nodempc(3,nodempc(3,index2))) 1092 cycle 1093 endif 1094! 1095 do j=1,3 1096 if(j.eq.1) then 1097 c4=aa(2)*pi(3)-aa(3)*pi(2) 1098 elseif(j.eq.2) then 1099 c4=aa(3)*pi(1)-aa(1)*pi(3) 1100 else 1101 c4=aa(1)*pi(2)-aa(2)*pi(1) 1102 endif 1103 c9=c4/c1 1104c write(*,*) 'couplings c4,c9',j,c4,c9 1105! 1106 index1=ipompc(nmpc) 1107 do 1108 node1=nodempc(1,index1) 1109 if(node1.eq.irotnode(k)) exit 1110 if(node1.eq.node) then 1111 c10=c9*(1.d0-1.d0/real(npt)) 1112 else 1113 c10=-c9/real(npt) 1114 endif 1115 if(j.eq.1) then 1116 coefmpc(index1)=coefmpc(index1)+c10 1117c write(*,*) 'couplings c10', 1118c & j,c10,coefmpc(index1) 1119 elseif(j.eq.2) then 1120 coefmpc(nodempc(3,index1))= 1121 & coefmpc(nodempc(3,index1))+c10 1122c write(*,*) 'couplings c10', 1123c & j,c10,coefmpc(nodempc(3,index1)) 1124 else 1125 coefmpc(nodempc(3,nodempc(3,index1)))= 1126 & coefmpc(nodempc(3,nodempc(3,index1)))+c10 1127c write(*,*) 'couplings c10', 1128c & j,c10,coefmpc(nodempc(3,nodempc(3,index1))) 1129 endif 1130 index1= 1131 & nodempc(3,nodempc(3,nodempc(3,index1))) 1132 enddo 1133 enddo 1134 index2=nodempc(3,nodempc(3,nodempc(3,index2))) 1135 enddo 1136 coefmpc(index2)=-npt 1137! 1138! assigning the degrees of freedom 1139! 1140 j=0 1141 index2=ipompc(nmpc) 1142 do 1143 j=j+1 1144 if(j.gt.3) j=1 1145 nodempc(2,index2)=j 1146c write(*,*) 'couplings a',(coefmpc(index2)) 1147 index2=nodempc(3,index2) 1148 if(nodempc(1,index2).eq.irotnode(k)) exit 1149 enddo 1150! 1151! look for biggest coefficient of all but the last 1152! regular node. The general form of the MPC is: 1153! regular nodes, rotationl dof and inhomogeneous dof 1154! 1155 indexmax=0 1156 index2=ipompc(nmpc) 1157 amax=1.d-5 1158! 1159! coefficients of the first terms in the translational 1160! dofs 1161! 1162 coef(1)=coefmpc(index2) 1163 coef(2)=coefmpc(nodempc(3,index2)) 1164 coef(3)=coefmpc(nodempc(3,nodempc(3,index2))) 1165! 1166! loop over all regular nodes 1167! 1168 do i=1,3*npt 1169 if(dabs(coefmpc(index2)).gt.amax) then 1170 idir=nodempc(2,index2) 1171! 1172! in the translational mpcs the coefficients of 1173! all npt terms are equal to 1 1174! in the rotational mpcs the coefficient of the 1175! dependent term should be different from the 1176! coefficient of the node corresponding to the 1177! translational dependent term with the same dof 1178! 1179c if(dabs(coefmpc(index2)-coef(idir)).lt.1.d-10) then 1180 if(dabs(coefmpc(index2)-coef(idir)).lt. 1181 & dabs(coefmpc(index2))/1000.)then 1182 index2=nodempc(3,index2) 1183 cycle 1184 endif 1185! 1186 node=nodempc(1,index2) 1187 idof=8*(node-1)+idir 1188! 1189! check whether the node was used in another MPC 1190! 1191 call nident(ikmpc,idof,nmpc-1,id) 1192 if(id.gt.0) then 1193 if(ikmpc(id).eq.idof) then 1194 index2=nodempc(3,index2) 1195 cycle 1196 endif 1197 endif 1198! 1199! check whether the node was used in a SPC 1200! 1201 call nident(ikboun,idof,nboun,id) 1202 if(id.gt.0) then 1203 if(ikboun(id).eq.idof) then 1204 index2=nodempc(3,index2) 1205 cycle 1206 endif 1207 endif 1208! 1209 amax=dabs(coefmpc(index2)) 1210 indexmax=index2 1211 endif 1212! 1213! after each node has been treated, a check is performed 1214! whether indexmax is already nonzero. Taking too many 1215! nodes into account may lead to dependent equations if 1216! all rotational dofs are constrained. 1217! 1218 if((i/3)*3.eq.i) then 1219 if(indexmax.ne.0) exit 1220 endif 1221! 1222 index2=nodempc(3,index2) 1223 enddo 1224! 1225 if(indexmax.eq.0) then 1226 write(*,*) 1227 & '*WARNING reading *DISTRIBUTING: no MPC is ' 1228 write(*,*) ' generated for the mean' 1229 write(*,*) ' rotation in node ',irefnode 1230 write(*,*) ' and direction ',kk 1231! 1232! removing the MPC 1233! 1234 mpcfreeold=mpcfree 1235 index2=ipompc(nmpc) 1236 ipompc(nmpc)=0 1237 mpcfree=index2 1238! 1239! removing the MPC from fields nodempc and coefmpc 1240! 1241 do 1242 nodempc(1,index2)=0 1243 nodempc(2,index2)=0 1244 coefmpc(index2)=0 1245 if(nodempc(3,index2).ne.0) then 1246 index2=nodempc(3,index2) 1247 else 1248 nodempc(3,index2)=mpcfreeold 1249 exit 1250 endif 1251 enddo 1252! 1253! decrementing nmpc 1254! 1255 nmpc=nmpc-1 1256! 1257 cycle 1258 endif 1259! 1260 nodemax=nodempc(1,indexmax) 1261 idirmax=nodempc(2,indexmax) 1262 index2=ipompc(nmpc) 1263! 1264 nodeold=nodempc(1,index2) 1265 idirold=nodempc(2,index2) 1266! 1267! exchange the node information in the MPC 1268! 1269 if(nodemax.ne.nodeold) then 1270! 1271 indexold(1)=index2 1272 indexold(2)=nodempc(3,index2) 1273 indexold(3)=nodempc(3,indexold(2)) 1274! 1275 index2=nodempc(3,indexold(3)) 1276 do 1277 if(nodempc(1,index2).eq.irotnode(k)) exit 1278 if(nodempc(1,index2).eq.nodemax) then 1279 if(nodempc(2,index2).eq.1) then 1280 indexnew(1)=index2 1281 elseif(nodempc(2,index2).eq.2) then 1282 indexnew(2)=index2 1283 else 1284 indexnew(3)=index2 1285 endif 1286 endif 1287 index2=nodempc(3,index2) 1288 enddo 1289! 1290 do j=1,3 1291 node=nodempc(1,indexold(j)) 1292 idir=nodempc(2,indexold(j)) 1293 xcoef=coefmpc(indexold(j)) 1294 nodempc(1,indexold(j))=nodempc(1,indexnew(j)) 1295 nodempc(2,indexold(j))=nodempc(2,indexnew(j)) 1296 coefmpc(indexold(j))=coefmpc(indexnew(j)) 1297 nodempc(1,indexnew(j))=node 1298 nodempc(2,indexnew(j))=idir 1299 coefmpc(indexnew(j))=xcoef 1300 enddo 1301 endif 1302! 1303! exchange the direction information in the MPC 1304! 1305 index2=ipompc(nmpc) 1306 if(idirmax.ne.1) then 1307 indexold(1)=index2 1308 if(idirmax.eq.2) then 1309 indexnew(1)=nodempc(3,index2) 1310 else 1311 indexnew(1)=nodempc(3,nodempc(3,index2)) 1312 endif 1313! 1314 do j=1,1 1315 idir=nodempc(2,indexold(j)) 1316 xcoef=coefmpc(indexold(j)) 1317 nodempc(2,indexold(j))=nodempc(2,indexnew(j)) 1318 coefmpc(indexold(j))=coefmpc(indexnew(j)) 1319 nodempc(2,indexnew(j))=idir 1320 coefmpc(indexnew(j))=xcoef 1321 enddo 1322 endif 1323! 1324! determining the dependent dof of the MPC 1325! 1326 idof=8*(nodemax-1)+idirmax 1327 call nident(ikmpc,idof,nmpc-1,id) 1328 do l=nmpc,id+2,-1 1329 ikmpc(l)=ikmpc(l-1) 1330 ilmpc(l)=ilmpc(l-1) 1331 enddo 1332 ikmpc(id+1)=idof 1333 ilmpc(id+1)=nmpc 1334! 1335 enddo 1336 enddo 1337! 1338! generating the translational MPC's (default) 1339! 1340c do m=1,3 1341c node=ics(1,1) 1342c idof=8*(node-1)+m 1343c call nident(ikmpc,idof,nmpc,id) 1344c! 1345c nmpc=nmpc+1 1346c if(nmpc.gt.nmpc_) then 1347c write(*,*) '*ERROR reading *COUPLING: increase nmpc_' 1348c ier=1 1349c return 1350c endif 1351c labmpc(nmpc)=' ' 1352c ipompc(nmpc)=mpcfree 1353c! 1354c! updating ikmpc and ilmpc 1355c! 1356c do j=nmpc,id+2,-1 1357c ikmpc(j)=ikmpc(j-1) 1358c ilmpc(j)=ilmpc(j-1) 1359c enddo 1360c ikmpc(id+1)=idof 1361c ilmpc(id+1)=nmpc 1362c! 1363c do j=1,npt 1364c node=ics(1,j) 1365c nodempc(1,mpcfree)=node 1366c nodempc(2,mpcfree)=m 1367c coefmpc(mpcfree)=1.d0 1368c mpcfree=nodempc(3,mpcfree) 1369c enddo 1370c nodempc(1,mpcfree)=irefnode 1371c nodempc(2,mpcfree)=m 1372c coefmpc(mpcfree)=-npt 1373c mpcfreeold=mpcfree 1374c mpcfree=nodempc(3,mpcfree) 1375c nodempc(3,mpcfreeold)=0 1376c enddo 1377 else 1378 write(*,*) 1379 & '*ERROR reading *COUPLING: the line following' 1380 write(*,*) ' *COUPLING must contain the' 1381 write(*,*) ' *KINEMATIC keyword or the' 1382 write(*,*) ' *DISTRIBUTING keyword' 1383 write(*,*) ' ' 1384 call inputerror(inpc,ipoinpc,iline, 1385 & "*COUPLING%",ier) 1386 return 1387 endif 1388! 1389 return 1390 end 1391 1392