FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
static_LIB_1d.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!-------------------------------------------------------------------------------
7 use hecmw, only : kint, kreal
9 implicit none
10
11contains
12 !
13 !=====================================================================*
14 ! STF_C1
15 !=====================================================================*
17 subroutine stf_c1( etype,nn,ecoord,area,gausses,stiff, u ,temperature )
18 use mmechgauss
19 integer(kind=kint), intent(in) :: etype
20 integer(kind=kint), intent(in) :: nn
21 real(kind=kreal), intent(in) :: ecoord(3,nn)
22 real(kind=kreal), intent(in) :: area
23 type(tgaussstatus), intent(in) :: gausses(:)
24 real(kind=kreal), intent(out) :: stiff(:,:)
25 real(kind=kreal), intent(in), optional :: u(:,:)
26 real(kind=kreal), intent(in), optional :: temperature(nn)
27
28 real(kind=kreal) llen, llen0, elem(3,nn)
29 logical :: ierr
30 real(kind=kreal) ina(1), outa(1), direc(3), direc0(3), coeff, strain
31 integer(kind=kint) :: i,j
32
33 ! we suppose the same material type in the element
34 if( present(u) ) then
35 elem(:,:) = ecoord(:,:) + u(:,:)
36 else
37 elem(:,:) = ecoord(:,:)
38 endif
39
40 direc = elem(:,2)-elem(:,1)
41 llen = dsqrt( dot_product(direc, direc) )
42 direc = direc/llen
43 direc0 = ecoord(:,2)-ecoord(:,1)
44 llen0 = dsqrt( dot_product(direc0, direc0) )
45
46 if( present(temperature) ) then
47 ina(1) = 0.5d0*(temperature(1)+temperature(2))
48 call fetch_tabledata( mc_isoelastic, gausses(1)%pMaterial%dict, outa, ierr, ina )
49 else
50 call fetch_tabledata( mc_isoelastic, gausses(1)%pMaterial%dict, outa, ierr )
51 endif
52 if( ierr ) outa(1) = gausses(1)%pMaterial%variables(m_youngs)
53 coeff = outa(1)*area*llen0/(llen*llen)
54 strain = gausses(1)%strain(1)
55
56 stiff(:,:) = 0.d0
57 do i=1,3
58 stiff(i,i) = coeff*strain
59 do j=1,3
60 stiff(i,j) = stiff(i,j) + coeff*(1.d0-2.d0*strain)*direc(i)*direc(j)
61 enddo
62 enddo
63
64 stiff(4:6,1:3) = -stiff(1:3,1:3)
65 stiff(1:3,4:6) = transpose(stiff(4:6,1:3))
66 stiff(4:6,4:6) = stiff(1:3,1:3)
67
68 end subroutine stf_c1
69
70 !
72 !---------------------------------------------------------------------*
73 subroutine update_c1( etype, nn, ecoord, area, u, du, qf ,gausses, TT, T0 )
74 !---------------------------------------------------------------------*
75 use m_fstr
76 use mmechgauss
77 ! I/F VARIAVLES
78 integer(kind=kint), intent(in) :: etype
79 integer(kind=kint), intent(in) :: nn
80 real(kind=kreal), intent(in) :: ecoord(3,nn)
81 real(kind=kreal), intent(in) :: area
82 real(kind=kreal), intent(in) :: u(3,nn)
83 real(kind=kreal), intent(in) :: du(3,nn)
84 real(kind=kreal), intent(out) :: qf(nn*3)
85 type(tgaussstatus), intent(inout) :: gausses(:)
86 real(kind=kreal), intent(in), optional :: tt(nn)
87 real(kind=kreal), intent(in), optional :: t0(nn)
88
89 ! LCOAL VARIAVLES
90 real(kind=kreal) :: direc(3), direc0(3)
91 real(kind=kreal) :: llen, llen0, ina(1), outa(1)
92 real(kind=kreal) :: elem(3,nn)
93 real(kind=kreal) :: young
94 real(kind=kreal) :: ttc, tt0, alp, alp0, epsth
95 logical :: ierr
96
97 qf(:) = 0.d0
98 ! we suppose the same material type in the element
99 elem(:,:) = ecoord(:,:) + u(:,:) + du(:,:)
100
101 direc = elem(:,2)-elem(:,1)
102 llen = dsqrt( dot_product(direc, direc) )
103 direc = direc/llen
104 direc0 = ecoord(:,2)-ecoord(:,1)
105 llen0 = dsqrt( dot_product(direc0, direc0) )
106
107 epsth = 0.d0
108 if( present(tt) .and. present(t0) ) then
109 ttc = 0.5d0*(tt(1)+tt(2))
110 tt0 = 0.5d0*(t0(1)+t0(2))
111
112 ina(1) = ttc
113 call fetch_tabledata( mc_isoelastic, gausses(1)%pMaterial%dict, outa, ierr, ina )
114 if( ierr ) outa(1) = gausses(1)%pMaterial%variables(m_youngs)
115 young = outa(1)
116
117 call fetch_tabledata( mc_themoexp, gausses(1)%pMaterial%dict, outa(:), ierr, ina )
118 if( ierr ) outa(1) = gausses(1)%pMaterial%variables(m_exapnsion)
119 alp = outa(1)
120
121 ina(1) = tt0
122 call fetch_tabledata( mc_themoexp, gausses(1)%pMaterial%dict, outa(:), ierr, ina )
123 if( ierr ) outa(1) = gausses(1)%pMaterial%variables(m_exapnsion)
124 alp0 = outa(1)
125
126 epsth=alp*(ttc-ref_temp)-alp0*(tt0-ref_temp)
127 else
128 call fetch_tabledata( mc_isoelastic, gausses(1)%pMaterial%dict, outa, ierr )
129 if( ierr ) outa(1) = gausses(1)%pMaterial%variables(m_youngs)
130 young = outa(1)
131 endif
132
133 gausses(1)%strain(1) = dlog(llen/llen0)
134 gausses(1)%stress(1) = young*(gausses(1)%strain(1)-epsth)
135
136 !set stress and strain for output
137 gausses(1)%strain_out(1) = gausses(1)%strain(1)
138 gausses(1)%stress_out(1) = gausses(1)%stress(1)
139
140 qf(1) = gausses(1)%stress(1)*area*llen0/llen
141 qf(1:3) = -qf(1)*direc
142 qf(4:6) = -qf(1:3)
143
144 end subroutine update_c1
145
146 !----------------------------------------------------------------------*
147 subroutine nodalstress_c1(ETYPE,NN,gausses,ndstrain,ndstress)
148 !----------------------------------------------------------------------*
149 !
150 ! Calculate Strain and Stress increment of solid elements
151 !
152 use mmechgauss
153 integer(kind=kint), intent(in) :: etype,nn
154 type(tgaussstatus), intent(in) :: gausses(:)
155 real(kind=kreal), intent(out) :: ndstrain(nn,6)
156 real(kind=kreal), intent(out) :: ndstress(nn,6)
157
158 ndstrain(1,1:6) = gausses(1)%strain(1:6)
159 ndstress(1,1:6) = gausses(1)%stress(1:6)
160 ndstrain(2,1:6) = gausses(1)%strain(1:6)
161 ndstress(2,1:6) = gausses(1)%stress(1:6)
162
163 end subroutine
164 !
165 !
166 !
167 !----------------------------------------------------------------------*
168 subroutine elementstress_c1(ETYPE,gausses,strain,stress)
169 !----------------------------------------------------------------------*
170 !
171 ! Calculate Strain and Stress increment of solid elements
172 !
173 use mmechgauss
174 integer(kind=kint), intent(in) :: ETYPE
175 type(tgaussstatus), intent(in) :: gausses(:)
176 real(kind=kreal), intent(out) :: strain(6)
177 real(kind=kreal), intent(out) :: stress(6)
178
179
180 strain(:) = gausses(1)%strain(1:6)
181 stress(:) = gausses(1)%stress(1:6)
182
183 end subroutine
184
185
186 !----------------------------------------------------------------------*
187 subroutine dl_c1(etype, nn, xx, yy, zz, rho, thick, ltype, params, &
188 vect, nsize)
189 !----------------------------------------------------------------------*
190 !** SET DLOAD
191 ! GRAV LTYPE=4 :GRAVITY FORCE
192 ! I/F VARIABLES
193 integer(kind = kint), intent(in) :: etype, nn
194 real(kind = kreal), intent(in) :: xx(:), yy(:), zz(:)
195 real(kind = kreal), intent(in) :: params(0:6)
196 real(kind = kreal), intent(inout) :: vect(:)
197 real(kind = kreal) :: rho, thick
198 integer(kind = kint) :: ltype, nsize
199 ! LOCAL VARIABLES
200 integer(kind = kint) :: ndof = 3
201 integer(kind = kint) :: ivol, isuf, i
202 real(kind = kreal) :: vx, vy, vz, val, a, aa
203 !--------------------------------------------------------------------
204 val = params(0)
205 !--------------------------------------------------------------
206
207 ivol = 0
208 isuf = 0
209 if( ltype .LT. 10 ) then
210 ivol = 1
211 else if( ltype .GE. 10 ) then
212 isuf = 1
213 end if
214
215 !--------------------------------------------------------------------
216 nsize = nn*ndof
217 !--------------------------------------------------------------------
218 vect(1:nsize) = 0.0d0
219
220 ! Volume force
221 if( ivol .EQ. 1 ) then
222 if( ltype .EQ. 4 ) then
223 aa = dsqrt( ( xx(2)-xx(1) )*( xx(2)-xx(1) ) &
224 +( yy(2)-yy(1) )*( yy(2)-yy(1) ) &
225 +( zz(2)-zz(1) )*( zz(2)-zz(1) ) )
226
227 a = thick
228 vx = params(1)
229 vy = params(2)
230 vz = params(3)
231 vx = vx/dsqrt( params(1)**2+params(2)**2+params(3)**2 )
232 vy = vy/dsqrt( params(1)**2+params(2)**2+params(3)**2 )
233 vz = vz/dsqrt( params(1)**2+params(2)**2+params(3)**2 )
234
235 do i = 1, 2
236 vect(3*i-2) = val*rho*a*0.5d0*aa*vx
237 vect(3*i-1) = val*rho*a*0.5d0*aa*vy
238 vect(3*i ) = val*rho*a*0.5d0*aa*vz
239 end do
240
241 end if
242 end if
243 !--------------------------------------------------------------------
244
245 return
246 end subroutine
247
248 !
249 !
250 !----------------------------------------------------------------------*
251 subroutine truss_diag_modify(hecMAT,hecMESH)
252 !----------------------------------------------------------------------*
253 !
254 use hecmw
255 type (hecmwST_matrix) :: hecMAT
256 type (hecmwST_local_mesh) :: hecMESH
257 integer(kind=kint) :: itype, is, iE, ic_type, icel, jS, j, n
258
259 do itype = 1, hecmesh%n_elem_type
260 ic_type = hecmesh%elem_type_item(itype)
261 if(ic_type == 301)then
262 is = hecmesh%elem_type_index(itype-1) + 1
263 ie = hecmesh%elem_type_index(itype )
264 do icel = is, ie
265 js = hecmesh%elem_node_index(icel-1)
266 do j=1,2
267 n = hecmesh%elem_node_item(js+j)
268 if( hecmat%D(9*n-8) == 0.0d0)then
269 hecmat%D(9*n-8) = 1.0d0
270 !call search_diag_modify(n,1,hecMAT,hecMESH)
271 endif
272 if( hecmat%D(9*n-4) == 0.0d0)then
273 hecmat%D(9*n-4) = 1.0d0
274 !call search_diag_modify(n,2,hecMAT,hecMESH)
275 endif
276 if( hecmat%D(9*n ) == 0.0d0)then
277 hecmat%D(9*n ) = 1.0d0
278 !call search_diag_modify(n,3,hecMAT,hecMESH)
279 endif
280 enddo
281 enddo
282 endif
283 enddo
284
285 end subroutine
286
287 !
288 !
289 !----------------------------------------------------------------------*
290 subroutine search_diag_modify(n,nn,hecMAT,hecMESH)
291 !----------------------------------------------------------------------*
292 !
293 use hecmw
294 type (hecmwST_matrix) :: hecMAT
295 type (hecmwST_local_mesh) :: hecMESH
296 integer :: n, nn, is, iE, i, a
297
298 if(nn == 1)then
299 a = 0
300 is = hecmat%IndexL(n-1)+1
301 ie = hecmat%IndexL(n )
302 do i=is,ie
303 if(hecmat%AL(9*i-8) /= 0.0d0) a = 1
304 if(hecmat%AL(9*i-7) /= 0.0d0) a = 1
305 if(hecmat%AL(9*i-6) /= 0.0d0) a = 1
306 enddo
307 is = hecmat%IndexU(n-1)+1
308 ie = hecmat%IndexU(n )
309 do i=is,ie
310 if(hecmat%AU(9*i-8) /= 0.0d0) a = 1
311 if(hecmat%AU(9*i-7) /= 0.0d0) a = 1
312 if(hecmat%AU(9*i-6) /= 0.0d0) a = 1
313 enddo
314 if(hecmat%D(9*n-7) /= 0.0d0) a = 1
315 if(hecmat%D(9*n-6) /= 0.0d0) a = 1
316 if(a == 0)then
317 hecmat%D(9*n-8) = 1.0d0
318 !write(*,"(a,i,a,i,a)")"### FIX DIAGONAL n:",n,", ID:",hecMESH%global_node_ID(n),", dof:1"
319 endif
320 endif
321 if(nn == 2)then
322 a = 0
323 is = hecmat%IndexL(n-1)+1
324 ie = hecmat%IndexL(n )
325 do i=is,ie
326 if(hecmat%AL(9*i-5) /= 0.0d0) a = 1
327 if(hecmat%AL(9*i-4) /= 0.0d0) a = 1
328 if(hecmat%AL(9*i-3) /= 0.0d0) a = 1
329 enddo
330 is = hecmat%IndexU(n-1)+1
331 ie = hecmat%IndexU(n )
332 do i=is,ie
333 if(hecmat%AU(9*i-5) /= 0.0d0) a = 1
334 if(hecmat%AU(9*i-4) /= 0.0d0) a = 1
335 if(hecmat%AU(9*i-3) /= 0.0d0) a = 1
336 enddo
337 if(hecmat%D(9*n-5) /= 0.0d0) a = 1
338 if(hecmat%D(9*n-3) /= 0.0d0) a = 1
339 if(a == 0)then
340 hecmat%D(9*n-4) = 1.0d0
341 !write(*,"(a,i,a,i,a)")"### FIX DIAGONAL n:",n,", ID:",hecMESH%global_node_ID(n),", dof:2"
342 endif
343 endif
344 if(nn == 3)then
345 a = 0
346 is = hecmat%IndexL(n-1)+1
347 ie = hecmat%IndexL(n )
348 do i=is,ie
349 if(hecmat%AL(9*i-2) /= 0.0d0) a = 1
350 if(hecmat%AL(9*i-1) /= 0.0d0) a = 1
351 if(hecmat%AL(9*i ) /= 0.0d0) a = 1
352 enddo
353 is = hecmat%IndexU(n-1)+1
354 ie = hecmat%IndexU(n )
355 do i=is,ie
356 if(hecmat%AU(9*i-2) /= 0.0d0) a = 1
357 if(hecmat%AU(9*i-1) /= 0.0d0) a = 1
358 if(hecmat%AU(9*i ) /= 0.0d0) a = 1
359 enddo
360 if(hecmat%D(9*n-2) /= 0.0d0) a = 1
361 if(hecmat%D(9*n-1) /= 0.0d0) a = 1
362 if(a == 0)then
363 hecmat%D(9*n ) = 1.0d0
364 !write(*,"(a,i,a,i,a)")"### FIX DIAGONAL n:",n,", ID:",hecMESH%global_node_ID(n),", dof:3"
365 endif
366 endif
367
368 end subroutine
369
370
371end module m_static_lib_1d
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
Definition: hecmw.f90:6
This module defined coomon data and basic structures for analysis.
Definition: m_fstr.f90:15
real(kind=kreal), pointer ref_temp
REFTEMP.
Definition: m_fstr.f90:120
This module provide common functions of 3D truss elements.
subroutine truss_diag_modify(hecmat, hecmesh)
subroutine search_diag_modify(n, nn, hecmat, hecmesh)
subroutine nodalstress_c1(etype, nn, gausses, ndstrain, ndstress)
subroutine update_c1(etype, nn, ecoord, area, u, du, qf, gausses, tt, t0)
Update strain and stress inside element.
subroutine stf_c1(etype, nn, ecoord, area, gausses, stiff, u, temperature)
This subroutine calculate stiff matrix of 2-nodes truss element.
subroutine elementstress_c1(etype, gausses, strain, stress)
subroutine dl_c1(etype, nn, xx, yy, zz, rho, thick, ltype, params, vect, nsize)
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13