1!$Id:$
2      subroutine fld2d1(d,ul,xl,ix,s,r,ndf,ndm,nst,isw)
3
4!      * * F E A P * * A Finite Element Analysis Program
5
6!....  Copyright (c) 1984-2017: Regents of the University of California
7!                               All rights reserved
8
9!-----[--.----+----.----+----.-----------------------------------------]
10!      Purpose:  2-D Finite Deformation Elasticity Routine
11!                Remark: This is a standard displacement model
12
13!      Inputs:
14!         d(*)      - Material set parameters
15!         ul(ndf,*) - Nodal solution parameters for element
16!         xl(ndm,*) - Nodal coordinates for element
17!         ix(*)     - Element nodal connection list
18!         ndf       - Number dof/node
19!         ndm       - Spatial dimension of mesh
20!         nst       - Dimension of element arrays
21!         isw       - Switch to control action
22
23!      Outputs:
24!         s(nst,*)  - Element matrix
25!         p(nst)    - Element vector
26!-----[--.----+----.----+----.-----------------------------------------]
27
28      implicit  none
29
30      include  'debugs.h'
31      include  'bdata.h'
32      include  'cdata.h'
33      include  'eldata.h'
34      include  'elengy.h'
35      include  'elplot.h'
36      include  'eltran.h'
37      include  'hdata.h'
38      include  'iofile.h'
39      include  'pmod2d.h'
40      include  'prstrs.h'
41      include  'ptdat6.h'
42      include  'qudshp.h'
43      include  'comblk.h'
44
45      integer   ndf,ndm,nst,isw
46      integer   i,i1, j,jj,j1, l, nn,nhv, istrt
47
48      real*8    bdb,bd3,dl,dc,di,dvol0,dvol,dmas0, jac0
49      real*8    cfac,lfac,xx1,xx2,xx3,yy1, tempi, ta
50
51      integer   ix(*)
52      real*8    d(*),ul(ndf,nen,*),xl(ndm,*),s(nst,*)
53      real*8    r(ndf,*),r1(3,64),al(2),ac(2),vl(2)
54      real*8    bbd(4,2),bb(6),shpr(64)
55      real*8    df(3,3,64),f(9,2,64),detf(2,64)
56      real*8    ds(6,6,5),dd(6,6),sigv(9),sigl(9,64)
57
58      save
59
60      data      lint / 0 /
61      data      f    /1152*0.0d0 /
62
63!     TEMPORARY SET OF TEMPERATURE VALUE
64
65      data      ta   /    0.0d0/
66
67!     No action for isw = 1
68
69      if(isw.eq.1) return
70
71!     Set quadrature and go to process isw
72
73      call quadr2d(d,.true.)
74
75      if(debug) call mprint(sg2,3,lint,3,'SG2')
76
77!     COMPUTE TANGENT STIFFNESS AND RESIDUAL FORCE VECTOR
78
79!     Compute shape functions and derivatives in reference configuration
80
81      do l = 1,lint
82        call interp2d(l, xl,ix,ndm,nel, .false.)
83
84        if(debug) call mprint(shp2(1,1,l),3,nel,3,'SHP')
85      end do
86      if(debug) call mprint(jac,1,lint,1,'JAC')
87
88!     Compute deformation gradient and determinant; transform shape
89!     functions to current configuration.
90
91      do l = 1,lint
92        call kine2d(shp2(1,1,l),xl,ul,f(1,1,l),df(1,1,l),
93     &              detf(1,l),ndm,ndf,nel,nen)
94      end do
95
96!     Integration order set to static
97
98      if(d(7).lt.0.0d0) then
99        cfac = 0.0d0
100        lfac = 0.0d0
101      else
102        cfac = d(7)
103        lfac = 1.0d0 - cfac
104      endif
105
106      nhv   = nint(d(15))
107      istrt = nint(d(84))
108
109      if(mod(isw,4).eq.0) go to 4
110
111!     LOOP OVER GAUSS POINTS
112
113      nn = 0
114      do l = 1,lint
115
116!       Set reference coordinates
117
118        xx1 = 0.0d0
119        xx2 = 0.0d0
120        do i = 1,nel
121          xx1 = xx1 + shp2(3,i,l)*xl(1,i)
122          xx2 = xx2 + shp2(3,i,l)*xl(2,i)
123        end do
124
125!       Check for axisymmetry
126
127        if(stype.eq.3) then
128          jac0  = jac(l)
129          dvol0 = jac0*xx1
130          do i = 1,nel
131            shpr(i) = shp2(3,i,l)/xx1
132          end do ! i
133        else
134          jac0  = 0.0d0
135          dvol0 = jac(l)
136          do i = 1,nel
137            shpr(i) = 0.0d0
138          end do ! i
139        end if
140        dvol  = dvol0*detf(1,l)
141        dmas0 = dvol0*d(4)
142
143!       Compute Cauchy stresses and spatial tangent tensor
144
145        call modlfd(l,d,f(1,1,l),df(1,1,l),detf(1,l),ta,
146     &             hr(nn+nh1),hr(nn+nh2),nhv,istrt,ds,sigv,bb,
147     &             .false.,isw)
148
149        if(isw.eq.13) then
150
151          epl(8) = epl(8) + estore*dvol0
152
153!         Compute velocity at point
154
155          vl(1) = shp2(3,1,l)*ul(1,1,4) + shp2(3,2,l)*ul(1,2,4)
156     &          + shp2(3,3,l)*ul(1,3,4) + shp2(3,4,l)*ul(1,4,4)
157
158          vl(2) = shp2(3,1,l)*ul(2,1,4) + shp2(3,2,l)*ul(2,2,4)
159     &          + shp2(3,3,l)*ul(2,3,4) + shp2(3,4,l)*ul(2,4,4)
160
161          tempi = 0.0d0
162          do i = 1,nel
163            tempi = tempi
164     &            + (ul(1,i,4)**2 + ul(2,i,4)**2)*shp2(3,i,l)
165          end do ! i
166
167!         Accumulate kinetic energy
168
169          epl(7) = epl(7) + 0.5d0*(lfac*tempi
170     &                    + cfac*(vl(1)**2 + vl(2)**2))*dmas0
171
172        elseif(isw.ne.14) then
173
174!       Store stress values for tplot
175
176        j1 = 6*(l-1)
177        do j = 1,4
178          tt(j+j1) = sigv(j)
179        end do
180
181!       Multiply tangent moduli and stresses by volume element.
182
183        do i = 1,4
184          sigv(i) = sigv(i)*dvol
185          do j = 1,4
186            dd(i,j) = ds(i,j,1)*dvol*ctan(1)
187          end do
188        end do
189
190!       Compute accelerations
191
192        al(1) = cfac*(shp2(3,1,l)*ul(1,1,5) + shp2(3,2,l)*ul(1,2,5)
193     &              + shp2(3,3,l)*ul(1,3,5) + shp2(3,4,l)*ul(1,4,5))
194
195        al(2) = cfac*(shp2(3,1,l)*ul(2,1,5) + shp2(3,2,l)*ul(2,2,5)
196     &              + shp2(3,3,l)*ul(2,3,5) + shp2(3,4,l)*ul(2,4,5))
197
198!       COMPUTE STRESS DIVERGENCE AND INERTIA TERMS
199
200        xx1 = xx1*d(65)**2
201        xx2 = xx2*d(65)**2
202
203        do i = 1,nel
204
205!         Compute inertial and body load effects
206
207          ac(1)   = (al(1) + lfac*ul(1,i,5))*dmas0
208     &            -  d(11)*dvol0 - xx1*dmas0
209          ac(2)   = (al(2) + lfac*ul(2,i,5))*dmas0
210     &            -  d(12)*dvol0 - xx2*dmas0
211
212!         Stress divergence term (used in geometric stiffness)
213
214          r1(1,i) = shp2(1,i,l)*sigv(1) + shp2(2,i,l)*sigv(4)
215          r1(2,i) = shp2(1,i,l)*sigv(4) + shp2(2,i,l)*sigv(2)
216
217!         Element residual
218
219          r(1,i) = r(1,i) - r1(1,i) - ac(1)*shp2(3,i,l)
220     &                    - shpr(i)*sigv(3)
221          r(2,i) = r(2,i) - r1(2,i) - ac(2)*shp2(3,i,l)
222
223        end do
224
225!       COMPUTE K (s(nst,nst) = K)
226
227        if(isw.eq.3) then
228
229!         PART 1. - Geometric and inertial part.
230
231          dc  = cfac*ctan(3)*dmas0
232          dl  = lfac*ctan(3)*dmas0
233          i1  = 0
234          do i = 1,nel
235
236            do jj = 1,2
237              s(i1+jj,i1+jj) = s(i1+jj,i1+jj) + shp2(3,i,l)*dl
238            end do
239
240            di  = dc*shp2(3,i,l)
241
242!           Include geometric stiffness
243
244            if(gflag) then
245              bd3 = shpr(i)*sigv(3)*ctan(1)
246              j1  = 0
247              do j = 1,i
248                bdb          = (r1(1,i)*shp2(1,j,l)
249     &                       +  r1(2,i)*shp2(2,j,l))*ctan(1)
250     &                       +       di*shp2(3,j,l)
251                s(i1+1,j1+1) = s(i1+1,j1+1) + bdb + bd3*shpr(j)
252                s(i1+2,j1+2) = s(i1+2,j1+2) + bdb
253                j1 = j1 + ndf
254              end do
255
256!           Include inertia only
257
258            else
259              j1  = 0
260              do j = 1,i
261                bdb          = di*shp2(3,j,l)
262                s(i1+1,j1+1) = s(i1+1,j1+1) + bdb
263                s(i1+2,j1+2) = s(i1+2,j1+2) + bdb
264                j1 = j1 + ndf
265              end do
266            endif
267
268            i1 = i1 + ndf
269          end do
270
271!         PART 2. - Tangent modulus part (based upon dd-array)
272
273          i1 = 0
274          do i  = 1,nel
275
276!           Compute bmat-t * dd * dvol
277
278            do jj = 1,4
279
280              bbd(jj,1) = shp2(1,i,l)*dd(1,jj)
281     &                  + shpr(  i  )*dd(3,jj)
282     &                  + shp2(2,i,l)*dd(4,jj)
283
284              bbd(jj,2) = shp2(1,i,l)*dd(4,jj)
285     &                  + shp2(2,i,l)*dd(2,jj)
286
287            end do ! jj
288
289!           Compute tangent stiffness
290
291            j1 = 0
292            do j  = 1,i
293
294              s(i1+1,j1+1) = s(i1+1,j1+1) + bbd(1,1)*shp2(1,j,l)
295     &                                    + bbd(3,1)*shpr(  j  )
296     &                                    + bbd(4,1)*shp2(2,j,l)
297
298              s(i1+2,j1+1) = s(i1+2,j1+1) + bbd(1,2)*shp2(1,j,l)
299     &                                    + bbd(3,2)*shpr(  j  )
300     &                                    + bbd(4,2)*shp2(2,j,l)
301
302              s(i1+1,j1+2) = s(i1+1,j1+2) + bbd(4,1)*shp2(1,j,l)
303     &                                    + bbd(2,1)*shp2(2,j,l)
304
305              s(i1+2,j1+2) = s(i1+2,j1+2) + bbd(4,2)*shp2(1,j,l)
306     &                                    + bbd(2,2)*shp2(2,j,l)
307
308              j1 = j1 + ndf
309            end do
310
311            i1 = i1 + ndf
312          end  do
313
314        endif ! end of tangent
315
316        endif ! end of isw options
317
318        nn = nn + nhv
319
320      end do
321
322      if(isw.eq.3) then
323
324!       Compute upper part by symmetry
325
326        do j = 1,nst
327          do i = 1,j
328            s(i,j) = s(j,i)
329          end do
330        end do
331
332        if(debug) then
333          call mprint(r,ndf,nel,ndf,'Resid')
334          call mprint(s,nst,nst,nst,'Stiff')
335        endif
336
337      endif
338
339      return
340
341!     OUTPUT STRESSES
342
343   4  xx1  = 0.d0
344      xx2  = 0.d0
345      xx3  = 0.d0
346      do i = 1,6
347        sigv(i) = 0.0d0
348      end do ! i
349
350!     LOOP OVER GAUSS POINTS
351
352      nn = 0
353
354      do l = 1,lint
355
356!       Compute Cauchy stresses and spatial tangent tensor at t-n+1
357
358        call modlfd(l,d,f(1,1,l),df(1,1,l),detf(1,l),ta,
359     &             hr(nn+nh1),hr(nn+nh2),nhv,istrt,ds,sigl(1,l),bb,
360     &             .false.,isw)
361
362        yy1 = 1.d0/dble(nel)
363        do i=1,nel
364          tempi = yy1 * shp2(3,i,l)
365          xx1   = xx1 + tempi*xl(1,i)
366          xx2   = xx2 + tempi*xl(2,i)
367        end do
368
369!       Compute average stresses and jacobian for printing
370
371        do i = 1,4
372          sigv(i) = sigv(i) + yy1 * sigl(i,l)
373        end do
374
375        nn = nn + nhv
376
377      end do
378
379!     OUTPUT STRESSES
380
381      if(isw.eq.4) then
382        mct = mct - 2
383        if(mct.le.0) then
384          write(iow,2001) o,head
385          if(ior.lt.0) write(*,2001) o,head
386          mct = 50
387        endif
388
389        write(iow,2002) n,ma,(sigv(jj),jj=1,6),xx1,xx2,xx3
390        if(ior.lt.0) then
391          write(*,2002) n,ma,(sigv(jj),jj=1,6),xx1,xx2,xx3
392        end if
393
394      elseif(isw.eq.8) then
395
396!       Project stress values to nodes
397
398        call stcn2d(ix,sigl,shp2,jac,hr(nph),hr(nph+numnp),hr(ner),
399     &              erav,lint,nel,9,numnp)
400
401      end if
402
403!     Format statements
404
4052001  format(a1,20a4//5x,'Element Stresses'//'  Elmt  Matl',
406     1   '  11-stress  22-stress  33-stress  12-stress',
407     2   '  23-stress  13-stress'/16x,'1-coord    2-coord    3-coord ')
408
4092002  format(2i6,1p,6e11.3/12x,1p,6e11.3/12x,0p,3f11.5)
410
411      end
412