mesh_mmg.F90

Go to the documentation of this file.

module mesh_mmg

  use choral_constants
  use choral_variables
  use basic_tools
  use real_type
  use graph_mod
  use cell_mod
  use mesh_mod

  implicit none
  private

#ifdef WMMG 
#include "mmg/mmg2d/libmmg2df.h"

  public :: mmg2d_retrievemesh, mmg2d_passmesh
#endif

  public :: remesh_levelset


contains

  !! Remesh respectively to a levelSet
  !!
  !!   m   = initial mesh / output mesh
  !!   val = level values at the vertexes of the mesh 'm'
  !!
  !! Output mesh has edges along the level set 0
  !!
  !! Optional parameters
  !!   verb  = verbisity
  !!   hausd = hausdorff distance
  !!   hMax  = maximal edge length
  !!
  subroutine remesh_levelset(m, val, verb, hausd, hMax)
    type(mesh), intent(inout) :: m
    real(DP)  , dimension(:)  :: val
    
    integer , optional, intent(in) :: verb
    real(DP), optional, intent(in) :: hausd
    real(DP), optional, intent(in) :: hMax

#ifdef WMMG 
    integer          :: ii, ier
    mmg5_data_ptr_t  :: mmgmesh, mmgsol

    if (choral_verb>0) write(*,*) &
         & "mesh_mmg        : reMesh_levelSet"

    if (.NOT.valid(m)) call quit("mesh_mmg: reMesh_levelSet:&
         & mesh 'm' not valid")

    if (size(val,1) /= m%nbNd) call quit("mesh_mmg: reMesh_levelSet:&
         & wrong size for 'val'")


    !! Adapt the code to the mesh dimension
    select case(m%dim)
    case(2)
       if ( m%cell_Count(cell_trg)==0 ) call quit(&
            & "mesh_mmg: reMesh_levelSet:&
            & dimension 2 mesh withour triangles")

       !! #####################
       !!
       !!  DATA TRANSFER TO MMG
       !!
       !! mesh
       call mmg2d_passmesh(mmgmesh, mmgsol, m)
       !!
       !! level values at the mesh vertexes
       call mmg2d_set_solsize(mmgmesh, mmgsol, mmg5_vertex, m%nbNd, mmg5_scalar, ier)
       do ii=1, m%nbNd
          call mmg2d_set_scalarsol(mmgsol, val(ii), ii, ier)
          if (ier/=1) exit
       end do
       !!
       !! clear the input mesh
       call clear(m)

       !! #####################
       !!
       !!  SET MMG PARAMETERS
       !!
       !!  verbosity, default = 0
       !!
       if ( present(verb) ) then
          call mmg2d_set_iparameter(mmgmesh,mmgsol,&
               & mmg2d_iparam_verbose, verb, ier)
       else
          call mmg2d_set_iparameter(mmgmesh,mmgsol,&
               & mmg2d_iparam_verbose, 0, ier)
          
       end if
       !!
       !!  Hausdorff distance, default = 0.01
       !!
       if ( present(hausd) ) then
          call mmg2d_set_dparameter(mmgmesh, mmgsol, & 
               & mmg2d_dparam_hausd, hausd, ier)          
       else
          call mmg2d_set_dparameter(mmgmesh, mmgsol, & 
               & mmg2d_dparam_hausd, 0.01_dp, ier)          
       end if
       !!
       !!  Hausdorff distance, default = 0.01
       !!
       if ( present(hmax) ) then
          call mmg2d_set_dparameter(mmgmesh, mmgsol, & 
               & mmg2d_dparam_hmax, hmax, ier)          
       end if
       !!
       !! ISOVAL
       !!
       call mmg2d_set_iparameter(mmgmesh,mmgsol,mmg2d_iparam_iso,1,ier)

       !! #####################
       !!
       !!  REMESH
       !!
       call mmg2d_mmg2dls(mmgmesh, mmgsol, ier)

       ! call mmg2d_saveMesh(mmgMesh,"fin.mesh",8, ier)
       ! call mmg2d_saveSol(mmgMesh,mmgSol,"fin.sol",7, ier)

       !! #####################
       !!
       !!  END
       !!
       !! retrieve the mesh from MMG
       call mmg2d_retrievemesh(m, mmgmesh, mmgsol)

       !! deallocate the mmg structures
       !!
       call mmg2d_free_all(mmg5_arg_start, &
            mmg5_arg_ppmesh,mmgmesh,mmg5_arg_ppmet,mmgsol, &
            mmg5_arg_end)


    case(3)
       if ( m%cell_Count(cell_tet)==0 ) call quit("mesh_mmg: reMesh_levelSet:&
         & dimension 3 mesh withour tetrahedra")

       call quit('mesh_mmg: reMesh_levelSet:&
         & dim 3 to be done')

    case default
       call quit("mesh_mmg: reMesh_levelSet: wrong dimension")

    end select

#else

    !! preproc variable WMMG not declared
    !!
    call quit("mesh_mmg: reMesh_levelSet: set 'WITH_MMG=1' in the&
         & file 'CONFIG.in', re-configure (mmg must be installed)")
#endif

  end subroutine remesh_levelset

  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!
  !!
  !!   HERE STARTS THE CODE INCLUDING MMG DATA / FUNCTIONS
  !!
  !!
  !!
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#ifdef WMMG 

  subroutine mmg2d_passmesh(mmgMesh, mmgSol, m1)  
    mmg5_data_ptr_t, intent(out) :: mmgmesh, mmgsol
    type(mesh)     , intent(in)  :: m1

    integer :: ier, nNd, nEd, nTr, cl, ii, j1

    if (choral_verb>0) write(*,*) &
         & "mesh_mmg        : mmg2D_passMesh"

    mmgmesh = 0
    mmgsol  = 0

    call mmg2d_init_mesh(mmg5_arg_start, &
         mmg5_arg_ppmesh, mmgmesh, mmg5_arg_ppmet, mmgsol, &
         mmg5_arg_end)

    nnd = m1%nbNd
    ned = m1%cell_Count(cell_edg)
    ntr = m1%cell_Count(cell_trg)

    call mmg2d_set_meshsize(mmgmesh, nnd, ntr, ned, ier)

    !! define mesh vertexes
    !!
    do ii=1, nnd
       call mmg2d_set_vertex(mmgmesh, &
            &m1%nd(1,ii), m1%nd(2,ii), 0,  ii, ier)
    end do

    !! vertex cells --> required vertexes
    !!
    do cl=1, m1%nbCl
       if (m1%clType(cl) == cell_vtx) then

          !! ii = index of the vertex in m%nd
          j1 = m1%clToNd%row(cl)
          ii = m1%clToNd%col(j1)

          call mmg2d_set_vertex(mmgmesh, &
               &m1%nd(1,ii), m1%nd(2,ii), m1%clTag(cl),  ii, ier)

          call mmg2d_set_requiredvertex(mmgmesh, ii, ier)

       end if
    end do


    !! define mesh edges
    !!
    ii = 0
    do cl=1, m1%nbCl
       if (m1%clType(cl) == cell_edg) then
          ii = ii + 1
          j1 = m1%clToNd%row(cl)

          call mmg2d_set_edge(mmgmesh, &
               & m1%clToNd%col(j1), m1%clToNd%col(j1+1), &
               & m1%clTag(cl),  ii, ier)

       end if
    end do

    !! define mesh triangles
    ii = 0
    do cl=1, m1%nbCl
       if (m1%clType(cl) == cell_trg) then
          ii = ii + 1
          j1 = m1%clToNd%row(cl)

          call mmg2d_set_triangle(mmgmesh, &
               & m1%clToNd%col(j1), m1%clToNd%col(j1+1), m1%clToNd%col(j1+2),&
               & m1%clTag(cl),  ii, ier)
       end if
    end do

  end subroutine mmg2d_passmesh

  subroutine mmg2d_retrievemesh(m2, mmgMesh, mmgSol)  
    type(mesh)     , intent(inout) :: m2
    mmg5_data_ptr_t, intent(inOut) :: mmgmesh, mmgsol

    logical, dimension(:)  , allocatable :: vtx_req
    integer, dimension(:)  , allocatable :: vtx_ref
    integer, dimension(:,:), allocatable :: clTab

    integer  :: ier, nNd, nEd, nTr, nb_vtx_cell, ii, jj, ll
    integer  :: v1, v2, v3, req, ref
    real(DP) :: x, y

    if (choral_verb>0) write(*,*) &
         & "mesh_mmg        : mmg2D_retrieveMesh"

    call clear(m2)

    call mmg2d_get_meshsize(mmgmesh, nnd, ntr, ned, ier)

    call allocmem(m2%nd, 3, nnd)
    call allocmem(vtx_req, nnd)
    call allocmem(vtx_ref, nnd)
    nb_vtx_cell = 0
    do ii=1, nnd
       call mmg2d_get_vertex(mmgmesh, x, y, vtx_ref(ii), ll, req, ier)
       m2%nd(1,ii) = x
       m2%nd(2,ii) = y       
       m2%nd(3,ii) = 0.0_rp       

       if (req/=0) then
          vtx_req(ii) = .true.
          nb_vtx_cell  = nb_vtx_cell + 1
       else
          vtx_req(ii) = .false.
       end if
    end do

    !! total number of cells = required vertexes + edges + triangles
    ii = nb_vtx_cell + ned + ntr

    !! allocations to build the mesh m2
    call allocmem(m2%clType, ii)
    call allocmem(m2%clTag , ii)
    call allocmem(cltab, 3, ii )

    !! set cell types
    ii = nb_vtx_cell
    m2%clType(1:ii)    = cell_vtx
    jj = ii + ned
    m2%clType(ii+1:jj) = cell_edg
    m2%clType(jj+1:)   = cell_trg

    !! set cell nodes
    !!
    !!   1- vertex cells
    !!
    jj = 0
    do ii=1, nnd
       if (vtx_req(ii)) then
          jj = jj + 1
          cltab(1, jj) = ii
          m2%clTag(jj) = vtx_ref(ii)
       end if
    end do
    call freemem(vtx_req)
    call freemem(vtx_ref)
    !!
    !!   2- edge cells
    !!
    jj = nb_vtx_cell
    do ii=1, ned

       !! index of the edge ii in the mesh m2
       jj = jj + 1

       call mmg2d_get_edge(mmgmesh, v1, v2, ref, ll, req, ier)
       cltab(1, jj) = v1
       cltab(2, jj) = v2
       m2%clTag(jj) = ref

    end do
    !!
    !!   3- triangle cells
    !!
    jj = nb_vtx_cell + ned 
    do ii=1, ntr

       !! index of the edge ii in the mesh m2
       jj = jj + 1

       call mmg2d_get_triangle(mmgmesh, v1, v2, v3, ref, req, ier)
       cltab(1, jj) = v1
       cltab(2, jj) = v2
       cltab(3, jj) = v3
       m2%clTag(jj) = ref
    end do

    !! end the mesh m2 construction
    !!
    call mesh_create_cltab(m2, cltab)
    call mesh_create_end(m2)

    if (.NOT.valid(m2)) call quit('mesh_mmg: &
         & mmg2d_retrieveMesh: construction not valid')

  end subroutine mmg2d_retrievemesh

#endif

end module mesh_mmg