FrontISTR  5.7.0
Large-scale structural analysis program with finit element method
dynamic_mat_ass_load.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 m_dynamic_mat_ass_load
8 contains
9 
10  !C
11  !C***
13  !C***
14  !C
15  subroutine dynamic_mat_ass_load(hecMESH, hecMAT, fstrSOLID, fstrDYNAMIC, fstrPARAM, iter )
16 
17  use m_fstr
18  use m_static_lib
19  use m_fstr_precheck
20  use m_table_dyn
21  use m_common_struct
22  use m_utilities
23 
24  implicit none
25  type(hecmwst_matrix) :: hecMAT
26  type(hecmwst_local_mesh) :: hecMESH
27  type(fstr_solid) :: fstrSOLID
28  type(fstr_dynamic) :: fstrDYNAMIC
29  type(fstr_param) :: fstrPARAM
30 
31  real(kind=kreal) :: xx(20), yy(20), zz(20)
32  real(kind=kreal) :: params(0:6)
33  real(kind=kreal) :: vect(60)
34  integer(kind=kint) :: iwk(60)
35  integer(kind=kint) :: nodLocal(20)
36  real(kind=kreal) :: tt(20), tt0(20), coords(3,3)
37  real(kind=kreal),pointer:: temp(:)
38  integer(kind=kint) :: ndof, ig0, ig, ityp, ltype, iS0, iE0, ik, in, i, j
39  integer(kind=kint) :: icel, ic_type, nn, is, isect, id, iset, nsize
40  integer(kind=kint) :: itype, iE, cdsys_ID
41  real(kind=kreal) :: val, rho, thick, pa1
42  logical :: fg_surf
43  logical, save :: isFirst = .true.
44 
45  integer(kind=kint) :: flag_u, ierror
46  integer(kind=kint), optional :: iter
47  real(kind=kreal) :: f_t
48 
49  integer(kind=kint) :: iiS, idofS, idofE
50  real(kind=kreal) :: ecoord(3, 20)
51  real(kind=kreal) :: v(6, 20), dv(6, 20), r(6*20)
52  real(kind=kreal) :: rhs
53  real(kind=kreal) :: unode_tmp(hecmat%NDOF*hecmesh%n_node)
54 
55  !for torque load
56  integer(kind=kint) :: n_rot, rid, n_nodes, idof
57  type(trotinfo) :: rinfo
58  real(kind=kreal) :: tval, normal(3), direc(3), ccoord(3), cdisp(3), cdiff(3)
59 
60  ndof = hecmat%NDOF
61  call hecmw_mat_clear_b( hecmat )
62  !C
63  !C CLOAD
64  !C
65  n_rot = fstrsolid%CLOAD_ngrp_rot
66  if( n_rot > 0 ) call fstr_rotinfo_init(n_rot, rinfo)
67 
68  do ig0 = 1, fstrsolid%CLOAD_ngrp_tot
69  ig = fstrsolid%CLOAD_ngrp_ID(ig0)
70  ityp = fstrsolid%CLOAD_ngrp_DOF(ig0)
71  val = fstrsolid%CLOAD_ngrp_val(ig0)
72 
73  flag_u = 0
74  call table_dyn(hecmesh, fstrsolid, fstrdynamic, ig0, f_t, flag_u)
75  val = val*f_t
76 
77  is0= hecmesh%node_group%grp_index(ig-1)+1
78  ie0= hecmesh%node_group%grp_index(ig )
79 
80  if( fstrsolid%CLOAD_ngrp_rotID(ig0) > 0 ) then ! setup torque load information
81  rid = fstrsolid%CLOAD_ngrp_rotID(ig0)
82  if( .not. rinfo%conds(rid)%active ) then
83  rinfo%conds(rid)%active = .true.
84  rinfo%conds(rid)%center_ngrp_id = fstrsolid%CLOAD_ngrp_centerID(ig0)
85  rinfo%conds(rid)%torque_ngrp_id = ig
86  endif
87  if( ityp>ndof ) ityp = ityp-ndof
88  rinfo%conds(rid)%vec(ityp) = val
89  cycle
90  endif
91 
92  do ik = is0, ie0
93  in = hecmesh%node_group%grp_item(ik)
94  hecmat%B( ndof*(in-1)+ityp ) = hecmat%B( ndof*(in-1)+ityp )+val
95  enddo
96  enddo
97 
98  !Add torque load to hecMAT%B
99  do rid = 1, n_rot
100  if( .not. rinfo%conds(rid)%active ) cycle
101  !get number of slave nodes
102  n_nodes = hecmw_ngrp_get_number(hecmesh, rinfo%conds(rid)%torque_ngrp_id)
103 
104  !get center node
105  ig = rinfo%conds(rid)%center_ngrp_id
106  do idof = 1, ndof
107  ccoord(idof) = hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, hecmesh%node)
108  cdisp(idof) = hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, fstrsolid%unode)
109  enddo
110  ccoord(1:ndof) = ccoord(1:ndof) + cdisp(1:ndof)
111 
112  tval = dsqrt(dot_product(rinfo%conds(rid)%vec(1:ndof),rinfo%conds(rid)%vec(1:ndof)))
113  if( tval < 1.d-16 ) then
114  write(*,*) '###ERROR### : norm of torque vector must be > 0.0'
115  call hecmw_abort( hecmw_comm_get_comm() )
116  endif
117  normal(1:ndof) = rinfo%conds(rid)%vec(1:ndof)/tval
118  tval = tval/dble(n_nodes)
119 
120  ig = rinfo%conds(rid)%torque_ngrp_id
121  is0 = hecmesh%node_group%grp_index(ig-1) + 1
122  ie0 = hecmesh%node_group%grp_index(ig )
123  do ik = is0, ie0
124  in = hecmesh%node_group%grp_item(ik)
125  cdiff(1:ndof) = hecmesh%node(ndof*(in-1)+1:ndof*in)+fstrsolid%unode(ndof*(in-1)+1:ndof*in)-ccoord(1:ndof)
126  call cross_product(normal,cdiff,vect(1:ndof))
127  val = dot_product(vect(1:ndof),vect(1:ndof))
128  if( val < 1.d-16 ) then
129  write(*,*) '###ERROR### : torque node is at the same position as that of center node in rotational surface.'
130  call hecmw_abort( hecmw_comm_get_comm() )
131  endif
132  vect(1:ndof) = (tval/val)*vect(1:ndof)
133  hecmat%B(ndof*(in-1)+1:ndof*in) = hecmat%B(ndof*(in-1)+1:ndof*in)+vect(1:ndof)
134  enddo
135  enddo
136  if( n_rot > 0 ) call fstr_rotinfo_finalize(rinfo)
137 
138  !C
139  !C DLOAD
140  !C
141  do ig0 = 1, fstrsolid%DLOAD_ngrp_tot
142  ig = fstrsolid%DLOAD_ngrp_ID(ig0)
143  ltype = fstrsolid%DLOAD_ngrp_LID(ig0)
144  do i = 0, 6
145  params(i) = fstrsolid%DLOAD_ngrp_params(i,ig0)
146  enddo
147  !C START & END
148  fg_surf = (ltype == 100)
149  if( fg_surf ) then ! surface group
150  is0 = hecmesh%surf_group%grp_index(ig-1) + 1
151  ie0 = hecmesh%surf_group%grp_index(ig )
152  else ! element group
153  is0 = hecmesh%elem_group%grp_index(ig-1) + 1
154  ie0 = hecmesh%elem_group%grp_index(ig )
155  endif
156 
157  do ik = is0, ie0
158  if( fg_surf ) then ! surface group
159  ltype = hecmesh%surf_group%grp_item(2*ik)*10
160  icel = hecmesh%surf_group%grp_item(2*ik-1)
161  ic_type = hecmesh%elem_type(icel)
162  else ! element group
163  icel = hecmesh%elem_group%grp_item(ik)
164  ic_type = hecmesh%elem_type(icel)
165  endif
166  !C** Create local stiffness
167  nn = hecmw_get_max_node(ic_type)
168  !C** node ID
169  is = hecmesh%elem_node_index(icel-1)
170  do j = 1, nn
171  nodlocal(j) = hecmesh%elem_node_item (is+j)
172  !C** nodal coordinate
173  xx(j) = hecmesh%node( 3*nodlocal(j)-2 )
174  yy(j) = hecmesh%node( 3*nodlocal(j)-1 )
175  zz(j) = hecmesh%node( 3*nodlocal(j) )
176  !C** create iwk array ***
177  do i = 1, ndof
178  iwk(ndof*(j-1)+i) = ndof*(nodlocal(j)-1)+i
179  enddo
180  enddo
181  !C** section ID
182  isect = hecmesh%section_ID(icel)
183  !C** Get Properties
184  rho = fstrsolid%elements(icel)%gausses(1)%pMaterial%variables(m_density)
185  call fstr_get_thickness(hecmesh,isect,thick)
186 
187  !C** Section Data
188  if( ndof == 2 ) then
189  id = hecmesh%section%sect_opt(isect)
190  if( id == 0 ) then
191  iset = 1
192  elseif( id == 1 ) then
193  iset = 0
194  elseif( id == 2 ) then
195  iset = 2
196  endif
197  pa1 = 1.0
198  endif
199 
200  !C** Create local stiffness
201  if( ic_type == 241 .or.ic_type == 242 .or. ic_type == 231 .or. ic_type == 232 ) then
202  call dl_c2(ic_type,nn,xx(1:nn),yy(1:nn),rho,pa1,ltype,params,vect(1:nn*ndof),nsize,iset)
203 
204  elseif( ic_type == 341 .or. ic_type == 351 .or. ic_type == 361 .or. &
205  ic_type == 342 .or. ic_type == 352 .or. ic_type == 362 ) then
206  call dl_c3(ic_type,nn,xx(1:nn),yy(1:nn),zz(1:nn),rho,ltype,params,vect(1:nn*ndof),nsize)
207 
208  elseif( ( ic_type == 741 ) .or. ( ic_type == 743 ) .or. ( ic_type == 731 ) ) then
209  call dl_shell(ic_type, nn, ndof, xx, yy, zz, rho, thick, ltype, params, vect, nsize, fstrsolid%elements(icel)%gausses)
210  elseif( ( ic_type==761 ) ) then
211  call dl_shell_33(ic_type, nn, ndof, xx, yy, zz, rho, thick, ltype, params, vect, nsize, &
212  fstrsolid%elements(icel)%gausses)
213  elseif( ( ic_type==781 ) ) then
214  call dl_shell_33(ic_type, nn, ndof, xx, yy, zz, rho, thick, ltype, params, vect, nsize, &
215  fstrsolid%elements(icel)%gausses)
216 
217  endif
218  !
219  !!!!!! time history
220 
221  flag_u = 10
222  call table_dyn(hecmesh, fstrsolid, fstrdynamic, ig0, f_t, flag_u)
223  do j=1,nsize
224  vect(j) = vect(j)*f_t
225  enddo
226  !
227  !C** Add vector
228  do j = 1, nsize
229  hecmat%B( iwk(j) )=hecmat%B( iwk(j) )+vect(j)
230  enddo
231  enddo
232  enddo
233 
234  if ( present(iter) ) then
235  if( iter == 1 ) then
236  do i = 1, ndof*hecmesh%n_node
237  unode_tmp(i) = fstrsolid%unode(i)
238  enddo
239 
240  do ig0 = 1, fstrsolid%BOUNDARY_ngrp_tot
241  ig = fstrsolid%BOUNDARY_ngrp_ID(ig0)
242  rhs = fstrsolid%BOUNDARY_ngrp_val(ig0)
243  ityp = fstrsolid%BOUNDARY_ngrp_type(ig0)
244  idofs = ityp/10
245  idofe = ityp-idofs*10
246  is0 = hecmesh%node_group%grp_index(ig-1) + 1
247  ie0 = hecmesh%node_group%grp_index(ig )
248 
249  do ik = is0, ie0
250  in = hecmesh%node_group%grp_item(ik)
251  do idof = idofs, idofe
252  unode_tmp( ndof*(in-1)+idof ) = rhs
253  enddo
254  enddo
255  enddo
256 
257  do itype = 1, hecmesh%n_elem_type
258  ic_type = hecmesh%elem_type_item(itype)
259  if( ic_type == 3414 ) then
260  nn = hecmw_get_max_node(ic_type)
261  if( nn > 20 ) stop "The number of elemental nodes > 20"
262 
263  is = hecmesh%elem_type_index(itype-1)+1
264  ie = hecmesh%elem_type_index(itype )
265  do icel = is, ie
266  if(fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype /= incomp_newtonian) then
267  write(*, *) '###ERROR### : This element is not supported for this material'
268  write(*, *) 'ic_type = ', ic_type, ', mtype = ', fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype
269  stop
270  call hecmw_abort(hecmw_comm_get_comm())
271  endif
272 
273  v = 0.0d0
274  dv= 0.0d0
275  iis = hecmesh%elem_node_index(icel-1)
276  do j = 1, nn
277  nodlocal(j) = hecmesh%elem_node_item(iis+j)
278  do i = 1, 3
279  ! nodal coordinates
280  ecoord(i,j) = hecmesh%node( 3*nodlocal(j)+i-3 )
281  ! nodal velocity
282  v(i,j) = unode_tmp( ndof*nodlocal(j)+i-ndof )
283  fstrsolid%unode( ndof*nodlocal(j)+i-ndof ) = v(i,j)
284  ! nodal velocity increment
285  dv(i,j) = fstrsolid%dunode( ndof*nodlocal(j)+i-ndof )
286  enddo
287  enddo
288 
289  call load_c3_vp &
290  ( ic_type, nn, ecoord(:,1:nn), v(1:4,1:nn), dv(1:4,1:nn), &
291  r(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), &
292  fstrdynamic%t_delta )
293 
294  do j = 1, nn
295  do i = 1, ndof
296  hecmat%B(ndof*(nodlocal(j)-1)+i) = hecmat%B(ndof*(nodlocal(j)-1)+i)+r(ndof*(j-1)+i)
297  enddo
298  enddo
299  enddo ! icel
300  endif
301  enddo ! itype
302 
303  else
304  do itype = 1, hecmesh%n_elem_type
305  ic_type = hecmesh%elem_type_item(itype)
306  if( ic_type == 3414 ) then
307  nn = hecmw_get_max_node(ic_type)
308  if( nn > 20 ) stop "The number of elemental nodes > 20"
309  is = hecmesh%elem_type_index(itype-1)+1
310  ie = hecmesh%elem_type_index(itype )
311  do icel = is, ie
312  iis = hecmesh%elem_node_index(icel-1)
313  do j = 1, nn
314  nodlocal(j) = hecmesh%elem_node_item(iis+j)
315  enddo
316  do j = 1, nn
317  do i = 1, ndof
318  hecmat%B(ndof*(nodlocal(j)-1)+i) = 0.0d0
319  enddo
320  enddo
321  enddo
322  endif
323  enddo
324  endif
325  endif
326 
327  !C
328  !C THERMAL LOAD USING TEMPERATURE
329  !C
330  !C Set Temperature
331  !C
332  if( fstrsolid%TEMP_ngrp_tot > 0 ) then
333 
334  if( hecmesh%my_rank .eq. 0 ) then
335  write(imsg,*) 'stop: THERMAL LOAD is not yet available in dynamic analysis!'
336  endif
337  call hecmw_abort( hecmw_comm_get_comm())
338 
339  allocate ( temp(hecmesh%n_node) )
340  temp=0
341  do ig0= 1, fstrsolid%TEMP_ngrp_tot
342  ig= fstrsolid%TEMP_ngrp_ID(ig0)
343  val=fstrsolid%TEMP_ngrp_val(ig0)
344  !C START & END
345  is0= hecmesh%node_group%grp_index(ig-1) + 1
346  ie0= hecmesh%node_group%grp_index(ig )
347  do ik= is0, ie0
348  in = hecmesh%node_group%grp_item(ik)
349  temp( in ) = val
350  enddo
351  enddo
352  !C
353  !C +-------------------------------+
354  !C | ELEMENT-by-ELEMENT ASSEMBLING |
355  !C | according to ELEMENT TYPE |
356  !C +-------------------------------+
357  !C===
358  do itype = 1, hecmesh%n_elem_type
359 
360  is = hecmesh%elem_type_index(itype-1)+1
361  ie = hecmesh%elem_type_index(itype )
362  ic_type = hecmesh%elem_type_item(itype)
363  if( hecmw_is_etype_link(ic_type) ) cycle
364  if( hecmw_is_etype_patch(ic_type) ) cycle
365  !C** Set number of nodes
366  nn = hecmw_get_max_node(ic_type)
367  !C element loop
368  do icel = is, ie
369  !C** node ID
370  is= hecmesh%elem_node_index(icel-1)
371  do j=1,nn
372  nodlocal(j)=hecmesh%elem_node_item(is+j)
373  !C** nodal coordinate
374  xx(j)=hecmesh%node(3*nodlocal(j)-2)
375  yy(j)=hecmesh%node(3*nodlocal(j)-1)
376  zz(j)=hecmesh%node(3*nodlocal(j) )
377  tt(j)=temp( nodlocal(j) )
378  tt0(j)=ref_temp
379  !C** create iwk array ***
380  do i=1,ndof
381  iwk(ndof*(j-1)+i)=ndof*(nodlocal(j)-1)+i
382  enddo
383  enddo
384 
385  !C** section ID
386  isect= hecmesh%section_ID(icel)
387  cdsys_id = fstrsolid%elements(icel)%gausses(1)%pMaterial%cdsys_ID
388  call get_coordsys( cdsys_id, hecmesh, fstrsolid, coords )
389 
390  !C** Section Data
391  if( ndof .eq. 2 ) then
392  id=hecmesh%section%sect_opt(isect)
393  if( id.eq.0 ) then
394  iset=1
395  elseif( id.eq.1) then
396  iset=0
397  elseif( id.eq.2) then
398  iset=2
399  endif
400  pa1=1.0
401  endif
402 
403  !C** Create local stiffness
404  if( ic_type == 241 .or. ic_type == 242 .or. ic_type == 231 .or. ic_type == 232 ) then
405  call tload_c2( ic_type, nn, xx(1:nn), yy(1:nn), tt(1:nn), tt0(1:nn), &
406  fstrsolid%elements(icel)%gausses,pa1,iset, vect(1:nn*2) )
407 
408  elseif( ic_type == 361 ) then
409  if( fstrsolid%sections(isect)%elemopt361 == kel361fi ) then
410  call tload_c3 &
411  ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
412  fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords )
413  else if( fstrsolid%sections(isect)%elemopt361 == kel361bbar ) then
414  call tload_c3d8bbar &
415  ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
416  fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords )
417  else if( fstrsolid%sections(isect)%elemopt361 == kel361ic ) then
418  call tload_c3d8ic &
419  ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
420  fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords )
421  else if( fstrsolid%sections(isect)%elemopt361 == kel361fbar ) then
422  call tload_c3d8fbar &
423  ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
424  fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords )
425  endif
426 
427  elseif (ic_type == 341 .or. ic_type == 351 .or. &
428  ic_type == 342 .or. ic_type == 352 .or. ic_type == 362 ) then
429  call tload_c3 &
430  ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
431  fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords )
432 
433  elseif ( ic_type == 741 .or. ic_type == 743 .or. ic_type == 731 ) then
434  if( myrank == 0 ) then
435  write(imsg,*) '*------------------------', &
436  '-------------------*'
437  write(imsg,*) ' Thermal loading option ', &
438  'not yet available.'
439  write(imsg,*) '*------------------------', &
440  '-------------------*'
441  call hecmw_abort( hecmw_comm_get_comm())
442  endif
443  endif
444  !C** Add vector
445  do j = 1, ndof*nn
446  hecmat%B( iwk(j) ) = hecmat%B( iwk(j) )+vect(j)
447  enddo
448  enddo
449  enddo
450  deallocate ( temp )
451  endif
452  end subroutine dynamic_mat_ass_load
453 
454 
455 end module m_dynamic_mat_ass_load
m_table_dyn
Table of lading step in dynamic analysis.
Definition: table_dyn.f90:6
m_utilities
This module provides aux functions.
Definition: utilities.f90:6
m_fstr::kel361fi
integer(kind=kint), parameter kel361fi
Definition: m_fstr.f90:77
m_fstr::myrank
integer(kind=kint) myrank
PARALLEL EXECUTION.
Definition: m_fstr.f90:96
m_fstr::kel361bbar
integer(kind=kint), parameter kel361bbar
Definition: m_fstr.f90:78
m_common_struct::fstr_rotinfo_init
subroutine fstr_rotinfo_init(n, rinfo)
Definition: m_common_struct.f90:99
m_fstr::fstr_solid
Definition: m_fstr.f90:238
m_fstr::fstr_dynamic
Data for DYNAMIC ANSLYSIS (fstrDYNAMIC)
Definition: m_fstr.f90:504
m_fstr::fstr_param
FSTR INNER CONTROL PARAMETERS (fstrPARAM)
Definition: m_fstr.f90:154
m_fstr_precheck::fstr_get_thickness
subroutine fstr_get_thickness(hecMESH, mid, thick)
Definition: fstr_precheck.f90:50
m_common_struct
This modules defines common structures for fem analysis.
Definition: m_common_struct.f90:6
m_dynamic_mat_ass_load::dynamic_mat_ass_load
subroutine dynamic_mat_ass_load(hecMESH, hecMAT, fstrSOLID, fstrDYNAMIC, fstrPARAM, iter)
This function sets boundary condition of external load.
Definition: dynamic_mat_ass_load.f90:16
m_fstr
This module defines common data and basic structures for analysis.
Definition: m_fstr.f90:15
m_fstr::ref_temp
real(kind=kreal), pointer ref_temp
REFTEMP.
Definition: m_fstr.f90:136
m_dynamic_mat_ass_load
This module contains function to set boundary condition of external load in dynamic analysis.
Definition: dynamic_mat_ass_load.f90:7
m_common_struct::fstr_rotinfo_finalize
subroutine fstr_rotinfo_finalize(rinfo)
Definition: m_common_struct.f90:117
m_fstr::get_coordsys
subroutine get_coordsys(cdsys_ID, hecMESH, fstrSOLID, coords)
This subroutine fetch coords defined by local coordinate system.
Definition: m_fstr.f90:1073
m_fstr::kel361fbar
integer(kind=kint), parameter kel361fbar
Definition: m_fstr.f90:80
m_fstr::kel361ic
integer(kind=kint), parameter kel361ic
Definition: m_fstr.f90:79
m_fstr_precheck
This module provides function to check input data of IFSTR solver.
Definition: fstr_precheck.f90:6
m_table_dyn::table_dyn
subroutine table_dyn(hecMESH, fstrSOLID, fstrDYNAMIC, ig0, f_t, flag_u)
Definition: table_dyn.f90:19
m_static_lib
This modules just summarizes all modules used in static analysis.
Definition: static_LIB.f90:6
m_common_struct::trotinfo
Definition: m_common_struct.f90:21
m_utilities::cross_product
subroutine cross_product(v1, v2, vn)
Definition: utilities.f90:330
m_fstr::imsg
integer(kind=kint), parameter imsg
Definition: m_fstr.f90:110