FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
fstr_Residual.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  use hecmw
9  implicit none
10 
11  public :: fstr_update_ndforce
12  public :: fstr_update_ndforce_spc
13  public :: fstr_get_residual
14  public :: fstr_get_norm_contact
16  public :: fstr_get_x_norm_contact
17 
18  private :: fstr_update_ndforce_mpc
19 
20 contains
21 
22  !C---------------------------------------------------------------------*
23  subroutine fstr_update_ndforce(cstep,hecMESH,hecMAT,fstrSOLID,conMAT)
24  !C---------------------------------------------------------------------*
30  use m_fstr
31  use muload
32  use m_fstr_spring
33  integer(kind=kint), intent(in) :: cstep
34  type(hecmwst_local_mesh), intent(in) :: hecmesh
35  type(hecmwst_matrix), intent(inout) :: hecmat
36  type(fstr_solid), intent(inout) :: fstrsolid
37  type(hecmwst_matrix), intent(inout), optional :: conmat
38  ! Local variables
39  integer(kind=kint) :: ndof, idof
40  real(kind=kreal) :: factor
41 
42  factor = fstrsolid%factor(2)
43 
44  ! Set residual load
45  do idof=1, hecmesh%n_node* hecmesh%n_dof
46  hecmat%B(idof)=fstrsolid%GL(idof)-fstrsolid%QFORCE(idof)
47  end do
48  ndof = hecmat%NDOF
49 
50  call fstr_update_ndforce_spring( cstep, hecmesh, fstrsolid, hecmat%B )
51 
52  ! Consider Uload
53  call uresidual( cstep, factor, hecmat%B )
54 
55  ! Consider EQUATION condition
56  call fstr_update_ndforce_mpc( hecmesh, hecmat%B )
57 
58  ! Consider SPC condition
59  call fstr_update_ndforce_spc( cstep, hecmesh, fstrsolid, hecmat%B )
60  if(present(conmat)) call fstr_update_ndforce_spc( cstep, hecmesh, fstrsolid, conmat%B )
61 
62  !
63  if( ndof==3 ) then
64  call hecmw_update_3_r(hecmesh,hecmat%B,hecmesh%n_node)
65  else if( ndof==2 ) then
66  call hecmw_update_2_r(hecmesh,hecmat%B,hecmesh%n_node)
67  else if( ndof==6 ) then
68  call hecmw_update_m_r(hecmesh,hecmat%B,hecmesh%n_node,6)
69  endif
70 
71  end subroutine fstr_update_ndforce
72 
73  subroutine fstr_update_ndforce_mpc( hecMESH, B )
74  use m_fstr
75  type(hecmwst_local_mesh), intent(in) :: hecmesh
76  real(kind=kreal), intent(inout) :: b(:)
77  ! Local variables
78  integer(kind=kint) ndof, ig0, is0, ie0, ik, in, idof
79  real(kind=kreal) :: rhs, lambda
80 
81  ndof = hecmesh%n_dof
82  outer: do ig0=1,hecmesh%mpc%n_mpc
83  is0= hecmesh%mpc%mpc_index(ig0-1)+1
84  ie0= hecmesh%mpc%mpc_index(ig0)
85  do ik= is0, ie0
86  if (hecmesh%mpc%mpc_dof(ik) > ndof) cycle outer
87  enddo
88  ! Suppose the lagrange multiplier= first dof of first node
89  in = hecmesh%mpc%mpc_item(is0)
90  idof = hecmesh%mpc%mpc_dof(is0)
91  rhs = hecmesh%mpc%mpc_val(is0)
92  lambda = b(ndof*(in-1)+idof)/rhs
93  ! update nodal residual
94  do ik= is0, ie0
95  in = hecmesh%mpc%mpc_item(ik)
96  idof = hecmesh%mpc%mpc_dof(ik)
97  rhs = hecmesh%mpc%mpc_val(ik)
98  b(ndof*(in-1)+idof) = b(ndof*(in-1)+idof) - rhs*lambda
99  enddo
100  enddo outer
101  end subroutine fstr_update_ndforce_mpc
102 
103  subroutine fstr_update_ndforce_spc( cstep, hecMESH, fstrSOLID, B )
104  use m_fstr
105  integer(kind=kint), intent(in) :: cstep
106  type(hecmwst_local_mesh), intent(in) :: hecmesh
107  type(fstr_solid), intent(in) :: fstrsolid
108  real(kind=kreal), intent(inout) :: b(:)
109  ! Local variables
110  integer(kind=kint) ndof, ig0, ig, ityp, is0, ie0, ik, in, idof1, idof2, idof
111  integer(kind=kint) :: grpid
112  real(kind=kreal) :: rhs
113 
114  ndof = hecmesh%n_dof
115  fstrsolid%REACTION = 0.d0
116 
117  do ig0= 1, fstrsolid%BOUNDARY_ngrp_tot
118  grpid = fstrsolid%BOUNDARY_ngrp_GRPID(ig0)
119  if( .not. fstr_isboundaryactive( fstrsolid, grpid, cstep ) ) cycle
120  ig= fstrsolid%BOUNDARY_ngrp_ID(ig0)
121  rhs= fstrsolid%BOUNDARY_ngrp_val(ig0)
122  ityp= fstrsolid%BOUNDARY_ngrp_type(ig0)
123  is0= hecmesh%node_group%grp_index(ig-1) + 1
124  ie0= hecmesh%node_group%grp_index(ig )
125  do ik= is0, ie0
126  in = hecmesh%node_group%grp_item(ik)
127  idof1 = ityp/10
128  idof2 = ityp - idof1*10
129  if( fstrsolid%BOUNDARY_ngrp_rotID(ig0) > 0 ) then
130  idof1 = 1
131  idof2 = ndof
132  end if
133  do idof=idof1,idof2
134  b( ndof*(in-1) + idof ) = 0.d0
135  !for output reaction force
136  fstrsolid%REACTION(ndof*(in-1)+idof) = fstrsolid%QFORCE(ndof*(in-1)+idof)
137  enddo
138  enddo
139  enddo
140  end subroutine fstr_update_ndforce_spc
141 
143  real(kind=kreal) function fstr_get_residual( force, hecMESH )
144  use m_fstr
145  real(kind=kreal), intent(in) :: force(:)
146  type(hecmwst_local_mesh), intent(in) :: hecmesh
147  integer :: ndof
148  ndof = hecmesh%n_dof
149  call hecmw_innerproduct_r(hecmesh,ndof,force,force,fstr_get_residual)
150  end function
151 
153  real(kind=kreal) function fstr_get_energy( force, displacement, hecMESH )
154  use m_fstr
155  real(kind=kreal), intent(in) :: force(:), displacement(:)
156  type(hecmwst_local_mesh), intent(in) :: hecmesh
157  integer :: ndof
158  ndof = hecmesh%n_dof
159  call hecmw_innerproduct_r(hecmesh, ndof, force, displacement, fstr_get_energy)
160  end function
161 
163  real(kind=kreal) function fstr_get_norm_contact(flag,hecMESH,hecMAT,fstrSOLID,fstrMAT)
164  use m_fstr
166  type(hecmwst_local_mesh), intent(in) :: hecmesh
167  type(hecmwst_matrix), intent(in) :: hecmat
168  type(fstr_solid), intent(in) :: fstrsolid
169  type(fstrst_matrix_contact_lagrange), intent(in) :: fstrmat
170  character(len=13) :: flag
171  real(kind=kreal) :: tmp1, tmp2, bi
172  integer :: i, i0, ndof
173  if( flag=='residualForce' )then
174  ndof = hecmesh%n_dof
175  call hecmw_innerproduct_r(hecmesh,ndof,hecmat%B,hecmat%B,tmp1)
176  tmp2 = 0.0d0
177  i0 = hecmesh%n_node*ndof
178  do i=1,fstrmat%num_lagrange
179  bi = hecmat%B(i0+i)
180  tmp2 = tmp2 + bi*bi
181  enddo
182  call hecmw_allreduce_r1(hecmesh,tmp2,hecmw_sum)
183  fstr_get_norm_contact = tmp1 + tmp2
184  elseif( flag==' force' )then
185  call hecmw_innerproduct_r(hecmesh,ndof,fstrsolid%QFORCE,fstrsolid%QFORCE,fstr_get_norm_contact)
186  endif
187  end function
188 
189  !
190  function fstr_get_norm_para_contact(hecMAT,fstrMAT,conMAT,hecMESH) result(rhsB)
191  use m_fstr
193  implicit none
194  type(hecmwst_matrix), intent(in) :: hecmat
195  type(fstrst_matrix_contact_lagrange), intent(in) :: fstrmat
196  type(hecmwst_matrix), intent(in) :: conmat
197  type(hecmwst_local_mesh), intent(in) :: hecmesh
198  !
199  real(kind=kreal) :: rhsb
200  integer(kind=kint) :: i,j,n,i0,ndof,n_loc,nndof
201  integer(kind=kint) :: offset, pid, lid
202  integer(kind=kint), allocatable :: displs(:)
203  real(kind=kreal), allocatable :: rhs_con_all(:), rhs_con(:)
204  !
205  ndof = hecmat%ndof
206  nndof = hecmat%N*ndof
207  n_loc = nndof + fstrmat%num_lagrange
208  allocate(displs(0:nprocs))
209  displs(:) = 0
210  displs(myrank+1) = n_loc
211  call hecmw_allreduce_i(hecmesh, displs, nprocs+1, hecmw_sum)
212  do i=1,nprocs
213  displs(i) = displs(i-1) + displs(i)
214  end do
215  offset = displs(myrank)
216  n = displs(nprocs)
217  allocate(rhs_con_all(n))
218  do i=1,n
219  rhs_con_all(i) = 0.0d0
220  end do
221  do i= hecmat%N+1,hecmat%NP
222  ! i0 = getExternalGlobalIndex(i,ndof,hecMAT%N)
223  pid = hecmesh%node_ID(i*2)
224  lid = hecmesh%node_ID(i*2-1)
225  i0 = displs(pid) + (lid-1)*ndof
226  if((i0 < 0.or.i0 > n)) then
227  ! print *,myrank,'ext dummy',i,i0/ndof
228  do j=1,ndof
229  if(conmat%b((i-1)*ndof+j) /= 0.0d0) then
230  print *,myrank,'i0',i,'conMAT%b',conmat%b((i-1)*ndof+j)
231  stop
232  endif
233  enddo
234  else
235  rhs_con_all(i0+1:i0+ndof) = conmat%b((i-1)*ndof+1:i*ndof)
236  endif
237  enddo
238  deallocate(displs)
239  call hecmw_allreduce_r(hecmesh, rhs_con_all, n, hecmw_sum)
240  allocate(rhs_con(n_loc))
241  do i=1,nndof
242  rhs_con(i) = rhs_con_all(offset+i) + hecmat%B(i) + conmat%B(i)
243  end do
244  deallocate(rhs_con_all)
245  i0 = hecmat%NP*ndof
246  do i=1,fstrmat%num_lagrange
247  rhs_con(nndof+i) = conmat%B(i0+i)
248  enddo
249  rhsb = 0.d0
250  rhsb = dot_product(rhs_con, rhs_con)
251  call hecmw_allreduce_r1(hecmesh, rhsb, hecmw_sum)
252  deallocate(rhs_con)
253 
254  end function fstr_get_norm_para_contact
255 
256  function fstr_get_x_norm_contact(hecMAT,fstrMAT,hecMESH) result(rhsX)
257  use m_fstr
259  implicit none
260  type(hecmwst_matrix), intent(in) :: hecmat
261  type(fstrst_matrix_contact_lagrange), intent(in) :: fstrmat
262  type(hecmwst_local_mesh), intent(in) :: hecmesh
263  real(kind=kreal) :: rhsx
264  integer(kind=kint) :: nndof, npndof, i
265 
266  nndof = hecmat%N * hecmat%NDOF
267  npndof = hecmat%NP * hecmat%NDOF
268  rhsx = 0.d0
269  do i=1,nndof
270  rhsx = rhsx + hecmat%X(i) ** 2
271  end do
272  do i=1,fstrmat%num_lagrange
273  rhsx = rhsx + hecmat%X(npndof+i) ** 2
274  end do
275  call hecmw_allreduce_r1(hecmesh, rhsx, hecmw_sum)
276 
277  end function fstr_get_x_norm_contact
278 
279 end module m_fstr_residual
This module provides functions of reconstructing.
Definition: hecmw.f90:6
This module provides function to calcualte residual of nodal force.
subroutine, public fstr_update_ndforce_spc(cstep, hecMESH, fstrSOLID, B)
subroutine, public fstr_update_ndforce(cstep, hecMESH, hecMAT, fstrSOLID, conMAT)
real(kind=kreal) function, public fstr_get_x_norm_contact(hecMAT, fstrMAT, hecMESH)
real(kind=kreal) function, public fstr_get_residual(force, hecMESH)
Calculate magnitude of a real vector.
real(kind=kreal) function, public fstr_get_norm_para_contact(hecMAT, fstrMAT, conMAT, hecMESH)
real(kind=kreal) function fstr_get_energy(force, displacement, hecMESH)
Calculate magnitude of a real vector.
real(kind=kreal) function, public fstr_get_norm_contact(flag, hecMESH, hecMAT, fstrSOLID, fstrMAT)
Calculate square norm.
This module provides functions to deal with spring force.
Definition: fstr_Spring.f90:7
subroutine fstr_update_ndforce_spring(cstep, hecMESH, fstrSOLID, B)
Definition: fstr_Spring.f90:46
This module defined coomon data and basic structures for analysis.
Definition: m_fstr.f90:15
integer(kind=kint) myrank
PARALLEL EXECUTION.
Definition: m_fstr.f90:80
integer(kind=kint) nprocs
Definition: m_fstr.f90:81
logical function fstr_isboundaryactive(fstrSOLID, nbc, cstep)
Definition: m_fstr.f90:997
This subroutine read in used-defined loading tangent.
Definition: uload.f90:7
subroutine uresidual(cstep, factor, residual)
This subroutine take consider of user-defined external loading.
Definition: uload.f90:39
Structure for Lagrange multiplier-related part of stiffness matrix (Lagrange multiplier-related matri...