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

Go to the source code of this file.

Functions/Subroutines

subroutine resultsmech_u1 (co, kon, ipkon, lakon, ne, v, stx, elcon, nelcon, rhcon, nrhcon, alcon, nalcon, alzero, ielmat, ielorien, norien, orab, ntmat_, t0, t1, ithermal, prestr, iprestr, eme, iperturb, fn, iout, qa, vold, nmethod, veold, dtime, time, ttime, plicon, nplicon, plkcon, nplkcon, xstateini, xstiff, xstate, npmat_, matname, mi, ielas, icmd, ncmat_, nstate_, stiini, vini, ener, eei, enerini, istep, iinc, reltime, calcul_fn, calcul_qa, calcul_cauchy, nener, ikin, nal, ne0, thicke, emeini, i, ielprop, prop)
 

Function/Subroutine Documentation

◆ resultsmech_u1()

subroutine resultsmech_u1 ( real*8, dimension(3,*), intent(in)  co,
integer, dimension(*), intent(in)  kon,
integer, dimension(*), intent(in)  ipkon,
character*8, dimension(*), intent(in)  lakon,
integer, intent(in)  ne,
real*8, dimension(0:mi(2),*), intent(in)  v,
real*8, dimension(6,mi(1),*), intent(inout)  stx,
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,
real*8, dimension(*), intent(in)  alzero,
integer, dimension(mi(3),*), intent(in)  ielmat,
integer, dimension(mi(3),*), intent(in)  ielorien,
integer, intent(in)  norien,
real*8, dimension(7,*), intent(in)  orab,
integer, intent(in)  ntmat_,
real*8, dimension(*), intent(in)  t0,
real*8, dimension(*), intent(in)  t1,
integer, dimension(2), intent(in)  ithermal,
real*8, dimension(6,mi(1),*), intent(in)  prestr,
integer, intent(in)  iprestr,
real*8, dimension(6,mi(1),*), intent(inout)  eme,
integer, dimension(*), intent(in)  iperturb,
real*8, dimension(0:mi(2),*), intent(inout)  fn,
integer, intent(in)  iout,
real*8, dimension(*), intent(inout)  qa,
real*8, dimension(0:mi(2),*), intent(in)  vold,
integer, intent(in)  nmethod,
real*8, dimension(0:mi(2),*), intent(in)  veold,
real*8, intent(in)  dtime,
real*8, intent(in)  time,
real*8, intent(in)  ttime,
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,
real*8, dimension(nstate_,mi(1),*), intent(in)  xstateini,
real*8, dimension(27,mi(1),*), intent(inout)  xstiff,
real*8, dimension(nstate_,mi(1),*), intent(in)  xstate,
integer, intent(in)  npmat_,
character*80, dimension(*), intent(in)  matname,
integer, dimension(*), intent(in)  mi,
integer, intent(in)  ielas,
integer, intent(in)  icmd,
integer, intent(in)  ncmat_,
integer, intent(in)  nstate_,
real*8, dimension(6,mi(1),*), intent(in)  stiini,
real*8, dimension(0:mi(2),*), intent(in)  vini,
real*8, dimension(mi(1),*), intent(inout)  ener,
real*8, dimension(6,mi(1),*), intent(inout)  eei,
real*8, dimension(mi(1),*), intent(in)  enerini,
integer, intent(in)  istep,
integer, intent(in)  iinc,
real*8, intent(in)  reltime,
integer, intent(in)  calcul_fn,
integer, intent(in)  calcul_qa,
integer, intent(in)  calcul_cauchy,
integer, intent(in)  nener,
integer, intent(in)  ikin,
integer, intent(inout)  nal,
integer, intent(in)  ne0,
real*8, dimension(mi(3),*), intent(in)  thicke,
real*8, dimension(6,mi(1),*), intent(in)  emeini,
integer, intent(in)  i,
integer, dimension(*), intent(in)  ielprop,
real*8, dimension(*), intent(in)  prop 
)
28 !
29 ! calculates nal,qa,fn,xstiff,ener,eme,eei,stx for user element 1
30 !
31 ! Yunhua Luo, An Efficient 3D Timoshenko Beam Element with
32 ! Consistent Shape Functions, Adv. Theor. Appl. Mech., Vol. 1,
33 ! 2008, no. 3, 95-106
34 !
35  implicit none
36 !
37  character*8 lakon(*)
38  character*80 amat,matname(*)
39 !
40  integer kon(*),konl(26),mi(*),
41  & nelcon(2,*),nrhcon(*),nalcon(2,*),ielmat(mi(3),*),
42  & ielorien(mi(3),*),ntmat_,ipkon(*),ne0,
43  & istep,iinc,ne,mattyp,ithermal(2),iprestr,i,j,k,m1,m2,jj,
44  & i1,kk,nener,indexe,nope,norien,iperturb(*),iout,
45  & nal,icmd,ihyper,nmethod,kode,imat,iorien,ielas,
46  & istiff,ncmat_,nstate_,ikin,ielprop(*),
47  & nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),npmat_,calcul_fn,
48  & calcul_cauchy,calcul_qa,index,node
49 !
50  real*8 co(3,*),v(0:mi(2),*),stiini(6,mi(1),*),
51  & stx(6,mi(1),*),xl(3,26),vl(0:mi(2),26),stre(6),prop(*),
52  & elcon(0:ncmat_,ntmat_,*),rhcon(0:1,ntmat_,*),
53  & alcon(0:6,ntmat_,*),vini(0:mi(2),*),
54  & alzero(*),orab(7,*),elas(21),rho,fn(0:mi(2),*),
55  & q(0:mi(2),26),t0(*),t1(*),prestr(6,mi(1),*),eme(6,mi(1),*),
56  & vold(0:mi(2),*),eloc(9),elconloc(21),eth(6),coords(3),
57  & ener(mi(1),*),emec(6),eei(6,mi(1),*),enerini(mi(1),*),
58  & veold(0:mi(2),*),e,un,um,tt,dl,qa(*),t0l,t1l,dtime,time,ttime,
59  & plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*),
60  & xstiff(27,mi(1),*),xstate(nstate_,mi(1),*),plconloc(802),
61  & xstateini(nstate_,mi(1),*),tm(3,3),a,reltime,
62  & thicke(mi(3),*),emeini(6,mi(1),*),aly,alz,bey,bez,xi(2),
63  & vlp(6,2),xi11,xi12,xi22,xk,offset1,offset2,e1(3),e2(3),e3(3)
64 !
65  intent(in) co,kon,ipkon,lakon,ne,v,
66  & elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
67  & ielmat,ielorien,norien,orab,ntmat_,t0,t1,ithermal,prestr,
68  & iprestr,iperturb,iout,vold,nmethod,
69  & veold,dtime,time,ttime,plicon,nplicon,plkcon,nplkcon,
70  & xstateini,xstate,npmat_,matname,mi,ielas,icmd,
71  & ncmat_,nstate_,stiini,vini,enerini,istep,iinc,
72  & reltime,calcul_fn,calcul_qa,calcul_cauchy,nener,
73  & ikin,ne0,thicke,emeini,i,ielprop,prop
74 !
75  intent(inout) nal,qa,fn,xstiff,ener,eme,eei,stx
76 !
77  nope=2
78 !
79 ! q contains the nodal forces per element; initialization of q
80 !
81  if((iperturb(1).ge.2).or.((iperturb(1).le.0).and.(iout.lt.1)))
82  & then
83  do m1=1,nope
84  do m2=0,mi(2)
85  q(m2,m1)=fn(m2,konl(m1))
86  enddo
87  enddo
88  endif
89 !
90  indexe=ipkon(i)
91 !
92 ! material and orientation
93 !
94  imat=ielmat(1,i)
95  amat=matname(imat)
96  if(norien.gt.0) then
97  iorien=ielorien(1,i)
98  else
99  iorien=0
100  endif
101  if(iorien.gt.0) then
102  write(*,*) '*ERROR in resultsmech_u1: no orientation'
103  write(*,*) ' calculation for this type of element'
104  call exit(201)
105  endif
106 !
107  if(nelcon(1,imat).lt.0) then
108  ihyper=1
109  else
110  ihyper=0
111  endif
112 !
113 ! properties of the cross section
114 !
115  index=ielprop(i)
116  a=prop(index+1)
117  xi11=prop(index+2)
118  xi12=prop(index+3)
119  xi22=prop(index+4)
120  xk=prop(index+5)
121  e2(1)=prop(index+6)
122  e2(2)=prop(index+7)
123  e2(3)=prop(index+8)
124  offset1=prop(index+9)
125  offset2=prop(index+10)
126 !
127  if(dabs(xi12).gt.0.d0) then
128  write(*,*) '*ERROR in resultsmech_u1: no nonzero cross moment'
129  write(*,*) ' of inertia for this type of element'
130  call exit(201)
131  endif
132 !
133  if(dabs(offset1).gt.0.d0) then
134  write(*,*) '*ERROR in resultsmech_u1: no offset in 1-direction'
135  write(*,*) ' for this type of element'
136  call exit(201)
137  endif
138 !
139  if(dabs(offset2).gt.0.d0) then
140  write(*,*) '*ERROR in resultsmech_u1: no offset in 2-direction'
141  write(*,*) ' for this type of element'
142  call exit(201)
143  endif
144 !
145 ! computation of the coordinates and displacements of the local nodes
146 !
147  do k=1,nope
148  konl(k)=kon(indexe+k)
149  do j=1,3
150  xl(j,k)=co(j,konl(k))
151  vl(j,k)=v(j,konl(k))
152  enddo
153  enddo
154 !
155 ! local axes e1-e2-e3 (e2 is given by the user (in the input deck
156 ! this is called e1), e1 is parallel to the beam axis)
157 !
158  do j=1,3
159  e1(j)=xl(j,2)-xl(j,1)
160  enddo
161 !
162  dl=dsqrt(e1(1)*e1(1)+e1(2)*e1(2)+e1(3)*e1(3))
163 !
164  do j=1,3
165  e1(j)=e1(j)/dl
166  enddo
167 !
168 ! e3 = e1 x e2
169 !
170  e3(1)=e1(2)*e2(3)-e1(3)*e2(2)
171  e3(2)=e1(3)*e2(1)-e1(1)*e2(3)
172  e3(3)=e1(1)*e2(2)-e1(2)*e2(1)
173 !
174 ! transformation matrix from the global into the local system
175 !
176  do j=1,3
177  tm(1,j)=e1(j)
178  tm(2,j)=e2(j)
179  tm(3,j)=e3(j)
180  enddo
181 !
182 ! calculating the temperature in the integration
183 ! point
184 !
185  t0l=0.d0
186  t1l=0.d0
187  if(ithermal(1).eq.1) then
188  do i1=1,nope
189  t0l=t0l+t0(konl(i1))/2.d0
190  t1l=t1l+t1(konl(i1))/2.d0
191  enddo
192  elseif(ithermal(1).ge.2) then
193  write(*,*) '*ERROR in resultsmech_u1: no thermal'
194  write(*,*) ' calculation for this type of element'
195  call exit(201)
196  endif
197  tt=t1l-t0l
198 !
199  kode=nelcon(1,imat)
200 !
201  if(kode.eq.2) then
202  mattyp=1
203  else
204  mattyp=0
205  endif
206 !
207 ! material data and local stiffness matrix
208 !
209  istiff=0
210  call materialdata_me(elcon,nelcon,rhcon,nrhcon,alcon,nalcon,
211  & imat,amat,iorien,coords,orab,ntmat_,elas,rho,
212  & i,ithermal,alzero,mattyp,t0l,t1l,
213  & ihyper,istiff,elconloc,eth,kode,plicon,
214  & nplicon,plkcon,nplkcon,npmat_,
215  & plconloc,mi(1),dtime,kk,
216  & xstiff,ncmat_)
217 !
218 ! correcting the thermal strains
219 !
220  do k=2,6
221  eth(k)=0.d0
222  enddo
223 !
224  if(mattyp.eq.1) then
225  e=elconloc(1)
226  un=elconloc(2)
227  um=e/(1.d0+un)
228  um=um/2.d0
229  else
230  write(*,*) '*ERROR in resultsmech_u1: no anisotropic material'
231  write(*,*) ' calculation for this type of element'
232  call exit(201)
233  endif
234 !
235 ! calculating the local displacement and rotation values
236 !
237  do k=1,nope
238  node=kon(indexe+k)
239  do j=1,3
240  vl(j,k)=tm(j,1)*v(1,node)
241  & +tm(j,2)*v(2,node)
242  & +tm(j,3)*v(3,node)
243  vl(j+3,k)=tm(j,1)*v(4,node)
244  & +tm(j,2)*v(5,node)
245  & +tm(j,3)*v(6,node)
246  enddo
247  enddo
248 !
249  xi(1)=0.d0
250  xi(2)=1.d0
251 !
252  aly=12.d0*e*xi11/(xk*um*a*dl*dl)
253  alz=12.d0*e*xi22/(xk*um*a*dl*dl)
254  bey=1.d0/(1.d0-aly)
255  bez=1.d0/(1.d0-alz)
256 !
257  do jj=1,nope
258 !
259 ! calculating the derivative of the local displacement and
260 ! rotation values
261 !
262  vlp(1,jj)=(-vl(1,1)+vl(1,2))/dl
263  vlp(2,jj)=
264  & bey*(6.d0*xi(jj)**2-3.d0*xi(jj)+aly*xi(jj))/dl*vl(2,1)+
265  & bey*(3.d0*xi(jj)**2+(aly-4.d0)*xi(jj)+(1.d0-aly/2.d0))*vl(6,1)
266  & +bey*(-6.d0*xi(jj)**2+6.d0*xi(jj)-aly)/dl*vl(2,2)
267  & +bey*(3.d0*xi(jj)**2-(2.d0+aly)*xi(jj)+aly/2.d0)*vl(6,2)
268  vlp(3,jj)=bez*(6.d0*xi(jj)*xi(jj)-6.d0*xi(jj)+alz)/dl*vl(3,1)+
269  & bez*(3.d0*xi(jj)**2+(alz-4.d0)*xi(jj)+(1.d0-alz/2.d0))*vl(5,1)
270  & +bez*(-6.d0*xi(jj)**2+6.d0*xi(jj)-alz)/dl*vl(3,2)
271  & +bez*(3.d0*xi(jj)**2-(2.d0+alz)*xi(jj)+alz/2.d0)*vl(5,2)
272  vlp(4,jj)=(vl(4,2)-vl(4,1))/dl
273  vlp(5,jj)=6.d0*bez*(2.d0*xi(jj)-1.d0)/(dl*dl)*vl(3,1)
274  & +bez*(6.d0*xi(jj)+(alz-4.d0))/dl*vl(5,1)
275  & +6.d0*bez*(-2.d0*xi(jj)+1.d0)/(dl*dl)*vl(3,2)
276  & +bez*(6.d0*xi(jj)-(alz+2))/dl*vl(5,2)
277  vlp(6,jj)=6.d0*bey*(2.d0*xi(jj)-1.d0)/(dl*dl)*vl(2,1)
278  & +bey*(6.d0*xi(jj)+(aly-4.d0))/dl*vl(6,1)
279  & +6.d0*bey*(-2.d0*xi(jj)+1.d0)/(dl*dl)*vl(2,2)
280  & +bey*(6.d0*xi(jj)-(aly+2.d0))/dl*vl(6,2)
281 !
282 ! calculation of the strains
283 !
284  eloc(1)=vlp(1,jj)
285  eloc(2)=-vlp(6,jj)
286  eloc(3)=vlp(5,jj)
287  eloc(4)=vlp(2,jj)-vl(6,jj)
288  eloc(5)=vlp(3,jj)+vl(5,jj)
289  eloc(6)=vlp(4,jj)
290 !
291 ! determining the mechanical strain
292 !
293  if(ithermal(1).ne.0) then
294  do m1=2,6
295  emec(m1)=eloc(m1)-eth(m1)
296  enddo
297  else
298  do m1=1,6
299  emec(m1)=eloc(m1)
300  enddo
301  endif
302 !
303 ! subtracting initial strains
304 !
305  if(iprestr.ne.0) then
306  write(*,*) '*ERROR in resultsmech_u1:'
307  write(*,*) ' no initial strains allowed for'
308  write(*,*) ' this user element'
309  call exit(201)
310  endif
311 !
312 ! calculating the section forces
313 !
314  stre(1)=e*a*emec(1)
315  stre(2)=e*xi11*emec(2)
316  stre(3)=e*xi22*emec(3)
317  stre(4)=xk*um*a*emec(4)
318  stre(5)=xk*um*a*emec(5)
319  stre(6)=xk*um*(xi11+xi22)*emec(6)
320 !
321 ! updating the internal energy and mechanical strain
322 !
323  if((iout.gt.0).or.(iout.eq.-2).or.(kode.le.-100).or.
324  & ((nmethod.eq.4).and.(iperturb(1).gt.1).and.
325  & (ithermal(1).le.1))) then
326  if(ithermal(1).eq.0) then
327  do m1=1,6
328  eth(m1)=0.d0
329  enddo
330  endif
331  if(nener.eq.1) then
332  ener(jj,i)=enerini(jj,i)+
333  & ((eloc(1)-eth(1)-emeini(1,jj,i))*
334  & (stre(1)+stiini(1,jj,i))+
335  & (eloc(2)-eth(2)-emeini(2,jj,i))*
336  & (stre(2)+stiini(2,jj,i))+
337  & (eloc(3)-eth(3)-emeini(3,jj,i))*
338  & (stre(3)+stiini(3,jj,i)))/2.d0+
339  & (eloc(4)-eth(4)-emeini(4,jj,i))*(stre(4)+stiini(4,jj,i))+
340  & (eloc(5)-eth(5)-emeini(5,jj,i))*(stre(5)+stiini(5,jj,i))+
341  & (eloc(6)-eth(6)-emeini(6,jj,i))*(stre(6)+stiini(6,jj,i))
342  endif
343  eme(1,jj,i)=eloc(1)-eth(1)
344  eme(2,jj,i)=eloc(2)-eth(2)
345  eme(3,jj,i)=eloc(3)-eth(3)
346  eme(4,jj,i)=eloc(4)-eth(4)
347  eme(5,jj,i)=eloc(5)-eth(5)
348  eme(6,jj,i)=eloc(6)-eth(6)
349  endif
350 !
351  if((iout.gt.0).or.(iout.eq.-2).or.(kode.le.-100)) then
352 !
353  eei(1,jj,i)=eloc(1)
354  eei(2,jj,i)=eloc(2)
355  eei(3,jj,i)=eloc(3)
356  eei(4,jj,i)=eloc(4)
357  eei(5,jj,i)=eloc(5)
358  eei(6,jj,i)=eloc(6)
359  endif
360 !
361 ! updating the kinetic energy
362 !
363  if(ikin.eq.1) then
364  write(*,*) '*ERROR in resultsmech_u1:'
365  write(*,*) ' no initial strains allowed for'
366  write(*,*) ' this user element'
367  call exit(201)
368  endif
369 !
370  stx(1,jj,i)=stre(1)
371  stx(2,jj,i)=stre(2)
372  stx(3,jj,i)=stre(3)
373  stx(4,jj,i)=stre(4)
374  stx(5,jj,i)=stre(5)
375  stx(6,jj,i)=stre(6)
376 !
377 ! calculation of the global nodal forces
378 !
379  if(calcul_fn.eq.1)then
380  node=kon(indexe+jj)
381  do j=1,3
382  fn(j,node)=fn(j,node)+tm(1,j)*stre(1)
383  & +tm(2,j)*stre(2)
384  & +tm(3,j)*stre(3)
385  fn(j+3,node)=fn(j+3,node)+tm(1,j)*stre(4)
386  & +tm(2,j)*stre(5)
387  & +tm(3,j)*stre(6)
388  enddo
389  endif
390  enddo
391 !
392 ! q contains the contributions to the nodal force in the nodes
393 ! belonging to the element at stake from other elements (elements
394 ! already treated). These contributions have to be
395 ! subtracted to get the contributions attributable to the element
396 ! at stake only
397 !
398  if(calcul_qa.eq.1) then
399  do m1=1,nope
400  do m2=1,3
401  qa(1)=qa(1)+dabs(fn(m2,konl(m1))-q(m2,m1))
402  enddo
403  enddo
404  nal=nal+3*nope
405  endif
406 !
407  return
materialdata_me
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_)
Definition: materialdata_me.f:23
nal
static ITG * nal
Definition: results.c:31
Hosted by OpenAircraft.com, (Michigan UAV, LLC)