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 normalsonsurface_robust(ipkon,kon,lakon,extnor,co,nk, 20 & ipoface,nodface,nactdof,mi,nodedesiinv,iregion, 21 & iponoelfa,ndesi,nodedesi,iponod2dto3d,ikboun,nboun, 22 & ne2d) 23! 24! calculating the normal direction onto the external surface; 25! the design variables are moved in this direction. 26! 27! if the design variables constitute a region, i.e. they 28! are constitute a set of faces, the normals at the boundary 29! of this set is determined based on the faces belonging to 30! this set only (so faces external to this set do not 31! contribute to the normal at the boundary) 32! 33! if it is not a region, the normal in a node is the mean 34! of the normal of all external faces to which this node 35! belongs, no matter how many other design variables 36! belong to these faces (if any) 37! 38 implicit none 39! 40 character*8 lakon(*) 41! 42 integer j,nelem,jface,indexe,ipkon(*),kon(*),nopem,node, 43 & ifaceq(8,6),ifacet(6,4),ifacew1(4,5),ifacew2(8,5),ne2d, 44 & konl(26),ipoface(*),nodface(5,*),mi(*),nodedesiinv(*), 45 & nactdof(0:mi(2),*),nopesurf(9),nnodes,iregion,nope, 46 & nopedesi,l,m,iflag,k,nk,iponoelfa(*),ndesi,nodedesi(*), 47 & nodemid,nodeboun1,nodeboun2,iponod2dto3d(2,*), 48 & ishift,expandhex(20),expandwed(15),konl2d(26),ikboun(*), 49 & idof,nboun,id,node2d 50! 51 real*8 extnor(3,*),xsj2(3),shp2(7,9),xs2(3,2),xi,et,dd, 52 & xquad(2,9),xtri(2,7),xl2m(3,9),co(3,*) 53! 54! local node numbers for relationship between 2D and 3D elements 55! 56 data expandhex /1,2,3,4, 57 & 1,2,3,4, 58 & 5,6,7,8, 59 & 5,6,7,8, 60 & 1,2,3,4/ 61 data expandwed /1,2,3, 62 & 1,2,3, 63 & 4,5,6, 64 & 4,5,6, 65 & 1,2,3/ 66! 67! nodes per face for hex elements 68! 69 data ifaceq /4,3,2,1,11,10,9,12, 70 & 5,6,7,8,13,14,15,16, 71 & 1,2,6,5,9,18,13,17, 72 & 2,3,7,6,10,19,14,18, 73 & 3,4,8,7,11,20,15,19, 74 & 4,1,5,8,12,17,16,20/ 75! 76! nodes per face for tet elements 77! 78 data ifacet /1,3,2,7,6,5, 79 & 1,2,4,5,9,8, 80 & 2,3,4,6,10,9, 81 & 1,4,3,8,10,7/ 82! 83! nodes per face for linear wedge elements 84! 85 data ifacew1 /1,3,2,0, 86 & 4,5,6,0, 87 & 1,2,5,4, 88 & 2,3,6,5, 89 & 3,1,4,6/ 90! 91! nodes per face for quadratic wedge elements 92! 93 data ifacew2 /1,3,2,9,8,7,0,0, 94 & 4,5,6,10,11,12,0,0, 95 & 1,2,5,4,7,14,10,13, 96 & 2,3,6,5,8,15,11,14, 97 & 3,1,4,6,9,13,12,15/ 98! 99! new added data for the local coodinates for nodes 100! 101 data xquad /-1.d0,-1.d0, 102 & 1.d0,-1.d0, 103 & 1.d0,1.d0, 104 & -1.d0,1.d0, 105 & 0.d0,-1.d0, 106 & 1.d0,0.d0, 107 & 0.d0,1.d0, 108 & -1.d0,0.d0, 109 & 0.d0,0.d0/ 110! 111 data xtri /0.d0,0.d0, 112 & 1.d0,0.d0, 113 & 0.d0,1.d0, 114 & .5d0,0.d0, 115 & .5d0,.5d0, 116 & 0.d0,.5d0, 117 & 0.333333333333333d0,0.333333333333333d0/ 118! 119 data iflag /2/ 120! 121! calculation of the normal to the external surface; 122! each external face contributes its normal 123! to each node belonging to the face providing that 124! 1) the node is not a design variable OR 125! 2) the node is a design variable and belongs to a design face 126! A design face is an external face for which more than nopedesi 127! nodes are design variables 128! 129! The appropriate normal component is set to zero for fixed dofs 130! 131 do j=1,nk 132! 133 if(ipoface(j).eq.0) cycle 134 indexe=ipoface(j) 135! 136 do 137! 138 nelem=nodface(3,indexe) 139 jface=nodface(4,indexe) 140c write(*,*) 'normalsonsurface_se ',j,nelem,jface 141! 142 if((lakon(nelem)(7:7).eq.'A').or. 143 & (lakon(nelem)(7:7).eq.'S').or. 144 & (lakon(nelem)(7:7).eq.'E')) then 145! 146! for plane stress/strain/axi only faces 147! 3 and higher are taken into account for the normal 148! 149 if(jface.le.2) then 150 indexe=nodface(5,indexe) 151 if(indexe.eq.0) then 152 exit 153 else 154 cycle 155 endif 156 endif 157 elseif(lakon(nelem)(7:7).eq.'L') then 158! 159! for shells only faces 2 and lower 160! taken into account for the normal 161! 162 if(jface.gt.2) then 163 indexe=nodface(5,indexe) 164 if(indexe.eq.0) then 165 exit 166 else 167 cycle 168 endif 169 endif 170 endif 171! 172! nopem: # of nodes in the surface 173! nope: # of nodes in the element 174! 175 if(lakon(nelem)(4:4).eq.'8') then 176 nopem=4 177 nope=8 178 nopedesi=3 179 ishift=8 180 elseif(lakon(nelem)(4:5).eq.'20') then 181 nopem=8 182 nope=20 183 nopedesi=5 184 ishift=20 185 elseif(lakon(nelem)(4:5).eq.'10') then 186 nopem=6 187 nope=10 188 nopedesi=4 189 elseif(lakon(nelem)(4:4).eq.'4') then 190 nopem=3 191 nope=4 192 nopedesi=3 193! 194! treatment of wedge faces 195! 196 elseif(lakon(nelem)(4:4).eq.'6') then 197 nope=6 198 if(jface.le.2) then 199 nopem=3 200 else 201 nopem=4 202 endif 203 nopedesi=3 204 ishift=6 205 elseif(lakon(nelem)(4:5).eq.'15') then 206 nope=15 207 if(jface.le.2) then 208 nopem=6 209 nopedesi=4 210 else 211 nopem=8 212 nopedesi=5 213 endif 214 ishift=15 215 endif 216 if(iregion.eq.0) then 217 nopedesi=0 218 endif 219! 220! actual position of the nodes belonging to the 221! surface 222! 223 if((lakon(nelem)(7:7).eq.'A').or. 224 & (lakon(nelem)(7:7).eq.'S').or. 225 & (lakon(nelem)(7:7).eq.'E')) then 226 if((lakon(nelem)(4:5).eq.'20').or. 227 & (lakon(nelem)(4:5).eq.'8 ')) then 228 do k=1,nope 229 konl(k)=kon(ipkon(nelem)+k) 230 konl2d(k)=kon(ipkon(nelem)+ishift+expandhex(k)) 231 enddo 232 elseif((lakon(nelem)(4:5).eq.'15').or. 233 & (lakon(nelem)(4:5).eq.'6 ')) then 234 do k=1,nope 235 konl(k)=kon(ipkon(nelem)+k) 236 konl2d(k)=kon(ipkon(nelem)+ishift+expandwed(k)) 237 enddo 238 endif 239 else 240 do k=1,nope 241 konl(k)=kon(ipkon(nelem)+k) 242 enddo 243 endif 244! 245c do k=1,nope 246c konl(k)=kon(ipkon(nelem)+k) 247c if((lakon(nelem)(7:7).eq.'A').or. 248c & (lakon(nelem)(7:7).eq.'S').or. 249c & (lakon(nelem)(7:7).eq.'E')) then 250c if((lakon(nelem)(4:5).eq.'20').or. 251c & (lakon(nelem)(4:5).eq.'8 ')) then 252c konl2d(k)=kon(ipkon(nelem)+ishift+expandhex(k)) 253c elseif((lakon(nelem)(4:5).eq.'15').or. 254c & (lakon(nelem)(4:5).eq.'6 ')) then 255c konl2d(k)=kon(ipkon(nelem)+ishift+expandwed(k)) 256c endif 257c endif 258c enddo 259! 260 if((nope.eq.20).or.(nope.eq.8)) then 261 do m=1,nopem 262 nopesurf(m)=konl(ifaceq(m,jface)) 263 do k=1,3 264 xl2m(k,m)=co(k,nopesurf(m)) 265 enddo 266 enddo 267 elseif((nope.eq.10).or.(nope.eq.4)) 268 & then 269 do m=1,nopem 270 nopesurf(m)=konl(ifacet(m,jface)) 271 do k=1,3 272 xl2m(k,m)=co(k,nopesurf(m)) 273 enddo 274 enddo 275 elseif(nope.eq.15) then 276 do m=1,nopem 277 nopesurf(m)=konl(ifacew2(m,jface)) 278 do k=1,3 279 xl2m(k,m)=co(k,nopesurf(m)) 280 enddo 281 enddo 282 else 283 do m=1,nopem 284 nopesurf(m)=konl(ifacew1(m,jface)) 285 do k=1,3 286 xl2m(k,m)=co(k,nopesurf(m)) 287 enddo 288 enddo 289 endif 290! 291! sum up how many designvariables are on that surface 292! 293 nnodes=0 294 do m=1,nopem 295 if(nodedesiinv(nopesurf(m)).eq.1) then 296 nnodes=nnodes+1 297 endif 298 enddo 299! 300! calculate the normal vector in the nodes belonging to the surface 301! 302 if(nopem.eq.8) then 303 do m=1,nopem 304 xi=xquad(1,m) 305 et=xquad(2,m) 306 call shape8q(xi,et,xl2m,xsj2,xs2,shp2,iflag) 307 dd=dsqrt(xsj2(1)*xsj2(1) + xsj2(2)*xsj2(2) 308 & + xsj2(3)*xsj2(3)) 309 xsj2(1)=xsj2(1)/dd 310 xsj2(2)=xsj2(2)/dd 311 xsj2(3)=xsj2(3)/dd 312! 313 if(nope.eq.20) then 314 node=konl(ifaceq(m,jface)) 315 elseif(nope.eq.15) then 316 node=konl(ifacew2(m,jface)) 317 endif 318c write(*,*) 'normalsonsurface_se',node,nelem,jface, 319c &xsj2(1), 320c &xsj2(2),xsj2(3),lakon(nelem) 321 if((nodedesiinv(node).eq.0).or. 322 & ((nodedesiinv(node).eq.1).and. 323 & (nnodes.ge.nopedesi))) then 324 extnor(1,node)=extnor(1,node) 325 & +xsj2(1) 326 extnor(2,node)=extnor(2,node) 327 & +xsj2(2) 328 extnor(3,node)=extnor(3,node) 329 & +xsj2(3) 330c write(*,*) 'normalsonsurface_se ',extnor(3,node) 331! 332! in case of plain strain/stress/axi elements 333! not considering the x3-direction and the 334! directions with fixed displacements 335! 336 if((lakon(nelem)(7:7).eq.'A').or. 337 & (lakon(nelem)(7:7).eq.'S').or. 338 & (lakon(nelem)(7:7).eq.'E')) then 339 if(nope.eq.20) then 340 node2d=konl2d(ifaceq(m,jface)) 341 elseif(nope.eq.15) then 342 node2d=konl2d(ifacew2(m,jface)) 343 endif 344 do l=1,2 345 idof=8*(node2d-1)+l 346 call nident(ikboun,idof,nboun,id) 347 if(id.gt.0) then 348 if(ikboun(id).eq.idof) then 349 extnor(l,node)=0.d0 350 endif 351 endif 352 enddo 353 extnor(3,node)=0.d0 354! 355! else, all the directions with fixed 356! displacements are not considered 357! 358! elseif(lakon(nelem)(7:7).eq.' ') then 359! do l=1,3 360! if(nactdof(l,node).le.0) then 361! extnor(l,node)=0.d0 362! endif 363! enddo 364 endif 365 endif 366 enddo 367 elseif(nopem.eq.4) then 368 do m=1,nopem 369 xi=xquad(1,m) 370 et=xquad(2,m) 371 call shape4q(xi,et,xl2m,xsj2,xs2,shp2,iflag) 372 dd=dsqrt(xsj2(1)*xsj2(1) + xsj2(2)*xsj2(2) 373 & + xsj2(3)*xsj2(3)) 374 xsj2(1)=xsj2(1)/dd 375 xsj2(2)=xsj2(2)/dd 376 xsj2(3)=xsj2(3)/dd 377! 378 if(nope.eq.8) then 379 node=konl(ifaceq(m,jface)) 380 elseif(nope.eq.6) then 381 node=konl(ifacew1(m,jface)) 382 endif 383c write(*,*) 'normalsonsurface_se',node,nelem,jface, 384c &xsj2(1), 385c &xsj2(2),xsj2(3),lakon(nelem) 386 if((nodedesiinv(node).eq.0).or. 387 & ((nodedesiinv(node).eq.1).and. 388 & (nnodes.ge.nopedesi))) then 389c write(*,*) 'normalsonsurface_se accepted' 390 extnor(1,node)=extnor(1,node) 391 & +xsj2(1) 392 extnor(2,node)=extnor(2,node) 393 & +xsj2(2) 394 extnor(3,node)=extnor(3,node) 395 & +xsj2(3) 396! 397! in case of plain strain/stress/axi elements 398! not considering the x3-direction and the 399! directions with fixed displacements 400! 401 if((lakon(nelem)(7:7).eq.'A').or. 402 & (lakon(nelem)(7:7).eq.'S').or. 403 & (lakon(nelem)(7:7).eq.'E')) then 404 if(nope.eq.8) then 405 node2d=konl2d(ifaceq(m,jface)) 406 elseif(nope.eq.6) then 407 node2d=konl2d(ifacew1(m,jface)) 408 endif 409 do l=1,2 410 idof=8*(node2d-1)+l 411 call nident(ikboun,idof,nboun,id) 412 if(id.gt.0) then 413 if(ikboun(id).eq.idof) then 414 extnor(l,node)=0.d0 415 endif 416 endif 417 enddo 418 extnor(3,node)=0.d0 419! 420! else, all the directions with fixed 421! displacements are not considered 422! 423! elseif(lakon(nelem)(7:7).eq.' ') then 424! do l=1,3 425! if(nactdof(l,node).le.0) then 426! extnor(l,node)=0.d0 427! endif 428! enddo 429 endif 430 endif 431 enddo 432 elseif(nopem.eq.6) then 433 do m=1,nopem 434 xi=xtri(1,m) 435 et=xtri(2,m) 436 call shape6tri(xi,et,xl2m,xsj2,xs2,shp2,iflag) 437 dd=dsqrt(xsj2(1)*xsj2(1) + xsj2(2)*xsj2(2) 438 & + xsj2(3)*xsj2(3)) 439 xsj2(1)=xsj2(1)/dd 440 xsj2(2)=xsj2(2)/dd 441 xsj2(3)=xsj2(3)/dd 442! 443 if(nope.eq.10) then 444 node=konl(ifacet(m,jface)) 445 elseif(nope.eq.15) then 446 node=konl(ifacew2(m,jface)) 447 endif 448 if((nodedesiinv(node).eq.0).or. 449 & ((nodedesiinv(node).eq.1).and. 450 & (nnodes.ge.nopedesi))) then 451 extnor(1,node)=extnor(1,node) 452 & +xsj2(1) 453 extnor(2,node)=extnor(2,node) 454 & +xsj2(2) 455 extnor(3,node)=extnor(3,node) 456 & +xsj2(3) 457! 458! in case of plain strain/stress/axi elements 459! not considering the x3-direction and the 460! directions with fixed displacements 461! 462 if((lakon(nelem)(7:7).eq.'A').or. 463 & (lakon(nelem)(7:7).eq.'S').or. 464 & (lakon(nelem)(7:7).eq.'E')) then 465 if(nope.eq.10) then 466 node2d=konl2d(ifacet(m,jface)) 467 elseif(nope.eq.15) then 468 node2d=konl2d(ifacew2(m,jface)) 469 endif 470 do l=1,2 471 idof=8*(node2d-1)+l 472 call nident(ikboun,idof,nboun,id) 473 if(id.gt.0) then 474 if(ikboun(id).eq.idof) then 475 extnor(l,node)=0.d0 476 endif 477 endif 478 enddo 479 extnor(3,node)=0.d0 480! 481! else, all the directions with fixed 482! displacements are not considered 483! 484! elseif(lakon(nelem)(7:7).eq.' ') then 485! do l=1,3 486! if(nactdof(l,node).le.0) then 487! extnor(l,node)=0.d0 488! endif 489! enddo 490 endif 491 endif 492 enddo 493 else 494 do m=1,nopem 495 xi=xtri(1,m) 496 et=xtri(2,m) 497 call shape3tri(xi,et,xl2m,xsj2,xs2,shp2,iflag) 498 dd=dsqrt(xsj2(1)*xsj2(1) + xsj2(2)*xsj2(2) 499 & + xsj2(3)*xsj2(3)) 500 xsj2(1)=xsj2(1)/dd 501 xsj2(2)=xsj2(2)/dd 502 xsj2(3)=xsj2(3)/dd 503! 504 if(nope.eq.6) then 505 node=konl(ifacew1(m,jface)) 506 elseif(nope.eq.4) then 507 node=konl(ifacet(m,jface)) 508 endif 509 if((nodedesiinv(node).eq.0).or. 510 & ((nodedesiinv(node).eq.1).and. 511 & (nnodes.ge.nopedesi))) then 512 extnor(1,node)=extnor(1,node) 513 & +xsj2(1) 514 extnor(2,node)=extnor(2,node) 515 & +xsj2(2) 516 extnor(3,node)=extnor(3,node) 517 & +xsj2(3) 518! 519! in case of plain strain/stress/axi elements 520! not considering the x3-direction and the 521! directions with fixed displacements 522! 523 if((lakon(nelem)(7:7).eq.'A').or. 524 & (lakon(nelem)(7:7).eq.'S').or. 525 & (lakon(nelem)(7:7).eq.'E')) then 526 if(nope.eq.6) then 527 node2d=konl2d(ifacew1(m,jface)) 528 elseif(nope.eq.4) then 529 node2d=konl2d(ifacet(m,jface)) 530 endif 531 do l=1,2 532 idof=8*(node2d-1)+l 533 call nident(ikboun,idof,nboun,id) 534 if(id.gt.0) then 535 if(ikboun(id).eq.idof) then 536 extnor(l,node)=0.d0 537 endif 538 endif 539 enddo 540 extnor(3,node)=0.d0 541! 542! else, all the directions with fixed 543! displacements are not considered 544! 545! elseif(lakon(nelem)(7:7).eq.' ') then 546! do l=1,3 547! if(nactdof(l,node).le.0) then 548! extnor(l,node)=0.d0 549! endif 550! enddo 551 endif 552 endif 553 enddo 554 endif 555! 556 indexe=nodface(5,indexe) 557 if(indexe.eq.0) exit 558! 559 enddo 560 enddo 561! 562! normalizing the normals 563! 564 do l=1,nk 565 dd=dsqrt(extnor(1,l)**2+extnor(2,l)**2+ 566 & extnor(3,l)**2) 567 if(dd.gt.0.d0) then 568 do m=1,3 569 extnor(m,l)=extnor(m,l)/dd 570 enddo 571 endif 572 enddo 573! 574! in case of 2D elements all expanded nodes have to have the same 575! normal direction to get the correct normal direction for the 2D model 576! 577 if(ne2d.ne.0) then 578 do l=1,ndesi 579 nodemid=nodedesi(l) 580c write(*,*) 'nodemid',nodemid 581 if(iponod2dto3d(1,nodemid).ne.0) then 582 nodeboun1=iponod2dto3d(1,nodemid) 583 nodeboun2=iponod2dto3d(2,nodemid) 584 do m=1,3 585 extnor(m,nodeboun1)=extnor(m,nodemid) 586 extnor(m,nodeboun2)=extnor(m,nodemid) 587 enddo 588c write(*,*) 'nodeboun1',nodeboun1,nodeboun2 589 endif 590 enddo 591 endif 592! 593 return 594 end 595