FrontISTR  5.7.0
Large-scale structural analysis program with finit element method
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 !-------------------------------------------------------------------------------
6 module m_dynamic_mat_ass_couple
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
14 contains
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 
388 end module m_dynamic_mat_ass_couple
389 
390 
391 
m_dynamic_mat_ass_couple::area_of_triangle2a
real(kind=kreal) function area_of_triangle2a(XX, YY, ZZ)
Definition: dynamic_mat_ass_couple.f90:188
m_dynamic_mat_ass_couple::dynamic_mat_ass_couple
subroutine dynamic_mat_ass_couple(hecMESH, hecMAT, fstrSOLID, fstrCPL)
Definition: dynamic_mat_ass_couple.f90:17
m_fstr::fstr_solid
Definition: m_fstr.f90:238
m_fstr::idbg
integer(kind=kint), parameter idbg
Definition: m_fstr.f90:111
m_dynamic_mat_ass_couple
This module contains functions relates to coupling analysis.
Definition: dynamic_mat_ass_couple.f90:6
m_fstr
This module defines common data and basic structures for analysis.
Definition: m_fstr.f90:15
m_fstr::fstr_couple
Data for coupling analysis.
Definition: m_fstr.f90:611