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 distattachline(xig,etg,pneigh,pnode,dist, 20 & nterms,xn,p) 21! 22! calculates the distance between a straight line through the node 23! with coordinates in "pnode" and direction vector "xn" and 24! the node with local coordinates xig and etg 25! in a face described by "nterms" nodes with coordinates 26! in "pneigh" 27! 28 implicit none 29! 30 integer nterms,i,j 31! 32 real*8 ratio(8),pneigh(3,*),pnode(3),dist,xi,et,xig,etg,p(3), 33 & xn(3),coeff,etm2,xim2,etm,xim,etp,xip,a2,et2,xi2,a 34! 35! 36! 37 if(nterms.eq.3) then 38 if(xig+etg.le.0.d0) then 39 ratio(2)=(xig+1.d0)/2.d0 40 ratio(3)=(etg+1.d0)/2.d0 41 else 42 ratio(2)=(1.d0-etg)/2.d0 43 ratio(3)=(1.d0-xig)/2.d0 44 endif 45 ratio(1)=1.d0-ratio(2)-ratio(3) 46 elseif(nterms.eq.4) then 47 xip=(1.d0+xig)/4.d0 48 xim=(1.d0-xig)/4.d0 49 etp=1.d0+etg 50 etm=1.d0-etg 51 ratio(1)=xim*etm 52 ratio(2)=xip*etm 53 ratio(3)=xip*etp 54 ratio(4)=xim*etp 55 elseif(nterms.eq.6) then 56 if(xig+etg.le.0.d0) then 57 xi=(xig+1.d0)/2.d0 58 et=(etg+1.d0)/2.d0 59 else 60 xi=(1.d0-etg)/2.d0 61 et=(1.d0-xig)/2.d0 62 endif 63 a=1.d0-xi-et 64 a2=2.d0*a 65 xi2=2.d0*xi 66 et2=2.d0*et 67 ratio(1)=a*(a2-1.d0) 68 ratio(2)=xi*(xi2-1.d0) 69 ratio(3)=et*(et2-1.d0) 70 ratio(4)=xi2*a2 71 ratio(5)=xi2*et2 72 ratio(6)=et2*a2 73 elseif(nterms.eq.8) then 74 xip=1.d0+xig 75 xim=1.d0-xig 76 xim2=xip*xim/2.d0 77 etp=1.d0+etg 78 etm=1.d0-etg 79 etm2=etp*etm/2.d0 80 ratio(5)=xim2*etm 81 ratio(6)=xip*etm2 82 ratio(7)=xim2*etp 83 ratio(8)=xim*etm2 84 xim=xim/4.d0 85 xip=xip/4.d0 86 ratio(1)=xim*etm*(-xig-etp) 87 ratio(2)=xip*etm*(xig-etp) 88 ratio(3)=xip*etp*(xig-etm) 89 ratio(4)=xim*etp*(-xig-etm) 90 else 91 write(*,*) '*ERROR in distattach: case with ',nterms 92 write(*,*) ' terms is not covered' 93 call exit(201) 94 endif 95! 96! calculating the position in the face 97! 98 do i=1,3 99 p(i)=0.d0 100 do j=1,nterms 101 p(i)=p(i)+ratio(j)*pneigh(i,j) 102 enddo 103 enddo 104! 105! calculating the distance 106! 107 coeff=0.d0 108 do i=1,3 109 coeff=coeff+xn(i)*(p(i)-pnode(i)) 110 enddo 111 dist=(p(1)-pnode(1)-coeff*xn(1))**2+(p(2)-pnode(2)- 112 & coeff*xn(2))**2+(p(3)-pnode(3)-coeff*xn(3))**2 113! 114 return 115 end 116 117