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 distributesens(istartdesi,ialdesi,ipkon,lakon, 20 & ipoface,ndesi,nodedesi,nodface,kon,co,dgdx,nobject, 21 & weightformgrad,nodedesiinv,iregion,objectset, 22 & dgdxglob,nk,physcon,nobjectstart) 23! 24! calculation of the relationship between facial distributions 25! and corresponding nodal values 26! 27! transforming the nodal sensitivities into facially 28! distributed sensitivities 29! 30 implicit none 31! 32 character*81 objectset(5,*) 33 character*8 lakon(*) 34! 35 integer idesvar,j,k,kk,l,m,m1,n,istartdesi(*), 36 & ielem,ipoface(*),nodface(5,*),ndesi,nodedesi(*),nope, 37 & ialdesi(*),nopes,indexe,ifour,indexel, 38 & mint3d,mint2d,ipkon(*),konl(26),iflag,ifaceq(8,6), 39 & ifacet(6,4),ifacew1(4,5),ifacew2(8,5),kon(*), 40 & nodes1(8),ifacel,iobject,nodes2(8),indexs,nk, 41 & ithree,i,node,ifacq(2,3,20),ifact(2,3,10),ifacw(2,3,15), 42 & nopedesi,nnodes,nodedesiinv(*),nobject,iregion,iactnode, 43 & nobjectstart 44! 45 real*8 xi,et,weight,xl(3,9),xs(3,2),xsj(3),shp(7,9), 46 & co(3,*),xsjj,weightformgrad(ndesi),dgdx(ndesi,nobject), 47 & dgdxglob(2,nk,nobject),physcon(*) 48! 49! 50! 51! flag for shape functions 52! 53 data iflag /3/ 54 data indexel /0/ 55! 56 save indexel 57! 58 include "gauss.f" 59! 60 data ifaceq /4,3,2,1,11,10,9,12, 61 & 5,6,7,8,13,14,15,16, 62 & 1,2,6,5,9,18,13,17, 63 & 2,3,7,6,10,19,14,18, 64 & 3,4,8,7,11,20,15,19, 65 & 4,1,5,8,12,17,16,20/ 66! 67! nodes per face for tet elements 68! 69 data ifacet /1,3,2,7,6,5, 70 & 1,2,4,5,9,8, 71 & 2,3,4,6,10,9, 72 & 1,4,3,8,10,7/ 73! 74! nodes per face for linear wedge elements 75! 76 data ifacew1 /1,3,2,0, 77 & 4,5,6,0, 78 & 1,2,5,4, 79 & 2,3,6,5, 80 & 3,1,4,6/ 81! 82! nodes per face for quadratic wedge elements 83! 84 data ifacew2 /1,3,2,9,8,7,0,0, 85 & 4,5,6,10,11,12,0,0, 86 & 1,2,5,4,7,14,10,13, 87 & 2,3,6,5,8,15,11,14, 88 & 3,1,4,6,9,13,12,15/ 89! 90! ifacq indicates for each local node number to which 91! faces this node belongs and the internal number 92! within this face. e.g., local node number 1 belongs to 93! face 1 with internal number 4, face 3 with internal number 94! 1 and face 6 with internal number 2 (first line underneath) 95! 96! ifacq, ifact and ifacw can be derived in a unique way from 97! ifaceq, ifacet and ifacew2 98! 99 data ifacq /1,4,3,1,6,2, 100 &1,3,3,2,4,1, 101 &1,2,4,2,5,1, 102 &1,1,5,2,6,1, 103 &2,1,3,4,6,3, 104 &2,2,3,3,4,4, 105 &2,3,4,3,5,4, 106 &2,4,5,3,6,4, 107 &1,7,3,5,0,0, 108 &1,6,4,5,0,0, 109 &1,5,5,5,0,0, 110 &1,8,6,5,0,0, 111 &2,5,3,7,0,0, 112 &2,6,4,7,0,0, 113 &2,7,5,7,0,0, 114 &2,8,6,7,0,0, 115 &3,8,6,6,0,0, 116 &3,6,4,8,0,0, 117 &4,6,5,8,0,0, 118 &5,6,6,8,0,0/ 119! 120 data ifact /1,1,2,1,4,1, 121 &1,3,2,2,3,1, 122 &1,2,3,2,4,3, 123 &2,3,3,3,4,2, 124 &1,6,2,4,0,0, 125 &1,5,3,4,0,0, 126 &1,4,4,6,0,0, 127 &2,6,4,4,0,0, 128 &2,5,3,6,0,0, 129 &3,5,4,5,0,0/ 130! 131 data ifacw /1,1,3,1,5,2, 132 &1,3,3,2,4,1, 133 &1,2,4,2,5,1, 134 &2,1,3,4,5,3, 135 &2,2,3,3,4,4, 136 &2,3,4,3,5,4, 137 &1,6,3,5,0,0, 138 &1,5,4,5,0,0, 139 &1,4,5,5,0,0, 140 &2,4,3,7,0,0, 141 &2,5,4,7,0,0, 142 &2,6,5,7,0,0, 143 &3,8,5,6,0,0, 144 &3,6,4,8,0,0, 145 &4,6,5,8,0,0/ 146! 147! flag for shape functions 148! 149 ifour=4 150 ithree=3 151! 152! 153! Loop over all designvariables 154! 155 do idesvar=1,ndesi 156! 157 node=nodedesi(idesvar) 158! 159! Loop over all elements belonging to the designvariable 160! 161 do j=istartdesi(idesvar),istartdesi(idesvar+1)-1 162 ielem=ialdesi(j) 163 164 indexe=ipkon(ielem) 165! 166! mint2d: # of integration points on the surface 167! nopes: # of nodes in the surface 168! nope: # of nodes in the element 169! 170 if(lakon(ielem)(4:5).eq.'8R') then 171 mint2d=1 172 nopes=4 173 nope=8 174 nopedesi=3 175 elseif(lakon(ielem)(4:4).eq.'8') then 176 mint2d=4 177 nopes=4 178 nope=8 179 nopedesi=3 180 elseif(lakon(ielem)(4:5).eq.'20') then 181 mint2d=9 182 nopes=8 183 nope=20 184 nopedesi=5 185 elseif(lakon(ielem)(4:5).eq.'10') then 186 mint2d=3 187 nopes=6 188 nope=10 189c nopedesi=3 190 nopedesi=4 191 elseif(lakon(ielem)(4:4).eq.'4') then 192 mint2d=1 193 nopes=3 194 nope=4 195 nopedesi=3 196 elseif(lakon(ielem)(4:4).eq.'6') then 197 mint2d=1 198 nopes=3 199 nope=6 200 nopedesi=3 201 elseif(lakon(ielem)(4:5).eq.'15') then 202 mint2d=3 203 nopes=6 204 nope=15 205c nopedesi=3 206 else 207 exit 208 endif 209 if(iregion.eq.0) nopedesi=0 210! 211! actual position of the nodes belonging to the 212! master surface 213! 214 do l=1,nope 215 konl(l)=kon(indexe+l) 216 enddo 217! 218 do i=1,nope 219 if(kon(indexe+i).eq.node) exit 220 enddo 221! 222! Loop over all faces of the actual element 223! 224 loop1: do m=1,3 225! 226! hexahedral elements 227! 228 if((nope.eq.20).or.(nope.eq.8)) then 229 k=ifacq(1,m,i) 230 m1=ifacq(2,m,i) 231 do l=1,nopes 232 nodes1(l)=konl(ifaceq(l,k)) 233 nodes2(l)=nodes1(l) 234 enddo 235 call insertsorti(nodes1,ifour) 236c call isortii(nodes1,iaux,ifour,kflag) 237! 238! tetrahedral element 239! 240 elseif((nope.eq.10).or.(nope.eq.4)) then 241 k=ifact(1,m,i) 242 m1=ifact(2,m,i) 243 do l=1,nopes 244 nodes1(l)=konl(ifacet(l,k)) 245 nodes2(l)=nodes1(l) 246 enddo 247 call insertsorti(nodes1,ithree) 248c call isortii(nodes1,iaux,ithree,kflag) 249 nodes1(4)=0 250! 251! wedge element 252! 253 elseif(nope.eq.15) then 254 k=ifacw(1,m,i) 255 m1=ifacw(2,m,i) 256 if(k.le.2) then 257 nopes=6 258 nopedesi=4 259 mint3d=3 260 do l=1,nopes 261 nodes1(l)=konl(ifacew2(l,k)) 262 nodes2(l)=nodes1(l) 263 enddo 264 call insertsorti(nodes1,ithree) 265c call isortii(nodes1,iaux,ithree,kflag) 266 nodes1(4)=0 267 else 268 nopes=8 269 nopedesi=5 270 mint3d=4 271 do l=1,nopes 272 nodes1(l)=konl(ifacew2(l,k)) 273 nodes2(l)=nodes1(l) 274 enddo 275 call insertsorti(nodes1,ithree) 276c call isortii(nodes1,iaux,ithree,kflag) 277 endif 278 else 279 k=ifacw(1,m,i) 280 m1=ifacw(2,m,i) 281 if(k.le.2) then 282 nopes=3 283 mint3d=1 284 do l=1,nopes 285 nodes1(l)=konl(ifacew1(l,k)) 286 nodes2(l)=nodes1(l) 287 enddo 288 call insertsorti(nodes1,ithree) 289c call isortii(nodes1,iaux,ithree,kflag) 290 nodes1(4)=0 291 else 292 nopes=4 293 mint3d=2 294 do l=1,nopes 295 nodes1(l)=konl(ifacew1(l,k)) 296 nodes2(l)=nodes1(l) 297 enddo 298 call insertsorti(nodes1,ithree) 299c call isortii(nodes1,iaux,ithree,kflag) 300 endif 301 endif 302 if(k.eq.0) exit 303! 304! Find external element face 305! 306 indexs=ipoface(nodes1(1)) 307 do 308 if(indexs.eq.0) cycle loop1 309 if((nodface(1,indexs).eq.nodes1(2)).and. 310 & (nodface(2,indexs).eq.nodes1(3))) then 311! 312! Check if sufficient designnodes on that surface 313! 314 nnodes=0 315 do n=1,nopes 316 if(nodedesiinv(nodes1(n)).eq.1) then 317 nnodes=nnodes+1 318 endif 319 enddo 320! 321! Assign coordinates to nodes of surface 322! 323c write(*,*) 'formgradient',nnodes,nopedesi 324 if(nnodes.ge.nopedesi) then 325 if((nope.eq.20).or.(nope.eq.8)) then 326 do l=1,nopes 327 do n=1,3 328 ifacel=ifaceq(l,nodface(4,indexs)) 329 xl(n,l)=co(n,konl(ifacel)) 330 enddo 331 enddo 332 elseif((nope.eq.10).or.(nope.eq.4)) then 333 do l=1,nopes 334 do n=1,3 335 ifacel=ifacet(l,nodface(4,indexs)) 336 xl(n,l)=co(n,konl(ifacel)) 337 enddo 338 enddo 339 elseif(nope.eq.15) then 340 do l=1,nopes 341 do n=1,3 342 ifacel=ifacew2(l,nodface(4,indexs)) 343 xl(n,l)=co(n,konl(ifacel)) 344 enddo 345 enddo 346 elseif(nope.eq.6) then 347 do l=1,nopes 348 do n=1,3 349 ifacel=ifacew1(l,nodface(4,indexs)) 350 xl(n,l)=co(n,konl(ifacel)) 351 enddo 352 enddo 353 endif 354! 355! Determine weighting factors 356! 357 do kk=1,mint2d 358 if(lakon(ielem)(4:5).eq.'20') then 359 xi=gauss2d3(1,kk) 360 et=gauss2d3(2,kk) 361 weight=weight2d3(kk) 362 call shape8q(xi,et,xl,xsj, 363 & xs,shp,iflag) 364 elseif(lakon(ielem)(4:5).eq.'10') then 365 xi=gauss2d5(1,kk) 366 et=gauss2d5(2,kk) 367 weight=weight2d5(kk) 368 call shape6tri(xi,et,xl,xsj, 369 & xs,shp,iflag) 370 elseif(lakon(ielem)(4:4).eq.'4') then 371 xi=gauss2d4(1,kk) 372 et=gauss2d4(2,kk) 373 weight=weight2d4(kk) 374 call shape3tri(xi,et,xl,xsj, 375 & xs,shp,iflag) 376 elseif(lakon(ielem)(4:5).eq.'15') then 377 if(k.le.2) then 378 xi=gauss2d5(1,kk) 379 et=gauss2d5(2,kk) 380 weight=weight2d5(kk) 381 call shape6tri(xi,et,xl,xsj, 382 & xs,shp,iflag) 383 else 384 xi=gauss2d3(1,kk) 385 et=gauss2d3(2,kk) 386 weight=weight2d3(kk) 387 call shape8q(xi,et,xl,xsj, 388 & xs,shp,iflag) 389 endif 390 elseif(lakon(ielem)(4:5).eq.'8R') then 391 xi=gauss2d1(1,kk) 392 et=gauss2d1(2,kk) 393 weight=weight2d1(kk) 394 call shape4q(xi,et,xl,xsj, 395 & xs,shp,iflag) 396 elseif(lakon(ielem)(4:4).eq.'8') then 397 xi=gauss2d2(1,kk) 398 et=gauss2d2(2,kk) 399 weight=weight2d2(kk) 400 call shape4q(xi,et,xl,xsj, 401 & xs,shp,iflag) 402 elseif(lakon(ielem)(4:4).eq.'6') then 403 if(k.le.2) then 404 xi=gauss2d4(1,kk) 405 et=gauss2d4(2,kk) 406 weight=weight2d4(kk) 407 call shape3tri(xi,et,xl,xsj, 408 & xs,shp,iflag) 409 else 410 xi=gauss2d2(1,kk) 411 et=gauss2d2(2,kk) 412 weight=weight2d2(kk) 413 call shape4q(xi,et,xl,xsj, 414 & xs,shp,iflag) 415 endif 416 endif 417! 418! Calculate Jacobian determinant 419! 420 xsjj=dsqrt(xsj(1)**2+xsj(2)**2+ 421 & xsj(3)**2) 422! 423! Evaluate Shape functions for weightmatrix 424! 425 weightformgrad(idesvar)=weightformgrad 426 & (idesvar)+weight*shp(4,m1)*xsjj 427! write(*,*) 'formgradient ',idesvar,j,m,kk, 428! & weightformgrad(idesvar) 429 enddo 430 endif 431 exit 432 endif 433 indexs=nodface(5,indexs) 434 enddo ! find external element faces 435! 436 enddo loop1 ! Loop over all faces of the actual element 437! 438 enddo ! Loop over all elements belonging to the designvariable 439! 440 enddo ! Loop over all designvariables 441! 442! Change sign of sensitivity of objective function in case of 443! a maximization problem 444! 445 if(physcon(11).le.0.d0) then 446 if(objectset(2,1)(17:19).eq.'MAX') then 447 do idesvar=1,ndesi 448 dgdx(idesvar,1)=-dgdx(idesvar,1) 449 enddo 450 endif 451 endif 452! 453! Scaling of sensitivities 454! 455 do idesvar=1,ndesi 456 do iobject=1+nobjectstart,nobject 457 if(objectset(1,iobject)(4:13).eq.'MEMBERSIZE') cycle 458 if(objectset(1,iobject)(1:9).eq.'FIXGROWTH') cycle 459 if(objectset(1,iobject)(1:12).eq.'FIXSHRINKAGE') cycle 460 if(weightformgrad(idesvar).gt.0.d0) then 461 dgdx(idesvar,iobject)=dgdx(idesvar,iobject) 462 & /weightformgrad(idesvar) 463 iactnode=nodedesi(idesvar) 464 dgdxglob(1,iactnode,iobject)=dgdx(idesvar,iobject) 465 else 466 iactnode=nodedesi(idesvar) 467 dgdxglob(1,iactnode,iobject)=dgdx(idesvar,iobject) 468 endif 469 enddo 470 enddo 471! 472 return 473 end 474 475