18 subroutine fstr_ass_load(cstep, ctime, hecMESH, hecMAT, fstrSOLID, fstrPARAM)
29 integer(kind=kint),
intent(in) :: cstep
30 real(kind=kreal),
intent(in) :: ctime
31 type(hecmwst_matrix),
intent(inout) :: hecmat
32 type(hecmwst_local_mesh),
intent(in) :: hecMESH
36 real(kind=kreal) :: xx(20), yy(20), zz(20)
37 real(kind=kreal) :: params(0:6)
38 real(kind=kreal) :: vect(60)
39 integer(kind=kint) :: iwk(60)
40 integer(kind=kint) :: nodLocal(20)
41 real(kind=kreal) :: tt(20), tt0(20), coords(3, 3), factor
42 integer(kind=kint) :: ndof, ig0, ig, ityp, ltype, iS0, iE0, ik, in, i, j
43 integer(kind=kint) :: icel, ic_type, nn, is, isect, id, iset, nsize
44 integer(kind=kint) :: itype, iE, ierror, grpid, cdsys_ID
45 real(kind=kreal) :: fval, rho, thick, pa1
47 integer(kind=kint) :: tstep
48 type(tmaterial),
pointer :: material
49 integer(kind=kint) :: ihead
53 integer(kind=kint) :: n_rot, rid, n_nodes, idof
55 real(kind=kreal) :: tval, normal(3), direc(3), ccoord(3), cdisp(3), cdiff(3)
62 n_rot = fstrsolid%CLOAD_ngrp_rot
65 fstrsolid%GL(:) = 0.0d0
66 fstrsolid%EFORCE(:) = 0.0d0
67 do ig0 = 1, fstrsolid%CLOAD_ngrp_tot
68 grpid = fstrsolid%CLOAD_ngrp_GRPID(ig0)
70 factor = fstrsolid%factor(2)
72 ig = fstrsolid%CLOAD_ngrp_ID(ig0)
73 ityp = fstrsolid%CLOAD_ngrp_DOF(ig0)
74 fval = fstrsolid%CLOAD_ngrp_val(ig0)
75 is0 = hecmesh%node_group%grp_index(ig-1) + 1
76 ie0 = hecmesh%node_group%grp_index(ig )
78 if( fstrsolid%CLOAD_ngrp_rotID(ig0) > 0 )
then
79 rid = fstrsolid%CLOAD_ngrp_rotID(ig0)
80 if( .not. rinfo%conds(rid)%active )
then
81 rinfo%conds(rid)%active = .true.
82 rinfo%conds(rid)%center_ngrp_id = fstrsolid%CLOAD_ngrp_centerID(ig0)
83 rinfo%conds(rid)%torque_ngrp_id = ig
85 if( ityp>ndof ) ityp = ityp-ndof
86 rinfo%conds(rid)%vec(ityp) = factor*fval
91 in = hecmesh%node_group%grp_item(ik)
92 fstrsolid%GL(ndof*(in-1)+ityp) = fstrsolid%GL(ndof*(in-1)+ityp)+factor*fval
98 if( .not. rinfo%conds(rid)%active ) cycle
100 n_nodes = hecmw_ngrp_get_number(hecmesh, rinfo%conds(rid)%torque_ngrp_id)
103 ig = rinfo%conds(rid)%center_ngrp_id
105 ccoord(idof) = hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, hecmesh%node)
106 cdisp(idof) = hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, fstrsolid%unode)
107 cdisp(idof) = cdisp(idof) + hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, fstrsolid%dunode)
109 ccoord(1:ndof) = ccoord(1:ndof) + cdisp(1:ndof)
111 tval = dsqrt(dot_product(rinfo%conds(rid)%vec(1:ndof),rinfo%conds(rid)%vec(1:ndof)))
112 if( tval < 1.d-16 )
then
113 write(*,*)
'###ERROR### : norm of torque vector must be > 0.0'
114 call hecmw_abort( hecmw_comm_get_comm() )
116 normal(1:ndof) = rinfo%conds(rid)%vec(1:ndof)/tval
117 tval = tval/dble(n_nodes)
119 ig = rinfo%conds(rid)%torque_ngrp_id
120 is0 = hecmesh%node_group%grp_index(ig-1) + 1
121 ie0 = hecmesh%node_group%grp_index(ig )
123 in = hecmesh%node_group%grp_item(ik)
124 cdiff(1:ndof) = hecmesh%node(ndof*(in-1)+1:ndof*in)+fstrsolid%unode(ndof*(in-1)+1:ndof*in) &
125 & +fstrsolid%dunode(ndof*(in-1)+1:ndof*in)-ccoord(1:ndof)
127 fval = dot_product(vect(1:ndof),vect(1:ndof))
128 if( fval < 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() )
132 vect(1:ndof) = (tval/fval)*vect(1:ndof)
133 fstrsolid%GL(ndof*(in-1)+1:ndof*in) = fstrsolid%GL(ndof*(in-1)+1:ndof*in)+vect(1:ndof)
141 do ig0 = 1, fstrsolid%DLOAD_ngrp_tot
142 grpid = fstrsolid%DLOAD_ngrp_GRPID(ig0)
144 factor = fstrsolid%factor(2)
146 ig = fstrsolid%DLOAD_ngrp_ID(ig0)
147 ltype = fstrsolid%DLOAD_ngrp_LID(ig0)
149 params(i)= fstrsolid%DLOAD_ngrp_params(i,ig0)
152 fg_surf = (ltype == 100)
154 is0 = hecmesh%surf_group%grp_index(ig-1) + 1
155 ie0 = hecmesh%surf_group%grp_index(ig )
157 is0 = hecmesh%elem_group%grp_index(ig-1) + 1
158 ie0 = hecmesh%elem_group%grp_index(ig )
162 ltype = hecmesh%surf_group%grp_item(2*ik)*10
163 icel = hecmesh%surf_group%grp_item(2*ik-1)
164 ic_type = hecmesh%elem_type(icel)
166 icel = hecmesh%elem_group%grp_item(ik)
167 ic_type = hecmesh%elem_type(icel)
169 if( hecmw_is_etype_link(ic_type) ) cycle
170 if( hecmw_is_etype_patch(ic_type) ) cycle
172 nn = hecmw_get_max_node(ic_type)
174 is = hecmesh%elem_node_index(icel-1)
175 if( fstrsolid%DLOAD_follow == 0 )
then
177 nodlocal(j) = hecmesh%elem_node_item (is+j)
179 xx(j) = hecmesh%node( 3*nodlocal(j)-2 )
180 yy(j) = hecmesh%node( 3*nodlocal(j)-1 )
181 zz(j) = hecmesh%node( 3*nodlocal(j) )
184 iwk( ndof*(j-1)+i ) = ndof*( nodlocal(j)-1 )+i
189 nodlocal(j) = hecmesh%elem_node_item (is+j)
192 xx(j) = hecmesh%node( 3*nodlocal(j)-2 )+fstrsolid%unode( 2*nodlocal(j)-1 )+fstrsolid%dunode( 2*nodlocal(j)-1 )
193 yy(j) = hecmesh%node( 3*nodlocal(j)-1 )+fstrsolid%unode( 2*nodlocal(j) )+fstrsolid%dunode( 2*nodlocal(j) )
194 else if (ndof==3)
then
195 xx(j) = hecmesh%node( 3*nodlocal(j)-2 )+fstrsolid%unode( 3*nodlocal(j)-2 )+fstrsolid%dunode( 3*nodlocal(j)-2 )
196 yy(j) = hecmesh%node( 3*nodlocal(j)-1 )+fstrsolid%unode( 3*nodlocal(j)-1 )+fstrsolid%dunode( 3*nodlocal(j)-1 )
197 zz(j) = hecmesh%node( 3*nodlocal(j) )+fstrsolid%unode( 3*nodlocal(j) )+fstrsolid%dunode( 3*nodlocal(j) )
198 else if (ndof==6)
then
199 xx(j) = hecmesh%node( 3*nodlocal(j)-2 )+fstrsolid%unode( 6*nodlocal(j)-5 )+fstrsolid%dunode( 6*nodlocal(j)-5 )
200 yy(j) = hecmesh%node( 3*nodlocal(j)-1 )+fstrsolid%unode( 6*nodlocal(j)-4 )+fstrsolid%dunode( 6*nodlocal(j)-4 )
201 zz(j) = hecmesh%node( 3*nodlocal(j) )+fstrsolid%unode( 6*nodlocal(j)-3 )+fstrsolid%dunode( 6*nodlocal(j)-3 )
205 iwk( ndof*(j-1)+i ) = ndof*( nodlocal(j)-1 )+i
210 isect = hecmesh%section_ID(icel)
212 material => fstrsolid%elements(icel)%gausses(1)%pMaterial
213 rho = material%variables(m_density)
217 id=hecmesh%section%sect_opt(isect)
220 else if( id == 1 )
then
222 else if( id == 2)
then
228 if (ic_type==301)
then
229 ihead = hecmesh%section%sect_R_index(isect-1)
230 call dl_c1(ic_type,nn,xx(1:nn),yy(1:nn),zz(1:nn),rho,thick,ltype,params,vect(1:nn*ndof),nsize)
232 elseif( ic_type == 241 .or. ic_type == 242 .or. ic_type == 231 .or. ic_type == 232 .or. ic_type == 2322 )
then
233 call dl_c2(ic_type,nn,xx(1:nn),yy(1:nn),rho,pa1,ltype,params,vect(1:nn*ndof),nsize,iset)
235 else if ( ic_type == 341 .or. ic_type == 351 .or. ic_type == 361 .or. &
236 ic_type == 342 .or. ic_type == 352 .or. ic_type == 362 )
then
237 call dl_c3(ic_type,nn,xx(1:nn),yy(1:nn),zz(1:nn),rho,ltype,params,vect(1:nn*ndof),nsize)
239 else if ( ic_type == 641 )
then
240 ihead = hecmesh%section%sect_R_index(isect-1)
241 call dl_beam_641(ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), rho, ltype, params, &
242 hecmesh%section%sect_R_item(ihead+1:), vect(1:nn*ndof), nsize)
244 else if( ( ic_type == 741 ) .or. ( ic_type == 743 ) .or. ( ic_type == 731 ) )
then
245 call dl_shell(ic_type, nn, ndof, xx, yy, zz, rho, thick, ltype, params, vect, nsize, fstrsolid%elements(icel)%gausses)
247 else if( ( ic_type==761 ) .or. ( ic_type==781 ) )
then
248 call dl_shell_33(ic_type, nn, ndof, xx, yy, zz, rho, thick, ltype, params, vect, nsize, &
249 fstrsolid%elements(icel)%gausses)
253 write(*,*)
"### WARNING: DLOAD",ic_type
258 fstrsolid%GL( iwk(j) )=fstrsolid%GL( iwk(j) )+factor*vect(j)
264 call uloading( cstep, factor, fstrsolid%GL )
269 if( hecmesh%n_dof == 3 )
then
270 call hecmw_update_3_r (hecmesh,fstrsolid%GL,hecmesh%n_node)
271 else if( hecmesh%n_dof == 2 )
then
272 call hecmw_update_2_r (hecmesh,fstrsolid%GL,hecmesh%n_node)
275 call hecmw_mat_clear_b( hecmat )
276 do i=1, hecmesh%n_node* hecmesh%n_dof
277 hecmat%B(i)=fstrsolid%GL(i)-fstrsolid%QFORCE(i)
280 do i=1, hecmat%NDOF*hecmat%NP
282 fstrsolid%EFORCE(i) = fstrsolid%GL(i)
291 if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0 )
then
292 do ig0 = 1, fstrsolid%TEMP_ngrp_tot
293 grpid = fstrsolid%TEMP_ngrp_GRPID(ig0)
295 factor = fstrsolid%factor(2)
297 ig = fstrsolid%TEMP_ngrp_ID(ig0)
298 fval =fstrsolid%TEMP_ngrp_val(ig0)
299 is0 = hecmesh%node_group%grp_index(ig-1)+1
300 ie0 = hecmesh%node_group%grp_index(ig )
302 in = hecmesh%node_group%grp_item(ik)
303 pa1 = fstrsolid%temp_bak( in )
304 fstrsolid%temperature( in ) = pa1+(fval-pa1)*factor
308 if( fstrsolid%TEMP_irres > 0 )
then
310 & fstrsolid%TEMP_rtype, fstrsolid%TEMP_interval, fstrsolid%TEMP_factor, ctime, &
311 & fstrsolid%temperature, fstrsolid%temp_bak)
315 do itype = 1, hecmesh%n_elem_type
317 is = hecmesh%elem_type_index(itype-1)+1
318 ie = hecmesh%elem_type_index(itype )
319 ic_type = hecmesh%elem_type_item(itype)
320 if( hecmw_is_etype_link(ic_type) ) cycle
321 if( hecmw_is_etype_patch(ic_type) ) cycle
323 nn = hecmw_get_max_node(ic_type)
329 is= hecmesh%elem_node_index(icel-1)
331 nodlocal(j)=hecmesh%elem_node_item(is+j)
334 xx(j)=hecmesh%node(3*nodlocal(j)-2)+fstrsolid%unode(ndof*nodlocal(j)-1)
335 yy(j)=hecmesh%node(3*nodlocal(j)-1)+fstrsolid%unode(ndof*nodlocal(j) )
336 else if (ndof==3)
then
337 xx(j)=hecmesh%node(3*nodlocal(j)-2)+fstrsolid%unode(ndof*nodlocal(j)-2)
338 yy(j)=hecmesh%node(3*nodlocal(j)-1)+fstrsolid%unode(ndof*nodlocal(j)-1)
339 zz(j)=hecmesh%node(3*nodlocal(j) )+fstrsolid%unode(ndof*nodlocal(j))
341 tt0(j)=fstrsolid%last_temp( nodlocal(j) )
342 tt(j) = fstrsolid%temperature( nodlocal(j) )
345 iwk(ndof*(j-1)+i)=ndof*(nodlocal(j)-1)+i
350 isect= hecmesh%section_ID(icel)
351 cdsys_id = hecmesh%section%sect_orien_ID(isect)
355 id=hecmesh%section%sect_opt(isect)
358 else if( id == 1 )
then
360 else if( id == 2 )
then
366 if( ic_type == 641 )
then
368 isect= hecmesh%section_ID(icel)
369 ihead = hecmesh%section%sect_R_index(isect-1)
371 call tload_beam_641( ic_type, nn, ndof, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
372 fstrsolid%elements(icel)%gausses, hecmesh%section%sect_R_item(ihead+1:), &
376 hecmat%B( iwk(j) ) = hecmat%B( iwk(j) )+vect(j)
382 if(ic_type == 241 .or. ic_type == 242 .or. ic_type == 231 .or. ic_type == 232 )
then
383 call tload_c2( ic_type, nn, xx(1:nn), yy(1:nn), tt(1:nn), tt0(1:nn), &
384 fstrsolid%elements(icel)%gausses,pa1, iset, vect(1:nn*2) )
386 else if( ic_type == 361 )
then
387 if( fstrsolid%sections(isect)%elemopt361 ==
kel361fi )
then
389 ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
390 fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords )
391 else if( fstrsolid%sections(isect)%elemopt361 ==
kel361bbar )
then
392 call tload_c3d8bbar &
393 ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
394 fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords )
395 else if( fstrsolid%sections(isect)%elemopt361 ==
kel361ic )
then
397 ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
398 fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords )
399 else if( fstrsolid%sections(isect)%elemopt361 ==
kel361fbar )
then
400 call tload_c3d8fbar &
401 ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
402 fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords )
405 else if( ic_type == 341 .or. ic_type == 351 .or. &
406 ic_type == 342 .or. ic_type == 352 .or. ic_type == 362 )
then
408 ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), &
409 fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords )
411 else if( ic_type == 741 .or. ic_type == 743 .or. ic_type == 731 )
then
413 write(
imsg,*)
'*------------------------', &
414 '-------------------*'
415 write(
imsg,*)
' Thermal loading option for shell elements', &
417 write(
imsg,*)
'*------------------------', &
418 '-------------------*'
419 call hecmw_abort( hecmw_comm_get_comm())
427 hecmat%B( iwk(j) ) = hecmat%B( iwk(j) )+vect(j)
437 if(
associated( fstrsolid%contacts ) .and. fstrparam%contact_algo ==
kcaalagrange )
then
438 do i = 1,
size(fstrsolid%contacts)
439 call ass_contact_force( fstrsolid%contacts(i), hecmesh%node, fstrsolid%unode, hecmat%B )
This modules defines common structures for fem analysis.
subroutine fstr_rotinfo_init(n, rinfo)
subroutine fstr_rotinfo_finalize(rinfo)
This module provides functions to take into acount external load.
subroutine fstr_ass_load(cstep, ctime, hecmesh, hecmat, fstrsolid, fstrparam)
This subroutine assmble following external force into fstrSOLIDGL and hecMATB afterwards.
This module provides function to check input data of IFSTR solver.
subroutine fstr_get_thickness(hecmesh, mid, thick)
This module provides functions to deal with spring force.
subroutine fstr_update_ndforce_spring(cstep, hecmesh, fstrsolid, b)
This module defined coomon data and basic structures for analysis.
integer(kind=kint), parameter kel361bbar
integer(kind=kint) myrank
PARALLEL EXECUTION.
integer(kind=kint), parameter imsg
logical function fstr_isloadactive(fstrsolid, nbc, cstep)
subroutine get_coordsys(cdsys_id, hecmesh, fstrsolid, coords)
This subroutine fetch coords defined by local coordinate system.
integer(kind=kint), parameter kel361fi
section control
integer(kind=kint), parameter kel361ic
integer(kind=kint), parameter kcaalagrange
integer(kind=kint), parameter kel361fbar
This modules just summarizes all modules used in static analysis.
This module provides aux functions.
subroutine cross_product(v1, v2, vn)
This modules defines a structure to record history dependent parameter in static analysis.
subroutine, public read_temperature_result(hecmesh, nstep, sstep, rtype, interval, factor, ctime, temp, temp_bak)
Read in temperature distribution from external file.
This subroutine read in used-defined loading tangent.
subroutine uloading(cstep, factor, exforce)
This subroutine take consider of user-defined external loading.
FSTR INNER CONTROL PARAMETERS (fstrPARAM)