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 cattet(kontet,netet_,ifac,ne,ipkon,kon,ifatet,ifreetet, 20 & bc,itetfa,ifreefa,planfa,ipofa,cotet,cg,ipoeln,ieln,ifreeln, 21 & lakon,kontetor,iquad,istartset,iendset,ialset,set,nset,filab, 22 & jfix,iparentel,jobnamec) 23! 24! catalogueing the tetrahedral elements of the mesh 25! 26 implicit none 27! 28 character*8 lakon(*) 29 character*81 set(*),elset 30 character*87 filab(*) 31 character*132 fnrfn,jobnamec 32! 33 integer kontet(4,*),netet_,i,ifac(4,*),ne,ipkon(*),kon(*),index,j, 34 & nodes(4),ifatet(4,*),ifreetet,itetfa(2,*),ifreefa,ipofa(*), 35 & ipoeln(*),ieln(2,*),node,ifreeln,kontetor(6,*),iquad,nset, 36 & istartset(*),iendset(*),ialset(*),indexe,k,nope,jfix(*), 37 & iparentel(*),id 38! 39 real*8 bc(4,*),planfa(4,*),cotet(3,*),cg(3,*) 40! 41! open a file for nodes which are not properly projected on the 42! free surface 43! 44 open(40,file='WarnNodeNotProjected.nam',status='unknown') 45 write(40,*) '*NSET,NSET=WarnNodeNotProjected' 46! 47! read the element set name to refine, if any 48! default: all tetrahedral elements are refined 49! 50 read(filab(48)(27:87),'(a61)') elset(1:61) 51! 52 do i=62,81 53 elset(i:i)=' ' 54 enddo 55! 56! determining the number of the set 57! 58c do i=1,nset 59c if(set(i).eq.elset) exit 60c enddo 61 call cident81(set,elset,nset,id) 62 i=nset+1 63 if(id.gt.0) then 64 if(elset.eq.set(id)) then 65 i=id 66 endif 67 endif 68! 69! if element set detected: 70! 71 if(i.le.nset) then 72! 73! replace 'C' by 'A' in all tet elements 74! 75 do j=1,ne 76 if((lakon(j)(1:4).eq.'C3D4').or. 77 & (lakon(j)(1:5).eq.'C3D10')) then 78 lakon(j)(1:1)='A' 79 endif 80 enddo 81! 82! replace 'A' by 'C' in all tet elements belonging to the set 83! to refine 84! 85 do j=istartset(i),iendset(i) 86 if(ialset(j).gt.0) then 87 k=ialset(j) 88 lakon(k)(1:1)='C' 89 else 90 k=ialset(j-2) 91 do 92 k=k-ialset(j) 93 if(k.ge.ialset(j-1)) exit 94 lakon(k)(1:1)='C' 95 enddo 96 endif 97 enddo 98! 99! setting field jfix in all VERTEX nodes belonging to 100! elements not to be refined to 1 (including the common 101! boundaries between the set to be refined and its complement) 102! 103 do i=1,ne 104 if(lakon(i)(1:5).eq.'C3D8I') then 105 nope=8 106 elseif(lakon(i)(4:5).eq.'20') then 107 nope=8 108 elseif(lakon(i)(4:4).eq.'8') then 109 nope=8 110 elseif(lakon(i)(1:5).eq.'C3D10') then 111 cycle 112 elseif(lakon(i)(1:4).eq.'C3D4') then 113 cycle 114 elseif(lakon(i)(4:5).eq.'15') then 115 nope=6 116 elseif(lakon(i)(4:4).eq.'6') then 117 nope=6 118 elseif(lakon(i)(1:2).eq.'ES') then 119 if(lakon(i)(7:7).eq.'C') then 120 cycle 121 else 122 nope=ichar(lakon(i)(8:8))-47 123 endif 124 elseif(lakon(i)(1:4).eq.'MASS') then 125 nope=1 126 elseif(lakon(i)(1:1).eq.'A') then 127 nope=4 128 else 129 cycle 130 endif 131! 132 indexe=ipkon(i) 133 do j=1,nope 134 jfix(kon(indexe+j))=1 135 enddo 136 enddo 137 endif 138! 139! change on 20200327: original element number should not be 140! overwritten: needed at reading time for ielmat, ielorien, 141! ielbody ... 142! 143c! 144c! determine the first tetrahedral element or first 145c! unused element 146c! 147c do i=1,ne 148c if((lakon(i)(1:4).eq.'C3D4').or. 149c & (lakon(i)(1:5).eq.'C3D10').or. 150c & (lakon(i)(1:1).eq.char(0))) exit 151c enddo 152! 153! determine the first unused element 154! 155 do i=1,ne 156 if(lakon(i)(1:1).eq.char(0)) exit 157 enddo 158c 159 ifreetet=i 160! 161 do 162 do j=i+1,netet_ 163! 164! element number supersedes largest one 165! 166 if(j.gt.ne) then 167 kontet(4,i)=j 168 i=j 169 exit 170 endif 171c! 172c! tetrahedral element 173c! 174c if((lakon(j)(1:4).eq.'C3D4').or. 175c & (lakon(j)(1:5).eq.'C3D10')) then 176c kontet(4,i)=j 177c i=j 178c exit 179c endif 180! 181! unused element 182! 183 if(lakon(j)(1:1).eq.char(0)) then 184 kontet(4,i)=j 185 i=j 186 exit 187 endif 188 enddo 189 if(j.eq.netet_) exit 190 enddo 191 kontet(4,netet_)=0 192! 193! initialization of ipofa and ifac 194! 195 do i=1,4*netet_ 196 ifac(4,i)=i+1 197 enddo 198 ifac(4,4*netet_)=0 199! 200! initialization of ieln 201! 202 do i=1,4*netet_ 203 ieln(2,i)=i+1 204 enddo 205 ieln(2,4*netet_)=0 206! 207! adding the tetrahedral elements one by one in kontet 208! tagging the original elements in kon to be removed 209! by using a model change in input deck jobname.rfn.inp 210! 211! creating the file name 212! 213 do i=1,132 214 if(ichar(jobnamec(i:i)).eq.0) exit 215 enddo 216 if(i.gt.125) then 217 write(*,*) '*ERROR in cattet' 218 write(*,*) ' jobname has more than 124 characters:' 219 write(*,*) jobnamec(1:132) 220 call exit(201) 221 endif 222 fnrfn(1:i+7)=jobnamec(1:i-1)//'.rfn.inp' 223! 224! opening the file 225! 226 open(2,file=fnrfn(1:i+7),status='unknown') 227! 228 do i=1,ne 229 if(ipkon(i).lt.0) cycle 230 if((lakon(i)(1:4).ne.'C3D4').and. 231 & (lakon(i)(1:5).ne.'C3D10')) cycle 232 index=ipkon(i) 233 do j=1,4 234 nodes(j)=kon(index+j) 235 enddo 236 write(2,101) 237 write(2,102) i 238! 239! if C3D10: store the middle nodes 240! If there is at least one C3D10 element iquad is set to 1 241! which means that the refined mesh will be fully quadratic 242! 243 if(lakon(i)(4:4).eq.'1') then 244 iquad=1 245 do j=1,6 246 kontetor(j,ifreetet)=kon(index+4+j) 247 enddo 248 endif 249! 250! defining the initial parent element (= element itself) 251! 252 iparentel(ifreetet)=i 253 call generatetet_refine(kontet,ifatet,ifreetet,bc,ifac,itetfa, 254 & ifreefa,planfa,ipofa,nodes,cotet,cg) 255 enddo 256! 257 101 format('*MODEL CHANGE,TYPE=ELEMENT,REMOVE') 258 102 format(i10) 259 close(2) 260! 261! generating the element per node relationship 262! 263 do j=1,netet_ 264 if(kontet(1,j).eq.0) cycle 265 do i=1,4 266 node=kontet(i,j) 267 index=ifreeln 268 ieln(1,index)=j 269 ifreeln=ieln(2,index) 270 if(ifreeln.eq.0) then 271 write(*,*) '*ERROR in cattet: increase the' 272 write(*,*) ' dimension of ieln' 273 call exit(201) 274 endif 275 ieln(2,index)=ipoeln(node) 276 ipoeln(node)=index 277 enddo 278 enddo 279! 280 return 281 end 282 283