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