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