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 gentiedmpc(tieset,ntie,itietri,ipkon,kon, 20 & lakon,set,istartset,iendset,ialset,cg,straight, 21 & koncont,co,xo,yo,zo,x,y,z,nx,ny,nz,nset, 22 & ifaceslave,istartfield,iendfield,ifield, 23 & ipompc,nodempc,coefmpc,nmpc,nmpctied,mpcfree,ikmpc,ilmpc, 24 & labmpc,ithermal,tietol,nef,ncont,imastop,ikboun,nboun,kind) 25! 26! generates MPC's for the slave tied contact nodes 27! 28 implicit none 29! 30 character*1 kind 31 character*8 lakon(*) 32 character*20 labmpc(*) 33 character*81 tieset(3,*),slavset,set(*) 34! 35 integer ntie,nset,istartset(*),iendset(*),ialset(*),iwrite, 36 & itietri(2,ntie),ipkon(*),kon(*),koncont(4,*),node, 37 & neigh(1),iflag,kneigh,i,j,k,l,isol,itri,ll,kflag,n,nx(*), 38 & ny(*),nz(*),nstart,ifaceq(8,6),ifacet(6,4),nboun, 39 & ifacew1(4,5),ifacew2(8,5),nelem,jface,indexe,imastop(3,*), 40 & nnodelem,nface,nope,nodef(8),idof,kstart,kend,jstart,id, 41 & jend,ifield(*),istartfield(*),iendfield(*),ifaceslave(*), 42 & ipompc(*),nodempc(3,*),nmpc,nmpctied,mpcfree,ikmpc(*), 43 & ilmpc(*),ithermal(*),nef,ncont,mpcfreeold,m,id1,ikboun(*), 44 & itriold,itrinew,ntriangle,ntriangle_,itriangle(100) 45! 46 real*8 cg(3,*),straight(16,*),co(3,*),p(3), 47 & dist,xo(*),yo(*),zo(*),x(*),y(*),z(*),pl(3,9), 48 & ratio(9),xi,et,coefmpc(*),tietol(3,*),tolloc 49! 50! nodes per face for hex elements 51! 52 data ifaceq /4,3,2,1,11,10,9,12, 53 & 5,6,7,8,13,14,15,16, 54 & 1,2,6,5,9,18,13,17, 55 & 2,3,7,6,10,19,14,18, 56 & 3,4,8,7,11,20,15,19, 57 & 4,1,5,8,12,17,16,20/ 58! 59! nodes per face for tet elements 60! 61 data ifacet /1,3,2,7,6,5, 62 & 1,2,4,5,9,8, 63 & 2,3,4,6,10,9, 64 & 1,4,3,8,10,7/ 65! 66! nodes per face for linear wedge elements 67! 68 data ifacew1 /1,3,2,0, 69 & 4,5,6,0, 70 & 1,2,5,4, 71 & 2,3,6,5, 72 & 3,1,4,6/ 73! 74! nodes per face for quadratic wedge elements 75! 76 data ifacew2 /1,3,2,9,8,7,0,0, 77 & 4,5,6,10,11,12,0,0, 78 & 1,2,5,4,7,14,10,13, 79 & 2,3,6,5,8,15,11,14, 80 & 3,1,4,6,9,13,12,15/ 81! 82! opening a file to store the nodes which are not connected 83! 84 iwrite=0 85 open(40,file='WarnNodeMissTiedContact.nam',status='unknown') 86 write(40,*) '*NSET,NSET=WarnNodeMissTiedContact' 87! 88 nmpctied=nmpc 89! 90! calculating a typical element size 91! 92 tolloc=0.d0 93 do i=1,ncont 94 tolloc=tolloc+dabs(straight(1,i)*cg(1,i)+ 95 & straight(2,i)*cg(2,i)+ 96 & straight(3,i)*cg(3,i)+ 97 & straight(4,i)) 98 enddo 99 tolloc=0.025d0*tolloc/ncont 100! 101! maximum number of neighboring master triangles for a slave node 102! 103 kflag=2 104! 105 do i=1,ntie 106 if(tieset(1,i)(81:81).ne.kind) cycle 107! 108! determining for which dofs MPC's have to be generated 109! 110 if(kind.eq.'T') then 111! 112! thermomechanical ties or CFD 113! 114 if(nef.gt.0) then 115 if(ithermal(2).le.1) then 116 kstart=1 117 kend=4 118 else 119 kstart=0 120 kend=4 121 endif 122 else 123 if(ithermal(2).le.1) then 124 kstart=1 125 kend=3 126 elseif(ithermal(2).eq.2) then 127 kstart=0 128 kend=0 129 else 130 kstart=0 131 kend=3 132 endif 133 endif 134 elseif(kind.eq.'E') then 135! 136! electromagnetic ties between the domains 137! 138 if((tieset(1,i)(1:2).eq.'12').or. 139 & (tieset(1,i)(1:2).eq.'13').or. 140 & (tieset(1,i)(1:2).eq.'23')) then 141 kstart=1 142 kend=3 143 elseif((tieset(1,i)(1:2).eq.'21').or. 144 & (tieset(1,i)(1:2).eq.'31')) then 145 kstart=5 146 kend=5 147 endif 148 endif 149! 150 iflag=0 151 kneigh=1 152 slavset=tieset(2,i) 153! 154! default tolerance if none is specified 155! 156 if(tietol(1,i).lt.1.d-10) tietol(1,i)=tolloc 157! 158! determining the slave set 159! 160 if(ifaceslave(i).eq.0) then 161c do j=1,nset 162c if(set(j).eq.slavset) then 163c exit 164c endif 165c enddo 166 call cident81(set,slavset,nset,id) 167 j=nset+1 168 if(id.gt.0) then 169 if(slavset.eq.set(id)) then 170 j=id 171 endif 172 endif 173 jstart=istartset(j) 174 jend=iendset(j) 175 else 176 jstart=istartfield(i) 177 jend=iendfield(i) 178 endif 179! 180 nstart=itietri(1,i)-1 181 n=itietri(2,i)-nstart 182 if(n.lt.kneigh) kneigh=n 183 do j=1,n 184 xo(j)=cg(1,nstart+j) 185 x(j)=xo(j) 186 nx(j)=j 187 yo(j)=cg(2,nstart+j) 188 y(j)=yo(j) 189 ny(j)=j 190 zo(j)=cg(3,nstart+j) 191 z(j)=zo(j) 192 nz(j)=j 193 enddo 194 call dsort(x,nx,n,kflag) 195 call dsort(y,ny,n,kflag) 196 call dsort(z,nz,n,kflag) 197! 198 do j=jstart,jend 199 if(((ifaceslave(i).eq.0).and.(ialset(j).gt.0)).or. 200 & (ifaceslave(i).eq.1)) then 201! 202 if(ifaceslave(i).eq.0) then 203 node=ialset(j) 204 else 205 node=ifield(j) 206 endif 207! 208 do k=1,3 209 p(k)=co(k,node) 210 enddo 211! 212! determining the kneigh neighboring master contact 213! triangle centers of gravity 214! 215 call near3d(xo,yo,zo,x,y,z,nx,ny,nz,p(1),p(2),p(3), 216 & n,neigh,kneigh) 217! 218 isol=0 219! 220 itriold=0 221 itri=neigh(1)+itietri(1,i)-1 222 ntriangle=0 223 ntriangle_=100 224! 225 loop1: do 226 do l=1,3 227 ll=4*l-3 228 dist=straight(ll,itri)*p(1)+ 229 & straight(ll+1,itri)*p(2)+ 230 & straight(ll+2,itri)*p(3)+ 231 & straight(ll+3,itri) 232c if(dist.gt.1.d-6) then 233c if(dist.gt.tietol(1,i)) then 234 if(dist.gt.tolloc) then 235 itrinew=imastop(l,itri) 236 if(itrinew.eq.0) then 237c write(*,*) '**border reached' 238 isol=-1 239 exit loop1 240 elseif(itrinew.eq.itriold) then 241c write(*,*) '**solution in between triangles' 242 isol=itri 243 exit loop1 244 else 245 call nident(itriangle,itrinew,ntriangle,id) 246 if(id.gt.0) then 247 if(itriangle(id).eq.itrinew) then 248c write(*,*) '**circular path; no solution' 249 isol=-2 250 exit loop1 251 endif 252 endif 253 ntriangle=ntriangle+1 254 if(ntriangle.gt.ntriangle_) then 255c write(*,*) '**too many iterations' 256 isol=-3 257 exit loop1 258 endif 259 do k=ntriangle,id+2,-1 260 itriangle(k)=itriangle(k-1) 261 enddo 262 itriangle(id+1)=itrinew 263 itriold=itri 264 itri=itrinew 265 cycle loop1 266 endif 267 elseif(l.eq.3) then 268c write(*,*) '**regular solution' 269 isol=itri 270 exit loop1 271 endif 272 enddo 273 enddo loop1 274! 275! if an opposite triangle is found: check the distance 276! perpendicular to the triangle 277! 278 if(isol.gt.0) then 279 dist=dabs(straight(13,itri)*p(1)+ 280 & straight(14,itri)*p(2)+ 281 & straight(15,itri)*p(3)+ 282 & straight(16,itri)) 283 if(dist.gt.tietol(1,i)) isol=0 284 endif 285! 286 if(isol.le.0) then 287! 288! no MPC is generated 289! 290 write(*,*) '*WARNING in gentiedmpc: no tied MPC' 291 write(*,*) ' generated for node ',node 292 if(isol.eq.0) then 293 write(*,*) 294 & ' opposite master face found but' 295 write(*,*) ' too far away; distance: ',dist 296 write(*,*) ' tolerance: ',tietol(1,i) 297 else 298 write(*,*) ' no opposite master face' 299 write(*,*) ' found; in-face tolerance: ', 300 & tolloc 301 endif 302 write(40,*) node 303 iwrite=1 304 else 305! 306 nelem=int(koncont(4,itri)/10.d0) 307 jface=koncont(4,itri)-10*nelem 308! 309 indexe=ipkon(nelem) 310 if(lakon(nelem)(4:4).eq.'2') then 311 nnodelem=8 312 nface=6 313 elseif(lakon(nelem)(4:4).eq.'8') then 314 nnodelem=4 315 nface=6 316 elseif(lakon(nelem)(4:5).eq.'10') then 317 nnodelem=6 318 nface=4 319 elseif(lakon(nelem)(4:4).eq.'4') then 320 nnodelem=3 321 nface=4 322 elseif(lakon(nelem)(4:5).eq.'15') then 323 if(jface.le.2) then 324 nnodelem=6 325 else 326 nnodelem=8 327 endif 328 nface=5 329 nope=15 330 elseif(lakon(nelem)(4:4).eq.'6') then 331 if(jface.le.2) then 332 nnodelem=3 333 else 334 nnodelem=4 335 endif 336 nface=5 337 nope=6 338 else 339 cycle 340 endif 341! 342! determining the nodes of the face 343! 344 if(nface.eq.4) then 345 do k=1,nnodelem 346 nodef(k)=kon(indexe+ifacet(k,jface)) 347 enddo 348 elseif(nface.eq.5) then 349 if(nope.eq.6) then 350 do k=1,nnodelem 351 nodef(k)=kon(indexe+ifacew1(k,jface)) 352 enddo 353 elseif(nope.eq.15) then 354 do k=1,nnodelem 355 nodef(k)=kon(indexe+ifacew2(k,jface)) 356 enddo 357 endif 358 elseif(nface.eq.6) then 359 do k=1,nnodelem 360 nodef(k)=kon(indexe+ifaceq(k,jface)) 361 enddo 362 endif 363! 364! attaching the node with coordinates in p 365! to the face 366! 367 do k=1,nnodelem 368 do l=1,3 369 pl(l,k)=co(l,nodef(k)) 370 enddo 371 enddo 372 call attach_2d(pl,p,nnodelem,ratio,dist,xi,et) 373! 374! adjusting the coordinates of the node (only 375! if the user did not specify ADJUST=NO on the 376! *TIE card) 377! 378 if(tietol(2,i).gt.0.d0) then 379 do k=1,3 380 co(k,node)=p(k) 381 enddo 382 endif 383! 384! generating MPC's 385! 386 do l=kstart,kend 387 idof=8*(node-1)+l 388 call nident(ikmpc,idof,nmpc,id) 389 if(id.gt.0) then 390 if(ikmpc(id).eq.idof) then 391 write(*,*) '*WARNING in gentiedmpc:' 392 write(*,*) ' DOF ',l,' of node ', 393 & node,' is not active;' 394 write(*,*) ' no tied constraint ', 395 & 'is generated' 396 write(40,*) node 397 iwrite=1 398 cycle 399 endif 400 endif 401! 402 call nident(ikboun,idof,nboun,id1) 403 if(id1.gt.0) then 404 if(ikboun(id1).eq.idof) then 405 write(*,*) '*WARNING in gentiedmpc:' 406 write(*,*) ' DOF ',l,' of node ', 407 & node,' is not active;' 408 write(*,*) ' no tied constraint ', 409 & 'is generated' 410 write(40,*) node 411 iwrite=1 412 cycle 413 endif 414 endif 415! 416 nmpc=nmpc+1 417 labmpc(nmpc)=' ' 418 ipompc(nmpc)=mpcfree 419! 420! updating ikmpc and ilmpc 421! 422 do m=nmpc,id+2,-1 423 ikmpc(m)=ikmpc(m-1) 424 ilmpc(m)=ilmpc(m-1) 425 enddo 426 ikmpc(id+1)=idof 427 ilmpc(id+1)=nmpc 428! 429 nodempc(1,mpcfree)=node 430 nodempc(2,mpcfree)=l 431 coefmpc(mpcfree)=1.d0 432 mpcfree=nodempc(3,mpcfree) 433 if(mpcfree.eq.0) then 434 write(*,*) 435 & '*ERROR in gentiedmpc: increase memmpc_' 436 call exit(201) 437 endif 438 do k=1,nnodelem 439 nodempc(1,mpcfree)=nodef(k) 440 nodempc(2,mpcfree)=l 441 coefmpc(mpcfree)=-ratio(k) 442 mpcfreeold=mpcfree 443 mpcfree=nodempc(3,mpcfree) 444 if(mpcfree.eq.0) then 445 write(*,*) 446 & '*ERROR in gentiedmpc: increase memmpc_' 447 call exit(201) 448 endif 449 enddo 450 nodempc(3,mpcfreeold)=0 451 enddo 452! 453 endif 454! 455 else 456 node=ialset(j-2) 457 do 458 node=node-ialset(j) 459 if(node.ge.ialset(j-1)) exit 460! 461 do k=1,3 462 p(k)=co(k,node) 463 enddo 464! 465! determining the kneigh neighboring master contact 466! triangle centers of gravity 467! 468 call near3d(xo,yo,zo,x,y,z,nx,ny,nz,p(1),p(2),p(3), 469 & n,neigh,kneigh) 470! 471 isol=0 472! 473 itriold=0 474 itri=neigh(1)+itietri(1,i)-1 475 ntriangle=0 476 ntriangle_=100 477! 478 loop2: do 479 do l=1,3 480 ll=4*l-3 481 dist=straight(ll,itri)*p(1)+ 482 & straight(ll+1,itri)*p(2)+ 483 & straight(ll+2,itri)*p(3)+ 484 & straight(ll+3,itri) 485c if(dist.gt.1.d-6) then 486c if(dist.gt.tietol(1,i)) then 487 if(dist.gt.tolloc) then 488 itrinew=imastop(l,itri) 489 if(itrinew.eq.0) then 490c write(*,*) '**border reached' 491 isol=-1 492 exit loop2 493 elseif(itrinew.eq.itriold) then 494c write(*,*) '**solution in between triangles' 495 isol=itri 496 exit loop2 497 else 498 call nident(itriangle,itrinew,ntriangle,id) 499 if(id.gt.0) then 500 if(itriangle(id).eq.itrinew) then 501c write(*,*) '**circular path; no solution' 502 isol=-2 503 exit loop2 504 endif 505 endif 506 ntriangle=ntriangle+1 507 if(ntriangle.gt.ntriangle_) then 508c write(*,*) '**too many iterations' 509 isol=-3 510 exit loop2 511 endif 512 do k=ntriangle,id+2,-1 513 itriangle(k)=itriangle(k-1) 514 enddo 515 itriangle(id+1)=itrinew 516 itriold=itri 517 itri=itrinew 518 cycle loop2 519 endif 520 elseif(l.eq.3) then 521c write(*,*) '**regular solution' 522 isol=itri 523 exit loop2 524 endif 525 enddo 526 enddo loop2 527! 528! if an opposite triangle is found: check the distance 529! perpendicular to the triangle 530! 531 if(isol.gt.0) then 532 dist=dabs(straight(13,itri)*p(1)+ 533 & straight(14,itri)*p(2)+ 534 & straight(15,itri)*p(3)+ 535 & straight(16,itri)) 536 if(dist.gt.tietol(1,i)) isol=0 537 endif 538! 539! check whether distance is larger than tietol(1,i): 540! no element is generated 541! 542 if(isol.le.0) then 543! 544! no MPC is generated 545! 546 write(*,*) '*WARNING in gentiedmpc: no tied MPC' 547 write(*,*) ' generated for node ',node 548 if(isol.eq.0) then 549 write(*,*) ' master face too far away' 550 write(*,*) ' distance: ',dist 551 write(*,*) ' tolerance: ',tietol(1,i) 552 else 553 write(*,*) ' no corresponding master face' 554 write(*,*) ' found; tolerance: ', 555 & tolloc 556c & tietol(1,i) 557 endif 558 write(40,*) node 559 iwrite=1 560 else 561! 562 nelem=int(koncont(4,itri)/10.d0) 563 jface=koncont(4,itri)-10*nelem 564! 565 indexe=ipkon(nelem) 566 if(lakon(nelem)(4:4).eq.'2') then 567 nnodelem=8 568 nface=6 569 elseif(lakon(nelem)(4:4).eq.'8') then 570 nnodelem=4 571 nface=6 572 elseif(lakon(nelem)(4:5).eq.'10') then 573 nnodelem=6 574 nface=4 575 elseif(lakon(nelem)(4:4).eq.'4') then 576 nnodelem=3 577 nface=4 578 elseif(lakon(nelem)(4:5).eq.'15') then 579 if(jface.le.2) then 580 nnodelem=6 581 else 582 nnodelem=8 583 endif 584 nface=5 585 nope=15 586 elseif(lakon(nelem)(4:4).eq.'6') then 587 if(jface.le.2) then 588 nnodelem=3 589 else 590 nnodelem=4 591 endif 592 nface=5 593 nope=6 594 else 595 cycle 596 endif 597! 598! determining the nodes of the face 599! 600 if(nface.eq.4) then 601 do k=1,nnodelem 602 nodef(k)=kon(indexe+ifacet(k,jface)) 603 enddo 604 elseif(nface.eq.5) then 605 if(nope.eq.6) then 606 do k=1,nnodelem 607 nodef(k)=kon(indexe+ifacew1(k,jface)) 608 enddo 609 elseif(nope.eq.15) then 610 do k=1,nnodelem 611 nodef(k)=kon(indexe+ifacew2(k,jface)) 612 enddo 613 endif 614 elseif(nface.eq.6) then 615 do k=1,nnodelem 616 nodef(k)=kon(indexe+ifaceq(k,jface)) 617 enddo 618 endif 619! 620! attaching the node with coordinates in p 621! to the face 622! 623 do k=1,nnodelem 624 do l=1,3 625 pl(l,k)=co(l,nodef(k)) 626 enddo 627 enddo 628 call attach_2d(pl,p,nnodelem,ratio,dist,xi,et) 629! 630! adjusting the coordinates of the node (only 631! if the user did not specify ADJUST=NO on the 632! *TIE card) 633! 634 if(tietol(2,i).gt.0.d0) then 635 do k=1,3 636 co(k,node)=p(k) 637 enddo 638 endif 639! 640! generating MPC's 641! 642 do l=kstart,kend 643 idof=8*(node-1)+l 644 call nident(ikmpc,idof,nmpc,id) 645 if(id.gt.0) then 646 if(ikmpc(id).eq.idof) then 647 write(*,*) '*WARNING in gentiedmpc:' 648 write(*,*) ' DOF ',l,' of node ', 649 & node,' is not active;' 650 write(*,*) ' no tied constraint ', 651 & 'is generated' 652 write(40,*) node 653 iwrite=1 654 cycle 655 endif 656 endif 657! 658 call nident(ikboun,idof,nboun,id1) 659 if(id1.gt.0) then 660 if(ikboun(id1).eq.idof) then 661 write(*,*) '*WARNING in gentiedmpc:' 662 write(*,*) ' DOF ',l,' of node ', 663 & node,' is not active;' 664 write(*,*) ' no tied constraint ', 665 & 'is generated' 666 write(40,*) node 667 iwrite=1 668 cycle 669 endif 670 endif 671! 672 nmpc=nmpc+1 673 labmpc(nmpc)=' ' 674 ipompc(nmpc)=mpcfree 675! 676! updating ikmpc and ilmpc 677! 678 do m=nmpc,id+2,-1 679 ikmpc(m)=ikmpc(m-1) 680 ilmpc(m)=ilmpc(m-1) 681 enddo 682 ikmpc(id+1)=idof 683 ilmpc(id+1)=nmpc 684! 685 nodempc(1,mpcfree)=node 686 nodempc(2,mpcfree)=l 687 coefmpc(mpcfree)=1.d0 688 mpcfree=nodempc(3,mpcfree) 689 if(mpcfree.eq.0) then 690 write(*,*) 691 & '*ERROR in gentiedmpc: increase memmpc_' 692 call exit(201) 693 endif 694 do k=1,nnodelem 695 nodempc(1,mpcfree)=nodef(k) 696 nodempc(2,mpcfree)=l 697 coefmpc(mpcfree)=-ratio(k) 698 mpcfreeold=mpcfree 699 mpcfree=nodempc(3,mpcfree) 700 if(mpcfree.eq.0) then 701 write(*,*) 702 & '*ERROR in gentiedmpc: increase memmpc_' 703 call exit(201) 704 endif 705 enddo 706 nodempc(3,mpcfreeold)=0 707 enddo 708 endif 709! 710 enddo 711 endif 712 enddo 713 enddo 714! 715! number of tied MPC's 716! 717 nmpctied=nmpc-nmpctied 718! 719 if(iwrite.eq.1) then 720 write(*,*) '*INFO in gentiedmpc:' 721 write(*,*) ' failed nodes are stored in file' 722 write(*,*) ' WarnNodeMissTiedContact.nam' 723 write(*,*) ' This file can be loaded into' 724 write(*,*) ' an active cgx-session by typing' 725 write(*,*) 726 & ' read WarnNodeMissTiedContact.nam inp' 727 write(*,*) 728 close(40) 729 else 730 close(40,status='delete') 731 endif 732! 733 return 734! 735 end 736 737