mesh_mod.F90

Go to the documentation of this file.

 

module mesh_mod


  use choral_constants
  use choral_variables
  use real_type
  use basic_tools
  use abstract_interfaces, only: r3tor
  use algebra_lin
  use cell_mod
  use graph_mod

  implicit none
  private

  public :: mesh
  public :: clear, print, valid
  public :: closest_node                 !! TESTED
  public :: write
  public :: getndcoord
  public :: flag_mesh_cells
  public :: cell_orientation
  public :: interface_orientation
  public :: mesh_create_cltab, mesh_create_end
  public :: mesh_clear_2

  !       %----------------------------------------%
  !       |                                        |
  !       |          DERIVED TYPE                  |
  !       |                                        |
  !       %----------------------------------------%
  type mesh

     integer :: nbcl=0
     integer :: nbnd=0
     integer :: dim=0
     integer :: nbitf = 0

     real(RP), dimension(:,:), allocatable :: nd

     integer , dimension(:)  , allocatable :: cltype

     integer , dimension(:)  , allocatable :: cltag

     integer, dimension(CELL_TOT_NB) :: cell_count = 0

     type(graph) :: cltond
     type(graph) :: ndtocl
     type(graph) :: itftocl
     type(graph) :: cltoitf

   contains

     final :: mesh_clear

  end type mesh

  !       %----------------------------------------%
  !       |                                        |
  !       |       GENERIc SUBROUTINES              |
  !       |                                        |
  !       %----------------------------------------%
  interface clear
     module procedure mesh_clear
  end interface clear

  interface mesh
     module procedure readmesh, mesh_create_1d
  end interface mesh

  interface write
     module procedure mesh_write
  end interface write

  interface valid
     module procedure mesh_valid
  end interface valid

  interface print
     module procedure mesh_print
  end interface print

contains


  subroutine mesh_clear(m)
    type(mesh), intent(inout) :: m

    !! clears the data tn the mesh built by 'call create(...)'
    !!
    m%nbCl  = 0
    m%nbNd  = 0
    call freemem(m%nd)
    call freemem(m%clType) 
    call freemem(m%clTag) 
    call clear(m%clToNd)
    

    !! clears the data tn the mesh built after 'call create(...)'
    !!
    call mesh_clear_2(m)

  end subroutine mesh_clear

  subroutine mesh_clear_2(m)
    type(mesh), intent(inout) :: m

    m%dim   = 0
    m%nbItf = 0

    call clear(m%ndToCl)

    call clear(m%itfToCl)
    call clear(m%clToItf)

    m%cell_count = 0

  end subroutine mesh_clear_2


  subroutine mesh_print(m) 
    type(mesh), intent(in) :: m

    integer :: ii

    write(*,*)"mesh_mod        : print"
    if (valid(m)) then
       write(*,*)"  Status                         = valid"
    else
       write(*,*)"  Status                         = not valid"
    end if

    write(*,*)"  Number of nodes                =", m%nbNd

    if (m%nbItf>0) then
       write(*,*)"  Interfaces assembled"
       write(*,*)"  Number of interfaces           =", m%nbItf
    else
       write(*,*)"  Interfaces not assembled"
    end if

    write(*,*)"  Number of cells                =", m%nbCl
    write(*,*)"  Number of cells per cell type"
    do ii=1, cell_tot_nb

       if ( m%cell_count(ii) == 0 ) cycle
       
       write(*,*)"                            "&
            &//cell_name(ii)//&
            &" =", m%cell_count(ii)
       
    end do

  end subroutine mesh_print


  function mesh_valid(m) result(b)
    type(mesh), intent(in) :: m
    logical :: b

    b = .false.
    if (.NOT. ( allocated(m%nd)     )) return
    if (.NOT. ( allocated(m%clType) )) return
    if (.NOT. ( allocated(m%clTag) )) return
    if ( valid(m%clToNd) < 0 ) return
    if ( valid(m%ndToCl) < 0 ) return

    b = ( m%nbCl>0 ) .AND. ( m%nbNd>0 )
    b = b .AND. ( m%nbCl == sum(m%cell_count) )
    b = b .AND. ( m%dim>0  ) .AND. ( m%dim<4 )
    b = b .AND. ( all( shape(m%nd)  == (/3,m%nbNd/)) )
    b = b .AND. ( size(m%clType,1) == m%nbCl )
    b = b .AND. ( size(m%clTag,1) == m%nbCl )
    b = b .AND. ( m%clToNd%nl == m%nbCl )
    b = b .AND. ( m%ndToCl%nl == m%nbNd )
    
    ! mesh with interfaces
    if (m%nbItf>0) then
       b =  b .AND. ( valid(m%clToItf) > 0 )
       b =  b .AND. ( valid(m%itfToCl) > 0 )
       if (.NOT.b) return

       b =  b .AND. ( m%clToItf%nl == m%nbCl )
       b =  b .AND. ( m%itfToCl%nl == m%nbItf)
    end if

  end function mesh_valid


  function readmesh(fileName, fileFormat) result(m)
    type(mesh)                    :: m
    character(len=*), intent(in)  :: fileName, fileFormat

    logical :: b
    
    if (choral_verb>0) write(*,*)&
         &"mesh_mod        : readMesh       = "//trim(filename)
    call clear(m)

    inquire(file=filename, exist=b)
    if (.NOT.b) call quit( 'mesh_mod: readMesh: uncorrect file ')

    open(unit=10, file=filename)
    select case(trim(fileformat))
    case("gmsh")
       call gmesh_readmesh(m)

    case default
       call quit( 'mesh_mod: readMesh: uncorrect file format' )
    end select
    close(10)

    call mesh_create_end(m)

    if (.NOT.valid(m)) call quit('mesh_mod: &
         & readMesh: construction not valid')

  end function readmesh

  subroutine gmesh_readmesh(m)
    type(mesh), intent(out) :: m

    integer :: ii, jj, kk, ll, nn
    integer, dimension(CELL_MAX_NBNODES + 8) :: tab  
    integer, dimension(:,:), allocatable     :: clTab

    !! reading the number of nodes
    do ii=1, 4
       read(10,*) 
    end do
    read(10,*) nn
    if (nn<=0) call quit( "mesh_mod: gmesh_readMesh: 1" )

    !! reading the nodes
    m%nbNd = nn
    call allocmem(m%nd, 3, nn)
    do ii=1, nn
       read (10, *) jj, m%nd(1:3,ii)
    end do

    !! reading the number of cells
    read(10,*) ; read(10,*) 
    read(10,*) nn
    if (nn<=0) call quit( "mesh_mod: gmesh_readMesh: 2" )
    m%nbCl = nn

    !! reading cell nodes
    call allocmem(m%clType, nn)
    call allocmem(m%clTag , nn)
    call allocmem(cltab, cell_max_nbnodes, nn )
    do ii=1, nn
       ! jj=cell index, kk=cel type, ll=cell number of tags
       read(10,*) jj, kk, ll, m%clTag(ii)

       ! cell type        
       kk = gmsh_to_cell(kk)
       m%clType(ii) = kk

       ! cell number of nodes
       jj = cell_nbnodes(kk)

       ! read cell nodes
       kk = 3 + ll + 1   ! first nd index 
       ll = 3 + ll + jj  ! last  nd index 
       backspace(10)
       read(10,*) tab(1:ll)   

       ! store cell nodes
       cltab(1:jj, ii) = tab(kk:ll)

    end do

    call mesh_create_cltab(m, cltab)

  end subroutine gmesh_readmesh

  subroutine mesh_create_cltab(m, clTab)
    type(mesh), intent(inout)  :: m
    integer   , dimension(:,:) :: clTab

    integer, dimension(:), allocatable :: nnz
    integer :: ii, jj, ll

    call allocmem(nnz, size(cltab, 2))

    ! cell number of nodes
    do ii=1, size(cltab, 2)

       ll = m%clType(ii) 
       jj = cell_nbnodes(ll)
       nnz(ii) = jj

    end do

    ! create cell to node graph
    m%clToNd = graph(nnz, m%nbCl, m%nbNd)

    ! fill in graph
    !
    jj = 1
    do ii=1, size(cltab, 2)

       ll = nnz(ii) 
       m%clToNd%col(jj: jj+ll-1) = cltab(1:ll, ii)   
       jj = jj + ll

    end do

  end subroutine mesh_create_cltab


  subroutine mesh_write(m, fileName, fileFormat, cell_tags)
    type(mesh)       , intent(in)    :: m
    character(len=*) , intent(in)    :: fileName, fileFormat
    logical, optional, intent(in)    :: cell_tags

    if (choral_verb>0) write(*,*)&
         &"mesh_mod        : writee         = "//trim(filename)

    if (.NOT.valid(m)) call quit('mesh_mod: mesh_write: mesh not valid')

    open(unit=10, file=filename)
    select case(trim(fileformat))
    case("gmsh")
       call gmesh_writemesh(m, cell_tags)

    case default
       call quit( 'mesh_mod: mesh_write: uncorrect file format' )
    end select
    close(10)

  end subroutine mesh_write

  subroutine gmesh_writemesh(m, cell_tags)
    type(mesh)       , intent(in) :: m
    logical, optional, intent(in) :: cell_tags

    integer :: ii, tt, j1, j2

    write(10,'(A11)') "$MeshFormat"
    write(10,'(A7)') "2.2 0 8"
    write(10,'(A14)') "$EndMeshFormat"
    write(10,'(A6)') "$Nodes"
    write(10,*) m%nbNd
    do ii=1, m%nbNd
       write(10, *) ii, m%nd(1:3,ii)
    end do
    write(10,'(A9)') "$EndNodes"
    write(10,'(A9)') "$Elements"
    write(10,*) m%nbCl
    j1 = 1
    do ii=1, m%nbCl

       j2 = m%clType(ii)
       tt = cell_to_gmsh(j2)

       j2 = m%clToNd%row(ii+1) - 1

       write(10, *) ii, tt, 2, m%clTag(ii), 1, m%clToNd%col(j1:j2)

       j1 = j2 + 1
    end do
    write(10,'(A12)') "$EndElements"

    if (present(cell_tags)) then
       if (cell_tags) then

          write(10,fmt="(A)") '$ElementData'
          write(10,fmt="(A)") '1'
          write(10,fmt="(A)") "Cell tags"
          write(10,fmt="(A)") '1'
          write(10,*) 0.0_rp
          write(10,fmt="(A)") '3'
          write(10,*) 1
          write(10,*) 1
          write(10,*) m%nbCl
          do ii=1, m%nbCl
             write(10,*) ii, m%clTag(ii)  
          end do
          write(10,fmt="(A)") '$EndElementData'

       end if
    end if
  end subroutine gmesh_writemesh

  function mesh_create_1d(a, b, n) result(m)
    type(mesh)           :: m
    real(RP), intent(in) :: a, b
    integer , intent(in) :: n

    real(RP) :: h, xi
    integer  :: ii, j1
    integer, dimension(:), allocatable :: p 

    if (choral_verb>0) write(*,*)"mesh_mod        : create_1D"
    call clear(m)
    
    h = (b-a)/real(n, RP) 
    call allocmem(m%nd, 3, n+1)

    m%nbCl = n + 2  ! n interval cells + 2 vertex cells = the boundary
    m%nbNd = n + 1

    do ii=0, n
       xi = a + real(ii, rp)*h
       m%nd(:,ii+1) = (/ xi, 0._rp, 0._rp/)
    end do

    call allocmem(m%clType, n+2)    
    m%clType(1:2) = cell_vtx
    m%clType(3: ) = cell_edg

    call allocmem(m%clTag , n+2)    
    m%clTag = 0

    !! cell to nodes graph
    !!
    call allocmem(p, n+2)
    p(1:2) = 1
    p(3: ) = 2
    m%clToNd = graph(p, m%nbCl, m%nbNd)
    m%clToNd%col(1) = 1
    m%clToNd%col(2) = n
    j1 = 3
    do ii=1, n
       m%clToNd%col(j1:j1+1) = (/ii, ii+1/)
       j1 = j1 + 2
    end do

    call mesh_create_end(m)

    if (.NOT.valid(m)) call quit('mesh_mod: &
         & mesh_create_1D: construction not valid')

  end function mesh_create_1d


  subroutine mesh_create_end(m)
    type(mesh), intent(inout) :: m

    integer :: cl, ct, dim
    
    if (choral_verb>0) write(*,*) "mesh_mod        : create_end     =", &
         & size(m%clType,1),"  cells "

    m%nbCl = size(m%clType,1)
    m%nbNd = size(m%nd,2)
    call transpose(m%ndToCl, m%clToNd)

    !! Set m%cell_count and m%dim
    !!
    m%cell_count = 0
    m%dim        = 0
    do cl=1, m%nbCl

       ct = m%clType(cl)
       m%cell_count(ct) = m%cell_count(ct) + 1

       dim = cell_dim(ct)
       if (dim>m%dim) m%dim = dim
       
    end do

  end subroutine mesh_create_end

  subroutine getndcoord(ndCoord, s, ndInd, n, m)
    integer                 , intent(in)  :: n, s
    real(RP), dimension(3,s), intent(out) :: ndCoord
    integer , dimension(n)  , intent(in)  :: ndInd
    type(mesh)              , intent(in)  :: m
    
    integer :: jj

    do jj = 1, n
       ndcoord(:,jj) = m%nd( :, ndind(jj) )          
    end do

  end subroutine getndcoord


  subroutine flag_mesh_cells(flag, m, dim, f)
    logical, dimension(:), allocatable :: flag
    type(mesh)           , intent(in)  :: m
    integer              , intent(in)  :: dim
    procedure(r3tor)     , optional    :: f

    integer , dimension(   CELL_MAX_NBNODES) :: p
    real(RP), dimension(3, CELL_MAX_NBNODES) :: X
    integer  :: cl, nbNd, ii, cpt, ct
    real(RP) :: f_x
    logical  :: bool

    if (choral_verb>1) write(*,*) &
         & "mesh_mod        : flag_mesh_cells"
    if (.NOT.valid(m)) call quit(&
         &'mesh_mod: flag_mesh_cells: mesh not valid')

    call allocmem(flag, m%nbCl)
    flag = .false.
    cpt  = 0

    do cl=1, m%nbCl
       ct = m%clType(cl)
       if (cell_dim(ct) == dim) then
          flag(cl) = .true.
          cpt = cpt + 1
       end if
    end do

    if (present(f)) then
       do cl=1, m%nbCl
          if (.NOT.flag(cl)) cycle
          
          call getrow(nbnd, p, cell_max_nbnodes, m%clToNd, cl)

          call getndcoord(x, cell_max_nbnodes, p(1:nbnd), nbnd, m)

          bool = .true.
          loop_nodes: do ii=1, nbnd

             f_x = f( x(:,ii) )

             bool = ( f_x >= -real_tol )
             if (.NOT.bool) exit loop_nodes

          end do loop_nodes
          flag(cl) = bool

          if (.NOT.bool) cpt = cpt - 1

       end do
    end if

    !! Display
    !!
    if (choral_verb>1) write(*,*) &
         & "  Number of flaged cells         =", cpt

  end subroutine flag_mesh_cells


  subroutine closest_node(V,  dist, pt, m, nDToNd, start)
    integer , dimension(:)  , allocatable :: V
    real(RP)                , intent(out) :: dist 
    real(RP), dimension(:,:), intent(in)  :: pt
    type(mesh)              , intent(in)  :: m
    type(graph)             , intent(in)  :: ndToNd   
    integer , optional      , intent(in)  :: start

    integer , dimension(  ndToNd%width) :: p
    real(RP), dimension(  ndToNd%width) :: d
    real(RP), dimension(3,ndToNd%width) :: Y
    real(RP), dimension(3)   :: x
    integer , dimension(1)   :: ind
    real(RP)  :: a, b
    integer   :: ii, jj, nd, cpt, wdt, sz
    logical   :: bool

    if (size(pt,1)/=3) call quit( "mesh_mod:&
         & closest_node: wrong argument 'pt'" )

    call allocmem(v, size(pt,2))
    cpt  = 0
    dist = 0._rp
    wdt = ndtond%width
    
    ! starting node
    nd = m%nbNd
    if (present(start))  nd = start

    do ii=1, size(pt,2)
       x = pt(:,ii)

       ! starting distance
       a = sum( (m%nd(:,nd) - x)**2)

       bool = .true.
       do while(bool)
          cpt = cpt + 1

          ! p(1:sz) = neighbour nodes of nd
          call getrow(sz, p, wdt, ndtond, nd)

          ! Y = corresponding coordinates
          ! b = minimal distance**2
          call getndcoord(y, wdt, p(1:sz), sz, m)
          do jj=1, sz
             y(:,jj) =  y(:,jj)-x(:)
          end do
          y(:,1:sz) = y(:,1:sz)**2
          do jj=1, sz
             d(jj) = sum( y(:,jj) )
          end do
          ind = minloc( d(1:sz) )
          b   = d(ind(1)) 

          ! one closer node found ?
          if (b < a - real_tol) then
             ! YES = restart search after "do while(bool)"        
             nd = p(ind(1))
             a   = b

          else
             ! NO = end the research
             v(ii) = nd
             bool=.false.
             if (a> dist)  dist=a
          end if
       end do
    end do

    a = re(cpt, size(pt,2))
    if (choral_verb>2) write(*,*)&
         &"mesh_mod        : closest_node   = &
         &averagge neighbour nodes used =", real(a,4)

  end subroutine closest_node

  function cell_orientation(m, cl) result(bool)
    logical                :: bool
    type(mesh), intent(in) :: m
    integer   , intent(in) :: cl

    integer, dimension(:), allocatable :: p
    integer :: wdt, sz
    real(RP), dimension(3) :: AB, AC, v

    wdt = m%clToNd%width
    call allocmem(p, wdt)

    select case(m%clType(cl))
    case(cell_edg, cell_edg_2)
       !!
       !! orientation of an edge in R

       call getrow(sz, p, wdt, m%clToNd, cl)
       ab = m%nd(:,p(2)) - m%nd(:,p(1))
       bool = (ab(1)>0._rp)

    case(cell_trg, cell_trg_2)
       !!
       !! orientation of a triangle in R^2

       call getrow(sz, p, wdt, m%clToNd, cl)
       ab = m%nd(:,p(2)) - m%nd(:,p(1))
       ac = m%nd(:,p(3)) - m%nd(:,p(1))

       v = crossprod_3d(ab, ac)
       bool = (v(3)>0._rp)

    case default
       bool = .false.
       call quit( 'mesh_mod: cell_orientation:&
               & cell type not supported' )

    end select

  end function cell_orientation


  subroutine interface_orientation(nb_itf, eps, n, m, cl)
    integer                 , intent(out) :: nb_itf    
    integer   , dimension(n), intent(out) :: eps
    type(mesh), intent(in)  :: m
    integer   , intent(in)  :: cl, n

    integer, dimension(2) :: p2
    integer, dimension(CELL_MAX_NBITF) :: p1

    integer :: ii, coBound

#if DBG>0
    if (m%nbItf<=0) call quit( &
         &   "mesh_mod: interface_orientation:&
         & try 'call define_interfaces(m)'" )
#endif

    ! interfaces of cell cl
    call getrow(nb_itf, p1, cell_max_nbitf, m%clToItf, cl)
    eps(1:nb_itf) = -1

    ! loop on the interfaces of cell cl
    do ii=1, nb_itf

       ! coboundary of the interface p1(ii)
       call getrow(cobound, p2, 2, m%itfToCl, p1(ii))

       ! tests if the edge p1(ii) of cell cl has 
       ! the local orientation associated with cell cl
       if (p2(1)==cl) eps(ii) = 1

    end do

  end subroutine interface_orientation

end module mesh_mod