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 beamextscheme(yil,ndim,nfield,lakonl,npropstart,prop, 20 & field,mi) 21! 22! provides the extrapolation scheme for beams with a cross section 23! which is not rectangular nor elliptical 24! 25 implicit none 26! 27 character*8 lakonl 28! 29 integer npropstart,mi(*),ndim,nfield,j,k,l 30! 31 real*8 prop(*),ratio,ratio2,yil(ndim,mi(1)),yig(nfield,mi(1)), 32 & field(999,20*mi(3)),r,scal,a8(8,8),pmean1(nfield), 33 & pmean2(nfield),t1,t2,t3,t4,a,b,dummy,shp(4,20),xis(8,3) 34! 35! extrapolation from a 2x2x2=8 integration point scheme in a hex to 36! the vertex nodes 37! 38 data a8 /2.549d0,-.683d0,.183d0,-.683d0,-.683d0,.183d0, 39 & -.04904d0,.183d0,-.683d0,2.549d0,-.683d0,.183d0, 40 & .183d0,-.683d0,.183d0,-.04904d0,-.683d0,.183d0, 41 & -.683d0,2.549d0,.183d0,-.04904d0,.183d0,-.683d0, 42 & .183d0,-.683d0,2.549d0,-.683d0,-.04904d0,.183d0, 43 & -.683d0,.183d0,-.683d0,.183d0,-.04904d0,.183d0, 44 & 2.549d0,-.683d0,.183d0,-.683d0,.183d0,-.683d0, 45 & .183d0,-.04904d0,-.683d0,2.549d0,-.683d0,.183d0, 46 & .183d0,-.04904d0,.183d0,-.683d0,-.683d0,.183d0, 47 & -.683d0,2.549d0,-.04904d0,.183d0,-.683d0,.183d0, 48 & .183d0,-.683d0,2.549d0,-.683d0/ 49! 50 if(lakonl(8:8).eq.'P') then 51! 52! pipe cross section 53! 54! the axis of the pipe is along the local xi-direction 55! the integration points are at positions +-0.57 along 56! the xi-axis. At each of these positions there are 8 57! integration points in the eta-zeta plane at one radial 58! position and 59! equally spaced along the circumferential direction 60! 61! ratio of inner radius to outer radius 62! 63 ratio=(prop(npropstart+1)-prop(npropstart+2))/ 64 & prop(npropstart+1) 65 ratio2=ratio*ratio 66! 67! radial location of integration points 68! 69 r=dsqrt((ratio2+1.d0)/2.d0) 70! 71! scaling factor between the radial location of the integration 72! points and the regular location of C3D20R integration points 73! 74 scal=dsqrt(2.d0/3.d0)/r 75! 76! calculating the mean values at each of the two xi-positions 77! 78 do k=1,nfield 79 pmean1(k)=(yil(k,2)+yil(k,4)+yil(k,6)+yil(k,8))/4.d0 80 pmean2(k)=(yil(k,10)+yil(k,12)+yil(k,14)+yil(k,16))/4.d0 81 enddo 82! 83! translating the results from the integration points of the 84! pipe section to the integration points of the C3D20R element 85! 86 do k=1,nfield 87 yig(k,1)=pmean1(k)+(yil(k,6)-pmean1(k))*scal 88 yig(k,2)=pmean2(k)+(yil(k,14)-pmean2(k))*scal 89 yig(k,3)=pmean1(k)+(yil(k,8)-pmean1(k))*scal 90 yig(k,4)=pmean2(k)+(yil(k,16)-pmean2(k))*scal 91 yig(k,5)=pmean1(k)+(yil(k,4)-pmean1(k))*scal 92 yig(k,6)=pmean2(k)+(yil(k,12)-pmean2(k))*scal 93 yig(k,7)=pmean1(k)+(yil(k,2)-pmean1(k))*scal 94 yig(k,8)=pmean2(k)+(yil(k,10)-pmean2(k))*scal 95 enddo 96! 97! standard extrapolation for the C3D20R element 98! 99 do j=1,8 100 do k=1,nfield 101 field(k,j)=0.d0 102 do l=1,8 103 field(k,j)=field(k,j)+a8(j,l)*yig(k,l) 104 enddo 105 enddo 106 enddo 107! 108 elseif(lakonl(8:8).eq.'B') then 109! 110! BOX cross section 111 a=prop(npropstart+1) 112 b=prop(npropstart+2) 113 t1=prop(npropstart+3) 114 t2=prop(npropstart+4) 115 t3=prop(npropstart+5) 116 t4=prop(npropstart+6) 117c 118c new local coordinates for node points of element 119c 120 xis(1,1) = -dsqrt(3.0d0) 121 xis(1,2) = -(t4-t2-2*b)/((-2*b)+t2+t4) 122 xis(1,3) = (t3-t1+2*a)/((-2*a)+t1+t3) 123 xis(2,1) = dsqrt(3.0d0) 124 xis(2,2) = -(t4-t2-2*b)/((-2*b)+t2+t4) 125 xis(2,3) = (t3-t1+2*a)/((-2*a)+t1+t3) 126 xis(3,1) = dsqrt(3.0d0) 127 xis(3,2) = -(t4-t2+2*b)/((-2*b)+t2+t4) 128 xis(3,3) = (t3-t1+2*a)/((-2*a)+t1+t3) 129 xis(4,1) = -dsqrt(3.0d0) 130 xis(4,2) = -(t4-t2+2*b)/((-2*b)+t2+t4) 131 xis(4,3) = (t3-t1+2*a)/((-2*a)+t1+t3) 132 xis(5,1) = -dsqrt(3.0d0) 133 xis(5,2) = -(t4-t2-2*b)/((-2*b)+t2+t4) 134 xis(5,3) = (t3-t1-2*a)/((-2*a)+t1+t3) 135 xis(6,1) = dsqrt(3.0d0) 136 xis(6,2) = -(t4-t2-2*b)/((-2*b)+t2+t4) 137 xis(6,3) = (t3-t1-2*a)/((-2*a)+t1+t3) 138 xis(7,1) = dsqrt(3.0d0) 139 xis(7,2) = -(t4-t2+2*b)/((-2*b)+t2+t4) 140 xis(7,3) = (t3-t1-2*a)/((-2*a)+t1+t3) 141 xis(8,1) = -dsqrt(3.0d0) 142 xis(8,2) = -(t4-t2+2*b)/((-2*b)+t2+t4) 143 xis(8,3) = (t3-t1-2*a)/((-2*a)+t1+t3) 144 xis(8,1) = -dsqrt(3.0d0) 145 xis(8,2) = -(t4-t2+2*b)/((-2*b)+t2+t4) 146 xis(8,3) = (t3-t1-2*a)/((-2*a)+t1+t3) 147! 148! extrapolate from int points to node points 149! 150 do k=1,nfield 151 yig(k,1)=yil(k,9) 152 yig(k,2)=yil(k,25) 153 yig(k,3)=yil(k,29) 154 yig(k,4)=yil(k,13) 155 yig(k,5)=yil(k,5) 156 yig(k,6)=yil(k,21) 157 yig(k,7)=yil(k,17) 158 yig(k,8)=yil(k,1) 159 enddo 160! 161 do j=1,8 162 call shape8h(xis(j,1),xis(j,2),xis(j,3),dummy,dummy,shp,1) 163 do k=1,nfield 164 field(k,j)=0.0d0 165 do l=1,8 166 field(k,j)=field(k,j)+shp(4,l)*yig(k,l) 167 enddo 168 enddo 169 enddo 170! 171 endif 172! 173 return 174 end 175