hecmw_util_f.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
 
module hecmw_util
#ifndef HECMW_SERIAL
  use mpi
#endif
  implicit none
  private :: hecmw_petot,hecmw_rank,hecmw_comm,hecmw_group
  public
 
  integer(kind=4),parameter:: kint  = 4
  integer(kind=4),parameter:: kreal = 8
 
  integer(kind=kint),parameter :: hecmw_name_len     =   63
  integer(kind=kint),parameter :: hecmw_header_len   =  127
  integer(kind=kint),parameter :: hecmw_msg_len      =  255
  integer(kind=kint),parameter :: hecmw_filename_len = 1023
 
  integer(kind=kint),parameter :: hecmw_exit_success      =  0
  integer(kind=kint),parameter :: hecmw_exit_usage        =  1
  integer(kind=kint),parameter :: hecmw_exit_input        =  2
  integer(kind=kint),parameter :: hecmw_exit_model        =  3
  integer(kind=kint),parameter :: hecmw_exit_solver_setup =  4
  integer(kind=kint),parameter :: hecmw_exit_noconv       = 10
  integer(kind=kint),parameter :: hecmw_exit_internal     = 20
 
  integer(kind=kint),parameter :: hecmw_sum              = 46801
  integer(kind=kint),parameter :: hecmw_prod             = 46802
  integer(kind=kint),parameter :: hecmw_max              = 46803
  integer(kind=kint),parameter :: hecmw_min              = 46804
  integer(kind=kint),parameter :: hecmw_lor              = 46805
  integer(kind=kint),parameter :: hecmw_land             = 46806
  integer(kind=kint),parameter :: hecmw_integer          = 53951
  integer(kind=kint),parameter :: hecmw_single_precision = 53952
  integer(kind=kint),parameter :: hecmw_double_precision = 53953
  integer(kind=kint),parameter :: hecmw_character        = 53954
 
#ifndef HECMW_SERIAL
  integer(kind=kint),parameter :: hecmw_status_size = mpi_status_size
#else
  integer(kind=kint),parameter :: hecmw_status_size = 1
#endif
 
  integer(kind=kint) :: hecmw_petot,hecmw_rank,hecmw_comm,hecmw_group
  !C
  !C +---------------+
  !C | SECTION info. |
  !C +---------------+
  !C===
  type hecmwst_section
    integer(kind=kint)         :: n_sect
    integer(kind=kint),pointer :: sect_type(:)
    integer(kind=kint),pointer :: sect_opt(:)
    integer(kind=kint),pointer :: sect_mat_id_index(:)
    integer(kind=kint),pointer :: sect_mat_id_item(:)
    integer(kind=kint),pointer :: sect_i_index(:)
    integer(kind=kint),pointer :: sect_i_item(:)
    integer(kind=kint),pointer :: sect_r_index(:)
    real(kind=kreal),pointer   :: sect_r_item(:)
    integer(kind=kint),pointer :: sect_orien_id(:) => null()
  end type hecmwst_section
 
  !C      for hecmwST_section%sect_type
  integer(kind=kint),parameter :: hecmw_sect_type_solid     = 1
  integer(kind=kint),parameter :: hecmw_sect_type_shell     = 2
  integer(kind=kint),parameter :: hecmw_sect_type_beam      = 3
  integer(kind=kint),parameter :: hecmw_sect_type_interface = 4
  !C      for hecmwST_section%sect_opt
  integer(kind=kint),parameter :: hecmw_sect_opt_pstress      =  0
  integer(kind=kint),parameter :: hecmw_sect_opt_pstrain      =  1
  integer(kind=kint),parameter :: hecmw_sect_opt_asymmetry    =  2
  integer(kind=kint),parameter :: hecmw_sect_opt_pstress_ri   = 10
  integer(kind=kint),parameter :: hecmw_sect_opt_pstrain_ri   = 11
  integer(kind=kint),parameter :: hecmw_sect_opt_asymmetry_ri = 12
  !C===
 
  !C
  !C +----------------+
  !C | MATERIAL info. |
  !C +----------------+
  !C===
  type hecmwst_material
    integer(kind=kint)         :: n_mat
    integer(kind=kint)         :: n_mat_item
    integer(kind=kint)         :: n_mat_subitem
    integer(kind=kint)         :: n_mat_table
    character(HECMW_NAME_LEN),pointer :: mat_name(:)
    integer(kind=kint),pointer :: mat_item_index(:)
    integer(kind=kint),pointer :: mat_subitem_index(:)
    integer(kind=kint),pointer :: mat_table_index(:)
    real(kind=kreal),pointer   :: mat_val(:)
    real(kind=kreal),pointer   :: mat_temp(:)
  end type hecmwst_material
  !C===
 
  !C
  !C +-----------+
  !C | MPC info. |
  !C +-----------+
  !C===
  type hecmwst_mpc
    integer(kind=kint)         :: n_mpc
    integer(kind=kint),pointer :: mpc_index(:)
    integer(kind=kint),pointer :: mpc_item(:)
    integer(kind=kint),pointer :: mpc_dof(:)
    real(kind=kreal),pointer   :: mpc_val(:)
    real(kind=kreal),pointer   :: mpc_const(:)
  end type hecmwst_mpc
  !C===
 
  !C
  !C +-----------+
  !C | AMPLITUDE |
  !C +-----------+
  !C===
  type hecmwst_amplitude
    integer(kind=kint)         :: n_amp
    character(len=HECMW_NAME_LEN),pointer :: amp_name(:)
    integer(kind=kint),pointer :: amp_type_definition(:)
    integer(kind=kint),pointer :: amp_type_time(:)
    integer(kind=kint),pointer :: amp_type_value(:)
    integer(kind=kint),pointer :: amp_index(:)
    real(kind=kreal),pointer   :: amp_val(:)
    real(kind=kreal),pointer   :: amp_table(:)
  end type hecmwst_amplitude
 
  !C      for hecmwST_amplitude%amp_type_definition
  integer(kind=kint),parameter :: hecmw_amp_typedef_tabular  = 1
  !C      for hecmwST_amplitude%amp_type_time
  integer(kind=kint),parameter :: hecmw_amp_typetime_step    = 1
  integer(kind=kint),parameter :: hecmw_amp_typetime_total   = 2
  !C      for hecmwST_amplitude%amp_type_value
  integer(kind=kint),parameter :: hecmw_amp_typeval_relative = 1
  integer(kind=kint),parameter :: hecmw_amp_typeval_absolute = 2
  !C===
 
  !C
  !C +-----------+
  !C | NODE grp. |
  !C +-----------+
  !C===
  type hecmwst_node_grp
    integer(kind=kint)         :: n_grp
    integer(kind=kint)         :: n_bc
    character(HECMW_NAME_LEN),pointer :: grp_name(:)
    integer(kind=kint),pointer :: grp_index(:)
    integer(kind=kint),pointer :: grp_item(:)
    integer(kind=kint),pointer :: bc_grp_id(:)
    integer(kind=kint),pointer :: bc_grp_type(:)
    integer(kind=kint),pointer :: bc_grp_index(:)
    integer(kind=kint),pointer :: bc_grp_dof(:)
    real(kind=kreal),pointer   :: bc_grp_val(:)
  end type hecmwst_node_grp
 
  !C      for hecmwST_node_grp%bc_grp_type
  integer(kind=kint),parameter :: hecmw_bcgrptype_desplacement = 1
  integer(kind=kint),parameter :: hecmw_bcgrptype_flux         = 2
  !C===
 
  !C
  !C +-----------+
  !C | ELEM grp. |
  !C +-----------+
  !C===
  type hecmwst_elem_grp
    integer(kind=kint)         :: n_grp
    integer(kind=kint)         :: n_bc
    character(HECMW_NAME_LEN),pointer :: grp_name(:)
    integer(kind=kint),pointer :: grp_index(:)
    integer(kind=kint),pointer :: grp_item(:)
    integer(kind=kint),pointer :: bc_grp_id(:)
    integer(kind=kint),pointer :: bc_grp_type(:)
    integer(kind=kint),pointer :: bc_grp_index(:)
    real(kind=kreal),pointer   :: bc_grp_val(:)
  end type hecmwst_elem_grp
 
  !C      for hecmwST_elem_grp%bc_grp_type
  !C        integer(kind=kint),parameter :: HECMW_BCGRPTYPE_DESPLACEMENT = 1
  !C        integer(kind=kint),parameter :: HECMW_BCGRPTYPE_FLUX         = 2
  !C===
 
  !C
  !C +-----------+
  !C | SURF grp. |
  !C +-----------+
  !C===
  type hecmwst_surf_grp
    integer(kind=kint)         :: n_grp
    integer(kind=kint)         :: n_bc
    character(HECMW_NAME_LEN),pointer:: grp_name(:)
    integer(kind=kint),pointer :: grp_index(:)
    integer(kind=kint),pointer :: grp_item(:)
    integer(kind=kint),pointer :: bc_grp_id(:)
    integer(kind=kint),pointer :: bc_grp_type(:)
    integer(kind=kint),pointer :: bc_grp_index(:)
    real(kind=kreal),pointer   :: bc_grp_val(:)
  end type hecmwst_surf_grp
 
  !C      for hecmwST_surf_grp%bc_grp_type
  !C        integer(kind=kint),parameter :: HECMW_BCGRPTYPE_DESPLACEMENT = 1
  !C        integer(kind=kint),parameter :: HECMW_BCGRPTYPE_FLUX         = 2
  !C
 
  !C
  !C +---------+
  !C | CONTACT |
  !C +---------+
  !C===
  type hecmwst_contact_pair
    integer(kind=kint)         :: n_pair
    character(HECMW_NAME_LEN),pointer:: name(:)
    integer(kind=kint),pointer :: type(:)
    integer(kind=kint),pointer :: slave_grp_id(:)
    integer(kind=kint),pointer :: master_grp_id(:)
    integer(kind=kint),pointer :: slave_orisgrp_id(:)
  end type hecmwst_contact_pair
 
  !C      for hecmwST_contact_pair%type
  integer(kind=kint),parameter :: hecmw_contact_type_node_surf = 1
  integer(kind=kint),parameter :: hecmw_contact_type_surf_surf = 2
  integer(kind=kint),parameter :: hecmw_contact_type_node_elem = 3
  !C===
 
  !C
  !C +----------------+
  !C | REFINE Origin. |
  !C +----------------+
  !C===
  type hecmwst_refine_origin
    integer(kind=kint),pointer :: index(:)
    integer(kind=kint),pointer :: item_index(:)
    integer(kind=kint),pointer :: item_item(:)
  end type hecmwst_refine_origin
  !C===
 
  !C
  !C +------------------+
  !C | LOCAL MESH info. |
  !C +------------------+
  !C===
  type hecmwst_local_mesh
 
    !C
    !C-- FILES, GENERAL
    !C
    character(HECMW_FILENAME_LEN)         :: gridfile
    character(HECMW_FILENAME_LEN),pointer :: files(:)
    character(HECMW_HEADER_LEN)           :: header
    integer(kind=kint) :: hecmw_flag_adapt
    integer(kind=kint) :: hecmw_flag_initcon
    integer(kind=kint) :: hecmw_n_file
    integer(kind=kint) :: hecmw_flag_parttype
    integer(kind=kint) :: hecmw_flag_partdepth
    integer(kind=kint) :: hecmw_flag_version
    integer(kind=kint) :: hecmw_flag_partcontact
    real(kind=kreal)   :: zero_temp
 
    !C
    !C-- NODE
    integer(kind=kint)         :: n_node
    integer(kind=kint)         :: n_node_gross
    !         For parallel contact with hanging slave nodes
    integer(kind=kint)         :: nn_middle
    integer(kind=kint)         :: nn_internal
    integer(kind=kint)         :: n_dof
    integer(kind=kint)         :: n_dof_grp
    integer(kind=kint)         :: n_dof_tot
    real(kind=kreal),pointer   :: node(:)
    integer(kind=kint),pointer :: node_id(:)
    integer(kind=kint),pointer :: global_node_id(:)
    integer(kind=kint),pointer :: node_val_index(:)
    real(kind=kreal),pointer   :: node_val_item(:)
    integer(kind=kint),pointer :: node_dof_index(:)
    integer(kind=kint),pointer :: node_dof_item(:)
    integer(kind=kint),pointer :: node_init_val_index(:)
    real(kind=kreal),pointer   :: node_init_val_item(:)
    integer(kind=kint),pointer :: node_internal_list(:)
    !C
    !C-- ELEMENT
    !C
    integer(kind=kint)         :: n_elem
    integer(kind=kint)         :: n_elem_gross
    integer(kind=kint)         :: ne_internal
    integer(kind=kint)         :: n_elem_type
    integer(kind=kint)         :: n_elem_mat_id
    integer(kind=kint),pointer :: elem_type_index(:)
    integer(kind=kint),pointer :: elem_type_item(:)
    integer(kind=kint),pointer :: elem_type(:)
    integer(kind=kint),pointer :: section_id(:)
    integer(kind=kint),pointer :: elem_mat_id_index(:)
    integer(kind=kint),pointer :: elem_mat_id_item(:)
    integer(kind=kint),pointer :: elem_node_index(:)
    integer(kind=kint),pointer :: elem_node_item(:)
    integer(kind=kint),pointer :: elem_id(:)
    integer(kind=kint),pointer :: global_elem_id(:)
    integer(kind=kint),pointer :: elem_internal_list(:)
    integer(kind=kint),pointer :: elem_mat_int_index(:)
    real(kind=kreal),pointer   :: elem_mat_int_val(:)
    integer(kind=kint),pointer :: elem_val_index(:)
    real(kind=kreal),pointer   :: elem_val_item(:)
    !integer(kind=kint)         :: is_33shell
    !integer(kind=kint)         :: is_33beam
    !integer(kind=kint)         :: is_heat
    !C
    !C-- COMMUNICATION
    !C
    integer(kind=kint)         :: zero
    integer(kind=kint)         :: mpi_comm
    integer(kind=kint)         :: petot
    integer(kind=kint)         :: pesmptot
    integer(kind=kint)         :: my_rank
    integer(kind=kint)         :: errnof
    integer(kind=kint)         :: n_subdomain
    integer(kind=kint)         :: n_neighbor_pe
    integer(kind=kint),pointer :: neighbor_pe(:)
    integer(kind=kint),pointer :: import_index(:)
    integer(kind=kint),pointer :: import_item(:)
    integer(kind=kint),pointer :: export_index(:)
    integer(kind=kint),pointer :: export_item(:)
    integer(kind=kint),pointer :: shared_index(:)
    integer(kind=kint),pointer :: shared_item(:)
 
    !C
    !C-- ADAPTATION
    !C
    integer(kind=kint)         :: coarse_grid_level
    integer(kind=kint)         :: n_adapt
    integer(kind=kint),pointer :: when_i_was_refined_node(:)
    integer(kind=kint),pointer :: when_i_was_refined_elem(:)
    integer(kind=kint),pointer :: adapt_parent_type(:)
    integer(kind=kint),pointer :: adapt_type (:)
    integer(kind=kint),pointer :: adapt_level(:)
    integer(kind=kint),pointer :: adapt_parent(:)
    integer(kind=kint),pointer :: adapt_children_index(:)
    integer(kind=kint),pointer :: adapt_children_item(:)
 
    integer(kind=kint) :: nn_array, ne_array, nx_array
    integer(kind=kint) :: n_adapt_edge, n_adapt_edge_global
    integer(kind=kint) :: n_adapt_act_node, n_adapt_act_edge
    integer(kind=kint) :: n_adapt_act_elem, n_adapt_act_elem_cur
    integer(kind=kint) :: n_adapt_elem_341, n_adapt_elem_351
    integer(kind=kint) :: n_adapt_elem_341_cur, n_adapt_elem_351_cur
    integer(kind=kint) :: n_adapt_act_elem_341, n_adapt_act_elem_351
 
    integer(kind=kint) :: n_adapt_node_cur, nn_adapt_internal_cur
    integer(kind=kint) :: n_adapt_node_old, nn_adapt_internal_old
    integer(kind=kint) :: n_adapt_elem_cur, n_adapt_elem_old
 
    integer(kind=kint), pointer :: adapt_edge_node(:), adapt_mid_edge (:)
    integer(kind=kint), pointer :: adapt_iemb     (:), adapt_edge_home(:)
    integer(kind=kint), pointer :: adapt_act_edge (:)
 
    integer(kind=kint), pointer ::                                           &
      &           adapt_import_edge_index(:), adapt_import_edge_item (:),&
      &           adapt_export_edge_index(:), adapt_export_edge_item (:),&
      &           adapt_import_elem_index(:), adapt_import_elem_item (:),&
      &           adapt_export_elem_index(:), adapt_export_elem_item (:),&
      &           adapt_import_new_index (:), adapt_import_new_item  (:),&
      &           adapt_export_new_index (:), adapt_export_new_item  (:)
    integer(kind=kint), pointer :: rev_neighbor_pe(:)
    integer(kind=kint), pointer :: adapt_act_elem_341(:)
    integer(kind=kint), pointer :: adapt_act_elem_351(:)
    integer(kind=kint), pointer :: adapt_oldtonew_node(:), adapt_newtoold_node(:)
    integer(kind=kint), pointer :: adapt_oldtonew_elem(:), adapt_newtoold_elem(:)
    integer(kind=kint), pointer :: adapt_iwk(:), adapt_children_local(:)
 
    !C
    !C-- REFINEMENT
    !C
    integer(kind=kint)         :: n_refine
    integer(kind=kint),pointer :: node_old2new(:)
    integer(kind=kint),pointer :: node_new2old(:)
    integer(kind=kint),pointer :: elem_old2new(:)
    integer(kind=kint),pointer :: elem_new2old(:)
    integer(kind=kint),pointer :: n_node_refine_hist(:)
 
    !C
    !C-- ETC.
    !C
    type (hecmwst_section)   :: section
    type (hecmwst_material)  :: material
    type (hecmwst_mpc)       :: mpc
    type (hecmwst_amplitude) :: amp
    type (hecmwst_node_grp)  :: node_group
    type (hecmwst_elem_grp)  :: elem_group
    type (hecmwst_surf_grp)  :: surf_group
    type (hecmwst_contact_pair):: contact_pair
    type (hecmwst_refine_origin):: refine_origin
 
  end type hecmwst_local_mesh
 
  !C      for hecmwST_local_mesh%hecmw_flag_parttype
  integer(kind=kint),parameter :: hecmw_flag_parttype_unknown   = 0
  integer(kind=kint),parameter :: hecmw_flag_parttype_nodebased = 1
  integer(kind=kint),parameter :: hecmw_flag_parttype_elembased = 2
 
  !C      for hecmwST_local_mesh%hecmw_flag_partcontact
  integer(kind=kint),parameter :: hecmw_flag_partcontact_unknown    = 0
  integer(kind=kint),parameter :: hecmw_flag_partcontact_aggregate  = 1
  integer(kind=kint),parameter :: hecmw_flag_partcontact_distribute = 2
  integer(kind=kint),parameter :: hecmw_flag_partcontact_simple     = 3
 
  !C      the mode above occupies bits 0-7, the ownership scheme bit 8.  OWNER_MASTER must stay zero so that a
  !C      distributed mesh holding a bare mode value reads back as the master-owner scheme.  Test with
  !C      hecmw_partcontact_get_mode()/_get_owner(), never for equality against the whole flag.
  integer(kind=kint),parameter :: hecmw_flag_partcontact_mode_mask    = 255
  integer(kind=kint),parameter :: hecmw_flag_partcontact_owner_mask   = 256
  integer(kind=kint),parameter :: hecmw_flag_partcontact_owner_master = 0
  integer(kind=kint),parameter :: hecmw_flag_partcontact_owner_slave  = 256
 
  !C
  !C +--------+
  !C | MATRIX |
  !C +--------+
  !C===
  type hecmwst_matrix_comm
    integer(kind=kint)                        :: zero
    integer(kind=kint)                        :: hecmw_comm
    integer(kind=kint)                        :: petot
    integer(kind=kint)                        :: pesmptot
    integer(kind=kint)                        :: my_rank
    integer(kind=kint)                        :: errnof
    integer(kind=kint)                        :: n_subdomain
    integer(kind=kint)                        :: n_neighbor_pe
    integer(kind=kint), dimension(:), pointer :: neighbor_pe
    integer(kind=kint), dimension(:), pointer :: import_index
    integer(kind=kint), dimension(:), pointer :: import_item
    integer(kind=kint), dimension(:), pointer :: export_index
    integer(kind=kint), dimension(:), pointer :: export_item
    integer(kind=kint), dimension(:), pointer :: shared_index
    integer(kind=kint), dimension(:), pointer :: shared_item
  end type hecmwst_matrix_comm
 
  type hecmwst_index_value_pair
    integer(kind=kint) :: i
    integer(kind=kint) :: j
    real(kind=kreal), dimension(3,3) :: val
  end type hecmwst_index_value_pair
 
  type hecmwst_matrix_lagrange
    integer(kind=kint) :: num_lagrange = 0
    integer(kind=kint) :: numl_lagrange = 0
    integer(kind=kint) :: numu_lagrange = 0
 
    integer(kind=kint), pointer  :: indexl_lagrange(:) => null() 
    integer(kind=kint), pointer  :: indexu_lagrange(:) => null() 
 
    integer(kind=kint), pointer  :: iteml_lagrange(:) => null() 
    integer(kind=kint), pointer  :: itemu_lagrange(:) => null() 
 
    real(kind=kreal),    pointer  :: al_lagrange(:) => null() 
    real(kind=kreal),    pointer  :: au_lagrange(:) => null() 
    real(kind=kreal),    pointer  :: d_lagrange(:) => null() 
 
    real(kind=kreal),    pointer  :: lagrange(:) => null() 
 
    integer(kind=kint), pointer  :: lag_node_table(:) => null() 
  end type hecmwst_matrix_lagrange
 
  type hecmwst_matrix
    integer(kind=kint) ::  n, np, npl, npu, npa, ndof
    real(kind=kreal), pointer :: d(:), b(:), x(:)
    real(kind=kreal), pointer :: al(:), au(:)
    real(kind=kreal), pointer :: a(:)
    integer(kind=kint), pointer :: indexl(:), indexu(:)
    integer(kind=kint), pointer :: indexa(:)
    integer(kind=kint), pointer ::  iteml(:),  itemu(:)
    integer(kind=kint), pointer ::  itema(:)
    integer(kind=kint ), dimension(100) :: iarray
    real   (kind=kreal), dimension(100) :: rarray
    logical :: symmetric = .true.
  end type hecmwst_matrix
contains
 
  !C
  !C***
  !C*** HECMW_PARTCONTACT_GET_MODE / HECMW_PARTCONTACT_GET_OWNER
  !C***
  !C
  !C    UNPACK hecmw_flag_partcontact, so that its bit layout stays internal to HEC-MW
  !C
  pure integer(kind=kint) function hecmw_partcontact_get_mode( flag_partcontact )
    integer(kind=kint), intent(in) :: flag_partcontact
    hecmw_partcontact_get_mode = iand( flag_partcontact, hecmw_flag_partcontact_mode_mask )
  end function hecmw_partcontact_get_mode
 
  pure integer(kind=kint) function hecmw_partcontact_get_owner( flag_partcontact )
    integer(kind=kint), intent(in) :: flag_partcontact
    hecmw_partcontact_get_owner = iand( flag_partcontact, hecmw_flag_partcontact_owner_mask )
  end function hecmw_partcontact_get_owner
 
  !C
  !C***
  !C*** HECMW_HAS_INTERNAL_NODE
  !C***
  !C
  !C    TEST WHETHER AN ELEMENT TOUCHES AN INTERNAL NODE OF THIS RANK
  !C
  !C    Local node numbering places the internal nodes first, so a local node id
  !C    greater than nn_internal belongs to another rank.
  !C
  pure logical function hecmw_has_internal_node( hecMESH, nn, nodLOCAL )
    type(hecmwst_local_mesh), intent(in) :: hecmesh
    integer(kind=kint), intent(in)       :: nn
    integer(kind=kint), intent(in)       :: nodlocal(:)
    integer(kind=kint) :: j
 
    hecmw_has_internal_node = .true.
    do j = 1, nn
      if( nodlocal(j) <= hecmesh%nn_internal ) return
    enddo
    hecmw_has_internal_node = .false.
  end function hecmw_has_internal_node
 
  !C
  !C***
  !C*** HECMW_INIT
  !C***
  !C
  !C    INIT. HECMW-FEM process's
  !C
  subroutine hecmw_init
    character(len=HECMW_FILENAME_LEN):: ctrlfile = "hecmw_ctrl.dat"
    call hecmw_init_ex(ctrlfile)
  end subroutine hecmw_init
 
  subroutine hecmw_comm_init
    integer(kind=kint) :: ierr
#ifndef HECMW_SERIAL
    !call MPI_INIT (ierr)
    call mpi_comm_size (mpi_comm_world, hecmw_petot, ierr)
    call mpi_comm_rank (mpi_comm_world, hecmw_rank,  ierr)
    call mpi_comm_dup  (mpi_comm_world, hecmw_comm,  ierr)
    call mpi_comm_group(mpi_comm_world, hecmw_group, ierr)
#else
    hecmw_petot=1
    hecmw_rank=0
    hecmw_comm=0
    hecmw_group=0
    ierr=0
#endif
    if(ierr /= 0) then
      call  hecmw_abort( hecmw_comm_get_comm( ) )
    endif
  end subroutine hecmw_comm_init
 
  subroutine hecmw_init_nocontrol()
    call hecmw_comm_init()
    call hecmw_comm_init_if(hecmw_comm, hecmw_petot, hecmw_rank, hecmw_group)
  end subroutine
 
  !C
  !C***
  !C*** HECMW_INIT_EX
  !C***
  !C
  !C    INIT. HECMW-FEM process's
  !C
  subroutine hecmw_init_ex(ctrlfile)
    character(len=HECMW_FILENAME_LEN):: ctrlfile
    integer(kind=kint) :: ierr
 
    call hecmw_comm_init()
    call hecmw_comm_init_if(hecmw_comm, hecmw_petot, hecmw_rank, hecmw_group)
 
    call hecmw_ctrl_init_ex_if(ctrlfile, ierr)
    if(ierr /= 0) then
      call  hecmw_abort( hecmw_comm_get_comm( ), hecmw_exit_input )
    endif
    !      call hecmw_couple_comm_init_if(ierr)
    !      if(ierr /= 0) then
    !          call  hecmw_abort( hecmw_comm_get_comm( ) )
    !      endif
 
  end subroutine hecmw_init_ex
 
 
  !C
  !C***
  !C*** HECMW_FINALIZE
  !C***
  !C
  !C    FINALIZE. HECMW-FEM process's
  !C
  subroutine hecmw_finalize
    integer(kind=kint) :: ierr
 
    call hecmw_ctrl_finalize_if()
 
#ifndef HECMW_SERIAL
    call mpi_finalize(ierr)
#endif
 
  end subroutine hecmw_finalize
 
 
  !C******************** MPI WRAPPER SUBROUTINES ************************
  !C
  !C***
  !C*** HECMW_ABORT
  !C***
  !C
  subroutine hecmw_abort(comm, code)
    integer(kind=kint) :: comm
    integer(kind=kint), intent(in), optional :: code
    integer(kind=kint) :: errorcode, ierror
 
    errorcode = hecmw_exit_internal
    if( present(code) ) errorcode = code
 
    ! MPI_ABORT does not flush the Fortran unit buffers; do it here so that
    ! the message written just before the abort survives.
    flush(6)
 
#ifndef HECMW_SERIAL
    call mpi_abort(comm, errorcode, ierror)
#else
    stop errorcode
#endif
  end subroutine hecmw_abort
 
  !C
  !C***
  !C*** HECMW_WTIME
  !C***
  !C
  function hecmw_wtime()
    real(kind=kreal) hecmw_wtime
    external  hecmw_wtime_fi
    real(kind=kreal) hecmw_wtime_fi
    hecmw_wtime = hecmw_wtime_fi()
  end function hecmw_wtime
  !C
  !C***
  !C*** HECMW_WTICK
  !C***
  !C
  function hecmw_wtick()
    real(kind=kreal) hecmw_wtick
    external hecmw_wtick_fi
    real(kind=kreal) hecmw_wtick_fi
    hecmw_wtick = hecmw_wtick_fi()
  end function hecmw_wtick
  !C
  !C***
  !C*** HECMW_COMM_GET_COMM
  !C***
  !C
  function hecmw_comm_get_comm() result(comm)
    integer(kind=kint) :: comm
 
    comm = hecmw_comm
  end function hecmw_comm_get_comm
 
  !C
  !C***
  !C*** HECMW_COMM_GET_RANK
  !C***
  !C
  function hecmw_comm_get_rank() result(rank)
    integer(kind=kint) :: rank
 
    rank = hecmw_rank
  end function hecmw_comm_get_rank
 
  !C
  !C***
  !C*** HECMW_COMM_GET_SIZE
  !C***
  !C
  function hecmw_comm_get_size() result(comm_size)
    integer(kind=kint) :: comm_size
 
    comm_size = hecmw_petot
  end function hecmw_comm_get_size
 
 
  !C***************  NULL POINTER SETTING UTILITY ****************
 
  subroutine hecmw_nullify_section( P )
    type( hecmwst_section ) :: P
    nullify( p%sect_type )
    nullify( p%sect_opt )
    nullify( p%sect_mat_ID_index )
    nullify( p%sect_mat_ID_item )
    nullify( p%sect_I_index )
    nullify( p%sect_I_item )
    nullify( p%sect_R_index )
    nullify( p%sect_R_item )
  end subroutine hecmw_nullify_section
 
  subroutine hecmw_nullify_material( P )
    type( hecmwst_material ) :: P
    nullify( p%mat_name )
    nullify( p%mat_item_index )
    nullify( p%mat_subitem_index )
    nullify( p%mat_table_index )
    nullify( p%mat_val )
    nullify( p%mat_temp )
  end subroutine hecmw_nullify_material
 
  subroutine hecmw_nullify_mpc( P )
    type( hecmwst_mpc ) :: P
    nullify( p%mpc_index )
    nullify( p%mpc_item )
    nullify( p%mpc_dof )
    nullify( p%mpc_val )
    nullify( p%mpc_const )
  end subroutine hecmw_nullify_mpc
 
  subroutine hecmw_initialize_mpc( mpc, n_mpc, n_item )
    type( hecmwst_mpc ), intent(inout) :: mpc
    integer(kind=kint), intent(in)    :: n_mpc
    integer(kind=kint), intent(in)    :: n_item
 
    mpc%n_mpc = n_mpc
    allocate( mpc%mpc_index(0:n_mpc) )
    allocate( mpc%mpc_item(n_item) )
    allocate( mpc%mpc_dof(n_item) )
    allocate( mpc%mpc_val(n_item) )
  end subroutine
 
  subroutine hecmw_finalize_mpc( P )
    type( hecmwst_mpc ) :: P
    if( associated(p%mpc_index) ) deallocate( p%mpc_index )
    if( associated(p%mpc_item) )  deallocate( p%mpc_item )
    if( associated(p%mpc_dof) )   deallocate( p%mpc_dof )
    if( associated(p%mpc_val) )   deallocate( p%mpc_val )
  end subroutine hecmw_finalize_mpc
 
  subroutine hecmw_nullify_amplitude( P )
    type( hecmwst_amplitude ) :: P
    nullify( p%amp_name )
    nullify( p%amp_type_definition )
    nullify( p%amp_type_time )
    nullify( p%amp_type_value )
    nullify( p%amp_index )
    nullify( p%amp_val )
    nullify( p%amp_table )
  end subroutine hecmw_nullify_amplitude
 
  subroutine hecmw_nullify_node_grp( P )
    type( hecmwst_node_grp ) :: P
    nullify( p%grp_name )
    nullify( p%grp_index )
    nullify( p%grp_item )
    nullify( p%bc_grp_ID )
    nullify( p%bc_grp_type )
    nullify( p%bc_grp_index )
    nullify( p%bc_grp_dof )
    nullify( p%bc_grp_val )
  end subroutine hecmw_nullify_node_grp
 
  subroutine hecmw_nullify_elem_grp( P )
    type( hecmwst_elem_grp ) :: P
    nullify( p%grp_name )
    nullify( p%grp_index )
    nullify( p%grp_item )
    nullify( p%bc_grp_ID )
    nullify( p%bc_grp_type )
    nullify( p%bc_grp_index )
    nullify( p%bc_grp_val )
  end subroutine hecmw_nullify_elem_grp
 
  subroutine hecmw_nullify_surf_grp( P )
    type( hecmwst_surf_grp ) :: P
    nullify( p%grp_name )
    nullify( p%grp_index )
    nullify( p%grp_item )
    nullify( p%bc_grp_ID )
    nullify( p%bc_grp_type )
    nullify( p%bc_grp_index )
    nullify( p%bc_grp_val )
  end subroutine hecmw_nullify_surf_grp
 
  subroutine hecmw_nullify_contact_pair( P )
    type( hecmwst_contact_pair ) :: P
    nullify( p%name )
    nullify( p%type )
    nullify( p%slave_grp_id )
    nullify( p%slave_orisgrp_id )
    nullify( p%master_grp_id )
  end subroutine hecmw_nullify_contact_pair
 
  subroutine hecmw_nullify_refine_origin( P )
    type( hecmwst_refine_origin ) :: P
    nullify( p%index )
    nullify( p%item_index )
    nullify( p%item_item )
  end subroutine hecmw_nullify_refine_origin
 
  subroutine hecmw_nullify_mesh( P )
    type( hecmwst_local_mesh ) :: P
    nullify( p%files )
    nullify( p%node )
    nullify( p%node_ID )
    nullify( p%global_node_ID )
    nullify( p%node_val_index )
    nullify( p%node_val_item )
    nullify( p%node_dof_index )
    nullify( p%node_dof_item )
    nullify( p%node_init_val_index )
    nullify( p%node_init_val_item )
    nullify( p%node_internal_list )
    nullify( p%elem_type_index )
    nullify( p%elem_type_item )
    nullify( p%elem_type )
    nullify( p%section_ID )
    nullify( p%elem_mat_ID_index )
    nullify( p%elem_mat_ID_item )
    nullify( p%elem_node_index )
    nullify( p%elem_node_item )
    nullify( p%elem_ID )
    nullify( p%global_elem_ID )
    nullify( p%elem_internal_list )
    nullify( p%elem_mat_int_index )
    nullify( p%elem_mat_int_val )
    nullify( p%elem_val_index )
    nullify( p%elem_val_item )
    nullify( p%neighbor_pe )
    nullify( p%import_index )
    nullify( p%import_item )
    nullify( p%export_index )
    nullify( p%export_item )
    nullify( p%shared_index )
    nullify( p%shared_item )
    nullify( p%when_i_was_refined_node )
    nullify( p%when_i_was_refined_elem )
    nullify( p%adapt_parent_type )
    nullify( p%adapt_type  )
    nullify( p%adapt_level )
    nullify( p%adapt_parent )
    nullify( p%adapt_children_index )
    nullify( p%adapt_children_item )
    nullify( p%adapt_edge_node )
    nullify( p%adapt_mid_edge  )
    nullify( p%adapt_iemb      )
    nullify( p%adapt_edge_home )
    nullify( p%adapt_act_edge  )
    nullify( p%adapt_import_edge_index )
    nullify( p%adapt_import_edge_item  )
    nullify( p%adapt_export_edge_index )
    nullify( p%adapt_export_edge_item  )
    nullify( p%adapt_import_elem_index )
    nullify( p%adapt_import_elem_item  )
    nullify( p%adapt_export_elem_index )
    nullify( p%adapt_export_elem_item  )
    nullify( p%adapt_import_new_index  )
    nullify( p%adapt_import_new_item   )
    nullify( p%adapt_export_new_index  )
    nullify( p%adapt_export_new_item   )
    nullify( p%rev_neighbor_pe )
    nullify( p%adapt_act_elem_341 )
    nullify( p%adapt_act_elem_351 )
    nullify( p%adapt_OLDtoNEW_node )
    nullify( p%adapt_NEWtoOLD_node )
    nullify( p%adapt_OLDtoNEW_elem )
    nullify( p%adapt_NEWtoOLD_elem )
    nullify( p%adapt_IWK )
    nullify( p%adapt_children_local )
    nullify( p%node_old2new )
    nullify( p%node_new2old )
    nullify( p%elem_old2new )
    nullify( p%elem_new2old )
    nullify( p%n_node_refine_hist )
 
    call hecmw_nullify_section( p%section )
    call hecmw_nullify_material( p%material )
    call hecmw_nullify_mpc( p%mpc )
    call hecmw_nullify_amplitude( p%amp )
    call hecmw_nullify_node_grp( p%node_group )
    call hecmw_nullify_elem_grp( p%elem_group )
    call hecmw_nullify_surf_grp( p%surf_group )
    call hecmw_nullify_contact_pair( p%contact_pair )
    call hecmw_nullify_refine_origin( p%refine_origin )
 
  end subroutine hecmw_nullify_mesh
 
 
  subroutine hecmw_nullify_matrix_comm( P )
    type( hecmwst_matrix_comm ) :: P
    nullify( p%neighbor_pe )
    nullify( p%import_index )
    nullify( p%import_item )
    nullify( p%export_index )
    nullify( p%export_item )
    nullify( p%shared_index )
    nullify( p%shared_item )
  end subroutine hecmw_nullify_matrix_comm
 
  subroutine hecmw_nullify_matrix( P )
    type( hecmwst_matrix ) :: P
    nullify( p%D )
    nullify( p%B )
    nullify( p%X )
    nullify( p%AL )
    nullify( p%AU )
    nullify( p%A )
    nullify( p%indexL )
    nullify( p%indexU )
    nullify( p%indexA )
    nullify( p%itemL )
    nullify( p%itemU )
    nullify( p%itemA )
  end subroutine hecmw_nullify_matrix
 
  subroutine hecmw_print_matrix( fname, P )
    character(len=*), intent(in)       :: fname
    type( hecmwst_matrix ), intent(in) :: P
 
    integer :: i, nf, nBlock
    integer :: io_stat
    nf = 777
    nblock = p%NDOF * p%NDOF
    open( unit=nf, file=fname, iostat=io_stat)
    if( io_stat /= 0 ) return
    write( nf, * ) p%N,p%NP,p%NPL,p%NPU,p%NDOF
    !---- index
    do i=0, p%NP
      write( nf,* )  p%indexL(i), p%indexU(i)
    enddo
    !---- itemL, AL
    do i=1, p%NPL
      write( nf,* ) p%itemL(i)
    enddo
    do i=1,nblock * p%NPL
      write( nf,* ) p%AL(i)
    enddo
    !---- itemU, AU
    do i=1,p%NPU
      write( nf,* ) p%itemU(i)
    enddo
    do i=1,nblock * p%NPU
      write( nf,* ) p%AU(i)
    enddo
    !---- D
    do i=1,nblock * p%NP
      write( nf, * ) p%D(i)
    enddo
    !---- B
    do i=1,p%NDOF * p%NP
      write( nf, * ) p%B(i)
    enddo
 
    close( nf )
  end subroutine
 
  subroutine hecmw_read_matrix( fname, P )
    character(len=*), intent(in)        :: fname
    type( hecmwst_matrix ), intent(out) :: P
 
    integer :: i, nf, nBlock, istat
    nf = 777
    open( unit=nf, file=fname, status='old', iostat= istat)
    if(istat /= 0) then
      print *, "cannot open file ",fname
      stop
    endif
    read( nf, * ) p%N,p%NP,p%NPL,p%NPU,p%NDOF
    nblock = p%NDOF * p%NDOF
 
    !----It is supposing of array are not allocated yet
    allocate( p%indexL(0:p%NP), p%indexU(0:p%NP) )
    allocate( p%itemL(p%NPL), p%itemU(p%NPU) )
    allocate( p%AL(p%NPL*nblock), p%AU(p%NPU*nblock) )
    allocate( p%D(p%NP*nblock) )
    allocate( p%B(p%NDOF*p%NP) )
    !---- index
    do i=0, p%NP
      read( nf,* )  p%indexL(i), p%indexU(i)
    enddo
    !---- itemL, AL
    do i=1, p%NPL
      read( nf,* ) p%itemL(i)
    enddo
    do i=1,nblock * p%NPL
      read( nf,* ) p%AL(i)
    enddo
    !---- itemU, AU
    do i=1,p%NPU
      read( nf,* ) p%itemU(i)
    enddo
    do i=1,nblock * p%NPU
      read( nf,* ) p%AU(i)
    enddo
    !---- D
    do i=1,nblock * p%NP
      read( nf, * ) p%D(i)
    enddo
    !---- B
    do i=1,p%NDOF * p%NP
      read( nf, * ) p%B(i)
    enddo
 
    close( nf )
  end subroutine
 
  subroutine hecmw_clone_matrix(hecMATorig,hecMAT)
    type (hecmwST_matrix    ) :: hecMATorig
    type (hecmwST_matrix    ),pointer :: hecMAT
    allocate(hecmat)
    call hecmw_nullify_matrix( hecmat )
 
    hecmat%B => hecmatorig%B
    hecmat%X => hecmatorig%X
    hecmat%D => hecmatorig%D
    hecmat%AL => hecmatorig%AL
    hecmat%AU => hecmatorig%AU
    hecmat%indexL => hecmatorig%indexL
    hecmat%indexU => hecmatorig%indexU
    hecmat%itemL  => hecmatorig%itemL
    hecmat%itemU  => hecmatorig%itemU
    hecmat%N    = hecmatorig%N
    hecmat%NP   = hecmatorig%NP
    hecmat%NPL  = hecmatorig%NPL
    hecmat%NPU  = hecmatorig%NPU
    hecmat%NDOF = hecmatorig%NDOF
    hecmat%Iarray = hecmatorig%Iarray
    hecmat%Rarray = hecmatorig%Rarray
  end subroutine hecmw_clone_matrix
 
  subroutine hecmw_copy_matrix(hecMATorig,hecMAT)
    type (hecmwST_matrix    ) :: hecMATorig
    type (hecmwST_matrix    ),pointer :: hecMAT
    integer(kind=kint) NDOF,NDOF2,N,NP,NPL,NPU
 
    allocate(hecmat)
    call hecmw_nullify_matrix( hecmat )
    n = hecmatorig%N
    ndof = hecmatorig%NDOF
    ndof2 = ndof*ndof
    np = hecmatorig%NP
    npl = hecmatorig%NPL
    npu = hecmatorig%NPU
    hecmat%N    = n
    hecmat%NP   = np
    hecmat%NPL  = npl
    hecmat%NPU  = npu
    hecmat%NDOF = ndof
    allocate(hecmat%B(ndof*np))
    allocate(hecmat%X(ndof*np))
    allocate(hecmat%D(ndof2*np))
    allocate(hecmat%AU(ndof2*npu))
    allocate(hecmat%AL(ndof2*npl))
    allocate(hecmat%indexL(0:np), hecmat%indexU(0:np), hecmat%itemL(npl), hecmat%itemU(npu))
    hecmat%B       = hecmatorig%B
    hecmat%X       = hecmatorig%X
    hecmat%D       = hecmatorig%D
    hecmat%AU      = hecmatorig%AU
    hecmat%AL      = hecmatorig%AL
    hecmat%indexL  = hecmatorig%indexL
    hecmat%indexU  = hecmatorig%indexU
    hecmat%itemL   = hecmatorig%itemL
    hecmat%itemU   = hecmatorig%itemU
    hecmat%Iarray  = hecmatorig%Iarray
    hecmat%Rarray  = hecmatorig%Rarray
  end subroutine hecmw_copy_matrix
 
  subroutine hecmw_blockmatrix_expand(hecMATorig,hecMAT,NDOF)
    type (hecmwST_matrix    ) :: hecMATorig
    type (hecmwST_matrix    ),pointer :: hecMAT
    integer(kind=kint) NDOF,NDOF2,oNDOF,oNDOF2,i,j,k
 
    ndof2 = ndof*ndof
    ondof = hecmatorig%NDOF
    ondof2 = ondof*ondof
    allocate(hecmat)
    call hecmw_nullify_matrix( hecmat )
 
    allocate(hecmat%B(ndof*hecmatorig%NP))
    allocate(hecmat%X(ndof*hecmatorig%NP))
    allocate(hecmat%D(ndof2*hecmatorig%NP))
    allocate(hecmat%AL(ndof2*hecmatorig%NPL))
    allocate(hecmat%AU(ndof2*hecmatorig%NPU))
    hecmat%indexL => hecmatorig%indexL
    hecmat%indexU => hecmatorig%indexU
    hecmat%itemL  => hecmatorig%itemL
    hecmat%itemU  => hecmatorig%itemU
    hecmat%N    = hecmatorig%N
    hecmat%NP   = hecmatorig%NP
    hecmat%NPL  = hecmatorig%NPL
    hecmat%NPU  = hecmatorig%NPU
    hecmat%NDOF = ndof
    hecmat%Iarray = hecmatorig%Iarray
    hecmat%Rarray = hecmatorig%Rarray
    hecmat%X = 0.0d0
    do i = 1, hecmatorig%NP
      do j = 1, ndof
        do k = 1, ndof
          if (j<=ondof .and. k<=ondof) then
            hecmat%D(ndof2*(i-1) + (j-1)*ndof + k) = hecmatorig%D(ondof2*(i-1) + (j-1)*ondof + k)
          else
            if(j==k) then
              hecmat%D(ndof2*(i-1) + (j-1)*ndof + k)=1
            else
              hecmat%D(ndof2*(i-1) + (j-1)*ndof + k)=0
            end if
          end if
        end do
        if (j<=ondof) then
          hecmat%B(ndof*(i-1) + j) = hecmatorig%B(ondof*(i-1) + j)
        else
          hecmat%B(ndof*(i-1) + j)=0
        end if
      end do
    end do
    do i = 1, hecmatorig%NPL
      do j = 1, ndof
        do k = 1, ndof
          if (j<=ondof .and. k<=ondof) then
            hecmat%AL(ndof2*(i-1) + (j-1)*ndof + k) = hecmatorig%AL(ondof2*(i-1) + (j-1)*ondof + k)
          else
            hecmat%AL(ndof2*(i-1) + (j-1)*ndof + k) = 0
          end if
        end do
      end do
    end do
    do i = 1, hecmatorig%NPU
      do j = 1, ndof
        do k = 1, ndof
          if (j<=ondof .and. k<=ondof) then
            hecmat%AU(ndof2*(i-1) + (j-1)*ndof + k) = hecmatorig%AU(ondof2*(i-1) + (j-1)*ondof + k)
          else
            hecmat%AU(ndof2*(i-1) + (j-1)*ndof + k) = 0
          end if
        end do
      end do
    end do
  end subroutine hecmw_blockmatrix_expand
  subroutine hecmw_vector_contract(hecMATorig,hecMAT,NDOF)
    type (hecmwST_matrix    ) :: hecMATorig
    type (hecmwST_matrix    ),pointer :: hecMAT
    integer(kind=kint) NDOF,NDOF2,oNDOF,i,j
    ndof2 = ndof*ndof
    ondof = hecmatorig%NDOF
    do i = 1, hecmatorig%NP
      do j = 1, ondof
        hecmatorig%X(ondof*(i-1) + j) = hecmat%X(ndof*(i-1) + j)
      end do
    end do
    hecmatorig%Iarray = hecmat%Iarray
    hecmatorig%Rarray = hecmat%Rarray
    deallocate(hecmat%B)
    deallocate(hecmat%D)
    deallocate(hecmat%X)
    deallocate(hecmat%AL)
    deallocate(hecmat%AU)
    deallocate(hecmat)
  end subroutine hecmw_vector_contract
end module hecmw_util