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 createtiedsurfs(nodface,ipoface,set,istartset, 20 & iendset,ialset,tieset,inomat,ne,ipkon,lakon,kon,ntie, 21 & tietol,nalset,nk,nset,iactive) 22! 23 implicit none 24! 25! creates ties for the surfaces in between domains of an 26! electromagnetic calculation 27! 28 character*8 lakon(*) 29 character*81 set(*),tieset(3,*),surfset 30! 31 integer ipkon(*),kon(*),ne,nodface(5,*),ipoface(*),istartset(*), 32 & iendset(*),ialset(*),inomat(*),ithree,ifour,ifaceq(8,6),iset, 33 & ifacet(6,4),ifacew(8,5),ifree,ifreenew,index,indexold,id, 34 & i,j,k,iactive(3),ntie,nodes(4),nalset,nk,nset,indexe 35! 36 real*8 tietol(3,*) 37! 38! nodes belonging to the element faces 39! 40 data ifaceq /4,3,2,1,11,10,9,12, 41 & 5,6,7,8,13,14,15,16, 42 & 1,2,6,5,9,18,13,17, 43 & 2,3,7,6,10,19,14,18, 44 & 3,4,8,7,11,20,15,19, 45 & 4,1,5,8,12,17,16,20/ 46 data ifacet /1,3,2,7,6,5, 47 & 1,2,4,5,9,8, 48 & 2,3,4,6,10,9, 49 & 1,4,3,8,10,7/ 50 data ifacew /1,3,2,9,8,7,0,0, 51 & 4,5,6,10,11,12,0,0, 52 & 1,2,5,4,7,14,10,13, 53 & 2,3,6,5,8,15,11,14, 54 & 4,6,3,1,12,15,9,13/ 55! 56 ithree=3 57 ifour=4 58! 59! determining the external element faces of the electromagnetic mesh 60! the faces are catalogued by the three lowes nodes numbers 61! in ascending order. ipoface(i) points to a face for which 62! node i is the lowest node and nodface(1,ipoface(i)) and 63! nodface(2,ipoface(i)) are the next lower ones. 64! nodface(3,ipoface(i)) contains the element number, 65! nodface(4,ipoface(i)) the face number and nodface(5,ipoface(i)) 66! is a pointer to the next surface for which node i is the 67! lowest node; if there are no more such surfaces the pointer 68! has the value zero 69! An external element face is one which belongs to one element 70! only 71! 72 ifree=1 73 do i=1,6*ne-1 74 nodface(5,i)=i+1 75 enddo 76 do i=1,ne 77 if(ipkon(i).lt.0) cycle 78 if(lakon(i)(1:3).ne.'C3D') cycle 79 indexe=ipkon(i) 80 if((lakon(i)(4:4).eq.'2').or.(lakon(i)(4:4).eq.'8')) then 81 do j=1,6 82 do k=1,4 83 nodes(k)=kon(indexe+ifaceq(k,j)) 84 enddo 85 call insertsorti(nodes,ifour) 86c call isortii(nodes,iaux,ifour,kflag) 87 indexold=0 88 index=ipoface(nodes(1)) 89 do 90! 91! adding a surface which has not been 92! catalogued so far 93! 94 if(index.eq.0) then 95 ifreenew=nodface(5,ifree) 96 nodface(1,ifree)=nodes(2) 97 nodface(2,ifree)=nodes(3) 98 nodface(3,ifree)=i 99 nodface(4,ifree)=j 100 nodface(5,ifree)=ipoface(nodes(1)) 101 ipoface(nodes(1))=ifree 102 ifree=ifreenew 103 exit 104 endif 105! 106! removing a surface which has already 107! been catalogued 108! 109 if((nodface(1,index).eq.nodes(2)).and. 110 & (nodface(2,index).eq.nodes(3))) then 111 if(indexold.eq.0) then 112 ipoface(nodes(1))=nodface(5,index) 113 else 114 nodface(5,indexold)=nodface(5,index) 115 endif 116 nodface(5,index)=ifree 117 ifree=index 118 exit 119 endif 120 indexold=index 121 index=nodface(5,index) 122 enddo 123 enddo 124 elseif((lakon(i)(4:4).eq.'4').or.(lakon(i)(4:5).eq.'10')) then 125 do j=1,4 126 do k=1,3 127 nodes(k)=kon(indexe+ifacet(k,j)) 128 enddo 129 call insertsorti(nodes,ithree) 130c call isortii(nodes,iaux,ithree,kflag) 131 indexold=0 132 index=ipoface(nodes(1)) 133 do 134! 135! adding a surface which has not been 136! catalogued so far 137! 138 if(index.eq.0) then 139 ifreenew=nodface(5,ifree) 140 nodface(1,ifree)=nodes(2) 141 nodface(2,ifree)=nodes(3) 142 nodface(3,ifree)=i 143 nodface(4,ifree)=j 144 nodface(5,ifree)=ipoface(nodes(1)) 145 ipoface(nodes(1))=ifree 146 ifree=ifreenew 147 exit 148 endif 149! 150! removing a surface which has already 151! been catalogued 152! 153 if((nodface(1,index).eq.nodes(2)).and. 154 & (nodface(2,index).eq.nodes(3))) then 155 if(indexold.eq.0) then 156 ipoface(nodes(1))=nodface(5,index) 157 else 158 nodface(5,indexold)=nodface(5,index) 159 endif 160 nodface(5,index)=ifree 161 ifree=index 162 exit 163 endif 164 indexold=index 165 index=nodface(5,index) 166 enddo 167 enddo 168 else 169 do j=1,5 170 if(j.le.2) then 171 do k=1,3 172 nodes(k)=kon(indexe+ifacew(k,j)) 173 enddo 174 call insertsorti(nodes,ithree) 175c call isortii(nodes,iaux,ithree,kflag) 176 else 177 do k=1,4 178 nodes(k)=kon(indexe+ifacew(k,j)) 179 enddo 180 call insertsorti(nodes,ifour) 181c call isortii(nodes,iaux,ifour,kflag) 182 endif 183 indexold=0 184 index=ipoface(nodes(1)) 185 do 186! 187! adding a surface which has not been 188! catalogued so far 189! 190 if(index.eq.0) then 191 ifreenew=nodface(5,ifree) 192 nodface(1,ifree)=nodes(2) 193 nodface(2,ifree)=nodes(3) 194 nodface(3,ifree)=i 195 nodface(4,ifree)=j 196 nodface(5,ifree)=ipoface(nodes(1)) 197 ipoface(nodes(1))=ifree 198 ifree=ifreenew 199 exit 200 endif 201! 202! removing a surface which has already 203! been catalogued 204! 205 if((nodface(1,index).eq.nodes(2)).and. 206 & (nodface(2,index).eq.nodes(3))) then 207 if(indexold.eq.0) then 208 ipoface(nodes(1))=nodface(5,index) 209 else 210 nodface(5,indexold)=nodface(5,index) 211 endif 212 nodface(5,index)=ifree 213 ifree=index 214 exit 215 endif 216 indexold=index 217 index=nodface(5,index) 218 enddo 219 enddo 220 endif 221 enddo 222! 223! boundary of domain 1 (phi-domain) 224! 225 surfset(1:22)='ELECTROMAGNETICSZONE1T' 226 do i=23,81 227 surfset(i:i)=' ' 228 enddo 229 call cident81(set,surfset,nset,id) 230 nset=nset+1 231 do j=nset,id+2,-1 232 istartset(j)=istartset(j-1) 233 iendset(j)=iendset(j-1) 234 set(j)=set(j-1) 235 enddo 236 iset=id+1 237ccc to remove start 238c iset=nset 239ccc to remove end 240 istartset(iset)=nalset+1 241 do i=1,nk 242 if(ipoface(i).eq.0) cycle 243 if(inomat(i).ne.1) cycle 244 index=ipoface(i) 245 do 246 if(index.eq.0) exit 247 nalset=nalset+1 248 ialset(nalset)=10*nodface(3,index)+nodface(4,index) 249 index=nodface(5,index) 250 enddo 251 enddo 252 if(istartset(iset).gt.nalset) then 253 iactive(1)=0 254 do j=iset+1,nset 255 istartset(j-1)=istartset(j) 256 iendset(j-1)=iendset(j) 257 set(j-1)=set(j) 258 enddo 259 nset=nset-1 260 else 261 iactive(1)=iset 262ccc to remove start 263c set(iset)(1:22)='ELECTROMAGNETICSZONE1T' 264c do i=23,81 265c set(iset)(i:i)=' ' 266c enddo 267ccc to remove end 268 set(iset)=surfset 269 iendset(iset)=nalset 270 endif 271! 272! boundary of domain 2 (A-V-domain) 273! 274 surfset(1:22)='ELECTROMAGNETICSZONE2T' 275 do i=23,81 276 surfset(i:i)=' ' 277 enddo 278 call cident81(set,surfset,nset,id) 279 nset=nset+1 280 do j=nset,id+2,-1 281 istartset(j)=istartset(j-1) 282 iendset(j)=iendset(j-1) 283 set(j)=set(j-1) 284 enddo 285 iset=id+1 286ccc to remove start 287c iset=nset 288ccc to remove end 289 istartset(iset)=nalset+1 290 do i=1,nk 291 if(ipoface(i).eq.0) cycle 292 if(inomat(i).ne.2) cycle 293 index=ipoface(i) 294 do 295 if(index.eq.0) exit 296 nalset=nalset+1 297 ialset(nalset)=10*nodface(3,index)+nodface(4,index) 298 index=nodface(5,index) 299 enddo 300 enddo 301 if(istartset(iset).gt.nalset) then 302 iactive(2)=0 303 do j=iset+1,nset 304 istartset(j-1)=istartset(j) 305 iendset(j-1)=iendset(j) 306 set(j-1)=set(j) 307 enddo 308 nset=nset-1 309 else 310 iactive(2)=iset 311ccc to remove start 312c set(iset)(1:22)='ELECTROMAGNETICSZONE2T' 313c do i=23,81 314c set(iset)(i:i)=' ' 315c enddo 316ccc to remove end 317 set(iset)=surfset 318 iendset(iset)=nalset 319 endif 320! 321! boundary of domain 3 (A-domain) 322! 323 surfset(1:22)='ELECTROMAGNETICSZONE3T' 324 do i=23,81 325 surfset(i:i)=' ' 326 enddo 327 call cident81(set,surfset,nset,id) 328 nset=nset+1 329 do j=nset,id+2,-1 330 istartset(j)=istartset(j-1) 331 iendset(j)=iendset(j-1) 332 set(j)=set(j-1) 333 enddo 334 iset=id+1 335ccc to remove start 336c iset=nset 337ccc to remove end 338 istartset(iset)=nalset+1 339 do i=1,nk 340 if(ipoface(i).eq.0) cycle 341 if(inomat(i).ne.3) cycle 342 index=ipoface(i) 343 do 344 if(index.eq.0) exit 345 nalset=nalset+1 346 ialset(nalset)=10*nodface(3,index)+nodface(4,index) 347 index=nodface(5,index) 348 enddo 349 enddo 350 if(istartset(iset).gt.nalset) then 351 iactive(3)=0 352 do j=iset+1,nset 353 istartset(j-1)=istartset(j) 354 iendset(j-1)=iendset(j) 355 set(j-1)=set(j) 356 enddo 357 nset=nset-1 358 else 359 iactive(3)=iset 360ccc to remove start 361c set(iset)(1:22)='ELECTROMAGNETICSZONE3T' 362c do i=23,81 363c set(iset)(i:i)=' ' 364c enddo 365ccc to remove end 366 set(iset)=surfset 367 iendset(iset)=nalset 368 endif 369! 370! create contact ties between the domains 371! 372 do i=1,3 373 do j=1,3 374 if((i.eq.j).or.((i.eq.3).and.(j.eq.2))) cycle 375 if(iactive(i)*iactive(j).gt.0) then 376 ntie=ntie+1 377 tietol(1,ntie)=0.d0 378 tietol(2,ntie)=0.d0 379 write(tieset(1,ntie)(1:1),'(i1)')i 380 write(tieset(1,ntie)(2:2),'(i1)')j 381 do k=3,80 382 tieset(1,ntie)(k:k)=' ' 383 enddo 384 tieset(1,ntie)(81:81)='E' 385! 386! slave set 387! this set does not have to be identified as nodal 388! or facial at this stage 389! 390 tieset(2,ntie)(1:20)='ELECTROMAGNETICSZONE' 391 write(tieset(2,ntie)(21:21),'(i1)')i 392 do k=22,81 393 tieset(2,ntie)(k:k)=' ' 394 enddo 395! 396! master set 397! 398 tieset(3,ntie)(1:22)='ELECTROMAGNETICSZONE T' 399 write(tieset(3,ntie)(21:21),'(i1)')j 400 do k=23,81 401 tieset(3,ntie)(k:k)=' ' 402 enddo 403 endif 404 enddo 405 enddo 406! 407 return 408 end 409