FrontISTR  5.7.0
Large-scale structural analysis program with finit element method
hecmw_Jacob_361.f90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! Copyright (c) 2019 FrontISTR Commons
3 ! This software is released under the MIT License, see LICENSE.txt
4 !-------------------------------------------------------------------------------
6 
7 module hecmw_jacob361
8 contains
9 
10  subroutine hecmw_jacob_361 ( hecMESH, iElem, DET, W, N, NX, NY, NZ)
11  use hecmw_util
12 
13  implicit none
14 
15  type(hecmwst_local_mesh):: hecmesh
16  integer(kind=kint):: ielem
17  real(kind=kreal):: det
18  real(kind=kreal):: w(8), n(8,8),nx(8,8),ny(8,8),nz(8,8)
19 
20  integer(kind=kint):: ilocal, j, jj
21  integer(kind=kint):: lx, ly, lz
22  real(kind=kreal):: dum
23  real(kind=kreal):: xx(8), yy(8), zz(8)
24  real(kind=kreal):: xg(2), wgt(2), hr(8), hs(8), ht(8)
25  real(kind=kreal):: h(2,2,2,8),bx(2,2,2,8),by(2,2,2,8),bz(2,2,2,8)
26  real(kind=kreal):: ri, si, ti, rp, sp, tp, rm, sm, tm
27  real(kind=kreal):: xj11,xj12,xj13,xj21,xj22,xj23,xj31,xj32,xj33
28  real(kind=kreal):: xji11,xji12,xji13,xji21,xji22,xji23,xji31,xji32,xji33
29 
30  data wgt/1.0d0,1.0d0/
31  data xg/-0.5773502691896d0, 0.5773502691896d0/
32  !C
33 
34  do j = 1, 8
35  jj = 8 - j
36  ilocal = hecmesh%elem_node_item ( 8*ielem -jj )
37  xx(j) = hecmesh%node( ilocal*3 -2 )
38  yy(j) = hecmesh%node( ilocal*3 -1 )
39  zz(j) = hecmesh%node( ilocal*3 )
40  w(j) = 1.0d0
41  end do
42 
43  do lx=1,2
44  ri=xg(lx)
45  do ly=1,2
46  si=xg(ly)
47  do lz=1,2
48  ti=xg(lz)
49  !C
50  rp=1.0+ri
51  sp=1.0+si
52  tp=1.0+ti
53  rm=1.0-ri
54  sm=1.0-si
55  tm=1.0-ti
56  !C
57  !C*INTERPOLATION FUNCTION
58  h(lx,ly,lz,1)=0.125*rm*sm*tm
59  h(lx,ly,lz,2)=0.125*rp*sm*tm
60  h(lx,ly,lz,3)=0.125*rp*sp*tm
61  h(lx,ly,lz,4)=0.125*rm*sp*tm
62  h(lx,ly,lz,5)=0.125*rm*sm*tp
63  h(lx,ly,lz,6)=0.125*rp*sm*tp
64  h(lx,ly,lz,7)=0.125*rp*sp*tp
65  h(lx,ly,lz,8)=0.125*rm*sp*tp
66  !C
67  !C*DERIVATIVE OF INTERPOLATION FUNCTION
68  !C* FOR R-COORDINATE
69  hr(1)=-.125*sm*tm
70  hr(2)= .125*sm*tm
71  hr(3)= .125*sp*tm
72  hr(4)=-.125*sp*tm
73  hr(5)=-.125*sm*tp
74  hr(6)= .125*sm*tp
75  hr(7)= .125*sp*tp
76  hr(8)=-.125*sp*tp
77  !C* FOR S-COORDINATE
78  hs(1)=-.125*rm*tm
79  hs(2)=-.125*rp*tm
80  hs(3)= .125*rp*tm
81  hs(4)= .125*rm*tm
82  hs(5)=-.125*rm*tp
83  hs(6)=-.125*rp*tp
84  hs(7)= .125*rp*tp
85  hs(8)= .125*rm*tp
86  !C* FOR T-COORDINATE
87  ht(1)=-.125*rm*sm
88  ht(2)=-.125*rp*sm
89  ht(3)=-.125*rp*sp
90  ht(4)=-.125*rm*sp
91  ht(5)= .125*rm*sm
92  ht(6)= .125*rp*sm
93  ht(7)= .125*rp*sp
94  ht(8)= .125*rm*sp
95  !C
96  !C*JACOBI MATRIX
97  xj11=hr(1)*xx(1)+hr(2)*xx(2)+hr(3)*xx(3)+hr(4)*xx(4) &
98  & +hr(5)*xx(5)+hr(6)*xx(6)+hr(7)*xx(7)+hr(8)*xx(8)
99  xj21=hs(1)*xx(1)+hs(2)*xx(2)+hs(3)*xx(3)+hs(4)*xx(4) &
100  & +hs(5)*xx(5)+hs(6)*xx(6)+hs(7)*xx(7)+hs(8)*xx(8)
101  xj31=ht(1)*xx(1)+ht(2)*xx(2)+ht(3)*xx(3)+ht(4)*xx(4) &
102  & +ht(5)*xx(5)+ht(6)*xx(6)+ht(7)*xx(7)+ht(8)*xx(8)
103  !C
104  xj12=hr(1)*yy(1)+hr(2)*yy(2)+hr(3)*yy(3)+hr(4)*yy(4) &
105  & +hr(5)*yy(5)+hr(6)*yy(6)+hr(7)*yy(7)+hr(8)*yy(8)
106  xj22=hs(1)*yy(1)+hs(2)*yy(2)+hs(3)*yy(3)+hs(4)*yy(4) &
107  & +hs(5)*yy(5)+hs(6)*yy(6)+hs(7)*yy(7)+hs(8)*yy(8)
108  xj32=ht(1)*yy(1)+ht(2)*yy(2)+ht(3)*yy(3)+ht(4)*yy(4) &
109  & +ht(5)*yy(5)+ht(6)*yy(6)+ht(7)*yy(7)+ht(8)*yy(8)
110  !C
111  xj13=hr(1)*zz(1)+hr(2)*zz(2)+hr(3)*zz(3)+hr(4)*zz(4) &
112  & +hr(5)*zz(5)+hr(6)*zz(6)+hr(7)*zz(7)+hr(8)*zz(8)
113  xj23=hs(1)*zz(1)+hs(2)*zz(2)+hs(3)*zz(3)+hs(4)*zz(4) &
114  & +hs(5)*zz(5)+hs(6)*zz(6)+hs(7)*zz(7)+hs(8)*zz(8)
115  xj33=ht(1)*zz(1)+ht(2)*zz(2)+ht(3)*zz(3)+ht(4)*zz(4) &
116  & +ht(5)*zz(5)+ht(6)*zz(6)+ht(7)*zz(7)+ht(8)*zz(8)
117  !C
118  !C*DETERMINANT OF JACOBIAN
119  !C
120  det=xj11*xj22*xj33 &
121  & +xj12*xj23*xj31 &
122  & +xj13*xj21*xj32 &
123  & -xj13*xj22*xj31 &
124  & -xj12*xj21*xj33 &
125  & -xj11*xj23*xj32
126  !C
127  !C* INVERSION OF JACOBIAN
128  !C
129  dum= -1.0/det
130  !C
131  xji11=dum*( xj22*xj33-xj23*xj32)
132  xji21=dum*(-xj21*xj33+xj23*xj31)
133  xji31=dum*( xj21*xj32-xj22*xj31)
134  xji12=dum*(-xj12*xj33+xj13*xj32)
135  xji22=dum*( xj11*xj33-xj13*xj31)
136  xji32=dum*(-xj11*xj32+xj12*xj31)
137  xji13=dum*( xj12*xj23-xj13*xj22)
138  xji23=dum*(-xj11*xj23+xj13*xj21)
139  xji33=dum*( xj11*xj22-xj12*xj21)
140  !C
141  do j=1, 8
142  bx(lx,ly,lz,j)=xji11*hr(j)+xji12*hs(j)+xji13*ht(j)
143  by(lx,ly,lz,j)=xji21*hr(j)+xji22*hs(j)+xji23*ht(j)
144  bz(lx,ly,lz,j)=xji31*hr(j)+xji32*hs(j)+xji33*ht(j)
145  end do
146  !C
147  !C
148  end do
149  end do
150  end do
151 
152  j = 1
153  do lx = 1, 2
154  do ly = 1, 2
155  do lz = 1, 2
156 
157  n(j,1) = h(lx,ly,lz,1)
158  n(j,2) = h(lx,ly,lz,2)
159  n(j,3) = h(lx,ly,lz,3)
160  n(j,4) = h(lx,ly,lz,4)
161  n(j,5) = h(lx,ly,lz,5)
162  n(j,6) = h(lx,ly,lz,6)
163  n(j,7) = h(lx,ly,lz,7)
164  n(j,8) = h(lx,ly,lz,8)
165 
166  nx(j,1) = bx(lx,ly,lz,1)
167  nx(j,2) = bx(lx,ly,lz,2)
168  nx(j,3) = bx(lx,ly,lz,3)
169  nx(j,4) = bx(lx,ly,lz,4)
170  nx(j,5) = bx(lx,ly,lz,5)
171  nx(j,6) = bx(lx,ly,lz,6)
172  nx(j,7) = bx(lx,ly,lz,7)
173  nx(j,8) = bx(lx,ly,lz,8)
174 
175  ny(j,1) = by(lx,ly,lz,1)
176  ny(j,2) = by(lx,ly,lz,2)
177  ny(j,3) = by(lx,ly,lz,3)
178  ny(j,4) = by(lx,ly,lz,4)
179  ny(j,5) = by(lx,ly,lz,5)
180  ny(j,6) = by(lx,ly,lz,6)
181  ny(j,7) = by(lx,ly,lz,7)
182  ny(j,8) = by(lx,ly,lz,8)
183 
184  nz(j,1) = bz(lx,ly,lz,1)
185  nz(j,2) = bz(lx,ly,lz,2)
186  nz(j,3) = bz(lx,ly,lz,3)
187  nz(j,4) = bz(lx,ly,lz,4)
188  nz(j,5) = bz(lx,ly,lz,5)
189  nz(j,6) = bz(lx,ly,lz,6)
190  nz(j,7) = bz(lx,ly,lz,7)
191  nz(j,8) = bz(lx,ly,lz,8)
192 
193  j = j + 1
194  end do
195  end do
196  end do
197 
198  !C
199  !C
200  end subroutine hecmw_jacob_361
201 
202  !*----------------------------------------------------------------------*
203  subroutine hecmw_jacob_361_surface( hecMESH, iElem, iSurface, NOD, DET, H )
204  !*----------------------------------------------------------------------*
205  use hecmw_util
206  implicit none
207 
208  type(hecmwst_local_mesh):: hecmesh
209  integer(kind=kint):: ielem, isurface, nod(4)
210  real(kind=kreal):: det
211 
212  real(kind=kreal):: xx(4),yy(4),zz(4)
213  real(kind=kreal):: h(2,2,4),hr(4),hs(4),pl(4)
214  real(kind=kreal):: xg(2),wgt(2)
215  real(kind=kreal):: ri,si,ti,rp,sp,tp,rm,sm,tm
216  real(kind=kreal):: xj11,xj21,xj31,xj12,xj22,xj32,xj13,xj23,xj33,wg
217  integer(kind=kint):: iNode
218  integer(kind=kint):: IG1,IG2,LX,LY,LZ,I
219  real(kind=kreal):: vx,vy,vz
220  real(kind=kreal):: g1x,g1y,g1z
221  real(kind=kreal):: g2x,g2y,g2z
222  real(kind=kreal):: g3x,g3y,g3z
223 
224  !**************************
225  !* GAUSS INTEGRATION POINT
226  !**************************
227  data xg/-0.5773502691896d0,0.5773502691896d0/
228  data wgt/1.0d0,1.0d0/
229  !*
230  !* Set node numbers for selected surface
231  !*
232  if( isurface.EQ.1 ) then
233  nod(1) = 1
234  nod(2) = 2
235  nod(3) = 3
236  nod(4) = 4
237  else if( isurface.EQ.2 ) then
238  nod(1) = 8
239  nod(2) = 7
240  nod(3) = 6
241  nod(4) = 5
242  else if( isurface.EQ.3 ) then
243  nod(1) = 5
244  nod(2) = 6
245  nod(3) = 2
246  nod(4) = 1
247  else if( isurface.EQ.4 ) then
248  nod(1) = 6
249  nod(2) = 7
250  nod(3) = 3
251  nod(4) = 2
252  else if( isurface.EQ.5 ) then
253  nod(1) = 7
254  nod(2) = 8
255  nod(3) = 4
256  nod(4) = 3
257  else if( isurface.EQ.6 ) then
258  nod(1) = 8
259  nod(2) = 5
260  nod(3) = 1
261  nod(4) = 4
262  endif
263 
264  do i = 1, 4
265  inode = hecmesh%elem_node_item( 8*(ielem-1) + nod(i) )
266  xx(i) = hecmesh%node( 3*inode -2 )
267  yy(i) = hecmesh%node( 3*inode -1 )
268  zz(i) = hecmesh%node( 3*inode )
269  end do
270 
271  !*
272  !*** SURFACE LOAD
273  !*
274  !* INTEGRATION OVER SURFACE
275  do ig2=1,2
276  si=xg(ig2)
277  do ig1=1,2
278  ri=xg(ig1)
279  h(ig1,ig2,1)=0.25*(1.0-ri)*(1.0-si)
280  h(ig1,ig2,2)=0.25*(1.0+ri)*(1.0-si)
281  h(ig1,ig2,3)=0.25*(1.0+ri)*(1.0+si)
282  h(ig1,ig2,4)=0.25*(1.0-ri)*(1.0+si)
283  hr(1)=-.25*(1.0-si)
284  hr(2)= .25*(1.0-si)
285  hr(3)= .25*(1.0+si)
286  hr(4)=-.25*(1.0+si)
287  hs(1)=-.25*(1.0-ri)
288  hs(2)=-.25*(1.0+ri)
289  hs(3)= .25*(1.0+ri)
290  hs(4)= .25*(1.0-ri)
291  g1x=0.0
292  g1y=0.0
293  g1z=0.0
294  g2x=0.0
295  g2y=0.0
296  g2z=0.0
297  do i=1,4
298  g1x=g1x+hr(i)*xx(i)
299  g1y=g1y+hr(i)*yy(i)
300  g1z=g1z+hr(i)*zz(i)
301  g2x=g2x+hs(i)*xx(i)
302  g2y=g2y+hs(i)*yy(i)
303  g2z=g2z+hs(i)*zz(i)
304  enddo
305  g3x=g1y*g2z-g1z*g2y
306  g3y=g1z*g2x-g1x*g2z
307  g3z=g1x*g2y-g1y*g2x
308  det=dsqrt(g3x**2+g3y**2+g3z**2)
309  enddo
310  enddo
311 
312  end subroutine hecmw_jacob_361_surface
313 
314 end module hecmw_jacob361
hecmw_util
I/O and Utility.
Definition: hecmw_util_f.F90:7
hecmw_util::hecmwst_local_mesh
Definition: hecmw_util_f.F90:234
hecmw_jacob361::hecmw_jacob_361_surface
subroutine hecmw_jacob_361_surface(hecMESH, iElem, iSurface, NOD, DET, H)
Definition: hecmw_Jacob_361.f90:204
hecmw_util::kreal
integer(kind=4), parameter kreal
Definition: hecmw_util_f.F90:16
hecmw_jacob361::hecmw_jacob_361
subroutine hecmw_jacob_361(hecMESH, iElem, DET, W, N, NX, NY, NZ)
Definition: hecmw_Jacob_361.f90:11
hecmw_jacob361
Jacobian calculation.
Definition: hecmw_Jacob_361.f90:7