precheck_mesh_quality.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 m_precheck_mesh_quality
 
  use hecmw
  use m_fstr
  use m_precheck_lib_elements
 
  implicit none
 
  private
  public :: precheck_mesh_quality
 
contains
 
  subroutine precheck_get_thickness(hecMESH, mid, thick)
    type(hecmwST_local_mesh) :: hecMESH
    integer(kind=kint) :: mid, ihead
    real(kind=kreal)   :: thick
 
    ihead = hecmesh%section%sect_R_index(mid-1)
    thick = hecmesh%section%sect_R_item(ihead+1)
  end subroutine precheck_get_thickness
 
  subroutine precheck_mesh_quality(hecMESH, hecMAT, elem_vol, elem_asp)
    implicit none
 
    type(hecmwst_matrix)     :: hecmat
    type(hecmwst_local_mesh) :: hecmesh
    real(kind=kreal), intent(out) :: elem_vol(:)
    real(kind=kreal), intent(out) :: elem_asp(:)
 
    !** Parameters
    real(kind=kreal), parameter :: aspect_warn_threshold = 50.0d0
 
    !** Local variables
    integer(kind=kint) :: nelem, mid, j, isect, icel, iis
    integer(kind=kint) :: ndof2, nelem_wo_mpc
    integer(kind=kint) :: ie, ia, jelem, ic_type, nn, js, je, itype
    integer(kind=kint) :: nodlocal(20), ntotsum(1)
    integer(kind=kint) :: nonzero
    integer(kind=kint) :: ntotbuf(2)
    !  Element type counts: 3-digit codes only (<=999). 4-digit codes (2322, 3414, 3422, 1031-1042) are internal/deprecated
    !  variants or patch elements not processed by precheck; warn and skip if encountered.
    integer(kind=kint), parameter :: etype_code_max = 999
    integer(kind=kint) :: etype_count(etype_code_max), i_etype
    real(kind=kreal)   :: ntdof2
    real(kind=kreal)   :: al, almin, almax, aa, thick, vol, avvol
    real(kind=kreal)   :: tvol, tvmax, tvmin, tlmax, tlmin, asp, aspmax
    real(kind=kreal)   :: xx(20), yy(20), zz(20)
    real(kind=kreal)   :: totsum(1), totmax(3), totmin(2)
 
    !C INIT
    elem_vol = 0.0d0
    elem_asp = 0.0d0
    nelem  = 0
    tvol   = 0.0
    tvmax  = 0.0
    tvmin  = 1.0e+20
    tlmax  = 0.0
    tlmin  = 1.0e+20
    aspmax = 0.0
 
    !C Mesh Summary (local values to ILOG, global values to stdout via rank 0)
    ndof2  = hecmesh%n_dof**2
 
    !C Local summary to log file
    write(ilog,"(a)") '###  Mesh Summary (local)  ###'
    write(ilog,"(a,i12)") '  Num of node:',hecmesh%n_node
    write(ilog,"(a,i12)") '  Num of DOF :',hecmesh%n_node*hecmesh%n_dof
    write(ilog,"(a,i12)") '  Num of elem:',hecmesh%n_elem
    nelem_wo_mpc = 0
    do itype = 1, hecmesh%n_elem_type
      js = hecmesh%elem_type_index(itype-1)
      je = hecmesh%elem_type_index(itype  )
      ic_type = hecmesh%elem_type_item(itype)
      if(.not. (hecmw_is_etype_link(ic_type) .or. hecmw_is_etype_patch(ic_type))) &
           nelem_wo_mpc = nelem_wo_mpc + je-js
    enddo
    write(ilog,"(a,i12)") '   ** w/o MPC/Patch:',nelem_wo_mpc
    do itype = 1, hecmesh%n_elem_type
      js = hecmesh%elem_type_index(itype-1)
      je = hecmesh%elem_type_index(itype  )
      ic_type = hecmesh%elem_type_item(itype)
      write(ilog,"(a,i4,a,i12)") '  Num of ',ic_type,':',je-js
    enddo
    nonzero = ndof2*(hecmat%NP + hecmat%NPU + hecmat%NPL)
    write(ilog,"(a,i12)") '  Num of NZ  :',nonzero
    ntdof2 = dble(hecmesh%n_node*hecmesh%n_dof)**2
    write(ilog,"(a,1pe12.5,a)") '  Sparsity   :',100.0d0*dble(nonzero)/dble(ntdof2),"[%]"
 
    !C Global summary to stdout (rank 0 only)
    ntotsum(1) = hecmesh%nn_internal
    call hecmw_allreduce_i(hecmesh, ntotsum, 1, hecmw_sum)
    ntotbuf(1) = hecmesh%n_elem
    ntotbuf(2) = nelem_wo_mpc
    call hecmw_allreduce_i(hecmesh, ntotbuf, 2, hecmw_sum)
    etype_count(:) = 0
    do itype = 1, hecmesh%n_elem_type
      ic_type = hecmesh%elem_type_item(itype)
      if( ic_type .ge. 1 .and. ic_type .le. etype_code_max ) then
        etype_count(ic_type) = etype_count(ic_type) + hecmesh%elem_type_index(itype) - hecmesh%elem_type_index(itype-1)
      else
        write(ilog,"(a,i6,a)") '  %%%  WARNING %%% Element type ', ic_type, ' is out of range for Mesh Summary count'
      endif
    enddo
    call hecmw_allreduce_i(hecmesh, etype_count, etype_code_max, hecmw_sum)
    if( hecmesh%my_rank .eq. 0 ) then
      write(*,"(a)") '###  Mesh Summary  ###'
      write(*,"(a,i12)") '  Num of node       :',ntotsum(1)
      write(*,"(a,i12)") '  Num of DOF        :',ntotsum(1)*hecmesh%n_dof
      write(*,"(a,i12)") '  Num of elem       :',ntotbuf(1)
      write(*,"(a,i12)") '   ** w/o MPC/Patch :',ntotbuf(2)
      do i_etype = 1, etype_code_max
        if( etype_count(i_etype) .gt. 0 ) then
          write(*,"(a,i5,a,i12)") '  Num of ',i_etype,':',etype_count(i_etype)
        endif
      enddo
    endif
 
    !C 3D
    if( hecmesh%n_dof .eq. 3 ) then
      do ie = 1, hecmesh%n_elem
        ia = hecmesh%elem_ID(ie*2)
        if( ia.ne.hecmesh%my_rank ) cycle
        jelem = hecmesh%global_elem_ID(ie)
        ic_type = hecmesh%elem_type(ie)
        if (.not. (hecmw_is_etype_rod(ic_type) .or. hecmw_is_etype_solid(ic_type) &
            & .or. hecmw_is_etype_beam(ic_type) .or. hecmw_is_etype_shell(ic_type))) then
          write(ilog,*) jelem, ' This Element cannot be checked. Type=',ic_type
          cycle
        endif
        vol = 0.0d0
        almax = 0.0d0
        almin = 0.0d0
        nn = hecmw_get_max_node(ic_type)
        js = hecmesh%elem_node_index(ie-1)
        je = hecmesh%elem_node_index(ie)
        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)  )
        enddo
 
        if    ( ic_type.eq.111 ) then
          isect = hecmesh%section_ID(ie)
          mid = hecmesh%section%sect_mat_ID_item(isect)
          call precheck_get_thickness( hecmesh,mid,aa )
          al = sqrt( (xx(2)-xx(1))**2+(yy(2)-yy(1))**2+(zz(2)-zz(1))**2 )
          vol = aa*al
          almax = al
          almin = al
        elseif( hecmw_is_etype_solid(ic_type) ) then
          call precheck_calc_vol_asp(ic_type, nn, xx, yy, zz, 0.0d0, vol, almax, almin)
        elseif( ic_type.eq.641 ) then
          vol = 1.0d-12
        elseif( ic_type.eq.761 .or. ic_type.eq.781 ) then
          ! 2nd-order shell: volume/edge calculation not yet implemented in elementInfo;
          ! use dummy value to avoid zero-volume error
          vol = 1.0d-12
        endif
 
        if( vol.le.0.0 ) then
          write(ilog,*) '  %%%  ERROR %%%  Volume of Element no.=',jelem,' is zero or negative.'
        endif
        nelem = nelem + 1
        elem_vol(ie) = vol
        tvol = tvol + vol
        if( vol.gt.tvmax ) tvmax = vol
        if( vol.lt.tvmin ) tvmin = vol
        if( almax.gt.tlmax ) tlmax = almax
        if( almin.lt.tlmin ) tlmin = almin
        asp = almax/almin
        elem_asp(ie) = asp
        if( asp.gt.aspmax ) aspmax = asp
        if( asp.gt.aspect_warn_threshold ) then
          write(ilog,*) '  %%%  WARNING %%% Aspect ratio of Element no.=',jelem, &
               &        ' exceeds ',aspect_warn_threshold
          write(ilog,*) '      Maximum length =',almax
          write(ilog,*) '      Minimum length =',almin
        endif
      enddo
 
    !C 2D
    elseif( hecmesh%n_dof .eq. 2 ) then
      do ie = 1, hecmesh%n_elem
        ia = hecmesh%elem_ID(ie*2)
        if( ia.ne.hecmesh%my_rank ) cycle
        jelem = hecmesh%global_elem_ID(ie)
        ic_type = hecmesh%elem_type(ie)
        if (.not. (hecmw_is_etype_rod(ic_type) .or. hecmw_is_etype_surface(ic_type))) then
          write(ilog,*) jelem, ' This Element cannot be checked. Type=',ic_type
          cycle
        endif
        vol = 0.0d0
        almax = 0.0d0
        almin = 0.0d0
        nn = hecmw_get_max_node(ic_type)
        js = hecmesh%elem_node_index(ie-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)
        enddo
 
        isect = hecmesh%section_ID(ie)
        mid = hecmesh%section%sect_mat_ID_item(isect)
        call precheck_get_thickness( hecmesh,mid,aa )
 
        if    ( ic_type.eq.111 ) then
          al = sqrt( (xx(2)-xx(1))**2+(yy(2)-yy(1))**2 )
          vol = aa*al
          if( al.gt.tlmax ) tlmax = al
          if( al.lt.tlmin ) tlmin = al
          aspmax = 1.0
        elseif( hecmw_is_etype_surface(ic_type) ) then
          call precheck_calc_vol_asp(ic_type, nn, xx, yy, zz, aa, vol, almax, almin)
        else
          vol = 0.0
        endif
 
        if( vol.le.0.0 ) then
          write(ilog,*) '  %%%  ERROR %%%  Volume of Element no.=',jelem,' is zero or negative.'
        endif
        nelem = nelem + 1
        elem_vol(ie) = vol
        tvol = tvol + vol
        if( vol.gt.tvmax ) tvmax = vol
        if( vol.lt.tvmin ) tvmin = vol
        if( almax.gt.tlmax ) tlmax = almax
        if( almin.lt.tlmin ) tlmin = almin
        asp = almax/almin
        elem_asp(ie) = asp
        if( asp.gt.aspmax ) aspmax = asp
        if( asp.gt.aspect_warn_threshold ) then
          write(ilog,*) '  %%%  WARNING %%% Aspect ratio of Element no.=',jelem, &
               &        ' exceeds ',aspect_warn_threshold
          write(ilog,*) '      Maximum length =',almax
          write(ilog,*) '      Minimum length =',almin
        endif
      enddo
 
    !C SHELL
    elseif( hecmesh%n_dof .eq. 6 ) then
      do ie = 1, hecmesh%n_elem
        ia = hecmesh%elem_ID(ie*2)
        if( ia.ne.hecmesh%my_rank ) cycle
        jelem = hecmesh%global_elem_ID(ie)
        ic_type = hecmesh%elem_type(ie)
        if (.not. (hecmw_is_etype_beam(ic_type) .or. hecmw_is_etype_shell(ic_type))) then
          write(ilog,*) jelem, ' This Element cannot be checked. Type=',ic_type
          cycle
        endif
        vol = 0.0d0
        almax = 0.0d0
        almin = 0.0d0
        nn = hecmw_get_max_node(ic_type)
        js = hecmesh%elem_node_index(ie-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)  )
        enddo
 
        isect = hecmesh%section_ID(ie)
        mid = hecmesh%section%sect_mat_ID_item(isect)
        call precheck_get_thickness( hecmesh,mid,aa )
 
        if    ( ic_type.eq.111 ) then
          al = sqrt( (xx(2)-xx(1))**2+(yy(2)-yy(1))**2+(zz(2)-zz(1))**2 )
          vol = aa*al
          if( al.gt.tlmax ) tlmax = al
          if( al.lt.tlmin ) tlmin = al
          aspmax = 1.0
        elseif( ic_type.eq.731 .or. ic_type.eq.741 ) then
          call precheck_calc_vol_asp(ic_type, nn, xx, yy, zz, aa, vol, almax, almin)
        endif
 
        if( vol.le.0.0 ) then
          write(ilog,*) '  %%%  ERROR %%%  Volume of Element no.=',jelem,' is zero or negative.'
        endif
        nelem = nelem + 1
        elem_vol(ie) = vol
        tvol = tvol + vol
        if( vol.gt.tvmax ) tvmax = vol
        if( vol.lt.tvmin ) tvmin = vol
        if( almax.gt.tlmax ) tlmax = almax
        if( almin.lt.tlmin ) tlmin = almin
        asp = almax/almin
        elem_asp(ie) = asp
        if( asp.gt.aspmax ) aspmax = asp
        if( asp.gt.aspect_warn_threshold ) then
          write(ilog,*) '  %%%  WARNING %%% Aspect ratio of Element no.=',jelem, &
               &        ' exceeds ',aspect_warn_threshold
          write(ilog,*) '      Maximum length =',almax
          write(ilog,*) '      Minimum length =',almin
        endif
      enddo
    endif
 
    !C Local Summary (per-rank, to ILOG only)
    avvol = tvol / nelem
    write(ilog,*) '###  Local Summary  ###'
    write(ilog,*) '  Total volume                = ',tvol
    write(ilog,*) '  Average volume of elements  = ',avvol
    write(ilog,*) '  Maximum volume of elements  = ',tvmax
    write(ilog,*) '  Minimum volume of elements  = ',tvmin
    write(ilog,*) '  Maximum edge length         = ',tlmax
    write(ilog,*) '  Minimum edge length         = ',tlmin
    write(ilog,*) '  Maximum aspect ratio        = ',aspmax
 
    !C Global Summary (allreduce, rank 0 outputs to both ILOG and stdout)
    totsum(1) = tvol
    call hecmw_allreduce_r(hecmesh,totsum,1,hecmw_sum)
    ntotsum(1) = nelem
    call hecmw_allreduce_i(hecmesh,ntotsum,1,hecmw_sum)
    totmax(1) = tvmax
    totmax(2) = tlmax
    totmax(3) = aspmax
    call hecmw_allreduce_r(hecmesh,totmax,3,hecmw_max)
    totmin(1) = tvmin
    totmin(2) = tlmin
    call hecmw_allreduce_r(hecmesh,totmin,2,hecmw_min)
    if( hecmesh%my_rank .eq. 0 ) then
      avvol = totsum(1) / ntotsum(1)
      write(ilog,*) '###  Global Summary  ###'
      write(ilog,*) '  Total volume                = ',totsum(1)
      write(*,*)    '  Total volume                = ',totsum(1)
      write(ilog,*) '  Average volume of elements  = ',avvol
      write(*,*)    '  Average volume of elements  = ',avvol
      write(ilog,*) '  Maximum volume of elements  = ',totmax(1)
      write(*,*)    '  Maximum volume of elements  = ',totmax(1)
      write(ilog,*) '  Minimum volume of elements  = ',totmin(1)
      write(*,*)    '  Minimum volume of elements  = ',totmin(1)
      write(ilog,*) '  Maximum edge length         = ',totmax(2)
      write(*,*)    '  Maximum edge length         = ',totmax(2)
      write(ilog,*) '  Minimum edge length         = ',totmin(2)
      write(*,*)    '  Minimum edge length         = ',totmin(2)
      write(ilog,*) '  Maximum aspect ratio        = ',totmax(3)
      write(*,*)    '  Maximum aspect ratio        = ',totmax(3)
    endif
 
  end subroutine precheck_mesh_quality
 
end module m_precheck_mesh_quality