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 objective_mass_dx(co,kon,ipkon,lakon,nelcon,rhcon, 20 & ielmat,ielorien,norien,ntmat_,matname,mi, 21 & thicke,mortar,nea,neb,ielprop,prop,distmin, 22 & ndesi,nodedesi,nobject,g0,dgdx,iobject,xmass, 23 & istartdesi,ialdesi,xdesi,idesvar) 24! 25! calculates the mass and its derivative w.r.t. the coordinates 26! of the mesh 27! 28 implicit none 29! 30 character*8 lakon(*),lakonl 31 character*80 amat,matname(*) 32! 33 integer kon(*),konl(26),nea,neb,mi(*),mint2d,nopes,nelcon(2,*), 34 & ielmat(mi(3),*),ielorien(mi(3),*),ntmat_,ipkon(*),iflag,null, 35 & mt,i,ii,j,k,jj,kk,indexe,nope,norien,ihyper,iactive, 36 & imat,mint3d,iorien,ilayer,nlayer,ki,kl,istartdesi(*), 37 & ielprop(*),mortar,idesvar,node,ndesi,nodedesi(*), 38 & nobject,iobject,ialdesi(*),ij,node1,node2, 39 & ifaceqexp(2,20),ifacewexp(2,15) 40! 41 real*8 co(3,*),shp(4,26),xl(3,26),vl(0:mi(2),26), 42 & prop(*),rhcon(0:1,ntmat_,*),xs2(3,7),thickness,rho,xl2(3,8), 43 & xsj2(3),shp2(7,8),xi,et,ze, 44 & xsj,weight,gs(8,4),a,tlayer(4),dlayer(4),xlayer(mi(3),4), 45 & thicke(mi(3),*),distmin,xmassel,g0(*),xmass(*),xdesi(3,*), 46 & dgdx(ndesi,nobject) 47! 48! 49! 50 include "gauss.f" 51! 52 iflag=3 53 null=0 54! 55 mt=mi(2)+1 56! 57! nodes in expansion direction of hex element 58! 59 data ifaceqexp /5,17, 60 & 6,18, 61 & 7,19, 62 & 8,20, 63 & 1,17, 64 & 2,18, 65 & 3,19, 66 & 4,20, 67 & 13,0, 68 & 14,0, 69 & 15,0, 70 & 16,0, 71 & 9,0, 72 & 10,0, 73 & 11,0, 74 & 12,0, 75 & 0,0, 76 & 0,0, 77 & 0,0, 78 & 0,0/ 79! 80! nodes in expansion direction of wedge element 81! 82 data ifacewexp /4,13, 83 & 5,14, 84 & 6,16, 85 & 1,13, 86 & 2,14, 87 & 3,15, 88 & 10,0, 89 & 11,0, 90 & 12,0, 91 & 7,0, 92 & 8,0, 93 & 9,0, 94 & 0,0, 95 & 0,0, 96 & 0,0/ 97! 98! ------------------------------------------------------------- 99! Initialisation of the loop for the variation of 100! the internal forces 101! ------------------------------------------------------------- 102! 103! 104! Loop over all elements in thread 105 do ij=nea,neb 106 if(idesvar.gt.0) then 107 i=ialdesi(ij) 108 else 109 i=ij 110 endif 111 112! initialisation of mass 113! 114 xmassel=0.d0 115 lakonl=lakon(i) 116! 117 if(ipkon(i).lt.0) cycle 118! 119! no 3D-fluid elements 120! 121 if(lakonl(1:1).eq.'F') cycle 122 if(lakonl(1:7).eq.'DCOUP3D') cycle 123! 124 if(lakonl(7:8).ne.'LC') then 125! 126 imat=ielmat(1,i) 127 amat=matname(imat) 128 if(norien.gt.0) then 129 iorien=max(0,ielorien(1,i)) 130 else 131 iorien=0 132 endif 133! 134 if(nelcon(1,imat).lt.0) then 135 ihyper=1 136 else 137 ihyper=0 138 endif 139 else 140! 141! composite materials 142! 143! determining the number of layers 144! 145 nlayer=0 146 do k=1,mi(3) 147 if(ielmat(k,i).ne.0) then 148 nlayer=nlayer+1 149 endif 150 enddo 151! 152! determining the layer thickness and global thickness 153! at the shell integration points 154! 155 if(lakonl(4:5).eq.'20') then 156! 157 mint2d=4 158 nopes=8 159! 160 iflag=1 161 indexe=ipkon(i) 162 do kk=1,mint2d 163 xi=gauss3d2(1,kk) 164 et=gauss3d2(2,kk) 165 call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag) 166 tlayer(kk)=0.d0 167 do k=1,nlayer 168 thickness=0.d0 169 do j=1,nopes 170 thickness=thickness+thicke(k,indexe+j)* 171 & shp2(4,j) 172 enddo 173 tlayer(kk)=tlayer(kk)+thickness 174 xlayer(k,kk)=thickness 175 enddo 176 enddo 177 iflag=3 178! 179 ilayer=0 180 do k=1,4 181 dlayer(k)=0.d0 182 enddo 183 elseif(lakonl(4:5).eq.'15') then 184! 185 mint2d=3 186 nopes=6 187! 188 iflag=1 189 indexe=ipkon(i) 190 do kk=1,mint2d 191 xi=gauss3d10(1,kk) 192 et=gauss3d10(2,kk) 193 call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 194 tlayer(kk)=0.d0 195 do k=1,nlayer 196 thickness=0.d0 197 do j=1,nopes 198 thickness=thickness+thicke(k,indexe+j)* 199 & shp2(4,j) 200 enddo 201 tlayer(kk)=tlayer(kk)+thickness 202 xlayer(k,kk)=thickness 203 enddo 204 enddo 205 iflag=3 206! 207 ilayer=0 208 do k=1,3 209 dlayer(k)=0.d0 210 enddo 211 endif 212! 213 endif 214! 215 indexe=ipkon(i) 216c Bernhardi start 217 if(lakonl(1:5).eq.'C3D8I') then 218 nope=11 219 elseif(lakonl(4:5).eq.'20') then 220c Bernhardi end 221 nope=20 222 elseif(lakonl(4:4).eq.'2') then 223 nope=26 224 elseif(lakonl(4:4).eq.'8') then 225 nope=8 226 elseif(lakonl(4:5).eq.'10') then 227 nope=10 228 elseif(lakonl(4:5).eq.'14') then 229 nope=14 230 elseif(lakonl(4:4).eq.'4') then 231 nope=4 232 elseif(lakonl(4:5).eq.'15') then 233 nope=15 234 elseif(lakonl(4:4).eq.'6') then 235 nope=6 236 elseif((lakonl(1:1).eq.'E').and.(lakonl(7:7).ne.'F')) then 237! 238! spring elements, contact spring elements and 239! dashpot elements 240! 241 if(lakonl(7:7).eq.'C') then 242! 243! contact spring elements 244! 245 if(mortar.eq.1) then 246! 247! face-to-face penalty 248! 249 nope=kon(ipkon(i)) 250 elseif(mortar.eq.0) then 251! 252! node-to-face penalty 253! 254 nope=ichar(lakonl(8:8))-47 255 konl(nope+1)=kon(indexe+nope+1) 256 endif 257 else 258! 259! genuine spring elements and dashpot elements 260! 261 nope=ichar(lakonl(8:8))-47 262 endif 263 else 264 cycle 265 endif 266! 267! check whether a design variable belongs to the element 268! 269 if(idesvar.gt.0) then 270 do ii=1,nope 271 node=kon(indexe+ii) 272 if(node.eq.nodedesi(idesvar)) then 273 iactive=ii 274 exit 275 endif 276 enddo 277 endif 278! 279 if(lakonl(4:5).eq.'8R') then 280 mint3d=1 281 elseif(lakonl(4:7).eq.'20RB') then 282 if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then 283 mint3d=50 284 else 285 call beamintscheme(lakonl,mint3d,ielprop(i),prop, 286 & null,xi,et,ze,weight) 287 endif 288 elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'26R').or. 289 & (lakonl(4:6).eq.'20R')) then 290 if(lakonl(7:8).eq.'LC') then 291 mint3d=8*nlayer 292 else 293 mint3d=8 294 endif 295 elseif(lakonl(4:4).eq.'2') then 296 mint3d=27 297 elseif((lakonl(4:5).eq.'10').or.(lakonl(4:5).eq.'14')) then 298 mint3d=4 299 elseif(lakonl(4:4).eq.'4') then 300 mint3d=1 301 elseif(lakonl(4:5).eq.'15') then 302 if(lakonl(7:8).eq.'LC') then 303 mint3d=6*nlayer 304 else 305 mint3d=9 306 endif 307 elseif(lakonl(4:4).eq.'6') then 308 mint3d=2 309 elseif(lakonl(1:1).eq.'E') then 310 mint3d=0 311 endif 312 313 do j=1,nope 314 konl(j)=kon(indexe+j) 315 do k=1,3 316 xl(k,j)=co(k,konl(j)) 317 enddo 318 enddo 319! 320! computation of the objective function and the derivate 321! if the designnode belongs to the considered element 322! if not, next element in loop will be taken 323! 324 if(idesvar.gt.0) then 325 do j=1,3 326 xl(j,iactive)=xl(j,iactive)+xdesi(j,idesvar) 327 enddo 328 if(lakonl(1:5).eq.'C3D20') then 329 if((lakonl(7:7).eq.'A').or. 330 & (lakonl(7:7).eq.'L').or. 331 & (lakonl(7:7).eq.'S').or. 332 & (lakonl(7:7).eq.'E')) then 333 node1=ifaceqexp(1,iactive) 334 node2=ifaceqexp(2,iactive) 335 do j=1,3 336 if(iactive.le.8) then 337 xl(j,node1)=xl(j,node1)+xdesi(j,idesvar) 338 xl(j,node2)=xl(j,node2)+xdesi(j,idesvar) 339 else 340 xl(j,node1)=xl(j,node1)+xdesi(j,idesvar) 341 endif 342 enddo 343 endif 344 elseif(lakonl(1:5).eq.'C3D15') then 345 if((lakonl(7:7).eq.'A').or. 346 & (lakonl(7:7).eq.'L').or. 347 & (lakonl(7:7).eq.'S').or. 348 & (lakonl(7:7).eq.'E')) then 349 node1=ifacewexp(1,iactive) 350 node2=ifacewexp(2,iactive) 351 do j=1,3 352 if(iactive.le.6) then 353 xl(j,node1)=xl(j,node1)+xdesi(j,idesvar) 354 xl(j,node2)=xl(j,node2)+xdesi(j,idesvar) 355 else 356 xl(j,node1)=xl(j,node1)+xdesi(j,idesvar) 357 endif 358 enddo 359 endif 360 elseif(lakonl(1:4).eq.'C3D8') then 361 if((lakonl(7:7).eq.'A').or. 362 & (lakonl(7:7).eq.'L').or. 363 & (lakonl(7:7).eq.'S').or. 364 & (lakonl(7:7).eq.'E')) then 365 node1=ifaceqexp(1,iactive) 366 do j=1,3 367 xl(j,node1)=xl(j,node1)+xdesi(j,idesvar) 368 enddo 369 endif 370 elseif(lakonl(1:4).eq.'C3D6') then 371 if((lakonl(7:7).eq.'A').or. 372 & (lakonl(7:7).eq.'L').or. 373 & (lakonl(7:7).eq.'S').or. 374 & (lakonl(7:7).eq.'E')) then 375 node1=ifaceqexp(1,iactive) 376 do j=1,3 377 xl(j,node1)=xl(j,node1)+xdesi(j,idesvar) 378 enddo 379 endif 380 endif 381 endif 382! 383 do jj=1,mint3d 384 if(lakonl(4:5).eq.'8R') then 385 xi=gauss3d1(1,jj) 386 et=gauss3d1(2,jj) 387 ze=gauss3d1(3,jj) 388 weight=weight3d1(jj) 389 elseif(lakonl(4:7).eq.'20RB') then 390 if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then 391 xi=gauss3d13(1,jj) 392 et=gauss3d13(2,jj) 393 ze=gauss3d13(3,jj) 394 weight=weight3d13(jj) 395 else 396 call beamintscheme(lakonl,mint3d,ielprop(i),prop, 397 & jj,xi,et,ze,weight) 398 endif 399 elseif((lakonl(4:4).eq.'8').or. 400 & (lakonl(4:6).eq.'20R').or.(lakonl(4:6).eq.'26R')) 401 & then 402 if(lakonl(7:8).ne.'LC') then 403 xi=gauss3d2(1,jj) 404 et=gauss3d2(2,jj) 405 ze=gauss3d2(3,jj) 406 weight=weight3d2(jj) 407 else 408 kl=mod(jj,8) 409 if(kl.eq.0) kl=8 410! 411 xi=gauss3d2(1,kl) 412 et=gauss3d2(2,kl) 413 ze=gauss3d2(3,kl) 414 weight=weight3d2(kl) 415! 416 ki=mod(jj,4) 417 if(ki.eq.0) ki=4 418! 419 if(kl.eq.1) then 420 ilayer=ilayer+1 421 if(ilayer.gt.1) then 422 do k=1,4 423 dlayer(k)=dlayer(k)+xlayer(ilayer-1,k) 424 enddo 425 endif 426 endif 427 ze=2.d0*(dlayer(ki)+(ze+1.d0)/ 428 & 2.d0*xlayer(ilayer,ki))/tlayer(ki)-1.d0 429 weight=weight*xlayer(ilayer,ki)/tlayer(ki) 430! 431! material and orientation 432! 433 imat=ielmat(ilayer,i) 434 amat=matname(imat) 435 if(norien.gt.0) then 436 iorien=max(0,ielorien(ilayer,i)) 437 else 438 iorien=0 439 endif 440! 441 if(nelcon(1,imat).lt.0) then 442 ihyper=1 443 else 444 ihyper=0 445 endif 446 endif 447 elseif(lakonl(4:4).eq.'2') then 448 xi=gauss3d3(1,jj) 449 et=gauss3d3(2,jj) 450 ze=gauss3d3(3,jj) 451 weight=weight3d3(jj) 452 elseif((lakonl(4:5).eq.'10').or.(lakonl(4:5).eq.'14')) 453 & then 454 xi=gauss3d5(1,jj) 455 et=gauss3d5(2,jj) 456 ze=gauss3d5(3,jj) 457 weight=weight3d5(jj) 458 elseif(lakonl(4:4).eq.'4') then 459 xi=gauss3d4(1,jj) 460 et=gauss3d4(2,jj) 461 ze=gauss3d4(3,jj) 462 weight=weight3d4(jj) 463 elseif(lakonl(4:5).eq.'15') then 464 if(lakonl(7:8).ne.'LC') then 465 xi=gauss3d8(1,jj) 466 et=gauss3d8(2,jj) 467 ze=gauss3d8(3,jj) 468 weight=weight3d8(jj) 469 else 470 kl=mod(jj,6) 471 if(kl.eq.0) kl=6 472! 473 xi=gauss3d10(1,kl) 474 et=gauss3d10(2,kl) 475 ze=gauss3d10(3,kl) 476 weight=weight3d10(kl) 477! 478 ki=mod(jj,3) 479 if(ki.eq.0) ki=3 480! 481 if(kl.eq.1) then 482 ilayer=ilayer+1 483 if(ilayer.gt.1) then 484 do k=1,3 485 dlayer(k)=dlayer(k)+xlayer(ilayer-1,k) 486 enddo 487 endif 488 endif 489 ze=2.d0*(dlayer(ki)+(ze+1.d0)/ 490 & 2.d0*xlayer(ilayer,ki))/tlayer(ki)-1.d0 491 weight=weight*xlayer(ilayer,ki)/tlayer(ki) 492! 493! material and orientation 494! 495 imat=ielmat(ilayer,i) 496 amat=matname(imat) 497 if(norien.gt.0) then 498 iorien=max(0,ielorien(ilayer,i)) 499 else 500 iorien=0 501 endif 502! 503 if(nelcon(1,imat).lt.0) then 504 ihyper=1 505 else 506 ihyper=0 507 endif 508 endif 509 elseif(lakonl(4:4).eq.'6') then 510 xi=gauss3d7(1,jj) 511 et=gauss3d7(2,jj) 512 ze=gauss3d7(3,jj) 513 weight=weight3d7(jj) 514 endif 515! 516c Bernhardi start 517 if(lakonl(1:5).eq.'C3D8R') then 518 call shape8hr(xl,xsj,shp,gs,a) 519 elseif(lakonl(1:5).eq.'C3D8I') then 520 call shape8hu(xi,et,ze,xl,xsj,shp,iflag) 521 elseif(nope.eq.20) then 522c Bernhardi end 523 if(lakonl(7:7).eq.'A') then 524 call shape20h_ax(xi,et,ze,xl,xsj,shp,iflag) 525 elseif((lakonl(7:7).eq.'E').or. 526 & (lakonl(7:7).eq.'S')) then 527 call shape20h_pl(xi,et,ze,xl,xsj,shp,iflag) 528 else 529 call shape20h(xi,et,ze,xl,xsj,shp,iflag) 530 endif 531c elseif(nope.eq.26) then 532c call shape26h(xi,et,ze,xl,xsj,shp,iflag,konl) 533 elseif(nope.eq.8) then 534 call shape8h(xi,et,ze,xl,xsj,shp,iflag) 535 elseif(nope.eq.10) then 536 call shape10tet(xi,et,ze,xl,xsj,shp,iflag) 537c elseif(nope.eq.14) then 538c call shape14tet(xi,et,ze,xl,xsj,shp,iflag,konl) 539 elseif(nope.eq.4) then 540 call shape4tet(xi,et,ze,xl,xsj,shp,iflag) 541 elseif(nope.eq.15) then 542 call shape15w(xi,et,ze,xl,xsj,shp,iflag) 543 else 544 call shape6w(xi,et,ze,xl,xsj,shp,iflag) 545 endif 546! 547! calculation of the objective function 548! 549 rho=rhcon(1,1,imat) 550 xmassel=xmassel+weight*xsj*rho 551! 552! write(5,'(i5,3x,i5,3x,i5,3x,e14.8,3x,e14.7,3x,e14.8, 553! & 3x,e14.8,3x,e14.8,3x,e20.14)') 554! & idesvar,i,jj,weight,xsj,xi,et,ze, 555! & xmassel 556! 557! end of loop over all integration points 558 enddo 559! 560! summation of the objective function over all elements 561! 562 if(idesvar.eq.0) then 563 xmass(i)=xmassel 564 g0(iobject)=g0(iobject)+xmassel 565 else 566 dgdx(idesvar,iobject)=dgdx(idesvar,iobject)+ 567 & (xmassel-xmass(i))/distmin 568 endif 569! 570! end of loop over all elements in thread 571! 572 enddo 573! 574! 575 return 576 end 577 578