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 getdesiinfo3d_robust(set,istartset,iendset,ialset, 20 & nset,mi,nactdof,ndesi,nodedesi,ntie,tieset,itmp,nmpc,nodempc, 21 & ipompc,nodedesiinv,iponoel,inoel,lakon,ipkon,kon,iregion, 22 & ipoface,nodface,nk,irandomtype) 23! 24! storing the design variables in nodedesi 25! marking which nodes are design variables in nodedesiinv 26! 27! a node is a design variable for the random field if: 28! 1) it belongs to the design variable set AND 29! 2) not all dofs in the node are defined by SPC's AND 30! 3) no MPC is applied to any of its dofs AND 31! 4) it belongs to at least one face whose number of 32! design variables exceeds half its nodes 33! 34 implicit none 35! 36 character*8 lakon(*) 37! 38 character*81 setname 39 character*81 set(*) 40 character*81 tieset(3,*) 41! 42 integer mi(*),istartset(*),iendset(*),ialset(*),ndesi, 43 & node,nodedesi(*),nset,ntie,i,j,k,l,m,nmpc,nodempc(3,*), 44 & nactdof(0:mi(2),*),itmp(*),ntmp,index,id,ipompc(*), 45 & nodedesiinv(*),iponoel(*),inoel(2,*),nelem,nope,nopedesi, 46 & ipkon(*),nnodes,kon(*),iregion,konl(26),iaux,kflag, 47 & ipoface(*),nodface(5,*),jfacem,nopesurf(9),ifaceq(8,6), 48 & ifacet(6,4),ifacew1(4,5),ifacew2(8,5),nopem,nk,ishift, 49 & indexe,jface,konl2d(26),expandhex(20),expandwed(15), 50 & idelta,idesi,irandomtype(*),iwrite 51! 52 setname(1:1)=' ' 53 ndesi=0 54! 55! nodes per face for hex elements 56! 57 data ifaceq /4,3,2,1,11,10,9,12, 58 & 5,6,7,8,13,14,15,16, 59 & 1,2,6,5,9,18,13,17, 60 & 2,3,7,6,10,19,14,18, 61 & 3,4,8,7,11,20,15,19, 62 & 4,1,5,8,12,17,16,20/ 63! 64! nodes per face for tet elements 65! 66 data ifacet /1,3,2,7,6,5, 67 & 1,2,4,5,9,8, 68 & 2,3,4,6,10,9, 69 & 1,4,3,8,10,7/ 70! 71! nodes per face for linear wedge elements 72! 73 data ifacew1 /1,3,2,0, 74 & 4,5,6,0, 75 & 1,2,5,4, 76 & 2,3,6,5, 77 & 3,1,4,6/ 78! 79! nodes per face for quadratic wedge elements 80! 81 data ifacew2 /1,3,2,9,8,7,0,0, 82 & 4,5,6,10,11,12,0,0, 83 & 1,2,5,4,7,14,10,13, 84 & 2,3,6,5,8,15,11,14, 85 & 3,1,4,6,9,13,12,15/ 86! 87! opening a file to store the nodes which are rejected as 88! design variables 89! 90 iwrite=0 91 open(40,file='WarnNodeDesignReject.nam',status='unknown') 92 write(40,*) '*NSET,NSET=WarnNodeDesignReject' 93! 94!------------------------------------------------------------------------ 95! assign the nodes of the set to the appropriate variables 96!------------------------------------------------------------------------ 97! 98 do i=1,nk 99 if(irandomtype(i).gt.0) then 100 ndesi=ndesi+1 101 nodedesi(ndesi)=i 102 nodedesiinv(i)=-2 103 endif 104 enddo 105! 106!------------------------------------------------------------------------ 107! Remove interior nodes from the set of design variables 108!------------------------------------------------------------------------ 109! 110 do j=1,nk 111! 112 if(ipoface(j).eq.0) cycle 113 indexe=ipoface(j) 114! 115 do 116! 117 nelem=nodface(3,indexe) 118 jface=nodface(4,indexe) 119! 120 if((lakon(nelem)(7:7).eq.'A').or. 121 & (lakon(nelem)(7:7).eq.'S').or. 122 & (lakon(nelem)(7:7).eq.'E')) then 123! 124! for plane stress/strain/axi only faces 125! 3 and higher are taken into account for the normal 126! 127 if(jface.le.2) then 128 indexe=nodface(5,indexe) 129 if(indexe.eq.0) then 130 exit 131 else 132 cycle 133 endif 134 endif 135 elseif(lakon(nelem)(7:7).eq.'L') then 136! 137! for shells only faces 2 and lower 138! taken into account for the normal 139! 140 if(jface.gt.2) then 141 indexe=nodface(5,indexe) 142 if(indexe.eq.0) then 143 exit 144 else 145 cycle 146 endif 147 endif 148 endif 149! 150! nopem: # of nodes in the surface 151! nope: # of nodes in the element 152! 153 if(lakon(nelem)(4:4).eq.'8') then 154 nopem=4 155 nope=8 156 nopedesi=3 157 ishift=8 158 elseif(lakon(nelem)(4:5).eq.'20') then 159 nopem=8 160 nope=20 161 nopedesi=5 162 ishift=20 163 elseif(lakon(nelem)(4:5).eq.'10') then 164 nopem=6 165 nope=10 166 nopedesi=4 167 elseif(lakon(nelem)(4:4).eq.'4') then 168 nopem=3 169 nope=4 170 nopedesi=3 171! 172! treatment of wedge faces 173! 174 elseif(lakon(nelem)(4:4).eq.'6') then 175 nope=6 176 if(jface.le.2) then 177 nopem=3 178 else 179 nopem=4 180 endif 181 nopedesi=3 182 ishift=6 183 elseif(lakon(nelem)(4:5).eq.'15') then 184 nope=15 185 if(jface.le.2) then 186 nopem=6 187 nopedesi=4 188 else 189 nopem=8 190 nopedesi=5 191 endif 192 ishift=15 193 endif 194 if(iregion.eq.0) then 195 nopedesi=0 196 endif 197! 198! actual position of the nodes belonging to the surface 199! 200 if((lakon(nelem)(7:7).eq.'A').or. 201 & (lakon(nelem)(7:7).eq.'S').or. 202 & (lakon(nelem)(7:7).eq.'E')) then 203 if((lakon(nelem)(4:5).eq.'20').or. 204 & (lakon(nelem)(4:5).eq.'8 ')) then 205 do k=1,nope 206 konl(k)=kon(ipkon(nelem)+k) 207 konl2d(k)=kon(ipkon(nelem)+ishift+expandhex(k)) 208 enddo 209 elseif((lakon(nelem)(4:5).eq.'15').or. 210 & (lakon(nelem)(4:5).eq.'6 ')) then 211 do k=1,nope 212 konl(k)=kon(ipkon(nelem)+k) 213 konl2d(k)=kon(ipkon(nelem)+ishift+expandwed(k)) 214 enddo 215 endif 216 else 217 do k=1,nope 218 konl(k)=kon(ipkon(nelem)+k) 219 enddo 220 endif 221! 222 if((nope.eq.20).or.(nope.eq.8)) then 223 do m=1,nopem 224 nopesurf(m)=konl(ifaceq(m,jface)) 225 enddo 226 elseif((nope.eq.10).or.(nope.eq.4)) then 227 do m=1,nopem 228 nopesurf(m)=konl(ifacet(m,jface)) 229 enddo 230 elseif(nope.eq.15) then 231 do m=1,nopem 232 nopesurf(m)=konl(ifacew2(m,jface)) 233 enddo 234 else 235 do m=1,nopem 236 nopesurf(m)=konl(ifacew1(m,jface)) 237 enddo 238 endif 239! 240! sum up how many designvariables are on that surface 241! 242 nnodes=0 243 do m=1,nopem 244 if(nodedesiinv(nopesurf(m)).ne.0) then 245 nnodes=nnodes+1 246 endif 247 enddo 248! 249! set design variables located on external surfaces 250! equal to 1 in nodedesiinv 251! if sufficient design variables are defined on this surface 252! 253 nnodes=0 254 do m=1,nopem 255 if((nodedesiinv(nopesurf(m)).lt.0).or. 256 & (nodedesiinv(nopesurf(m)).eq.1)) then 257 nnodes=nnodes+1 258 endif 259 enddo 260! 261 if(nnodes.ge.nopedesi) then 262 do m=1,nopem 263 if(nodedesiinv(nopesurf(m)).lt.0) then 264 nodedesiinv(nopesurf(m))=1 265 endif 266 enddo 267 else 268 do m=1,nopem 269 if(nodedesiinv(nopesurf(m)).eq.-2) then 270 nodedesiinv(nopesurf(m))=-1 271 endif 272 enddo 273 endif 274 275 indexe=nodface(5,indexe) 276 if(indexe.eq.0) exit 277! 278 enddo 279 enddo 280! 281 kflag=1 282 call isortii(nodedesi,iaux,ndesi,kflag) 283! 284!------------------------------------------------------------------------ 285! if node i in nodedesiinv(i) is -2 or -1 --> delete node i from 286! set of designvariables 287!------------------------------------------------------------------------ 288! 289 do i=1,nk 290 if(nodedesiinv(i).eq.-2) then 291! 292 write(*,*) '*WARNING in getdesiinfo:' 293 write(*,*) ' node ',i,' is removed' 294 write(*,*) ' from the set of design' 295 write(*,*) ' variables as it is an' 296 write(*,*) ' internal node' 297 write(40,*) i 298 iwrite=1 299! 300 elseif(nodedesiinv(i).eq.-1) then 301! 302 write(*,*) '*WARNING in getdesiinfo:' 303 write(*,*) ' node ',i,' is removed' 304 write(*,*) ' from the set of design' 305 write(*,*) ' variables as not sufficient ' 306 write(*,*) ' other variables are on the ' 307 write(*,*) ' surrounding element faces ' 308 write(40,*) i 309 iwrite=1 310! 311 nodedesiinv(i)=0 312 call nident(nodedesi,i,ndesi,id) 313 do k=id+1,ndesi 314 nodedesi(k-1)=nodedesi(k) 315 enddo 316 ndesi=ndesi-1 317 endif 318 enddo 319! 320!------------------------------------------------------------------------ 321! write remaining design variables in the dat file 322!------------------------------------------------------------------------ 323! 324 write(5,*)'NUMBER OF DESIGN VARIABLES TAKEN INTO ACCOUNT ' 325 write(5,*)'FOR RANDOM FIELD ASSESSMENT:' 326 write(5,*) ndesi 327 write(5,*)'' 328 write(5,*)'NODE NUMBERS OF DESIGN VARIABLES: ' 329 do i=1,ndesi 330 write(5,*) nodedesi(i),',' 331 enddo 332! 333 if(iwrite.eq.1) then 334 write(*,*) '*INFO in getdesiinfo3d_robust:' 335 write(*,*) ' rejected design nodes are stored in' 336 write(*,*) ' file WarnNodeDesignReject.nam' 337 write(*,*) ' This file can be loaded into' 338 write(*,*) ' an active cgx-session by typing' 339 write(*,*) 340 & ' read WarnNodeDesignReject.nam inp' 341 write(*,*) 342 close(40) 343 else 344 close(40,status='delete') 345 endif 346! 347 return 348 end 349 350