FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
dynamic_mat_ass_couple.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!-------------------------------------------------------------------------------
7 use m_fstr
8 private :: area_of_triangle
9 private :: area_of_triangle2
10 private :: area_of_squre
11 private :: node_on_surface
12
13 integer(kind=kint), private :: icall = 0
14contains
15
16 subroutine dynamic_mat_ass_couple( hecMESH, hecMAT, fstrSOLID, fstrCPL )
17 use m_fstr
18 implicit none
19 type(hecmwst_matrix) :: hecMAT
20 type(hecmwst_local_mesh) :: hecMESH
21 type(fstr_solid) :: fstrSOLID
22 type(fstr_couple) :: fstrCPL
23 ! local
24 integer(kind=kint) :: ig0, ig, is, ie, ik
25 integer(kind=kint) :: i, j, count
26 integer(kind=kint) :: eid, sid, etype
27 integer(kind=kint) :: node(20)
28 integer(kind=kint) :: node_n
29 real(kind=kreal) :: px, py, pz ! traction on surface
30 real(kind=kreal) :: vx, vy, vz ! force on vertex
31
32 ! ================== modified by K. Tagami ============= 2010/02/25 ====
33 !!!! real(kind=kreal) :: xx(4), yy(4), zz(4)
34 integer, parameter :: MaxNumNodesOnSurf = 8
35 !parameter(MaxNumNodesOnSurf=8)
36 real(kind=kreal) :: xx(maxnumnodesonsurf)
37 real(kind=kreal) :: yy(maxnumnodesonsurf)
38 real(kind=kreal) :: zz(maxnumnodesonsurf)
39 ! ================================================================
40
41 real(kind=kreal) :: area, wg
42 integer(kind=kint) :: ierr
43
44 ! ============= added by K. Tagami =========for debug ====== 2010/03/02 ====
45 integer( kind=kint ) :: node_global
46 ! ==========================================================
47
48 icall = icall + 1
49 ierr = 0
50
51 do ig0= 1, fstrsolid%COUPLE_ngrp_tot
52 ig= fstrsolid%COUPLE_ngrp_ID(ig0)
53 is= hecmesh%surf_group%grp_index(ig-1) + 1
54 ie= hecmesh%surf_group%grp_index(ig )
55 do ik= is, ie
56 sid = hecmesh%surf_group%grp_item(2*ik)
57 eid = hecmesh%surf_group%grp_item(2*ik-1)
58 etype = hecmesh%elem_type(eid)
59 call node_on_surface( hecmesh, etype, eid, sid, node, node_n )
60
61 !-------------------------------------------------------------
62 ! traction on surface
63 !-------------------------------------------------------------
64 px = 0
65 py = 0
66 pz = 0
67 count=0
68 do i=1, node_n
69 j=3*fstrcpl%index( node(i) )
70
71 ! ============= added by K. Tagami =========for debug ====== 2010/03/02 ====
72 ! node_global = hecMESH%global_node_ID( node(i) )
73 ! ======================================================================
74 !
75 if( icall == 1 ) then
76 if( j<=0 ) then
77 write(idbg,'(a,i0,a)') "dynamic_mat_ass_couple: traction for node ", &
78 & hecmesh%global_node_ID(node(i)), " on coupling surface not found"
79 ! ============ added by K. Tagami ============== for debug ==== 2010/03/02
80 ! write(IDBG,*) "local : ", node(i), " Global :", node_global, " not found"
81 ! ==============================================================
82 ierr = 1
83 cycle
84 else
85 write(idbg,'(a,i0,a)') "dynamic_mat_ass_couple: traction for node ", &
86 & hecmesh%global_node_ID(node(i)), " on coupling surface OK"
87 ! ============ added by K. Tagami ============== for debug ==== 2010/03/02
88 ! write(IDBG,*) "local : ", node(i), " Global :", node_global, " OK"
89 ! ==============================================================
90 endif
91 endif
92 px = px + fstrcpl%trac(j-2)
93 py = py + fstrcpl%trac(j-1)
94 pz = pz + fstrcpl%trac(j )
95 count = count +1
96 end do
97 if( count == 0 ) cycle
98 px = px / count
99 py = py / count
100 pz = pz / count
101 !-------------------------------------------------------------
102 ! force on vertex
103 !-------------------------------------------------------------
104 do i=1, node_n
105 j=3*node(i)
106 xx(i)=hecmesh%node(j-2)
107 yy(i)=hecmesh%node(j-1)
108 zz(i)=hecmesh%node(j )
109 end do
110 area = 0.0d0
111 if( node_n == 3 ) then
112 area = area_of_triangle( xx,yy,zz )
113 else if( node_n == 4 ) then
114 area = area_of_squre( xx,yy,zz )
115 else if( node_n == 6 ) then
116 area = area_of_triangle2(xx,yy,zz)
117 else
118 write(*,*) "#Error : in FSTR_MAT_ASS_COUPLE "
119 call hecmw_abort( hecmw_comm_get_comm())
120 end if
121 wg = area / node_n
122
123 vx = px * wg
124 vy = py * wg
125 vz = pz * wg
126 !-------------------------------------------------------------
127 ! add in B
128 !-------------------------------------------------------------
129 do i=1, node_n
130 j=3*node(i)
131
132 hecmat%B(j-2)=hecmat%B(j-2) + vx
133 hecmat%B(j-1)=hecmat%B(j-1) + vy
134 hecmat%B(j )=hecmat%B(j ) + vz
135 end do
136 end do
137 end do
138 call hecmw_barrier(hecmesh)
139 if( ierr == 1 ) then
140 write(*,*) "#Error : in FSTR_MAT_ASS_COUPLE"
141 call hecmw_abort( hecmw_comm_get_comm())
142 endif
143 end subroutine dynamic_mat_ass_couple
144
145 !==============================================================================
146 ! CALC AREA
147 !==============================================================================
148
149 function area_of_triangle( XX,YY,ZZ )
150 implicit none
151 real(kind=kreal) :: xx(*), yy(*), zz(*)
152 real(kind=kreal) :: v1x, v1y, v1z
153 real(kind=kreal) :: v2x, v2y, v2z
154 real(kind=kreal) :: v3x, v3y, v3z
155 real(kind=kreal) :: area_of_triangle
156
157 v1x=xx(2)-xx(1)
158 v1y=yy(2)-yy(1)
159 v1z=zz(2)-zz(1)
160 v2x=xx(3)-xx(1)
161 v2y=yy(3)-yy(1)
162 v2z=zz(3)-zz(1)
163 v3x= v1y*v2z-v1z*v2y
164 v3y=-v1x*v2z+v1z*v2x
165 v3z= v1x*v2y-v1y*v2x
166 area_of_triangle = sqrt( v3x*v3x + v3y*v3y + v3z*v3z )*0.5
167 end function area_of_triangle
168
169 function area_of_triangle2( XX,YY,ZZ )
170 implicit none
171 real(kind=kreal) :: xx(*), yy(*), zz(*)
172 real(kind=kreal) :: area_of_triangle2
173 real(kind=kreal) :: x(3), y(3), z(3)
174
175 x(1)=xx(1)
176 x(2)=xx(3)
177 x(3)=xx(5)
178 y(1)=yy(1)
179 y(2)=yy(3)
180 y(3)=yy(5)
181 z(1)=zz(1)
182 z(2)=zz(3)
183 z(3)=zz(5)
184 area_of_triangle2 = area_of_triangle(x,y,z)
185 end function area_of_triangle2
186
187 function area_of_triangle2a( XX,YY,ZZ )
188 implicit none
189 real(kind=kreal) xx(*),yy(*),zz(*)
190 real(kind=kreal) area_of_triangle2a
191 real(kind=kreal) x(3),y(3),z(3)
192 x(1)=xx(1)
193 x(2)=xx(2)
194 x(3)=xx(3)
195 y(1)=yy(1)
196 y(2)=yy(2)
197 y(3)=yy(3)
198 z(1)=zz(1)
199 z(2)=zz(2)
200 z(3)=zz(3)
201 area_of_triangle2a = area_of_triangle(x,y,z)
202 end function area_of_triangle2a
203
204
205 function area_of_squre ( XX,YY,ZZ )
206 implicit none
207 ! I/F VARIABLES
208 real(kind=kreal) :: xx(*), yy(*), zz(*)
209 real(kind=kreal) :: area_of_squre
210 ! LOCAL VARIABLES
211 integer(kind=kint), parameter :: nn = 8
212 integer(kind=kint), parameter :: ng = 2
213 !parameter(NN=8, NG=2)
214 real(kind=kreal) :: h(nn), hr(nn), hs(nn), ht(nn)
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, det, wg
217 integer(kind=kint) :: ig1, ig2, lx, ly, lz, i
218 real(kind=kreal) :: vx, vy, vz, xcod, ycod, zcod
219 real(kind=kreal) :: ax, ay, az, rx, ry, rz, hx, hy, hz, val
220 real(kind=kreal) :: phx, phy, phz
221 real(kind=kreal) :: g1x, g1y, g1z
222 real(kind=kreal) :: g2x, g2y, g2z
223 real(kind=kreal) :: g3x, g3y, g3z
224 real(kind=kreal) :: xsum, coefx, coefy, coefz
225 real(kind=kreal) :: area, xg(2), wgt(2)
226 data wgt/1.0, 1.0/
227 data xg/-0.5773502691896, 0.5773502691896/
228
229 area = 0.0
230 ! INTEGRATION OVER SURFACE
231 do ig2=1,ng
232 si=xg(ig2)
233 do ig1=1,ng
234 ri=xg(ig1)
235 h(1)=0.25*(1.0-ri)*(1.0-si)
236 h(2)=0.25*(1.0+ri)*(1.0-si)
237 h(3)=0.25*(1.0+ri)*(1.0+si)
238 h(4)=0.25*(1.0-ri)*(1.0+si)
239 hr(1)=-.25*(1.0-si)
240 hr(2)= .25*(1.0-si)
241 hr(3)= .25*(1.0+si)
242 hr(4)=-.25*(1.0+si)
243 hs(1)=-.25*(1.0-ri)
244 hs(2)=-.25*(1.0+ri)
245 hs(3)= .25*(1.0+ri)
246 hs(4)= .25*(1.0-ri)
247 g1x=0.0
248 g1y=0.0
249 g1z=0.0
250 g2x=0.0
251 g2y=0.0
252 g2z=0.0
253 do i=1,nn
254 g1x=g1x+hr(i)*xx(i)
255 g1y=g1y+hr(i)*yy(i)
256 g1z=g1z+hr(i)*zz(i)
257 g2x=g2x+hs(i)*xx(i)
258 g2y=g2y+hs(i)*yy(i)
259 g2z=g2z+hs(i)*zz(i)
260 enddo
261 g3x=g1y*g2z-g1z*g2y
262 g3y=g1z*g2x-g1x*g2z
263 g3z=g1x*g2y-g1y*g2x
264 xsum=dsqrt(g3x**2+g3y**2+g3z**2)
265 g3x=g3x/xsum
266 g3y=g3y/xsum
267 g3z=g3z/xsum
268 !JACOBI MATRIX
269 xj11=g1x
270 xj12=g1y
271 xj13=g1z
272 xj21=g2x
273 xj22=g2y
274 xj23=g2z
275 xj31=g3x
276 xj32=g3y
277 xj33=g3z
278 !DETERMINANT OF JACOBIAN
279 det=xj11*xj22*xj33 &
280 +xj12*xj23*xj31 &
281 +xj13*xj21*xj32 &
282 -xj13*xj22*xj31 &
283 -xj12*xj21*xj33 &
284 -xj11*xj23*xj32
285 wg=wgt(ig1)*wgt(ig2)*det
286 do i = 1, nn
287 area = area + h(i)*wg
288 enddo
289 enddo
290 enddo
291 area_of_squre = area;
292 end function area_of_squre
293
294 !==============================================================================
295 ! GET NODES ON SURFACE
296 !==============================================================================
297
298 subroutine node_on_surface( hecMESH, etype, eid, sid, node, node_n )
299 implicit none
300 type(hecmwst_local_mesh) :: hecmesh
301 integer(kind=kint) :: eid
302 integer(kind=kint) :: sid
303 integer(kind=kint) :: node(*)
304 integer(kind=kint) :: node_n
305 ! local parameters
306 integer(kind=kint) :: etype
307 integer(kind=kint) :: is
308 integer(kind=kint) :: tbl341(3, 4)
309 integer(kind=kint) :: tbl342(6, 4)
310 integer(kind=kint) :: tbl361(4, 6)
311 !! vertex id tables by definition of fstr
312 data tbl341 / 1,2,3, 1,2,4, 2,3,4, 3,1,4 /
313 data tbl342 / 1,5,2,6,3,7, 1,5,2,9,4,8, 2,6,3,10,4,9, 3,7,1,10,4,8 /
314 data tbl361 / 1,2,3,4, 5,6,7,8, 1,2,6,5, 2,3,7,6, 3,4,8,7, 4,1,5,8 /
315
316 ! =============== added by K. Tagami ======== experimental ==== 2010/02/25 ==
317 integer(kind=kint) :: tbl351(4,5)
318 data tbl351 / 1,2,3,1, 4,5,6,4, 1,2,5,4, 2,3,6,5, 3,1,4,6 /
319
320 integer(kind=kint) :: tbl362(8,6)
321 data tbl362 / 1, 9, 2,10, 3,11, 4,12, &
322 & 5,13, 6,14, 7,15, 8,16, &
323 & 1, 9, 2,18, 6,13, 5,17, &
324 & 2,10, 3,19, 7,14, 6,18, &
325 & 3,11, 4,20, 8,15, 7,19, &
326 & 4,12, 1,17, 5,16, 8,20 /
327 ! =================================================================
328
329 etype = hecmesh%elem_type( eid )
330 is = hecmesh%elem_node_index( eid-1 )
331 select case( etype )
332 case ( 341 )
333 node_n = 3
334 node(1) = hecmesh%elem_node_item (is + tbl341(1,sid) )
335 node(2) = hecmesh%elem_node_item (is + tbl341(2,sid) )
336 node(3) = hecmesh%elem_node_item (is + tbl341(3,sid) )
337 case ( 342 )
338 node_n = 6
339 node(1) = hecmesh%elem_node_item (is + tbl342(1,sid) )
340 node(2) = hecmesh%elem_node_item (is + tbl342(2,sid) )
341 node(3) = hecmesh%elem_node_item (is + tbl342(3,sid) )
342 node(4) = hecmesh%elem_node_item (is + tbl342(4,sid) )
343 node(5) = hecmesh%elem_node_item (is + tbl342(5,sid) )
344 node(6) = hecmesh%elem_node_item (is + tbl342(6,sid) )
345 case ( 361 )
346 node_n = 4
347 node(1) = hecmesh%elem_node_item (is + tbl361(1,sid) )
348 node(2) = hecmesh%elem_node_item (is + tbl361(2,sid) )
349 node(3) = hecmesh%elem_node_item (is + tbl361(3,sid) )
350 node(4) = hecmesh%elem_node_item (is + tbl361(4,sid) )
351
352 ! ==================== Added by K. Tagami ========== experimental == 2010/02/25
353 case ( 351 )
354 node_n = 4
355 node(1) = hecmesh%elem_node_item (is + tbl351(1,sid) )
356 node(2) = hecmesh%elem_node_item (is + tbl351(2,sid) )
357 node(3) = hecmesh%elem_node_item (is + tbl351(3,sid) )
358 node(4) = hecmesh%elem_node_item (is + tbl351(4,sid) )
359 if ( node(1) == node(4) ) then
360 node_n = 3
361 endif
362 ! --------------- K. Tagami -------------- under construction ---
363 ! case ( 362 )
364 ! node_n = 8
365 ! node(1) = hecMESH%elem_node_item (is + tbl362(1,sid) )
366 ! node(2) = hecMESH%elem_node_item (is + tbl362(2,sid) )
367 ! node(3) = hecMESH%elem_node_item (is + tbl362(3,sid) )
368 ! node(4) = hecMESH%elem_node_item (is + tbl362(4,sid) )
369 ! node(5) = hecMESH%elem_node_item (is + tbl362(5,sid) )
370 ! node(6) = hecMESH%elem_node_item (is + tbl362(6,sid) )
371 ! node(7) = hecMESH%elem_node_item (is + tbl362(7,sid) )
372 ! node(8) = hecMESH%elem_node_item (is + tbl362(8,sid) )
373 ! --------------------------------------------------------------
374 ! ================================================================
375 case default
376 if( hecmesh%my_rank==0 ) then
377 write(*,*) '##Error: not supported element type for coupling analysis ', etype
378 ! =================== modified by K. Tagami ===== 2010/02/25
379 ! write(*,*) ' This version supports element type 341,342 and 361 only.'
380 write(*,*) ' This version supports element type 341,342, 351, and 361 only.'
381 ! ----------------------------------------------------------
382 endif
383 call hecmw_abort( hecmw_comm_get_comm() )
384 end select
385
386 end subroutine node_on_surface
387
389
390
391
This module contains functions relates to coupling analysis.
real(kind=kreal) function area_of_triangle2a(xx, yy, zz)
subroutine dynamic_mat_ass_couple(hecmesh, hecmat, fstrsolid, fstrcpl)
This module defined coomon data and basic structures for analysis.
Definition: m_fstr.f90:15
integer(kind=kint), parameter idbg
Definition: m_fstr.f90:95
Data for coupling analysis.
Definition: m_fstr.f90:580