FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
dynamic_mat_ass_bc.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
8
9contains
10
11
13 subroutine dynamic_mat_ass_bc(hecMESH, hecMAT, fstrSOLID ,fstrDYNAMIC, fstrPARAM, fstrMAT, iter, conMAT)
14 use m_fstr
15 use m_table_dyn
18 use mcontact
19 use m_utilities
20
21 implicit none
22 type(hecmwst_matrix) :: hecMAT
23 type(hecmwst_local_mesh) :: hecMESH
24 type(fstr_solid) :: fstrSOLID
25 type(fstr_dynamic) :: fstrDYNAMIC
26 type(fstr_param) :: fstrPARAM
27 type(fstrst_matrix_contact_lagrange) :: fstrMAT
28 integer, optional :: iter
29 type(hecmwst_matrix), optional :: conMAT
30
31 integer(kind=kint) :: ig0, ig, ityp, NDOF, iS0, iE0, ik, in, idofS, idofE, idof
32
33 integer(kind=kint) :: flag_u
34 real(kind=kreal) :: rhs, f_t, f_t1
35
36 !for rotation
37 integer(kind=kint) :: n_rot, rid, n_nodes
38 type(trotinfo) :: rinfo
39 real(kind=kreal) :: theta, normal(3), direc(3), ccoord(3), cdiff(3), cdiff0(3)
40 real(kind=kreal) :: cdisp(3), cddisp(3)
41 !
42 ndof = hecmat%NDOF
43 n_rot = fstrsolid%BOUNDARY_ngrp_rot
44 if( n_rot > 0 ) call fstr_rotinfo_init(n_rot, rinfo)
45 fstrsolid%REACTION = 0.d0
46
47 flag_u = 1
48 !C=============================C
49 !C-- implicit dynamic analysis
50 !C=============================C
51 if( fstrdynamic%idx_eqa == 1 ) then
52
53 do ig0 = 1, fstrsolid%BOUNDARY_ngrp_tot
54 ig = fstrsolid%BOUNDARY_ngrp_ID(ig0)
55 rhs = fstrsolid%BOUNDARY_ngrp_val(ig0)
56
57 if( present(iter) ) then
58 if( iter>1 ) then
59 rhs=0.d0
60 else
61 fstrdynamic%i_step = fstrdynamic%i_step-1
62 fstrdynamic%t_curr = fstrdynamic%t_curr - fstrdynamic%t_delta
63 call table_dyn(hecmesh, fstrsolid, fstrdynamic, ig0, f_t1, flag_u)
64 fstrdynamic%i_step = fstrdynamic%i_step+1
65 fstrdynamic%t_curr = fstrdynamic%t_curr + fstrdynamic%t_delta
66 call table_dyn(hecmesh, fstrsolid, fstrdynamic, ig0, f_t, flag_u)
67 rhs = rhs * (f_t-f_t1)
68 endif
69 else
70 call table_dyn(hecmesh, fstrsolid, fstrdynamic, ig0, f_t, flag_u)
71 rhs = rhs * f_t
72 endif
73
74 ityp = fstrsolid%BOUNDARY_ngrp_type(ig0)
75 idofs = ityp/10
76 idofe = ityp - idofs*10
77
78 is0 = hecmesh%node_group%grp_index(ig-1) + 1
79 ie0 = hecmesh%node_group%grp_index(ig )
80
81 if( fstrsolid%BOUNDARY_ngrp_rotID(ig0) > 0 ) then ! setup rotation information
82 rid = fstrsolid%BOUNDARY_ngrp_rotID(ig0)
83 if( .not. rinfo%conds(rid)%active ) then
84 rinfo%conds(rid)%active = .true.
85 rinfo%conds(rid)%center_ngrp_id = fstrsolid%BOUNDARY_ngrp_centerID(ig0)
86 rinfo%conds(rid)%torque_ngrp_id = ig
87 endif
88 do idof=idofs,idofe
89 if( idof>ndof ) then
90 rinfo%conds(rid)%vec(idof-ndof) = rhs
91 else
92 rinfo%conds(rid)%vec(idof) = rhs
93 endif
94 enddo
95 cycle
96 endif
97
98 do ik = is0, ie0
99 in = hecmesh%node_group%grp_item(ik)
100
101 do idof = idofs, idofe
102 if(present(conmat)) then
103 call hecmw_mat_ass_bc(hecmat, in, idof, rhs, conmat)
104 else
105 call hecmw_mat_ass_bc(hecmat, in, idof, rhs)
106 endif
107 if( fstr_is_contact_active() .and. fstrparam%contact_algo == kcaslagrange &
108 .and. fstrparam%nlgeom .and. fstrdynamic%idx_resp == 1 ) then
109 if(present(conmat)) then
110 call fstr_mat_ass_bc_contact(conmat,fstrmat,in,idof,rhs)
111 else
112 call fstr_mat_ass_bc_contact(hecmat,fstrmat,in,idof,rhs)
113 endif
114 endif
115
116 !for output reaction force
117 fstrsolid%REACTION(ndof*(in-1)+idof) = fstrsolid%QFORCE(ndof*(in-1)+idof)
118 enddo
119 enddo
120
121 enddo
122
123 !Apply rotational boundary condition
124 do rid = 1, n_rot
125 if( .not. rinfo%conds(rid)%active ) cycle
126 cdiff = 0.d0
127 cdiff0 = 0.d0
128 cddisp = 0.d0
129
130 if( f_t > 0.d0 ) then
131 ig = rinfo%conds(rid)%center_ngrp_id
132 do idof = 1, ndof
133 ccoord(idof) = hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, hecmesh%node)
134 cdisp(idof) = hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, fstrsolid%unode)
135 cddisp(idof) = hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, hecmat%B)
136 enddo
137 ccoord(1:ndof) = ccoord(1:ndof) + cdisp(1:ndof)
138 endif
139
140 ig = rinfo%conds(rid)%torque_ngrp_id
141 is0 = hecmesh%node_group%grp_index(ig-1) + 1
142 ie0 = hecmesh%node_group%grp_index(ig )
143 do ik = is0, ie0
144 in = hecmesh%node_group%grp_item(ik)
145 if( f_t > 0.d0 ) then
146 cdiff0(1:ndof) = hecmesh%node(ndof*(in-1)+1:ndof*in)+fstrsolid%unode(ndof*(in-1)+1:ndof*in)-ccoord(1:ndof)
147 cdiff(1:ndof) = cdiff0(1:ndof)
148 call rotate_3dvector_by_rodrigues_formula(rinfo%conds(rid)%vec(1:ndof),cdiff(1:ndof))
149 endif
150 do idof = 1, ndof
151 rhs = cdiff(idof)-cdiff0(idof)+cddisp(idof)
152 if(present(conmat)) then
153 call hecmw_mat_ass_bc(hecmat, in, idof, rhs, conmat)
154 else
155 call hecmw_mat_ass_bc(hecmat, in, idof, rhs)
156 endif
157 if( fstr_is_contact_active() .and. fstrparam%solution_type == kststatic &
158 .and. fstrparam%contact_algo == kcaslagrange ) then
159 if(present(conmat)) then
160 call fstr_mat_ass_bc_contact(conmat,fstrmat,in,idof,rhs)
161 else
162 call fstr_mat_ass_bc_contact(hecmat,fstrmat,in,idof,rhs)
163 endif
164 endif
165
166 !for output reaction force
167 fstrsolid%REACTION(ndof*(in-1)+idof) = fstrsolid%QFORCE(ndof*(in-1)+idof)
168 enddo
169 enddo
170 enddo
171 !C
172 !C-- end of implicit dynamic analysis
173 !C
174
175 !C=============================C
176 !C-- explicit dynamic analysis
177 !C=============================C
178 else if( fstrdynamic%idx_eqa == 11 ) then
179 !C
180 ndof = hecmat%NDOF
181 do ig0 = 1, fstrsolid%BOUNDARY_ngrp_tot
182 ig = fstrsolid%BOUNDARY_ngrp_ID(ig0)
183 rhs = fstrsolid%BOUNDARY_ngrp_val(ig0)
184
185 call table_dyn(hecmesh, fstrsolid, fstrdynamic, ig0, f_t, flag_u)
186 rhs = rhs * f_t
187
188 ityp = fstrsolid%BOUNDARY_ngrp_type(ig0)
189
190 is0 = hecmesh%node_group%grp_index(ig-1) + 1
191 ie0 = hecmesh%node_group%grp_index(ig )
192 idofs = ityp/10
193 idofe = ityp - idofs*10
194
195 do ik = is0, ie0
196 in = hecmesh%node_group%grp_item(ik)
197
198 do idof = idofs, idofe
199 hecmat%B (ndof*in-(ndof-idof)) = rhs
200 fstrdynamic%VEC1(ndof*in-(ndof-idof)) = 1.0d0
201
202 !for output reaction force
203 fstrsolid%REACTION(ndof*(in-1)+idof) = fstrsolid%QFORCE(ndof*(in-1)+idof)
204 end do
205 enddo
206 enddo
207 !C
208 !C-- end of explicit dynamic analysis
209 !C
210 end if
211
212 if( n_rot > 0 ) call fstr_rotinfo_finalize(rinfo)
213
214 end subroutine dynamic_mat_ass_bc
215
216
217 !C***
219 !C***
220 subroutine dynamic_bc_init(hecMESH, hecMAT, fstrSOLID ,fstrDYNAMIC)
221 use m_fstr
222 use m_table_dyn
223
224 implicit none
225 type(hecmwst_matrix) :: hecmat
226 type(hecmwst_local_mesh) :: hecMESH
227 type(fstr_solid) :: fstrSOLID
228 type(fstr_dynamic) :: fstrDYNAMIC
229
230 integer(kind=kint) :: NDOF, ig0, ig, ityp, iS0, iE0, ik, in, idofS, idofE, idof
231 integer(kind=kint) :: flag_u
232 real(kind=kreal) :: rhs, f_t
233
234 flag_u = 1
235 ndof = hecmat%NDOF
236
237 do ig0 = 1, fstrsolid%BOUNDARY_ngrp_tot
238 ig = fstrsolid%BOUNDARY_ngrp_ID(ig0)
239 rhs = fstrsolid%BOUNDARY_ngrp_val(ig0)
240
241
242 call table_dyn(hecmesh, fstrsolid, fstrdynamic, ig0, f_t, flag_u)
243 rhs = rhs * f_t
244
245 ityp = fstrsolid%BOUNDARY_ngrp_type(ig0)
246
247 is0 = hecmesh%node_group%grp_index(ig-1) + 1
248 ie0 = hecmesh%node_group%grp_index(ig )
249 idofs = ityp/10
250 idofe = ityp - idofs*10
251
252 do ik = is0, ie0
253 in = hecmesh%node_group%grp_item(ik)
254
255 do idof = idofs, idofe
256 fstrdynamic%DISP(ndof*in-(ndof-idof),1) = rhs
257 end do
258 enddo
259 enddo
260
261 return
262 end subroutine dynamic_bc_init
263
265 subroutine dynamic_explicit_ass_bc(hecMESH, hecMAT, fstrSOLID ,fstrDYNAMIC, iter)
266 use m_fstr
267 use m_table_dyn
270 use mcontact
271 use m_utilities
272
273 implicit none
274 type(hecmwst_matrix) :: hecmat
275 type(hecmwst_local_mesh) :: hecMESH
276 type(fstr_solid) :: fstrSOLID
277 type(fstr_dynamic) :: fstrDYNAMIC
278 integer, optional :: iter
279
280 integer(kind=kint) :: ig0, ig, ityp, NDOF, iS0, iE0, ik, in, idofS, idofE, idof
281
282 integer(kind=kint) :: flag_u
283 real(kind=kreal) :: rhs, f_t, f_t1
284
285 !for rotation
286 integer(kind=kint) :: n_rot, rid, n_nodes
287 type(trotinfo) :: rinfo
288 real(kind=kreal) :: theta, normal(3), direc(3), ccoord(3), cdiff(3), cdiff0(3)
289 real(kind=kreal) :: cdisp(3), cddisp(3)
290 !
291 ndof = hecmat%NDOF
292 n_rot = fstrsolid%BOUNDARY_ngrp_rot
293 if( n_rot > 0 ) call fstr_rotinfo_init(n_rot, rinfo)
294 fstrsolid%REACTION = 0.d0
295
296 flag_u = 1
297
298 !C
299 ndof = hecmat%NDOF
300 do ig0 = 1, fstrsolid%BOUNDARY_ngrp_tot
301 ig = fstrsolid%BOUNDARY_ngrp_ID(ig0)
302 rhs = fstrsolid%BOUNDARY_ngrp_val(ig0)
303
304 call table_dyn(hecmesh, fstrsolid, fstrdynamic, ig0, f_t, flag_u)
305 rhs = rhs * f_t
306
307 ityp = fstrsolid%BOUNDARY_ngrp_type(ig0)
308
309 is0 = hecmesh%node_group%grp_index(ig-1) + 1
310 ie0 = hecmesh%node_group%grp_index(ig )
311 idofs = ityp/10
312 idofe = ityp - idofs*10
313
314 do ik = is0, ie0
315 in = hecmesh%node_group%grp_item(ik)
316
317 do idof = idofs, idofe
318 hecmat%B(ndof*in-(ndof-idof)) = rhs*fstrdynamic%VEC1(ndof*in-(ndof-idof))
319 ! fstrDYNAMIC%VEC1(NDOF*in-(NDOF-idof)) = 1.0d0
320
321 !for output reaction force
322 fstrsolid%REACTION(ndof*(in-1)+idof) = fstrsolid%QFORCE(ndof*(in-1)+idof)
323 end do
324 enddo
325 enddo
326
327 if( n_rot > 0 ) call fstr_rotinfo_finalize(rinfo)
328
329 end subroutine dynamic_explicit_ass_bc
330
331end module m_dynamic_mat_ass_bc
This module provides functions of reconstructing.
This module provides functions: 1) obtain contact stiffness matrix of each contact pair and assemble ...
subroutine, public fstr_mat_ass_bc_contact(hecmat, fstrmat, inode, idof, rhs)
Modify Lagrange multiplier-related part of stiffness matrix and right-hand side vector for dealing wi...
This module contains functions to set displacement boundary condition in dynamic analysis.
subroutine dynamic_mat_ass_bc(hecmesh, hecmat, fstrsolid, fstrdynamic, fstrparam, fstrmat, iter, conmat)
This subroutine setup disp bundary condition.
subroutine dynamic_explicit_ass_bc(hecmesh, hecmat, fstrsolid, fstrdynamic, iter)
This subroutine setup disp boundary condition.
subroutine dynamic_bc_init(hecmesh, hecmat, fstrsolid, fstrdynamic)
This subroutine setup initial condition of displacement.
This module defined coomon data and basic structures for analysis.
Definition: m_fstr.f90:15
integer(kind=kint), parameter kcaslagrange
contact analysis algorithm
Definition: m_fstr.f90:53
integer(kind=kint), parameter kststatic
Definition: m_fstr.f90:37
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 rotate_3dvector_by_rodrigues_formula(r, v)
Definition: utilities.f90:515
This module provides functions to calcualte contact stiff matrix.
Definition: fstr_contact.f90:6
logical function fstr_is_contact_active()
Structure for Lagrange multiplier-related part of stiffness matrix (Lagrange multiplier-related matri...
Data for DYNAMIC ANSLYSIS (fstrDYNAMIC)
Definition: m_fstr.f90:473
FSTR INNER CONTROL PARAMETERS (fstrPARAM)
Definition: m_fstr.f90:138