26 type(hecmwst_matrix) :: hecMAT
27 type(hecmwst_local_mesh) :: hecMESH
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
44 logical,
save :: isFirst = .true.
46 integer(kind=kint) :: flag_u, ierror
47 integer(kind=kint),
optional :: iter
48 real(kind=kreal) :: f_t, t_t
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)
57 integer(kind=kint) :: n_rot, rid, n_nodes, idof
59 real(kind=kreal) :: tval, normal(3), direc(3), ccoord(3), cdisp(3), cdiff(3)
62 call hecmw_mat_clear_b( hecmat )
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
74 n_rot = fstrsolid%CLOAD_ngrp_rot
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)
83 call table_dyn(hecmesh, fstrsolid, fstrdynamic, ig0, f_t, flag_u)
86 is0= hecmesh%node_group%grp_index(ig-1)+1
87 ie0= hecmesh%node_group%grp_index(ig )
89 if( fstrsolid%CLOAD_ngrp_rotID(ig0) > 0 )
then
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
96 if( ityp>ndof ) ityp = ityp-ndof
97 rinfo%conds(rid)%vec(ityp) = val
102 in = hecmesh%node_group%grp_item(ik)
103 hecmat%B( ndof*(in-1)+ityp ) = hecmat%B( ndof*(in-1)+ityp )+val
109 if( .not. rinfo%conds(rid)%active ) cycle
111 n_nodes = hecmw_ngrp_get_number(hecmesh, rinfo%conds(rid)%torque_ngrp_id)
114 ig = rinfo%conds(rid)%center_ngrp_id
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)
119 ccoord(1:ndof) = ccoord(1:ndof) + cdisp(1:ndof)
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() )
126 normal(1:ndof) = rinfo%conds(rid)%vec(1:ndof)/tval
127 tval = tval/dble(n_nodes)
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 )
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)
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() )
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)
150 do ig0 = 1, fstrsolid%DLOAD_ngrp_tot
151 ig = fstrsolid%DLOAD_ngrp_ID(ig0)
152 ltype = fstrsolid%DLOAD_ngrp_LID(ig0)
154 params(i) = fstrsolid%DLOAD_ngrp_params(i,ig0)
157 fg_surf = (ltype == 100)
159 is0 = hecmesh%surf_group%grp_index(ig-1) + 1
160 ie0 = hecmesh%surf_group%grp_index(ig )
162 is0 = hecmesh%elem_group%grp_index(ig-1) + 1
163 ie0 = hecmesh%elem_group%grp_index(ig )
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)
172 icel = hecmesh%elem_group%grp_item(ik)
173 ic_type = hecmesh%elem_type(icel)
177 if( fstrsolid%elements(icel)%elemact_flag == kelact_inactive ) cycle
180 nn = hecmw_get_max_node(ic_type)
182 is = hecmesh%elem_node_index(icel-1)
184 nodlocal(j) = hecmesh%elem_node_item (is+j)
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) )
191 iwk(ndof*(j-1)+i) = ndof*(nodlocal(j)-1)+i
195 isect = hecmesh%section_ID(icel)
197 rho = fstrsolid%elements(icel)%gausses(1)%pMaterial%variables(m_density)
202 id = hecmesh%section%sect_opt(isect)
205 elseif( id == 1 )
then
207 elseif( id == 2 )
then
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)
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)
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)
235 call table_dyn(hecmesh, fstrsolid, fstrdynamic, ig0, f_t, flag_u)
237 vect(j) = vect(j)*f_t
242 hecmat%B( iwk(j) )=hecmat%B( iwk(j) )+vect(j)
247 if (
present(iter) )
then
249 do i = 1, ndof*hecmesh%n_node
250 unode_tmp(i) = fstrsolid%unode(i)
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)
258 idofe = ityp-idofs*10
259 is0 = hecmesh%node_group%grp_index(ig-1) + 1
260 ie0 = hecmesh%node_group%grp_index(ig )
263 in = hecmesh%node_group%grp_item(ik)
264 do idof = idofs, idofe
265 unode_tmp( ndof*(in-1)+idof ) = rhs
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"
276 is = hecmesh%elem_type_index(itype-1)+1
277 ie = hecmesh%elem_type_index(itype )
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
283 call hecmw_abort(hecmw_comm_get_comm())
288 iis = hecmesh%elem_node_index(icel-1)
290 nodlocal(j) = hecmesh%elem_node_item(iis+j)
293 ecoord(i,j) = hecmesh%node( 3*nodlocal(j)+i-3 )
295 v(i,j) = unode_tmp( ndof*nodlocal(j)+i-ndof )
296 fstrsolid%unode( ndof*nodlocal(j)+i-ndof ) = v(i,j)
298 dv(i,j) = fstrsolid%dunode( ndof*nodlocal(j)+i-ndof )
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 )
309 hecmat%B(ndof*(nodlocal(j)-1)+i) = hecmat%B(ndof*(nodlocal(j)-1)+i)+r(ndof*(j-1)+i)
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 )
325 iis = hecmesh%elem_node_index(icel-1)
327 nodlocal(j) = hecmesh%elem_node_item(iis+j)
331 hecmat%B(ndof*(nodlocal(j)-1)+i) = 0.0d0
345 if( fstrsolid%TEMP_ngrp_tot > 0 )
then
347 if( hecmesh%my_rank .eq. 0 )
then
348 write(
imsg,*)
'stop: THERMAL LOAD is not yet available in dynamic analysis!'
350 call hecmw_abort( hecmw_comm_get_comm())
352 allocate ( temp(hecmesh%n_node) )
354 do ig0= 1, fstrsolid%TEMP_ngrp_tot
355 ig= fstrsolid%TEMP_ngrp_ID(ig0)
356 val=fstrsolid%TEMP_ngrp_val(ig0)
358 is0= hecmesh%node_group%grp_index(ig-1) + 1
359 ie0= hecmesh%node_group%grp_index(ig )
361 in = hecmesh%node_group%grp_item(ik)
371 do itype = 1, hecmesh%n_elem_type
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
379 nn = hecmw_get_max_node(ic_type)
384 if( fstrsolid%elements(icel)%elemact_flag == kelact_inactive ) cycle
387 is= hecmesh%elem_node_index(icel-1)
389 nodlocal(j)=hecmesh%elem_node_item(is+j)
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) )
398 iwk(ndof*(j-1)+i)=ndof*(nodlocal(j)-1)+i
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 )
408 if( ndof .eq. 2 )
then
409 id=hecmesh%section%sect_opt(isect)
412 elseif( id.eq.1)
then
414 elseif( id.eq.2)
then
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) )
425 elseif( ic_type == 361 )
then
426 if( fstrsolid%sections(isect)%elemopt361 ==
kel361fi )
then
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
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 )
444 elseif (ic_type == 341 .or. ic_type == 351 .or. &
445 ic_type == 342 .or. ic_type == 352 .or. ic_type == 362 )
then
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 )
450 elseif ( ic_type == 741 .or. ic_type == 743 .or. ic_type == 731 )
then
452 write(
imsg,*)
'*------------------------', &
453 '-------------------*'
454 write(
imsg,*)
' Thermal loading option ', &
456 write(
imsg,*)
'*------------------------', &
457 '-------------------*'
458 call hecmw_abort( hecmw_comm_get_comm())
463 hecmat%B( iwk(j) ) = hecmat%B( iwk(j) )+vect(j)
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.
integer(kind=kint), parameter kel361bbar
subroutine get_coordsys(cdsys_ID, hecMESH, fstrSOLID, coords)
This subroutine fetch coords defined by local coordinate system.
integer(kind=kint) myrank
PARALLEL EXECUTION.
integer(kind=kint), parameter imsg
integer(kind=kint), parameter kel361fi
integer(kind=kint), parameter kel361ic
real(kind=kreal), pointer ref_temp
REFTEMP.
integer(kind=kint), parameter kel361fbar
This modules just summarizes all modules used in static analysis.
Table of lading step in dynamic analysis.
subroutine table_dyn(hecMESH, fstrSOLID, fstrDYNAMIC, ig0, f_t, flag_u)
This module provides aux functions.
subroutine cross_product(v1, v2, vn)
Data for DYNAMIC ANSLYSIS (fstrDYNAMIC)
FSTR INNER CONTROL PARAMETERS (fstrPARAM)