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