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

Go to the source code of this file.

Functions/Subroutines

subroutine materialdata_me (elcon, nelcon, rhcon, nrhcon, alcon, nalcon, imat, amat, iorien, pgauss, orab, ntmat_, elas, rho, iel, ithermal, alzero, mattyp, t0l, t1l, ihyper, istiff, elconloc, eth, kode, plicon, nplicon, plkcon, nplkcon, npmat_, plconloc, mi, dtime, iint, xstiff, ncmat_)
 

Function/Subroutine Documentation

◆ materialdata_me()

subroutine materialdata_me ( real*8, dimension(0:ncmat_,ntmat_,*), intent(in)  elcon,
integer, dimension(2,*), intent(in)  nelcon,
real*8, dimension(0:1,ntmat_,*), intent(in)  rhcon,
integer, dimension(*), intent(in)  nrhcon,
real*8, dimension(0:6,ntmat_,*), intent(in)  alcon,
integer, dimension(2,*), intent(in)  nalcon,
integer, intent(in)  imat,
character*80, intent(in)  amat,
integer, intent(in)  iorien,
real*8, dimension(3), intent(in)  pgauss,
real*8, dimension(7,*), intent(in)  orab,
integer, intent(in)  ntmat_,
real*8, dimension(21), intent(inout)  elas,
real*8, intent(inout)  rho,
integer, intent(in)  iel,
integer, intent(in)  ithermal,
real*8, dimension(*), intent(in)  alzero,
integer, intent(inout)  mattyp,
real*8, intent(in)  t0l,
real*8, intent(in)  t1l,
integer, intent(in)  ihyper,
integer, intent(in)  istiff,
real*8, dimension(21), intent(inout)  elconloc,
real*8, dimension(6), intent(inout)  eth,
integer, intent(in)  kode,
real*8, dimension(0:2*npmat_,ntmat_,*), intent(in)  plicon,
integer, dimension(0:ntmat_,*), intent(in)  nplicon,
real*8, dimension(0:2*npmat_,ntmat_,*), intent(in)  plkcon,
integer, dimension(0:ntmat_,*), intent(in)  nplkcon,
integer, intent(in)  npmat_,
real*8, dimension(802), intent(inout)  plconloc,
integer, dimension(*), intent(in)  mi,
real*8, intent(in)  dtime,
integer, intent(in)  iint,
real*8, dimension(27,mi(1),*), intent(in)  xstiff,
integer, intent(in)  ncmat_ 
)
23 !
24  implicit none
25 !
26 ! determines the material data for element iel
27 !
28 ! istiff=0: only interpolation of material data
29 ! istiff=1: copy the consistent tangent matrix from the field
30 ! xstiff and check for zero entries
31 !
32  character*80 amat
33 !
34  integer nelcon(2,*),nrhcon(*),nalcon(2,*),
35  & imat,iorien,ithermal,j,k,mattyp,kal(2,6),j1,j2,j3,j4,
36  & jj,ntmat_,istiff,nelconst,ihyper,kode,itemp,kin,nelas,
37  & iel,iint,mi(*),ncmat_,id,two,seven,
38  & nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),npmat_
39 !
40  real*8 elcon(0:ncmat_,ntmat_,*),rhcon(0:1,ntmat_,*),
41  & alcon(0:6,ntmat_,*),eth(6),xstiff(27,mi(1),*),
42  & orab(7,*),elas(21),alph(6),alzero(*),rho,t0l,t1l,
43  & skl(3,3),xa(3,3),elconloc(21),emax,pgauss(3),
44  & plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*),
45  & plconloc(802),dtime
46 !
47  intent(in) elcon,nelcon,rhcon,nrhcon,alcon,nalcon,
48  & imat,amat,iorien,pgauss,orab,ntmat_,iel,ithermal,
49  & alzero,t0l,t1l,ihyper,istiff,kode,plicon,
50  & nplicon,plkcon,nplkcon,npmat_,mi,dtime,iint,
51  & xstiff,ncmat_
52 !
53  intent(inout) plconloc,eth,elconloc,elas,mattyp,rho
54 !
55  kal=reshape((/1,1,2,2,3,3,1,2,1,3,2,3/),(/2,6/))
56 !
57  two=2
58  seven=7
59 !
60 ! nelconst: # constants read from file
61 ! nelas: # constants in the local tangent stiffness matrix
62 !
63  if(istiff.eq.1) then
64 
65  nelas=nelcon(1,imat)
66  if((nelas.lt.0).or.((nelas.ne.2).and.(iorien.ne.0))) nelas=21
67 !
68 ! calculating the density (needed for the mass matrix and
69 ! gravity or centrifugal loading)
70 !
71  if(ithermal.eq.0) then
72  rho=rhcon(1,1,imat)
73  else
74  call ident2(rhcon(0,1,imat),t1l,nrhcon(imat),two,id)
75  if(nrhcon(imat).eq.0) then
76  continue
77  elseif(nrhcon(imat).eq.1) then
78  rho=rhcon(1,1,imat)
79  elseif(id.eq.0) then
80  rho=rhcon(1,1,imat)
81  elseif(id.eq.nrhcon(imat)) then
82  rho=rhcon(1,id,imat)
83  else
84  rho=rhcon(1,id,imat)+
85  & (rhcon(1,id+1,imat)-rhcon(1,id,imat))*
86  & (t1l-rhcon(0,id,imat))/
87  & (rhcon(0,id+1,imat)-rhcon(0,id,imat))
88  endif
89  endif
90 !
91 ! for nonlinear behavior (nonlinear geometric or
92 ! nonlinear material behavior): copy the stiffness matrix
93 ! from the last stress calculation
94 !
95  do j=1,21
96  elas(j)=xstiff(j,iint,iel)
97  enddo
98 !
99 ! check whether the fully anisotropic case can be
100 ! considered as orthotropic
101 !
102  if(nelas.eq.21) then
103  emax=0.d0
104  do j=1,9
105  emax=max(emax,dabs(elas(j)))
106  enddo
107  do j=10,21
108  if(dabs(elas(j)).gt.emax*1.d-10) then
109  emax=-1.d0
110  exit
111  endif
112  enddo
113  if(emax.gt.0.d0) nelas=9
114  endif
115 !
116 ! determining the type: isotropic, orthotropic or anisotropic
117 !
118  if(nelas.le.2) then
119  mattyp=1
120  elseif(nelas.le.9) then
121  mattyp=2
122  else
123  mattyp=3
124  endif
125 !
126  else
127 !
128  nelconst=nelcon(1,imat)
129 !
130  if(nelconst.lt.0) then
131 !
132 ! inelastic material or user material
133 !
134  if(nelconst.eq.-1) then
135  nelconst=3
136  elseif(nelconst.eq.-2) then
137  nelconst=3
138  elseif(nelconst.eq.-3) then
139  nelconst=2
140  elseif(nelconst.eq.-4) then
141  nelconst=3
142  elseif(nelconst.eq.-5) then
143  nelconst=6
144  elseif(nelconst.eq.-6) then
145  nelconst=9
146  elseif(nelconst.eq.-7) then
147  nelconst=3
148  elseif(nelconst.eq.-8) then
149  nelconst=7
150  elseif(nelconst.eq.-9) then
151  nelconst=12
152  elseif(nelconst.eq.-10) then
153  nelconst=2
154  elseif(nelconst.eq.-11) then
155  nelconst=4
156  elseif(nelconst.eq.-12) then
157  nelconst=6
158  elseif(nelconst.eq.-13) then
159  nelconst=5
160  elseif(nelconst.eq.-14) then
161  nelconst=6
162  elseif(nelconst.eq.-15) then
163  nelconst=3
164  elseif(nelconst.eq.-16) then
165  nelconst=6
166  elseif(nelconst.eq.-17) then
167  nelconst=9
168  elseif(nelconst.eq.-50) then
169  nelconst=5
170  elseif(nelconst.eq.-51) then
171  nelconst=2
172  elseif(nelconst.eq.-52) then
173  nelconst=5
174  elseif(nelconst.le.-100) then
175  nelconst=-nelconst-100
176  endif
177 !
178  endif
179 !
180 ! in case no initial temperatures are defined, the calculation
181 ! is assumed athermal, and the first available set material
182 ! constants are used
183 !
184  if(ithermal.eq.0) then
185  if(ihyper.ne.1) then
186  do k=1,nelconst
187  elconloc(k)=elcon(k,1,imat)
188  enddo
189  else
190  do k=1,nelconst
191  elconloc(k)=elcon(k,1,imat)
192  enddo
193 !
194  itemp=1
195 !
196  if((kode.lt.-50).and.(kode.gt.-100)) then
197  plconloc(1)=0.d0
198  plconloc(2)=0.d0
199  plconloc(3)=0.d0
200  plconloc(801)=nplicon(1,imat)+0.5d0
201  plconloc(802)=nplkcon(1,imat)+0.5d0
202 !
203 ! isotropic hardening
204 !
205  if(nplicon(1,imat).ne.0) then
206  kin=0
207  call plcopy(plicon,nplicon,plconloc,npmat_,ntmat_,
208  & imat,itemp,iel,kin)
209  endif
210 !
211 ! kinematic hardening
212 !
213  if(nplkcon(1,imat).ne.0) then
214  kin=1
215  call plcopy(plkcon,nplkcon,plconloc,npmat_,ntmat_,
216  & imat,itemp,iel,kin)
217  endif
218 !
219  endif
220 !
221  endif
222  else
223 !
224 ! calculating the expansion coefficients
225 !
226  call ident2(alcon(0,1,imat),t1l,nalcon(2,imat),seven,id)
227  if(nalcon(2,imat).eq.0) then
228  do k=1,6
229  alph(k)=0.d0
230  enddo
231  continue
232  elseif(nalcon(2,imat).eq.1) then
233  do k=1,nalcon(1,imat)
234  alph(k)=alcon(k,1,imat)*(t1l-alzero(imat))
235  enddo
236  elseif(id.eq.0) then
237  do k=1,nalcon(1,imat)
238  alph(k)=alcon(k,1,imat)*(t1l-alzero(imat))
239  enddo
240  elseif(id.eq.nalcon(2,imat)) then
241  do k=1,nalcon(1,imat)
242  alph(k)=alcon(k,id,imat)*(t1l-alzero(imat))
243  enddo
244  else
245  do k=1,nalcon(1,imat)
246  alph(k)=(alcon(k,id,imat)+
247  & (alcon(k,id+1,imat)-alcon(k,id,imat))*
248  & (t1l-alcon(0,id,imat))/
249  & (alcon(0,id+1,imat)-alcon(0,id,imat)))
250  & *(t1l-alzero(imat))
251  enddo
252  endif
253 !
254 ! subtracting the initial temperature influence
255 !
256  call ident2(alcon(0,1,imat),t0l,nalcon(2,imat),seven,id)
257  if(nalcon(2,imat).eq.0) then
258  continue
259  elseif(nalcon(2,imat).eq.1) then
260  do k=1,nalcon(1,imat)
261  alph(k)=alph(k)-alcon(k,1,imat)*(t0l-alzero(imat))
262  enddo
263  elseif(id.eq.0) then
264  do k=1,nalcon(1,imat)
265  alph(k)=alph(k)-alcon(k,1,imat)*(t0l-alzero(imat))
266  enddo
267  elseif(id.eq.nalcon(2,imat)) then
268  do k=1,nalcon(1,imat)
269  alph(k)=alph(k)-alcon(k,id,imat)*(t0l-alzero(imat))
270  enddo
271  else
272  do k=1,nalcon(1,imat)
273  alph(k)=alph(k)-(alcon(k,id,imat)+
274  & (alcon(k,id+1,imat)-alcon(k,id,imat))*
275  & (t0l-alcon(0,id,imat))/
276  & (alcon(0,id+1,imat)-alcon(0,id,imat)))
277  & *(t0l-alzero(imat))
278  enddo
279  endif
280 !
281 ! storing the thermal strains
282 !
283  if(nalcon(1,imat).eq.1) then
284  do k=1,3
285  eth(k)=alph(1)
286  enddo
287  do k=4,6
288  eth(k)=0.d0
289  enddo
290  elseif(nalcon(1,imat).eq.3) then
291  do k=1,3
292  eth(k)=alph(k)
293  enddo
294  do k=4,6
295  eth(k)=0.d0
296  enddo
297  else
298  do k=1,6
299  eth(k)=alph(k)
300  enddo
301  endif
302 !
303 ! calculating the hardening coefficients
304 !
305 ! for the calculation of the stresses, the whole curve
306 ! has to be stored:
307 ! plconloc(2*k-1), k=1...20: equivalent plastic strain values (iso)
308 ! plconloc(2*k),k=1...20: corresponding stresses (iso)
309 ! plconloc(39+2*k),k=1...20: equivalent plastic strain values (kin)
310 ! plconloc(40+2*k),k=1...20: corresponding stresses (kin)
311 !
312 ! initialization
313 !
314  if((kode.lt.-50).and.(kode.gt.-100)) then
315  if(npmat_.eq.0) then
316  plconloc(801)=0.5d0
317  plconloc(802)=0.5d0
318  else
319  plconloc(1)=0.d0
320  plconloc(2)=0.d0
321  plconloc(3)=0.d0
322  plconloc(801)=nplicon(1,imat)+0.5d0
323  plconloc(802)=nplkcon(1,imat)+0.5d0
324 !
325 ! isotropic hardening
326 !
327  if(nplicon(1,imat).ne.0) then
328 !
329  if(nplicon(0,imat).eq.1) then
330  id=-1
331  else
332  call ident2(plicon(0,1,imat),t1l,
333  & nplicon(0,imat),2*npmat_+1,id)
334  endif
335 !
336  if(nplicon(0,imat).eq.0) then
337  continue
338  elseif((nplicon(0,imat).eq.1).or.(id.eq.0).or.
339  & (id.eq.nplicon(0,imat))) then
340  if(id.le.0) then
341  itemp=1
342  else
343  itemp=id
344  endif
345  kin=0
346  call plcopy(plicon,nplicon,plconloc,npmat_,
347  & ntmat_,imat,itemp,iel,kin)
348  if((id.eq.0).or.(id.eq.nplicon(0,imat))) then
349  endif
350  else
351  kin=0
352  call plmix(plicon,nplicon,plconloc,npmat_,
353  & ntmat_,imat,id+1,t1l,iel,kin)
354  endif
355  endif
356 !
357 ! kinematic hardening
358 !
359  if(nplkcon(1,imat).ne.0) then
360 !
361  if(nplkcon(0,imat).eq.1) then
362  id=-1
363  else
364  call ident2(plkcon(0,1,imat),t1l,
365  & nplkcon(0,imat),2*npmat_+1,id)
366  endif
367 !
368  if(nplkcon(0,imat).eq.0) then
369  continue
370  elseif((nplkcon(0,imat).eq.1).or.(id.eq.0).or.
371  & (id.eq.nplkcon(0,imat))) then
372  if(id.le.0)then
373  itemp=1
374  else
375  itemp=id
376  endif
377  kin=1
378  call plcopy(plkcon,nplkcon,plconloc,npmat_,
379  & ntmat_,imat,itemp,iel,kin)
380  if((id.eq.0).or.(id.eq.nplkcon(0,imat))) then
381  endif
382  else
383  kin=1
384  call plmix(plkcon,nplkcon,plconloc,npmat_,
385  & ntmat_,imat,id+1,t1l,iel,kin)
386  endif
387  endif
388  endif
389  endif
390 !
391 ! calculating the elastic constants
392 !
393  call ident2(elcon(0,1,imat),t1l,nelcon(2,imat),ncmat_+1,id)
394  if(nelcon(2,imat).eq.0) then
395  continue
396  elseif(nelcon(2,imat).eq.1) then
397  do k=1,nelconst
398  elconloc(k)=elcon(k,1,imat)
399  enddo
400  elseif(id.eq.0) then
401  do k=1,nelconst
402  elconloc(k)=elcon(k,1,imat)
403  enddo
404  elseif(id.eq.nelcon(2,imat)) then
405  do k=1,nelconst
406  elconloc(k)=elcon(k,id,imat)
407  enddo
408  else
409  do k=1,nelconst
410  elconloc(k)=elcon(k,id,imat)+
411  & (elcon(k,id+1,imat)-elcon(k,id,imat))*
412  & (t1l-elcon(0,id,imat))/
413  & (elcon(0,id+1,imat)-elcon(0,id,imat))
414  enddo
415  endif
416 !
417 ! modifying the thermal constants if anisotropic and
418 ! a transformation was defined
419 !
420  if((iorien.ne.0).and.(nalcon(1,imat).gt.1)) then
421 !
422 ! calculating the transformation matrix
423 !
424  call transformatrix(orab(1,iorien),pgauss,skl)
425 !
426 ! transforming the thermal strain
427 !
428  xa(1,1)=eth(1)
429  xa(1,2)=eth(4)
430  xa(1,3)=eth(5)
431  xa(2,1)=eth(4)
432  xa(2,2)=eth(2)
433  xa(2,3)=eth(6)
434  xa(3,1)=eth(5)
435  xa(3,2)=eth(6)
436  xa(3,3)=eth(3)
437 !
438  do jj=1,6
439  eth(jj)=0.d0
440  j1=kal(1,jj)
441  j2=kal(2,jj)
442  do j3=1,3
443  do j4=1,3
444  eth(jj)=eth(jj)+
445  & xa(j3,j4)*skl(j1,j3)*skl(j2,j4)
446  enddo
447  enddo
448  enddo
449  endif
450  endif
451  endif
452 !
453  return
plcopy
subroutine plcopy(plcon, nplcon, plconloc, npmat_, ntmat_, imat, itemp, nelem, kin)
Definition: plcopy.f:21
max
#define max(a, b)
Definition: cascade.c:32
ident2
subroutine ident2(x, px, n, ninc, id)
Definition: ident2.f:27
plmix
subroutine plmix(plcon, nplcon, plconloc, npmat_, ntmat_, imat, j, temp, nelem, kin)
Definition: plmix.f:21
transformatrix
subroutine transformatrix(xab, p, a)
Definition: transformatrix.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)