FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
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
8  use elementinfo
9  implicit none
10 
11 contains
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 
371 end 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 elementstress_c1(ETYPE, gausses, strain, stress)
subroutine nodalstress_c1(ETYPE, NN, gausses, ndstrain, ndstress)
subroutine search_diag_modify(n, nn, hecMAT, hecMESH)
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 truss_diag_modify(hecMAT, hecMESH)
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