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