xref: /phasta/phSolver/common/e3metric.f (revision 595995161822a203c8467e0e4a253d7bd7d6df32)
1c-----------------------------------------------------------------------
2c
3c  compute the metrics of the mapping from global to local
4c  coordinates and the jacobian of the mapping (weighted by
5c  the quadrature weight
6c
7c-----------------------------------------------------------------------
8      subroutine e3metric(  xl,      shgl,     dxidx,
9     &                      shg,     WdetJ)
10
11      include "common.h"
12
13      real*8     xl(npro,nenl,nsd),    shgl(npro,nsd,nshl),
14     &           dxidx(npro,nsd,nsd),  shg(npro,nshl,nsd),
15     &           WdetJ(npro)
16
17      real*8     dxdxi(npro,nsd,nsd),  tmp(npro)
18
19c
20c.... compute the deformation gradient
21c
22      dxdxi = zero
23c
24       do n = 1, nenl
25          dxdxi(:,1,1) = dxdxi(:,1,1) + xl(:,n,1) * shgl(:,1,n)
26          dxdxi(:,1,2) = dxdxi(:,1,2) + xl(:,n,1) * shgl(:,2,n)
27          dxdxi(:,1,3) = dxdxi(:,1,3) + xl(:,n,1) * shgl(:,3,n)
28          dxdxi(:,2,1) = dxdxi(:,2,1) + xl(:,n,2) * shgl(:,1,n)
29          dxdxi(:,2,2) = dxdxi(:,2,2) + xl(:,n,2) * shgl(:,2,n)
30          dxdxi(:,2,3) = dxdxi(:,2,3) + xl(:,n,2) * shgl(:,3,n)
31          dxdxi(:,3,1) = dxdxi(:,3,1) + xl(:,n,3) * shgl(:,1,n)
32          dxdxi(:,3,2) = dxdxi(:,3,2) + xl(:,n,3) * shgl(:,2,n)
33          dxdxi(:,3,3) = dxdxi(:,3,3) + xl(:,n,3) * shgl(:,3,n)
34       enddo
35c
36c.... compute the inverse of deformation gradient
37c
38       dxidx(:,1,1) =   dxdxi(:,2,2) * dxdxi(:,3,3)
39     &                - dxdxi(:,3,2) * dxdxi(:,2,3)
40       dxidx(:,1,2) =   dxdxi(:,3,2) * dxdxi(:,1,3)
41     &                - dxdxi(:,1,2) * dxdxi(:,3,3)
42       dxidx(:,1,3) =  dxdxi(:,1,2) * dxdxi(:,2,3)
43     &                - dxdxi(:,1,3) * dxdxi(:,2,2)
44       tmp          = one / ( dxidx(:,1,1) * dxdxi(:,1,1)
45     &                       + dxidx(:,1,2) * dxdxi(:,2,1)
46     &                       + dxidx(:,1,3) * dxdxi(:,3,1) )
47       dxidx(:,1,1) = dxidx(:,1,1) * tmp
48       dxidx(:,1,2) = dxidx(:,1,2) * tmp
49       dxidx(:,1,3) = dxidx(:,1,3) * tmp
50       dxidx(:,2,1) = (dxdxi(:,2,3) * dxdxi(:,3,1)
51     &                - dxdxi(:,2,1) * dxdxi(:,3,3)) * tmp
52       dxidx(:,2,2) = (dxdxi(:,1,1) * dxdxi(:,3,3)
53     &                - dxdxi(:,3,1) * dxdxi(:,1,3)) * tmp
54       dxidx(:,2,3) = (dxdxi(:,2,1) * dxdxi(:,1,3)
55     &                - dxdxi(:,1,1) * dxdxi(:,2,3)) * tmp
56       dxidx(:,3,1) = (dxdxi(:,2,1) * dxdxi(:,3,2)
57     &                - dxdxi(:,2,2) * dxdxi(:,3,1)) * tmp
58       dxidx(:,3,2) = (dxdxi(:,3,1) * dxdxi(:,1,2)
59     &                - dxdxi(:,1,1) * dxdxi(:,3,2)) * tmp
60       dxidx(:,3,3) = (dxdxi(:,1,1) * dxdxi(:,2,2)
61     &                - dxdxi(:,1,2) * dxdxi(:,2,1)) * tmp
62c
63       WdetJ = Qwt(lcsyst,intp) / tmp
64c
65c.... compute the global gradient of shape-functions
66c
67       do n = 1, nshl
68          shg(:,n,1) = shgl(:,1,n) * dxidx(:,1,1) +
69     &                 shgl(:,2,n) * dxidx(:,2,1) +
70     &                 shgl(:,3,n) * dxidx(:,3,1)
71          shg(:,n,2) = shgl(:,1,n) * dxidx(:,1,2) +
72     &                 shgl(:,2,n) * dxidx(:,2,2) +
73     &                 shgl(:,3,n) * dxidx(:,3,2)
74          shg(:,n,3) = shgl(:,1,n) * dxidx(:,1,3) +
75     &                 shgl(:,2,n) * dxidx(:,2,3) +
76     &                 shgl(:,3,n) * dxidx(:,3,3)
77       enddo
78
79       return
80       end
81
82
83
84