FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
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
13 public :: fstr_get_residual
14 public :: fstr_get_norm_contact
17
18 private :: fstr_update_ndforce_mpc
19
20contains
21
22 !C---------------------------------------------------------------------*
23 subroutine fstr_update_ndforce(cstep,hecMESH,hecMAT,fstrSOLID,conMAT)
24 !C---------------------------------------------------------------------*
30 use m_fstr
31 use muload
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
279end module m_fstr_residual
This module provides functions of reconstructing.
Definition: hecmw.f90:6
This module provides function to calcualte residual of nodal force.
real(kind=kreal) function, public fstr_get_x_norm_contact(hecmat, fstrmat, 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.
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)
subroutine, public fstr_update_ndforce_spc(cstep, hecmesh, fstrsolid, b)
subroutine, public fstr_update_ndforce(cstep, hecmesh, hecmat, fstrsolid, conmat)
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
logical function fstr_isboundaryactive(fstrsolid, nbc, cstep)
Definition: m_fstr.f90:997
integer(kind=kint) myrank
PARALLEL EXECUTION.
Definition: m_fstr.f90:80
integer(kind=kint) nprocs
Definition: m_fstr.f90:81
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...