mesh_tools.f90

Go to the documentation of this file.

module mesh_tools

  use choral_constants
  use choral_variables
  use real_type
  use basic_tools
  use algebra_set, only: cap_sorted_set
  use graph_mod
  use cell_mod
  use mesh_mod

  implicit none
  private

  public :: mesh_analyse                   !! TESTED
  public :: maxedgelength, minedgelength

contains

  function maxedgelength(m) result(h)
    type(mesh), intent(in) :: m
    real(RP)             :: h

    integer, dimension(:), allocatable :: row
    real(RP), dimension(3) :: x
    integer  :: ii, sz, wdt
    real(RP) :: l

    h = 0._rp
    wdt = m%clToNd%width
    call allocmem(row, wdt)

    do ii=1, m%nbCl

       call getrow(sz, row, wdt, m%clToNd, ii)

       select case(m%clType(ii))          
       case(cell_vtx)
          cycle

       case(cell_edg, cell_edg_2)   ! indicqtive if order = 2
          x = m%nd(:,row(1)) - m%nd(:,row(2))
          l = sqrt(sum(x*x))
          if (l>h) h=l


       case(cell_trg, cell_trg_2)   ! indicqtive if order = 2
          x = m%nd(:,row(1)) - m%nd(:,row(2))
          l = sqrt(sum(x*x))
          if (l>h) h=l

          x = m%nd(:,row(2)) - m%nd(:,row(3))
          l = sqrt(sum(x*x))
          if (l>h) h=l

          x = m%nd(:,row(3)) - m%nd(:,row(1))
          l = sqrt(sum(x*x))
          if (l>h) h=l

       case(cell_tet, cell_tet_2)
          x = m%nd(:,row(1)) - m%nd(:,row(2))
          l = sqrt(sum(x*x))
          if (l>h) h=l

          x = m%nd(:,row(1)) - m%nd(:,row(3))
          l = sqrt(sum(x*x))
          if (l>h) h=l

          x = m%nd(:,row(1)) - m%nd(:,row(4))
          l = sqrt(sum(x*x))
          if (l>h) h=l

          x = m%nd(:,row(2)) - m%nd(:,row(3))
          l = sqrt(sum(x*x))
          if (l>h) h=l

          x = m%nd(:,row(2)) - m%nd(:,row(4))
          l = sqrt(sum(x*x))
          if (l>h) h=l

          x = m%nd(:,row(3)) - m%nd(:,row(4))
          l = sqrt(sum(x*x))
          if (l>h) h=l

       case default
          call quit ("mesh_tools: maxEdgeLength:&
               & cell type not supported")

       end select
    end do

  end function maxedgelength

  function minedgelength(m) result(h)
    type(mesh), intent(in) :: m
    real(RP)               :: h

    integer, dimension(:), allocatable :: row
    real(RP), dimension(3) :: x
    integer  :: ii, sz, wdt
    real(RP) :: l

    h = 1e10_rp
    wdt = m%clToNd%width
    call allocmem(row, wdt)

    do ii=1, m%nbCl

       call getrow(sz, row, wdt, m%clToNd, ii)

       select case(m%clType(ii))          
       case(cell_vtx)
          cycle

       case(cell_edg, cell_edg_2)   ! indicqtive if order = 2
          x = m%nd(:,row(1)) - m%nd(:,row(2))
          l = sqrt(sum(x*x))
          if (l<h) h=l

       case(cell_trg, cell_trg_2)   ! indicqtive if order = 2
          x = m%nd(:,row(1)) - m%nd(:,row(2))
          l = sqrt(sum(x*x))
          if (l<h) h=l

          x = m%nd(:,row(2)) - m%nd(:,row(3))
          l = sqrt(sum(x*x))
          if (l<h) h=l

          x = m%nd(:,row(3)) - m%nd(:,row(1))
          l = sqrt(sum(x*x))
          if (l<h) h=l

       case(cell_tet, cell_tet_2)
          x = m%nd(:,row(1)) - m%nd(:,row(2))
          l = sqrt(sum(x*x))
          if (l<h) h=l

          x = m%nd(:,row(1)) - m%nd(:,row(3))
          l = sqrt(sum(x*x))
          if (l<h) h=l

          x = m%nd(:,row(1)) - m%nd(:,row(4))
          l = sqrt(sum(x*x))
          if (l<h) h=l

          x = m%nd(:,row(2)) - m%nd(:,row(3))
          l = sqrt(sum(x*x))
          if (l<h) h=l

          x = m%nd(:,row(2)) - m%nd(:,row(4))
          l = sqrt(sum(x*x))
          if (l<h) h=l

          x = m%nd(:,row(3)) - m%nd(:,row(4))
          l = sqrt(sum(x*x))
          if (l<h) h=l

       case default
          call quit ("mesh_tools: minEdgeLength:&
               & cell type not supported")

       end select
    end do

  end function minedgelength


  function mesh_analyse(m, verb) result(mesh_desc)
    integer   , dimension(4,6) :: mesh_desc
    type(mesh), intent(in)     :: m
    integer   , intent(in)     :: verb

    integer :: ii, cl_t, dim, nbVtx, nbNodes
    integer, dimension(CELL_MAX_NBNODES) :: row
    integer, dimension(4, m%nbNd) :: p, q

    if (choral_verb>0) write(*,*) "mesh_tools      : mesh_analyse"
    mesh_desc = 0

    
    if(.NOT.valid(m)) then
       call quit( "mesh_tools: mesh_analyse: mesh 'm' not valid" )
    end if

    p = 0; q = 0
    do ii=1, m%nbCl

       cl_t    = m%clType(ii)
       dim     = cell_dim(cl_t)
       nbvtx   = cell_nbvtx(cl_t)
       
       mesh_desc( dim + 1, 1) = mesh_desc( dim + 1, 1) + 1

       call getrow(nbnodes, row, cell_max_nbnodes, m%clToNd, ii)
       p(dim + 1, row(1:nbvtx   )) = 1
       q(dim + 1, row(1:nbnodes )) = 1
       
    end do

    do ii=1, 4
       mesh_desc(ii,2) = sum( p(ii,:) )
       mesh_desc(ii,3) = sum( q(ii,:) )
    end do

    mesh_desc(2,6) = mesh_desc(2,1) - mesh_desc(2,2)
    if (m%dim==1) goto 10

    mesh_desc(3,4) = count_2dcell_edges(m)
    mesh_desc(3,6) = mesh_desc(3,1) - mesh_desc(3,4) + mesh_desc(3,2)
    if (m%dim==2) goto 10

    mesh_desc(4,4) = count_3dcell_edges(m)
    mesh_desc(4,5) = count_3dcell_faces(m)
    mesh_desc(4,6) = mesh_desc(4,1) - mesh_desc(4,5) &
         &         + mesh_desc(4,4) - mesh_desc(4,2)

10  continue

    if (verb>0) call display_mesh_desc()  

    if (mesh_desc(1,1)/=mesh_desc(1,2)) call quit( &
         & "mesh_tools: mesh_analyse: 2" )

    if (mesh_desc(1,1)/=mesh_desc(1,3)) call quit( &
         & "mesh_tools: mesh_analyse: 3" )

  contains

    subroutine display_mesh_desc()  

      call print(m)
      write(*,*)"  0D cells"
      write(*,*)"    nb cells              =", mesh_desc(1,1)
      write(*,*)"    nb vertexes           =", mesh_desc(1,2)
      write(*,*)"    nb nodes              =", mesh_desc(1,3)

      write(*,*)"  1D cells"
      write(*,*)"    nb cells              =", mesh_desc(2,1)
      write(*,*)"    nb vertexes           =", mesh_desc(2,2)
      write(*,*)"    Euler charac          =", mesh_desc(2,6)
      write(*,*)"    nb nodes              =", mesh_desc(2,3)
      if (m%dim==1) return

      write(*,*)"  2D cells"
      write(*,*)"    nb cells              =", mesh_desc(3,1)
      write(*,*)"    nb edges              =", mesh_desc(3,4)
      write(*,*)"    nb vertexes           =", mesh_desc(3,2)
      write(*,*)"    Euler charac          =", mesh_desc(3,6)
      write(*,*)"    nb nodes              =", mesh_desc(3,3)
      if (m%dim==2) return

      write(*,*)"  3D cells"
      write(*,*)"    nb cells              =", mesh_desc(4,1)
      write(*,*)"    nb faces              =", mesh_desc(4,5)
      write(*,*)"    nb edges              =", mesh_desc(4,4)
      write(*,*)"    nb vertexes           =", mesh_desc(4,2)
      write(*,*)"    Euler charac          =", mesh_desc(4,6)
      write(*,*)"    nb nodes              =", mesh_desc(4,2)

    end subroutine display_mesh_desc

  end function mesh_analyse

  function count_3dcell_edges(m) result(nEd)
    type(mesh), intent(in)  :: m
    integer                 :: nEd
    
    integer :: wdt
    
    ned = 0
    
    if (m%dim/=3) return
    
    ! Maximal number of neiuhbour cells
    ! of the nodes of the mesh
    ! 
    wdt = m%ndToCl%width

    call cell_loop()

  contains

    subroutine cell_loop()
    
      integer, dimension(wdt, CELL_MAX_NBVTX):: vtx_toCl
      integer, dimension(CELL_MAX_NBVTX)     :: vtx_nbCl 
      integer, dimension(CELL_MAX_NBNODES)   :: cl_vtx
      integer, dimension(wdt)                :: ed_toCl
      
      integer :: ii, jj, ll, v1, v2
      integer :: cl, cl_nbVtx, cl_t
      integer :: ed_nbCl, cl2, cl2_t
      
      do cl=1, m%nbCl
         
         cl_t = m%clType(cl)
         if ( cell_dim(cl_t) /= 3 ) cycle
         
         !! cl_nbVtx  = number of vertexes of cell cl
         !! cl_vtx(:) = vertexes of cl (global indexzs)
         !!
         call getrow(ii, cl_vtx, cell_max_nbnodes, &
              &         m%clToNd, cl)
         cl_nbvtx = cell_nbvtx( cl_t )
         
         !! For the vertex ii of cell cl :
         !!   vtx_nbCl(ii)   = number of neighbour cells 
         !!   vtx_toCl(:,ii) = neighbour cells
         !!
         do ii=1, cl_nbvtx
            call getrow(vtx_nbcl(ii), vtx_tocl(:,ii), wdt, &
                 & m%ndToCl, cl_vtx(ii))
         end do
         
         ! Loop on cell cl edges
         do ii=1, cell_nbed(cl_t)
            
            !! edge ii vertexes = local indexing 
            !!                    in the cell cl
            v1 = cell_ed_vtx(1, ii, cl_t)
            v2 = cell_ed_vtx(2, ii, cl_t)
            
            
            ! Neighbour cells to the edge ii
            !
            jj = vtx_nbcl(v1)
            ll = vtx_nbcl(v2)
            call cap_sorted_set(ed_nbcl, ed_tocl, wdt, &
                 & vtx_tocl( 1: jj, v1), jj,      &
                 & vtx_tocl( 1: ll, v2), ll       )
            
            ! Test if that edge is a new face
            !
            loop: do jj=1, ed_nbcl
               cl2   = ed_tocl(jj)
               cl2_t = m%clType(cl2)
               if (cell_dim(cl2_t)==3) exit loop          
            end do loop
            if (cl2==cl) ned = ned + 1   ! new edge
            
         end do
      end do
    end subroutine cell_loop
    
  end function count_3dcell_edges


  
  function count_2dcell_edges(m) result(n)
    integer :: n    
    type(mesh), intent(in)  :: m

    integer, dimension(:), allocatable :: nnz
    type(graph) :: clToEd

    integer  :: cl
    integer  :: cl_t, nb_vtx, nb_ed, width

    if (choral_verb>1) write(*,*) &
         & "mesh_tools      : count_2DCell_edges"

    !! STEP 1: build clToEd pattern
    !!
    call allocmem(nnz, m%nbCl)
    do cl=1, m%nbCl
       cl_t = m%clType(cl)
       if ( cell_dim(cl_t) ==2 ) then
          nnz(cl) = cell_nbed(cl_t)
       else
          nnz(cl) = 0
       end if
    end do
    if (sum(nnz)==0) then
       n = 0
       return
    end if
    cltoed = graph(nnz, m%nbCl, 0)
    

    !! STEP 2 : count edges
    !!
    n     = 0
    width = m%ndToCl%width
    do cl_t=1, cell_tot_nb
       if ( cell_dim(cl_t)     /= 2 ) cycle
       if ( m%cell_count(cl_t) == 0 ) cycle

       nb_vtx = cell_nbvtx(cl_t)
       nb_ed  = cell_nbed(cl_t)
       call cell_loop()

    end do

  contains

    !! count the edges of the cells with type cl_t
    !!
    subroutine cell_loop()

      integer, dimension(nb_vtx)         :: vtx_nbCl
      integer, dimension(width, nb_vtx)  :: vtx_cl
      integer, dimension(CELL_MAX_NBVTX) :: cl_vtx
      integer, dimension(width)          :: coBord
      integer :: ed, v1, v2
      integer :: j1, j2, l1, l2

      do cl=1, m%nbCl
         l1 = m%clType(cl) 
         if (l1 /= cl_t ) cycle

         !! cell cl vertexes
         !!
         j1 = m%clToNd%row(cl)
         j2 = j1 + nb_vtx - 1
         cl_vtx(1:nb_vtx) =  m%clToNd%col( j1:j2 )

         !! neighbour cells for each vertex
         !! of the cell cl
         !!
         do v1=1, nb_vtx
            l1 = cl_vtx(v1)
            j1 = m%ndToCl%row(l1)
            j2 = m%ndToCl%row(l1+1)

            l1 = j2-j1
            vtx_nbcl(v1) = l1

            !! vtx_cl is sorted !
            vtx_cl(1:l1, v1) = m%ndToCl%col(j1:j2-1)

         end do

         !! get the co-boundary of the edges of cl
         !! and update the total number of edges
         !!
         j1 = cltoed%row(cl) 
         do ed=1, nb_ed

            !! v1, v2 = vertexes of the current edge, local  indexes
            v1 = cell_ed_vtx(1, ed, cl_t)
            v2 = cell_ed_vtx(2, ed, cl_t)

            !! list the cells sharing the vertexes v1 and v2
            !! coBord is sorted !
            l1 = vtx_nbcl(v1)
            l2 = vtx_nbcl(v2)

            call cap_sorted_set(j2, cobord, width, &
                 & vtx_cl(1:l1, v1), l1, vtx_cl(1:l2, v2), l2)

            !! determine if the current edge is a new edge
            !!
            do l1=1, j2      
               if ( nnz( cobord(l1) ) == 0 ) cycle
               exit
            end do 
            if (cobord(l1)==cl) n = n + 1

            !! increment edge index
            j1 = j1 + 1

         end do
      end do

    end subroutine cell_loop
    
  end function count_2dcell_edges

  function count_3dcell_faces(m) result(n)
    integer :: n    
    type(mesh), intent(in)  :: m

    integer, dimension(:), allocatable :: nnz
    type(graph) :: clToFc

    integer  :: cl
    integer  :: cl_t, nb_vtx, nb_fc, width

    if (choral_verb>1) write(*,*) &
         & "mesh_tools      : count_3DCell_faces"

    !! STEP 1: build clToFc pattern
    !!
    call allocmem(nnz, m%nbCl)
    do cl=1, m%nbCl
       cl_t = m%clType(cl)
       if ( cell_dim(cl_t) == 3 ) then
          nnz(cl) = cell_nbfc(cl_t)
       else
          nnz(cl) = 0
       end if
    end do
    cltofc = graph(nnz, m%nbCl, 0)


    !! STEP 2 : count faces
    !!
    n     = 0
    width = m%ndToCl%width
    do cl_t=1, cell_tot_nb
       if ( cell_dim(cl_t)     /= 3 ) cycle
       if ( m%cell_count(cl_t) == 0 ) cycle

       nb_vtx = cell_nbvtx(cl_t)
       nb_fc  = cell_nbfc(cl_t)
       call cell_loop()

    end do

  contains

    !! count the faces of the cells with type cl_t
    !!
    subroutine cell_loop()

      integer, dimension(nb_vtx)         :: vtx_nbCl
      integer, dimension(width, nb_vtx)  :: vtx_cl
      integer, dimension(CELL_MAX_NBVTX) :: cl_vtx
      integer, dimension(width)          :: coBord, aux
      integer :: fc, v1, v2
      integer :: j1, j2, l1, l2, ii

      do cl=1, m%nbCl
         l1 = m%clType(cl) 
         if (l1 /= cl_t ) cycle

         !! cell cl vertexes
         !!
         j1 = m%clToNd%row(cl)
         j2 = j1 + nb_vtx - 1
         cl_vtx(1:nb_vtx) =  m%clToNd%col( j1:j2 )

         !! neighbour cells for each vertex
         !! of the cell cl
         !!
         do v1=1, nb_vtx
            l1 = cl_vtx(v1)
            j1 = m%ndToCl%row(l1)
            j2 = m%ndToCl%row(l1+1)

            l1 = j2-j1
            vtx_nbcl(v1) = l1

            !! vtx_cl is sorted !
            vtx_cl(1:l1, v1) = m%ndToCl%col(j1:j2-1)

         end do

         !! get the co-boundary of the faces of cl
         !! and update the total number of faces
         !!
         j1 = cltofc%row(cl) 
         do fc=1, nb_fc

            !! v1, v2 = vertexes of the current face, local  indexes
            v1 = cell_fc_vtx(1, fc, cl_t)
            v2 = cell_fc_vtx(2, fc, cl_t)

            !! list the cells sharing the vertexes v1 and v2
            !! coBord is sorted !
            l1 = vtx_nbcl(v1)
            l2 = vtx_nbcl(v2)

            call cap_sorted_set(j2, cobord, width, &
                 & vtx_cl(1:l1, v1), l1, vtx_cl(1:l2, v2), l2)

            do ii=3, cell_fc_nbvtx(fc, cl_t)
               v2 = cell_fc_vtx(ii, fc, cl_t)
               l2 = vtx_nbcl(v2)

               l1 = j2
               aux(1:l1) = cobord(1:l1)

               call cap_sorted_set(j2, cobord(1:l1), l1, &
                 & aux(1:l1), l1, vtx_cl(1:l2, v2), l2)

            end do

            !! determine if the current face is a new face
            !!
            do l1=1, j2      
               if ( nnz( cobord(l1) ) == 0 ) cycle
               exit
            end do 
            if (cobord(l1)==cl) n = n + 1

            !! increment face index
            j1 = j1 + 1

         end do
      end do

    end subroutine cell_loop
    
  end function count_3dcell_faces


end module mesh_tools