hecmw_precond_rbm.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2026 FrontISTR Commons
! This software is released under the MIT License, see License.txt
!-------------------------------------------------------------------------------
module hecmw_precond_rbm
  use hecmw_util
  use hecmw_etype
  use m_hecmw_comm_f, only: hecmw_allreduce_r
  implicit none
 
  private
  public :: hecmw_precond_rbm_from_mesh
 
contains
 
  subroutine hecmw_precond_rbm_from_mesh(hecMESH, ndof, rbm)
    implicit none
    type(hecmwst_local_mesh), intent(in)  :: hecmesh
    integer(kind=kint),       intent(in)  :: ndof
    real(kind=kreal),         intent(out) :: rbm(*)   
    integer(kind=kint), allocatable :: mark(:)
    integer(kind=kint) :: nint, vl, nmode, node, base
    integer(kind=kint) :: itype, ic_type, nn, is, ie, icel, iis, j, nod
    real(kind=kreal) :: cx, cy, cz, x, y, z, csum(4)
 
    nint = hecmesh%nn_internal
    vl   = nint * ndof
    select case (ndof)
    case (1);    nmode = 1
    case (2);    nmode = 3
    case default; nmode = 6
    end select
 
    ! mark 33struct shell rotation nodes (latter-half element nodes)
    allocate(mark(hecmesh%n_node)); mark = 0
    if (ndof == 3) then
      do itype = 1, hecmesh%n_elem_type
        ic_type = hecmesh%elem_type_item(itype)
        if (hecmw_is_etype_33struct(ic_type)) then
          nn = hecmw_get_max_node(ic_type)
          is = hecmesh%elem_type_index(itype-1)+1
          ie = hecmesh%elem_type_index(itype)
          do icel = is, ie
            iis = hecmesh%elem_node_index(icel-1)
            do j = nn/2+1, nn
              nod = hecmesh%elem_node_item(iis+j); mark(nod) = 1
            end do
          end do
        end if
      end do
    end if
 
    ! global centroid, for well-conditioned rotation modes far from the origin
    ! (sum of coordinates and node count reduced together in one allreduce)
    csum = 0.0d0
    do node = 1, nint
      csum(1) = csum(1) + hecmesh%node(3*node-2)
      csum(2) = csum(2) + hecmesh%node(3*node-1)
      csum(3) = csum(3) + hecmesh%node(3*node)
    end do
    csum(4) = real(nint, kreal)
    call hecmw_allreduce_r(hecmesh, csum, 4, hecmw_sum)
    cx = 0.0d0; cy = 0.0d0; cz = 0.0d0
    if (csum(4) > 0.0d0) then
      cx = csum(1)/csum(4); cy = csum(2)/csum(4); cz = csum(3)/csum(4)
    end if
 
    rbm(1:vl*nmode) = 0.0d0
    if (ndof == 1) then
      do node = 1, nint
        rbm(node) = 1.0d0                                   ! constant
      end do
    else if (ndof == 2) then
      do node = 1, nint
        x = hecmesh%node(3*node-2) - cx; y = hecmesh%node(3*node-1) - cy
        base = 2*(node-1)
        rbm(        base+1) = 1.0d0                          ! trans-x
        rbm(   vl + base+2) = 1.0d0                          ! trans-y
        rbm( 2*vl + base+1) = -y                            ! rot-z
        rbm( 2*vl + base+2) =  x
      end do
    else   ! ndof = 3 or 6
      do node = 1, nint
        base = ndof*(node-1)
        if (ndof == 3 .and. mark(node) == 1) then
          ! 33struct rotation node: rotational rigid-body modes only [0 | I3]
          rbm( 3*vl + base+1) = 1.0d0
          rbm( 4*vl + base+2) = 1.0d0
          rbm( 5*vl + base+3) = 1.0d0
        else
          x = hecmesh%node(3*node-2) - cx
          y = hecmesh%node(3*node-1) - cy
          z = hecmesh%node(3*node)   - cz
          rbm(        base+1) = 1.0d0                        ! trans-x
          rbm(   vl + base+2) = 1.0d0                        ! trans-y
          rbm( 2*vl + base+3) = 1.0d0                        ! trans-z
          rbm( 3*vl + base+2) = -z                          ! rot-x
          rbm( 3*vl + base+3) =  y
          rbm( 4*vl + base+1) =  z                          ! rot-y
          rbm( 4*vl + base+3) = -x
          rbm( 5*vl + base+1) = -y                          ! rot-z
          rbm( 5*vl + base+2) =  x
          if (ndof == 6) then
            rbm( 3*vl + base+4) = 1.0d0                      ! shell rotational DOF identity
            rbm( 4*vl + base+5) = 1.0d0
            rbm( 5*vl + base+6) = 1.0d0
          end if
        end if
      end do
    end if
 
    deallocate(mark)
  end subroutine hecmw_precond_rbm_from_mesh
 
end module hecmw_precond_rbm