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 normalsforequ_se(nk,co,iponoelfa,inoelfa,konfa, 20 & ipkonfa,lakonfa,ne,iponor,xnor,nodedesiinv,jobnamef, 21 & iponexp,nmpc,labmpc,ipompc,nodempc,ipretinfo,kon,ipkon,lakon, 22 & iponoel,inoel,iponor2d,knor2d,iponod2dto3d,ipoface,nodface) 23! 24! calculates normals on surface for mesh modification 25! purposes in an optimization loop 26! 27! during optimization the coordinates of the design variables 28! are changed leading to a changed geometry. In order to keep 29! a good quality mesh the other nodes may have to be moved as 30! well. The external shape in these nodes has to be kept, which 31! can be guaranteed by MPC's. These MPC's are based on the local 32! normal(s) in a node. At sharp corners more than one normal 33! may be necessary. 34! 35! the equations are stored in file jobname.equ 36! 37! the user can use this file for the appropriate mesh 38! modifications. 39! 40 implicit none 41! 42 character*132 jobnamef,fnequ 43 character*8 lakonfa(*),lakonfaloc 44 character*20 labmpc(*) 45 character*8 lakon(*) 46! 47 integer nk,iponoelfa(*),inoelfa(3,*),konfa(*),ipkonfa(*),ne, 48 & i,ndepnodes,index,nexp,nel,ielem,indexe,j,iel(100), 49 & jl(100),ial(100),ifi(100),indexx,k,l,ifix,nemin,jact,ixfree, 50 & node,iponor(*),nodedesiinv(*),len,ndet(3),nsort(3),two, 51 & three,iponexp(2,*),nmpc,ipompc(*),nodempc(3,*), 52 & node1,node2,node3,ipretinfo(*),ieq,pretflag,inoel(2,*),nope, 53 & nodepret,ixfreei,ixfreej,kon(*),ipkon(*),iponoel(*),iface, 54 & inode,ifaceq(8,6),ifacew(8,5),iposn,iponor2d(2,*),flag2d, 55 & knor2d(*),node2d,iponod2dto3d(2,*),nopesurf(8),ipoface(*), 56 & nodface(5,*),konl(20),nopem,ifaceqmid(6),ifacewmid(5),node3d 57! 58 real*8 co(3,*),xnor(*),xno(3,100),xi,et,coloc6(2,6),coloc8(2,8), 59 & xl(3,8),dd,xnoref(3),dot,xnorloc(3,3),det(3),sort(3),xdir, 60 & ydir,zdir 61! 62! 63! 64! In this routine the faces at the free surface play an 65! important role. They are considered to be like a layer of 66! shell elements. Therefore, the term "shell elements" in this 67! routine is basically equivalent to "external faces" 68! 69 70 data coloc6 /0.d0,0.d0,1.d0,0.d0,0.d0,1.d0,0.5d0,0.d0, 71 & 0.5d0,0.5d0,0.d0,0.5d0/ 72 data coloc8 /-1.d0,-1.d0,1.d0,-1.d0,1.d0,1.d0,-1.d0,1.d0, 73 & 0.d0,-1.d0,1.d0,0.d0,0.d0,1.d0,-1.d0,0.d0/ 74! 75 data ifaceqmid /0, 76 & 0, 77 & 5, 78 & 6, 79 & 7, 80 & 8/ 81! 82! 83! 84 data ifacewmid /0, 85 & 0, 86 & 4, 87 & 5, 88 & 6/ 89! 90! nodes per face for hex elements 91! 92 data ifaceq /4,3,2,1,11,10,9,12, 93 & 5,6,7,8,13,14,15,16, 94 & 1,2,6,5,9,18,13,17, 95 & 2,3,7,6,10,19,14,18, 96 & 3,4,8,7,11,20,15,19, 97 & 4,1,5,8,12,17,16,20/ 98! 99! nodes per face for quadratic wedge elements 100! 101 data ifacew /1,3,2,9,8,7,0,0, 102 & 4,5,6,10,11,12,0,0, 103 & 1,2,5,4,7,14,10,13, 104 & 2,3,6,5,8,15,11,14, 105 & 3,1,4,6,9,13,12,15/ 106! 107 two=2 108 three=3 109! 110 do len=1,132 111 if(jobnamef(len:len).eq.' ') exit 112 enddo 113 len=len-1 114 115 fnequ=jobnamef(1:len)//'.equ' 116 open(20,file=fnequ(1:len+4),status='unknown',err=100) 117 close(20,status='delete',err=101) 118 open(20,file=fnequ(1:len+4),status='unknown',err=100) 119 write(20,102) 120! write(20,103) 121 102 format('**SUMMARY OF EQUATIONS FOR MESH-UPDATE') 122 103 format('*EQUATION') 123! 124 ixfree=0 125! 126 do i=1,nk 127 ndepnodes=0 128 index=iponoelfa(i) 129 if(index.eq.0) cycle 130! 131! nexp indicates how many times the node was expanded 132! 133 nexp=0 134! 135! locating the shell elements to which node i belongs 136! 137 nel=0 138 do 139 if(index.eq.0) exit 140 ielem=inoelfa(1,index) 141 indexe=ipkonfa(ielem) 142 nel=nel+1 143 if(nel.gt.100) then 144 write(*,*) '*ERROR in normalsforequ_se: more ' 145 write(*,*) ' than 100 shell elements ' 146 write(*,*) ' share the same node' 147 call exit(201) 148 endif 149 j=inoelfa(2,index) 150 jl(nel)=j 151 iel(nel)=ielem 152 index=inoelfa(3,index) 153 enddo 154! 155 if(nel.gt.0) then 156 do j=1,nel 157 ial(j)=0 158 enddo 159! 160! estimate the normal 161! 162 do j=1,nel 163 indexe=ipkonfa(iel(j)) 164 indexx=iponor(indexe+jl(j)) 165 if(indexx.ge.0) then 166 do k=1,3 167 xno(k,j)=xnor(indexx+k) 168 enddo 169 ifi(j)=1 170 cycle 171 else 172 ifi(j)=0 173 endif 174! 175! local normal on the element (Jacobian) 176! 177 lakonfaloc=lakonfa(iel(j)) 178 if(lakonfa(iel(j))(2:2).eq.'3') then 179 xi=coloc6(1,jl(j)) 180 et=coloc6(2,jl(j)) 181 do k=1,3 182 node=konfa(indexe+k) 183 do l=1,3 184 xl(l,k)=co(l,node) 185 enddo 186 enddo 187 call norshell3(xi,et,xl,xno(1,j)) 188 elseif(lakonfa(iel(j))(2:2).eq.'4') then 189 xi=coloc8(1,jl(j)) 190 et=coloc8(2,jl(j)) 191 do k=1,4 192 node=konfa(indexe+k) 193 do l=1,3 194 xl(l,k)=co(l,node) 195 enddo 196 enddo 197 call norshell4(xi,et,xl,xno(1,j)) 198 elseif(lakonfa(iel(j))(2:2).eq.'6') then 199 xi=coloc6(1,jl(j)) 200 et=coloc6(2,jl(j)) 201 do k=1,6 202 node=konfa(indexe+k) 203 do l=1,3 204 xl(l,k)=co(l,node) 205 enddo 206 enddo 207 call norshell6(xi,et,xl,xno(1,j)) 208 elseif(lakonfa(iel(j))(2:2).eq.'8') then 209 xi=coloc8(1,jl(j)) 210 et=coloc8(2,jl(j)) 211 do k=1,8 212 node=konfa(indexe+k) 213 do l=1,3 214 xl(l,k)=co(l,node) 215 enddo 216 enddo 217 call norshell8(xi,et,xl,xno(1,j)) 218 endif 219! 220 dd=dsqrt(xno(1,j)**2+xno(2,j)**2+xno(3,j)**2) 221 if(dd.lt.1.d-10) then 222 write(*,*) '*ERROR in normalsforequ_se: size ' 223 write(*,*) ' of estimatedshell normal in 224 &node ',i,' element ',iel(j) 225 write(*,*) ' is smaller than 1.e-10' 226 call exit(201) 227 endif 228 do k=1,3 229 xno(k,j)=xno(k,j)/dd 230 enddo 231 enddo 232! 233 do 234! 235! determining a fixed normal which was not treated yet, 236! or, if none is left, the minimum element number of all 237! elements containing node i and for which no normal was 238! determined yet 239! 240! if ial(j)=0: the normal on this element has not been 241! treated yet 242! if ial(j)=2: normal has been treated 243! 244 ifix=0 245 nemin=ne+1 246 do j=1,nel 247 if(ial(j).ne.0) cycle 248 if(ifi(j).eq.1) then 249 jact=j 250 ifix=1 251 exit 252 endif 253 enddo 254 if(ifix.eq.0) then 255 do j=1,nel 256 if(ial(j).eq.0) then 257 if(iel(j).lt.nemin) then 258 nemin=iel(j) 259 jact=j 260 endif 261 endif 262 enddo 263 if(nemin.eq.ne+1) exit 264 endif 265! 266 do j=1,3 267 xnoref(j)=xno(j,jact) 268 enddo 269! 270! determining all elements whose normal in node i makes an 271! angle smaller than 0.5 or 20 degrees with the reference normal, 272! depending whether the reference normal was given by the 273! user or is being calculated; the thickness and offset must 274! also fit. 275! 276! if ial(j)=1: normal on element is being treated now 277! 278 do j=1,nel 279 if(ial(j).eq.2) cycle 280 if(j.eq.jact) then 281 ial(jact)=1 282 else 283 dot=xno(1,j)*xnoref(1)+xno(2,j)*xnoref(2)+ 284 & xno(3,j)*xnoref(3) 285 if(ifix.eq.0) then 286 if(dot.gt.0.939693d0)then 287 if((lakonfa(iel(j))(1:3).eq. 288 & lakonfa(iel(jact))(1:3)) 289 & .or. 290 & ((lakonfa(iel(j))(1:1).eq.'S').and. 291 & (lakonfa(iel(jact))(1:1).eq.'S'))) 292 & ial(j)=1 293 else 294 if((lakonfa(iel(j))(1:1).eq.'S').and. 295 & (lakonfa(iel(jact))(1:1).eq.'S')) then 296! 297! if the normals have the opposite 298! direction, the expanded nodes are on a 299! straight line 300! 301 if(dot.le.-0.999962d0) then 302 write(*,*) '*INFO in gen3dnor: in some 303 &nodes opposite normals are defined' 304 endif 305 endif 306 endif 307 else 308 if(dot.gt.0.999962d0) then 309 if((lakonfa(iel(j))(1:3).eq. 310 & lakonfa(iel(jact))(1:3)) 311 & .or. 312 & ((lakonfa(iel(j))(1:1).eq.'S').and. 313 & (lakonfa(iel(jact))(1:1).eq.'S'))) 314 & ial(j)=1 315 else 316 if((lakonfa(iel(j))(1:1).eq.'S').and. 317 & (lakonfa(iel(jact))(1:1).eq.'S')) then 318! 319! if the normals have the opposite 320! direction, the expanded nodes are on a 321! straight line 322! 323 if(dot.le.-0.999962d0) then 324 write(*,*) '*INFO in gen3dnor: in some 325 &nodes opposite normals are defined' 326 endif 327 endif 328 endif 329 endif 330 endif 331 enddo 332! 333! determining the mean normal for the selected elements 334! 335 if(ifix.eq.0) then 336 do j=1,3 337 xnoref(j)=0.d0 338 enddo 339 do j=1,nel 340 if(ial(j).eq.1) then 341 do k=1,3 342 xnoref(k)=xnoref(k)+xno(k,j) 343 enddo 344 endif 345 enddo 346 dd=dsqrt(xnoref(1)**2+xnoref(2)**2+xnoref(3)**2) 347 if(dd.lt.1.d-10) then 348 write(*,*) '*ERROR in gen3dnor: size of' 349 write(*,*) ' estimated shell normal is' 350 write(*,*) ' smaller than 1.e-10' 351 call exit(201) 352 endif 353 do j=1,3 354 xnoref(j)=xnoref(j)/dd 355 enddo 356 endif 357! 358! updating the pointers iponor 359! 360 nexp=nexp+1 361 do j=1,nel 362 if(ial(j).eq.1) then 363 ial(j)=2 364 if(ifix.eq.0) then 365 iponor(ipkonfa(iel(j))+jl(j))=ixfree 366 elseif(j.ne.jact) then 367 iponor(ipkonfa(iel(j))+jl(j))= 368 & iponor(ipkonfa(iel(jact))+jl(jact)) 369 endif 370 endif 371 enddo 372! 373! storing the normal in xnor and generating 3 nodes 374! for knor 375! 376 if(ifix.eq.0) then 377 do j=1,3 378 xnor(ixfree+j)=xnoref(j) 379 enddo 380 ixfree=ixfree+3 381 endif 382! 383 enddo 384 endif 385! 386! save nexp and ixfree 387! 388 iponexp(1,i)=nexp 389 iponexp(2,i)=ixfree 390! 391 enddo 392! 393! find nodes created by "*PRETENSION SECTION" 394! 395! find pretension node if existing 396! 397 pretflag=0 398 do i=1,nmpc 399 if(labmpc(i)(1:10).eq.'PRETENSION') then 400 pretflag=1 401 index=ipompc(i) 402 index=nodempc(3,index) 403 index=nodempc(3,index) 404 nodepret=nodempc(1,index) 405 pretflag=1 406 exit 407 endif 408 enddo 409! 410 if(pretflag.eq.1) then 411 do i=1,nmpc 412 if(labmpc(i)(1:11).eq.'THERMALPRET') cycle 413! 414 ieq=0 415 index=ipompc(i) 416 if(index.eq.0) cycle 417 node1=nodempc(1,index) 418 index=nodempc(3,index) 419 node2=nodempc(1,index) 420 index=nodempc(3,index) 421 node3=nodempc(1,index) 422 if(node3.eq.nodepret) then 423 ipretinfo(node2)=node1 424 ipretinfo(node1)=-1 425 endif 426 enddo 427 endif 428! 429! correct nodes on free pretension surface" 430! 431 do i=1,nk 432 if(ipretinfo(i).le.0) cycle 433! 434 nexp=iponexp(1,i) 435 ixfreei=iponexp(2,i) 436 ixfreej=iponexp(2,ipretinfo(i)) 437! 438 do j=1,nexp 439 k=j*3-3 440 zdir=xnor(ixfreei+1-1-k)+xnor(ixfreej+1-1-k) 441 ydir=xnor(ixfreei+1-2-k)+xnor(ixfreej+1-2-k) 442 xdir=xnor(ixfreei+1-3-k)+xnor(ixfreej+1-3-k) 443 dd=(xdir)**2+(ydir)**2+(zdir)**2 444! 445 if(dd.gt.1.0e-12) then 446 ipretinfo(i)=0 447 endif 448! 449 enddo 450! 451 enddo 452! 453!--------------------------------------------------------------------------- 454! 455! write equations in file "jobname.equ" 456! in case of a 2D model just write the node numbers in the file 457! 458 do i=1,nk 459 flag2d=0 460! 461! check for additional pretension nodes 462! 463 if(ipretinfo(i).ne.0) cycle 464! 465! check if node is a designvariable 466! 467 if(nodedesiinv(i).eq.0) then 468! 469! consideration of plain stress/strain 2d-elements 470! and rotational symmetry elements 471! 472 if(iponoel(i).eq.0) cycle 473 ielem=inoel(1,iponoel(i)) 474 if((lakon(ielem)(7:7).eq.'A').or. 475 & (lakon(ielem)(7:7).eq.'S').or. 476 & (lakon(ielem)(7:7).eq.'E')) then 477! 478 if(lakon(ielem)(4:5).eq.'20') then 479 nope=20 480 elseif (lakon(ielem)(4:4).eq.'8') then 481 nope=8 482 elseif (lakon(ielem)(4:5).eq.'15') then 483 nope=15 484 else 485 cycle 486 endif 487! 488 indexe=ipkon(ielem) 489 do inode=1,nope 490 if(i.eq.kon(indexe+inode)) then 491 exit 492 endif 493 enddo 494 if(lakon(ielem)(4:5).eq.'20') then 495 if((inode.ne.ifaceq(6,3)).and. 496 & (inode.ne.ifaceq(8,3)).and. 497 & (inode.ne.ifaceq(6,5)).and. 498 & (inode.ne.ifaceq(8,5))) cycle 499! 500! replace 3D node number by 2D node number 501! 502c write(*,*) 'normalsforequ_se ',i,iponoelfa(i) 503c write(*,*) 'normalsforequ_se ',i,inoelfa(1,iponoelfa(i)) 504c write(*,*) 'normalsforequ_se ',(konfa(ipkonfa(iface)+j),j=1,4) 505c iface=inoelfa(1,iponoelfa(i)) 506c if(iface.le.2) cycle 507 node=kon(indexe+inode+4) 508 flag2d=1 509 elseif(lakon(ielem)(4:5).eq.'15') then 510 if((inode.ne.ifacew(6,3)).and. 511 & (inode.ne.ifacew(8,3)).and. 512 & (inode.ne.ifacew(6,4))) cycle 513! 514! replace 3D node number by 2D node number 515! 516c iface=inoelfa(1,iponoelfa(i)) 517c if(iface.le.2) cycle 518 node=kon(indexe+inode+3) 519 flag2d=1 520 else 521 cycle 522 endif 523 elseif(lakon(ielem)(7:7).eq.'L') then 524! 525! no output for shell elements necessary 526! 527 cycle 528 else 529! 530! in case of a 3D model no change of node number 531 node=i 532 endif 533! 534! write equations in case nexp is greater or equal 3 535! 536 nexp=iponexp(1,i) 537 ixfree=iponexp(2,i) 538! 539 if((nexp.ge.3).and.(flag2d.eq.0)) then 540 do j=1,3 541 write(20,106) 1 542 write(20,105) node,j,1 543 enddo 544! 545! write equations in case nexp is 1 546! 547 elseif((nexp.eq.1).and.(flag2d.eq.0)) then 548 j=1 549 do l=1,3 550 xnorloc(4-l,j)=xnor(ixfree+1-l) 551 sort(4-l)=dabs(xnor(ixfree+1-l)) 552 nsort(4-l)=4-l 553 enddo 554 call dsort(sort,nsort,three,two) 555 write(20,106) 3 556 write(20,104) node,nsort(3),xnorloc(nsort(3),1), 557 & node,nsort(2),xnorloc(nsort(2),1), 558 & node,nsort(1),xnorloc(nsort(1),1) 559! 560! write equations in case nexp is 2 561! 562 elseif((nexp.eq.2).and.(flag2d.eq.0)) then 563 do j=1,nexp 564 k=j*3-3 565 do l=1,3 566 xnorloc(4-l,j)=xnor(ixfree+1-l-k) 567 enddo 568 enddo 569 ndet(1)=1 570 ndet(2)=2 571 ndet(3)=3 572 det(1)=dabs(xnorloc(1,1)*xnorloc(2,2)- 573 & xnorloc(1,2)*xnorloc(2,1)) 574 det(2)=dabs(xnorloc(1,1)*xnorloc(3,2)- 575 & xnorloc(1,2)*xnorloc(3,1)) 576 det(3)=dabs(xnorloc(2,1)*xnorloc(3,2)- 577 & xnorloc(2,2)*xnorloc(3,1)) 578 call dsort(det,ndet,three,two) 579 580 if(ndet(3).eq.1) then 581! if((dabs(xnorloc(1,1)).ge.dabs(xnorloc(2,1))).and. 582! & (dabs(xnorloc(2,2)).gt.1.d-5)) then 583 if((dabs(xnorloc(1,1)).gt.1.d-5).and. 584 & (dabs(xnorloc(2,2)).gt.1.d-5)) then 585 write(20,106) 3 586 write(20,104) node,1,xnorloc(1,1), 587 & node,2,xnorloc(2,1),node,3,xnorloc(3,1) 588 write(20,106) 3 589 write(20,104) node,2,xnorloc(2,2), 590 & node,1,xnorloc(1,2),node,3,xnorloc(3,2) 591 else 592 write(20,106) 3 593 write(20,104) node,2,xnorloc(2,1), 594 & node,1,xnorloc(1,1),node,3,xnorloc(3,1) 595 write(20,106) 3 596 write(20,104) node,1,xnorloc(1,2), 597 & node,2,xnorloc(2,2),node,3,xnorloc(3,2) 598 endif 599 elseif(ndet(3).eq.2) then 600! if((dabs(xnorloc(1,1)).ge.dabs(xnorloc(3,1))).and. 601! & (dabs(xnorloc(3,2)).gt.1.d-5)) then 602 if((dabs(xnorloc(1,1)).gt.1.d-5).and. 603 & (dabs(xnorloc(3,2)).gt.1.d-5)) then 604 write(20,106) 3 605 write(20,104) node,1,xnorloc(1,1), 606 & node,3,xnorloc(3,1),node,2,xnorloc(2,1) 607 write(20,106) 3 608 write(20,104) node,3,xnorloc(3,2), 609 & node,1,xnorloc(1,2),node,2,xnorloc(2,2) 610 else 611 write(20,106) 3 612 write(20,104) node,3,xnorloc(3,1), 613 & node,1,xnorloc(1,1),node,2,xnorloc(2,1) 614 write(20,106) 3 615 write(20,104) node,1,xnorloc(1,2), 616 & node,3,xnorloc(3,2),node,2,xnorloc(2,2) 617 endif 618 elseif(ndet(3).eq.3) then 619! if((dabs(xnorloc(2,1)).ge.dabs(xnorloc(3,1))).and. 620! & (dabs(xnorloc(3,2)).gt.1.d-5)) then 621 if((dabs(xnorloc(2,1)).gt.1.d-5).and. 622 & (dabs(xnorloc(3,2)).gt.1.d-5)) then 623 write(20,106) 3 624 write(20,104) node,2,xnorloc(2,1), 625 & node,3,xnorloc(3,1),node,1,xnorloc(1,1) 626 write(20,106) 3 627 write(20,104) node,3,xnorloc(3,2), 628 & node,2,xnorloc(2,2),node,1,xnorloc(1,2) 629 else 630 write(20,106) 3 631 write(20,104) node,3,xnorloc(3,1), 632 & node,2,xnorloc(2,1),node,1,xnorloc(1,1) 633 write(20,106) 3 634 write(20,104) node,2,xnorloc(2,2), 635 & node,3,xnorloc(3,2),node,1,xnorloc(1,2) 636 endif 637 endif 638! 639! WORKAROUND: MPC's in combination with expanded 2D models does not work 640! in case of expanded 2D models create a set with all surface nodes 641! which are not in the designvariables set. These nodes are fully 642! constrained 643! 644 elseif(flag2d.eq.1) then 645 write(20,'(i10,a1)') node,',' 646 endif 647 elseif(nodedesiinv(i).eq.1) then 648 if(iponoel(i).eq.0) cycle 649 ielem=inoel(1,iponoel(i)) 650 if((lakon(ielem)(7:7).eq.'A').or. 651 & (lakon(ielem)(7:7).eq.'S').or. 652 & (lakon(ielem)(7:7).eq.'L').or. 653 & (lakon(ielem)(7:7).eq.'E')) then 654! 655 nodedesiinv(iponod2dto3d(1,i))=-1 656 nodedesiinv(iponod2dto3d(2,i))=-1 657 endif 658 endif 659! 660 enddo 661! 662!------------------------------------------------------------------------------- 663! 664! in case of plain strain/stress/axi 2D models write midnodes belonging 665! to these 2D elements to file. This naturally only applies to 666! quadratic elements 667! 668 do i=1,nk 669 670 if(ipoface(i).eq.0) cycle 671 indexe=ipoface(i) 672! 673 do 674 ielem=nodface(3,indexe) 675 iface=nodface(4,indexe) 676! 677 if((lakon(ielem)(7:7).eq.'A').or. 678 & (lakon(ielem)(7:7).eq.'S').or. 679 & (lakon(ielem)(7:7).eq.'E')) then 680! 681! faces in z-direction (expansion direction) do not play 682! a role in the optimization of plane stress/strain/axi 683! elements (corresponds to iface=1 and iface=2) 684! 685 if(iface.gt.2) then 686! 687 if(lakon(ielem)(4:5).eq.'20') then 688 nope=20 689 nopem=8 690 elseif(lakon(ielem)(4:5).eq.'15') then 691 nope=15 692 nopem=8 693 else 694 indexe=nodface(5,indexe) 695 if(indexe.eq.0) then 696 exit 697 else 698 cycle 699 endif 700 endif 701c! 702c! actual position of the nodes belonging to the 703c! surface and surface normal 704c do j=1,nope 705c konl(j)=kon(ipkon(ielem)+j) 706c enddo 707c! 708c do j=1,nopem 709c nopesurf(j)=konl(ifaceq(j,iface)) 710c do k=1,3 711c xl(k,j)=co(k,nopesurf(j)) 712c enddo 713c enddo 714c xi=0.d0 715c et=0.d0 716c call norshell8(xi,et,xl,xno(1,1)) 717c dd=dsqrt(xno(1,1)**2+xno(2,1)**2+ 718c & xno(3,1)**2) 719c if(dd.gt.0.d0) then 720c do j=1,3 721c xno(j,1)=xno(j,1)/dd 722c enddo 723c endif 724! 725! node number and equation of the 2D node 726! 727 if(nope.eq.20) then 728 node2d=kon(ipkon(ielem)+nope+ifaceqmid(iface)) 729 iposn=iponor2d(2,ipkon(ielem)+ifaceqmid(iface)) 730 elseif(nope.eq.15) then 731 node2d=kon(ipkon(ielem)+nope+ifacewmid(iface)) 732 iposn=iponor2d(2,ipkon(ielem)+ifacewmid(iface)) 733 endif 734! 735! 3D-equivalent of the 2D-design variable 736! The user defines 2D nodes of plane stress/strain/axi 737! elements as design variables. Internally, these are 738! replaced by the first node in the 3D expansion 739! 740 node3d=knor2d(iposn+1) 741! 742! write the 2D node to file if it is not a design variable 743! 744 if(nodedesiinv(node3d).eq.0) then 745 write(20,'(i10,a1)') node2d,',' 746! write(20,106) 3 747! write(20,104) node2d,1,xno(1,1), 748! & node2d,2,xno(2,1), 749! & node2d,3,xno(3,1) 750 endif 751 endif 752 endif 753 indexe=nodface(5,indexe) 754 if(indexe.eq.0) exit 755! 756 enddo 757 enddo 758! 759!------------------------------------------------------------------------------- 760! 761! in case of 2D shell models fix all nodes which are no design variable 762! 763! do i=1,nk 764! if(ipoface(i).eq.0) cycle 765! if(nodedesiinv(i).eq.0) then 766! write(20,'(i10,a1)') i,',' 767! endif 768! enddo 769! 770 do i=1,nk 771 if(nodedesiinv(i).eq.-1) then 772 nodedesiinv(i)=0 773 endif 774 enddo 775! 776 close(20) 777 return 778! 779 104 format(3(i10,",",i1,",",e20.13,",")) 780 105 format(1(i10,",",i1,",",i1,",")) 781 106 format(i1) 782 107 format(2(i10,",",i1,",",e20.13,",")) 783! 784 100 write(*,*) '*ERROR in openfile: could not open file ', 785 & fnequ(1:len+4) 786 call exit(201) 787 101 write(*,*) '*ERROR in openfile: could not delete file ', 788 & fnequ(1:len+4) 789 call exit(201) 790! 791 end 792