CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
extrapolate.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine extrapolate (yi, yn, ipkon, inum, kon, lakon, nfield, nk, ne, mi, ndim, orab, ielorien, co, iorienloc, cflag, vold, force, ielmat, thicke, ielprop, prop)
 

Function/Subroutine Documentation

◆ extrapolate()

subroutine extrapolate ( real*8, dimension(ndim,mi(1),*)  yi,
real*8, dimension(nfield,*)  yn,
integer, dimension(*)  ipkon,
integer, dimension(*)  inum,
integer, dimension(*)  kon,
character*8, dimension(*)  lakon,
integer  nfield,
integer  nk,
integer  ne,
integer, dimension(*)  mi,
integer  ndim,
real*8, dimension(7,*)  orab,
integer, dimension(mi(3),*)  ielorien,
real*8, dimension(3,*)  co,
integer  iorienloc,
character*1  cflag,
real*8, dimension(0:mi(2),*)  vold,
logical  force,
integer, dimension(mi(3),*)  ielmat,
real*8, dimension(mi(3),*)  thicke,
integer, dimension(*)  ielprop,
real*8, dimension(*)  prop 
)
22 !
23 ! extrapolates field values at the integration points to the
24 ! nodes
25 !
26 ! the C3D20RB element has 50 integration points, however, the
27 ! first 8 integration points coincide with those of a C3D20R
28 ! element. In this routine the C3D20RBR and C3D20RBC
29 ! elements are treated as an ordinary C3D20R element
30 !
31 ! the number of internal state variables is limited to 999
32 ! (cfr. array field)
33 !
34  implicit none
35 !
36  logical force
37 !
38  character*1 cflag
39  character*8 lakon(*),lakonl
40 !
41  integer ipkon(*),inum(*),kon(*),mi(*),ne,indexe,nope,
42  & nonei20(3,12),nfield,nonei10(3,6),nk,i,j,k,l,ndim,
43  & nonei15(3,9),iorienloc,iorien,ielorien(mi(3),*),konl,
44  & mint3d,m,iflag,jj,ll,ielmat(mi(3),*),ielprop(*),
45  & nlayer,nopeexp,ilayer,kk,mint2d,nopes,kl,ki,null,
46  & itet(4),iwedge(2,9)
47 !
48  real*8 yi(ndim,mi(1),*),yn(nfield,*),field(999,20*mi(3)),a8(8,8),
49  & a4(4,4),a27(20,27),a9(6,9),a2(6,2),orab(7,*),co(3,*),prop(*),
50  & coords(3,mi(1)),xi,et,ze,xl(3,20),xsj,shp(4,20),weight,
51  & yiloc(6,mi(1)),a(3,3),b(3,3),c(3,3),vold(0:mi(2),*),tlayer(4),
52  & dlayer(4),xlayer(mi(3),4),thickness,xs2(3,7),xl2(3,8),
53  & xsj2(3),shp2(7,8),thicke(mi(3),*),coloc(3,8),
54  & xwedge(2,2,9),a14(8,14),a6(6,6)
55 !
56  include "gauss.f"
57 !
58  data coloc /-1.d0,-1.d0,-1.d0,
59  & 1.d0,-1.d0,-1.d0,
60  & -1.d0,1.d0,-1.d0,
61  & 1.d0,1.d0,-1.d0,
62  & -1.d0,-1.d0,1.d0,
63  & 1.d0,-1.d0,1.d0,
64  & -1.d0,1.d0,1.d0,
65  & 1.d0,1.d0,1.d0/
66  data null /0/
67 !
68 ! a 10-node tet is remeshed into 10 4-node tets at contact
69 ! interfaces; itet contains the linear tet elements to which
70 ! the integration points of the parent 10-node tet belong
71 !
72  data itet /1,2,3,10/
73 !
74  data iwedge /1,0,2,0,3,0,1,5,2,6,3,7,5,0,6,0,7,0/
75 !
76  data xwedge /0.975615382715435242d0,0.243846172845647580d-1,
77  & 0.d0,0.d0,
78  & 0.975615382715435242d0,0.243846172845647580d-1,
79  & 0.d0,0.d0,
80  & 0.975615382715435242d0,0.243846172845647580d-1,
81  & 0.d0,0.d0,
82  & 0.d0,.5d0,.5d0,0.d0,
83  & 0.d0,.5d0,.5d0,0.d0,
84  & 0.d0,.5d0,.5d0,0.d0,
85  & 0.243846172845647580d-1,0.975615382715435242d0,
86  & 0.d0,0.d0,
87  & 0.243846172845647580d-1,0.975615382715435242d0,
88  & 0.d0,0.d0,
89  & 0.243846172845647580d-1,0.975615382715435242d0,
90  & 0.d0,0.d0/
91 !
92  data nonei10 /5,1,2,6,2,3,7,3,1,8,1,4,9,2,4,10,3,4/
93 !
94  data nonei15 /7,1,2,8,2,3,9,3,1,10,4,5,11,5,6,12,6,4,
95  & 13,1,4,14,2,5,15,3,6/
96 !
97  data nonei20 /9,1,2,10,2,3,11,3,4,12,4,1,
98  & 13,5,6,14,6,7,15,7,8,16,8,5,
99  & 17,1,5,18,2,6,19,3,7,20,4,8/
100 !
101  data a2 / 1.1455,-0.1455,1.1455,-0.1455,1.1455,-0.1455,
102  & -0.1455,1.1455,-0.1455,1.1455,-0.1455,1.1455/
103  data a4 / 1.92705, -0.30902, -0.30902, -0.30902,
104  & -0.30902, 1.92705, -0.30902, -0.30902,
105  & -0.30902, -0.30902, 1.92705, -0.30902,
106  & -0.30902, -0.30902, -0.30902, 1.92705/
107 !
108 ! extrapolation from a 6 integration point scheme in a wedge to
109 ! the vertex nodes
110 !
111  data a6 / 2.04904, 0.00000, 0.00000,-0.54904, 0.00000, 0.00000,
112  & -0.34151, 1.70753,-0.34151, 0.09151,-0.45753, 0.09151,
113  & -0.34151,-0.34151, 1.70753, 0.09151, 0.09151,-0.45753,
114  & -0.54904, 0.00000, 0.00000, 2.04904, 0.00000, 0.00000,
115  & 0.09151,-0.45753, 0.09151,-0.34151, 1.70753,-0.34151,
116  & 0.09151, 0.09151,-0.45753,-0.34151,-0.34151, 1.70753/
117 !
118  data a9 / 1.63138,-0.32628,-0.32628,-0.52027, 0.10405, 0.10405,
119  & -0.32628, 1.63138,-0.32628, 0.10405,-0.52027, 0.10405,
120  & -0.32628,-0.32628, 1.63138, 0.10405, 0.10405,-0.52027,
121  & 0.55556,-0.11111,-0.11111, 0.55556,-0.11111,-0.11111,
122  & -0.11111, 0.55556,-0.11111,-0.11111, 0.55556,-0.11111,
123  & -0.11111,-0.11111, 0.55556,-0.11111,-0.11111, 0.55556,
124  & -0.52027, 0.10405, 0.10405, 1.63138,-0.32628,-0.32628,
125  & 0.10405,-0.52027, 0.10405,-0.32628, 1.63138,-0.32628,
126  & 0.10405, 0.10405,-0.52027,-0.32628,-0.32628, 1.63138/
127 !
128 ! extrapolation from a 2x2x2=8 integration point scheme in a hex to
129 ! the vertex nodes
130 !
131  data a8 /2.549,-.683,.183,-.683,-.683,.183,
132  & -.04904,.183,-.683,2.549,-.683,.183,
133  & .183,-.683,.183,-.04904,-.683,.183,
134  & -.683,2.549,.183,-.04904,.183,-.683,
135  & .183,-.683,2.549,-.683,-.04904,.183,
136  & -.683,.183,-.683,.183,-.04904,.183,
137  & 2.549,-.683,.183,-.683,.183,-.683,
138  & .183,-.04904,-.683,2.549,-.683,.183,
139  & .183,-.04904,.183,-.683,-.683,.183,
140  & -.683,2.549,-.04904,.183,-.683,.183,
141  & .183,-.683,2.549,-.683/
142 !
143 ! extrapolation from a 14 integration point scheme in a hex to
144 ! the vertex nodes
145 !
146  data a14 /
147  & 0.1396e+01,-0.3026e+00,0.1124e-01,-0.3026e+00,
148  & -0.3026e+00,0.1124e-01,0.4901e-01,0.1124e-01,
149  & -0.3026e+00,0.1396e+01,-0.3026e+00,0.1124e-01,
150  & 0.1124e-01,-0.3026e+00,0.1124e-01,0.4901e-01,
151  & 0.1124e-01,-0.3026e+00,0.1396e+01,-0.3026e+00,
152  & 0.4901e-01,0.1124e-01,-0.3026e+00,0.1124e-01,
153  & -0.3026e+00,0.1124e-01,-0.3026e+00,0.1396e+01,
154  & 0.1124e-01,0.4901e-01,0.1124e-01,-0.3026e+00,
155  & -0.3026e+00,0.1124e-01,0.4901e-01,0.1124e-01,
156  & 0.1396e+01,-0.3026e+00,0.1124e-01,-0.3026e+00,
157  & 0.1124e-01,-0.3026e+00,0.1124e-01,0.4901e-01,
158  & -0.3026e+00,0.1396e+01,-0.3026e+00,0.1124e-01,
159  & 0.4901e-01,0.1124e-01,-0.3026e+00,0.1124e-01,
160  & 0.1124e-01,-0.3026e+00,0.1396e+01,-0.3026e+00,
161  & 0.1124e-01,0.4901e-01,0.1124e-01,-0.3026e+00,
162  & -0.3026e+00,0.1124e-01,-0.3026e+00,0.1396e+01,
163  & 0.2069e+00,0.2069e+00,-0.6408e-01,-0.6408e-01,
164  & 0.2069e+00,0.2069e+00,-0.6408e-01,-0.6408e-01,
165  & -0.6408e-01,0.2069e+00,0.2069e+00,-0.6408e-01,
166  & -0.6408e-01,0.2069e+00,0.2069e+00,-0.6408e-01,
167  & -0.6408e-01,-0.6408e-01,0.2069e+00,0.2069e+00,
168  & -0.6408e-01,-0.6408e-01,0.2069e+00,0.2069e+00,
169  & 0.2069e+00,-0.6408e-01,-0.6408e-01,0.2069e+00,
170  & 0.2069e+00,-0.6408e-01,-0.6408e-01,0.2069e+00,
171  & 0.2069e+00,0.2069e+00,0.2069e+00,0.2069e+00,
172  & -0.6408e-01,-0.6408e-01,-0.6408e-01,-0.6408e-01,
173  & -0.6408e-01,-0.6408e-01,-0.6408e-01,-0.6408e-01,
174  & 0.2069e+00,0.2069e+00,0.2069e+00,0.2069e+00/
175 !
176 ! extrapolation from a 3x3x3=27 integration point scheme in a hex to
177 ! the all nodes in a 20-node element
178 !
179  data a27 /
180  & 2.37499,-0.12559,-0.16145,-0.12559,-0.12559,-0.16145, 0.11575,
181  & -0.16145, 0.32628, 0.11111, 0.11111, 0.32628, 0.11111,-0.10405,
182  & -0.10405, 0.11111, 0.32628, 0.11111,-0.10405, 0.11111,-0.31246,
183  & -0.31246, 0.31481, 0.31481, 0.31481, 0.31481,-0.16902,-0.16902,
184  & 1.28439,-0.27072,-0.19444,-0.27072,-0.19444, 0.15961,-0.00661,
185  & 0.15961,-0.27072,-0.27072, 0.15961, 0.15961,-0.12559, 2.37499,
186  & -0.12559,-0.16145,-0.16145,-0.12559,-0.16145, 0.11575, 0.32628,
187  & 0.32628, 0.11111, 0.11111, 0.11111, 0.11111,-0.10405,-0.10405,
188  & 0.11111, 0.32628, 0.11111,-0.10405,-0.31246, 0.31481, 0.31481,
189  & -0.31246, 0.31481,-0.16902,-0.16902, 0.31481,-0.27072,-0.19444,
190  & -0.27072, 1.28439, 0.15961,-0.00661, 0.15961,-0.19444,-0.27072,
191  & 0.15961, 0.15961,-0.27072,-0.48824,-0.48824,-0.48824,-0.48824,
192  & 0.22898, 0.22898, 0.22898, 0.22898, 0.05556, 0.05556, 0.05556,
193  & 0.05556, 0.05556, 0.05556, 0.05556, 0.05556,-0.22222,-0.22222,
194  & -0.22222,-0.22222, 0.31481,-0.31246,-0.31246, 0.31481,-0.16902,
195  & 0.31481, 0.31481,-0.16902,-0.27072, 1.28439,-0.27072,-0.19444,
196  & 0.15961,-0.19444, 0.15961,-0.00661, 0.15961,-0.27072,-0.27072,
197  & 0.15961,-0.12559,-0.16145,-0.12559, 2.37499,-0.16145, 0.11575,
198  & -0.16145,-0.12559, 0.11111, 0.11111, 0.32628, 0.32628,-0.10405,
199  & -0.10405, 0.11111, 0.11111, 0.11111,-0.10405, 0.11111, 0.32628,
200  & 0.31481, 0.31481,-0.31246,-0.31246,-0.16902,-0.16902, 0.31481,
201  & 0.31481,-0.19444,-0.27072, 1.28439,-0.27072,-0.00661, 0.15961,
202  & -0.19444, 0.15961, 0.15961, 0.15961,-0.27072,-0.27072,-0.16145,
203  & -0.12559, 2.37499,-0.12559, 0.11575,-0.16145,-0.12559,-0.16145,
204  & 0.11111, 0.32628, 0.32628, 0.11111,-0.10405, 0.11111, 0.11111,
205  & -0.10405,-0.10405, 0.11111, 0.32628, 0.11111,-0.31246, 0.31481,
206  & -0.16902, 0.31481,-0.31246, 0.31481,-0.16902, 0.31481,-0.27072,
207  & 0.15961, 0.15961,-0.27072,-0.27072, 0.15961, 0.15961,-0.27072,
208  & 1.28439,-0.19444,-0.00661,-0.19444,-0.48824,-0.48824, 0.22898,
209  & 0.22898,-0.48824,-0.48824, 0.22898, 0.22898, 0.05556,-0.22222,
210  & 0.05556,-0.22222, 0.05556,-0.22222, 0.05556,-0.22222, 0.05556,
211  & 0.05556, 0.05556, 0.05556, 0.31481,-0.31246, 0.31481,-0.16902,
212  & 0.31481,-0.31246, 0.31481,-0.16902,-0.27072,-0.27072, 0.15961,
213  & 0.15961,-0.27072,-0.27072, 0.15961, 0.15961,-0.19444, 1.28439,
214  & -0.19444,-0.00661,-0.48824, 0.22898, 0.22898,-0.48824,-0.48824,
215  & 0.22898, 0.22898,-0.48824,-0.22222, 0.05556,-0.22222, 0.05556,
216  & -0.22222, 0.05556,-0.22222, 0.05556, 0.05556, 0.05556, 0.05556,
217  & 0.05556,-0.29630,-0.29630,-0.29630,-0.29630,-0.29630,-0.29630,
218  & -0.29630,-0.29630,-0.11111,-0.11111,-0.11111,-0.11111,-0.11111,
219  & -0.11111,-0.11111,-0.11111,-0.11111,-0.11111,-0.11111,-0.11111,
220  & 0.22898,-0.48824,-0.48824, 0.22898, 0.22898,-0.48824,-0.48824,
221  & 0.22898,-0.22222, 0.05556,-0.22222, 0.05556,-0.22222, 0.05556,
222  & -0.22222, 0.05556, 0.05556, 0.05556, 0.05556, 0.05556, 0.31481,
223  & -0.16902, 0.31481,-0.31246, 0.31481,-0.16902, 0.31481,-0.31246,
224  & 0.15961, 0.15961,-0.27072,-0.27072, 0.15961, 0.15961,-0.27072,
225  & -0.27072,-0.19444,-0.00661,-0.19444, 1.28439, 0.22898, 0.22898,
226  & -0.48824,-0.48824, 0.22898, 0.22898,-0.48824,-0.48824, 0.05556,
227  & -0.22222, 0.05556,-0.22222, 0.05556,-0.22222, 0.05556,-0.22222,
228  & 0.05556, 0.05556, 0.05556, 0.05556,-0.16902, 0.31481,-0.31246,
229  & 0.31481,-0.16902, 0.31481,-0.31246, 0.31481, 0.15961,-0.27072,
230  & -0.27072, 0.15961, 0.15961,-0.27072,-0.27072, 0.15961,-0.00661,
231  & -0.19444, 1.28439,-0.19444,-0.12559,-0.16145, 0.11575,-0.16145,
232  & 2.37499,-0.12559,-0.16145,-0.12559, 0.11111,-0.10405,-0.10405,
233  & 0.11111, 0.32628, 0.11111, 0.11111, 0.32628, 0.32628, 0.11111,
234  & -0.10405, 0.11111, 0.31481, 0.31481,-0.16902,-0.16902,-0.31246,
235  & -0.31246, 0.31481, 0.31481,-0.19444, 0.15961,-0.00661, 0.15961,
236  & 1.28439,-0.27072,-0.19444,-0.27072,-0.27072,-0.27072, 0.15961,
237  & 0.15961,-0.16145,-0.12559,-0.16145, 0.11575,-0.12559, 2.37499,
238  & -0.12559,-0.16145, 0.11111, 0.11111,-0.10405,-0.10405, 0.32628,
239  & 0.32628, 0.11111, 0.11111, 0.11111, 0.32628, 0.11111,-0.10405,
240  & 0.31481,-0.16902,-0.16902, 0.31481,-0.31246, 0.31481, 0.31481,
241  & -0.31246, 0.15961,-0.00661, 0.15961,-0.19444,-0.27072,-0.19444,
242  & -0.27072, 1.28439,-0.27072, 0.15961, 0.15961,-0.27072, 0.22898,
243  & 0.22898, 0.22898, 0.22898,-0.48824,-0.48824,-0.48824,-0.48824,
244  & 0.05556, 0.05556, 0.05556, 0.05556, 0.05556, 0.05556, 0.05556,
245  & 0.05556,-0.22222,-0.22222,-0.22222,-0.22222,-0.16902, 0.31481,
246  & 0.31481,-0.16902, 0.31481,-0.31246,-0.31246, 0.31481, 0.15961,
247  & -0.19444, 0.15961,-0.00661,-0.27072, 1.28439,-0.27072,-0.19444,
248  & 0.15961,-0.27072,-0.27072, 0.15961,-0.16145, 0.11575,-0.16145,
249  & -0.12559,-0.12559,-0.16145,-0.12559, 2.37499,-0.10405,-0.10405,
250  & 0.11111, 0.11111, 0.11111, 0.11111, 0.32628, 0.32628, 0.11111,
251  & -0.10405, 0.11111, 0.32628,-0.16902,-0.16902, 0.31481, 0.31481,
252  & 0.31481, 0.31481,-0.31246,-0.31246,-0.00661, 0.15961,-0.19444,
253  & 0.15961,-0.19444,-0.27072, 1.28439,-0.27072, 0.15961, 0.15961,
254  & -0.27072,-0.27072, 0.11575,-0.16145,-0.12559,-0.16145,-0.16145,
255  & -0.12559, 2.37499,-0.12559,-0.10405, 0.11111, 0.11111,-0.10405,
256  & 0.11111, 0.32628, 0.32628, 0.11111,-0.10405, 0.11111, 0.32628,
257  & 0.11111/
258 !
259  data iflag /1/
260 !
261  do i=1,nk
262  inum(i)=0
263  enddo
264 !
265  do i=1,nk
266  do j=1,nfield
267  yn(j,i)=0.d0
268  enddo
269  enddo
270 !
271  do i=1,ne
272 !
273  if(ipkon(i).lt.0) cycle
274 !
275  indexe=ipkon(i)
276  lakonl=lakon(i)
277 !
278  if(lakonl(7:8).eq.'LC') then
279  nlayer=0
280  do j=1,mi(3)
281  if(ielmat(j,i).gt.0) then
282  nlayer=nlayer+1
283  else
284  exit
285  endif
286  enddo
287 !
288  if(lakonl(4:4).eq.'2') then
289  nopeexp=28
290  elseif(lakonl(4:5).eq.'15') then
291  nopeexp=21
292  endif
293  endif
294 !
295  if(lakonl(4:4).eq.'2') then
296  nope=20
297  elseif(lakonl(4:4).eq.'8') then
298  nope=8
299  elseif(lakonl(4:5).eq.'10') then
300  nope=10
301  elseif(lakonl(4:4).eq.'4') then
302  nope=4
303  elseif(lakonl(4:5).eq.'15') then
304  nope=15
305  elseif(lakonl(4:4).eq.'6') then
306  nope=6
307  elseif((lakonl(1:1).eq.'E').and.(lakonl(7:7).eq.'A'))then
308 !
309  inum(kon(indexe+1))=inum(kon(indexe+1))+1
310  inum(kon(indexe+2))=inum(kon(indexe+2))+1
311  cycle
312  elseif(lakonl(1:7).eq.'ESPRNGF') then
313  read(lakonl(8:8),'(i1)') nope
314  nope=nope+1
315  inum(kon(indexe+nope))=-1
316  cycle
317  elseif(lakonl(1:1).eq.'U') then
318  call extrapolate_u(yi,yn,ipkon,inum,kon,lakon,nfield,nk,
319  & ne,mi,ndim,orab,ielorien,co,iorienloc,cflag,
320  & vold,force,ielmat,thicke,ielprop,prop,i)
321  return
322  else
323  cycle
324  endif
325 !
326 ! storage in local coordinates
327 !
328 ! calculation of the integration point coordinates for
329 ! output in the local system
330 !
331  if((iorienloc.ne.0).and.
332  & ((lakonl(7:8).eq.'LC').or.(ielorien(1,i).ne.0))) then
333 !
334  if(lakonl(7:8).ne.'LC') then
335  iorien=ielorien(1,i)
336 !david -start
337  elseif(lakonl(4:5).eq.'20') then
338 !
339 ! composite materials
340 !
341  mint2d=4
342  nopes=8
343 !
344 ! determining the layer thickness and global thickness
345 ! at the shell integration points
346 !
347 ! xlayer: actual layer thickness
348 ! tlayer: total thickness of all layers
349 ! dlayer: total thickness of all layers up to the actual one
350 !
351  indexe=ipkon(i)
352  do kk=1,mint2d
353  xi=gauss3d2(1,kk)
354  et=gauss3d2(2,kk)
355  call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
356  tlayer(kk)=0.d0
357  do k=1,nlayer
358  thickness=0.d0
359  do j=1,nopes
360  thickness=thickness+thicke(k,indexe+j)*shp2(4,j)
361  enddo
362  tlayer(kk)=tlayer(kk)+thickness
363  xlayer(k,kk)=thickness
364  enddo
365  enddo
366 !
367  ilayer=0
368  do k=1,mint2d
369  dlayer(k)=0.d0
370  enddo
371 !
372 ! S6-composite shell
373 !
374  elseif(lakonl(4:5).eq.'15') then
375 !
376 ! composite materials
377 !
378  mint2d=3
379  nopes=6
380 !
381 ! determining the layer thickness and global thickness
382 ! at the shell integration points
383 !
384 ! xlayer: actual layer thickness
385 ! tlayer: total thickness of all layers
386 ! dlayer: total thickness of all layers up to the actual one
387 !
388  indexe=ipkon(i)
389  do kk=1,mint2d
390  xi=gauss3d10(1,kk)
391  et=gauss3d10(2,kk)
392  call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
393  tlayer(kk)=0.d0
394  do k=1,nlayer
395  thickness=0.d0
396  do j=1,nopes
397  thickness=thickness+thicke(k,indexe+j)*shp2(4,j)
398  enddo
399  tlayer(kk)=tlayer(kk)+thickness
400  xlayer(k,kk)=thickness
401  enddo
402  enddo
403 !
404  ilayer=0
405  do k=1,mint2d
406  dlayer(k)=0.d0
407  enddo
408  endif
409 !
410  if((lakon(i)(4:5).eq.'8R').or.(lakon(i)(1:1).eq.'F')) then
411  mint3d=1
412  elseif((lakon(i)(4:7).eq.'20RB').and.
413  & (lakon(i)(8:8).ne.'R').and.
414  & (lakon(i)(8:8).ne.'C')) then
415  call beamintscheme(lakon(i),mint3d,ielprop(i),prop,
416  & null,xi,et,ze,weight)
417  elseif((lakon(i)(4:4).eq.'8').or.
418  & (lakon(i)(4:6).eq.'20R')) then
419  if(lakonl(7:8).eq.'LC') then
420  mint3d=8*nlayer
421  else
422  mint3d=8
423  endif
424  elseif(lakon(i)(4:4).eq.'2') then
425  mint3d=27
426  elseif((lakon(i)(4:5).eq.'10')) then
427  mint3d=4
428  elseif(lakon(i)(4:4).eq.'4') then
429  mint3d=1
430  elseif(lakon(i)(4:5).eq.'15') then
431  if(lakonl(7:8).eq.'LC') then
432  mint3d=6*nlayer
433  else
434  mint3d=9
435  endif
436  elseif(lakon(i)(4:4).eq.'6') then
437  mint3d=2
438  endif
439 !
440  do j=1,nope
441  konl=kon(indexe+j)
442  do k=1,3
443  xl(k,j)=co(k,konl)
444  enddo
445  enddo
446 !
447  do j=1,mint3d
448  if((lakon(i)(4:5).eq.'8R').or.
449  & (lakon(i)(1:4).eq.'F3D8')) then
450  xi=gauss3d1(1,j)
451  et=gauss3d1(2,j)
452  ze=gauss3d1(3,j)
453  weight=weight3d1(j)
454  elseif((lakon(i)(4:7).eq.'20RB').and.
455  & (lakon(i)(8:8).ne.'R').and.
456  & (lakon(i)(8:8).ne.'C')) then
457  call beamintscheme(lakon(i),mint3d,ielprop(i),prop,
458  & j,xi,et,ze,weight)
459  elseif((lakon(i)(4:4).eq.'8').or.
460  & (lakon(i)(4:6).eq.'20R'))
461  & then
462  if(lakonl(7:8).ne.'LC') then
463  xi=gauss3d2(1,j)
464  et=gauss3d2(2,j)
465  ze=gauss3d2(3,j)
466  weight=weight3d2(j)
467  else
468 !
469 ! kl: number of the integration point within the layer
470 !
471  kl=mod(j,8)
472  if(kl.eq.0) kl=8
473 !
474  xi=gauss3d2(1,kl)
475  et=gauss3d2(2,kl)
476  ze=gauss3d2(3,kl)
477  weight=weight3d2(kl)
478 !
479 ! ki: position of the integration point (1...4)
480 !
481  ki=mod(j,4)
482  if(ki.eq.0) ki=4
483 !
484  if(kl.eq.1) then
485  ilayer=ilayer+1
486  if(ilayer.gt.1) then
487  do k=1,4
488  dlayer(k)=dlayer(k)+xlayer(ilayer-1,k)
489  enddo
490  endif
491  endif
492  ze=2.d0*(dlayer(ki)+(ze+1.d0)/2.d0*
493  & xlayer(ilayer,ki))/tlayer(ki)-1.d0
494  weight=weight*xlayer(ilayer,ki)/tlayer(ki)
495 !
496 ! material and orientation
497 !
498  iorien=ielorien(ilayer,i)
499  endif
500  elseif(lakon(i)(4:4).eq.'2') then
501  xi=gauss3d3(1,j)
502  et=gauss3d3(2,j)
503  ze=gauss3d3(3,j)
504  weight=weight3d3(j)
505  elseif(lakon(i)(4:5).eq.'10') then
506  xi=gauss3d5(1,j)
507  et=gauss3d5(2,j)
508  ze=gauss3d5(3,j)
509  weight=weight3d5(j)
510  elseif(lakon(i)(4:4).eq.'4') then
511  xi=gauss3d4(1,j)
512  et=gauss3d4(2,j)
513  ze=gauss3d4(3,j)
514  weight=weight3d4(j)
515  elseif(lakon(i)(4:5).eq.'15') then
516  if(lakonl(7:8).ne.'LC') then
517  xi=gauss3d8(1,j)
518  et=gauss3d8(2,j)
519  ze=gauss3d8(3,j)
520  weight=weight3d8(j)
521  else
522 !
523 ! kl: number of the integration point within the layer
524 !
525  kl=mod(j,6)
526  if(kl.eq.0) kl=6
527 !
528  xi=gauss3d10(1,kl)
529  et=gauss3d10(2,kl)
530  ze=gauss3d10(3,kl)
531  weight=weight3d10(kl)
532 !
533 ! ki: position of the integration point (1...3)
534 !
535  ki=mod(j,3)
536  if(ki.eq.0) ki=3
537 !
538  if(kl.eq.1) then
539  ilayer=ilayer+1
540  if(ilayer.gt.1) then
541  do k=1,3
542  dlayer(k)=dlayer(k)+xlayer(ilayer-1,k)
543  enddo
544  endif
545  endif
546  ze=2.d0*(dlayer(ki)+(ze+1.d0)/2.d0*
547  & xlayer(ilayer,ki))/tlayer(ki)-1.d0
548  weight=weight*xlayer(ilayer,ki)/tlayer(ki)
549 !
550 ! material and orientation
551 !
552  iorien=ielorien(ilayer,i)
553  endif
554  elseif(lakon(i)(1:4).eq.'C3D6') then
555  xi=gauss3d7(1,j)
556  et=gauss3d7(2,j)
557  ze=gauss3d7(3,j)
558  weight=weight3d7(j)
559  elseif(lakon(i)(1:4).eq.'F3D6') then
560  xi=gauss3d14(1,j)
561  et=gauss3d14(2,j)
562  ze=gauss3d14(3,j)
563  weight=weight3d14(j)
564  endif
565 !
566  if(nope.eq.20) then
567  call shape20h(xi,et,ze,xl,xsj,shp,iflag)
568  elseif(nope.eq.8) then
569  call shape8h(xi,et,ze,xl,xsj,shp,iflag)
570  elseif(nope.eq.10) then
571  call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
572  elseif(nope.eq.4) then
573  call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
574  elseif(nope.eq.15) then
575  call shape15w(xi,et,ze,xl,xsj,shp,iflag)
576  else
577  call shape6w(xi,et,ze,xl,xsj,shp,iflag)
578  endif
579 !
580 ! layer without orientation in a composite
581 !
582  if(iorien.eq.0) then
583  if(nfield.eq.3) then
584  do k=1,3
585  yiloc(k,j)=yi(k,j,i)
586  enddo
587  elseif(nfield.eq.6) then
588  do k=1,6
589  yiloc(k,j)=yi(k,j,i)
590  enddo
591  endif
592  cycle
593  endif
594 !
595 ! coordinates of the integration point
596 !
597  do k=1,3
598  coords(k,j)=0.d0
599  do l=1,nope
600  coords(k,j)=coords(k,j)+xl(k,l)*shp(4,l)
601  enddo
602  enddo
603 !
604 ! transforming the vector or tensor field
605 !
606  if(nfield.eq.3) then
607  call transformatrix(orab(1,iorien),coords(1,j),a)
608  yiloc(1,j)=yi(1,j,i)*a(1,1)+yi(2,j,i)*a(2,1)+
609  & yi(3,j,i)*a(3,1)
610  yiloc(2,j)=yi(1,j,i)*a(1,2)+yi(2,j,i)*a(2,2)+
611  & yi(3,j,i)*a(3,2)
612  yiloc(3,j)=yi(1,j,i)*a(1,3)+yi(2,j,i)*a(2,3)+
613  & yi(3,j,i)*a(3,3)
614  elseif(nfield.eq.6) then
615  call transformatrix(orab(1,iorien),coords(1,j),a)
616  b(1,1)=yi(1,j,i)
617  b(2,2)=yi(2,j,i)
618  b(3,3)=yi(3,j,i)
619  b(1,2)=yi(4,j,i)
620  b(1,3)=yi(5,j,i)
621  b(2,3)=yi(6,j,i)
622  b(2,1)=b(1,2)
623  b(3,1)=b(1,3)
624  b(3,2)=b(2,3)
625  do k=1,3
626  do l=1,3
627  c(k,l)=0.d0
628  do m=1,3
629  c(k,l)=c(k,l)+b(k,m)*a(m,l)
630  enddo
631  enddo
632  enddo
633  do k=1,3
634  do l=k,3
635  b(k,l)=0.d0
636  do m=1,3
637  b(k,l)=b(k,l)+a(m,k)*c(m,l)
638  enddo
639  enddo
640  enddo
641  yiloc(1,j)=b(1,1)
642  yiloc(2,j)=b(2,2)
643  yiloc(3,j)=b(3,3)
644  yiloc(4,j)=b(1,2)
645  yiloc(5,j)=b(1,3)
646  yiloc(6,j)=b(2,3)
647  endif
648  enddo
649 !
650  if(lakonl(1:1).eq.'F') then
651  do j=1,8
652  do k=1,nfield
653  field(k,j)=yiloc(k,1)
654  enddo
655  enddo
656  elseif((lakonl(4:7).eq.'20RB').and.
657  & (lakonl(8:8).ne.'R').and.
658  & (lakonl(8:8).ne.'C')) then
659  call beamextscheme(yi(1,1,i),ndim,nfield,lakonl,
660  & ielprop(i),prop,field,mi)
661 c Bernhardi start
662  elseif((lakonl(4:6).eq.'20R').or.
663  & (lakonl(4:5).eq.'8 ').or.(lakonl(4:5).eq.'8I')) then
664 c Bernhardi end
665  if(lakonl(7:8).ne.'LC') then
666  do j=1,8
667  do k=1,nfield
668  field(k,j)=0.d0
669  do l=1,8
670  field(k,j)=field(k,j)+a8(j,l)*yiloc(k,l)
671  enddo
672  enddo
673  enddo
674  else
675  do m=1,nlayer
676  jj=20*(m-1)
677  ll=8*(m-1)
678  do j=1,8
679  do k=1,nfield
680  field(k,jj+j)=0.d0
681  do l=1,8
682  field(k,jj+j)=
683  & field(k,jj+j)+a8(j,l)*yiloc(k,ll+l)
684  enddo
685  enddo
686  enddo
687  enddo
688  endif
689  elseif(lakonl(4:4).eq.'8') then
690  do j=1,8
691  do k=1,nfield
692  field(k,j)=yiloc(k,1)
693  enddo
694  enddo
695  elseif(lakonl(4:5).eq.'10') then
696  do j=1,4
697  do k=1,nfield
698  field(k,j)=0.d0
699  do l=1,4
700  field(k,j)=field(k,j)+a4(j,l)*yiloc(k,l)
701  enddo
702  enddo
703  enddo
704  elseif(lakonl(4:4).eq.'2') then
705  do j=1,20
706  do k=1,nfield
707  field(k,j)=0.d0
708  do l=1,27
709  field(k,j)=field(k,j)+a27(j,l)*yiloc(k,l)
710  enddo
711  enddo
712  enddo
713  elseif(lakonl(4:4).eq.'4') then
714  do j=1,4
715  do k=1,nfield
716  field(k,j)=yiloc(k,1)
717  enddo
718  enddo
719  elseif(lakonl(4:4).eq.'1') then
720  if(lakonl(7:8).ne.'LC') then
721  do j=1,6
722  do k=1,nfield
723  field(k,j)=0.d0
724  do l=1,9
725  field(k,j)=field(k,j)+a9(j,l)*yiloc(k,l)
726  enddo
727  enddo
728  enddo
729  else
730  do m=1,nlayer
731  jj=15*(m-1)
732  ll=6*(m-1)
733  do j=1,6
734  do k=1,nfield
735  field(k,jj+j)=0.d0
736  do l=1,6
737  field(k,jj+j)=
738  & field(k,jj+j)+a6(j,l)*yiloc(k,ll+l)
739  enddo
740  enddo
741  enddo
742  enddo
743  endif
744  else
745  do j=1,6
746  do k=1,nfield
747  field(k,j)=0.d0
748  do l=1,2
749  field(k,j)=field(k,j)+a2(j,l)*yiloc(k,l)
750  enddo
751  enddo
752  enddo
753  endif
754  else
755 !
756 ! storage in global coordinates
757 !
758 ! determining the field values in the vertex nodes
759 ! for C3D20R and C3D8: trilinear extrapolation (= use of the
760 ! C3D8 linear interpolation functions)
761 ! for C3D8R: constant field value in each element
762 ! for C3D10: use of the C3D4 linear interpolation functions
763 ! for C3D4: constant field value in each element
764 ! for C3D15: use of the C3D6 linear interpolation functions
765 ! for C3D6: use of a linear interpolation function
766 !
767  if(lakonl(1:1).eq.'F') then
768  do j=1,8
769  do k=1,nfield
770  field(k,j)=yi(k,1,i)
771  enddo
772  enddo
773  elseif((lakonl(4:7).eq.'20RB').and.
774  & (lakonl(8:8).ne.'R').and.
775  & (lakonl(8:8).ne.'C')) then
776  call beamextscheme(yi(1,1,i),ndim,nfield,lakonl,
777  & ielprop(i),prop,field,mi)
778 !
779 c Bernhardi start
780  elseif((lakonl(4:6).eq.'20R').or.
781  & (lakonl(4:5).eq.'8 ').or.(lakonl(4:5).eq.'8I')) then
782 c Bernhardi end
783  if(lakonl(7:8).ne.'LC') then
784  do j=1,8
785  do k=1,nfield
786  field(k,j)=0.d0
787  do l=1,8
788  field(k,j)=field(k,j)+a8(j,l)*yi(k,l,i)
789  enddo
790  enddo
791  enddo
792  else
793  do m=1,nlayer
794  jj=20*(m-1)
795  ll=8*(m-1)
796  do j=1,8
797  do k=1,nfield
798  field(k,jj+j)=0.d0
799  do l=1,8
800  field(k,jj+j)=
801  & field(k,jj+j)+a8(j,l)*yi(k,ll+l,i)
802  enddo
803  enddo
804  enddo
805  enddo
806  endif
807  elseif(lakonl(4:4).eq.'8') then
808  do j=1,8
809  do k=1,nfield
810  field(k,j)=yi(k,1,i)
811  enddo
812  enddo
813  elseif(lakonl(4:5).eq.'10') then
814  do j=1,4
815  do k=1,nfield
816  field(k,j)=0.d0
817  do l=1,4
818  field(k,j)=field(k,j)+a4(j,l)*yi(k,l,i)
819  enddo
820  enddo
821  enddo
822  elseif(lakonl(4:4).eq.'2') then
823  do j=1,20
824  do k=1,nfield
825  field(k,j)=0.d0
826  do l=1,27
827  field(k,j)=field(k,j)+a27(j,l)*yi(k,l,i)
828  enddo
829  enddo
830  enddo
831  elseif(lakonl(4:4).eq.'4') then
832  do j=1,4
833  do k=1,nfield
834  field(k,j)=yi(k,1,i)
835  enddo
836  enddo
837  elseif(lakonl(4:4).eq.'1') then
838  if(lakonl(7:8).ne.'LC') then
839  do j=1,6
840  do k=1,nfield
841  field(k,j)=0.d0
842  do l=1,9
843  field(k,j)=field(k,j)+a9(j,l)*yi(k,l,i)
844  enddo
845  enddo
846  enddo
847  else
848  do m=1,nlayer
849  jj=15*(m-1)
850  ll=6*(m-1)
851  do j=1,6
852  do k=1,nfield
853  field(k,jj+j)=0.d0
854  do l=1,6
855  field(k,jj+j)=
856  & field(k,jj+j)+a6(j,l)*yi(k,ll+l,i)
857  enddo
858  enddo
859  enddo
860  enddo
861  endif
862  else
863  do j=1,6
864  do k=1,nfield
865  field(k,j)=0.d0
866  do l=1,2
867  field(k,j)=field(k,j)+a2(j,l)*yi(k,l,i)
868  enddo
869  enddo
870  enddo
871  endif
872  endif
873 !
874 ! determining the field values in the midside nodes
875 !
876  if(lakonl(4:6).eq.'20R') then
877  if(lakonl(7:8).ne.'LC') then
878  do j=9,20
879  do k=1,nfield
880  field(k,j)=(field(k,nonei20(2,j-8))+
881  & field(k,nonei20(3,j-8)))/2.d0
882  enddo
883  enddo
884  else
885  do m=1,nlayer
886  jj=20*(m-1)
887  do j=9,20
888  do k=1,nfield
889  field(k,jj+j)=(field(k,jj+nonei20(2,j-8))
890  & +field(k,jj+nonei20(3,j-8)))/2.d0
891  enddo
892  enddo
893  enddo
894  endif
895  elseif(lakonl(4:5).eq.'10') then
896  do j=5,10
897  do k=1,nfield
898  field(k,j)=(field(k,nonei10(2,j-4))+
899  & field(k,nonei10(3,j-4)))/2.d0
900  enddo
901  enddo
902  elseif(lakonl(4:5).eq.'15') then
903  if(lakonl(7:8).ne.'LC') then
904  do j=7,15
905  do k=1,nfield
906  field(k,j)=(field(k,nonei15(2,j-6))+
907  & field(k,nonei15(3,j-6)))/2.d0
908  enddo
909  enddo
910  else
911  do m=1,nlayer
912  jj=15*(m-1)
913  do j=7,15
914  do k=1,nfield
915  field(k,jj+j)=(field(k,jj+nonei15(2,j-6))
916  & +field(k,jj+nonei15(3,j-6)))/2.d0
917  enddo
918  enddo
919  enddo
920  endif
921  endif
922 !
923 ! transferring the field values into yn
924 !
925  if(lakonl(7:8).ne.'LC') then
926  do j=1,nope
927  do k=1,nfield
928  yn(k,kon(indexe+j))=yn(k,kon(indexe+j))+
929  & field(k,j)
930  enddo
931  inum(kon(indexe+j))=inum(kon(indexe+j))+1
932  enddo
933  else
934  do j=1,nope*nlayer
935  do k=1,nfield
936  yn(k,kon(indexe+nopeexp+j))=
937  & yn(k,kon(indexe+nopeexp+j))+field(k,j)
938  enddo
939  inum(kon(indexe+nopeexp+j))=
940  & inum(kon(indexe+nopeexp+j))+1
941  enddo
942  endif
943 !
944 c Bernhardi start
945 c incompatible modes elements
946  if(lakonl(1:5).eq.'C3D8I') then
947  do j=1,3
948  do k=1,nfield
949  yn(k,kon(indexe+nope+j))=0.0d0
950  enddo
951  enddo
952  endif
953 c Bernhardi end
954 !
955  enddo
956 !
957 ! taking the mean of nodal contributions coming from different
958 ! elements having the node in common
959 !
960  do i=1,nk
961  if(inum(i).gt.0) then
962  do j=1,nfield
963  yn(j,i)=yn(j,i)/inum(i)
964  enddo
965  endif
966  enddo
967 !
968 ! for 1d and 2d elements only:
969 ! finding the solution in the original nodes
970 !
971  if((cflag.ne.' ').and.(cflag.ne.'E')) then
972  call map3dto1d2d(yn,ipkon,inum,kon,lakon,nfield,nk,ne,cflag,
973  & co,vold,force,mi)
974  endif
975 !
976  return
shape6w
subroutine shape6w(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape6w.f:20
shape8q
subroutine shape8q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape8q.f:20
beamextscheme
subroutine beamextscheme(yil, ndim, nfield, lakonl, npropstart, prop, field, mi)
Definition: beamextscheme.f:21
shape10tet
subroutine shape10tet(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape10tet.f:20
shape8h
subroutine shape8h(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape8h.f:20
extrapolate_u
subroutine extrapolate_u(yi, yn, ipkon, inum, kon, lakon, nfield, nk, ne, mi, ndim, orab, ielorien, co, iorienloc, cflag, vold, force, ielmat, thicke, ielprop, prop, i)
Definition: extrapolate_u.f:22
shape15w
subroutine shape15w(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape15w.f:20
transformatrix
subroutine transformatrix(xab, p, a)
Definition: transformatrix.f:20
shape20h
subroutine shape20h(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape20h.f:20
beamintscheme
subroutine beamintscheme(lakonl, mint3d, npropstart, prop, kk, xi, et, ze, weight)
Definition: beamintscheme.f:21
map3dto1d2d
subroutine map3dto1d2d(yn, ipkon, inum, kon, lakon, nfield, nk, ne, cflag, co, vold, force, mi)
Definition: map3dto1d2d.f:21
shape4tet
subroutine shape4tet(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape4tet.f:20
shape6tri
subroutine shape6tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape6tri.f:20
thickness
subroutine thickness(dgdx, nobject, nodedesiboun, ndesiboun, objectset, xo, yo, zo, x, y, z, nx, ny, nz, co, ifree, ndesia, ndesib, iobject, ndesi, dgdxglob, nk)
Definition: thickness.f:22
Hosted by OpenAircraft.com, (Michigan UAV, LLC)