1! 2! CalculiX - A 3-dimensional finite element program 3! Copyright (C) 1998 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 checkjac(cotet,node,pnew,kontet,c1,jflag, 20 & iedtet,iedgmid,ipoeled,ieled,iedge,ibadnodes,nbadnodes, 21 & iwrite) 22! 23! check the Jacobian in all integration points of all elements 24! belonging to midnode "node" if the coordinates of "node" are 25! changed from cotet(1..3,node) to pnew(1..3); if the Jacobian 26! is strictly positive, the coordinates are changed, if not, they 27! are left unchanged 28! "node" belongs to edge "iedge" 29! 30 implicit none 31! 32 integer node,kontet(4,*),i,ielem,indexe,iflag,jflag,iwrite, 33 & j,n,iedtet(6,*),iedgmid(*),ipoeled(*),ieled(2,*), 34 & iedge,k,m,ibadnodes(*),nbadnodes,id,listed,neighedge 35! 36 real*8 cotet(3,*),pnew(3),pold(3),xsj,c1,c2,xi,et,ze,xl(3,10), 37 & shp(4,10) 38! 39 include "gauss.f" 40! 41 iflag=2 42! 43! backup of the old coordinates 44! 45 do i=1,3 46 pold(i)=cotet(i,node) 47 enddo 48! 49! storing the new coordinates in cotet 50! 51 c2=c1 52 do i=1,3 53 cotet(i,node)=c2*pnew(i)+(1.d0-c2)*pold(i) 54 enddo 55! 56! trial loops; if bad elements the projection is reduced and the 57! elements are rechecked 58! 59 n=3 60 loop1: do j=1,n 61 indexe=ipoeled(iedge) 62! 63! loop over all elements belonging to edge "iedge" 64! 65 loop2: do 66 if(indexe.eq.0) exit loop1 67 ielem=ieled(1,indexe) 68! 69! vertex nodes 70! 71 do k=1,4 72 do m=1,3 73 xl(m,k)=cotet(m,kontet(k,ielem)) 74 enddo 75 enddo 76! 77! middle nodes 78! 79 do k=1,6 80 do m=1,3 81 xl(m,k+4)=cotet(m,iedgmid(iedtet(k,ielem))) 82 enddo 83 enddo 84! 85 do k=1,4 86 xi=gauss3d5(1,k) 87 et=gauss3d5(2,k) 88 ze=gauss3d5(3,k) 89 call shape10tet(xi,et,ze,xl,xsj,shp,iflag) 90 if(xsj.le.0.d0) then 91c write(*,*) 'checkjac bad element ',ielem,' node ',node 92 exit loop2 93 endif 94 enddo 95 indexe=ieled(2,indexe) 96 enddo loop2 97! 98! last loop: revert to non-projected coordinates 99! 100 if(j.eq.n) then 101 do i=1,3 102 cotet(i,node)=pold(i) 103 enddo 104 exit 105 endif 106! 107! not last loop: reduce the projection and retry 108! 109 c2=c2/2.d0 110 do i=1,3 111 cotet(i,node)=c2*pnew(i)+(1.d0-c2)*pold(i) 112 enddo 113 enddo loop1 114! 115! if node was not completely projected: list all edges of 116! neighboring elements (corresponds 1-to-1 to midnodes) 117! -> application of fminsirefine for all neighboring edges/midnodes 118! which are not on the free surface 119! 120 if(j.ne.1) then 121 indexe=ipoeled(iedge) 122 do 123 if(indexe.eq.0) exit 124! 125! adjacent element 126! 127 ielem=ieled(1,indexe) 128 do k=1,6 129 neighedge=iedtet(k,ielem) 130! 131! edge of adjacent element 132! 133 listed=0 134 call nident(ibadnodes,neighedge,nbadnodes,id) 135 if(id.gt.0) then 136 if(ibadnodes(id).eq.neighedge) then 137 listed=1 138 endif 139 endif 140 if(listed.eq.0) then 141 nbadnodes=nbadnodes+1 142 do m=nbadnodes,id+2,-1 143 ibadnodes(m)=ibadnodes(m-1) 144 enddo 145 ibadnodes(id+1)=neighedge 146 endif 147 enddo 148 indexe=ieled(2,indexe) 149 enddo 150 endif 151! 152! if jflag=1 and j>1: last chance for a complete projection 153! did not work 154! 155 if((jflag.eq.1).and.(j.ne.1)) then 156 write(*,*) '*WARNING in checkjac: projection of midnode ',node 157 write(*,*) ' had to be reduced to keep the adjacent' 158 write(*,*) ' elements regular' 159 write(*,*) 160 write(40,*) node 161 iwrite=1 162 endif 163! 164 return 165 end 166