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(set,istartset,iendset,ialset,nset, 20 & mi,nactdof,ndesi,nodedesi,ntie,tieset,itmp,nmpc,nodempc, 21 & ipompc,nodedesiinv,iponoel,inoel,lakon,ipkon,kon,iregion, 22 & ipoface,nodface,nk) 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 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,iwrite 49! 50 setname(1:1)=' ' 51 ndesi=0 52! 53! nodes per face for hex elements 54! 55 data ifaceq /4,3,2,1,11,10,9,12, 56 & 5,6,7,8,13,14,15,16, 57 & 1,2,6,5,9,18,13,17, 58 & 2,3,7,6,10,19,14,18, 59 & 3,4,8,7,11,20,15,19, 60 & 4,1,5,8,12,17,16,20/ 61! 62! nodes per face for tet elements 63! 64 data ifacet /1,3,2,7,6,5, 65 & 1,2,4,5,9,8, 66 & 2,3,4,6,10,9, 67 & 1,4,3,8,10,7/ 68! 69! nodes per face for linear wedge elements 70! 71 data ifacew1 /1,3,2,0, 72 & 4,5,6,0, 73 & 1,2,5,4, 74 & 2,3,6,5, 75 & 3,1,4,6/ 76! 77! nodes per face for quadratic wedge elements 78! 79 data ifacew2 /1,3,2,9,8,7,0,0, 80 & 4,5,6,10,11,12,0,0, 81 & 1,2,5,4,7,14,10,13, 82 & 2,3,6,5,8,15,11,14, 83 & 3,1,4,6,9,13,12,15/ 84! 85! Search for the set name of the set with the design variables 86! 87 do i=1,ntie 88 if(tieset(1,i)(81:81).eq.'D') then 89 setname=tieset(2,i) 90 endif 91 enddo 92! 93! Check for the existence of the name 94! 95 if(setname(1:1).eq.' ') then 96 write(*,*) '*ERROR in getdesiinfo: name of node set ' 97 write(*,*) ' has not yet been defined. ' 98 call exit(201) 99 endif 100! 101! catalogue all nodes (dependent and independent) which 102! belong to MPC's and sort them in increasing order 103! 104 ntmp=0 105 do i=1,nmpc 106 index=ipompc(i) 107 do 108 if(index.eq.0) exit 109 node=nodempc(1,index) 110 call nident(itmp,node,ntmp,id) 111 if(id.gt.0) then 112 if(itmp(id).eq.node) then 113 index=nodempc(3,index) 114 cycle 115 endif 116 endif 117 ntmp=ntmp+1 118 do j=ntmp,id+2,-1 119 itmp(j)=itmp(j-1) 120 enddo 121 itmp(id+1)=node 122 index=nodempc(3,index) 123 enddo 124 enddo 125! 126! opening a file to store the nodes which are rejected as 127! design variables 128! 129 iwrite=0 130 open(40,file='WarnNodeDesignReject.nam',status='unknown') 131 write(40,*) '*NSET,NSET=WarnNodeDesignReject' 132! 133! Search the name of the node set in "set(i)" and 134! assign the nodes of the set to the appropriate variables 135! 136c do i=1,nset 137c if(setname.eq.set(i)) then 138 call cident81(set,setname,nset,i) 139 if(i.gt.0) then 140 if(setname.eq.set(i)) then 141 loop1: do j=istartset(i),iendset(i) 142 if(ialset(j).gt.0) then 143 node=ialset(j) 144! 145! check for SPC-constraints: if a node is constrained in 146! all dofs it is removed from the design node set 147! 148 do l=1,3 149 if(nactdof(l,node).gt.0) exit 150 if(l.eq.3) then 151 write(*,*) '*WARNING in getdesiinfo:' 152 write(*,*) ' node ',node,' has no' 153 write(*,*) ' active dofs and' 154 write(*,*) ' is removed from the set' 155 write(*,*) ' of design variables' 156 write(40,*) node 157 iwrite=1 158 cycle loop1 159 endif 160 enddo 161! 162! check for MPC-constraints 163! 164 call nident(itmp,node,ntmp,id) 165 if(id.gt.0) then 166 if(itmp(id).eq.node) then 167 write(*,*) '*WARNING in getdesiinfo:' 168 write(*,*) ' node ',node,' is subject' 169 write(*,*) ' to MPC-constraints and' 170 write(*,*) ' is removed from the set' 171 write(*,*) ' of design variables' 172 write(40,*) node 173 iwrite=1 174 cycle loop1 175 endif 176 endif 177! 178 ndesi=ndesi+1 179 nodedesi(ndesi)=node 180 else 181 k=ialset(j-2) 182 loop2: do 183 k=k-ialset(j) 184 if(k.ge.ialset(j-1)) exit 185! 186! check for SPC-constraints: if a node is constrained in 187! all dofs it is removed from the design node set 188! 189 do l=1,3 190 if(nactdof(l,k).gt.0) exit 191 if(l.eq.3) then 192 write(*,*) '*WARNING in getdesiinfo:' 193 write(*,*) ' node ',k,' has no' 194 write(*,*) ' active dofs and' 195 write(*,*) ' is removed from the set' 196 write(*,*) ' of design variables' 197 write(40,*) k 198 iwrite=1 199 cycle loop2 200 endif 201 enddo 202! 203! check for MPC-constraints 204! 205 call nident(itmp,k,ntmp,id) 206 if(id.gt.0) then 207 if(itmp(id).eq.k) then 208 write(*,*) '*WARNING in getdesiinfo:' 209 write(*,*) ' node ',k,' is subject' 210 write(*,*) ' to MPC-constraints and' 211 write(*,*) ' is removed from the set' 212 write(*,*) ' of design variables' 213 write(40,*) k 214 iwrite=1 215 cycle loop2 216 endif 217 endif 218! 219 ndesi=ndesi+1 220 nodedesi(ndesi)=k 221 enddo loop2 222 endif 223 enddo loop1 224 endif 225 endif 226c endif 227c enddo 228! 229! creating field nodedesiinv indicating for each node whether 230! it is a design variable or not 231! 232 do i=1,ndesi 233 index=nodedesi(i) 234 nodedesiinv(index)=-1 235 enddo 236! 237 kflag=1 238 call isortii(nodedesi,iaux,ndesi,kflag) 239! 240! A design node is also removed from nodedesi if it does not 241! belong to a face whose number of design variables exceeds half 242! of its nodes 243! 244 do i=1,nk 245 node=i 246 if(ipoface(node).eq.0) cycle 247 index=ipoface(node) 248 do 249 nelem=nodface(3,index) 250 jfacem=nodface(4,index) 251! 252 if(lakon(nelem)(4:4).eq.'8') then 253 nope=8 254 nopedesi=3 255 nopem=4 256 elseif(lakon(nelem)(4:5).eq.'20') then 257 nope=20 258 nopedesi=5 259 nopem=8 260 elseif(lakon(nelem)(4:5).eq.'10') then 261 nope=10 262 nopedesi=4 263 nopem=6 264 elseif(lakon(nelem)(4:4).eq.'4') then 265 nope=4 266 nopedesi=3 267 nopem=3 268 elseif(lakon(nelem)(4:4).eq.'6') then 269 nope=6 270 if(jfacem.le.2) then 271 nopem=3 272 nopedesi=3 273 else 274 nopem=4 275 nopedesi=3 276 endif 277 elseif(lakon(nelem)(4:5).eq.'15') then 278 nope=15 279 if(jfacem.le.2) then 280 nopem=6 281 nopedesi=4 282 else 283 nopem=8 284 nopedesi=5 285 endif 286 endif 287 if(iregion.eq.0) nopedesi=0 288! 289! actual position of the nodes belonging to the 290! master surface 291! 292 do k=1,nope 293 konl(k)=kon(ipkon(nelem)+k) 294 enddo 295! 296 if((nope.eq.20).or.(nope.eq.8)) then 297 do m=1,nopem 298 nopesurf(m)=konl(ifaceq(m,jfacem)) 299 enddo 300 elseif((nope.eq.10).or.(nope.eq.4)) then 301 do m=1,nopem 302 nopesurf(m)=konl(ifacet(m,jfacem)) 303 enddo 304 elseif(nope.eq.15) then 305 do m=1,nopem 306 nopesurf(m)=konl(ifacew2(m,jfacem)) 307 enddo 308 else 309 do m=1,nopem 310 nopesurf(m)=konl(ifacew1(m,jfacem)) 311 enddo 312 endif 313! 314! sum up how many designvariables are on that surface 315! 316 nnodes=0 317 do m=1,nopem 318 if(nodedesiinv(nopesurf(m)).ne.0) then 319 nnodes=nnodes+1 320 endif 321 enddo 322! 323 if(nnodes.ge.nopedesi) then 324 do m=1,nopem 325 if(nodedesiinv(nopesurf(m)).eq.-1) then 326 nodedesiinv(nopesurf(m))=1 327 endif 328 enddo 329 endif 330 index=nodface(5,index) 331 if(index.eq.0) exit 332 enddo 333 enddo 334! 335! if node i in nodedesi(i) is -1 --> delete node i from 336! set of designvariables 337! 338 do i=1,nk 339 if(nodedesiinv(i).eq.-1) then 340! 341 write(*,*) '*WARNING in getdesiinfo:' 342 write(*,*) ' node ',i,' is removed' 343 write(*,*) ' from the set of design' 344 write(*,*) ' variables as not sufficient ' 345 write(*,*) ' other variables are on the ' 346 write(*,*) ' surrounding element faces ' 347 write(40,*) i 348 iwrite=1 349! 350 nodedesiinv(i)=0 351 call nident(nodedesi,i,ndesi,id) 352 do k=id+1,ndesi 353 nodedesi(k-1)=nodedesi(k) 354 enddo 355 ndesi=ndesi-1 356 endif 357 enddo 358! 359 if(iwrite.eq.1) then 360 write(*,*) '*INFO in getdesiinfo:' 361 write(*,*) ' rejected design nodes are stored in' 362 write(*,*) ' file WarnNodeDesignReject.nam' 363 write(*,*) ' This file can be loaded into' 364 write(*,*) ' an active cgx-session by typing' 365 write(*,*) 366 & ' read WarnNodeDesignReject.nam inp' 367 write(*,*) 368 close(40) 369 else 370 close(40,status='delete') 371 endif 372! 373 return 374 end 375