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 genmidnodes(nktet_,ipoed,iedgmid,iexternedg, 20 & iedgext,cotet,nktet,iedg,jfix,ipoeled,ieled,kontet, 21 & iedtet,iwrite) 22! 23! generates midnodes in the middle between the neighboring vertex 24! nodes 25! 26 implicit none 27! 28 integer nktet_,ipoed(*),index,i,k,iedgmid(*),iexternedg(*), 29 & iedgext(3,*),node1,node2,iedg(3,*),nktet,jfix(*),ichange, 30 & ielem,iflag,indexe,iwrite,m,ipoeled(*),ieled(2,*), 31 & kontet(4,*),iedtet(6,*) 32! 33 real*8 cotet(3,*),xi,et,ze,shp(4,10),xsj,xl(3,10) 34! 35 include "gauss.f" 36! 37 iflag=2 38! 39! generate middle nodes 40! 41 do i=1,nktet_ 42 index=ipoed(i) 43! 44 do 45 if(index.eq.0) exit 46! 47 node1=iedg(1,index) 48 node2=iedg(2,index) 49! 50 if((jfix(node1).eq.1).and.(jfix(node2).eq.1).and. 51 & (iexternedg(index).gt.0)) then 52 iedgmid(index)=iedgext(2,iexternedg(index)) 53 else 54 nktet=nktet+1 55 iedgmid(index)=nktet 56 do k=1,3 57 cotet(k,nktet)= 58 & (cotet(k,node1)+cotet(k,node2))/2.d0 59 enddo 60 endif 61! 62 index=iedg(3,index) 63 enddo 64 enddo 65! 66! check quality of elements (may be bad due to the midnodes on 67! the fixed edges; these are not necessarily in the middle between 68! the vertex nodes) 69! 70 do 71 ichange=0 72 loop: do i=1,nktet_ 73 index=ipoed(i) 74! 75 do 76 if(index.eq.0) exit 77! 78 node1=iedg(1,index) 79 node2=iedg(2,index) 80! 81 if((jfix(node1).eq.1).and.(jfix(node2).eq.1).and. 82 & (iexternedg(index).gt.0)) then 83! 84! check whether all elements in the shell of the edge 85! have positive Jacobians 86! 87 indexe=ipoeled(index) 88 do 89 if(indexe.eq.0) exit 90 ielem=ieled(1,indexe) 91! 92! vertex nodes 93! 94 do k=1,4 95 do m=1,3 96 xl(m,k)=cotet(m,kontet(k,ielem)) 97 enddo 98 enddo 99! 100! middle nodes 101! 102 do k=1,6 103 do m=1,3 104 xl(m,k+4)=cotet(m,iedgmid(iedtet(k,ielem))) 105 enddo 106 enddo 107! 108 do k=1,4 109 xi=gauss3d5(1,k) 110 et=gauss3d5(2,k) 111 ze=gauss3d5(3,k) 112 call shape10tet(xi,et,ze,xl,xsj,shp,iflag) 113 if(xsj.le.0.d0) then 114 do m=1,3 115 cotet(m,iedgmid(index))= 116 & (cotet(m,node1)+cotet(m,node2))/2.d0 117 enddo 118 write(*,*) '*WARNING in genmidnodes: ' 119 write(*,*) ' fixed midnode ',iedgmid(index) 120 write(*,*) 121 & ' had to be moved into the middle' 122 write(*,*) 123 & ' of its neighboring vertex nodes' 124 write(*,*) ' to keep the adjacent' 125 write(*,*) ' elements regular' 126 write(*,*) 127 write(40,*) iedgmid(index) 128 iwrite=1 129 ichange=1 130 cycle loop 131 endif 132 enddo 133 indexe=ieled(2,indexe) 134 enddo 135 endif 136! 137 index=iedg(3,index) 138 enddo 139 enddo loop 140 if(ichange.eq.0) exit 141 enddo 142! 143 return 144 end 145