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