FrontISTR  6.0.0
Large-scale structural analysis program with finit element method
fstr_CreateMatrix.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 !-------------------------------------------------------------------------------
9 
10 module m_fstr_creatematrix_and_dampingforce
11  use m_fstr
12  implicit none
13 
14  private
15  public :: fstr_creatematrix_and_dampingforce
16  public :: fstr_calc_kinetic_energy
17  public :: fstr_calc_mass_product
18  public :: fstr_report_lumped_mass_fallback
19 
21  integer(kind=kint), parameter :: kLumpedMassFallbackTypes(5) = (/ 301, 611, 641, 761, 781 /)
22 
23 contains
24 
25  !---------------------------------------------------------------------*
36  subroutine fstr_creatematrix_and_dampingforce( hecMESH, hecMAT, fstrSOLID, time, tincr, fstrDYNAMIC, coef )
37  !---------------------------------------------------------------------*
38  use m_static_lib
39  use mmechgauss
40  use m_dynamic_mass
41  use m_elemact
42  use m_fstr_nodalkinematics, only: fstr_ensure_finite_rotation_state
43  use m_fstr_finiterotationkinematics, only: fstr_uses_finite_rotation_kinematics, &
44  shellcomposenodaldisplacement
45 
46  type(hecmwst_local_mesh) :: hecmesh
47  type(hecmwst_matrix) :: hecmat
48  type(fstr_solid) :: fstrsolid
49  real(kind=kreal), intent(in) :: time
50  real(kind=kreal), intent(in) :: tincr
51  type(fstr_dynamic), intent(in), optional :: fstrdynamic
52  real(kind=kreal), intent(in), optional :: coef(6)
53 
54  integer(kind=kint) :: ndof, itype, is, ie, ic_type, nn, icel, iis, i, j, in, jn
55  integer(kind=kint) :: nodlocal(fstrsolid%max_ncon)
56  real(kind=kreal) :: stiff_mat(fstrsolid%max_ncon_stf*6, fstrsolid%max_ncon_stf*6)
57  real(kind=kreal) :: mass_mat(fstrsolid%max_ncon_stf*6, fstrsolid%max_ncon_stf*6)
58  real(kind=kreal) :: damp_mat(fstrsolid%max_ncon_stf*6, fstrsolid%max_ncon_stf*6)
59  real(kind=kreal) :: mat(fstrsolid%max_ncon_stf*6, fstrsolid%max_ncon_stf*6), lumped(fstrsolid%max_ncon_stf*6)
60  real(kind=kreal) :: tt(fstrsolid%max_ncon), ecoord(3,fstrsolid%max_ncon)
61  real(kind=kreal) :: u(6,fstrsolid%max_ncon), u_inc(6,fstrsolid%max_ncon)
62  real(kind=kreal) :: du(fstrsolid%max_ncon_stf*6), u_prev(6,fstrsolid%max_ncon)
63  real(kind=kreal) :: triad_tri(9,fstrsolid%max_ncon)
64  real(kind=kreal) :: triad_cur(9,fstrsolid%max_ncon)
65  real(kind=kreal) :: triad_ref(9,fstrsolid%max_ncon)
66  real(kind=kreal) :: shell_drill(fstrsolid%max_ncon)
67  real(kind=kreal) :: veca(fstrsolid%max_ncon_stf*6), vecb(fstrsolid%max_ncon_stf*6)
68  real(kind=kreal) :: acc(fstrsolid%max_ncon_stf*6), vec(fstrsolid%max_ncon_stf*6), df(fstrsolid%max_ncon_stf*6)
69  real(kind=kreal) :: a1, a2, a3, b1, b2, b3, wm, thick
70  type(tmaterial), pointer :: material
71  logical :: is_dynamic, use_lumped_mass
72 
73  is_dynamic = present(fstrdynamic) .and. present(coef)
74  use_lumped_mass = .true.
75  if( is_dynamic ) use_lumped_mass = (fstrdynamic%idx_mas == kmasslumped)
76  a1 = 0.0d0; a2 = 0.0d0; a3 = 0.0d0
77  b1 = 0.0d0; b2 = 0.0d0; b3 = 0.0d0
78  wm = 1.0d0
79 
80  call hecmw_mat_clear( hecmat )
81  if( is_dynamic ) then
82  fstrsolid%DFORCE = 0.0d0
83  a1 = coef(1); a2 = coef(2); a3 = coef(3)
84  b1 = coef(4); b2 = coef(5); b3 = coef(6)
85  ! HHT: the equation of motion is divided by 1+alpha, which leaves only the inertia scaled
86  if( fstrdynamic%idx_eqa == 2 ) wm = 1.0d0/(1.0d0+fstrdynamic%alpha)
87  endif
88 
89  ndof = hecmat%NDOF
90  call fstr_ensure_finite_rotation_state( hecmesh, fstrsolid, ndof )
91 
92  do itype = 1, hecmesh%n_elem_type
93  is = hecmesh%elem_type_index(itype-1) + 1
94  ie = hecmesh%elem_type_index(itype )
95  ic_type = hecmesh%elem_type_item(itype)
96 
97  ! ----- Ignore link and patch elements
98  if (hecmw_is_etype_link(ic_type)) cycle
99  if (hecmw_is_etype_patch(ic_type)) cycle
100 
101  nn = hecmw_get_max_node(ic_type)
102 
103  !$omp parallel default(none), &
104  !$omp& private(icel,iiS,nn,j,nodLOCAL,i,in,jn,ecoord,du,u,u_inc,u_prev,tt, &
105  !$omp& triad_tri,triad_cur,triad_ref,shell_drill, &
106  !$omp& material,thick,stiff_mat,mass_mat,damp_mat, &
107  !$omp& lumped,mat,df,vec,acc,vecA,vecB), &
108  !$omp& shared(iS,iE,hecMESH,ndof,fstrSOLID,ic_type,hecMAT,time,tincr,fstrDYNAMIC, &
109  !$omp& is_dynamic,use_lumped_mass,a1,a2,a3,b1,b2,b3,wm)
110  !$omp do
111  do icel = is, ie
112 
113  ! ----- nodal coordinate & displacement
114  iis = hecmesh%elem_node_index(icel-1)
115  nn = hecmesh%elem_node_index(icel) - iis
116 
117  ! ----- rows beyond nn_internal are not equations, so an element with no internal node assembles nothing used here
118  if( .not. hecmw_has_internal_node( hecmesh, nn, hecmesh%elem_node_item(iis+1:iis+nn) ) ) cycle
119 
120  tt(:) = 0.d0
121  do j = 1, nn
122  in = hecmesh%elem_node_item(iis+j)
123  nodlocal(j) = in
124  do i = 1, 3
125  ecoord(i,j) = hecmesh%node(3*(in-1)+i)
126  enddo
127  do i = 1, ndof
128  u_inc(i,j) = fstrsolid%dunode(ndof*(in-1)+i)
129  u_prev(i,j) = fstrsolid%unode(ndof*(in-1)+i)
130  u(i,j) = u_prev(i,j)+u_inc(i,j)
131  enddo
132  if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0 ) then
133  tt(j) = fstrsolid%temperature( nodlocal(j) )
134  endif
135  if( ic_type == 741 .or. ic_type == 743 .or. ic_type == 731 ) then
136  triad_tri(1:9,j) = 0.0d0
137  triad_cur(1:9,j) = 0.0d0
138  triad_ref(1:9,j) = 0.0d0
139  shell_drill(j) = 0.0d0
140  if( associated(fstrsolid%shell_dtriad) ) triad_tri(1:9,j) = fstrsolid%shell_dtriad(9*(in-1)+1:9*(in-1)+9)
141  if( associated(fstrsolid%shell_triad) ) triad_cur(1:9,j) = fstrsolid%shell_triad(9*(in-1)+1:9*(in-1)+9)
142  if( associated(fstrsolid%shell_ref_triad) ) triad_ref(1:9,j) = fstrsolid%shell_ref_triad(9*(in-1)+1:9*(in-1)+9)
143  if( associated(fstrsolid%shell_ddrill) ) shell_drill(j) = fstrsolid%shell_ddrill(in)
144  endif
145  enddo
146 
147  if( ic_type == 761 .or. ic_type == 781 ) then
148  do j = 1, nn/2
149  in = nodlocal(j+nn/2)
150  triad_tri(1:9,j) = 0.0d0
151  triad_cur(1:9,j) = 0.0d0
152  triad_ref(1:9,j) = 0.0d0
153  shell_drill(j) = 0.0d0
154  if( associated(fstrsolid%shell_dtriad) ) triad_tri(1:9,j) = fstrsolid%shell_dtriad(9*(in-1)+1:9*(in-1)+9)
155  if( associated(fstrsolid%shell_triad) ) triad_cur(1:9,j) = fstrsolid%shell_triad(9*(in-1)+1:9*(in-1)+9)
156  if( associated(fstrsolid%shell_ref_triad) ) triad_ref(1:9,j) = fstrsolid%shell_ref_triad(9*(in-1)+1:9*(in-1)+9)
157  if( associated(fstrsolid%shell_ddrill) ) shell_drill(j) = fstrsolid%shell_ddrill(in)
158  enddo
159  endif
160 
161  if( ic_type == 741 .or. ic_type == 743 .or. ic_type == 731 ) then
162  material => fstrsolid%elements(icel)%gausses(1)%pMaterial
163  if( fstr_uses_finite_rotation_kinematics(ic_type, nn, material) ) then
164  call shellcomposenodaldisplacement(ndof, nn, u_prev(1:6,1:nn), u_inc(1:6,1:nn), u(1:6,1:nn))
165  endif
166  endif
167 
168  ! ----- inactive element : assemble dummy stiffness and skip the rest
169  if( fstrsolid%elements(icel)%elemact_flag == kelact_inactive ) then
170  stiff_mat = 0.0d0
171  call stf_dummy( ndof, nn, ecoord(:,1:nn), &
172  & stiff_mat(1:nn*ndof, 1:nn*ndof), fstrsolid%elements(icel) )
173  call hecmw_mat_ass_elem(hecmat, nn, nodlocal, stiff_mat)
174  cycle
175  endif
176 
177  ! ----- element K, M, C (section / material context resolved inside)
178  stiff_mat = 0.0d0
179  mass_mat = 0.0d0
180  damp_mat = 0.0d0
181 
182  call calc_stiff_and_mass_elem( ic_type, nn, ndof, ecoord, u, u_prev, tt, &
183  time, tincr, is_dynamic, fstrsolid, hecmesh, icel, nodlocal, &
184  triad_tri, triad_cur, triad_ref, shell_drill, &
185  stiff_mat, mass_mat, damp_mat, lumped, use_lumped_mass )
186 
187  ! ----- assemble system matrix (and the dynamic damping force)
188  if( is_dynamic ) then
189  ! gathered with nn and nodLOCAL returned by calc_stiff_and_mass_elem, which compresses them for smoothed elements
190  do j = 1, nn
191  in = nodlocal(j)
192  do i = 1, ndof
193  du(ndof*(j-1)+i) = fstrsolid%dunode(ndof*(in-1)+i)
194  vec(ndof*(j-1)+i) = fstrdynamic%VEL(ndof*(in-1)+i,1)
195  acc(ndof*(j-1)+i) = fstrdynamic%ACC(ndof*(in-1)+i,1)
196  enddo
197  enddo
198 
199  ! Newmark intermediate vectors
200  do i = 1, nn*ndof
201  veca(i) = -a3*du(i) + a2*vec(i) + a1*acc(i)
202  vecb(i) = -b3*du(i) + b2*vec(i) + b1*acc(i)
203  enddo
204 
205  ! effective dynamic system matrix A = c1*K + c2*M (+ b3*C) and damping force df
206  material => fstrsolid%elements(icel)%gausses(1)%pMaterial
207  call calc_damping_mat_and_force_elem( ic_type, nn, ndof, material, fstrdynamic, &
208  a3, b3, wm, veca, vecb, stiff_mat, mass_mat, damp_mat, mat, df )
209 
210  call hecmw_mat_ass_elem(hecmat, nn, nodlocal, mat)
211 
212  do j = 1, nn
213  do i = 1, ndof
214  !$omp atomic
215  fstrsolid%DFORCE(ndof*(nodlocal(j)-1)+i) = fstrsolid%DFORCE(ndof*(nodlocal(j)-1)+i)+df(ndof*(j-1)+i)
216  enddo
217  enddo
218  else
219  ! static / eigen / frequency : assemble K only
220  call hecmw_mat_ass_elem(hecmat, nn, nodlocal, stiff_mat)
221  endif
222 
223  enddo ! icel
224  !$omp end do
225  !$omp end parallel
226  enddo ! itype
227 
228  end subroutine fstr_creatematrix_and_dampingforce
229 
232  real(kind=kreal) function fstr_calc_kinetic_energy( &
233  hecMESH, fstrSOLID, fstrDYNAMIC, velocity) result(kinetic_energy)
234  use m_elemact
235 
236  type(hecmwst_local_mesh), intent(in) :: hecmesh
237  type(fstr_solid), intent(in) :: fstrsolid
238  type(fstr_dynamic), intent(in) :: fstrdynamic
239  real(kind=kreal), intent(in) :: velocity(:)
240 
241  integer(kind=kint) :: ndof, itype, is, ie, ic_type, icel, iis, nn
242  integer(kind=kint) :: i, j, row, in, nsize
243  integer(kind=kint) :: nodlocal(fstrsolid%max_ncon)
244  real(kind=kreal) :: ecoord(3,fstrsolid%max_ncon)
245  real(kind=kreal) :: mass_mat(20*6,20*6), lumped(20*6), elem_velocity(20*6)
246  logical :: use_lumped_mass
247 
248  ndof = hecmesh%n_dof
249  use_lumped_mass = (fstrdynamic%idx_mas == kmasslumped)
250  kinetic_energy = 0.0d0
251 
252  do itype = 1, hecmesh%n_elem_type
253  is = hecmesh%elem_type_index(itype-1) + 1
254  ie = hecmesh%elem_type_index(itype)
255  ic_type = hecmesh%elem_type_item(itype)
256 
257  if( hecmw_is_etype_link(ic_type) ) cycle
258  if( hecmw_is_etype_patch(ic_type) ) cycle
259  if( hecmw_is_etype_smoothing(ic_type) ) cycle
260 
261  !$omp parallel default(none), &
262  !$omp& private(icel,iiS,nn,i,j,row,in,nsize,nodLOCAL,ecoord,mass_mat,lumped,elem_velocity), &
263  !$omp& shared(iS,iE,ic_type,ndof,hecMESH,fstrSOLID,velocity,use_lumped_mass) &
264  !$omp& reduction(+:kinetic_energy)
265  !$omp do
266  do icel = is, ie
267  iis = hecmesh%elem_node_index(icel-1)
268  nn = hecmesh%elem_node_index(icel) - iis
269  if( .not.hecmw_has_internal_node(hecmesh, nn, hecmesh%elem_node_item(iis+1:iis+nn)) ) cycle
270  if( fstrsolid%elements(icel)%elemact_flag == kelact_inactive ) cycle
271 
272  ecoord = 0.0d0
273  elem_velocity = 0.0d0
274  do j = 1, nn
275  in = hecmesh%elem_node_item(iis+j)
276  nodlocal(j) = in
277  ecoord(1:3,j) = hecmesh%node(3*(in-1)+1:3*in)
278  elem_velocity(ndof*(j-1)+1:ndof*j) = velocity(ndof*(in-1)+1:ndof*in)
279  enddo
280 
281  mass_mat = 0.0d0
282  lumped = 0.0d0
283  call calc_mass_elem(ic_type, nn, ndof, ecoord, fstrsolid, hecmesh, icel, &
284  use_lumped_mass, mass_mat, lumped)
285 
286  nsize = nn*ndof
287  do j = 1, nn
288  if( nodlocal(j) > hecmesh%nn_internal ) cycle
289  do i = 1, ndof
290  row = ndof*(j-1)+i
291  kinetic_energy = kinetic_energy + &
292  elem_velocity(row)*dot_product(mass_mat(row,1:nsize), elem_velocity(1:nsize))
293  enddo
294  enddo
295  enddo
296  !$omp end do
297  !$omp end parallel
298  enddo
299 
300  call hecmw_allreduce_r1(hecmesh, kinetic_energy, hecmw_sum)
301  kinetic_energy = 0.5d0*kinetic_energy
302  end function fstr_calc_kinetic_energy
303 
306  subroutine fstr_calc_mass_product( hecMESH, fstrSOLID, fstrDYNAMIC, vec, mass_vec )
307  use m_elemact
308 
309  type(hecmwst_local_mesh), intent(in) :: hecmesh
310  type(fstr_solid), intent(in) :: fstrsolid
311  type(fstr_dynamic), intent(in) :: fstrdynamic
312  real(kind=kreal), intent(in) :: vec(:)
313  real(kind=kreal), intent(out) :: mass_vec(:)
314 
315  integer(kind=kint) :: ndof, itype, is, ie, ic_type, icel, iis, nn
316  integer(kind=kint) :: i, j, row, in, nsize
317  integer(kind=kint) :: nodlocal(fstrsolid%max_ncon)
318  real(kind=kreal) :: ecoord(3,fstrsolid%max_ncon)
319  real(kind=kreal) :: mass_mat(20*6,20*6), lumped(20*6), elem_vec(20*6), prod
320  logical :: use_lumped_mass
321 
322  ndof = hecmesh%n_dof
323  use_lumped_mass = (fstrdynamic%idx_mas == kmasslumped)
324  mass_vec(:) = 0.0d0
325 
326  do itype = 1, hecmesh%n_elem_type
327  is = hecmesh%elem_type_index(itype-1) + 1
328  ie = hecmesh%elem_type_index(itype)
329  ic_type = hecmesh%elem_type_item(itype)
330 
331  if( hecmw_is_etype_link(ic_type) ) cycle
332  if( hecmw_is_etype_patch(ic_type) ) cycle
333  if( hecmw_is_etype_smoothing(ic_type) ) cycle
334 
335  !$omp parallel default(none), &
336  !$omp& private(icel,iiS,nn,i,j,row,in,nsize,nodLOCAL,ecoord,mass_mat,lumped,elem_vec,prod), &
337  !$omp& shared(iS,iE,ic_type,ndof,hecMESH,fstrSOLID,vec,mass_vec,use_lumped_mass)
338  !$omp do
339  do icel = is, ie
340  iis = hecmesh%elem_node_index(icel-1)
341  nn = hecmesh%elem_node_index(icel) - iis
342  if( .not.hecmw_has_internal_node(hecmesh, nn, hecmesh%elem_node_item(iis+1:iis+nn)) ) cycle
343  if( fstrsolid%elements(icel)%elemact_flag == kelact_inactive ) cycle
344 
345  ecoord = 0.0d0
346  elem_vec = 0.0d0
347  do j = 1, nn
348  in = hecmesh%elem_node_item(iis+j)
349  nodlocal(j) = in
350  ecoord(1:3,j) = hecmesh%node(3*(in-1)+1:3*in)
351  elem_vec(ndof*(j-1)+1:ndof*j) = vec(ndof*(in-1)+1:ndof*in)
352  enddo
353 
354  mass_mat = 0.0d0
355  lumped = 0.0d0
356  call calc_mass_elem(ic_type, nn, ndof, ecoord, fstrsolid, hecmesh, icel, &
357  use_lumped_mass, mass_mat, lumped)
358 
359  nsize = nn*ndof
360  do j = 1, nn
361  in = nodlocal(j)
362  if( in > hecmesh%nn_internal ) cycle
363  do i = 1, ndof
364  row = ndof*(j-1)+i
365  prod = dot_product(mass_mat(row,1:nsize), elem_vec(1:nsize))
366  !$omp atomic
367  mass_vec(ndof*(in-1)+i) = mass_vec(ndof*(in-1)+i) + prod
368  enddo
369  enddo
370  enddo
371  !$omp end do
372  !$omp end parallel
373  enddo
374  end subroutine fstr_calc_mass_product
375 
378  subroutine fstr_report_lumped_mass_fallback( hecMESH, fstrDYNAMIC )
379  type(hecmwst_local_mesh), intent(in) :: hecmesh
380  type(fstr_dynamic), intent(in) :: fstrdynamic
381 
382  integer(kind=kint) :: itype, i
383  integer(kind=kint) :: found(size(klumpedmassfallbacktypes))
384 
385  if( fstrdynamic%idx_mas /= kmassconsistent ) return
386 
387  found = 0
388  do itype = 1, hecmesh%n_elem_type
389  ! a distributed mesh may list element types that have no element in this subdomain
390  if( hecmesh%elem_type_index(itype) == hecmesh%elem_type_index(itype-1) ) cycle
391  do i = 1, size(klumpedmassfallbacktypes)
392  if( hecmesh%elem_type_item(itype) == klumpedmassfallbacktypes(i) ) found(i) = 1
393  enddo
394  enddo
395  call hecmw_allreduce_i(hecmesh, found, size(found), hecmw_max)
396 
397  if( hecmesh%my_rank /= 0 ) return
398  do i = 1, size(klumpedmassfallbacktypes)
399  if( found(i) == 0 ) cycle
400  write(imsg,'(a,i0,a)') 'INFO: consistent mass is unavailable for element type ', &
401  klumpedmassfallbacktypes(i), '; lumped mass is used'
402  enddo
403  end subroutine fstr_report_lumped_mass_fallback
404 
405  !---------------------------------------------------------------------*
414  subroutine calc_stiff_and_mass_elem( ic_type, nn, ndof, ecoord, u, u_prev, tt, &
415  time, tincr, is_dynamic, fstrSOLID, hecMESH, icel, nodLOCAL, &
416  triad_tri, triad_cur, triad_ref, shell_drill, &
417  stiff_mat, mass_mat, damp_mat, lumped, use_lumped_mass )
418  !---------------------------------------------------------------------*
419  use m_static_lib
420  use mmechgauss
421 
422  integer(kind=kint), intent(in) :: ic_type, ndof
423  integer(kind=kint), intent(inout) :: nn ! some legacy STF_* routines mutate this
424  real(kind=kreal), intent(in) :: ecoord(:,:), u(:,:), u_prev(:,:), tt(:)
425  real(kind=kreal), intent(in) :: time, tincr
426  logical, intent(in) :: is_dynamic
427  type(fstr_solid), intent(inout) :: fstrsolid
428  type(hecmwst_local_mesh), intent(in) :: hecmesh
429  integer(kind=kint), intent(in) :: icel
430  integer(kind=kint), intent(inout) :: nodlocal(:)
431  real(kind=kreal), intent(in) :: triad_tri(:,:), triad_cur(:,:)
432  real(kind=kreal), intent(in) :: triad_ref(:,:), shell_drill(:)
433  real(kind=kreal), intent(inout) :: stiff_mat(:,:), mass_mat(:,:), damp_mat(:,:)
434  real(kind=kreal), intent(inout) :: lumped(:)
435  logical, intent(in) :: use_lumped_mass
436 
437  type(tmaterial), pointer :: material
438  integer(kind=kint) :: isect, ihead, cdsys_id
439  real(kind=kreal) :: coords(3,3), thick
440 
441  ! ----- section / material context
442  isect = hecmesh%section_ID(icel)
443  ihead = hecmesh%section%sect_R_index(isect-1)
444  cdsys_id = hecmesh%section%sect_orien_ID(isect)
445  coords = 0.0d0
446  if( cdsys_id > 0 ) call get_coordsys(cdsys_id, hecmesh, fstrsolid, coords, icel)
447 
448  material => fstrsolid%elements(icel)%gausses(1)%pMaterial
449  thick = hecmesh%section%sect_R_item(ihead+1)
450 
451  if( ic_type==241 .or. ic_type==242 .or. ic_type==231 .or. ic_type==232 .or. ic_type==2322) then
452  if( material%nlgeom_flag /= infinitesimal ) call createmat_abort( ic_type, 2 )
453  call stf_c2( ic_type, nn, ecoord(1:2,1:nn), fstrsolid%elements(icel)%gausses(:), thick, &
454  stiff_mat(1:nn*ndof,1:nn*ndof), fstrsolid%elements(icel)%iset, u(1:2,1:nn) )
455 
456  elseif( ic_type==301 ) then
457  call stf_c1( nn, ecoord(:,1:nn), thick, fstrsolid%elements(icel)%gausses(:), &
458  stiff_mat(1:nn*ndof,1:nn*ndof), u(1:3,1:nn) )
459 
460  elseif( ic_type==361 ) then
461  if( fstrsolid%sections(isect)%elemopt361 == kel361fi ) then ! full integration element
462  call stf_c3( ic_type, nn, ecoord(:,1:nn), fstrsolid%elements(icel)%gausses(:), &
463  stiff_mat(1:nn*ndof,1:nn*ndof), cdsys_id, coords, time, tincr, u(1:3,1:nn), tt(1:nn) )
464  else if( fstrsolid%sections(isect)%elemopt361 == kel361bbar ) then ! B-bar element
465  call stf_c3d8bbar( ic_type, nn, ecoord(:,1:nn), fstrsolid%elements(icel)%gausses(:), &
466  stiff_mat(1:nn*ndof,1:nn*ndof), cdsys_id, coords, time, tincr, u(1:3,1:nn), tt(1:nn) )
467  else if( fstrsolid%sections(isect)%elemopt361 == kel361ic ) then ! incompatible element
468  call stf_c3d8ic( ic_type, nn, ecoord(:,1:nn), fstrsolid%elements(icel)%gausses(:), &
469  stiff_mat(1:nn*ndof,1:nn*ndof), cdsys_id, coords, time, tincr, u(1:3,1:nn), &
470  fstrsolid%elements(icel)%aux, tt(1:nn) )
471  else if( fstrsolid%sections(isect)%elemopt361 == kel361fbar ) then ! F-bar element
472  call stf_c3d8fbar( ic_type, nn, ecoord(:,1:nn), fstrsolid%elements(icel)%gausses(:), &
473  stiff_mat(1:nn*ndof,1:nn*ndof), cdsys_id, coords, time, tincr, u(1:3,1:nn), tt(1:nn) )
474  else if( fstrsolid%sections(isect)%elemopt361 == kel361up ) then ! UP (u-p mixed) element
475  call stf_c3_up( ic_type, nn, ecoord(:,1:nn), fstrsolid%elements(icel)%gausses(:), &
476  stiff_mat(1:nn*ndof,1:nn*ndof), cdsys_id, coords, time, tincr, &
477  1, u(1:3,1:nn), tt(1:nn) )
478  endif
479 
480  elseif( ic_type==341 .or. ic_type==351 .or. ic_type==342 .or. ic_type==352 .or. ic_type==362 ) then
481  ! SESNS option: stiffness is assembled separately via ic_type==881/891;
482  ! return stiff_mat = 0 here and still compute the element mass below.
483  if( .not. (ic_type==341 .and. fstrsolid%sections(isect)%elemopt341 == kel341sesns) ) then
484  call stf_c3( ic_type, nn, ecoord(:,1:nn), fstrsolid%elements(icel)%gausses(:), &
485  stiff_mat(1:nn*ndof,1:nn*ndof), cdsys_id, coords, time, tincr, u(1:3,1:nn), tt(1:nn) )
486  endif
487 
488  else if( ic_type == 511 ) then
489  call stf_connector( nn, ecoord(:,1:nn), fstrsolid%elements(icel)%gausses(:), &
490  stiff_mat(1:nn*ndof,1:nn*ndof), u(1:3,1:nn))
491 
492  if( is_dynamic ) call dmp_connector( nn, ecoord(:,1:nn), fstrsolid%elements(icel)%gausses(:), &
493  damp_mat(1:nn*ndof,1:nn*ndof), u(1:3,1:nn))
494 
495  else if( ic_type == 611 ) then
496  if( material%nlgeom_flag /= infinitesimal ) call createmat_abort( ic_type, 2 )
497  call stf_beam(nn, ecoord, hecmesh%section%sect_R_item(ihead+1:), &
498  material%variables(m_youngs), material%variables(m_poisson), stiff_mat(1:nn*ndof,1:nn*ndof), &
499  fstrsolid%sections(isect)%elemopt611)
500 
501  else if( ic_type == 641 ) then
502  if( material%nlgeom_flag /= infinitesimal ) call createmat_abort( ic_type, 2 )
503  call stf_beam_641_from_611(ecoord(1:3,1:4), fstrsolid%elements(icel)%gausses(:), &
504  hecmesh%section%sect_R_item(ihead+1:), stiff_mat(1:nn*ndof,1:nn*ndof), &
505  fstrsolid%sections(isect)%elemopt611)
506 
507  else if( ( ic_type == 741 ) .or. ( ic_type == 743 ) .or. ( ic_type == 731 ) ) then
508  call stf_shell_mitc(ic_type, nn, ndof, ecoord(1:3,1:nn), fstrsolid%elements(icel)%gausses(:), &
509  stiff_mat(1:nn*ndof,1:nn*ndof), thick, 0, nddisp=u(1:ndof,1:nn), &
510  element=fstrsolid%elements(icel), ndtriad=triad_tri(1:9,1:nn), &
511  ndreftriad=triad_ref(1:9,1:nn), &
512  ndcurtriad=triad_cur(1:9,1:nn), nddrill=shell_drill(1:nn), cdsys_id=cdsys_id, coords=coords)
513 
514  else if( ic_type == 761 ) then ! for shell-solid mixed analysis
515  if( material%nlgeom_flag == infinitesimal ) then
516  call stf_shell_mitc(731, 3, 6, ecoord(1:3,1:3), fstrsolid%elements(icel)%gausses(:), &
517  stiff_mat(1:nn*ndof,1:nn*ndof), thick, 2, cdsys_id=cdsys_id, coords=coords)
518  else
519  call stf_shell_mitc33(731, 3, 6, ecoord(1:3,1:3), u(1:3,1:6), u_prev(1:3,1:6), &
520  fstrsolid%elements(icel)%gausses(:), stiff_mat(1:nn*ndof,1:nn*ndof), thick, 2, &
521  fstrsolid%elements(icel), triad_tri(1:9,1:3), triad_ref(1:9,1:3), &
522  triad_cur(1:9,1:3), shell_drill(1:3), cdsys_id=cdsys_id, coords=coords)
523  endif
524 
525  else if( ic_type == 781 ) then ! for shell-solid mixed analysis
526  call stf_shell_mitc33(741, 4, 6, ecoord(1:3,1:4), u(1:3,1:8), u_prev(1:3,1:8), &
527  fstrsolid%elements(icel)%gausses(:), stiff_mat(1:nn*ndof,1:nn*ndof), thick, 1, &
528  fstrsolid%elements(icel), triad_tri(1:9,1:4), triad_ref(1:9,1:4), &
529  triad_cur(1:9,1:4), shell_drill(1:4), cdsys_id=cdsys_id, coords=coords)
530 
531  elseif( ic_type==3414 ) then
532  if( material%mtype /= incomp_newtonian ) call createmat_abort( ic_type, 3, material%mtype )
533  call stf_c3_vp( ic_type, nn, ecoord(:,1:nn), fstrsolid%elements(icel)%gausses(:), &
534  stiff_mat(1:nn*ndof,1:nn*ndof), tincr, u_prev(1:4,1:nn) )
535 
536  else if( ic_type == 881 .or. ic_type == 891 ) then ! for selective es/ns smoothed fem
537  call stf_c3d4_sesns( ic_type, nn, nodlocal, ecoord(:,1:nn), fstrsolid%elements(icel)%gausses(:), &
538  stiff_mat, cdsys_id, coords, time, tincr, u(1:3,1:nn), tt(1:nn) )
539 
540  else
541  call createmat_abort( ic_type, 1 )
542  endif
543 
544  if( is_dynamic ) call calc_mass_elem( ic_type, nn, ndof, ecoord, fstrsolid, hecmesh, icel, &
545  use_lumped_mass, mass_mat, lumped )
546 
547  end subroutine calc_stiff_and_mass_elem
548 
549  !---------------------------------------------------------------------*
557  subroutine calc_mass_elem( ic_type, nn, ndof, ecoord, fstrSOLID, hecMESH, icel, &
558  use_lumped_mass, mass_mat, lumped )
559  !---------------------------------------------------------------------*
560  use m_static_lib
561  use mmechgauss
562  use m_dynamic_mass
563 
564  integer(kind=kint), intent(in) :: ic_type, nn, ndof, icel
565  real(kind=kreal), intent(in) :: ecoord(:,:)
566  type(fstr_solid), intent(in) :: fstrsolid
567  type(hecmwst_local_mesh), intent(in) :: hecMESH
568  logical, intent(in) :: use_lumped_mass
569  real(kind=kreal), intent(inout) :: mass_mat(:,:), lumped(:)
570 
571  type(tmaterial), pointer :: material
572  integer(kind=kint) :: i, isect, ihead, sec_opt
573  real(kind=kreal) :: thick, rho, length, surf
574 
575  isect = hecmesh%section_ID(icel)
576  ihead = hecmesh%section%sect_R_index(isect-1)
577  sec_opt = hecmesh%section%sect_opt(isect)
578 
579  material => fstrsolid%elements(icel)%gausses(1)%pMaterial
580  thick = hecmesh%section%sect_R_item(ihead+1)
581 
582  if( ic_type==241 .or. ic_type==242 .or. ic_type==231 .or. ic_type==232 .or. ic_type==2322) then
583  call mass_c2(ic_type, nn, ecoord(1:2,1:nn), fstrsolid%elements(icel)%gausses, &
584  sec_opt, thick, mass_mat, lumped, is_lumped=use_lumped_mass)
585 
586  elseif( ic_type==361 .or. ic_type==341 .or. ic_type==351 .or. ic_type==342 .or. ic_type==352 &
587  .or. ic_type==362 ) then
588  call mass_c3(ic_type, nn, ecoord(1:3,1:nn), fstrsolid%elements(icel)%gausses, &
589  mass_mat, lumped, is_lumped=use_lumped_mass)
590 
591  else if( ( ic_type == 741 ) .or. ( ic_type == 743 ) .or. ( ic_type == 731 ) ) then
592  rho = material%variables(m_density)
593  call mass_shell(ic_type, nn, ecoord(1:3,1:nn), rho, thick, fstrsolid%elements(icel)%gausses, mass_mat, lumped)
594  if( use_lumped_mass ) call set_diagonal_mass(nn*ndof, lumped, mass_mat)
595 
596  elseif( ic_type==301 ) then
597  surf = hecmesh%section%sect_R_item(ihead+1)
598  length = get_length(ecoord(1:3,1:nn))
599  rho = material%variables(m_density)
600  mass_mat = 0.0d0
601  do i = 1, nn*ndof
602  mass_mat(i,i) = 0.5d0*surf*length*rho
603  enddo
604 
605  else if( ic_type == 611 ) then
606  surf = hecmesh%section%sect_R_item(ihead+4)
607  length = get_length(ecoord(1:3,1:nn))
608  rho = material%variables(m_density)
609  call mass_beam(surf, length, rho, mass_mat)
610 
611  else if( ic_type == 641 ) then
612  surf = hecmesh%section%sect_R_item(ihead+4)
613  length = get_length(ecoord(1:3,1:nn))
614  rho = material%variables(m_density)
615  call mass_beam_33(surf, length, rho, mass_mat)
616 
617  else if( ic_type == 761 ) then ! for shell-solid mixed analysis
618  surf = get_face3(ecoord(1:3,1:nn))
619  rho = material%variables(m_density)
620  call mass_s3(surf, thick, rho, mass_mat)
621 
622  else if( ic_type == 781 ) then ! for shell-solid mixed analysis
623  surf = get_face4(ecoord(1:3,1:nn))
624  rho = material%variables(m_density)
625  call mass_s4(surf, thick, rho, mass_mat)
626 
627  else if( ic_type == 511 .or. ic_type == 3414 .or. ic_type == 881 .or. ic_type == 891 ) then
628  ! no element mass
629 
630  else
631  call createmat_abort( ic_type, 4 )
632  endif
633 
634  end subroutine calc_mass_elem
635 
636  subroutine set_diagonal_mass(nsize, lumped, mass)
637  integer(kind=kint), intent(in) :: nsize
638  real(kind=kreal), intent(in) :: lumped(:)
639  real(kind=kreal), intent(inout) :: mass(:,:)
640  integer(kind=kint) :: i
641 
642  mass = 0.0d0
643  do i = 1, nsize
644  mass(i,i) = lumped(i)
645  enddo
646  end subroutine set_diagonal_mass
647 
648  !---------------------------------------------------------------------*
656  subroutine calc_damping_mat_and_force_elem( ic_type, nn, ndof, material, fstrDYNAMIC, &
657  a3, b3, wm, vecA, vecB, stiff_mat, mass_mat, damp_mat, mat, df )
658  !---------------------------------------------------------------------*
659  integer(kind=kint), intent(in) :: ic_type, nn, ndof
660  type(tmaterial), pointer, intent(in) :: material
661  type(fstr_dynamic), intent(in) :: fstrDYNAMIC
662  real(kind=kreal), intent(in) :: a3, b3, wm
663  real(kind=kreal), intent(in) :: veca(:), vecb(:)
664  real(kind=kreal), intent(in) :: stiff_mat(:,:), mass_mat(:,:), damp_mat(:,:)
665  real(kind=kreal), intent(out) :: mat(:,:), df(:)
666 
667  integer(kind=kint) :: i, j
668  real(kind=kreal) :: ray_m, ray_k, c1, c2
669  real(kind=kreal) :: vecc(nn*ndof), kb(nn*ndof)
670 
671  if( material%is_elem_Rayleigh_damping ) then
672  ray_m = material%variables(m_damping_rm)
673  ray_k = material%variables(m_damping_rk)
674  else
675  ray_m = fstrdynamic%ray_m
676  ray_k = fstrdynamic%ray_k
677  endif
678 
679  ! ----- LHS: A = c1*K + c2*M (+ b3*C for connector)
680  c1 = 1.d0 + ray_k*b3
681  c2 = wm*a3 + ray_m*b3
682  do i = 1, nn*ndof
683  do j = 1, nn*ndof
684  mat(j,i) = c1*stiff_mat(j,i) + c2*mass_mat(j,i)
685  enddo
686  enddo
687  if( ic_type == 511 ) then
688  do i = 1, nn*ndof
689  do j = 1, nn*ndof
690  mat(j,i) = mat(j,i) + b3*damp_mat(j,i)
691  enddo
692  enddo
693  endif
694 
695  ! ----- RHS: df = M*(wm*vecA + ray_m*vecB) + ray_k*K*vecB (+ C*vecB for connector)
696  do i = 1, nn*ndof
697  vecc(i) = wm*veca(i) + ray_m*vecb(i)
698  enddo
699 
700  kb = 0.0d0
701  do i = 1, nn*ndof
702  do j = 1, nn*ndof
703  kb(i) = kb(i) + stiff_mat(i,j)*vecb(j)
704  enddo
705  enddo
706 
707  df = 0.0d0
708  df(1:nn*ndof) = matmul(mass_mat(1:nn*ndof,1:nn*ndof), vecc(1:nn*ndof))
709  do i = 1, nn*ndof
710  df(i) = df(i) + ray_k*kb(i)
711  enddo
712  if( ic_type == 511 ) then
713  do i = 1, nn*ndof
714  do j = 1, nn*ndof
715  df(i) = df(i) + damp_mat(i,j)*vecb(j)
716  enddo
717  enddo
718  endif
719 
720  end subroutine calc_damping_mat_and_force_elem
721 
722  subroutine createmat_abort( ic_type, flag, mtype )
723  integer(kind=kint), intent(in) :: ic_type
724  integer(kind=kint), intent(in) :: flag
725  integer(kind=kint), intent(in), optional :: mtype
726 
727  !$omp critical
728  if( flag == 1 ) then
729  write(*,*) '###ERROR### : Element type not supported for static analysis'
730  else if( flag == 2 ) then
731  write(*,*) '###ERROR### : Element type not supported for nonlinear static analysis'
732  else if( flag == 3 ) then
733  write(*,*) '###ERROR### : This element is not supported for this material'
734  else if( flag == 4 ) then
735  write(*,*) '###ERROR### : Element type not supported for dynamic analysis'
736  endif
737  write(*,*) ' ic_type = ', ic_type
738  if( present(mtype) ) write(*,*) ' mtype = ', mtype
739  !$omp end critical
740  call hecmw_abort(hecmw_comm_get_comm())
741  end subroutine createmat_abort
742 
743 end module m_fstr_creatematrix_and_dampingforce
m_dynamic_mass
This module contains subroutines used in 3d eigen analysis for.
Definition: dynamic_mass.f90:6
m_dynamic_mass::get_length
real(kind=kreal) function get_length(ecoord)
Definition: dynamic_mass.f90:246
m_dynamic_mass::mass_beam_33
subroutine mass_beam_33(surf, length, rho, mass)
Definition: dynamic_mass.f90:203
m_dynamic_mass::get_face4
real(kind=kreal) function get_face4(ecoord)
Definition: dynamic_mass.f90:277
m_dynamic_mass::mass_s3
subroutine mass_s3(surf, thick, rho, mass)
Definition: dynamic_mass.f90:164
m_dynamic_mass::get_face3
real(kind=kreal) function get_face3(ecoord)
Definition: dynamic_mass.f90:257
m_dynamic_mass::mass_beam
subroutine mass_beam(surf, length, rho, mass)
Definition: dynamic_mass.f90:190
m_dynamic_mass::mass_c3
subroutine mass_c3(etype, nn, ecoord, gausses, mass, lumped, temperature, is_lumped)
Definition: dynamic_mass.f90:88
m_dynamic_mass::mass_c2
subroutine mass_c2(etype, nn, ecoord, gausses, sec_opt, thick, mass, lumped, temperature, is_lumped)
Definition: dynamic_mass.f90:12
m_dynamic_mass::mass_s4
subroutine mass_s4(surf, thick, rho, mass)
Definition: dynamic_mass.f90:177
m_elemact
This module defined elemact data and function.
Definition: m_element_activation.f90:6
m_elemact::stf_dummy
subroutine stf_dummy(ndof, nn, ecoord, stiff, element)
Definition: m_element_activation.f90:41
m_elemact::kelact_inactive
integer, parameter kelact_inactive
Definition: m_element_activation.f90:30
m_fstr_creatematrix_and_dampingforce
This module assembles the tangent stiffness matrix and, in the implicit dynamic case,...
Definition: fstr_CreateMatrix.f90:10
m_fstr_creatematrix_and_dampingforce::fstr_calc_kinetic_energy
real(kind=kreal) function, public fstr_calc_kinetic_energy(hecMESH, fstrSOLID, fstrDYNAMIC, velocity)
Recompute each element mass matrix and evaluate 0.5*v^T*M*v. Only rows owned by this rank are accumul...
Definition: fstr_CreateMatrix.f90:234
m_fstr_creatematrix_and_dampingforce::fstr_creatematrix_and_dampingforce
subroutine, public fstr_creatematrix_and_dampingforce(hecMESH, hecMAT, fstrSOLID, time, tincr, fstrDYNAMIC, coef)
Assemble the system matrix and, optionally, the dynamic damping force.
Definition: fstr_CreateMatrix.f90:37
m_fstr_creatematrix_and_dampingforce::calc_mass_elem
subroutine calc_mass_elem(ic_type, nn, ndof, ecoord, fstrSOLID, hecMESH, icel, use_lumped_mass, mass_mat, lumped)
Compute the element mass matrix for the given element type.
Definition: fstr_CreateMatrix.f90:559
m_fstr_creatematrix_and_dampingforce::fstr_calc_mass_product
subroutine, public fstr_calc_mass_product(hecMESH, fstrSOLID, fstrDYNAMIC, vec, mass_vec)
Recompute each element mass matrix and evaluate M*vec. Only rows owned by this rank are computed; the...
Definition: fstr_CreateMatrix.f90:307
m_fstr_creatematrix_and_dampingforce::set_diagonal_mass
subroutine set_diagonal_mass(nsize, lumped, mass)
Definition: fstr_CreateMatrix.f90:637
m_fstr_creatematrix_and_dampingforce::fstr_report_lumped_mass_fallback
subroutine, public fstr_report_lumped_mass_fallback(hecMESH, fstrDYNAMIC)
Report the element types that use the lumped mass although the consistent mass is selected....
Definition: fstr_CreateMatrix.f90:379
m_fstr_finiterotationkinematics
Shared finite-rotation nodal kinematics and rotation algebra.
Definition: fstr_FiniteRotationKinematics.f90:7
m_fstr_finiterotationkinematics::fstr_uses_finite_rotation_kinematics
logical function, public fstr_uses_finite_rotation_kinematics(etype, nn, material)
Definition: fstr_FiniteRotationKinematics.f90:42
m_fstr_finiterotationkinematics::shellcomposenodaldisplacement
pure subroutine, public shellcomposenodaldisplacement(ndof, nn, disp_old, disp_inc, disp_new)
Definition: fstr_FiniteRotationKinematics.f90:299
m_fstr_nodalkinematics
Finite-rotation nodal kinematics for NLGEOM.
Definition: fstr_NodalKinematics.f90:14
m_fstr_nodalkinematics::fstr_ensure_finite_rotation_state
subroutine, public fstr_ensure_finite_rotation_state(hecMESH, fstrSOLID, ndof)
Build the per-node reference frames once, by averaging the triads contributed by adjacent local shell...
Definition: fstr_NodalKinematics.f90:34
m_fstr
This module defines common data and basic structures for analysis.
Definition: m_fstr.F90:15
m_fstr::kel361bbar
integer(kind=kint), parameter kel361bbar
Definition: m_fstr.F90:85
m_fstr::kel341sesns
integer(kind=kint), parameter kel341sesns
Definition: m_fstr.F90:82
m_fstr::imsg
integer(kind=kint), parameter imsg
Definition: m_fstr.F90:121
m_fstr::kel361up
integer(kind=kint), parameter kel361up
Definition: m_fstr.F90:88
m_fstr::kel361fi
integer(kind=kint), parameter kel361fi
Definition: m_fstr.F90:84
m_fstr::kmassconsistent
integer(kind=kint), parameter kmassconsistent
Definition: m_fstr.F90:74
m_fstr::get_coordsys
subroutine get_coordsys(cdsys_ID, hecMESH, fstrSOLID, coords, icel)
This subroutine fetch coords defined by local coordinate system.
Definition: m_fstr.F90:1307
m_fstr::kel361ic
integer(kind=kint), parameter kel361ic
Definition: m_fstr.F90:86
m_fstr::kel361fbar
integer(kind=kint), parameter kel361fbar
Definition: m_fstr.F90:87
m_fstr::kmasslumped
integer(kind=kint), parameter kmasslumped
mass matrix type
Definition: m_fstr.F90:73
m_static_lib
This modules just summarizes all modules used in static analysis.
Definition: static_LIB.f90:6
mmechgauss
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
m_fstr::fstr_dynamic
Data for DYNAMIC ANSLYSIS (fstrDYNAMIC)
Definition: m_fstr.F90:525
m_fstr::fstr_solid
Definition: m_fstr.F90:230