1!$Id:$
2      subroutine tjac3d ( ss , xl, ndm, nel, shp, detj )
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: Compute jacobian determinant of tetrahedral element
11
12!      Inputs:
13!         ss(3)     - Natural coordinate point
14!         xl(ndm,*) - Nodal coordinates for element
15!         ndm       - Spatial dimension of mesh
16!         nel       - Number of nodes on element
17
18!      Outputs:
19!         shp(3,*)  - Derivatives of shape functions for element
20!                     in natural coordinates.
21!         detj      - Jacobian determinant at point
22!-----[--.----+----.----+----.-----------------------------------------]
23
24      implicit  none
25
26      integer   ndm, nel, i, j, k
27      real*8    detj, ss4
28
29      real*8    ss(3), xl(ndm,*), shp(3,11), xs(3,3)
30
31      save
32
33!     Compute shape functions and their derivatives
34
35      ss4 = 1.d0 - ss(1) - ss(2) - ss(3)
36      if(nel.eq.4) then
37        shp(1, 1) = 1.d0
38        shp(1, 2) = 0.d0
39        shp(1, 3) = 0.d0
40        shp(1, 4) =-1.d0
41        shp(2, 1) = 0.d0
42        shp(2, 2) = 1.d0
43        shp(2, 3) = 0.d0
44        shp(2, 4) =-1.d0
45        shp(3, 1) = 0.d0
46        shp(3, 2) = 0.d0
47        shp(3, 3) = 1.d0
48        shp(3, 4) =-1.d0
49      elseif(nel.ge.10) then
50
51        shp(1, 1) = 4.d0*ss(1) - 1.d0
52        shp(1, 2) = 0.d0
53        shp(1, 3) = 0.d0
54        shp(1, 4) =-4.d0*ss4 + 1.d0
55        shp(1, 5) = 4.d0*ss(2)
56        shp(1, 6) = 0.d0
57        shp(1, 7) = 4.d0*ss(3)
58        shp(1, 8) = 4.d0*(ss4 - ss(1))
59        shp(1, 9) =-4.d0*ss(2)
60        shp(1,10) =-4.d0*ss(3)
61
62        shp(2, 1) = 0.d0
63        shp(2, 2) = 4.d0*ss(2) - 1.d0
64        shp(2, 3) = 0.d0
65        shp(2, 4) =-4.d0*ss4 + 1.d0
66        shp(2, 5) = 4.d0*ss(1)
67        shp(2, 6) = 4.d0*ss(3)
68        shp(2, 7) = 0.d0
69        shp(2, 8) =-4.d0*ss(1)
70        shp(2, 9) = 4.d0*(ss4 - ss(2))
71        shp(2,10) =-4.d0*ss(3)
72
73        shp(3, 1) = 0.d0
74        shp(3, 2) = 0.d0
75        shp(3, 3) = 4.d0*ss(3) - 1.d0
76        shp(3, 4) =-4.d0*ss4 + 1.d0
77        shp(3, 5) = 0.d0
78        shp(3, 6) = 4.d0*ss(2)
79        shp(3, 7) = 4.d0*ss(1)
80        shp(3, 8) =-4.d0*ss(1)
81        shp(3, 9) =-4.d0*ss(2)
82        shp(3,10) = 4.d0*(ss4 - ss(3))
83        if(nel.eq.11) then
84          shp(1,11) = 256.d0*ss(2)*ss(3)*(ss4 - ss(1))
85          shp(2,11) = 256.d0*ss(1)*ss(3)*(ss4 - ss(2))
86          shp(3,11) = 256.d0*ss(1)*ss(2)*(ss4 - ss(3))
87        endif
88
89      endif
90
91!     Compute jacobian matrix
92
93      do i = 1,3
94        do j = 1,3
95          xs(i,j) = 0.0d0
96          do k = 1,nel
97            xs(i,j) = xs(i,j) - shp(i,k)*xl(j,k)
98          end do
99        end do
100      end do
101
102!     Compute jacobian determinant
103
104      detj = xs(1,1)*(xs(2,2)*xs(3,3) - xs(2,3)*xs(3,2))
105     &     + xs(1,2)*(xs(2,3)*xs(3,1) - xs(2,1)*xs(3,3))
106     &     + xs(1,3)*(xs(2,1)*xs(3,2) - xs(2,2)*xs(3,1))
107
108      end
109