FrontISTR  5.7.0
Large-scale structural analysis program with finit element method
ElasticLinear.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 module m_elasticlinear
7  use hecmw_util
8  use mmaterial
9 
10  implicit none
11 
12 contains
13 
15  subroutine calelasticmatrix( matl, sectType, D, temp, hdflag )
16  type( tmaterial ), intent(in) :: matl
17  integer, intent(in) :: sectType
18  real(kind=kreal), intent(out) :: d(:,:)
19  real(kind=kreal), intent(in) :: temp
20  real(kind=kreal) :: ee, pp, coef1, coef2, ina(1), outa(2)
21  integer(kind=kint), intent(in), optional :: hdflag
22  logical :: ierr
23  real(kind=kreal) :: k, g
24 
25  d(:,:)=0.d0
26 
27  ina(1) = temp
28  call fetch_tabledata( mc_isoelastic, matl%dict, outa, ierr, ina )
29  if( ierr ) then
30  ee = matl%variables(m_youngs)
31  pp = matl%variables(m_poisson)
32  else
33  ee = outa(1)
34  pp = outa(2)
35  endif
36 
37  select case (secttype)
38  case (d3)
39  k = ee/(1.d0-2.d0*pp)/3.d0
40  g = ee/(1.d0+pp)*0.5d0
41  if( present(hdflag) ) then
42  if( hdflag == 1 ) k = 0.d0
43  if( hdflag == 2 ) g = 0.d0
44  end if
45  d(1,1)=k+(4.d0/3.d0)*g
46  d(1,2)=k-(2.d0/3.d0)*g
47  d(4,4)=g
48  d(1,3)=d(1,2)
49  d(2,1)=d(1,2)
50  d(2,2)=d(1,1)
51  d(2,3)=d(1,2)
52  d(3,1)=d(1,3)
53  d(3,2)=d(2,3)
54  d(3,3)=d(1,1)
55  d(5,5)=d(4,4)
56  d(6,6)=d(4,4)
57  case (planestress)
58  coef1=ee/(1.d0-pp*pp)
59  coef2=0.5d0*(1.d0-pp)
60  d(1,1)=coef1
61  d(1,2)=coef1*pp
62  d(1,3)=0.d0
63  d(2,1)=d(1,2)
64  d(2,2)=d(1,1)
65  d(2,3)=0.d0
66  d(3,1)=0.d0
67  d(3,2)=0.d0
68  d(3,3)=coef1*coef2
69  case (planestrain)
70  coef1=ee/((1.d0+pp)*(1.d0-2.d0*pp))
71  coef2=ee/(2.d0*(1.d0+pp))
72  d(1,1)=coef1*(1.d0-pp)
73  d(1,2)=coef1*pp
74  d(1,3)=0.d0
75  d(2,1)=d(1,2)
76  d(2,2)=d(1,1)
77  d(2,3)=0.d0
78  d(3,1)=0.d0
79  d(3,2)=0.d0
80  d(3,3)=coef2
81  case (axissymetric)
82  coef1=ee*(1.d0-pp)/((1.d0+pp)*(1.d0-2.d0*pp))
83  coef2=(1.d0-2.d0*pp)/(2.d0*(1.d0-pp))
84  d(1,1)=coef1
85  d(1,2)=coef1*pp/(1.d0-pp)
86  d(1,3)=0.d0
87  d(1,4)=d(1,2)
88  d(2,1)=d(1,2)
89  d(2,2)=d(1,1)
90  d(2,3)=0.d0
91  d(2,4)=d(1,2)
92  d(3,1)=0.d0
93  d(3,2)=0.d0
94  d(3,3)=coef1*coef2
95  d(3,4)=0.d0
96  d(4,1)=d(1,4)
97  d(4,2)=d(2,4)
98  d(4,3)=0.d0
99  d(4,4)=d(1,1)
100  case default
101  stop "Section type not defined"
102  end select
103 
104  end subroutine
105 
106 
108  subroutine calelasticmatrix_ortho( matl, sectType, bij, DMAT, temp )
109  use m_utilities
110  type( tmaterial ), intent(in) :: matl
111  integer, intent(in) :: secttype
112  real(kind=kreal), intent(in) :: bij(3,3)
113  real(kind=kreal), intent(out) :: dmat(:,:)
114  real(kind=kreal), intent(in) :: temp
115  real(kind=kreal) :: e1, e2, e3, g12, g23, g13, nyu12, nyu23,nyu13
116  real(kind=kreal) :: nyu21,nyu32,nyu31, delta1, ina(1), outa(9)
117  real(kind=kreal) :: tm(6,6)
118  logical :: ierr
119 
120  ina(1) = temp
121  call fetch_tabledata( mc_orthoelastic, matl%dict, outa, ierr, ina )
122  if( ierr ) then
123  stop "Fails in fetching orthotropic elastic constants!"
124  endif
125 
126  e1 = outa(1)
127  e2 = outa(2)
128  e3 = outa(3)
129  nyu12 = outa(4)
130  nyu13 = outa(5)
131  nyu23 = outa(6)
132  g12 = outa(7)
133  g13 = outa(8)
134  g23 = outa(9)
135  nyu21 = e2/e1*nyu12
136  nyu32 = e3/e2*nyu23
137  nyu31 = e3/e1*nyu13
138  delta1 = 1.d0/(1.d0 -nyu12*nyu21 -nyu23*nyu32 -nyu31*nyu13 -2.d0*nyu21*nyu32*nyu13)
139 
140  dmat(:,:)=0.d0
141  dmat(1,1) = e1*(1.d0-nyu23*nyu32)*delta1
142  dmat(2,2) = e2*(1.d0-nyu13*nyu31)*delta1
143  dmat(3,3) = e3*(1.d0-nyu12*nyu21)*delta1
144  dmat(1,2) = e1*(nyu21+nyu31*nyu23)*delta1
145  dmat(1,3) = e1*(nyu31+nyu21*nyu32)*delta1
146  dmat(2,3) = e2*(nyu32+nyu12*nyu31)*delta1
147  dmat(4,4) = g12
148  dmat(5,5) = g23
149  dmat(6,6) = g13
150 
151  dmat(2,1) = dmat(1,2)
152  dmat(3,2) = dmat(2,3)
153  dmat(3,1) = dmat(1,3)
154 
155  call transformation(bij, tm)
156 
157  dmat = matmul( transpose(tm), dmat)
158  dmat = matmul( dmat, (tm) )
159 
160  end subroutine
161 
162 
163  !####################################################################
164  subroutine linearelastic_shell &
165  (matl, secttype, c, &
166  e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, &
167  alpha, n_layer)
168  !####################################################################
169 
170  type( tmaterial ), intent(in) :: matl
171  integer, intent(in) :: sectType, n_layer
172  real(kind = kreal), intent(out) :: c(:, :, :, :)
173  real(kind = kreal), intent(in) :: e1_hat(3), e2_hat(3), e3_hat(3)
174  real(kind = kreal), intent(in) :: cg1(3), cg2(3), cg3(3)
175  real(kind = kreal), intent(out) :: alpha
176 
177  !--------------------------------------------------------------------
178 
179  real(kind = kreal) :: ee, pp, ee2, g12, g23, g31, theta, pp2
180  real(kind = kreal) :: outa(2)
181  real(kind = kreal) :: lambda1, lambda2, mu, k_correction
182  real(kind = kreal) :: c_hat(3, 3, 3, 3), d(5,5), d_hat(5,5), d_temp(5,5), t(5,5)
183  real(kind = kreal) :: e_hat_dot_cg(3, 3)
184  real(kind = kreal) :: alpha_over_mu
185 
186  integer :: index_i(5), index_j(5), index_k(5), index_l(5)
187  integer :: is, js, ii, ij, ik, il
188  integer :: i, j, k, l, total_layer, matl_type
189 
190  logical :: ierr
191 
192  !--------------------------------------------------------------------
193 
194  call fetch_tabledata(mc_isoelastic, matl%dict, outa, ierr)
195 
196  !--------------------------------------------------------------------
197 
198  matl_type = matl%shell_var(n_layer)%ortho
199  if(matl_type == 0)then
200  if( ierr ) then
201 
202  ee = matl%shell_var(n_layer)%ee
203  pp = matl%shell_var(n_layer)%pp
204  alpha_over_mu = matl%variables(m_alpha_over_mu)
205 
206  else
207 
208  ee = outa(1)
209  pp = outa(2)
210  alpha_over_mu = matl%variables(m_alpha_over_mu)
211 
212  end if
213 
214  !--------------------------------------------------------------------
215 
216  ! Elastic constant
217  lambda1 = ee/( 1.0d0-pp*pp )
218  lambda2 = pp*lambda1
219  mu = 0.5d0*ee/( 1.0d0+pp )
220 
221  !--------------------------------------------------------------------
222 
223  ! Shear correction factor
224  k_correction = 5.0d0/6.0d0
225  !k_correction = 1.0D0
226 
227  !--------------------------------------------------------------------
228 
229  alpha = alpha_over_mu*mu
230 
231  !--------------------------------------------------------------------
232 
233  ! Constitutive tensor
234  c_hat(:, :, :, :) = 0.0d0
235 
236  !--------------------------------------------------------
237 
238  c_hat(1, 1, 1, 1) = lambda1
239  c_hat(1, 1, 2, 2) = lambda2
240  c_hat(2, 2, 1, 1) = lambda2
241  c_hat(2, 2, 2, 2) = lambda1
242  c_hat(1, 2, 1, 2) = mu
243  c_hat(1, 2, 2, 1) = mu
244  c_hat(2, 1, 1, 2) = mu
245  c_hat(2, 1, 2, 1) = mu
246  c_hat(1, 3, 1, 3) = k_correction*mu
247  c_hat(1, 3, 3, 1) = k_correction*mu
248  c_hat(2, 3, 2, 3) = k_correction*mu
249  c_hat(2, 3, 3, 2) = k_correction*mu
250  c_hat(3, 1, 3, 1) = k_correction*mu
251  c_hat(3, 1, 1, 3) = k_correction*mu
252  c_hat(3, 2, 3, 2) = k_correction*mu
253  c_hat(3, 2, 2, 3) = k_correction*mu
254 
255  !--------------------------------------------------------
256 
257  elseif(matl_type == 1)then
258  total_layer = matl%totallyr
259 
260  ee = matl%shell_var(n_layer)%ee
261  pp = matl%shell_var(n_layer)%pp
262  ee2 = matl%shell_var(n_layer)%ee2
263  g12 = matl%shell_var(n_layer)%g12
264  g23 = matl%shell_var(n_layer)%g23
265  g31 = matl%shell_var(n_layer)%g31
266  theta = matl%shell_var(n_layer)%angle
267 
268  alpha_over_mu = matl%variables(m_alpha_over_mu)
269 
270  !--------------------------------------------------------------------
271 
272  ! Shear correction factor
273  k_correction = 5.0d0/6.0d0
274 
275  !--------------------------------------------------------------------
276 
277  alpha = alpha_over_mu * 0.5d0 * ee / ( 1.0d0+pp )
278 
279  !--------------------------------------------------------------------
280 
281  ! write(*,*) ee,pp,ee2,g12,g23,g31,theta,alpha_over_mu,n_layer
282 
283  d(:,:) = 0.0d0
284  d_hat(:,:) = 0.0d0
285  d_temp(:,:) = 0.0d0
286  t(:,:) = 0.0d0
287 
288  pp2 = pp * ee2 / ee
289 
290  d(1,1) = ee / (1.0d0 - pp * pp2)
291  d(1,2) = pp2 * ee / (1.0d0- pp * pp2)
292  d(2,1) = pp2 * ee / (1.0d0- pp * pp2)
293  d(2,2) = ee2 / (1.0d0 - pp * pp2)
294  d(3,3) = g12
295  d(4,4) = g23
296  d(5,5) = g31
297 
298  t(1,1) = cos(theta) * cos(theta)
299  t(1,2) = sin(theta) * sin(theta)
300  t(2,1) = sin(theta) * sin(theta)
301  t(2,2) = cos(theta) * cos(theta)
302  t(3,3) = cos(theta) * cos(theta) - sin(theta) * sin(theta)
303  t(1,3) = sin(theta) * cos(theta)
304  t(2,3) = -sin(theta) * cos(theta)
305  t(3,1) = -2.0d0 * sin(theta) * cos(theta)
306  t(3,2) = 2.0d0 * sin(theta) * cos(theta)
307  t(4,4) = cos(theta)
308  t(4,5) = sin(theta)
309  t(5,4) = -sin(theta)
310  t(5,5) = cos(theta)
311 
312  !--------------------- D_temp = [D]*[T]
313 
314  d_temp(1,1) = d(1,1)*t(1,1)+d(1,2)*t(2,1)
315  d_temp(1,2) = d(1,1)*t(1,2)+d(1,2)*t(2,2)
316  d_temp(2,1) = d(2,1)*t(1,1)+d(2,2)*t(2,1)
317  d_temp(2,2) = d(2,1)*t(1,2)+d(2,2)*t(2,2)
318  d_temp(3,1) = d(3,3)*t(3,1)
319  d_temp(3,2) = d(3,3)*t(3,2)
320  d_temp(1,3) = d(1,1)*t(1,3)+d(1,2)*t(2,3)
321  d_temp(2,3) = d(2,1)*t(1,3)+d(2,2)*t(2,3)
322  d_temp(3,3) = d(3,3)*t(3,3)
323  d_temp(4,4) = d(4,4)*t(4,4)
324  d_temp(4,5) = d(4,4)*t(4,5)
325  d_temp(5,4) = d(5,5)*t(5,4)
326  d_temp(5,5) = d(5,5)*t(5,5)
327 
328  !--------------------- D_hat = [trans_T]*[D_temp]
329 
330  d_hat(1,1) = t(1,1)*d_temp(1,1)+t(1,2)*d_temp(2,1)+t(3,1)*d_temp(3,1)
331  d_hat(1,2) = t(1,1)*d_temp(1,2)+t(1,2)*d_temp(2,2)+t(3,1)*d_temp(3,2)
332  d_hat(2,1) = t(2,1)*d_temp(1,1)+t(2,2)*d_temp(2,1)+t(3,2)*d_temp(3,1)
333  d_hat(2,2) = t(2,1)*d_temp(1,2)+t(2,2)*d_temp(2,2)+t(3,2)*d_temp(3,2)
334  d_hat(3,1) = t(1,3)*d_temp(1,1)+t(2,3)*d_temp(2,1)+t(3,3)*d_temp(3,1)
335  d_hat(3,2) = t(1,3)*d_temp(1,2)+t(2,3)*d_temp(2,2)+t(3,3)*d_temp(3,2)
336  d_hat(1,3) = t(1,1)*d_temp(1,3)+t(1,2)*d_temp(2,3)+t(3,1)*d_temp(3,3)
337  d_hat(2,3) = t(2,1)*d_temp(1,3)+t(2,2)*d_temp(2,3)+t(3,2)*d_temp(3,3)
338  d_hat(3,3) = t(1,3)*d_temp(1,3)+t(2,3)*d_temp(2,3)+t(3,3)*d_temp(3,3)
339  d_hat(4,4) = t(4,4)*d_temp(4,4)+t(5,4)*d_temp(5,4)
340  d_hat(4,5) = t(4,4)*d_temp(4,5)+t(5,4)*d_temp(5,5)
341  d_hat(5,4) = t(4,5)*d_temp(4,4)+t(5,5)*d_temp(5,4)
342  d_hat(5,5) = t(4,5)*d_temp(4,5)+t(5,5)*d_temp(5,5)
343 
344  !--------------------------------------------------------------------
345 
346  ! Constitutive tensor
347  c_hat(:, :, :, :) = 0.0d0
348  !write(*,*) 'Elastic linear.f90 make c_hat'
349 
350  !-------D_hat to c_hat
351 
352  index_i(1) = 1
353  index_i(2) = 2
354  index_i(3) = 1
355  index_i(4) = 2
356  index_i(5) = 3
357 
358  index_j(1) = 1
359  index_j(2) = 2
360  index_j(3) = 2
361  index_j(4) = 3
362  index_j(5) = 1
363 
364  index_k(1) = 1
365  index_k(2) = 2
366  index_k(3) = 1
367  index_k(4) = 2
368  index_k(5) = 3
369 
370  index_l(1) = 1
371  index_l(2) = 2
372  index_l(3) = 2
373  index_l(4) = 3
374  index_l(5) = 1
375 
376  !--------------------------------------------------------------------
377 
378  do js = 1, 5
379  do is = 1, 5
380 
381  ii = index_i(is)
382  ij = index_j(is)
383  ik = index_k(js)
384  il = index_l(js)
385 
386  c_hat(ii, ij, ik, il) = d_hat(is, js)
387 
388  end do
389  end do
390 
391  !--------------------------------------------------------
392 
393  else
394  write(*,*) 'shell matl type isnot collect'
395  stop
396  endif
397 
398  select case( secttype )
399  case( shell )
400 
401  e_hat_dot_cg(1, 1) &
402  = e1_hat(1)*cg1(1)+e1_hat(2)*cg1(2) &
403  +e1_hat(3)*cg1(3)
404  e_hat_dot_cg(2, 1) &
405  = e2_hat(1)*cg1(1)+e2_hat(2)*cg1(2) &
406  +e2_hat(3)*cg1(3)
407  e_hat_dot_cg(3, 1) = 0.0d0
408  e_hat_dot_cg(1, 2) &
409  = e1_hat(1)*cg2(1)+e1_hat(2)*cg2(2) &
410  +e1_hat(3)*cg2(3)
411  e_hat_dot_cg(2, 2) &
412  = e2_hat(1)*cg2(1)+e2_hat(2)*cg2(2) &
413  +e2_hat(3)*cg2(3)
414  e_hat_dot_cg(3, 2) = 0.0d0
415  e_hat_dot_cg(1, 3) &
416  = e1_hat(1)*cg3(1)+e1_hat(2)*cg3(2) &
417  +e1_hat(3)*cg3(3)
418  e_hat_dot_cg(2, 3) &
419  = e2_hat(1)*cg3(1)+e2_hat(2)*cg3(2) &
420  +e2_hat(3)*cg3(3)
421  e_hat_dot_cg(3, 3) &
422  = e3_hat(1)*cg3(1)+e3_hat(2)*cg3(2) &
423  +e3_hat(3)*cg3(3)
424 
425  !--------------------------------------------------------
426 
427  ! Constitutive tensor
428 
429  c(1, 1, 1, 1) = 0.0d0
430  c(2, 2, 1, 1) = 0.0d0
431  c(1, 2, 1, 1) = 0.0d0
432  c(2, 2, 2, 2) = 0.0d0
433  c(1, 2, 2, 2) = 0.0d0
434  c(1, 2, 1, 2) = 0.0d0
435  c(3, 1, 1, 1) = 0.0d0
436  c(3, 1, 2, 2) = 0.0d0
437  c(3, 1, 1, 2) = 0.0d0
438  c(2, 3, 1, 1) = 0.0d0
439  c(2, 3, 2, 2) = 0.0d0
440  c(2, 3, 1, 2) = 0.0d0
441  c(3, 1, 3, 1) = 0.0d0
442  c(3, 1, 2, 3) = 0.0d0
443  c(2, 3, 2, 3) = 0.0d0
444 
445  do l = 1, 2
446 
447  do k = 1, 2
448 
449  do j = 1, 2
450 
451  do i = 1, 2
452 
453  c(1, 1, 1, 1) &
454  = c(1, 1, 1, 1) &
455  +c_hat(i, j, k, l) &
456  *e_hat_dot_cg(i, 1)*e_hat_dot_cg(j ,1) &
457  *e_hat_dot_cg(k, 1)*e_hat_dot_cg(l, 1)
458  c(2, 2, 1, 1) &
459  = c(2, 2, 1, 1) &
460  +c_hat(i, j, k, l) &
461  *e_hat_dot_cg(i, 2)*e_hat_dot_cg(j, 2) &
462  *e_hat_dot_cg(k, 1)*e_hat_dot_cg(l, 1)
463  c(1, 2, 1, 1) &
464  = c(1, 2, 1, 1) &
465  +c_hat(i, j, k, l) &
466  *e_hat_dot_cg(i, 1)*e_hat_dot_cg(j, 2) &
467  *e_hat_dot_cg(k, 1)*e_hat_dot_cg(l, 1)
468  c(2, 2, 2, 2) &
469  = c(2, 2, 2, 2) &
470  +c_hat(i, j, k, l) &
471  *e_hat_dot_cg(i, 2)*e_hat_dot_cg(j, 2) &
472  *e_hat_dot_cg(k, 2)*e_hat_dot_cg(l, 2)
473 
474  c(1, 2, 2, 2) &
475  = c(1, 2, 2, 2) &
476  +c_hat(i, j, k, l) &
477  *e_hat_dot_cg(i, 1)*e_hat_dot_cg(j, 2) &
478  *e_hat_dot_cg(k, 2)*e_hat_dot_cg(l, 2)
479  c(1, 2, 1, 2) &
480  = c(1, 2, 1, 2) &
481  +c_hat(i, j, k, l) &
482  *e_hat_dot_cg(i, 1)*e_hat_dot_cg(j, 2) &
483  *e_hat_dot_cg(k, 1)*e_hat_dot_cg(l, 2)
484 
485  end do
486 
487  do i = 1, 3
488 
489  c(3, 1, 1, 1) &
490  = c(3, 1, 1, 1) &
491  +c_hat(i, j, k, l) &
492  *e_hat_dot_cg(i, 3)*e_hat_dot_cg(j, 1) &
493  *e_hat_dot_cg(k, 1)*e_hat_dot_cg(l, 1)
494  c(3, 1, 2, 2) &
495  = c(3, 1, 2, 2) &
496  +c_hat(i, j, k, l) &
497  *e_hat_dot_cg(i, 3)*e_hat_dot_cg(j, 1) &
498  *e_hat_dot_cg(k, 2)*e_hat_dot_cg(l, 2)
499  c(3, 1, 1, 2) &
500  = c(3, 1, 1, 2) &
501  +c_hat(i, j, k, l) &
502  *e_hat_dot_cg(i, 3)*e_hat_dot_cg(j, 1) &
503  *e_hat_dot_cg(k, 1)*e_hat_dot_cg(l, 2)
504 
505  end do
506 
507  end do
508 
509  do j = 1, 3
510 
511  do i = 1, 2
512 
513  c(2, 3, 1, 1) &
514  = c(2, 3, 1, 1) &
515  +c_hat(i, j, k, l) &
516  *e_hat_dot_cg(i, 2)*e_hat_dot_cg(j, 3) &
517  *e_hat_dot_cg(k, 1)*e_hat_dot_cg(l, 1)
518  c(2, 3, 2, 2) &
519  = c(2, 3, 2, 2) &
520  +c_hat(i, j, k, l) &
521  *e_hat_dot_cg(i, 2)*e_hat_dot_cg(j, 3) &
522  *e_hat_dot_cg(k, 2)*e_hat_dot_cg(l, 2)
523  c(2, 3, 1, 2) &
524  = c(2, 3, 1, 2) &
525  +c_hat(i, j, k, l) &
526  *e_hat_dot_cg(i, 2)*e_hat_dot_cg(j, 3) &
527  *e_hat_dot_cg(k, 1)*e_hat_dot_cg(l, 2)
528 
529  end do
530 
531  end do
532 
533  end do
534 
535  do k = 1, 3
536 
537  do j = 1, 2
538 
539  do i = 1, 3
540 
541  c(3, 1, 3, 1) &
542  = c(3, 1, 3, 1) &
543  +c_hat(i, j, k, l) &
544  *e_hat_dot_cg(i, 3)*e_hat_dot_cg(j, 1) &
545  *e_hat_dot_cg(k, 3)*e_hat_dot_cg(l, 1)
546 
547  end do
548 
549  end do
550 
551  end do
552 
553  end do
554 
555  do l = 1, 3
556 
557  do k = 1, 2
558 
559  do j = 1, 2
560 
561  do i = 1, 3
562 
563  c(3, 1, 2, 3) &
564  = c(3, 1, 2, 3) &
565  +c_hat(i, j, k, l) &
566  *e_hat_dot_cg(i, 3)*e_hat_dot_cg(j, 1) &
567  *e_hat_dot_cg(k, 2)*e_hat_dot_cg(l, 3)
568 
569  end do
570 
571  end do
572 
573  do j = 1, 3
574 
575  do i = 1, 2
576 
577  c(2, 3, 2, 3) &
578  = c(2, 3, 2, 3) &
579  +c_hat(i, j, k, l) &
580  *e_hat_dot_cg(i, 2)*e_hat_dot_cg(j, 3) &
581  *e_hat_dot_cg(k, 2)*e_hat_dot_cg(l, 3)
582 
583  end do
584 
585  end do
586 
587  end do
588 
589  end do
590 
591  c(1, 1, 2, 2) = c(2, 2, 1, 1)
592  c(1, 1, 1, 2) = c(1, 2, 1, 1)
593  c(1, 1, 2, 3) = c(2, 3, 1, 1)
594  c(1, 1, 3, 1) = c(3, 1, 1, 1)
595  c(2, 2, 1, 2) = c(1, 2, 2, 2)
596  c(2, 2, 2, 3) = c(2, 3, 2, 2)
597  c(2, 2, 3, 1) = c(3, 1, 2, 2)
598  c(1, 2, 2, 3) = c(2, 3, 1, 2)
599  c(1, 2, 3, 1) = c(3, 1, 1, 2)
600  c(2, 3, 3, 1) = c(3, 1, 2, 3)
601 
602  !--------------------------------------------------------
603 
604  ! DO l = 1, 3
605  !
606  ! DO k = 1, 3
607  !
608  ! DO j = 1, 3
609  !
610  ! DO i = 1, 3
611  !
612  ! c(i, j, k, l) = 0.0D0
613  !
614  ! DO ll = 1, 3
615  !
616  ! DO kk = 1, 3
617  !
618  ! DO jj = 1, 3
619  !
620  ! DO ii = 1, 3
621  !
622  ! c(i, j, k, l) &
623  ! = c(i, j, k, l) &
624  ! +c_hat(ii, jj, kk, ll) &
625  ! *e_hat_dot_cg(ii, i)*e_hat_dot_cg(jj, j) &
626  ! *e_hat_dot_cg(kk, k)*e_hat_dot_cg(ll, l)
627  !
628  ! END DO
629  !
630  ! END DO
631  !
632  ! END DO
633  !
634  ! END DO
635  !
636  ! END DO
637  !
638  ! END DO
639  !
640  ! END DO
641  !
642  ! END DO
643 
644  !--------------------------------------------------------
645 
646  end select
647 
648  !--------------------------------------------------------------------
649 
650  return
651 
652  !####################################################################
653  end subroutine linearelastic_shell
654  !####################################################################
655  ! > (Gaku Hashimoto, The University of Tokyo, 2012/11/15)
656 
657 
658 
659 end module m_elasticlinear
mcontact::mu
real(kind=kreal), save mu
penalty, default value
Definition: fstr_contact.f90:17
m_elasticlinear::calelasticmatrix_ortho
subroutine calelasticmatrix_ortho(matl, sectType, bij, DMAT, temp)
Calculate orthotropic elastic matrix.
Definition: ElasticLinear.f90:109
mmaterial::mc_isoelastic
character(len=dict_key_length) mc_isoelastic
Definition: material.f90:140
m_utilities
This module provides aux functions.
Definition: utilities.f90:6
mmaterial::shell
integer(kind=kint), parameter shell
Definition: material.f90:88
mmaterial::axissymetric
integer(kind=kint), parameter axissymetric
Definition: material.f90:87
mmaterial::m_youngs
integer(kind=kint), parameter m_youngs
Definition: material.f90:92
mmaterial::m_poisson
integer(kind=kint), parameter m_poisson
Definition: material.f90:93
mmaterial::m_alpha_over_mu
integer(kind=kint), parameter m_alpha_over_mu
Definition: material.f90:107
hecmw_util
I/O and Utility.
Definition: hecmw_util_f.F90:7
m_utilities::transformation
subroutine transformation(jacob, tm)
Definition: utilities.f90:339
mmaterial::mc_orthoelastic
character(len=dict_key_length) mc_orthoelastic
Definition: material.f90:141
mmaterial::d3
integer(kind=kint), parameter d3
Definition: material.f90:84
hecmw_util::kreal
integer(kind=4), parameter kreal
Definition: hecmw_util_f.F90:16
m_elasticlinear
This module provides functions for elastic material.
Definition: ElasticLinear.f90:6
mmaterial::tmaterial
Structure to manage all material related data.
Definition: material.f90:167
mmaterial::planestrain
integer(kind=kint), parameter planestrain
Definition: material.f90:86
mmaterial
This module summarizes all information of material properties.
Definition: material.f90:6
m_elasticlinear::linearelastic_shell
subroutine linearelastic_shell(matl, sectType, c, e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, alpha, n_layer)
Definition: ElasticLinear.f90:168
m_elasticlinear::calelasticmatrix
subroutine calelasticmatrix(matl, sectType, D, temp, hdflag)
Calculate isotropic elastic matrix.
Definition: ElasticLinear.f90:16
mmaterial::planestress
integer(kind=kint), parameter planestress
Definition: material.f90:85