!======================================================================! ! ! ! Software Name : FrontISTR Ver. 3.6 ! ! ! ! Module Name : Static Analysis ! ! ! ! Written by X. YUAN(AdavanceSoft), K. Sato(Advancesoft) ! ! ! ! ! ! Contact address : IIS,The University of Tokyo, CISS ! ! ! ! "Structural Analysis for Large Scale Assembly" ! ! ! !======================================================================! !======================================================================! ! !> \brief This module provides functions to take into acount external load !! !> \author date version !> X.Yuan(Advancesoft) 2009/08/26 original !> X.Yuan 2013/03/18 consider anisotropic expansion ! !======================================================================! !======================================================================! module m_fstr_ass_load implicit none contains ! !======================================================================! !> This subroutine assmble following external force into fstrSOLID%GL and hecMAT%B afterwards !> -# concentrated nodal force !> -# surface pressure !> -# volume force !> -# thermal force subroutine fstr_ass_load(cstep, hecMESH, hecMAT, fstrSOLID, fstrPARAM) !======================================================================! use m_fstr use m_static_lib use m_fstr_precheck use mMechGauss use mReadTemp use mULoad !#ifdef PARA_CONTACT use m_fstr_para_contact !#endif integer, intent(in) :: cstep !< current step type (hecmwST_matrix),intent(inout) :: hecMAT !< hecmw matrix type (hecmwST_local_mesh),intent(in) :: hecMESH !< hecmw mesh type (fstr_solid),intent(inout) :: fstrSOLID !< fstr_solid type (fstr_param),intent(inout) :: fstrPARAM !< analysis control parameters ! Gaku Hashimoto (The University of Tokyo, 2015/03/23) < ! real(kind=kreal) :: xx(20), yy(20), zz(20) real(kind=kreal) :: xx(30), yy(30), zz(30) ! > Gaku Hashimoto (The University of Tokyo, 2015/03/23) real(kind=kreal) :: params(0:6) ! Gaku Hashimoto (The University of Tokyo, 2015/03/23) < ! real(kind=kreal) :: vect(60) ! integer(kind=kint) :: iwk(60) ! integer(kind=kint) :: nodLocal(20) ! real(kind=kreal) :: tt(20), tt0(20), coords(3,3), factor real(kind=kreal) :: vect(30*3) integer(kind=kint) :: iwk(30*3) integer(kind=kint) :: nodLocal(30) real(kind=kreal) :: tt(30), tt0(30), coords(3,3), factor ! > Gaku Hashimoto (The University of Tokyo, 2015/03/23) integer(kind=kint) :: ndof, ig0, ig, ityp, ltype, iS0, iE0, ik, in, i, j integer(kind=kint) :: icel, ic_type, nn, iS, isect, id, iset, nsize integer(kind=kint) :: itype, iE, ierror, grpid, cdsys_ID real(kind=kreal) :: fval, rho, thick, pa1 logical :: fg_surf integer(kind=kint) :: tstep type( tMaterial ), pointer :: material !< material information integer(kind=kint) :: ihead real(kind=kreal) :: a ndof = hecMAT%NDOF ! ------------------------------------------------------------------- ! CLOAD ! ------------------------------------------------------------------- fstrSOLID%GL(:) = 0.0d0 do ig0 = 1, fstrSOLID%CLOAD_ngrp_tot grpid = fstrSOLID%CLOAD_ngrp_GRPID(ig0) if( .not. fstr_isLoadActive( fstrSOLID, grpid, cstep ) ) cycle factor = fstrSOLID%factor(2) if( cstep > 1 ) then if( fstr_isLoadActive( fstrSOLID, grpid, cstep-1 ) ) factor = 1.0d0 endif ig = fstrSOLID%CLOAD_ngrp_ID(ig0) ityp = fstrSOLID%CLOAD_ngrp_DOF(ig0) fval = fstrSOLID%CLOAD_ngrp_val(ig0) iS0 = hecMESH%node_group%grp_index(ig-1) + 1 iE0 = hecMESH%node_group%grp_index(ig ) do ik = iS0, iE0 in = hecMESH%node_group%grp_item(ik) fstrSOLID%GL(ndof*(in-1)+ityp) = fstrSOLID%GL(ndof*(in-1)+ityp)+factor*fval enddo enddo ! ! ! ------------------------------------------------------------------- ! DLOAD ! ------------------------------------------------------------------- do ig0 = 1, fstrSOLID%DLOAD_ngrp_tot grpid = fstrSOLID%DLOAD_ngrp_GRPID(ig0) if( .not. fstr_isLoadActive( fstrSOLID, grpid, cstep ) ) cycle factor = fstrSOLID%factor(2) if( cstep > 1 ) then if( fstr_isLoadActive( fstrSOLID, grpid, cstep-1 ) ) factor = 1.0d0 endif ig = fstrSOLID%DLOAD_ngrp_ID(ig0) ltype = fstrSOLID%DLOAD_ngrp_LID(ig0) do i = 0, 6 params(i)= fstrSOLID%DLOAD_ngrp_params(i,ig0) enddo ! ----- START & END fg_surf = (ltype == 100) if( fg_surf ) then ! surface group iS0 = hecMESH%surf_group%grp_index(ig-1) + 1 iE0 = hecMESH%surf_group%grp_index(ig ) else ! element group iS0 = hecMESH%elem_group%grp_index(ig-1) + 1 iE0 = hecMESH%elem_group%grp_index(ig ) endif do ik = iS0, iE0 if( fg_surf ) then ! surface group ltype = hecMESH%surf_group%grp_item(2*ik)*10 icel = hecMESH%surf_group%grp_item(2*ik-1) ic_type = hecMESH%elem_type(icel) else ! element group icel = hecMESH%elem_group%grp_item(ik) ic_type = hecMESH%elem_type(icel) endif if( hecmw_is_etype_link(ic_type) ) cycle ! if( ic_type==3422 ) ic_type=342 nn = hecmw_get_max_node(ic_type) ! ----- node ID iS = hecMESH%elem_node_index(icel-1) if( fstrSOLID%DLOAD_follow == 0 ) then do j = 1, nn nodLOCAL(j) = hecMESH%elem_node_item (iS+j) ! ----- nodal coordinate xx(j) = hecMESH%node( 3*nodLOCAL(j)-2 ) yy(j) = hecMESH%node( 3*nodLOCAL(j)-1 ) zz(j) = hecMESH%node( 3*nodLOCAL(j) ) ! ----- create iwk array *** do i = 1, ndof iwk( ndof*(j-1)+i ) = ndof*( nodLOCAL(j)-1 )+i enddo enddo else do j = 1, nn nodLOCAL(j) = hecMESH%elem_node_item (iS+j) ! ----- nodal coordinate if (ndof==2) then xx(j) = hecMESH%node( 3*nodLOCAL(j)-2 )+fstrSOLID%unode( 2*nodLOCAL(j)-1 )+fstrSOLID%dunode( 2*nodLOCAL(j)-1 ) yy(j) = hecMESH%node( 3*nodLOCAL(j)-1 )+fstrSOLID%unode( 2*nodLOCAL(j) )+fstrSOLID%dunode( 2*nodLOCAL(j) ) else if (ndof==3) then xx(j) = hecMESH%node( 3*nodLOCAL(j)-2 )+fstrSOLID%unode( 3*nodLOCAL(j)-2 )+fstrSOLID%dunode( 3*nodLOCAL(j)-2 ) yy(j) = hecMESH%node( 3*nodLOCAL(j)-1 )+fstrSOLID%unode( 3*nodLOCAL(j)-1 )+fstrSOLID%dunode( 3*nodLOCAL(j)-1 ) zz(j) = hecMESH%node( 3*nodLOCAL(j) )+fstrSOLID%unode( 3*nodLOCAL(j) )+fstrSOLID%dunode( 3*nodLOCAL(j) ) else if (ndof==6) then xx(j) = hecMESH%node( 3*nodLOCAL(j)-2 )+fstrSOLID%unode( 6*nodLOCAL(j)-5 )+fstrSOLID%dunode( 6*nodLOCAL(j)-5 ) yy(j) = hecMESH%node( 3*nodLOCAL(j)-1 )+fstrSOLID%unode( 6*nodLOCAL(j)-4 )+fstrSOLID%dunode( 6*nodLOCAL(j)-4 ) zz(j) = hecMESH%node( 3*nodLOCAL(j) )+fstrSOLID%unode( 6*nodLOCAL(j)-3 )+fstrSOLID%dunode( 6*nodLOCAL(j)-3 ) endif ! ----- create iwk array *** do i = 1, ndof iwk( ndof*(j-1)+i ) = ndof*( nodLOCAL(j)-1 )+i enddo enddo end if ! ----- section ID isect = hecMESH%section_ID(icel) ! ----- Get Properties material => fstrSOLID%elements(icel)%gausses(1)%pMaterial rho = material%variables(M_DENSITY) call fstr_get_thickness(hecMESH,isect,thick) ! ----- Section Data if( ndof == 2 ) then id=hecMESH%section%sect_opt(isect) if( id == 0 ) then iset = 1 else if( id == 1 ) then iset = 0 else if( id == 2) then iset = 2 endif pa1=1.d0 endif ! ----- Create local stiffness if (ic_type==301)then ihead = hecMESH%section%sect_R_index(isect-1) call DL_C1(ic_type,nn,xx(1:nn),yy(1:nn),zz(1:nn),rho,thick,ltype,params,vect(1:nn*ndof),nsize) elseif( ic_type == 241 .or. ic_type == 242 .or. ic_type == 231 .or. ic_type == 232 .or. ic_type == 2322 ) then call DL_C2(ic_type,nn,xx(1:nn),yy(1:nn),rho,pa1,ltype,params,vect(1:nn*ndof),nsize,iset) ! Gaku Hashimoto (The University of Tokyo, 2015/03/23) < ! else if ( ic_type == 341 .or. ic_type == 351 .or. ic_type == 361 .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 == 361 .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) call DL_C3(ic_type,nn,xx(1:nn),yy(1:nn),zz(1:nn),rho,ltype,params,vect(1:nn*ndof),nsize) else if ( ic_type == 641 ) then ihead = hecMESH%section%sect_R_index(isect-1) call DL_Beam_641(ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), rho, ltype, params, & hecMESH%section%sect_R_item(ihead+1:), vect(1:nn*ndof), nsize) else if( ( ic_type == 741 ) .or. ( ic_type == 743 ) .or. ( ic_type == 731 ) ) then call DL_Shell(ic_type, nn, ndof, xx, yy, zz, rho, thick, ltype, params, vect, nsize, fstrSOLID%elements(icel)%gausses) else if( ( ic_type==761 ) .or. ( ic_type==781 ) ) then call DL_Shell_33(ic_type, nn, ndof, xx, yy, zz, rho, thick, ltype, params, vect, nsize, & fstrSOLID%elements(icel)%gausses) else nsize = 0 write(*,*)"### WARNING: DLOAD",ic_type endif ! ----- Add vector do j=1,nsize fstrSOLID%GL( iwk(j) )=fstrSOLID%GL( iwk(j) )+factor*vect(j) enddo enddo enddo ! -----Uload call uloading( cstep, factor, fstrSOLID%GL ) !C !C Update for fstrSOLID%GL !C if( hecMESH%n_dof == 3 ) then if(paraContactFlag) then call paraContact_update_3_R (hecMESH,fstrSOLID%GL) else call hecmw_update_3_R (hecMESH,fstrSOLID%GL,hecMESH%n_node) endif else if( hecMESH%n_dof == 2 ) then call hecmw_update_2_R (hecMESH,fstrSOLID%GL,hecMESH%n_node) endif call hecmw_mat_clear_b( hecMAT ) do i=1, hecMESH%n_node* hecMESH%n_dof hecMAT%B(i)=fstrSOLID%GL(i)-fstrSOLID%QFORCE(i) enddo ! ! ! ------------------------------------------------------------------- ! TLOAD : THERMAL LOAD USING TEMPERATURE ! ------------------------------------------------------------------- !C !C Set Temperature !C if( fstrSOLID%TEMP_ngrp_tot > 0 .or. fstrSOLID%TEMP_irres > 0 ) then do ig0 = 1, fstrSOLID%TEMP_ngrp_tot grpid = fstrSOLID%TEMP_ngrp_GRPID(ig0) if( .not. fstr_isLoadActive( fstrSOLID, grpid, cstep ) ) cycle factor = fstrSOLID%factor(2) if( cstep > 1 ) then if( fstr_isLoadActive( fstrSOLID, grpid, cstep-1 ) ) factor = 1.d0 endif ig = fstrSOLID%TEMP_ngrp_ID(ig0) fval =fstrSOLID%TEMP_ngrp_val(ig0) iS0 = hecMESH%node_group%grp_index(ig-1)+1 iE0 = hecMESH%node_group%grp_index(ig ) do ik = iS0, iE0 in = hecMESH%node_group%grp_item(ik) pa1 = fstrSOLID%temp_bak( in ) fstrSOLID%temperature( in ) = pa1+(fval-pa1)*factor enddo enddo if( fstrSOLID%TEMP_irres > 0 ) then call read_temperature_result(hecMESH, fstrSOLID%TEMP_irres, fstrSOLID%TEMP_tstep, fstrSOLID%temperature) endif ! ----- element TYPE loop. 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 ! ----- Set number of nodes nn = hecmw_get_max_node(ic_type) ! ----- element loop do icel = iS, iE ! ----- node ID iS= hecMESH%elem_node_index(icel-1) do j=1,nn nodLOCAL(j)=hecMESH%elem_node_item(iS+j) ! ----- nodal coordinate if (ndof==2) then xx(j)=hecMESH%node(3*nodLOCAL(j)-2)+fstrSOLID%unode(ndof*nodLOCAL(j)-1) yy(j)=hecMESH%node(3*nodLOCAL(j)-1)+fstrSOLID%unode(ndof*nodLOCAL(j) ) else if (ndof==3) then xx(j)=hecMESH%node(3*nodLOCAL(j)-2)+fstrSOLID%unode(ndof*nodLOCAL(j)-2) yy(j)=hecMESH%node(3*nodLOCAL(j)-1)+fstrSOLID%unode(ndof*nodLOCAL(j)-1) zz(j)=hecMESH%node(3*nodLOCAL(j) )+fstrSOLID%unode(ndof*nodLOCAL(j)) endif tt0(j)=fstrSOLID%last_temp( nodLOCAL(j) ) tt(j) = fstrSOLID%temperature( nodLOCAL(j) ) ! ----- create iwk array *** do i=1,ndof iwk(ndof*(j-1)+i)=ndof*(nodLOCAL(j)-1)+i enddo enddo ! ----- section Data isect= hecMESH%section_ID(icel) cdsys_ID = hecMESH%section%sect_orien_ID(isect) call get_coordsys(cdsys_ID, hecMESH, fstrSOLID, coords) if( ndof == 2 ) then id=hecMESH%section%sect_opt(isect) if( id==0 ) then iset = 1 else if( id == 1 ) then iset = 0 else if( id == 2 ) then iset = 2 endif pa1 = 1.0d0 endif IF( ic_type == 641 ) THEN isect= hecMESH%section_ID(icel) ihead = hecMESH%section%sect_R_index(isect-1) CALL TLOAD_Beam_641( ic_type, nn, ndof, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), & fstrSOLID%elements(icel)%gausses, hecMESH%section%sect_R_item(ihead+1:), & vect(1:nn*ndof) ) DO j = 1, ndof*nn hecMAT%B( iwk(j) ) = hecMAT%B( iwk(j) )+vect(j) END DO CYCLE END IF ! ----- Create local stiffness if(ic_type == 241 .or. ic_type == 242 .or. ic_type == 231 .or. ic_type == 232 ) then call TLOAD_C2( ic_type, nn, xx(1:nn), yy(1:nn), tt(1:nn), tt0(1:nn), & fstrSOLID%elements(icel)%gausses,pa1, iset, vect(1:nn*2) ) else if( ic_type == 361 ) then call TLOAD_C3D8Bbar & ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), & fstrSOLID%elements(icel)%gausses, vect(1:nn*ndof), cdsys_ID, coords ) ! 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) call TLOAD_C3 & ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), & fstrSOLID%elements(icel)%gausses, vect(1:nn*ndof), cdsys_ID, coords ) else if( ic_type == 741 .or. ic_type == 743 .or. ic_type == 731 ) then if( myrank == 0 ) then WRITE(IMSG,*) '*------------------------', & '-------------------*' WRITE(IMSG,*) ' Thermal loading option for shell elements', & 'not yet available.' WRITE(IMSG,*) '*------------------------', & '-------------------*' call hecmw_abort( hecmw_comm_get_comm()) endif endif ! ----- Add vector do j = 1, ndof*nn ! fstrSOLID%GL( iwk(j) )=fstrSOLID%GL( iwk(j) )+vect(j) hecMAT%B( iwk(j) ) = hecMAT%B( iwk(j) )+vect(j) enddo enddo enddo endif if( associated( fstrSOLID%contacts ) .and. fstrPARAM%contact_algo == kcaALagrange ) then do i = 1, size(fstrSOLID%contacts) call ass_contact_force( fstrSOLID%contacts(i), hecMESH%node, fstrSOLID%unode, hecMAT%B ) enddo endif end subroutine fstr_ass_load subroutine fstr_AddSPRING(cstep, sub_step, hecMESH, hecMAT, fstrSOLID, fstrPARAM) use m_fstr use m_static_lib integer, intent(in) :: cstep !< current step integer, intent(in) :: sub_step !< current sub_step type (hecmwST_matrix),intent(inout) :: hecMAT !< hecmw matrix type (hecmwST_local_mesh),intent(in) :: hecMESH !< hecmw mesh type (fstr_solid),intent(inout) :: fstrSOLID !< fstr_solid type (fstr_param),intent(inout) :: fstrPARAM !< analysis control parameters integer(kind=kint) :: grpid, ndof, ig0, ig, ityp, iS0, iE0, ik, in, idx, num real(kind=kreal) :: fval, fact ndof = hecMAT%NDOF do ig0= 1, fstrSOLID%SPRING_ngrp_tot grpid = fstrSOLID%SPRING_ngrp_GRPID(ig0) if( .not. fstr_isLoadActive( fstrSOLID, grpid, cstep ) ) cycle ig= fstrSOLID%SPRING_ngrp_ID(ig0) ityp= fstrSOLID%SPRING_ngrp_DOF(ig0) fval= fstrSOLID%SPRING_ngrp_val(ig0) if( fval < 0.d0 ) then num = fstrSOLID%step_ctrl(cstep)%num_substep fact = DFLOAT( num - sub_step ) / DFLOAT( num ) fval = - fval * fact endif iS0= hecMESH%node_group%grp_index(ig-1) + 1 iE0= hecMESH%node_group%grp_index(ig ) do ik= iS0, iE0 in = hecMESH%node_group%grp_item(ik) idx = ndof**2 * (in - 1) + ndof * (ityp - 1) + ityp hecMAT%D(idx) = hecMAT%D(idx) + fval enddo enddo end subroutine fstr_AddSPRING end module m_fstr_ass_load