FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
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
14contains
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
236end module m_fstr_stiffmatrix
This module provides function to calcualte tangent stiffness matrix.
subroutine, public fstr_stiffmatrix(hecmesh, hecmat, fstrsolid, time, tincr)
接線剛性マトリックスを作成するサブルーチン
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