FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
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
8contains
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
314end module hecmw_jacob361
Jacobian calculation.
subroutine hecmw_jacob_361(hecmesh, ielem, det, w, n, nx, ny, nz)
subroutine hecmw_jacob_361_surface(hecmesh, ielem, isurface, nod, det, h)
I/O and Utility.
Definition: hecmw_util_f.F90:7
integer(kind=4), parameter kreal