FrontISTR  5.9.0
Large-scale structural analysis program with finit element method
fstr_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_fstr_ass_load
8  use m_fstr
9  use m_static_lib
10  use m_fstr_precheck
11  use m_fstr_elemact
12  use mmechgauss
13  use mreadtemp
14  use muload
15  use m_fstr_spring
16  use m_common_struct
17  use m_utilities
18 
19  implicit none
20 
21 contains
22  !
23  !======================================================================!
29 
30  subroutine fstr_ass_load(cstep, ctime, hecMESH, hecMAT, fstrSOLID, fstrPARAM)
31  !======================================================================!
32  integer(kind=kint), intent(in) :: cstep
33  real(kind=kreal), intent(in) :: ctime
34  type(hecmwst_matrix), intent(inout) :: hecmat
35  type(hecmwst_local_mesh), intent(in) :: hecMESH
36  type(fstr_solid), intent(inout) :: fstrSOLID
37  type(fstr_param), intent(inout) :: fstrPARAM
38 
39  ! Initialize the global load vector
40  fstrsolid%GL(:) = 0.0d0
41  fstrsolid%EFORCE(:) = 0.0d0
42 
43  ! Process concentrated nodal forces (CLOAD)
44  call process_concentrated_loads(cstep, ctime, hecmesh, fstrsolid)
45 
46  ! Process distributed loads (DLOAD) - surface pressure and volume force
47  call process_distributed_loads(cstep, ctime, hecmesh, fstrsolid)
48 
49  ! Process user-defined loads
50  call process_user_loads(cstep, fstrsolid)
51 
52  ! Update global load vector
53  call hecmw_update_r(hecmesh, fstrsolid%GL, hecmesh%n_node, hecmesh%n_dof)
54 
55  ! Update right-hand side vector
56  call hecmw_mat_clear_b(hecmat)
57  call update_rhs_vector(hecmesh, hecmat, fstrsolid)
58 
59  ! Process thermal loads (TLOAD)
60  call process_thermal_loads(cstep, ctime, hecmesh, hecmat, fstrsolid)
61 
62  ! Process spring forces
63  call fstr_update_ndforce_spring(cstep, hecmesh, fstrsolid, hecmat%B)
64 
65  end subroutine fstr_ass_load
66 
67  !======================================================================!
69  !======================================================================!
70  subroutine process_concentrated_loads(cstep, ctime, hecMESH, fstrSOLID)
71  integer(kind=kint), intent(in) :: cstep
72  real(kind=kreal), intent(in) :: ctime
73  type(hecmwst_local_mesh), intent(in) :: hecmesh
74  type(fstr_solid), intent(inout) :: fstrSOLID
75 
76  integer(kind=kint) :: n_rot, rid, n_nodes, idof, ndof
77  integer(kind=kint) :: ig0, ig, ityp, iS0, iE0, ik, in, grpid, jj_n_amp
78  real(kind=kreal) :: aval, fval, tval
79  real(kind=kreal) :: normal(3), direc(3), ccoord(3), cdisp(3), cdiff(3)
80  real(kind=kreal) :: vect(60)
81  type(trotinfo) :: rinfo
82 
83  ndof = hecmesh%n_dof
84 
85  ! Initialize rotation information for torque loads
86  n_rot = fstrsolid%CLOAD_ngrp_rot
87  if (n_rot > 0) call fstr_rotinfo_init(n_rot, rinfo)
88 
89  ! Process all concentrated loads
90  do ig0 = 1, fstrsolid%CLOAD_ngrp_tot
91  grpid = fstrsolid%CLOAD_ngrp_GRPID(ig0)
92  if (.not. fstr_isloadactive(fstrsolid, grpid, cstep)) cycle
93  ig = fstrsolid%CLOAD_ngrp_ID(ig0)
94  ityp = fstrsolid%CLOAD_ngrp_DOF(ig0)
95  fval = fstrsolid%CLOAD_ngrp_val(ig0)
96  jj_n_amp = fstrsolid%CLOAD_ngrp_amp(ig0)
97  aval = fval
98  if (jj_n_amp <= 0) then ! Amplitude not defined
99  aval = fstrsolid%FACTOR(2) * fval
100  else
101  call fstr_get_amplitude(hecmesh, fstrsolid, cstep, jj_n_amp, ctime, aval)
102  endif
103 
104  ! A load continuing from a previous step is held at its nominal value,
105  ! BUT only when no amplitude is given. An explicit AMP= prescribes the
106  ! load's time history, so it is honored in every active step (consistent
107  ! with the dynamic load path, which has no such override).
108  if (jj_n_amp <= 0 .and. fstr_isloadactive(fstrsolid, grpid, cstep-1)) aval = fval
109  is0 = hecmesh%node_group%grp_index(ig-1) + 1
110  ie0 = hecmesh%node_group%grp_index(ig)
111 
112  if( fstrsolid%CLOAD_ngrp_rotID(ig0) > 0 ) then ! setup torque load information
113  rid = fstrsolid%CLOAD_ngrp_rotID(ig0)
114  if (.not. rinfo%conds(rid)%active) then
115  rinfo%conds(rid)%active = .true.
116  rinfo%conds(rid)%center_ngrp_id = fstrsolid%CLOAD_ngrp_centerID(ig0)
117  rinfo%conds(rid)%torque_ngrp_id = ig
118  endif
119  if (ityp > ndof) ityp = ityp - ndof
120  rinfo%conds(rid)%vec(ityp) = aval
121  cycle
122  endif
123 
124  do ik = is0, ie0
125  in = hecmesh%node_group%grp_item(ik)
126  fstrsolid%GL(ndof*(in-1)+ityp) = fstrsolid%GL(ndof*(in-1)+ityp) + aval
127  enddo
128  enddo
129 
130  !Add torque load to fstrSOLID%GL
131  do rid = 1, n_rot
132  if (.not. rinfo%conds(rid)%active) cycle
133  ! Get number of slave nodes
134  n_nodes = hecmw_ngrp_get_number(hecmesh, rinfo%conds(rid)%torque_ngrp_id)
135 
136  ! Get center node
137  ig = rinfo%conds(rid)%center_ngrp_id
138  do idof = 1, ndof
139  ccoord(idof) = hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, hecmesh%node)
140  cdisp(idof) = hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, fstrsolid%unode)
141  cdisp(idof) = cdisp(idof) + hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, fstrsolid%dunode)
142  enddo
143  ccoord(1:ndof) = ccoord(1:ndof) + cdisp(1:ndof)
144 
145  tval = dsqrt(dot_product(rinfo%conds(rid)%vec(1:ndof), rinfo%conds(rid)%vec(1:ndof)))
146  if (tval < 1.d-16) then
147  write(*,*) '###ERROR### : norm of torque vector must be > 0.0'
148  call hecmw_abort(hecmw_comm_get_comm())
149  endif
150  normal(1:ndof) = rinfo%conds(rid)%vec(1:ndof) / tval
151  tval = tval / dble(n_nodes)
152 
153  ig = rinfo%conds(rid)%torque_ngrp_id
154  is0 = hecmesh%node_group%grp_index(ig-1) + 1
155  ie0 = hecmesh%node_group%grp_index(ig)
156  do ik = is0, ie0
157  in = hecmesh%node_group%grp_item(ik)
158  cdiff(1:ndof) = hecmesh%node(ndof*(in-1)+1:ndof*in) + fstrsolid%unode(ndof*(in-1)+1:ndof*in) &
159  & + fstrsolid%dunode(ndof*(in-1)+1:ndof*in) - ccoord(1:ndof)
160  call cross_product(normal,cdiff,vect(1:ndof))
161  fval = dot_product(vect(1:ndof), vect(1:ndof))
162  if (fval < 1.d-16) then
163  write(*,*) '###ERROR### : torque node is at the same position as that of center node in rotational surface.'
164  call hecmw_abort(hecmw_comm_get_comm())
165  endif
166  vect(1:ndof) = (tval/fval) * vect(1:ndof)
167  fstrsolid%GL(ndof*(in-1)+1:ndof*in) = fstrsolid%GL(ndof*(in-1)+1:ndof*in) + vect(1:ndof)
168  enddo
169  enddo
170  if (n_rot > 0) call fstr_rotinfo_finalize(rinfo)
171 
172  end subroutine process_concentrated_loads
173  !
174  ! -------------------------------------------------------------------
175  ! DLOAD
176  ! -------------------------------------------------------------------
177  subroutine process_distributed_loads(cstep, ctime, hecMESH, fstrSOLID)
178  integer(kind=kint), intent(in) :: cstep
179  real(kind=kreal), intent(in) :: ctime
180  type(hecmwst_local_mesh), intent(in) :: hecmesh
181  type(fstr_solid), intent(inout) :: fstrSOLID
182 
183  integer(kind=kint) :: ndof, ig0, ig, ltype, iS0, iE0, ik, icel, ic_type, nn, is
184  integer(kind=kint) :: isect, id, iset, ihead, nsize, grpid, i, j, jj_n_amp
185  integer(kind=kint) :: iwk(60), nodLocal(20)
186  real(kind=kreal) :: xx(20), yy(20), zz(20), vect(60), params(0:6)
187  real(kind=kreal) :: factor, aval, rho, thick, pa1
188  logical :: fg_surf
189  type(tmaterial), pointer :: material
190 
191  ndof = hecmesh%n_dof
192 
193  do ig0 = 1, fstrsolid%DLOAD_ngrp_tot
194  grpid = fstrsolid%DLOAD_ngrp_GRPID(ig0)
195  if (.not. fstr_isloadactive(fstrsolid, grpid, cstep)) cycle
196  jj_n_amp = fstrsolid%DLOAD_ngrp_amp(ig0)
197  ig = fstrsolid%DLOAD_ngrp_ID(ig0)
198  ltype = fstrsolid%DLOAD_ngrp_LID(ig0)
199  do i = 0, 6
200  params(i) = fstrsolid%DLOAD_ngrp_params(i, ig0)
201  enddo
202  if (jj_n_amp <= 0) then ! Amplitude not defined
203  factor = fstrsolid%factor(2)
204  else
205  ! params(0) is the load magnitude (pressure / body force value,
206  ! gravity acceleration, or angular velocity). Passing it to the
207  ! evaluator lets VALUE=ABSOLUTE replace the magnitude with the
208  ! amplitude value, while VALUE=RELATIVE reduces the factor to a(t).
209  aval = params(0)
210  call fstr_get_amplitude(hecmesh, fstrsolid, cstep, jj_n_amp, ctime, aval)
211  if (dabs(params(0)) > 1.d-30) then
212  factor = aval / params(0)
213  else
214  factor = 1.0d0
215  endif
216  endif
217 
218  ! Continuing load held at nominal only when no amplitude is given;
219  ! an explicit AMP= is honored in every active step (see CLOAD above).
220  if (jj_n_amp <= 0 .and. fstr_isloadactive(fstrsolid, grpid, cstep-1)) factor = 1.0d0
221  ! ----- START & END
222  fg_surf = (ltype == 100)
223  if( fg_surf ) then ! surface group
224  is0 = hecmesh%surf_group%grp_index(ig-1) + 1
225  ie0 = hecmesh%surf_group%grp_index(ig)
226  else ! element group
227  is0 = hecmesh%elem_group%grp_index(ig-1) + 1
228  ie0 = hecmesh%elem_group%grp_index(ig)
229  endif
230  do ik = is0, ie0
231  if( fg_surf ) then ! surface group
232  ltype = hecmesh%surf_group%grp_item(2*ik) * 10
233  icel = hecmesh%surf_group%grp_item(2*ik-1)
234  ic_type = hecmesh%elem_type(icel)
235  else ! element group
236  icel = hecmesh%elem_group%grp_item(ik)
237  ic_type = hecmesh%elem_type(icel)
238  endif
239 
240  !ELEMENT ACTIVATION
241  if( fstrsolid%elements(icel)%elemact_flag == kelact_inactive ) cycle
242 
243  if (hecmw_is_etype_link(ic_type)) cycle
244  if (hecmw_is_etype_patch(ic_type)) cycle
245  ! if( ic_type==3422 ) ic_type=342
246  nn = hecmw_get_max_node(ic_type)
247  ! ----- node ID
248  is = hecmesh%elem_node_index(icel-1)
249  if (fstrsolid%DLOAD_follow == 0) then
250  do j = 1, nn
251  nodlocal(j) = hecmesh%elem_node_item (is+j)
252  ! ----- nodal coordinate
253  xx(j) = hecmesh%node( 3*nodlocal(j)-2 )
254  yy(j) = hecmesh%node( 3*nodlocal(j)-1 )
255  zz(j) = hecmesh%node( 3*nodlocal(j) )
256  ! ----- create iwk array ***
257  do i = 1, ndof
258  iwk( ndof*(j-1)+i ) = ndof*( nodlocal(j)-1 )+i
259  enddo
260  enddo
261  else
262  do j = 1, nn
263  nodlocal(j) = hecmesh%elem_node_item (is+j)
264  ! ----- nodal coordinate
265  if (ndof==2) then
266  xx(j) = hecmesh%node( 3*nodlocal(j)-2 )+fstrsolid%unode( 2*nodlocal(j)-1 )+fstrsolid%dunode( 2*nodlocal(j)-1 )
267  yy(j) = hecmesh%node( 3*nodlocal(j)-1 )+fstrsolid%unode( 2*nodlocal(j) )+fstrsolid%dunode( 2*nodlocal(j) )
268  else if (ndof==3) then
269  xx(j) = hecmesh%node( 3*nodlocal(j)-2 )+fstrsolid%unode( 3*nodlocal(j)-2 )+fstrsolid%dunode( 3*nodlocal(j)-2 )
270  yy(j) = hecmesh%node( 3*nodlocal(j)-1 )+fstrsolid%unode( 3*nodlocal(j)-1 )+fstrsolid%dunode( 3*nodlocal(j)-1 )
271  zz(j) = hecmesh%node( 3*nodlocal(j) )+fstrsolid%unode( 3*nodlocal(j) )+fstrsolid%dunode( 3*nodlocal(j) )
272  else if (ndof==6) then
273  xx(j) = hecmesh%node( 3*nodlocal(j)-2 )+fstrsolid%unode( 6*nodlocal(j)-5 )+fstrsolid%dunode( 6*nodlocal(j)-5 )
274  yy(j) = hecmesh%node( 3*nodlocal(j)-1 )+fstrsolid%unode( 6*nodlocal(j)-4 )+fstrsolid%dunode( 6*nodlocal(j)-4 )
275  zz(j) = hecmesh%node( 3*nodlocal(j) )+fstrsolid%unode( 6*nodlocal(j)-3 )+fstrsolid%dunode( 6*nodlocal(j)-3 )
276  endif
277  ! ----- create iwk array ***
278  do i = 1, ndof
279  iwk( ndof*(j-1)+i ) = ndof*( nodlocal(j)-1 )+i
280  enddo
281  enddo
282  endif
283  ! ----- section ID
284  isect = hecmesh%section_ID(icel)
285  ! ----- Get Properties
286  material => fstrsolid%elements(icel)%gausses(1)%pMaterial
287  rho = material%variables(m_density)
288  call fstr_get_thickness(hecmesh, isect, thick)
289  ! ----- Section Data
290  if (ndof == 2) then
291  id = hecmesh%section%sect_opt(isect)
292  if (id == 0) then
293  iset = 1
294  else if (id == 1) then
295  iset = 0
296  else if (id == 2) then
297  iset = 2
298  endif
299  pa1 = 1.d0
300  endif
301  ! ----- Create local stiffness
302  if (ic_type==301)then
303  ihead = hecmesh%section%sect_R_index(isect-1)
304  call dl_c1(ic_type,nn,xx(1:nn),yy(1:nn),zz(1:nn),rho,thick,ltype,params,vect(1:nn*ndof),nsize)
305 
306  elseif( ic_type == 241 .or. ic_type == 242 .or. ic_type == 231 .or. ic_type == 232 .or. ic_type == 2322 ) then
307  call dl_c2(ic_type,nn,xx(1:nn),yy(1:nn),rho,pa1,ltype,params,vect(1:nn*ndof),nsize,iset)
308 
309  else if ( ic_type == 341 .or. ic_type == 351 .or. ic_type == 361 .or. &
310  ic_type == 342 .or. ic_type == 352 .or. ic_type == 362 ) then
311  call dl_c3(ic_type,nn,xx(1:nn),yy(1:nn),zz(1:nn),rho,ltype,params,vect(1:nn*ndof),nsize)
312 
313  else if ( ic_type == 641 ) then
314  ihead = hecmesh%section%sect_R_index(isect-1)
315  call dl_beam_641(ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), rho, ltype, params, &
316  hecmesh%section%sect_R_item(ihead+1:), vect(1:nn*ndof), nsize)
317 
318  else if( ( ic_type == 741 ) .or. ( ic_type == 743 ) .or. ( ic_type == 731 ) ) then
319  call dl_shell(ic_type, nn, ndof, xx, yy, zz, rho, thick, ltype, params, vect, nsize, fstrsolid%elements(icel)%gausses)
320 
321  else if( ( ic_type==761 ) .or. ( ic_type==781 ) ) then
322  call dl_shell_33(ic_type, nn, ndof, xx, yy, zz, rho, thick, ltype, params, vect, nsize, &
323  fstrsolid%elements(icel)%gausses)
324 
325  else
326  nsize = 0
327  write(*,*)"### WARNING: DLOAD",ic_type
328 
329  endif
330  ! ----- Add vector
331  do j = 1, nsize
332  fstrsolid%GL(iwk(j)) = fstrsolid%GL(iwk(j)) + factor * vect(j)
333  enddo
334  enddo
335  enddo
336 
337  end subroutine process_distributed_loads
338 
339  ! -----Uload
340  subroutine process_user_loads(cstep, fstrSOLID)
341  integer(kind=kint), intent(in) :: cstep
342  type(fstr_solid), intent(inout) :: fstrSOLID
343 
344  real(kind=kreal) :: factor
345 
346  factor = fstrsolid%factor(2)
347  call uloading(cstep, factor, fstrsolid%GL)
348 
349  end subroutine process_user_loads
350 
351  !======================================================================!
353  !======================================================================!
354  subroutine update_rhs_vector(hecMESH, hecMAT, fstrSOLID)
355  type(hecmwst_local_mesh), intent(in) :: hecmesh
356  type(hecmwst_matrix), intent(inout) :: hecMAT
357  type(fstr_solid), intent(inout) :: fstrSOLID
358 
359  integer(kind=kint) :: i
360 
361  do i = 1, hecmesh%n_node * hecmesh%n_dof
362  hecmat%B(i) = fstrsolid%GL(i) - fstrsolid%QFORCE(i)
363  enddo
364 
365  do i = 1, hecmat%NDOF * hecmat%NP
366  !thermal load is not considered
367  fstrsolid%EFORCE(i) = fstrsolid%GL(i)
368  enddo
369 
370  end subroutine update_rhs_vector
371 
372  ! -------------------------------------------------------------------
373  ! TLOAD : THERMAL LOAD USING TEMPERATURE
374  ! -------------------------------------------------------------------
375  subroutine process_thermal_loads(cstep, ctime, hecMESH, hecMAT, fstrSOLID)
376  integer(kind=kint), intent(in) :: cstep
377  real(kind=kreal), intent(in) :: ctime
378  type(hecmwst_local_mesh), intent(in) :: hecmesh
379  type(hecmwst_matrix), intent(inout) :: hecMAT
380  type(fstr_solid), intent(inout) :: fstrSOLID
381 
382  integer(kind=kint) :: ndof, ig0, ig, iS0, iE0, ik, in, grpid
383  integer(kind=kint) :: itype, is, iE, icel, ic_type, nn, isect, cdsys_ID, id, iset
384  integer(kind=kint) :: i, j, ihead, tstep, nodLocal(20), iwk(60)
385  real(kind=kreal) :: factor, fval, pa1
386  real(kind=kreal) :: xx(20), yy(20), zz(20), tt(20), tt0(20), coords(3,3), vect(60)
387  real(kind=kreal) :: local_coords(3,3) ! Local copy for coordinate transformation
388 
389  ndof = hecmesh%n_dof
390 
391  if (fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0) then
392  do ig0 = 1, fstrsolid%TEMP_ngrp_tot
393  grpid = fstrsolid%TEMP_ngrp_GRPID(ig0)
394  if (.not. fstr_isloadactive(fstrsolid, grpid, cstep)) cycle
395  factor = fstrsolid%factor(2)
396  if (fstr_isloadactive(fstrsolid, grpid, cstep-1)) factor = 1.0d0
397  ig = fstrsolid%TEMP_ngrp_ID(ig0)
398  fval = fstrsolid%TEMP_ngrp_val(ig0)
399  is0 = hecmesh%node_group%grp_index(ig-1) + 1
400  ie0 = hecmesh%node_group%grp_index(ig)
401  do ik = is0, ie0
402  in = hecmesh%node_group%grp_item(ik)
403  pa1 = fstrsolid%temp_bak(in)
404  fstrsolid%temperature(in) = pa1 + (fval - pa1) * factor
405  enddo
406  enddo
407 
408  if (fstrsolid%TEMP_irres > 0) then
409  call read_temperature_result(hecmesh, fstrsolid%TEMP_irres, fstrsolid%TEMP_tstep, &
410  & fstrsolid%TEMP_rtype, fstrsolid%TEMP_interval, fstrsolid%TEMP_factor, ctime, &
411  & fstrsolid%temperature, fstrsolid%temp_bak)
412  endif
413  endif
414 
415  ! ----- elemact element
416  if( fstrsolid%elemact%ELEMACT_egrp_tot > 0 ) &
417  & call fstr_update_elemact_solid( hecmesh, fstrsolid, cstep, ctime )
418 
419  if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0 ) then
420  ! ----- element TYPE loop.
421  do itype = 1, hecmesh%n_elem_type
422  is = hecmesh%elem_type_index(itype-1) + 1
423  ie = hecmesh%elem_type_index(itype)
424  ic_type = hecmesh%elem_type_item(itype)
425  if (hecmw_is_etype_link(ic_type)) cycle
426  if (hecmw_is_etype_patch(ic_type)) cycle
427  ! ----- Set number of nodes
428  nn = hecmw_get_max_node(ic_type)
429 
430  ! ----- element loop
431  do icel = is, ie
432 
433  !ELEMENT ACTIVATION
434  if( fstrsolid%elements(icel)%elemact_flag == kelact_inactive ) cycle
435 
436  ! ----- node ID
437  is = hecmesh%elem_node_index(icel-1)
438  do j = 1, nn
439  nodlocal(j) = hecmesh%elem_node_item(is+j)
440  ! ----- nodal coordinate
441  if (ndof == 2) then
442  xx(j) = hecmesh%node(3*nodlocal(j)-2) + fstrsolid%unode(ndof*nodlocal(j)-1)
443  yy(j) = hecmesh%node(3*nodlocal(j)-1) + fstrsolid%unode(ndof*nodlocal(j))
444  else if (ndof == 3) then
445  xx(j) = hecmesh%node(3*nodlocal(j)-2) + fstrsolid%unode(ndof*nodlocal(j)-2)
446  yy(j) = hecmesh%node(3*nodlocal(j)-1) + fstrsolid%unode(ndof*nodlocal(j)-1)
447  zz(j) = hecmesh%node(3*nodlocal(j)) + fstrsolid%unode(ndof*nodlocal(j))
448  endif
449  tt0(j) = fstrsolid%last_temp(nodlocal(j))
450  tt(j) = fstrsolid%temperature(nodlocal(j))
451  ! ----- create iwk array ***
452  do i = 1, ndof
453  iwk(ndof*(j-1)+i) = ndof*(nodlocal(j)-1)+i
454  enddo
455  enddo
456 
457  ! ----- section Data
458  isect = hecmesh%section_ID(icel)
459  cdsys_id = hecmesh%section%sect_orien_ID(isect)
460  call get_coordsys(cdsys_id, hecmesh, fstrsolid, coords, icel)
461 
462  if( ndof == 2 ) then
463  id=hecmesh%section%sect_opt(isect)
464  if( id==0 ) then
465  iset = 1
466  else if (id == 1) then
467  iset = 0
468  else if (id == 2) then
469  iset = 2
470  endif
471  pa1 = 1.0d0
472  endif
473 
474  if (ic_type == 641) then
475  isect = hecmesh%section_ID(icel)
476  ihead = hecmesh%section%sect_R_index(isect-1)
477 
478  call tload_beam_641( ic_type, nn, ndof, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
479  fstrsolid%elements(icel)%gausses, hecmesh%section%sect_R_item(ihead+1:), &
480  vect(1:nn*ndof))
481 
482  do j = 1, ndof*nn
483  hecmat%B(iwk(j)) = hecmat%B(iwk(j)) + vect(j)
484  enddo
485  cycle
486  endif
487 
488  ! Local copy of coordinate data
489  local_coords = coords
490 
491  ! Calculate thermal load based on element type
492  call calculate_thermal_load(ic_type, nn, xx, yy, zz, tt, tt0, isect, ndof, &
493  hecmesh, fstrsolid, icel, vect, cdsys_id, local_coords, &
494  iset, pa1, iwk, hecmat%B)
495  enddo
496  enddo
497  endif
498 
499  end subroutine process_thermal_loads
500 
501  !======================================================================!
503  !======================================================================!
504  subroutine calculate_thermal_load(ic_type, nn, xx, yy, zz, tt, tt0, isect, ndof, &
505  hecMESH, fstrSOLID, icel, vect, cdsys_ID, coords, &
506  iset, pa1, iwk, B)
507  integer(kind=kint), intent(in) :: ic_type, nn, isect, ndof, cdsys_ID, iset
508  real(kind=kreal), intent(in) :: xx(*), yy(*), zz(*), tt(*), tt0(*), pa1
509  real(kind=kreal), intent(inout) :: coords(3,3) ! Changed INTENT from IN to INOUT
510  type(hecmwst_local_mesh), intent(in) :: hecmesh
511  type(fstr_solid), intent(in) :: fstrSOLID
512  integer(kind=kint), intent(in) :: icel
513  real(kind=kreal), intent(out) :: vect(*)
514  integer(kind=kint), intent(in) :: iwk(*)
515  real(kind=kreal), intent(inout) :: b(*)
516 
517  integer(kind=kint) :: j, myrank
518 
519  myrank = 0
520 
521  ! 2D elements
522  if (ic_type == 241 .or. ic_type == 242 .or. ic_type == 231 .or. ic_type == 232) then
523  call tload_c2(ic_type, nn, xx(1:nn), yy(1:nn), tt(1:nn), tt0(1:nn), &
524  fstrsolid%elements(icel)%gausses, pa1, iset, vect(1:nn*2))
525 
526  else if (ic_type == 361) then
527  if (fstrsolid%sections(isect)%elemopt361 == kel361fi) then
528  call tload_c3 &
529  ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
530  fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords)
531  else if (fstrsolid%sections(isect)%elemopt361 == kel361bbar) then
532  call tload_c3d8bbar &
533  ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
534  fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords)
535  else if (fstrsolid%sections(isect)%elemopt361 == kel361ic) then
536  call tload_c3d8ic &
537  ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
538  fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords)
539  else if (fstrsolid%sections(isect)%elemopt361 == kel361fbar) then
540  call tload_c3d8fbar &
541  ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
542  fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords)
543  else if (fstrsolid%sections(isect)%elemopt361 == kel361up) then
544  ! UP element : reuse the standard C3 thermal load. Note the thermal
545  ! strain is volumetric and TLOAD_C3 uses the full elastic matrix, so
546  ! this is not strictly consistent with the u-p deviatoric/pressure
547  ! split (see report); provided for basic thermal loading support.
548  call tload_c3 &
549  ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
550  fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords )
551  endif
552 
553  else if (ic_type == 341 .or. ic_type == 351 .or. &
554  ic_type == 342 .or. ic_type == 352 .or. ic_type == 362) then
555  call tload_c3 &
556  ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
557  fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords)
558 
559  else if (ic_type == 741 .or. ic_type == 743 .or. ic_type == 731) then
560  if (myrank == 0) then
561  write(imsg,*) '*------------------------', &
562  '-------------------*'
563  write(imsg,*) ' Thermal loading option for shell elements', &
564  'not yet available.'
565  write(imsg,*) '*------------------------', &
566  '-------------------*'
567  call hecmw_abort(hecmw_comm_get_comm())
568  endif
569  endif
570 
571  ! ----- Add vector
572  do j = 1, ndof*nn
573  b(iwk(j)) = b(iwk(j)) + vect(j)
574  enddo
575 
576  end subroutine calculate_thermal_load
577 
578 end module m_fstr_ass_load
m_common_struct
This modules defines common structures for fem analysis.
Definition: m_common_struct.f90:6
m_common_struct::fstr_rotinfo_init
subroutine fstr_rotinfo_init(n, rinfo)
Definition: m_common_struct.f90:74
m_common_struct::fstr_rotinfo_finalize
subroutine fstr_rotinfo_finalize(rinfo)
Definition: m_common_struct.f90:92
m_fstr_ass_load
This module provides functions to take into account external load.
Definition: fstr_ass_load.f90:7
m_fstr_ass_load::process_concentrated_loads
subroutine process_concentrated_loads(cstep, ctime, hecMESH, fstrSOLID)
Process concentrated nodal forces (CLOAD)
Definition: fstr_ass_load.f90:71
m_fstr_ass_load::fstr_ass_load
subroutine fstr_ass_load(cstep, ctime, hecMESH, hecMAT, fstrSOLID, fstrPARAM)
This subroutine assmble following external force into fstrSOLIDGL and hecMATB afterwards.
Definition: fstr_ass_load.f90:31
m_fstr_ass_load::update_rhs_vector
subroutine update_rhs_vector(hecMESH, hecMAT, fstrSOLID)
Update right-hand side vector.
Definition: fstr_ass_load.f90:355
m_fstr_ass_load::process_user_loads
subroutine process_user_loads(cstep, fstrSOLID)
Definition: fstr_ass_load.f90:341
m_fstr_ass_load::process_distributed_loads
subroutine process_distributed_loads(cstep, ctime, hecMESH, fstrSOLID)
Definition: fstr_ass_load.f90:178
m_fstr_ass_load::calculate_thermal_load
subroutine calculate_thermal_load(ic_type, nn, xx, yy, zz, tt, tt0, isect, ndof, hecMESH, fstrSOLID, icel, vect, cdsys_ID, coords, iset, pa1, iwk, B)
Calculate thermal load based on element type.
Definition: fstr_ass_load.f90:507
m_fstr_ass_load::process_thermal_loads
subroutine process_thermal_loads(cstep, ctime, hecMESH, hecMAT, fstrSOLID)
Definition: fstr_ass_load.f90:376
m_fstr_elemact
This module provide a function to elemact elements.
Definition: fstr_element_activation.f90:6
m_fstr_elemact::fstr_update_elemact_solid
subroutine fstr_update_elemact_solid(hecMESH, fstrSOLID, cstep, ctime)
Definition: fstr_element_activation.f90:40
m_fstr_precheck
This module provides the entry point for ELEMCHECK (pre-analysis input validation)
Definition: fstr_precheck.f90:6
m_fstr_precheck::fstr_get_thickness
subroutine fstr_get_thickness(hecMESH, mid, thick)
Definition: fstr_precheck.f90:81
m_fstr_spring
This module provides functions to deal with spring force.
Definition: fstr_Spring.f90:7
m_fstr_spring::fstr_update_ndforce_spring
subroutine fstr_update_ndforce_spring(cstep, hecMESH, fstrSOLID, B)
Definition: fstr_Spring.f90:58
m_fstr
This module defines common data and basic structures for analysis.
Definition: m_fstr.F90:15
m_fstr::kel361bbar
integer(kind=kint), parameter kel361bbar
Definition: m_fstr.F90:81
m_fstr::imsg
integer(kind=kint), parameter imsg
Definition: m_fstr.F90:117
m_fstr::kel361up
integer(kind=kint), parameter kel361up
Definition: m_fstr.F90:84
m_fstr::fstr_isloadactive
logical function fstr_isloadactive(fstrSOLID, nbc, cstep)
Definition: m_fstr.F90:1102
m_fstr::kel361fi
integer(kind=kint), parameter kel361fi
Definition: m_fstr.F90:80
m_fstr::fstr_get_amplitude
subroutine fstr_get_amplitude(hecMESH, fstrSOLID, cstep, jj_n_amp, time, value)
Evaluate the amplitude-scaled target value a(t) for static analysis.
Definition: m_fstr.F90:1254
m_fstr::get_coordsys
subroutine get_coordsys(cdsys_ID, hecMESH, fstrSOLID, coords, icel)
This subroutine fetch coords defined by local coordinate system.
Definition: m_fstr.F90:1147
m_fstr::kel361ic
integer(kind=kint), parameter kel361ic
Definition: m_fstr.F90:82
m_fstr::kel361fbar
integer(kind=kint), parameter kel361fbar
Definition: m_fstr.F90:83
m_static_lib
This modules just summarizes all modules used in static analysis.
Definition: static_LIB.f90:6
m_utilities
This module provides aux functions.
Definition: utilities.f90:6
m_utilities::cross_product
subroutine cross_product(v1, v2, vn)
Definition: utilities.f90:406
mmechgauss
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
mreadtemp::read_temperature_result
subroutine, public read_temperature_result(hecMESH, nstep, sstep, rtype, interval, factor, ctime, temp, temp_bak)
Read in temperature distribution from external file.
Definition: readtemp.f90:14
muload
This subroutine read in used-defined loading tangent.
Definition: uload.f90:7
muload::uloading
subroutine uloading(cstep, factor, exForce)
This subroutine take consider of user-defined external loading.
Definition: uload.f90:31
m_common_struct::trotinfo
Definition: m_common_struct.f90:21
m_fstr::fstr_param
FSTR INNER CONTROL PARAMETERS (fstrPARAM)
Definition: m_fstr.F90:161
m_fstr::fstr_solid
Definition: m_fstr.F90:252