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