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