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