FrontISTR  5.9.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(cstep, hecMESH, hecMAT, fstrSOLID ,fstrDYNAMIC, fstrPARAM, hecLagMAT, t_curr, iter, conMAT)
14  use m_fstr
15  use m_table_dyn
16  use mcontact
17  use m_utilities
18 
19  implicit none
20  integer, intent(in) :: cstep
21  type(hecmwst_matrix) :: hecMAT
22  type(hecmwst_local_mesh) :: hecMESH
23  type(fstr_solid) :: fstrSOLID
24  type(fstr_dynamic) :: fstrDYNAMIC
25  type(fstr_param) :: fstrPARAM
26  type(hecmwst_matrix_lagrange) :: hecLagMAT
27  real(kind=kreal) :: t_curr
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, grpid
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  real(kind=kreal) :: rotation_factor
42  !
43  ndof = hecmat%NDOF
44  n_rot = fstrsolid%BOUNDARY_ngrp_rot
45  if( n_rot > 0 ) call fstr_rotinfo_init(n_rot, rinfo)
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  grpid = fstrsolid%BOUNDARY_ngrp_GRPID(ig0)
56  if( .not. fstr_isboundaryactive( fstrsolid, grpid, cstep ) ) cycle
57  rhs = fstrsolid%BOUNDARY_ngrp_val(ig0)
58 
59  if( present(iter) ) then
60  if( iter>1 ) then
61  rhs=0.d0
62  else
63  call table_dyn(hecmesh, fstrsolid, fstrdynamic, ig0, t_curr-fstrdynamic%t_delta, f_t1, flag_u)
64  call table_dyn(hecmesh, fstrsolid, fstrdynamic, ig0, t_curr, f_t, flag_u)
65  rhs = rhs * (f_t-f_t1)
66  endif
67  else
68  call table_dyn(hecmesh, fstrsolid, fstrdynamic, ig0, t_curr, 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  enddo
114  enddo
115 
116  enddo
117 
118  !Apply rotational boundary condition
119  do rid = 1, n_rot
120  if( .not. rinfo%conds(rid)%active ) cycle
121  cdiff = 0.d0
122  cdiff0 = 0.d0
123  cddisp = 0.d0
124 
125  if( present(iter) ) then
126  if( iter > 1 ) then
127  rotation_factor = 0.d0 ! No additional rotation for subsequent iterations
128  else
129  rotation_factor = 1.0d0 ! Use the rotation vector as-is (already contains increment)
130  endif
131  else
132  rotation_factor = 1.0d0 ! Use the rotation vector as-is
133  endif
134 
135  if( rotation_factor > 0.d0 ) then
136  ig = rinfo%conds(rid)%center_ngrp_id
137  do idof = 1, ndof
138  ccoord(idof) = hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, hecmesh%node)
139  cdisp(idof) = hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, fstrsolid%unode)
140  cddisp(idof) = hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, hecmat%B)
141  enddo
142  ccoord(1:ndof) = ccoord(1:ndof) + cdisp(1:ndof)
143  endif
144 
145  ig = rinfo%conds(rid)%torque_ngrp_id
146  is0 = hecmesh%node_group%grp_index(ig-1) + 1
147  ie0 = hecmesh%node_group%grp_index(ig )
148  do ik = is0, ie0
149  in = hecmesh%node_group%grp_item(ik)
150  if( rotation_factor > 0.d0 ) then
151  cdiff0(1:ndof) = hecmesh%node(ndof*(in-1)+1:ndof*in)+fstrsolid%unode(ndof*(in-1)+1:ndof*in)-ccoord(1:ndof)
152  cdiff(1:ndof) = cdiff0(1:ndof)
153  call rotate_3dvector_by_rodrigues_formula(rinfo%conds(rid)%vec(1:ndof),cdiff(1:ndof))
154  endif
155  do idof = 1, ndof
156  rhs = cdiff(idof)-cdiff0(idof)+cddisp(idof)
157  if(present(conmat)) then
158  call hecmw_mat_ass_bc(hecmat, in, idof, rhs, conmat)
159  else
160  call hecmw_mat_ass_bc(hecmat, in, idof, rhs)
161  endif
162  if( fstr_is_contact_active() .and. fstrparam%solution_type == kststatic &
163  .and. fstrparam%contact_algo == kcaslagrange ) then
164  if(present(conmat)) then
165  call hecmw_mat_ass_bc_contactlag(conmat,heclagmat,in,idof,rhs)
166  else
167  call hecmw_mat_ass_bc_contactlag(hecmat,heclagmat,in,idof,rhs)
168  endif
169  endif
170  enddo
171  enddo
172  enddo
173  !C
174  !C-- end of implicit dynamic analysis
175  !C
176 
177  !C=============================C
178  !C-- explicit dynamic analysis
179  !C=============================C
180  else if( fstrdynamic%idx_eqa == 11 ) then
181  !C
182  ndof = hecmat%NDOF
183  do ig0 = 1, fstrsolid%BOUNDARY_ngrp_tot
184  ig = fstrsolid%BOUNDARY_ngrp_ID(ig0)
185  rhs = fstrsolid%BOUNDARY_ngrp_val(ig0)
186 
187  call table_dyn(hecmesh, fstrsolid, fstrdynamic, ig0, t_curr, f_t, flag_u)
188  rhs = rhs * f_t
189 
190  ityp = fstrsolid%BOUNDARY_ngrp_type(ig0)
191 
192  is0 = hecmesh%node_group%grp_index(ig-1) + 1
193  ie0 = hecmesh%node_group%grp_index(ig )
194  idofs = ityp/10
195  idofe = ityp - idofs*10
196 
197  do ik = is0, ie0
198  in = hecmesh%node_group%grp_item(ik)
199 
200  do idof = idofs, idofe
201  hecmat%B (ndof*in-(ndof-idof)) = rhs
202  fstrdynamic%VEC1(ndof*in-(ndof-idof)) = 1.0d0
203  end do
204  enddo
205  enddo
206  !C
207  !C-- end of explicit dynamic analysis
208  !C
209  end if
210 
211  if( n_rot > 0 ) call fstr_rotinfo_finalize(rinfo)
212 
213  end subroutine dynamic_mat_ass_bc
214 
215 
216  !C***
218  !C***
219  subroutine dynamic_bc_init(hecMESH, hecMAT, fstrSOLID ,fstrDYNAMIC, t_curr)
220  use m_fstr
221  use m_table_dyn
222 
223  implicit none
224  type(hecmwst_matrix) :: hecmat
225  type(hecmwst_local_mesh) :: hecMESH
226  type(fstr_solid) :: fstrSOLID
227  type(fstr_dynamic) :: fstrDYNAMIC
228  real(kind=kreal) :: t_curr
229 
230  integer(kind=kint) :: NDOF, ig0, ig, ityp, iS0, iE0, ik, in, idofS, idofE, idof
231  integer(kind=kint) :: flag_u, grpid
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  grpid = fstrsolid%BOUNDARY_ngrp_GRPID(ig0)
241  if( .not. fstr_isboundaryactive( fstrsolid, grpid, 1 ) ) cycle
242 
243  call table_dyn(hecmesh, fstrsolid, fstrdynamic, ig0, t_curr, f_t, flag_u)
244  rhs = rhs * f_t
245 
246  ityp = fstrsolid%BOUNDARY_ngrp_type(ig0)
247 
248  is0 = hecmesh%node_group%grp_index(ig-1) + 1
249  ie0 = hecmesh%node_group%grp_index(ig )
250  idofs = ityp/10
251  idofe = ityp - idofs*10
252 
253  do ik = is0, ie0
254  in = hecmesh%node_group%grp_item(ik)
255 
256  do idof = idofs, idofe
257  fstrdynamic%DISP(ndof*in-(ndof-idof),1) = rhs
258  end do
259  enddo
260  enddo
261 
262  return
263  end subroutine dynamic_bc_init
264 
266  subroutine dynamic_explicit_ass_bc(hecMESH, hecMAT, fstrSOLID ,fstrDYNAMIC, t_curr, iter)
267  use m_fstr
268  use m_table_dyn
269  use mcontact
270  use m_utilities
271 
272  implicit none
273  type(hecmwst_matrix) :: hecmat
274  type(hecmwst_local_mesh) :: hecMESH
275  type(fstr_solid) :: fstrSOLID
276  type(fstr_dynamic) :: fstrDYNAMIC
277  real(kind=kreal) :: t_curr
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 
295  flag_u = 1
296 
297  !C
298  ndof = hecmat%NDOF
299  do ig0 = 1, fstrsolid%BOUNDARY_ngrp_tot
300  ig = fstrsolid%BOUNDARY_ngrp_ID(ig0)
301  rhs = fstrsolid%BOUNDARY_ngrp_val(ig0)
302 
303  call table_dyn(hecmesh, fstrsolid, fstrdynamic, ig0, t_curr, f_t, flag_u)
304  rhs = rhs * f_t
305 
306  ityp = fstrsolid%BOUNDARY_ngrp_type(ig0)
307 
308  is0 = hecmesh%node_group%grp_index(ig-1) + 1
309  ie0 = hecmesh%node_group%grp_index(ig )
310  idofs = ityp/10
311  idofe = ityp - idofs*10
312 
313  do ik = is0, ie0
314  in = hecmesh%node_group%grp_item(ik)
315 
316  do idof = idofs, idofe
317  hecmat%B(ndof*in-(ndof-idof)) = rhs*fstrdynamic%VEC1(ndof*in-(ndof-idof))
318  ! fstrDYNAMIC%VEC1(NDOF*in-(NDOF-idof)) = 1.0d0
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_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_mat_ass_bc
subroutine dynamic_mat_ass_bc(cstep, hecMESH, hecMAT, fstrSOLID, fstrDYNAMIC, fstrPARAM, hecLagMAT, t_curr, 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, t_curr, iter)
This subroutine setup disp boundary condition.
Definition: dynamic_mat_ass_bc.f90:267
m_dynamic_mat_ass_bc::dynamic_bc_init
subroutine dynamic_bc_init(hecMESH, hecMAT, fstrSOLID, fstrDYNAMIC, t_curr)
This subroutine setup initial condition of displacement.
Definition: dynamic_mat_ass_bc.f90:220
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:61
m_fstr::kststatic
integer(kind=kint), parameter kststatic
Definition: m_fstr.F90:39
m_fstr::fstr_isboundaryactive
logical function fstr_isboundaryactive(fstrSOLID, nbc, cstep)
Definition: m_fstr.F90:1053
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, t_curr, 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:591
mcontact
Top-level contact analysis module (System level)
Definition: fstr_contact.f90:33
mcontact::fstr_is_contact_active
logical function fstr_is_contact_active()
Definition: fstr_contact.f90:349
m_fstr::fstr_dynamic
Data for DYNAMIC ANSLYSIS (fstrDYNAMIC)
Definition: m_fstr.F90:530
m_fstr::fstr_param
FSTR INNER CONTROL PARAMETERS (fstrPARAM)
Definition: m_fstr.F90:156
m_fstr::fstr_solid
Definition: m_fstr.F90:242