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 findextsurface(nodface,ipoface,ne,ipkon,lakon, 20 & kon,konfa,ipkonfa,nk,lakonfa,nsurfs,ifreemax) 21! 22 implicit none 23! 24! 1) catalogues the external surface of the structure: see 25! comments further down 26! 2) creates topology fields for the external faces: 27! ipkonfa(1...nsurfs): pointer into konfa 28! konfa(..): konfa(ipkonfa(i)+1,....ipkonfa(i+1)) contains 29! the nodes belonging to face i 30! lakonfa(1...nsurfs): label (S3, S4, S6 or S8) 31! 32 character*8 lakon(*),lakonfa(*) 33! 34 integer ipkon(*),kon(*),ne,nodface(5,*),ipoface(*),nk,nopem,m, 35 & ithree,ifour,ifaceq(8,6),konfa(*),ipkonfa(*),nelemm,jfacem, 36 & ifacet(6,4),ifacew2(8,5),ifree,ifreenew,index,indexold, 37 & i,j,k,nodes(4),indexe,konl(26),nope,nsurfs, 38 & ifacew1(4,5),ifreemax 39! 40! 41! 42! nodes belonging to the element faces 43! 44 data ifaceq /4,3,2,1,11,10,9,12, 45 & 5,6,7,8,13,14,15,16, 46 & 1,2,6,5,9,18,13,17, 47 & 2,3,7,6,10,19,14,18, 48 & 3,4,8,7,11,20,15,19, 49 & 4,1,5,8,12,17,16,20/ 50 data ifacet /1,3,2,7,6,5, 51 & 1,2,4,5,9,8, 52 & 2,3,4,6,10,9, 53 & 1,4,3,8,10,7/ 54 data ifacew1 /1,3,2,0, 55 & 4,5,6,0, 56 & 1,2,5,4, 57 & 2,3,6,5, 58 & 3,1,4,6/ 59 data ifacew2 /1,3,2,9,8,7,0,0, 60 & 4,5,6,10,11,12,0,0, 61 & 1,2,5,4,7,14,10,13, 62 & 2,3,6,5,8,15,11,14, 63 & 3,1,4,6,9,13,12,15/ 64! 65 ithree=3 66 ifour=4 67! 68! determining the external element faces of the mesh; 69! the faces are catalogued by the three lowest nodes numbers 70! in ascending order. ipoface(i) points to a face for which 71! node i is the lowest node and nodface(1,ipoface(i)) and 72! nodface(2,ipoface(i)) are the next lower ones. 73! nodface(3,ipoface(i)) contains the element number, 74! nodface(4,ipoface(i)) the face number and nodface(5,ipoface(i)) 75! is a pointer to the next surface for which node i is the 76! lowest node; if there are no more such surfaces the pointer 77! has the value zero 78! An external element face is one which belongs to one element 79! only 80! 81 ifree=1 82 ifreemax=1 83 do i=1,6*ne-1 84 nodface(5,i)=i+1 85 enddo 86 do i=1,ne 87 if(ipkon(i).lt.0) cycle 88 if(lakon(i)(1:3).ne.'C3D') cycle 89 indexe=ipkon(i) 90 if((lakon(i)(4:4).eq.'2').or.(lakon(i)(4:4).eq.'8')) then 91 do j=1,6 92 do k=1,4 93 nodes(k)=kon(indexe+ifaceq(k,j)) 94 enddo 95 call insertsorti(nodes,ifour) 96c call isortii(nodes,iaux,ifour,kflag) 97 indexold=0 98 index=ipoface(nodes(1)) 99 do 100! 101! adding a surface which has not been 102! catalogued so far 103! 104 if(index.eq.0) then 105 ifreemax=max(ifreemax,ifree) 106 ifreenew=nodface(5,ifree) 107 nodface(1,ifree)=nodes(2) 108 nodface(2,ifree)=nodes(3) 109 nodface(3,ifree)=i 110 nodface(4,ifree)=j 111 nodface(5,ifree)=ipoface(nodes(1)) 112 ipoface(nodes(1))=ifree 113 ifree=ifreenew 114 exit 115 endif 116! 117! removing a surface which has already 118! been catalogued 119! 120 if((nodface(1,index).eq.nodes(2)).and. 121 & (nodface(2,index).eq.nodes(3))) then 122 if(indexold.eq.0) then 123 ipoface(nodes(1))=nodface(5,index) 124 else 125 nodface(5,indexold)=nodface(5,index) 126 endif 127 nodface(5,index)=ifree 128 ifree=index 129 exit 130 endif 131 indexold=index 132 index=nodface(5,index) 133 enddo 134 enddo 135 elseif((lakon(i)(4:4).eq.'4').or.(lakon(i)(4:5).eq.'10')) then 136 do j=1,4 137 do k=1,3 138 nodes(k)=kon(indexe+ifacet(k,j)) 139 enddo 140 call insertsorti(nodes,ithree) 141c call isortii(nodes,iaux,ithree,kflag) 142 indexold=0 143 index=ipoface(nodes(1)) 144 do 145! 146! adding a surface which has not been 147! catalogued so far 148! 149 if(index.eq.0) then 150 ifreemax=max(ifreemax,ifree) 151 ifreenew=nodface(5,ifree) 152 nodface(1,ifree)=nodes(2) 153 nodface(2,ifree)=nodes(3) 154 nodface(3,ifree)=i 155 nodface(4,ifree)=j 156 nodface(5,ifree)=ipoface(nodes(1)) 157 ipoface(nodes(1))=ifree 158 ifree=ifreenew 159 exit 160 endif 161! 162! removing a surface which has already 163! been catalogued 164! 165 if((nodface(1,index).eq.nodes(2)).and. 166 & (nodface(2,index).eq.nodes(3))) then 167 if(indexold.eq.0) then 168 ipoface(nodes(1))=nodface(5,index) 169 else 170 nodface(5,indexold)=nodface(5,index) 171 endif 172 nodface(5,index)=ifree 173 ifree=index 174 exit 175 endif 176 indexold=index 177 index=nodface(5,index) 178 enddo 179 enddo 180 else 181 do j=1,5 182 if(j.le.2) then 183 do k=1,3 184 nodes(k)=kon(indexe+ifacew2(k,j)) 185 enddo 186 call insertsorti(nodes,ithree) 187c call isortii(nodes,iaux,ithree,kflag) 188 else 189 do k=1,4 190 nodes(k)=kon(indexe+ifacew2(k,j)) 191 enddo 192 call insertsorti(nodes,ifour) 193c call isortii(nodes,iaux,ifour,kflag) 194 endif 195 indexold=0 196 index=ipoface(nodes(1)) 197 do 198! 199! adding a surface which has not been 200! catalogued so far 201! 202 if(index.eq.0) then 203 ifreemax=max(ifreemax,ifree) 204 ifreenew=nodface(5,ifree) 205 nodface(1,ifree)=nodes(2) 206 nodface(2,ifree)=nodes(3) 207 nodface(3,ifree)=i 208 nodface(4,ifree)=j 209 nodface(5,ifree)=ipoface(nodes(1)) 210 ipoface(nodes(1))=ifree 211 ifree=ifreenew 212 exit 213 endif 214! 215! removing a surface which has already 216! been catalogued 217! 218 if((nodface(1,index).eq.nodes(2)).and. 219 & (nodface(2,index).eq.nodes(3))) then 220 if(indexold.eq.0) then 221 ipoface(nodes(1))=nodface(5,index) 222 else 223 nodface(5,indexold)=nodface(5,index) 224 endif 225 nodface(5,index)=ifree 226 ifree=index 227 exit 228 endif 229 indexold=index 230 index=nodface(5,index) 231 enddo 232 enddo 233 endif 234 enddo 235! 236! Create fields konfa and ipkonfa for calculation of normals of shell 237! elements 238! 239 nsurfs=0 240 ifree=0 241 do j=1,nk 242! 243 if(ipoface(j).eq.0) cycle 244 indexe=ipoface(j) 245! 246 do 247 nsurfs=nsurfs+1 248 nelemm=nodface(3,indexe) 249 jfacem=nodface(4,indexe) 250! 251! treatment of hexahedral elements 252! 253 if(lakon(nelemm)(4:4).eq.'8') then 254 nopem=4 255 nope=8 256 elseif(lakon(nelemm)(4:5).eq.'20') then 257 nopem=8 258 nope=20 259! 260! treatment of tethrahedral elements 261! 262 elseif(lakon(nelemm)(4:5).eq.'10') then 263 nopem=6 264 nope=10 265 elseif(lakon(nelemm)(4:4).eq.'4') then 266 nopem=3 267 nope=4 268! 269! treatment of wedge faces 270! 271 elseif(lakon(nelemm)(4:4).eq.'6') then 272 nope=6 273 if(jfacem.le.2) then 274 nopem=3 275 else 276 nopem=4 277 endif 278 elseif(lakon(nelemm)(4:5).eq.'15') then 279 nope=15 280 if(jfacem.le.2) then 281 nopem=6 282 else 283 nopem=8 284 endif 285 endif 286! 287! actual position of the nodes 288! 289 do k=1,nope 290 konl(k)=kon(ipkon(nelemm)+k) 291 enddo 292! 293! quadratic quad shell element 294! 295 if((nope.eq.20).or. 296 & ((nope.eq.15).and.(jfacem.gt.2))) then 297 lakonfa(nsurfs)='S8' 298 ipkonfa(nsurfs)=ifree 299 do m=1,nopem 300 if(nope.eq.20) then 301 konfa(ifree+m)=konl(ifaceq(m,jfacem)) 302 else 303 konfa(ifree+m)=konl(ifacew2(m,jfacem)) 304 endif 305 enddo 306 ifree=ifree+nopem 307! 308! linear quad shell element 309! 310 elseif((nope.eq.8).or. 311 & ((nope.eq.6).and.(jfacem.gt.2))) then 312 lakonfa(nsurfs)='S4' 313 ipkonfa(nsurfs)=ifree 314 do m=1,nopem 315 if(nope.eq.8) then 316 konfa(ifree+m)=konl(ifaceq(m,jfacem)) 317 else 318 konfa(ifree+m)=konl(ifacew1(m,jfacem)) 319 endif 320 enddo 321 ifree=ifree+nopem 322! 323! quadratic tri shell element 324! 325 elseif((nope.eq.10).or. 326 & ((nope.eq.15).and.(jfacem.le.2))) then 327 lakonfa(nsurfs)='S6' 328 ipkonfa(nsurfs)=ifree 329 do m=1,nopem 330 if(nope.eq.10) then 331 konfa(ifree+m)=konl(ifacet(m,jfacem)) 332 else 333 konfa(ifree+m)=konl(ifacew2(m,jfacem)) 334 endif 335 enddo 336 ifree=ifree+nopem 337! 338! linear tri shell element 339! 340 elseif((nope.eq.4).or. 341 & ((nope.eq.6).and.(jfacem.le.2))) then 342 lakonfa(nsurfs)='S3' 343 ipkonfa(nsurfs)=ifree 344 do m=1,nopem 345 if(nope.eq.4) then 346 konfa(ifree+m)=konl(ifacet(m,jfacem)) 347 else 348 konfa(ifree+m)=konl(ifacew1(m,jfacem)) 349 endif 350 enddo 351 ifree=ifree+nopem 352 endif 353 indexe=nodface(5,indexe) 354 if(indexe.eq.0) exit 355 enddo 356 enddo 357! 358 return 359 end 360