FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
fstr_StiffMatrix.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 m_fstr
9  implicit none
10 
11  private
12  public :: fstr_stiffmatrix
13 
14 contains
15 
16  !---------------------------------------------------------------------*
18  subroutine fstr_stiffmatrix( hecMESH, hecMAT, fstrSOLID, time, tincr)
19  !---------------------------------------------------------------------*
20  use m_static_lib
21  use mmechgauss
22 
23  type (hecmwst_local_mesh) :: hecmesh
24  type (hecmwst_matrix) :: hecmat
25  type (fstr_solid) :: fstrsolid
26  real(kind=kreal),intent(in) :: time
27  real(kind=kreal),intent(in) :: tincr
28 
29  type( tmaterial ), pointer :: material
30 
31  real(kind=kreal) :: stiffness(20*6, 20*6)
32  integer(kind=kint) :: nodlocal(20)
33  real(kind=kreal) :: tt(20), ecoord(3,20)
34  real(kind=kreal) :: thick
35  integer(kind=kint) :: ndof, itype, is, ie, ic_type, nn, icel, iis, i, j
36  real(kind=kreal) :: u(6,20), du(6,20), coords(3,3), u_prev(6,20)
37  integer :: isect, ihead, cdsys_id
38 
39  ! ----- initialize
40  call hecmw_mat_clear( hecmat )
41 
42  ndof = hecmat%NDOF
43  do itype= 1, hecmesh%n_elem_type
44  is= hecmesh%elem_type_index(itype-1) + 1
45  ie= hecmesh%elem_type_index(itype )
46  ic_type= hecmesh%elem_type_item(itype)
47  ! ----- Ignore link and patch elements
48  if (hecmw_is_etype_link(ic_type)) cycle
49  if (hecmw_is_etype_patch(ic_type)) cycle
50  ! ----- Set number of nodes
51  nn = hecmw_get_max_node(ic_type)
52 
53  ! ----- element loop
54  !$omp parallel default(none), &
55  !$omp& private(icel,iiS,j,nodLOCAL,i,ecoord,du,u,u_prev,tt,cdsys_ID,coords, &
56  !$omp& material,thick,stiffness,isect,ihead), &
57  !$omp& shared(iS,iE,hecMESH,nn,ndof,fstrSOLID,ic_type,hecMAT,time,tincr)
58  !$omp do
59  do icel= is, ie
60 
61  ! ----- nodal coordinate & displacement
62  iis= hecmesh%elem_node_index(icel-1)
63  do j=1,nn
64  nodlocal(j)= hecmesh%elem_node_item (iis+j)
65  do i=1, 3
66  ecoord(i,j) = hecmesh%node(3*nodlocal(j)+i-3)
67  enddo
68  do i=1,ndof
69  du(i,j) = fstrsolid%dunode(ndof*nodlocal(j)+i-ndof)
70  u(i,j) = fstrsolid%unode(ndof*nodlocal(j)+i-ndof) + du(i,j)
71  u_prev(i,j) = fstrsolid%unode(ndof*nodlocal(j)+i-ndof)
72  enddo
73  if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres >0 ) &
74  tt(j)=fstrsolid%temperature( nodlocal(j) )
75  enddo
76 
77  isect = hecmesh%section_ID(icel)
78  cdsys_id = hecmesh%section%sect_orien_ID(isect)
79  if( cdsys_id > 0 ) call get_coordsys(cdsys_id, hecmesh, fstrsolid, coords)
80 
81  material => fstrsolid%elements(icel)%gausses(1)%pMaterial
82  thick = material%variables(m_thick)
83  if( getspacedimension( ic_type )==2 ) thick =1.d0
84  if( ic_type==241 .or. ic_type==242 .or. ic_type==231 .or. ic_type==232 .or. ic_type==2322) then
85  if( material%nlgeom_flag /= infinitesimal ) call stiffmat_abort( ic_type, 2 )
86  call stf_c2( ic_type,nn,ecoord(1:2,1:nn),fstrsolid%elements(icel)%gausses(:),thick, &
87  stiffness(1:nn*ndof,1:nn*ndof), fstrsolid%elements(icel)%iset, &
88  u(1:2,1:nn) )
89 
90  elseif ( ic_type==301 ) then
91  isect= hecmesh%section_ID(icel)
92  ihead = hecmesh%section%sect_R_index(isect-1)
93  thick = hecmesh%section%sect_R_item(ihead+1)
94  call stf_c1( ic_type,nn,ecoord(:,1:nn),thick,fstrsolid%elements(icel)%gausses(:), &
95  stiffness(1:nn*ndof,1:nn*ndof), u(1:3,1:nn) )
96 
97  elseif ( ic_type==361 ) then
98 
99  if( fstrsolid%sections(isect)%elemopt361 == kel361fi ) then ! full integration element
100  if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres >0 ) then
101  call stf_c3 &
102  ( ic_type, nn, ecoord(:, 1:nn), fstrsolid%elements(icel)%gausses(:), &
103  stiffness(1:nn*ndof, 1:nn*ndof), cdsys_id, coords, time, tincr, u(1:3,1:nn), tt(1:nn) )
104  else
105  call stf_c3 &
106  ( ic_type,nn,ecoord(:, 1:nn),fstrsolid%elements(icel)%gausses(:), &
107  stiffness(1:nn*ndof, 1:nn*ndof), cdsys_id, coords, time, tincr, u(1:3, 1:nn) )
108  endif
109  else if( fstrsolid%sections(isect)%elemopt361 == kel361bbar ) then ! B-bar element
110  if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres >0 ) then
111  call stf_c3d8bbar &
112  ( ic_type, nn, ecoord(:, 1:nn), fstrsolid%elements(icel)%gausses(:), &
113  stiffness(1:nn*ndof,1:nn*ndof), cdsys_id, coords, time, tincr, u(1:3, 1:nn), tt(1:nn) )
114  else
115  call stf_c3d8bbar &
116  ( ic_type, nn, ecoord(:, 1:nn),fstrsolid%elements(icel)%gausses(:), &
117  stiffness(1:nn*ndof, 1:nn*ndof), cdsys_id, coords, time, tincr, u(1:3, 1:nn) )
118  endif
119  else if( fstrsolid%sections(isect)%elemopt361 == kel361ic ) then ! incompatible element
120  if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres >0 ) then
121  CALL stf_c3d8ic &
122  ( ic_type, nn, ecoord(:,1:nn), fstrsolid%elements(icel)%gausses(:), &
123  stiffness(1:nn*ndof, 1:nn*ndof), cdsys_id, coords, time, tincr, u(1:3,1:nn), &
124  fstrsolid%elements(icel)%aux, tt(1:nn) )
125  else
126  CALL stf_c3d8ic &
127  ( ic_type, nn, ecoord(:,1:nn), fstrsolid%elements(icel)%gausses(:), &
128  stiffness(1:nn*ndof, 1:nn*ndof), cdsys_id, coords, time, tincr, u(1:3,1:nn), &
129  fstrsolid%elements(icel)%aux )
130  endif
131  else if( fstrsolid%sections(isect)%elemopt361 == kel361fbar ) then ! F-bar element
132  if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres >0 ) then
133  call stf_c3d8fbar &
134  ( ic_type, nn, ecoord(:, 1:nn), fstrsolid%elements(icel)%gausses(:), &
135  stiffness(1:nn*ndof,1:nn*ndof), cdsys_id, coords, time, tincr, u(1:3, 1:nn), tt(1:nn) )
136  else
137  call stf_c3d8fbar &
138  ( ic_type, nn, ecoord(:, 1:nn),fstrsolid%elements(icel)%gausses(:), &
139  stiffness(1:nn*ndof, 1:nn*ndof), cdsys_id, coords, time, tincr, u(1:3, 1:nn) )
140  endif
141  endif
142 
143  elseif (ic_type==341 .or. ic_type==351 .or. &
144  ic_type==342 .or. ic_type==352 .or. ic_type==362 ) then
145  if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres >0 ) then
146  call stf_c3 &
147  ( ic_type, nn, ecoord(:, 1:nn), fstrsolid%elements(icel)%gausses(:), &
148  stiffness(1:nn*ndof, 1:nn*ndof), cdsys_id, coords, time, tincr, u(1:3,1:nn), tt(1:nn) )
149  else
150  call stf_c3 &
151  ( ic_type,nn,ecoord(:, 1:nn),fstrsolid%elements(icel)%gausses(:), &
152  stiffness(1:nn*ndof, 1:nn*ndof), cdsys_id, coords, time, tincr, u(1:3, 1:nn) )
153  endif
154 
155  else if( ic_type == 611) then
156  if( material%nlgeom_flag /= infinitesimal ) call stiffmat_abort( ic_type, 2 )
157  isect = hecmesh%section_ID(icel)
158  ihead = hecmesh%section%sect_R_index(isect-1)
159  call stf_beam(ic_type, nn, ecoord, hecmesh%section%sect_R_item(ihead+1:), &
160  & material%variables(m_youngs), material%variables(m_poisson), stiffness(1:nn*ndof,1:nn*ndof))
161 
162  else if( ic_type == 641 ) then
163  if( material%nlgeom_flag /= infinitesimal ) call stiffmat_abort( ic_type, 2 )
164  isect = hecmesh%section_ID(icel)
165  ihead = hecmesh%section%sect_R_index(isect-1)
166  call stf_beam_641(ic_type, nn, ecoord, fstrsolid%elements(icel)%gausses(:), &
167  & hecmesh%section%sect_R_item(ihead+1:), stiffness(1:nn*ndof,1:nn*ndof))
168 
169  else if( ( ic_type == 741 ) .or. ( ic_type == 743 ) .or. ( ic_type == 731 ) ) then
170  if( material%nlgeom_flag /= infinitesimal ) call stiffmat_abort( ic_type, 2 )
171  isect = hecmesh%section_ID(icel)
172  ihead = hecmesh%section%sect_R_index(isect-1)
173  thick = hecmesh%section%sect_R_item(ihead+1)
174  call stf_shell_mitc(ic_type, nn, ndof, ecoord(1:3, 1:nn), fstrsolid%elements(icel)%gausses(:), &
175  & stiffness(1:nn*ndof, 1:nn*ndof), thick, 0)
176 
177  else if( ic_type == 761 ) then !for shell-solid mixed analysis
178  if( material%nlgeom_flag /= infinitesimal ) call stiffmat_abort( ic_type, 2 )
179  isect = hecmesh%section_ID(icel)
180  ihead = hecmesh%section%sect_R_index(isect-1)
181  thick = hecmesh%section%sect_R_item(ihead+1)
182  call stf_shell_mitc(731, 3, 6, ecoord(1:3, 1:3), fstrsolid%elements(icel)%gausses(:), &
183  & stiffness(1:nn*ndof, 1:nn*ndof), thick, 2)
184 
185  else if( ic_type == 781 ) then !for shell-solid mixed analysis
186  if( material%nlgeom_flag /= infinitesimal ) call stiffmat_abort( ic_type, 2 )
187  isect = hecmesh%section_ID(icel)
188  ihead = hecmesh%section%sect_R_index(isect-1)
189  thick = hecmesh%section%sect_R_item(ihead+1)
190  call stf_shell_mitc(741, 4, 6, ecoord(1:3, 1:4), fstrsolid%elements(icel)%gausses(:), &
191  & stiffness(1:nn*ndof, 1:nn*ndof), thick, 1)
192 
193  elseif ( ic_type==3414 ) then
194  if(fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype /= incomp_newtonian) then
195  write(*, *) '###ERROR### : This element is not supported for this material'
196  write(*, *) 'ic_type = ', ic_type, ', mtype = ', fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype
197  call hecmw_abort(hecmw_comm_get_comm())
198  endif
199  call stf_c3_vp &
200  ( ic_type, nn, ecoord(:, 1:nn),fstrsolid%elements(icel)%gausses(:), &
201  stiffness(1:nn*ndof, 1:nn*ndof), tincr, u_prev(1:4, 1:nn) )
202  !
203  ! elseif ( ic_type==731) then
204  ! call STF_S3(xx,yy,zz,ee,pp,thick,local_stf)
205  ! call fstr_local_stf_restore_temp(local_stf, nn*ndof, stiffness)
206  ! elseif ( ic_type==741) then
207  ! call STF_S4(xx,yy,zz,ee,pp,thick,local_stf)
208  ! call fstr_local_stf_restore_temp(local_stf, nn*ndof, stiffness)
209  else
210  call stiffmat_abort( ic_type, 1 )
211  endif
212  !
213  ! ----- CONSTRUCT the GLOBAL MATRIX STARTED
214  call hecmw_mat_ass_elem(hecmat, nn, nodlocal, stiffness)
215 
216  enddo ! icel
217  !$omp end do
218  !$omp end parallel
219  enddo ! itype
220 
221  end subroutine fstr_stiffmatrix
222 
223  subroutine stiffmat_abort( ic_type, flag )
224  integer(kind=kint), intent(in) :: ic_type
225  integer(kind=kint), intent(in) :: flag
226 
227  if( flag == 1 ) then
228  write(*,*) '###ERROR### : Element type not supported for static analysis'
229  else if( flag == 2 ) then
230  write(*,*) '###ERROR### : Element type not supported for nonlinear static analysis'
231  endif
232  write(*,*) ' ic_type = ', ic_type
233  call hecmw_abort(hecmw_comm_get_comm())
234  end subroutine
235 
236 end module m_fstr_stiffmatrix
This module provides function to calcualte tangent stiffness matrix.
subroutine, public fstr_stiffmatrix(hecMESH, hecMAT, fstrSOLID, time, tincr)
接線剛性マトリックスを作成するサブルーチン
subroutine stiffmat_abort(ic_type, flag)
This module defined coomon data and basic structures for analysis.
Definition: m_fstr.f90:15
integer(kind=kint), parameter kel361bbar
Definition: m_fstr.f90:69
subroutine get_coordsys(cdsys_ID, hecMESH, fstrSOLID, coords)
This subroutine fetch coords defined by local coordinate system.
Definition: m_fstr.f90:1032
integer(kind=kint), parameter kel361fi
section control
Definition: m_fstr.f90:68
integer(kind=kint), parameter kel361ic
Definition: m_fstr.f90:70
integer(kind=kint), parameter kel361fbar
Definition: m_fstr.f90:71
This modules just summarizes all modules used in static analysis.
Definition: static_LIB.f90:6
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6