FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
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 
9 contains
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 
331 end 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_bc_init(hecMESH, hecMAT, fstrSOLID, fstrDYNAMIC)
This subroutine setup initial condition of displacement.
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.
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