hecmw_ebc_defer.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_ebc_defer
  use hecmw_util
  use m_hecmw_comm_f
  implicit none
 
  private
  public :: hecmwst_ebc
  public :: hecmw_ebc_init
  public :: hecmw_ebc_set
  public :: hecmw_ebc_apply
  public :: hecmw_ebc_finalize
 
  type hecmwst_ebc
    integer(kind=kint) :: ndof = 0
    integer(kind=kint), pointer :: mark(:) => null()
    real(kind=kreal), pointer :: val(:) => null()
  end type hecmwst_ebc
 
contains
 
  !C
  !C***
  !C*** hecmw_ebc_init
  !C***
  !C
  subroutine hecmw_ebc_init(hecMAT, hecEBC)
    implicit none
    type (hecmwst_matrix), intent(in) :: hecmat
    type (hecmwst_ebc), intent(inout) :: hecebc
    integer(kind=kint) :: npndof
 
    call hecmw_ebc_finalize(hecebc)
    hecebc%ndof = hecmat%NDOF
    npndof = hecmat%NP * hecmat%NDOF
    allocate(hecebc%mark(npndof))
    allocate(hecebc%val(npndof))
    hecebc%mark(:) = 0
    hecebc%val(:) = 0.d0
  end subroutine hecmw_ebc_init
 
  !C
  !C***
  !C*** hecmw_ebc_set
  !C***
  !C
  subroutine hecmw_ebc_set(hecEBC, inode, idof, val)
    implicit none
    type (hecmwst_ebc), intent(inout) :: hecebc
    integer(kind=kint), intent(in) :: inode, idof
    real(kind=kreal), intent(in) :: val
    integer(kind=kint) :: k
 
    if (idof > hecebc%ndof) return
    k = hecebc%ndof * (inode - 1) + idof
    if (hecebc%mark(k) /= 0 .and. hecebc%val(k) /= val) then
      write(*,'(a,i0,a,i0,a,i0,a,2(1pe14.6))') 'WARNING: rank ', hecmw_comm_get_rank(), &
          ': boundary value of local node ', inode, ' dof ', idof, ' overwritten:', hecebc%val(k), val
    endif
    hecebc%mark(k) = 1
    hecebc%val(k) = val
  end subroutine hecmw_ebc_set
 
  !C
  !C***
  !C*** hecmw_ebc_apply
  !C***
  !C
  subroutine hecmw_ebc_apply(hecMESH, hecMAT, hecEBC, conMAT)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_matrix), intent(inout) :: hecmat
    type (hecmwst_ebc), intent(inout) :: hecebc
    type (hecmwst_matrix), intent(inout), optional :: conmat
 
    call hecmw_ebc_extend(hecmesh, hecmat, hecebc)
    call hecmw_ebc_convert_slave(hecmesh, hecmat, hecebc)
    call hecmw_ebc_impose(hecmat, hecebc, 1.d0)
    if (present(conmat)) call hecmw_ebc_impose(conmat, hecebc, 0.d0)
  end subroutine hecmw_ebc_apply
 
  !C
  !C***
  !C*** hecmw_ebc_extend
  !C***
  !C
  subroutine hecmw_ebc_extend(hecMESH, hecMAT, hecEBC)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_matrix), intent(in) :: hecmat
    type (hecmwst_ebc), intent(inout) :: hecebc
    integer(kind=kint), pointer :: mark(:)
    real(kind=kreal), pointer :: val(:)
    integer(kind=kint) :: ndof, npndof, npndof_old, i
 
    ndof = hecmat%NDOF
    npndof = hecmat%NP * ndof
    npndof_old = size(hecebc%mark)
 
    if (npndof > npndof_old) then
      allocate(mark(npndof))
      allocate(val(npndof))
      mark(:) = 0
      val(:) = 0.d0
      do i = 1, npndof_old
        mark(i) = hecebc%mark(i)
        val(i) = hecebc%val(i)
      enddo
      deallocate(hecebc%mark)
      deallocate(hecebc%val)
      hecebc%mark => mark
      hecebc%val => val
    endif
 
    call hecmw_assemble_i(hecmesh, hecebc%mark, hecmat%NP, ndof)
    call hecmw_assemble_r(hecmesh, hecebc%val, hecmat%NP, ndof)
    call hecmw_update_i(hecmesh, hecebc%mark, hecmat%NP, ndof)
    call hecmw_update_r(hecmesh, hecebc%val, hecmat%NP, ndof)
 
    do i = 1, npndof
      if (hecebc%mark(i) == 0) cycle
      hecebc%val(i) = hecebc%val(i) / hecebc%mark(i)
      hecebc%mark(i) = 1
    enddo
  end subroutine hecmw_ebc_extend
 
  !C
  !C***
  !C*** hecmw_ebc_convert_slave
  !C***
  !C
  subroutine hecmw_ebc_convert_slave(hecMESH, hecMAT, hecEBC)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_matrix), intent(in) :: hecmat
    type (hecmwst_ebc), intent(inout) :: hecebc
    integer(kind=kint), allocatable :: dmark(:)
    real(kind=kreal), allocatable :: dval(:)
    integer(kind=kint) :: ndof, npndof, i, j, k, km, ks, kk, nchange
    real(kind=kreal) :: um, tol
 
    ndof = hecmat%NDOF
    npndof = hecmat%NP * ndof
    allocate(dmark(npndof))
    allocate(dval(npndof))
    dmark(:) = 0
    dval(:) = 0.d0
    nchange = 0
 
    outer: do i = 1, hecmesh%mpc%n_mpc
      do j = hecmesh%mpc%mpc_index(i-1)+1, hecmesh%mpc%mpc_index(i)
        if (hecmesh%mpc%mpc_dof(j) > ndof) cycle outer
      enddo
      k = hecmesh%mpc%mpc_index(i-1) + 1
      ks = ndof * (hecmesh%mpc%mpc_item(k) - 1) + hecmesh%mpc%mpc_dof(k)
      if (hecebc%mark(ks) == 0) cycle
      if (hecmesh%mpc%mpc_index(i) - hecmesh%mpc%mpc_index(i-1) /= 2) then
        write(*,'(a,i0,a,i0,a)') 'ERROR: a boundary condition is given to node ', &
            hecmesh%global_node_ID(hecmesh%mpc%mpc_item(k)), ' dof ', hecmesh%mpc%mpc_dof(k), &
            ', the slave of an !EQUATION with more than one master; it cannot be imposed as a boundary condition'
        call hecmw_abort(hecmw_comm_get_comm())
      endif
      km = k + 1
      kk = ndof * (hecmesh%mpc%mpc_item(km) - 1) + hecmesh%mpc%mpc_dof(km)
      um = (hecmesh%mpc%mpc_const(i) - hecmesh%mpc%mpc_val(k) * hecebc%val(ks)) / hecmesh%mpc%mpc_val(km)
      tol = 1.d-10 * max(1.d0, abs(um))
      if ((hecebc%mark(kk) /= 0 .and. abs(hecebc%val(kk) - um) > tol) .or. &
          (dmark(kk) /= 0 .and. abs(dval(kk) - um) > tol)) then
        write(*,'(a,i0,a,i0,a,i0,a,i0,a)') 'ERROR: the boundary condition on node ', &
            hecmesh%global_node_ID(hecmesh%mpc%mpc_item(k)), ' dof ', hecmesh%mpc%mpc_dof(k), &
            ', the slave of an !EQUATION, conflicts with the one on its master node ', &
            hecmesh%global_node_ID(hecmesh%mpc%mpc_item(km)), ' dof ', hecmesh%mpc%mpc_dof(km), ''
        call hecmw_abort(hecmw_comm_get_comm())
      endif
      if (hecebc%mark(kk) == 0) then
        dmark(kk) = 1
        dval(kk) = um
      endif
      dmark(ks) = -1
      nchange = nchange + 1
    enddo outer
 
    call hecmw_allreduce_i1(hecmesh, nchange, hecmw_sum)
    if (nchange > 0) then
      call hecmw_assemble_i(hecmesh, dmark, hecmat%NP, ndof)
      call hecmw_assemble_r(hecmesh, dval, hecmat%NP, ndof)
      do i = 1, npndof
        if (dmark(i) > 0) then
          hecebc%mark(i) = 1
          hecebc%val(i) = dval(i) / dmark(i)
        else if (dmark(i) < 0) then
          hecebc%mark(i) = 0
          hecebc%val(i) = 0.d0
        endif
      enddo
      call hecmw_update_i(hecmesh, hecebc%mark, hecmat%NP, ndof)
      call hecmw_update_r(hecmesh, hecebc%val, hecmat%NP, ndof)
    endif
 
    deallocate(dmark)
    deallocate(dval)
  end subroutine hecmw_ebc_convert_slave
 
  !C
  !C***
  !C*** hecmw_ebc_impose
  !C***
  !C
  subroutine hecmw_ebc_impose(hecMAT, hecEBC, diag)
    implicit none
    type (hecmwst_matrix), intent(inout) :: hecmat
    type (hecmwst_ebc), intent(in) :: hecebc
    real(kind=kreal), intent(in) :: diag
    integer(kind=kint) :: ndof, ndof2, i, j, k, idof, jdof, ir, jc, idx
 
    ndof = hecmat%NDOF
    ndof2 = ndof * ndof
 
    do i = 1, hecmat%NP
      do idof = 1, ndof
        ir = ndof * (i - 1) + idof
        do jdof = 1, ndof
          jc = ndof * (i - 1) + jdof
          idx = ndof2 * (i - 1) + ndof * (idof - 1) + jdof
          if (hecebc%mark(ir) /= 0) then
            hecmat%D(idx) = 0.d0
          else if (hecebc%mark(jc) /= 0) then
            hecmat%B(ir) = hecmat%B(ir) - hecmat%D(idx) * hecebc%val(jc)
            hecmat%D(idx) = 0.d0
          endif
        enddo
      enddo
      do k = hecmat%indexL(i-1) + 1, hecmat%indexL(i)
        j = hecmat%itemL(k)
        do idof = 1, ndof
          ir = ndof * (i - 1) + idof
          do jdof = 1, ndof
            jc = ndof * (j - 1) + jdof
            idx = ndof2 * (k - 1) + ndof * (idof - 1) + jdof
            if (hecebc%mark(ir) /= 0) then
              hecmat%AL(idx) = 0.d0
            else if (hecebc%mark(jc) /= 0) then
              hecmat%B(ir) = hecmat%B(ir) - hecmat%AL(idx) * hecebc%val(jc)
              hecmat%AL(idx) = 0.d0
            endif
          enddo
        enddo
      enddo
      do k = hecmat%indexU(i-1) + 1, hecmat%indexU(i)
        j = hecmat%itemU(k)
        do idof = 1, ndof
          ir = ndof * (i - 1) + idof
          do jdof = 1, ndof
            jc = ndof * (j - 1) + jdof
            idx = ndof2 * (k - 1) + ndof * (idof - 1) + jdof
            if (hecebc%mark(ir) /= 0) then
              hecmat%AU(idx) = 0.d0
            else if (hecebc%mark(jc) /= 0) then
              hecmat%B(ir) = hecmat%B(ir) - hecmat%AU(idx) * hecebc%val(jc)
              hecmat%AU(idx) = 0.d0
            endif
          enddo
        enddo
      enddo
    enddo
 
    do i = 1, hecmat%NP
      do idof = 1, ndof
        ir = ndof * (i - 1) + idof
        if (hecebc%mark(ir) == 0) cycle
        hecmat%D(ndof2 * (i - 1) + ndof * (idof - 1) + idof) = diag
        hecmat%B(ir) = diag * hecebc%val(ir)
      enddo
    enddo
  end subroutine hecmw_ebc_impose
 
  !C
  !C***
  !C*** hecmw_ebc_finalize
  !C***
  !C
  subroutine hecmw_ebc_finalize(hecEBC)
    implicit none
    type (hecmwst_ebc), intent(inout) :: hecebc
 
    hecebc%ndof = 0
    if (associated(hecebc%mark)) deallocate(hecebc%mark)
    if (associated(hecebc%val)) deallocate(hecebc%val)
    nullify(hecebc%mark)
    nullify(hecebc%val)
  end subroutine hecmw_ebc_finalize
 
end module hecmw_ebc_defer