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 adjacentbounodes(ifront,ifrontrel,nfront,iedno,iedg, 20 & nnfront,ifront2,ifrontrel2,ibounnod,nbounnod,istartfront, 21 & iendfront,ibounedg,istartcrackfro,iendcrackfro,ncrack, 22 & ibounnod2,istartcrackbou,iendcrackbou,isubsurffront,iedno2, 23 & stress,stress2,iresort,ieled,kontri,costruc,costruc2,temp, 24 & temp2,nstep,ier) 25! 26! sorting the crack boundary nodes according to adjacency 27! 28 implicit none 29! 30 integer ifront(*),ifrontrel(*),nfront,iedno(2,*),iedg(3,*), 31 & ichange,nnfront,nfront2,node,noderel,iedge,ibounnod(*), 32 & nbounnod,i,j,k,istartfront(*),iendfront(*),ifront2(*), 33 & ifrontrel2(*),id,ibounedg(*),nodestart,istartcrackfro(*), 34 & iendcrackfro(*),ncrack,iact,nodeprev,noderelprev,kflag, 35 & ibounnod2(*),istartcrackbou(*),iendcrackbou(*),nbounnod2, 36 & isubsurffront(*),isubsurf,iedno2(2,*),iresort(*),itri, 37 & nodenext,ieled(2,*),kontri(3,*),nstep,ier 38! 39 real*8 stress(6,nstep,*),stress2(6,nstep,*),costruc(3,*), 40 & costruc2(3,*),temp(nstep,*),temp2(nstep,*) 41! 42! nnfront is the number of distinct fronts 43! 44 nnfront=0 45 nfront2=0 46 nbounnod2=0 47 ncrack=0 48 kflag=2 49c do i=1,nfront 50c write(*,*) 'adjacentbounodes ',i,ifront(i) 51c enddo 52! 53 do 54! 55! for all recorded front nodes i the entry in ifrontrel(i) is 56! set to 0. If no zero is left, the job is done. 57! 58 ichange=0 59! 60 loop: do i=1,nfront 61 if(ifrontrel(i).gt.0) then 62 ichange=1 63! 64! find the end of a front: start with an arbitrary node 65! 66c write(*,*) 'look for new front ' 67 node=ifront(i) 68c write(*,*) node 69 nodestart=node 70 noderel=ifrontrel(i) 71! 72! take one of the edges to which the node belongs 73! check whether the nodes belonging to the edge are in the 74! same order as in the element topology of the crack 75! 76 iedge=ibounedg(iedno(1,noderel)) 77! 78! crack triangle belonging to the edge 79! 80 itri=ieled(1,iedge) 81! 82! finding the other node of the edge 83! 84 if(iedg(1,iedge).ne.node) then 85 nodenext=iedg(1,iedge) 86 else 87 nodenext=iedg(2,iedge) 88 endif 89! 90! node and nodenext also belong to itri; check whether the 91! order in the element topology is nodenext -> node 92! (because the order is reverted as soon as an external node 93! is found (subsurface crack) or the start node is retreated 94! (surface crack) 95! 96 do j=1,3 97 if(kontri(j,itri).eq.node) exit 98 enddo 99 k=j+1 100 if(k.gt.3) k=1 101c k=mod(j+1,3) 102! 103! take the other edge is the order is not the same 104! 105 if(kontri(k,itri).eq.nodenext) then 106 iedge=ibounedg(iedno(2,noderel)) 107 endif 108! 109 do 110! 111! find the other node on the edge 112! 113 if(iedg(1,iedge).ne.node) then 114 node=iedg(1,iedge) 115 else 116 node=iedg(2,iedge) 117 endif 118c write(*,*) node 119 if(node.eq.nodestart) then 120! 121! subsurface crack 122! 123 call nident(ifront,node,nfront,id) 124 ifrontrel(id)=0 125 isubsurf=1 126 exit 127 endif 128! 129! check whether this is a front node 130! 131 call nident(ifront,node,nfront,id) 132 if(id.gt.0) then 133 if(ifront(id).eq.node) then 134! 135! front node! search for other edge adjacent to the node 136! 137 noderel=ifrontrel(id) 138 if(ibounedg(iedno(1,noderel)).eq.iedge) then 139 iedge=ibounedg(iedno(2,noderel)) 140 else 141 iedge=ibounedg(iedno(1,noderel)) 142 endif 143 cycle 144 endif 145 endif 146! 147 isubsurf=0 148 exit 149 enddo 150! 151! if surface crack: 152! node is a boundary node which is external to the structure 153! but adjacent to a front node (i.e. a boundary node internal 154! to the structure) 155! 156! if subsurface crack: 157! node is an arbitary node belonging to the crack front 158! 159! finding the relative node number 160! 161 call nident(ibounnod,node,nbounnod,noderel) 162c write(*,*) 'adjacentbounodes start of new front' 163! 164! defining a new crack 165! 166 ncrack=ncrack+1 167 nodestart=node 168! 169! defining a new front segment 170! 171 nnfront=nnfront+1 172 iact=1 173! 174! new node belonging to a segment (a segment contains the internal 175! nodes plus the immediately adjacent external nodes 176! 177 nfront2=nfront2+1 178 ifront2(nfront2)=node 179! 180 nbounnod2=nbounnod2+1 181 ibounnod2(nbounnod2)=node 182! 183c write(*,*) node 184 ifrontrel2(nfront2)=noderel 185 istartfront(nnfront)=nfront2 186 isubsurffront(nnfront)=isubsurf 187 istartcrackfro(ncrack)=nfront2 188 istartcrackbou(ncrack)=nbounnod2 189! 190! looking for other nodes belonging to the segment 191! 192 do 193! 194! find the other node on the edge 195! 196 if(iedg(1,iedge).ne.node) then 197 nodeprev=node 198 noderelprev=noderel 199 node=iedg(1,iedge) 200 else 201 nodeprev=node 202 noderelprev=noderel 203 node=iedg(2,iedge) 204 endif 205! 206! reached starting node 207! 208 if(node.eq.nodestart) then 209 iendcrackbou(ncrack)=nbounnod2 210 iendfront(nnfront)=nfront2 211 exit 212 endif 213! 214! check whether this is a front node 215! 216 call nident(ifront,node,nfront,id) 217 if(id.gt.0) then 218 if(ifront(id).eq.node) then 219! 220! front node ! 221! 222! check whether new front is to be activated 223! 224 if(iact.eq.0) then 225! 226! new front 227! 228 iact=1 229 nnfront=nnfront+1 230 isubsurffront(nnfront)=isubsurf 231! 232! first node of new front 233! 234 nfront2=nfront2+1 235 ifront2(nfront2)=nodeprev 236 ifrontrel2(nfront2)=noderelprev 237 istartfront(nnfront)=nfront2 238 endif 239! 240 noderel=ifrontrel(id) 241 if(ibounedg(iedno(1,noderel)).eq.iedge) then 242 iedge=ibounedg(iedno(2,noderel)) 243 else 244 iedge=ibounedg(iedno(1,noderel)) 245 endif 246! 247 nfront2=nfront2+1 248 ifront2(nfront2)=node 249! 250 nbounnod2=nbounnod2+1 251 ibounnod2(nbounnod2)=node 252! 253c write(*,*) node 254 ifrontrel2(nfront2)=noderel 255! 256! setting the entry in frontrel for a node which was treated 257! to zero 258! 259 ifrontrel(id)=0 260 cycle 261 endif 262 endif 263! 264! no front node! 265! 266 call nident(ibounnod,node,nbounnod,noderel) 267 if(ibounedg(iedno(1,noderel)).eq.iedge) then 268 iedge=ibounedg(iedno(2,noderel)) 269 else 270 iedge=ibounedg(iedno(1,noderel)) 271 endif 272! 273 nbounnod2=nbounnod2+1 274 ibounnod2(nbounnod2)=node 275! 276 if(iact.eq.1) then 277! 278! last node of the segment (external boundary node immediately 279! adjacent to an internal node) 280! 281 nfront2=nfront2+1 282 ifront2(nfront2)=node 283c write(*,*) node 284 ifrontrel2(nfront2)=noderel 285! 286! finishing the actual segment 287! 288 iendfront(nnfront)=nfront2 289 iact=0 290 endif 291 enddo 292! 293! end of crack 294! 295 iendcrackfro(ncrack)=iendfront(nnfront) 296 endif 297 enddo loop ! end loop over the front nodes 298! 299! if no front node is left untreated, leave 300! 301 if(ichange.eq.0) exit 302 enddo 303! 304! adjust the dependent fields of ibounnod 305! 306 do i=1,nbounnod2 307 node=ibounnod2(i) 308 call nident(ibounnod,node,nbounnod,id) 309 do j=1,2 310 iedno2(j,i)=iedno(j,id) 311 enddo 312 do k=1,nstep 313 temp2(k,i)=temp(k,id) 314 enddo 315 do k=1,nstep 316 do j=1,6 317 stress2(j,k,i)=stress(j,k,id) 318 enddo 319 enddo 320 do j=1,3 321 costruc2(j,i)=costruc(j,id) 322 enddo 323 enddo 324! 325! copy the ibounnod2 field and its dependent fields 326! 327 do i=1,nbounnod2 328 ibounnod(i)=ibounnod2(i) 329 do j=1,2 330 iedno(j,i)=iedno2(j,i) 331 enddo 332 do k=1,nstep 333 temp(k,i)=temp2(k,i) 334 enddo 335 do k=1,nstep 336 do j=1,6 337 stress(j,k,i)=stress2(j,k,i) 338 enddo 339 enddo 340 do j=1,3 341 costruc(j,i)=costruc2(j,i) 342 enddo 343 enddo 344! 345! calculating in which order ibounnod was resorted 346! 347 do i=1,nbounnod2 348 iresort(i)=i 349 enddo 350! 351 nfront=nfront 352 call isortii(ibounnod2,iresort,nbounnod2,kflag) 353! 354! copy the temporary fields (ifront2 and ifrontrel2) 355! 356 do i=1,nfront2 357 ifront(i)=ifront2(i) 358 ifrontrel(i)=iresort(ifrontrel2(i)) 359 enddo 360 nfront=nfront2 361! 362c write(*,*) 363c write(*,*) 'adjacentbounodes: front nodes' 364c write(*,*) 365c do i=1,nfront 366c write(*,*) i,ifront(i) 367c enddo 368! 369 if(nfront.eq.0) then 370 write(*,*) '*ERROR in adjacentbounodes: no front node found' 371 ier=1 372c call exit(201) 373 endif 374! 375 return 376 end 377 378