!======================================================================! ! ! ! Software Name : FrontISTR Ver. 3.6 ! ! ! ! Module Name : Static Analysis ! ! ! ! Written by K. Sato(Advancesoft), X. YUAN(AdavanceSoft) ! ! ! ! ! ! Contact address : IIS,The University of Tokyo, CISS ! ! ! ! "Structural Analysis for Large Scale Assembly" ! ! ! !======================================================================! !======================================================================! ! !> \brief This module provides function to calcualte to do updates !! !> \author date version !> X.Yuan(Advancesoft) 2009/08/28 original !> X.Yuan 2013/03/18 consider anisotropic materials ! !======================================================================! module m_fstr_Update implicit none contains !=====================================================================* ! UPDATE_C3 !> \brief 変位／応力・ひずみ／内力のアップデート ! !> \author K. Sato(Advancesoft), X. YUAN(AdavanceSoft) !> \version 0.00 !! !> \if AS !> \par サブルーチン構成 !> -# 変位の更新 \f$ u_{n+1}^{(k)} = u_{n+1}^{(k-1)} + \delta u^{(k)} \f$ !> -# ひずみ・応力の更新 \f$ \varepsilon_{n+1}^{(k)} = \varepsilon_{n+1}^{(k-1)} + \delta \varepsilon^{(k)} \f$, \f$ \sigma_{n+1}^{(k)} = \sigma_{n+1}^{(k-1)} + \delta \sigma^{(k)} \f$ !> -# 内力（等価節点力）の計算 \f$ Q_{n+1}^{(k-1)} ( u_{n+1}^{(k-1)} ) \f$ !> \endif ! subroutine fstr_UpdateNewton ( hecMESH, hecMAT, fstrSOLID, tincr,iter, strainEnergy) !=====================================================================* use m_fstr use m_static_lib !#ifdef PARA_CONTACT use m_fstr_para_contact !#endif type (hecmwST_matrix) :: hecMAT !< linear equation, its right side modified here type (hecmwST_local_mesh) :: hecMESH !< mesh information type (fstr_solid) :: fstrSOLID !< we need boundary conditions of curr step real(kind=kreal),intent(in) :: tincr !< time increment integer, intent(in) :: iter !< NR iterations ! Gaku Hashimoto (The University of Tokyo, 2015/03/23) < ! integer(kind=kint) :: nodLOCAL(20) ! real(kind=kreal) :: ecoord(3, 20) integer(kind=kint) :: nodLOCAL(30) real(kind=kreal) :: ecoord(3, 30) ! > Gaku Hashimoto (The University of Tokyo, 2015/03/23) real(kind=kreal) :: thick integer(kind=kint) :: ndof, itype, iS, iE, ic_type, nn, icel, iiS, i, j ! Gaku Hashimoto (The University of Tokyo, 2015/03/23) < ! real(kind=kreal) :: total_disp(6, 20), du(6, 20), ddu(6, 20) ! real(kind=kreal) :: tt(20), tt0(20), ttn(20), qf(20*6), coords(3, 3) real(kind=kreal) :: total_disp(6, 30), du(6, 30), ddu(6, 30) real(kind=kreal) :: tt(30), tt0(30), ttn(30), qf(30*6), coords(3, 3) ! > Gaku Hashimoto (The University of Tokyo, 2015/03/23) integer :: ig0, ig, ik, in, ierror, isect, ihead, cdsys_ID real(kind=kreal), optional :: strainEnergy real(kind=kreal) :: tmp ! Gaku Hashimoto (The University of Tokyo, 2015/03/23) < INTEGER( KIND = kint ) nn_p REAL( KIND = kreal ) :: ddlambda(30) REAL( KIND = kreal ) :: lambda(30) ! > Gaku Hashimoto (The University of Tokyo, 2015/03/23) ndof = hecMAT%NDOF fstrSOLID%QFORCE=0.0d0 tt0 = 0.d0 ttn = 0.d0 tt = 0.d0 ! -------------------------------------------------------------------- ! updated ! 1. stress and strain : ep^(k) = ep^(k-1)+dep^(k) ! sgm^(k) = sgm^(k-1)+dsgm^(k) ! 2. Internal Force : Q^(k-1) ( u^(k-1) ) ! -------------------------------------------------------------------- ! ! ---------------------------------------------------------------------------------- ! calculate the Strain and Stress and Internal Force ( Equivalent Nodal Force ) ! ---------------------------------------------------------------------------------- ! do itype = 1, hecMESH%n_elem_type iS = hecMESH%elem_type_index(itype-1)+1 iE = hecMESH%elem_type_index(itype ) ic_type= hecMESH%elem_type_item(itype) if (hecmw_is_etype_link(ic_type)) cycle nn = hecmw_get_max_node(ic_type) if( nn>20 ) stop "Elemental nodes>20!" ! element loop !$omp parallel default(none), & !$omp& private(icel,iiS,j,nodLOCAL,i,ecoord,ddu,du,total_disp, & !$omp& cdsys_ID,coords,thick,qf,isect,ihead,tmp), & !$omp& shared(iS,iE,hecMESH,nn,fstrSOLID,ndof,hecMAT,ic_type,fstrPR, & !$omp& strainEnergy,iter,tincr), & !$omp& firstprivate(tt0,ttn,tt) !$omp do do icel = iS, iE ! ----- nodal coordinate iiS = hecMESH%elem_node_index(icel-1) do j = 1, nn nodLOCAL(j) = hecMESH%elem_node_item (iiS+j) do i = 1, 3 ecoord(i,j) = hecMESH%node(3*nodLOCAL(j)+i-3) enddo if( fstrSOLID%TEMP_ngrp_tot > 0 .or. fstrSOLID%TEMP_irres > 0 ) then if( isElastoplastic(fstrSOLID%elements(icel)%gausses(1)%pMaterial%mtype) .or. & fstrSOLID%elements(icel)%gausses(1)%pMaterial%mtype == NORTON ) then tt0(j)=fstrSOLID%last_temp( nodLOCAL(j) ) else tt0(j) = 0.d0 if( hecMESH%hecmw_flag_initcon == 1 ) tt0(j) = hecMESH%node_init_val_item(nodLOCAL(j)) endif ttn(j) = fstrSOLID%last_temp( nodLOCAL(j) ) tt(j) = fstrSOLID%temperature( nodLOCAL(j) ) endif enddo ! nodal displacement do j = 1, nn nodLOCAL(j) = hecMESH%elem_node_item (iiS+j) do i = 1, ndof ddu(i,j) = hecMAT%X(ndof*nodLOCAL(j)+i-ndof) du(i,j) = fstrSOLID%dunode(ndof*nodLOCAL(j)+i-ndof) total_disp(i,j) = fstrSOLID%unode(ndof*nodLOCAL(j)+i-ndof) enddo enddo isect = hecMESH%section_ID(icel) cdsys_ID = hecMESH%section%sect_orien_ID(isect) if( cdsys_ID > 0 ) call get_coordsys(cdsys_ID, hecMESH, fstrSOLID, coords) ! ===== calculate the Internal Force if( getSpaceDimension( ic_type ) == 2 ) thick = 1.0d0 if( ic_type == 241 .or. ic_type == 242 .or. ic_type == 231 .or. ic_type == 232 .or. ic_type == 2322 ) then call UPDATE_C2( ic_type,nn,ecoord(1:3,1:nn),fstrSOLID%elements(icel)%gausses(:), & thick,fstrSOLID%elements(icel)%iset, & total_disp(1:2,1:nn), ddu(1:2,1:nn), qf(1:nn*ndof) ) else if( ic_type == 301 ) then isect= hecMESH%section_ID(icel) ihead = hecMESH%section%sect_R_index(isect-1) thick = hecMESH%section%sect_R_item(ihead+1) call UPDATE_C1( ic_type,nn,ecoord(:,1:nn), thick, total_disp(1:3,1:nn), du(1:3,1:nn), & qf(1:nn*ndof),fstrSOLID%elements(icel)%gausses(:) ) else if( ic_type == 361 ) then if( fstrPR%solution_type == kstSTATIC ) then if( fstrSOLID%TEMP_ngrp_tot > 0 .or. fstrSOLID%TEMP_irres > 0 ) then call UpdateST_C3D8Bbar & ( ic_type, nn, ecoord(1, 1:nn), ecoord(2, 1:nn), ecoord(3, 1:nn), tt(1:nn), tt0(1:nn), & ddu(1:3, 1:nn), fstrSOLID%elements(icel)%gausses(:), cdsys_ID, coords ) else call UpdateST_C3D8IC & ( ic_type, nn, ecoord(1, 1:nn), ecoord(2, 1:nn), ecoord(3, 1:nn), tt(1:nn), tt0(1:nn), & ddu(1:3, 1:nn), fstrSOLID%elements(icel)%gausses(:), cdsys_ID, coords ) endif else if( fstrSOLID%TEMP_ngrp_tot > 0 .or. fstrSOLID%TEMP_irres > 0 ) then call UPDATE_C3D8Bbar & ( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_ID, coords, & qf(1:nn*ndof),fstrSOLID%elements(icel)%gausses(:), iter, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn) ) else call Update_C3D8Bbar & ( ic_type,nn,ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_ID, coords, & qf(1:nn*ndof),fstrSOLID%elements(icel)%gausses(:), iter, tincr ) endif endif ! Gaku Hashimoto (The University of Tokyo, 2015/03/23) < ! else if (ic_type == 341 .or. ic_type == 351 .or. & ! ic_type == 342 .or. ic_type == 352 .or. ic_type == 362 ) then else if (ic_type == 341 .or. ic_type == 351 .or. & ic_type == 342 .or. ic_type == 352 .or. ic_type == 362 .or. ic_type == 363) then ! > Gaku Hashimoto (The University of Tokyo, 2015/03/23) if( fstrSOLID%TEMP_ngrp_tot > 0 .or. fstrSOLID%TEMP_irres > 0 ) then call UPDATE_C3 & ( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_ID, coords, & qf(1:nn*ndof), fstrSOLID%elements(icel)%gausses(:), iter, tincr, & tt(1:nn), tt0(1:nn), ttn(1:nn) ) else call UPDATE_C3 & ( ic_type,nn,ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_ID, coords, & qf(1:nn*ndof),fstrSOLID%elements(icel)%gausses(:), iter, tincr ) endif ! Gaku Hashimoto (The University of Tokyo, 2015/03/23) < else if( ic_type == 3610 ) then nn_p = SIZE( fstrSOLID%elements(icel)%p ) DO j = 1, nn_p lambda(j) = -0.5D0*fstrSOLID%elements(icel)%p(j) END DO if( fstrSOLID%TEMP_ngrp_tot > 0 .or. fstrSOLID%TEMP_irres > 0 ) then CALL Update_C3_up & ( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), ddu(1:3,1:nn), & qf(1:nn*ndof), fstrSOLID%elements(icel)%gausses(:), iter, tincr, & nn_p, lambda(1:nn_p), ddlambda(1:nn_p), tt(1:nn), tt0(1:nn) ) else CALL Update_C3_up & ( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), ddu(1:3,1:nn), & qf(1:nn*ndof), fstrSOLID%elements(icel)%gausses(:), iter, tincr, & nn_p, lambda(1:nn_p), ddlambda(1:nn_p) ) endif DO j = 1, nn_p lambda(j) = lambda(j)+ddlambda(j) fstrSOLID%elements(icel)%p(j) = -2.0D0*lambda(j) END DO ! > Gaku Hashimoto (The University of Tokyo, 2015/03/23) ! else if ( ic_type==731) then ! call UPDATE_S3(xx,yy,zz,ee,pp,thick,local_stf) ! call fstr_local_stf_restore_temp(local_stf, nn*ndof, stiffness) ! else if ( ic_type==741) then ! call UPDATE_S4(xx,yy,zz,ee,pp,thick,local_stf) ! call fstr_local_stf_restore_temp(local_stf, nn*ndof, stiffness) else write(*, *) '###ERROR### : Element type not supported for nonlinear static analysis' write(*, *) ' ic_type = ', ic_type call hecmw_abort(hecmw_comm_get_comm()) endif ! ----- calculate the global internal force ( Q(u_{n+1}^{k-1}) ) do j = 1, nn do i = 1, ndof !$omp atomic fstrSOLID%QFORCE(ndof*(nodLOCAL(j)-1)+i) & = fstrSOLID%QFORCE(ndof*(nodLOCAL(j)-1)+i)+qf(ndof*(j-1)+i) enddo enddo ! ----- calculate strain energy if(present(strainEnergy))then do j = 1, nn do i = 1, ndof tmp = 0.5d0*( fstrSOLID%elements(icel)%equiForces(ndof*(j-1)+i) & +qf(ndof*(j-1)+i) )*ddu(i,j) !$omp atomic strainEnergy = strainEnergy+tmp fstrSOLID%elements(icel)%equiForces(ndof*(j-1)+i) = qf(ndof*(j-1)+i) enddo enddo endif enddo ! icel !$omp end do !$omp end parallel enddo ! itype !C !C Update for fstrSOLID%QFORCE !C if( ndof==3 ) then if(paraContactFlag) then call paraContact_update_3_R(hecMESH,fstrSOLID%QFORCE) else call hecmw_update_3_R(hecMESH,fstrSOLID%QFORCE,hecMESH%n_node) endif else if( ndof==2 ) then call hecmw_update_2_R(hecMESH,fstrSOLID%QFORCE,hecMESH%n_node) else if( ndof==6 ) then call hecmw_update_m_R(hecMESH,fstrSOLID%QFORCE,hecMESH%n_node,6) endif end subroutine fstr_UpdateNewton !> Update elastiplastic status subroutine fstr_UpdateState( hecMESH, fstrSOLID, tincr) use m_fstr use m_static_lib use m_ElastoPlastic use mCreep use mViscoElastic type (hecmwST_local_mesh) :: hecMESH !< hecmw mesh type (fstr_solid) :: fstrSOLID !< fstr_solid real(kind=kreal) :: tincr integer(kind=kint) :: itype, iS, iE, ic_type, icel, ngauss, i if( associated( fstrSOLID%temperature ) ) then do i = 1, hecMESH%n_node fstrSOLID%last_temp(i) = fstrSOLID%temperature(i) end do endif do itype = 1, hecMESH%n_elem_type iS = hecMESH%elem_type_index(itype-1) + 1 iE = hecMESH%elem_type_index(itype ) ic_type= hecMESH%elem_type_item(itype) if( ic_type == 301 ) ic_type = 111 if( hecmw_is_etype_link(ic_type) ) cycle ngauss = NumOfQuadPoints( ic_type ) do icel = iS, iE if( isElastoplastic( fstrSOLID%elements(icel)%gausses(1)%pMaterial%mtype ) ) then do i = 1, ngauss call updateEPState( fstrSOLID%elements(icel)%gausses(i) ) enddo elseif( fstrSOLID%elements(icel)%gausses(1)%pMaterial%mtype == NORTON ) then if( tincr>0.d0 ) then do i = 1, ngauss fstrSOLID%elements(icel)%gausses(i)%ttime = & fstrSOLID%elements(icel)%gausses(i)%ttime+tincr call updateViscoState( fstrSOLID%elements(icel)%gausses(i) ) enddo endif elseif( isViscoelastic( fstrSOLID%elements(icel)%gausses(1)%pMaterial%mtype ) ) then if( tincr > 0.d0 ) then do i = 1, ngauss fstrSOLID%elements(icel)%gausses(i)%ttime = & fstrSOLID%elements(icel)%gausses(i)%ttime+tincr call updateViscoElasticState( fstrSOLID%elements(icel)%gausses(i) ) enddo endif endif do i = 1, ngauss fstrSOLID%elements(icel)%gausses(i)%strain_bak = fstrSOLID%elements(icel)%gausses(i)%strain fstrSOLID%elements(icel)%gausses(i)%stress_bak = fstrSOLID%elements(icel)%gausses(i)%stress enddo enddo enddo end subroutine fstr_UpdateState !> Update at linear static analysis ---------------------------------------------------------------------- subroutine fstr_Update3D( hecMESH, fstrSOLID ) use m_fstr use m_static_lib type (hecmwST_local_mesh) :: hecMESH type (fstr_solid) :: fstrSOLID !C** local variables integer(kind=kint) :: itype, icel, iS, iE, jS, i, j, ic_type, ig, nn, isect, ihead, iflag, mixflag ! Gaku Hashimoto (The University of Tokyo, 2015/03/23) < ! integer(kind=kint) :: nodLOCAL(20), flag_dof ! real(kind=kreal) :: xx(20), yy(20), zz(20), tt(20), tt0(20), edisp(60), force(60) ! real(kind=kreal) :: ecoord(3, 20), stiff(60, 60) integer(kind=kint) :: nodLOCAL(30), flag_dof real(kind=kreal) :: xx(30), yy(30), zz(30), tt(30), tt0(30), edisp(3*30), force(3*30) real(kind=kreal) :: ecoord(3, 30), stiff(3*30, 3*30) ! > Gaku Hashimoto (The University of Tokyo, 2015/03/23) real(kind=kreal) :: thick, coords(3, 3) real(kind=kreal), allocatable :: temp(:) integer(kind=kint), allocatable :: id_spc(:) real(kind=kreal) :: a integer(kind=kint) :: cdsys_ID tt = 0.d0 tt0 = 0.d0 !C !C set temperature !C allocate( temp(hecMESH%n_node) ) temp = 0.0d0 do i = 1, fstrSOLID%TEMP_ngrp_tot ig = fstrSOLID%TEMP_ngrp_ID(i) iS = hecMESH%node_group%grp_index(ig-1) + 1 iE = hecMESH%node_group%grp_index(ig ) do j = iS, iE temp( hecMESH%node_group%grp_item(j) ) = fstrSOLID%TEMP_ngrp_val(i) enddo enddo !C !C set boundary force !C fstrSOLID%QFORCE = 0.0d0 allocate ( id_spc(hecMESH%n_node) ) id_spc = 0 do i = 1, fstrSOLID%BOUNDARY_ngrp_tot ig = fstrSOLID%BOUNDARY_ngrp_ID(i) iS = hecMESH%node_group%grp_index(ig-1) + 1 iE = hecMESH%node_group%grp_index(ig ) do j = iS, iE id_spc( hecMESH%node_group%grp_item(j) ) = 1 enddo enddo !C +-------------------------------+ !C | according to ELEMENT TYPE | !C +-------------------------------+ do itype = 1, hecMESH%n_elem_type iS = hecMESH%elem_type_index(itype-1) + 1 iE = hecMESH%elem_type_index(itype ) ic_type = hecMESH%elem_type_item(itype) if( ic_type == fe_tet10nc ) ic_type = fe_tet10n if( .not. (hecmw_is_etype_solid(ic_type) .or. ic_type == 781 .or. ic_type == 761 .or. ic_type == 641 )) cycle nn = hecmw_get_max_node( ic_type ) !C element loop !$omp parallel default(none), & !$omp& private(icel,js,j,nodLOCAL,xx,yy,zz,ecoord,edisp,isect,cdsys_ID,ihead,thick,mixflag,stiff,iflag,coords), & !$omp& shared(iS,iE,hecMESH,temp,ref_temp,fstrSOLID,id_spc,force), & !$omp& firstprivate(nn,ic_type,tt,tt0) !$omp do do icel = iS, iE jS = hecMESH%elem_node_index(icel-1) do j = 1, nn nodLOCAL(j) = hecMESH%elem_node_item(jS+j) xx(j) = hecMESH%node(3*nodLOCAL(j)-2) yy(j) = hecMESH%node(3*nodLOCAL(j)-1) zz(j) = hecMESH%node(3*nodLOCAL(j) ) if( fstrSOLID%TEMP_ngrp_tot > 0 .or. fstrSOLID%TEMP_irres > 0 ) then if( isElastoplastic(fstrSOLID%elements(icel)%gausses(1)%pMaterial%mtype) .or. & fstrSOLID%elements(icel)%gausses(1)%pMaterial%mtype==NORTON ) then tt0(j)=fstrSOLID%last_temp( nodLOCAL(j) ) else tt0(j) = 0.d0 if( hecMESH%hecmw_flag_initcon==1 ) tt0(j) = hecMESH%node_init_val_item(nodLOCAL(j)) endif tt(j) = fstrSOLID%temperature( nodLOCAL(j) ) endif ecoord(1,j) = hecMESH%node(3*nodLOCAL(j)-2) ecoord(2,j) = hecMESH%node(3*nodLOCAL(j)-1) ecoord(3,j) = hecMESH%node(3*nodLOCAL(j) ) edisp(3*j-2) = fstrSOLID%unode(3*nodLOCAL(j)-2) edisp(3*j-1) = fstrSOLID%unode(3*nodLOCAL(j)-1) edisp(3*j ) = fstrSOLID%unode(3*nodLOCAL(j) ) enddo isect = hecMESH%section_ID(icel) cdsys_ID = hecMESH%section%sect_orien_ID(isect) if( cdsys_ID > 0 ) call get_coordsys(cdsys_ID, hecMESH, fstrSOLID, coords) !--- calculate stress and strain of gauss point if( ic_type == 781 ) then !for shell-solid mixed analysis isect= hecMESH%section_ID(icel) ihead = hecMESH%section%sect_R_index(isect-1) thick = hecMESH%section%sect_R_item(ihead+1) nn = 4; ic_type = 741; mixflag = 1 call STF_Shell_MITC(ic_type, nn, 6, ecoord(1:3, 1:8), fstrSOLID%elements(icel)%gausses, & stiff, thick, mixflag) ic_type = 781; nn = 8 CYCLE else if( ic_type == 761 ) then !for shell-solid mixed analysis isect= hecMESH%section_ID(icel) ihead = hecMESH%section%sect_R_index(isect-1) thick = hecMESH%section%sect_R_item(ihead+1) nn = 3; ic_type = 731; mixflag = 2 call STF_Shell_MITC(ic_type, nn, 6, ecoord(1:3, 1:8), fstrSOLID%elements(icel)%gausses, & stiff, thick, mixflag) ic_type = 761; nn = 6 CYCLE else if( ic_type == 641 ) THEN isect = hecMESH%section_ID(icel) ihead = hecMESH%section%sect_R_index(isect-1) CALL STF_Beam_641( ic_type, nn, ecoord, fstrSOLID%elements(icel)%gausses(:), & hecMESH%section%sect_R_item(ihead+1:), stiff ) iflag = 0 DO j = 1, nn IF( id_spc( nodLOCAL(j) ) == 1 ) iflag = 1 END DO if( iflag == 1 ) THEN IF( ic_type == 641 ) THEN isect = hecMESH%section_ID(icel) ihead = hecMESH%section%sect_R_index(isect-1) CALL STF_Beam_641( ic_type, nn, ecoord, fstrSOLID%elements(icel)%gausses(:), & hecMESH%section%sect_R_item(ihead+1:), stiff ) END IF force(1:nn*3) = MATMUL( stiff(1:nn*3,1:nn*3), edisp(1:nn*3) ) DO j = 1, nn IF( id_spc( nodLOCAL(j) ) == 1 ) THEN !$omp atomic fstrSOLID%QFORCE( 3*nodLOCAL(j)-2 ) = & fstrSOLID%QFORCE( 3*nodLOCAL(j)-2 )+force(3*j-2) !$omp atomic fstrSOLID%QFORCE( 3*nodLOCAL(j)-1 ) = & fstrSOLID%QFORCE( 3*nodLOCAL(j)-1 )+force(3*j-1) !$omp atomic fstrSOLID%QFORCE( 3*nodLOCAL(j) ) = & fstrSOLID%QFORCE( 3*nodLOCAL(j) )+force(3*j ) END IF END DO ENDIF CYCLE else if( ic_type == 361 ) then if( fstrSOLID%TEMP_ngrp_tot > 0 .or. fstrSOLID%TEMP_irres > 0 ) then call UpdateST_C3D8Bbar & ( ic_type, nn, xx, yy, zz, tt, tt0, edisp, & fstrSOLID%elements(icel)%gausses, cdsys_ID, coords ) else call UpdateST_C3D8IC & ( ic_type, nn, xx, yy, zz, tt, tt0, edisp, & fstrSOLID%elements(icel)%gausses, cdsys_ID, coords ) endif else if( ic_type == 301 ) then isect = hecMESH%section_ID(icel) ihead = hecMESH%section%sect_R_index(isect-1) thick = hecMESH%section%sect_R_item(ihead+1) call UpdateST_C1( ic_type, nn, xx, yy, zz, tt, tt0, thick, edisp, fstrSOLID%elements(icel)%gausses ) else call UpdateST_C3 & ( ic_type, nn, xx, yy, zz, tt, tt0, edisp, & fstrSOLID%elements(icel)%gausses, cdsys_ID, coords ) endif !--- calculate reaction force iflag = 0 do j = 1, nn if( id_spc( nodLOCAL(j) ) == 1 ) iflag = 1 enddo if( iflag == 1 ) then if( ic_type == 361 ) then if( fstrSOLID%TEMP_ngrp_tot > 0 .or. fstrSOLID%TEMP_irres > 0 ) then call STF_C3D8Bbar & ( ic_type, nn, ecoord, fstrSOLID%elements(icel)%gausses, stiff, cdsys_ID, coords, 1.0D0 ) else call STF_C3D8IC & ( ic_type, nn, ecoord, fstrSOLID%elements(icel)%gausses, stiff, cdsys_ID, coords ) endif else if( ic_type == 301 ) then isect = hecMESH%section_ID(icel) ihead = hecMESH%section%sect_R_index(isect-1) thick = hecMESH%section%sect_R_item(ihead+1) call STF_C1( ic_type, nn, ecoord, thick, fstrSOLID%elements(icel)%gausses, stiff ) else call STF_C3 & (ic_type, nn, ecoord, fstrSOLID%elements(icel)%gausses, stiff, cdsys_ID, coords, 1.0D0) endif force(1:nn*3) = matmul( stiff(1:nn*3,1:nn*3), edisp(1:nn*3) ) do j = 1, nn if( id_spc( nodLOCAL(j) ) == 1 ) then !$omp atomic fstrSOLID%QFORCE(3*nodLOCAL(j)-2) = fstrSOLID%QFORCE(3*nodLOCAL(j)-2) + force(3*j-2) !$omp atomic fstrSOLID%QFORCE(3*nodLOCAL(j)-1) = fstrSOLID%QFORCE(3*nodLOCAL(j)-1) + force(3*j-1) !$omp atomic fstrSOLID%QFORCE(3*nodLOCAL(j) ) = fstrSOLID%QFORCE(3*nodLOCAL(j) ) + force(3*j ) endif enddo endif enddo !$omp end do !$omp end parallel enddo call hecmw_update_3_R( hecMESH, fstrSOLID%QFORCE, hecMESH%n_node ) deallocate( temp ) deallocate( id_spc ) end subroutine fstr_Update3D !> Update at linear static analysis subroutine fstr_Update2D( hecMESH, fstrSOLID ) use m_fstr use m_static_lib type (hecmwST_local_mesh) :: hecMESH type (fstr_solid) :: fstrSOLID !C** local variables integer(kind=kint) :: itype, icel, iS, iE, jS, i, j, ic_type, ig, nn, iflag integer(kind=kint) :: nodLOCAL(8) real(kind=kreal) :: xx(8), yy(8), tt(8), tt0(8), edisp(16), force(16) real(kind=kreal) :: ecoord(2,8), stiff(16,16) real(kind=kreal), allocatable :: temp(:) integer(kind=kint), allocatable :: id_spc(:) !C !C set temperature !C allocate( temp(hecMESH%n_node) ) temp = 0.0d0 do i = 1, fstrSOLID%TEMP_ngrp_tot ig = fstrSOLID%TEMP_ngrp_ID(i) iS = hecMESH%node_group%grp_index(ig-1)+1 iE = hecMESH%node_group%grp_index(ig ) do j = iS, iE temp( hecMESH%node_group%grp_item(j) ) = fstrSOLID%TEMP_ngrp_val(i) enddo enddo !C !C set boundary force !C fstrSOLID%QFORCE = 0.0d0 allocate ( id_spc(hecMESH%n_node) ) id_spc = 0 do i = 1, fstrSOLID%BOUNDARY_ngrp_tot ig = fstrSOLID%BOUNDARY_ngrp_ID(i) iS = hecMESH%node_group%grp_index(ig-1) + 1 iE = hecMESH%node_group%grp_index(ig ) do j = iS, iE id_spc( hecMESH%node_group%grp_item(j) ) = 1 enddo enddo !C +-------------------------------+ !C | according to ELEMENT TYPE | !C +-------------------------------+ do itype = 1, hecMESH%n_elem_type iS = hecMESH%elem_type_index(itype-1) + 1 iE = hecMESH%elem_type_index(itype ) ic_type = hecMESH%elem_type_item(itype) if( .not. hecmw_is_etype_surface(ic_type) ) cycle nn = hecmw_get_max_node( ic_type ) !C element loop do icel = iS, iE jS = hecMESH%elem_node_index(icel-1) do j = 1, nn nodLOCAL(j) = hecMESH%elem_node_item(jS+j) xx(j) = hecMESH%node(3*nodLOCAL(j)-2) yy(j) = hecMESH%node(3*nodLOCAL(j)-1) tt(j) = temp( nodLOCAL(j) ) tt0(j)= ref_temp ecoord(1,j) = hecMESH%node(3*nodLOCAL(j)-2) ecoord(2,j) = hecMESH%node(3*nodLOCAL(j)-1) edisp(2*j-1) = fstrSOLID%unode(2*nodLOCAL(j)-1) edisp(2*j ) = fstrSOLID%unode(2*nodLOCAL(j) ) enddo !--- calculate stress and strain of gauss points call UpdateST_C2( ic_type, nn, xx, yy, tt, tt0, 1.0d0, fstrSOLID%elements(icel)%iset, & edisp, fstrSOLID%elements(icel)%gausses ) !--- calculate reaction force iflag = 0 do j = 1, nn if( id_spc( nodLOCAL(j) ) == 1 ) iflag = 1 enddo if( iflag == 1 ) then call STF_C2( ic_type, nn, ecoord, fstrSOLID%elements(icel)%gausses, 1.0d0, & stiff, fstrSOLID%elements(icel)%iset ) force(1:nn*2) = matmul( stiff(1:nn*2,1:nn*2), edisp(1:nn*2) ) do j = 1, nn if( id_spc( nodLOCAL(j) ) == 1 ) then fstrSOLID%QFORCE(2*nodLOCAL(j)-1) = fstrSOLID%QFORCE(2*nodLOCAL(j)-1) + force(2*j-1) fstrSOLID%QFORCE(2*nodLOCAL(j) ) = fstrSOLID%QFORCE(2*nodLOCAL(j) ) + force(2*j ) endif enddo endif enddo enddo call hecmw_update_2_R( hecMESH, fstrSOLID%QFORCE, hecMESH%n_node ) deallocate( temp ) deallocate( id_spc ) end subroutine fstr_Update2D !> Update at linear static analysis subroutine fstr_Update6D( hecMESH, fstrSOLID ) use m_fstr use m_static_lib type (hecmwST_local_mesh) :: hecMESH type (fstr_solid) :: fstrSOLID !C** local variables integer(kind=kint) :: itype, icel, iS, iE, jS, i, j, ic_type, ig, nn, isect, ihead, iflag, mixflag integer(kind=kint) :: nodLOCAL(9) real(kind=kreal) :: ecoord(3,9), edisp(54), force(54), stiff(54,54) real(kind=kreal) :: thick, logr integer(kind=kint), allocatable :: id_spc(:) !C !C set boundary force !C fstrSOLID%QFORCE = 0.0d0 allocate ( id_spc(hecMESH%n_node) ) id_spc = 0 do i = 1, fstrSOLID%BOUNDARY_ngrp_tot ig = fstrSOLID%BOUNDARY_ngrp_ID(i) iS = hecMESH%node_group%grp_index(ig-1) + 1 iE = hecMESH%node_group%grp_index(ig ) do j = iS, iE id_spc( hecMESH%node_group%grp_item(j) ) = 1 enddo enddo !C +-------------------------------+ !C | according to ELEMENT TYPE | !C +-------------------------------+ do itype = 1, hecMESH%n_elem_type iS = hecMESH%elem_type_index(itype-1) + 1 iE = hecMESH%elem_type_index(itype ) ic_type = hecMESH%elem_type_item(itype) if( .not. hecmw_is_etype_shell(ic_type) ) cycle nn = hecmw_get_max_node( ic_type ) !C element loop do icel = iS, iE jS = hecMESH%elem_node_index(icel-1) do j = 1, nn nodLOCAL(j) = hecMESH%elem_node_item(jS+j) ecoord(1,j) = hecMESH%node(3*nodLOCAL(j)-2) ecoord(2,j) = hecMESH%node(3*nodLOCAL(j)-1) ecoord(3,j) = hecMESH%node(3*nodLOCAL(j) ) edisp(6*j-5) = fstrSOLID%unode(6*nodLOCAL(j)-5) edisp(6*j-4) = fstrSOLID%unode(6*nodLOCAL(j)-4) edisp(6*j-3) = fstrSOLID%unode(6*nodLOCAL(j)-3) edisp(6*j-2) = fstrSOLID%unode(6*nodLOCAL(j)-2) edisp(6*j-1) = fstrSOLID%unode(6*nodLOCAL(j)-1) edisp(6*j ) = fstrSOLID%unode(6*nodLOCAL(j) ) enddo isect = hecMESH%section_ID(icel) ihead = hecMESH%section%sect_R_index(isect-1) thick = hecMESH%section%sect_R_item(ihead+1) !--- calculate reaction force if( ic_type == 731 .or. ic_type == 741 .or. ic_type == 743 ) then iflag = 0 do j = 1, nn if( id_spc( nodLOCAL(j) ) == 1 ) iflag = 1 enddo if( iflag == 1 ) then mixflag = 0 call STF_Shell_MITC( ic_type, nn, 6, ecoord(1:3,1:nn), fstrSOLID%elements(icel)%gausses, & stiff(1:nn*6,1:nn*6), thick, mixflag ) force(1:nn*6) = matmul( stiff(1:nn*6,1:nn*6), edisp(1:nn*6) ) do j = 1, nn if( id_spc( nodLOCAL(j) ) == 1 ) then fstrSOLID%QFORCE(6*nodLOCAL(j)-5) = fstrSOLID%QFORCE(6*nodLOCAL(j)-5) + force(6*j-5) fstrSOLID%QFORCE(6*nodLOCAL(j)-4) = fstrSOLID%QFORCE(6*nodLOCAL(j)-4) + force(6*j-4) fstrSOLID%QFORCE(6*nodLOCAL(j)-3) = fstrSOLID%QFORCE(6*nodLOCAL(j)-3) + force(6*j-3) fstrSOLID%QFORCE(6*nodLOCAL(j)-2) = fstrSOLID%QFORCE(6*nodLOCAL(j)-2) + force(6*j-2) fstrSOLID%QFORCE(6*nodLOCAL(j)-1) = fstrSOLID%QFORCE(6*nodLOCAL(j)-1) + force(6*j-1) fstrSOLID%QFORCE(6*nodLOCAL(j) ) = fstrSOLID%QFORCE(6*nodLOCAL(j) ) + force(6*j ) endif enddo endif endif enddo enddo call hecmw_update_m_R( hecMESH, fstrSOLID%QFORCE, hecMESH%n_node, 6 ) deallocate( id_spc ) end subroutine fstr_Update6D end module m_fstr_Update