FrontISTR  5.8.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) :: factor, 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  jj_n_amp = fstrsolid%CLOAD_ngrp_amp(ig0)
94  if (jj_n_amp <= 0) then ! Amplitude not defined
95  factor = fstrsolid%FACTOR(2)
96  else
97  call table_amp(hecmesh, fstrsolid, cstep, jj_n_amp, ctime, factor)
98  endif
99 
100  if (fstr_isloadactive(fstrsolid, grpid, cstep-1)) factor = 1.0d0
101  ig = fstrsolid%CLOAD_ngrp_ID(ig0)
102  ityp = fstrsolid%CLOAD_ngrp_DOF(ig0)
103  fval = fstrsolid%CLOAD_ngrp_val(ig0)
104  is0 = hecmesh%node_group%grp_index(ig-1) + 1
105  ie0 = hecmesh%node_group%grp_index(ig)
106 
107  if( fstrsolid%CLOAD_ngrp_rotID(ig0) > 0 ) then ! setup torque load information
108  rid = fstrsolid%CLOAD_ngrp_rotID(ig0)
109  if (.not. rinfo%conds(rid)%active) then
110  rinfo%conds(rid)%active = .true.
111  rinfo%conds(rid)%center_ngrp_id = fstrsolid%CLOAD_ngrp_centerID(ig0)
112  rinfo%conds(rid)%torque_ngrp_id = ig
113  endif
114  if (ityp > ndof) ityp = ityp - ndof
115  rinfo%conds(rid)%vec(ityp) = factor * fval
116  cycle
117  endif
118 
119  do ik = is0, ie0
120  in = hecmesh%node_group%grp_item(ik)
121  fstrsolid%GL(ndof*(in-1)+ityp) = fstrsolid%GL(ndof*(in-1)+ityp) + factor*fval
122  enddo
123  enddo
124 
125  !Add torque load to fstrSOLID%GL
126  do rid = 1, n_rot
127  if (.not. rinfo%conds(rid)%active) cycle
128  ! Get number of slave nodes
129  n_nodes = hecmw_ngrp_get_number(hecmesh, rinfo%conds(rid)%torque_ngrp_id)
130 
131  ! Get center node
132  ig = rinfo%conds(rid)%center_ngrp_id
133  do idof = 1, ndof
134  ccoord(idof) = hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, hecmesh%node)
135  cdisp(idof) = hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, fstrsolid%unode)
136  cdisp(idof) = cdisp(idof) + hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, fstrsolid%dunode)
137  enddo
138  ccoord(1:ndof) = ccoord(1:ndof) + cdisp(1:ndof)
139 
140  tval = dsqrt(dot_product(rinfo%conds(rid)%vec(1:ndof), rinfo%conds(rid)%vec(1:ndof)))
141  if (tval < 1.d-16) then
142  write(*,*) '###ERROR### : norm of torque vector must be > 0.0'
143  call hecmw_abort(hecmw_comm_get_comm())
144  endif
145  normal(1:ndof) = rinfo%conds(rid)%vec(1:ndof) / tval
146  tval = tval / dble(n_nodes)
147 
148  ig = rinfo%conds(rid)%torque_ngrp_id
149  is0 = hecmesh%node_group%grp_index(ig-1) + 1
150  ie0 = hecmesh%node_group%grp_index(ig)
151  do ik = is0, ie0
152  in = hecmesh%node_group%grp_item(ik)
153  cdiff(1:ndof) = hecmesh%node(ndof*(in-1)+1:ndof*in) + fstrsolid%unode(ndof*(in-1)+1:ndof*in) &
154  & + fstrsolid%dunode(ndof*(in-1)+1:ndof*in) - ccoord(1:ndof)
155  call cross_product(normal,cdiff,vect(1:ndof))
156  fval = dot_product(vect(1:ndof), vect(1:ndof))
157  if (fval < 1.d-16) then
158  write(*,*) '###ERROR### : torque node is at the same position as that of center node in rotational surface.'
159  call hecmw_abort(hecmw_comm_get_comm())
160  endif
161  vect(1:ndof) = (tval/fval) * vect(1:ndof)
162  fstrsolid%GL(ndof*(in-1)+1:ndof*in) = fstrsolid%GL(ndof*(in-1)+1:ndof*in) + vect(1:ndof)
163  enddo
164  enddo
165  if (n_rot > 0) call fstr_rotinfo_finalize(rinfo)
166 
167  end subroutine process_concentrated_loads
168  !
169  ! -------------------------------------------------------------------
170  ! DLOAD
171  ! -------------------------------------------------------------------
172  subroutine process_distributed_loads(cstep, ctime, hecMESH, fstrSOLID)
173  integer(kind=kint), intent(in) :: cstep
174  real(kind=kreal), intent(in) :: ctime
175  type(hecmwst_local_mesh), intent(in) :: hecmesh
176  type(fstr_solid), intent(inout) :: fstrSOLID
177 
178  integer(kind=kint) :: ndof, ig0, ig, ltype, iS0, iE0, ik, icel, ic_type, nn, is
179  integer(kind=kint) :: isect, id, iset, ihead, nsize, grpid, i, j, jj_n_amp
180  integer(kind=kint) :: iwk(60), nodLocal(20)
181  real(kind=kreal) :: xx(20), yy(20), zz(20), vect(60), params(0:6)
182  real(kind=kreal) :: factor, rho, thick, pa1
183  logical :: fg_surf
184  type(tmaterial), pointer :: material
185 
186  ndof = hecmesh%n_dof
187 
188  do ig0 = 1, fstrsolid%DLOAD_ngrp_tot
189  grpid = fstrsolid%DLOAD_ngrp_GRPID(ig0)
190  if (.not. fstr_isloadactive(fstrsolid, grpid, cstep)) cycle
191  jj_n_amp = fstrsolid%DLOAD_ngrp_amp(ig0)
192  if (jj_n_amp <= 0) then ! Amplitude not defined
193  factor = fstrsolid%factor(2)
194  else
195  call table_amp(hecmesh, fstrsolid, cstep, jj_n_amp, ctime, factor)
196  endif
197 
198  if (fstr_isloadactive(fstrsolid, grpid, cstep-1)) factor = 1.0d0
199  ig = fstrsolid%DLOAD_ngrp_ID(ig0)
200  ltype = fstrsolid%DLOAD_ngrp_LID(ig0)
201  do i = 0, 6
202  params(i) = fstrsolid%DLOAD_ngrp_params(i, ig0)
203  enddo
204  ! ----- START & END
205  fg_surf = (ltype == 100)
206  if( fg_surf ) then ! surface group
207  is0 = hecmesh%surf_group%grp_index(ig-1) + 1
208  ie0 = hecmesh%surf_group%grp_index(ig)
209  else ! element group
210  is0 = hecmesh%elem_group%grp_index(ig-1) + 1
211  ie0 = hecmesh%elem_group%grp_index(ig)
212  endif
213  do ik = is0, ie0
214  if( fg_surf ) then ! surface group
215  ltype = hecmesh%surf_group%grp_item(2*ik) * 10
216  icel = hecmesh%surf_group%grp_item(2*ik-1)
217  ic_type = hecmesh%elem_type(icel)
218  else ! element group
219  icel = hecmesh%elem_group%grp_item(ik)
220  ic_type = hecmesh%elem_type(icel)
221  endif
222 
223  !ELEMENT ACTIVATION
224  if( fstrsolid%elements(icel)%elemact_flag == kelact_inactive ) cycle
225 
226  if (hecmw_is_etype_link(ic_type)) cycle
227  if (hecmw_is_etype_patch(ic_type)) cycle
228  ! if( ic_type==3422 ) ic_type=342
229  nn = hecmw_get_max_node(ic_type)
230  ! ----- node ID
231  is = hecmesh%elem_node_index(icel-1)
232  if (fstrsolid%DLOAD_follow == 0) then
233  do j = 1, nn
234  nodlocal(j) = hecmesh%elem_node_item (is+j)
235  ! ----- nodal coordinate
236  xx(j) = hecmesh%node( 3*nodlocal(j)-2 )
237  yy(j) = hecmesh%node( 3*nodlocal(j)-1 )
238  zz(j) = hecmesh%node( 3*nodlocal(j) )
239  ! ----- create iwk array ***
240  do i = 1, ndof
241  iwk( ndof*(j-1)+i ) = ndof*( nodlocal(j)-1 )+i
242  enddo
243  enddo
244  else
245  do j = 1, nn
246  nodlocal(j) = hecmesh%elem_node_item (is+j)
247  ! ----- nodal coordinate
248  if (ndof==2) then
249  xx(j) = hecmesh%node( 3*nodlocal(j)-2 )+fstrsolid%unode( 2*nodlocal(j)-1 )+fstrsolid%dunode( 2*nodlocal(j)-1 )
250  yy(j) = hecmesh%node( 3*nodlocal(j)-1 )+fstrsolid%unode( 2*nodlocal(j) )+fstrsolid%dunode( 2*nodlocal(j) )
251  else if (ndof==3) then
252  xx(j) = hecmesh%node( 3*nodlocal(j)-2 )+fstrsolid%unode( 3*nodlocal(j)-2 )+fstrsolid%dunode( 3*nodlocal(j)-2 )
253  yy(j) = hecmesh%node( 3*nodlocal(j)-1 )+fstrsolid%unode( 3*nodlocal(j)-1 )+fstrsolid%dunode( 3*nodlocal(j)-1 )
254  zz(j) = hecmesh%node( 3*nodlocal(j) )+fstrsolid%unode( 3*nodlocal(j) )+fstrsolid%dunode( 3*nodlocal(j) )
255  else if (ndof==6) then
256  xx(j) = hecmesh%node( 3*nodlocal(j)-2 )+fstrsolid%unode( 6*nodlocal(j)-5 )+fstrsolid%dunode( 6*nodlocal(j)-5 )
257  yy(j) = hecmesh%node( 3*nodlocal(j)-1 )+fstrsolid%unode( 6*nodlocal(j)-4 )+fstrsolid%dunode( 6*nodlocal(j)-4 )
258  zz(j) = hecmesh%node( 3*nodlocal(j) )+fstrsolid%unode( 6*nodlocal(j)-3 )+fstrsolid%dunode( 6*nodlocal(j)-3 )
259  endif
260  ! ----- create iwk array ***
261  do i = 1, ndof
262  iwk( ndof*(j-1)+i ) = ndof*( nodlocal(j)-1 )+i
263  enddo
264  enddo
265  endif
266  ! ----- section ID
267  isect = hecmesh%section_ID(icel)
268  ! ----- Get Properties
269  material => fstrsolid%elements(icel)%gausses(1)%pMaterial
270  rho = material%variables(m_density)
271  call fstr_get_thickness(hecmesh, isect, thick)
272  ! ----- Section Data
273  if (ndof == 2) then
274  id = hecmesh%section%sect_opt(isect)
275  if (id == 0) then
276  iset = 1
277  else if (id == 1) then
278  iset = 0
279  else if (id == 2) then
280  iset = 2
281  endif
282  pa1 = 1.d0
283  endif
284  ! ----- Create local stiffness
285  if (ic_type==301)then
286  ihead = hecmesh%section%sect_R_index(isect-1)
287  call dl_c1(ic_type,nn,xx(1:nn),yy(1:nn),zz(1:nn),rho,thick,ltype,params,vect(1:nn*ndof),nsize)
288 
289  elseif( ic_type == 241 .or. ic_type == 242 .or. ic_type == 231 .or. ic_type == 232 .or. ic_type == 2322 ) then
290  call dl_c2(ic_type,nn,xx(1:nn),yy(1:nn),rho,pa1,ltype,params,vect(1:nn*ndof),nsize,iset)
291 
292  else if ( ic_type == 341 .or. ic_type == 351 .or. ic_type == 361 .or. &
293  ic_type == 342 .or. ic_type == 352 .or. ic_type == 362 ) then
294  call dl_c3(ic_type,nn,xx(1:nn),yy(1:nn),zz(1:nn),rho,ltype,params,vect(1:nn*ndof),nsize)
295 
296  else if ( ic_type == 641 ) then
297  ihead = hecmesh%section%sect_R_index(isect-1)
298  call dl_beam_641(ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), rho, ltype, params, &
299  hecmesh%section%sect_R_item(ihead+1:), vect(1:nn*ndof), nsize)
300 
301  else if( ( ic_type == 741 ) .or. ( ic_type == 743 ) .or. ( ic_type == 731 ) ) then
302  call dl_shell(ic_type, nn, ndof, xx, yy, zz, rho, thick, ltype, params, vect, nsize, fstrsolid%elements(icel)%gausses)
303 
304  else if( ( ic_type==761 ) .or. ( ic_type==781 ) ) then
305  call dl_shell_33(ic_type, nn, ndof, xx, yy, zz, rho, thick, ltype, params, vect, nsize, &
306  fstrsolid%elements(icel)%gausses)
307 
308  else
309  nsize = 0
310  write(*,*)"### WARNING: DLOAD",ic_type
311 
312  endif
313  ! ----- Add vector
314  do j = 1, nsize
315  fstrsolid%GL(iwk(j)) = fstrsolid%GL(iwk(j)) + factor * vect(j)
316  enddo
317  enddo
318  enddo
319 
320  end subroutine process_distributed_loads
321 
322  ! -----Uload
323  subroutine process_user_loads(cstep, fstrSOLID)
324  integer(kind=kint), intent(in) :: cstep
325  type(fstr_solid), intent(inout) :: fstrSOLID
326 
327  real(kind=kreal) :: factor
328 
329  factor = fstrsolid%factor(2)
330  call uloading(cstep, factor, fstrsolid%GL)
331 
332  end subroutine process_user_loads
333 
334  !======================================================================!
336  !======================================================================!
337  subroutine update_rhs_vector(hecMESH, hecMAT, fstrSOLID)
338  type(hecmwst_local_mesh), intent(in) :: hecmesh
339  type(hecmwst_matrix), intent(inout) :: hecMAT
340  type(fstr_solid), intent(inout) :: fstrSOLID
341 
342  integer(kind=kint) :: i
343 
344  do i = 1, hecmesh%n_node * hecmesh%n_dof
345  hecmat%B(i) = fstrsolid%GL(i) - fstrsolid%QFORCE(i)
346  enddo
347 
348  do i = 1, hecmat%NDOF * hecmat%NP
349  !thermal load is not considered
350  fstrsolid%EFORCE(i) = fstrsolid%GL(i)
351  enddo
352 
353  end subroutine update_rhs_vector
354 
355  ! -------------------------------------------------------------------
356  ! TLOAD : THERMAL LOAD USING TEMPERATURE
357  ! -------------------------------------------------------------------
358  subroutine process_thermal_loads(cstep, ctime, hecMESH, hecMAT, fstrSOLID)
359  integer(kind=kint), intent(in) :: cstep
360  real(kind=kreal), intent(in) :: ctime
361  type(hecmwst_local_mesh), intent(in) :: hecmesh
362  type(hecmwst_matrix), intent(inout) :: hecMAT
363  type(fstr_solid), intent(inout) :: fstrSOLID
364 
365  integer(kind=kint) :: ndof, ig0, ig, iS0, iE0, ik, in, grpid
366  integer(kind=kint) :: itype, is, iE, icel, ic_type, nn, isect, cdsys_ID, id, iset
367  integer(kind=kint) :: i, j, ihead, tstep, nodLocal(20), iwk(60)
368  real(kind=kreal) :: factor, fval, pa1
369  real(kind=kreal) :: xx(20), yy(20), zz(20), tt(20), tt0(20), coords(3,3), vect(60)
370  real(kind=kreal) :: local_coords(3,3) ! Local copy for coordinate transformation
371 
372  ndof = hecmesh%n_dof
373 
374  if (fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0) then
375  do ig0 = 1, fstrsolid%TEMP_ngrp_tot
376  grpid = fstrsolid%TEMP_ngrp_GRPID(ig0)
377  if (.not. fstr_isloadactive(fstrsolid, grpid, cstep)) cycle
378  factor = fstrsolid%factor(2)
379  if (fstr_isloadactive(fstrsolid, grpid, cstep-1)) factor = 1.0d0
380  ig = fstrsolid%TEMP_ngrp_ID(ig0)
381  fval = fstrsolid%TEMP_ngrp_val(ig0)
382  is0 = hecmesh%node_group%grp_index(ig-1) + 1
383  ie0 = hecmesh%node_group%grp_index(ig)
384  do ik = is0, ie0
385  in = hecmesh%node_group%grp_item(ik)
386  pa1 = fstrsolid%temp_bak(in)
387  fstrsolid%temperature(in) = pa1 + (fval - pa1) * factor
388  enddo
389  enddo
390 
391  if (fstrsolid%TEMP_irres > 0) then
392  call read_temperature_result(hecmesh, fstrsolid%TEMP_irres, fstrsolid%TEMP_tstep, &
393  & fstrsolid%TEMP_rtype, fstrsolid%TEMP_interval, fstrsolid%TEMP_factor, ctime, &
394  & fstrsolid%temperature, fstrsolid%temp_bak)
395  endif
396  endif
397 
398  ! ----- elemact element
399  if( fstrsolid%elemact%ELEMACT_egrp_tot > 0 ) &
400  & call fstr_update_elemact_solid( hecmesh, fstrsolid, cstep, ctime )
401 
402  if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0 ) then
403  ! ----- element TYPE loop.
404  do itype = 1, hecmesh%n_elem_type
405  is = hecmesh%elem_type_index(itype-1) + 1
406  ie = hecmesh%elem_type_index(itype)
407  ic_type = hecmesh%elem_type_item(itype)
408  if (hecmw_is_etype_link(ic_type)) cycle
409  if (hecmw_is_etype_patch(ic_type)) cycle
410  ! ----- Set number of nodes
411  nn = hecmw_get_max_node(ic_type)
412 
413  ! ----- element loop
414  do icel = is, ie
415 
416  !ELEMENT ACTIVATION
417  if( fstrsolid%elements(icel)%elemact_flag == kelact_inactive ) cycle
418 
419  ! ----- node ID
420  is = hecmesh%elem_node_index(icel-1)
421  do j = 1, nn
422  nodlocal(j) = hecmesh%elem_node_item(is+j)
423  ! ----- nodal coordinate
424  if (ndof == 2) then
425  xx(j) = hecmesh%node(3*nodlocal(j)-2) + fstrsolid%unode(ndof*nodlocal(j)-1)
426  yy(j) = hecmesh%node(3*nodlocal(j)-1) + fstrsolid%unode(ndof*nodlocal(j))
427  else if (ndof == 3) then
428  xx(j) = hecmesh%node(3*nodlocal(j)-2) + fstrsolid%unode(ndof*nodlocal(j)-2)
429  yy(j) = hecmesh%node(3*nodlocal(j)-1) + fstrsolid%unode(ndof*nodlocal(j)-1)
430  zz(j) = hecmesh%node(3*nodlocal(j)) + fstrsolid%unode(ndof*nodlocal(j))
431  endif
432  tt0(j) = fstrsolid%last_temp(nodlocal(j))
433  tt(j) = fstrsolid%temperature(nodlocal(j))
434  ! ----- create iwk array ***
435  do i = 1, ndof
436  iwk(ndof*(j-1)+i) = ndof*(nodlocal(j)-1)+i
437  enddo
438  enddo
439 
440  ! ----- section Data
441  isect = hecmesh%section_ID(icel)
442  cdsys_id = hecmesh%section%sect_orien_ID(isect)
443  call get_coordsys(cdsys_id, hecmesh, fstrsolid, coords)
444 
445  if (ndof == 2) then
446  id = hecmesh%section%sect_opt(isect)
447  if (id == 0) then
448  iset = 1
449  else if (id == 1) then
450  iset = 0
451  else if (id == 2) then
452  iset = 2
453  endif
454  pa1 = 1.0d0
455  endif
456 
457  if (ic_type == 641) then
458  isect = hecmesh%section_ID(icel)
459  ihead = hecmesh%section%sect_R_index(isect-1)
460 
461  call tload_beam_641( ic_type, nn, ndof, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
462  fstrsolid%elements(icel)%gausses, hecmesh%section%sect_R_item(ihead+1:), &
463  vect(1:nn*ndof))
464 
465  do j = 1, ndof*nn
466  hecmat%B(iwk(j)) = hecmat%B(iwk(j)) + vect(j)
467  enddo
468  cycle
469  endif
470 
471  ! Local copy of coordinate data
472  local_coords = coords
473 
474  ! Calculate thermal load based on element type
475  call calculate_thermal_load(ic_type, nn, xx, yy, zz, tt, tt0, isect, ndof, &
476  hecmesh, fstrsolid, icel, vect, cdsys_id, local_coords, &
477  iset, pa1, iwk, hecmat%B)
478  enddo
479  enddo
480  endif
481 
482  end subroutine process_thermal_loads
483 
484  !======================================================================!
486  !======================================================================!
487  subroutine calculate_thermal_load(ic_type, nn, xx, yy, zz, tt, tt0, isect, ndof, &
488  hecMESH, fstrSOLID, icel, vect, cdsys_ID, coords, &
489  iset, pa1, iwk, B)
490  integer(kind=kint), intent(in) :: ic_type, nn, isect, ndof, cdsys_ID, iset
491  real(kind=kreal), intent(in) :: xx(*), yy(*), zz(*), tt(*), tt0(*), pa1
492  real(kind=kreal), intent(inout) :: coords(3,3) ! Changed INTENT from IN to INOUT
493  type(hecmwst_local_mesh), intent(in) :: hecmesh
494  type(fstr_solid), intent(in) :: fstrSOLID
495  integer(kind=kint), intent(in) :: icel
496  real(kind=kreal), intent(out) :: vect(*)
497  integer(kind=kint), intent(in) :: iwk(*)
498  real(kind=kreal), intent(inout) :: b(*)
499 
500  integer(kind=kint) :: j, myrank
501 
502  myrank = 0
503 
504  ! 2D elements
505  if (ic_type == 241 .or. ic_type == 242 .or. ic_type == 231 .or. ic_type == 232) then
506  call tload_c2(ic_type, nn, xx(1:nn), yy(1:nn), tt(1:nn), tt0(1:nn), &
507  fstrsolid%elements(icel)%gausses, pa1, iset, vect(1:nn*2))
508 
509  else if (ic_type == 361) then
510  if (fstrsolid%sections(isect)%elemopt361 == kel361fi) then
511  call tload_c3 &
512  ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
513  fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords)
514  else if (fstrsolid%sections(isect)%elemopt361 == kel361bbar) then
515  call tload_c3d8bbar &
516  ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
517  fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords)
518  else if (fstrsolid%sections(isect)%elemopt361 == kel361ic) then
519  call tload_c3d8ic &
520  ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
521  fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords)
522  else if (fstrsolid%sections(isect)%elemopt361 == kel361fbar) then
523  call tload_c3d8fbar &
524  ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
525  fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords)
526  endif
527 
528  else if (ic_type == 341 .or. ic_type == 351 .or. &
529  ic_type == 342 .or. ic_type == 352 .or. ic_type == 362) then
530  call tload_c3 &
531  ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
532  fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords)
533 
534  else if (ic_type == 741 .or. ic_type == 743 .or. ic_type == 731) then
535  if (myrank == 0) then
536  write(imsg,*) '*------------------------', &
537  '-------------------*'
538  write(imsg,*) ' Thermal loading option for shell elements', &
539  'not yet available.'
540  write(imsg,*) '*------------------------', &
541  '-------------------*'
542  call hecmw_abort(hecmw_comm_get_comm())
543  endif
544  endif
545 
546  ! ----- Add vector
547  do j = 1, ndof*nn
548  b(iwk(j)) = b(iwk(j)) + vect(j)
549  enddo
550 
551  end subroutine calculate_thermal_load
552 
553 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:99
m_common_struct::fstr_rotinfo_finalize
subroutine fstr_rotinfo_finalize(rinfo)
Definition: m_common_struct.f90:117
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:338
m_fstr_ass_load::process_user_loads
subroutine process_user_loads(cstep, fstrSOLID)
Definition: fstr_ass_load.f90:324
m_fstr_ass_load::process_distributed_loads
subroutine process_distributed_loads(cstep, ctime, hecMESH, fstrSOLID)
Definition: fstr_ass_load.f90:173
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:490
m_fstr_ass_load::process_thermal_loads
subroutine process_thermal_loads(cstep, ctime, hecMESH, hecMAT, fstrSOLID)
Definition: fstr_ass_load.f90:359
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:14
m_fstr_precheck
This module provides function to check input data of IFSTR solver.
Definition: fstr_precheck.f90:6
m_fstr_precheck::fstr_get_thickness
subroutine fstr_get_thickness(hecMESH, mid, thick)
Definition: fstr_precheck.f90:50
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:46
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:79
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:1093
m_fstr::imsg
integer(kind=kint), parameter imsg
Definition: m_fstr.f90:111
m_fstr::fstr_isloadactive
logical function fstr_isloadactive(fstrSOLID, nbc, cstep)
Definition: m_fstr.f90:1048
m_fstr::kel361fi
integer(kind=kint), parameter kel361fi
Definition: m_fstr.f90:78
m_fstr::kel361ic
integer(kind=kint), parameter kel361ic
Definition: m_fstr.f90:80
m_fstr::table_amp
subroutine table_amp(hecMESH, fstrSOLID, cstep, jj_n_amp, time, f_t)
Definition: m_fstr.f90:1176
m_fstr::kel361fbar
integer(kind=kint), parameter kel361fbar
Definition: m_fstr.f90:81
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:367
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:155
m_fstr::fstr_solid
Definition: m_fstr.f90:240