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

Go to the source code of this file.

Functions/Subroutines

subroutine e_c3d_u1 (co, kon, lakonl, p1, p2, omx, bodyfx, nbody, s, sm, ff, nelem, nmethod, elcon, nelcon, rhcon, nrhcon, alcon, nalcon, alzero, ielmat, ielorien, norien, orab, ntmat_, t0, t1, ithermal, vold, iperturb, nelemload, sideload, xload, nload, idist, sti, stx, iexpl, plicon, nplicon, plkcon, nplkcon, xstiff, npmat_, dtime, matname, mi, ncmat_, mass, stiffness, buckling, rhsi, intscheme, ttime, time, istep, iinc, coriolis, xloadold, reltime, ipompc, nodempc, coefmpc, nmpc, ikmpc, ilmpc, veold, ne0, ipkon, thicke, integerglob, doubleglob, tieset, istartset, iendset, ialset, ntie, nasym, ielprop, prop)
 

Function/Subroutine Documentation

◆ e_c3d_u1()

subroutine e_c3d_u1 ( real*8, dimension(3,*), intent(in)  co,
integer, dimension(*), intent(in)  kon,
character*8, intent(in)  lakonl,
real*8, dimension(3), intent(in)  p1,
real*8, dimension(3), intent(in)  p2,
real*8, intent(in)  omx,
real*8, dimension(3), intent(in)  bodyfx,
integer, intent(in)  nbody,
real*8, dimension(60,60), intent(inout)  s,
real*8, dimension(60,60), intent(inout)  sm,
real*8, dimension(60), intent(inout)  ff,
integer, intent(in)  nelem,
integer, intent(inout)  nmethod,
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, intent(in)  ithermal,
real*8, dimension(0:mi(2),*), intent(in)  vold,
integer, dimension(*), intent(in)  iperturb,
integer, dimension(2,*), intent(in)  nelemload,
character*20, dimension(*), intent(in)  sideload,
real*8, dimension(2,*), intent(inout)  xload,
integer, intent(in)  nload,
integer, intent(in)  idist,
real*8, dimension(6,mi(1),*), intent(in)  sti,
real*8, dimension(6,mi(1),*), intent(in)  stx,
integer, intent(in)  iexpl,
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(27,mi(1),*), intent(in)  xstiff,
integer, intent(in)  npmat_,
real*8, intent(in)  dtime,
character*80, dimension(*), intent(in)  matname,
integer, dimension(*), intent(in)  mi,
integer, intent(in)  ncmat_,
integer, intent(in)  mass,
integer, intent(in)  stiffness,
integer, intent(in)  buckling,
integer, intent(in)  rhsi,
integer, intent(in)  intscheme,
real*8, intent(in)  ttime,
real*8, intent(in)  time,
integer, intent(in)  istep,
integer, intent(in)  iinc,
integer, intent(in)  coriolis,
real*8, dimension(2,*), intent(in)  xloadold,
real*8, intent(in)  reltime,
integer, dimension(*), intent(in)  ipompc,
integer, dimension(3,*), intent(in)  nodempc,
real*8, dimension(*), intent(in)  coefmpc,
integer, intent(in)  nmpc,
integer, dimension(*), intent(in)  ikmpc,
integer, dimension(*), intent(in)  ilmpc,
real*8, dimension(0:mi(2),*), intent(in)  veold,
integer, intent(in)  ne0,
integer, dimension(*), intent(in)  ipkon,
real*8, dimension(mi(3),*), intent(in)  thicke,
integer, dimension(*), intent(in)  integerglob,
real*8, dimension(*), intent(in)  doubleglob,
character*81, dimension(3,*), intent(in)  tieset,
integer, dimension(*), intent(in)  istartset,
integer, dimension(*), intent(in)  iendset,
integer, dimension(*), intent(in)  ialset,
integer, intent(in)  ntie,
integer, intent(in)  nasym,
integer, dimension(*), intent(in)  ielprop,
real*8, dimension(*), intent(in)  prop 
)
31 !
32 ! computation of the element matrix and rhs for the user element
33 ! of type u1
34 !
35 ! This is a beam type element. Reference:
36 ! Yunhua Luo, An Efficient 3D Timoshenko Beam Element with
37 ! Consistent Shape Functions, Adv. Theor. Appl. Mech., Vol. 1,
38 ! 2008, no. 3, 95-106
39 !
40 ! special case for which the beam axis goes through the
41 ! center of gravity of the cross section and the 1-direction
42 ! corresponds with a principal axis
43 !
44 !
45 ! INPUT:
46 !
47 ! co(1..3,i) coordinates of node i
48 ! kon(*) contains the topology of all elements. The
49 ! topology of element i starts at kon(ipkon(i)+1)
50 ! and continues until all nodes are covered. The
51 ! number of nodes depends on the element label
52 ! lakonl label of current element nelem (character*8)
53 ! p1(1..3) coordinates of a point on the rotation axis
54 ! (if applicable)
55 ! p2(1..3) unit vector on rotation axis (if applicable)
56 ! omx rotational speed square
57 ! bodyfx(1..3) acceleration vector
58 ! nbody number of body loads
59 ! nelem element number
60 ! nmethod procedure:
61 ! 1: static analysis
62 ! 2: frequency analysis
63 ! 3: buckling analysis
64 ! 4: (linear or nonlinear) dynamic analysis
65 ! 5: steady state dynamics analysis
66 ! 6: Coriolis frequency calculation
67 ! 7: flutter frequency calculation
68 ! 8: magnetostatics
69 ! 9: magnetodynamics
70 ! 10: electromagnetic eigenvalue problems
71 ! 11: superelement creation or Green function
72 ! calculation
73 ! 12: sensitivity analysis
74 ! elcon elastic constants (cf. List of variables and
75 ! their meaning in the User's Manual)
76 ! nelcon integers describing the elastic constant fields
77 ! (cf. User's Manual)
78 ! rhcon density constants (cf. User's Manual)
79 ! nrhcon integers describing the density constant fields
80 ! (cf. User's Manual)
81 ! alcon thermal expansion constants (cf. User's Manual)
82 ! nalcon integers describing the thermal expansion
83 ! constants (cf. User's Manual)
84 ! alzero thermal expansion reference values (cf. User's Manual)
85 ! ielmat(i) material for element i
86 ! ielorien(i) orientation for element i
87 ! norien number of orientations
88 ! orab(7,*) description of all local coordinate systems.
89 ! (cf. List of variables and their meaning in the
90 ! User's manual)
91 ! ntmat_ maximum number of material temperature points
92 ! t0(i) temperature in node i at start of calculation
93 ! t1(i) temperature in node i at the end of the current
94 ! increment
95 ! ithermal cf. ithermal(1) in the List of variables and
96 ! their meaning in the User's Manual
97 ! vold(0..mi(2),i) value of the variables in node i at the end
98 ! of the previous iteration
99 ! iperturb(*) describes the kind of nonlinearity of the
100 ! calculation, cf. User's Manual
101 ! nelemload field describing facial distributed load (cf.
102 ! User's Manual)
103 ! sideload field describing facial distributed load (cf.
104 ! User's Manual)
105 ! xload facial distributed load values (cf.
106 ! User's Manual)
107 ! nload number of facial distributed loads
108 ! idist 0: no distributed forces in the model
109 ! 1: distributed forces are present in the model
110 ! (body forces or thermal loads or residual
111 ! stresses or distributed facial loads)
112 ! sti(i,j,k) stress component i in integration point j
113 ! of element k at the end of the previous step
114 ! stx(i,j,k) stress component i in integration point j
115 ! of element k at the end of a static step
116 ! describing a reference buckling load
117 ! iexpl 0: structure: implicit dynamics,
118 ! fluid: incompressible
119 ! iexpl 1: structure: implicit dynamics,
120 ! fluid: compressible
121 ! iexpl 2: structure: explicit dynamics,
122 ! fluid: incompressible
123 ! iexpl 3: structure: explicit dynamics,
124 ! fluid: compressible
125 ! plicon,nplicon fields describing isotropic hardening of
126 ! a plastic material or spring constants of
127 ! a nonlinear spring (cf. User's Manual)
128 ! plkcon,nplkcon fields describing kinematic hardening of
129 ! a plastic material or gap conductance
130 ! constants (cf. User's Manual)
131 ! xstiff(i,j,k) stiffness (i=1...21) and conductivity
132 ! (i=22..27) constants in integration point j
133 ! of element k
134 ! npmat_ maximum number of plastic constants
135 ! dtime length of present time increment
136 ! matname(i) name of material i
137 ! mi(1) max # of integration points per element (max
138 ! over all elements)
139 ! mi(2) max degree of freedom per node (max over all
140 ! nodes) in fields like v(0:mi(2))...
141 ! mi(3) max number of layers in the structure
142 ! ncmat_ max number of elastic constants
143 ! mass(1) if 1: mass matrix needed, else not
144 ! mass(2) if 1: heat capacity matrix needed, else not
145 ! stiffness if 1: stiffness matrix needed, else not
146 ! buckling if 1: linear buckling calculation, else not
147 ! rhsi if 1: right hand side needed, else not
148 ! intscheme 1: integration point scheme for the calculation
149 ! of the initial acceleration in a dynamic
150 ! step
151 ! 0: any other case
152 ! ttime total time at the start of the present increment
153 ! time step time at the end of the present increment
154 ! istep current step number
155 ! iinc current increment number within the actual step
156 ! coriolis if 1: Coriolis forces are taken into account,
157 ! else not
158 ! xloadold load at the end of the previous increment
159 ! (cf. User's Manual)
160 ! reltime relative step time (between 0 and 1)
161 ! ipompc(i) pointer into field nodempc and coefmpc for
162 ! MPC i
163 ! nodempc field containing the nodes and dof's of all
164 ! MPC's (cf. User's Manual)
165 ! coefmpc field containing the coefficients of all
166 ! MPC's (cf. User's Manual)
167 ! nmpc number of MPC's in the model
168 ! ikmpc(i) field containing an integer code number
169 ! describing the dependent dof of some MPC:
170 ! ikmpc(i)=8*(node-1)+dir, where the dependent
171 ! dof is applied in direction dir of node "node"
172 ! ikmpc is sorted in ascending order
173 ! ilmpc(i) MPC number corresponding to ikmpc(i) (cf.
174 ! User's manual)
175 ! veold(j,i) time rate of variable j in node i at the end
176 ! of the previous iteration
177 ! ne0 element number before the first contact element;
178 ! contact elements are appended after ne0
179 ! ipkon(i) points to the location in field kon preceding
180 ! the topology of element i
181 ! thicke(j,i) thickness of layer j in node i
182 ! integerglob integer field needed for determining the
183 ! interface loading of a submodel
184 ! doubleglob real field needed for determining the
185 ! interface loading of a submodel
186 ! tieset(1..3,*) tie character information for tie i (cf.
187 ! User's Manual)
188 ! istartset integer describing set information (cf.
189 ! User's Manual)
190 ! iendset integer describing set information (cf.
191 ! User's Manual)
192 ! ialset integer describing set information (cf.
193 ! User's Manual)
194 ! ntie number of ties in the model
195 ! nasym 1: asymmetric contributions, else purely
196 ! symmetric
197 ! ielprop(i) points to the location in field prop preceding
198 ! the properties of element i
199 ! prop(*) contains the properties and some beam
200 ! elements (cf. User's Manual)
201 !
202 !
203 ! OUTPUT:
204 !
205 ! s element stiffness matrix
206 ! sm element mass matrix
207 ! ff element load vector
208 ! xload facial distributed load values (cf.
209 ! User's Manual)
210 ! nmethod may be changed by the user, e.g. a
211 ! change to 0 for negative Jacobians leads
212 ! to a program exit in the calling program
213 ! for the other value: cf. above
214 !
215  implicit none
216 !
217  integer mass,stiffness,buckling,rhsi,coriolis
218 !
219  character*8 lakonl
220  character*20 sideload(*)
221  character*80 matname(*),amat
222  character*81 tieset(3,*)
223 !
224  integer konl(26),nelemload(2,*),nbody,nelem,mi(*),kon(*),
225  & ielprop(*),index,mattyp,ithermal,iperturb(*),nload,idist,
226  & i,j,k,i1,nmethod,kk,nelcon(2,*),nrhcon(*),nalcon(2,*),
227  & ielmat(mi(3),*),ielorien(mi(3),*),ipkon(*),indexe,
228  & ntmat_,nope,norien,ihyper,iexpl,kode,imat,iorien,istiff,
229  & ncmat_,intscheme,istep,iinc,ipompc(*),nodempc(3,*),
230  & nmpc,ikmpc(*),ilmpc(*),ne0,ndof,istartset(*),iendset(*),
231  & ialset(*),ntie,integerglob(*),nasym,nplicon(0:ntmat_,*),
232  & nplkcon(0:ntmat_,*),npmat_
233 !
234  real*8 co(3,*),xl(3,26),veold(0:mi(2),*),rho,s(60,60),bodyfx(3),
235  & ff(60),elconloc(21),coords(3),p1(3),elcon(0:ncmat_,ntmat_,*),
236  & p2(3),eth(6),rhcon(0:1,ntmat_,*),reltime,prop(*),tm(3,3),
237  & alcon(0:6,ntmat_,*),alzero(*),orab(7,*),t0(*),t1(*),
238  & xloadold(2,*),vold(0:mi(2),*),xload(2,*),omx,e,un,um,tt,
239  & sm(60,60),sti(6,mi(1),*),stx(6,mi(1),*),t0l,t1l,coefmpc(*),
240  & elas(21),thicke(mi(3),*),doubleglob(*),dl,e2(3),e3(3),
241  & plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*),
242  & xstiff(27,mi(1),*),plconloc(802),dtime,ttime,time,tmg(12,12),
243  & a,xi11,xi12,xi22,xk,e1(3),offset1,offset2,y1,y2,y3,z1,z2,z3,
244  & sg(12,12)
245 !
246  intent(in) co,kon,lakonl,p1,p2,omx,bodyfx,nbody,
247  & nelem,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
248  & ielmat,ielorien,norien,orab,ntmat_,
249  & t0,t1,ithermal,vold,iperturb,nelemload,
250  & sideload,nload,idist,sti,stx,iexpl,plicon,
251  & nplicon,plkcon,nplkcon,xstiff,npmat_,dtime,
252  & matname,mi,ncmat_,mass,stiffness,buckling,rhsi,intscheme,
253  & ttime,time,istep,iinc,coriolis,xloadold,reltime,
254  & ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,veold,
255  & ne0,ipkon,thicke,
256  & integerglob,doubleglob,tieset,istartset,iendset,ialset,ntie,
257  & nasym,ielprop,prop
258 !
259  intent(inout) s,sm,ff,xload,nmethod
260 !
261  indexe=ipkon(nelem)
262 !
263 ! material and orientation
264 !
265  imat=ielmat(1,nelem)
266  amat=matname(imat)
267  if(norien.gt.0) then
268  iorien=ielorien(1,nelem)
269  else
270  iorien=0
271  endif
272 !
273  if(nelcon(1,imat).lt.0) then
274  ihyper=1
275  else
276  ihyper=0
277  endif
278 !
279 ! properties of the cross section
280 !
281  index=ielprop(nelem)
282  a=prop(index+1)
283  xi11=prop(index+2)
284  xi12=prop(index+3)
285  xi22=prop(index+4)
286  xk=prop(index+5)
287  e2(1)=prop(index+6)
288  e2(2)=prop(index+7)
289  e2(3)=prop(index+8)
290  offset1=prop(index+9)
291  offset2=prop(index+10)
292 !
293  nope=2
294  ndof=6
295 !
296  if(dabs(xi12).gt.0.d0) then
297  write(*,*) '*ERROR in e_c3d_u1: no nonzero cross moment'
298  write(*,*) ' of inertia for this type of element'
299  call exit(201)
300  endif
301 !
302  if(dabs(offset1).gt.0.d0) then
303  write(*,*) '*ERROR in e_c3d_u1: no offset in 1-direction'
304  write(*,*) ' for this type of element'
305  call exit(201)
306  endif
307 !
308  if(dabs(offset2).gt.0.d0) then
309  write(*,*) '*ERROR in e_c3d_u1: no offset in 2-direction'
310  write(*,*) ' for this type of element'
311  call exit(201)
312  endif
313 !
314 ! computation of the coordinates of the local nodes
315 !
316  do i=1,nope
317  konl(i)=kon(indexe+i)
318  do j=1,3
319  xl(j,i)=co(j,konl(i))
320  enddo
321  enddo
322 !
323 ! local axes e1-e2-e3 (e2 is given by the user (in the input deck
324 ! this is called e1), e1 is parallel to the beam axis)
325 !
326  do j=1,3
327  e1(j)=xl(j,2)-xl(j,1)
328  enddo
329 !
330  dl=dsqrt(e1(1)*e1(1)+e1(2)*e1(2)+e1(3)*e1(3))
331 !
332  do j=1,3
333  e1(j)=e1(j)/dl
334  enddo
335 !
336 ! e3 = e1 x e2
337 !
338  e3(1)=e1(2)*e2(3)-e1(3)*e2(2)
339  e3(2)=e1(3)*e2(1)-e1(1)*e2(3)
340  e3(3)=e1(1)*e2(2)-e1(2)*e2(1)
341 !
342 ! transformation matrix from the global into the local system
343 !
344  do j=1,3
345  tm(1,j)=e1(j)
346  tm(2,j)=e2(j)
347  tm(3,j)=e3(j)
348  enddo
349 !
350 ! initialisation for distributed forces
351 !
352  if(rhsi.eq.1) then
353  if(idist.ne.0) then
354  do i=1,ndof*nope
355  ff(i)=0.d0
356  enddo
357  endif
358  endif
359 !
360 ! displacements for 2nd order static and modal theory
361 !
362  if(((iperturb(1).eq.1).or.(iperturb(2).eq.1)).and.
363  & (stiffness.eq.1).and.(buckling.eq.0)) then
364  write(*,*) '*ERROR in e_c3d_u1: no second order'
365  write(*,*) ' calculation for this type of element'
366  call exit(201)
367  endif
368 !
369 ! initialisation of sm
370 !
371  if((mass.eq.1).or.(buckling.eq.1).or.(coriolis.eq.1)) then
372  write(*,*) '*ERROR in e_c3d_u1: no dynamic or buckling'
373  write(*,*) ' calculation for this type of element'
374  call exit(201)
375  endif
376 !
377 ! initialisation of s
378 !
379  do i=1,ndof*nope
380  do j=1,ndof*nope
381  s(i,j)=0.d0
382  enddo
383  enddo
384 !
385 ! computation of the matrix
386 !
387 ! calculating the temperature in the integration
388 ! point
389 !
390  t0l=0.d0
391  t1l=0.d0
392  if(ithermal.eq.1) then
393  do i1=1,nope
394  t0l=t0l+t0(konl(i1))/2.d0
395  t1l=t1l+t1(konl(i1))/2.d0
396  enddo
397  elseif(ithermal.ge.2) then
398  write(*,*) '*ERROR in e_c3d_u1: no thermal'
399  write(*,*) ' calculation for this type of element'
400  call exit(201)
401  endif
402  tt=t1l-t0l
403 !
404 ! calculating the coordinates of the integration point
405 ! for material orientation purposes (for cylindrical
406 ! coordinate systems)
407 !
408  if(iorien.gt.0) then
409  write(*,*) '*ERROR in e_c3d_u1: no orientation'
410  write(*,*) ' calculation for this type of element'
411  call exit(201)
412  endif
413 !
414 ! for deformation plasticity: calculating the Jacobian
415 ! and the inverse of the deformation gradient
416 ! needed to convert the stiffness matrix in the spatial
417 ! frame of reference to the material frame
418 !
419  kode=nelcon(1,imat)
420 !
421  if(kode.eq.2) then
422  mattyp=1
423  else
424  mattyp=0
425  endif
426 !
427 ! material data and local stiffness matrix
428 !
429  istiff=0
430  call materialdata_me(elcon,nelcon,rhcon,nrhcon,alcon,nalcon,
431  & imat,amat,iorien,coords,orab,ntmat_,elas,rho,
432  & nelem,ithermal,alzero,mattyp,t0l,t1l,
433  & ihyper,istiff,elconloc,eth,kode,plicon,
434  & nplicon,plkcon,nplkcon,npmat_,
435  & plconloc,mi(1),dtime,kk,
436  & xstiff,ncmat_)
437 !
438  if(mattyp.eq.1) then
439  e=elconloc(1)
440  un=elconloc(2)
441  um=e/(1.d0+un)
442  um=um/2.d0
443  elseif(mattyp.eq.2) then
444  write(*,*) '*ERROR in e_c3d_u1: no orthotropic material'
445  write(*,*) ' calculation for this type of element'
446  call exit(201)
447  else
448  write(*,*) '*ERROR in e_c3d_u1: no anisotropic material'
449  write(*,*) ' calculation for this type of element'
450  call exit(201)
451  endif
452 !
453 ! initialization for the body forces
454 !
455  if(rhsi.eq.1) then
456  if(nbody.ne.0) then
457  write(*,*) '*ERROR in e_c3d_u1: no body forces'
458  write(*,*) ' for this type of element'
459  call exit(201)
460  endif
461  endif
462 !
463  if(buckling.eq.1) then
464  write(*,*) '*ERROR in e_c3d_u1: no buckling '
465  write(*,*) ' calculation for this type of element'
466  call exit(201)
467  endif
468 !
469 ! determination of the stiffness, and/or mass and/or
470 ! buckling matrix
471 !
472  if((stiffness.eq.1).or.(mass.eq.1).or.(buckling.eq.1).or.
473  & (coriolis.eq.1)) then
474 !
475  if(((iperturb(1).ne.1).and.(iperturb(2).ne.1)).or.
476  & (buckling.eq.1)) then
477 !
478 ! stiffness matrix
479 !
480  y1=xk*um*a*e*xi11*(12.d0*e*xi11+xk*um*a*dl*dl)
481  y2=(12.d0*e*xi11-xk*um*a*dl*dl)**2
482  y3=4.d0*e*xi11*((xk*um*a)**2*dl**4+
483  & 3.d0*xk*um*a*dl*dl*e*xi11+
484  & 36.d0*(e*xi11)**2)
485  z1=xk*um*a*e*xi22*(12.d0*e*xi22+xk*um*a*dl*dl)
486  z2=(12.d0*e*xi22-xk*um*a*dl*dl)**2
487  z3=4.d0*e*xi22*((xk*um*a)**2*dl**4+
488  & 3.d0*xk*um*a*dl*dl*e*xi22+
489  & 36.d0*(e*xi22)**2)
490 !
491 ! stiffness matrix S' in local coordinates
492 !
493  s(1,1)=e*a/dl
494  s(1,7)=-s(1,1)
495  s(2,2)=12.d0*y1/(dl*y2)
496  s(2,6)=6.d0*y1/y2
497  s(2,8)=-s(2,2)
498  s(2,12)=s(2,6)
499  s(3,3)=12.d0*z1/(dl*z2)
500  s(3,5)=-6.d0*z1/z2
501  s(3,9)=s(3,3)
502  s(3,11)=s(3,5)
503  s(4,4)=um*(xi11+xi22)/dl
504  s(4,10)=-s(4,4)
505  s(5,5)=z3/(dl*z2)
506  s(5,9)=6.d0*z1/z2
507  s(5,11)=-2.d0*e*xi22*(72.d0*(e*xi22)**2-
508  & (xk*um*a)**2*dl**4-
509  & 30.d0*xk*um*a*dl*dl*e*xi22)/(dl*z2)
510  s(6,6)=y3/(dl*y2)
511  s(6,8)=-6.d0*y1/y2
512  s(6,12)=-2.d0*e*xi11*(-(xk*um*a)**2*dl**4-
513  & 30.d0*xk*um*a*dl*dl*e*xi11+
514  & 72.d0*(e*xi11)**2)/(dl*y2)
515  s(1,7)=-s(1,1)
516  s(7,7)=s(1,1)
517  s(8,8)=12.d0*y1/(dl*y2)
518  s(8,12)=-6.d0*y1/y2
519  s(9,9)=12.d0*z1/(dl*z2)
520  s(9,11)=6.d0*z1/z2
521  s(10,10)=s(4,4)
522  s(11,11)=z3/(dl*z2)
523  s(12,12)=y3/(dl*y2)
524 !
525 ! completing the symmetric part
526 !
527  do i=1,12
528  do j=1,i
529  s(i,j)=s(j,i)
530  enddo
531  enddo
532 !
533 ! 12 x 12 transformation matrix
534 !
535  do i=1,12
536  do j=1,12
537  tmg(i,j)=0.d0
538  enddo
539  enddo
540  do i=1,3
541  do j=1,3
542  tmg(i,j)=tm(i,j)
543  tmg(i+3,j+3)=tm(i,j)
544  tmg(i+6,j+6)=tm(i,j)
545  tmg(i+9,j+9)=tm(i,j)
546  enddo
547  enddo
548 !
549 ! stiffness matrix in global coordinates: S = T^T*S'*T
550 !
551  do i=1,12
552  do j=1,12
553  sg(i,j)=0.d0
554  do k=1,12
555  sg(i,j)=sg(i,j)+s(i,k)*tmg(k,j)
556  enddo
557  enddo
558  enddo
559 !
560 ! only upper triangular matrix
561 !
562  do i=1,12
563  do j=i,12
564  s(i,j)=0.d0
565  do k=1,12
566  s(i,j)=s(i,j)+tmg(k,i)*sg(k,j)
567  enddo
568  enddo
569  enddo
570 !
571  endif
572 !
573  endif
574 !
575  return
z1
static double * z1
Definition: filtermain.c:48
idist
static ITG * idist
Definition: radflowload.c:39
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
Hosted by OpenAircraft.com, (Michigan UAV, LLC)