t13.f90

! ------------------------------------------------------------------------------
!
!  Gmsh Fortran tutorial 13
!
!  Remeshing an STL file without an underlying CAD model
!
! ------------------------------------------------------------------------------

program t13

use, intrinsic :: iso_c_binding
use gmsh

implicit none
type(gmsh_t) :: gmsh
integer(c_int) :: ret
real(c_double), parameter :: lc = .1
character(len=GMSH_API_MAX_STR_LEN) :: cmd

call gmsh%initialize()
! Create ONELAB parameters with remeshing options:
call gmsh%onelab%set('['//&
  '{'//&
    '"type":"number",'//&
    '"name":"Parameters/Angle for surface detection",'//&
    '"values":[40],'//&
    '"min":20,'//&
    '"max":120,'//&
    '"step":1'//&
  '},'//&
  '{'//&
    '"type":"number",'//&
    '"name":"Parameters/Create surfaces guaranteed to be parametrizable",'//&
    '"values":[0],'//&
    '"choices":[0, 1]'//&
  '},'//&
  '{'//&
    '"type":"number",'//&
    '"name":"Parameters/Apply funny mesh size field?",'//&
    '"values":[0],'//&
    '"choices":[0, 1]'//&
  '}'//&
']')

! Create the geometry and mesh it:
call createGeometryAndMesh()

call get_command(cmd)
if (index(cmd, "-nopopup") == 0) then
    call gmsh%fltk%initialize()
    do while (transfer(gmsh%fltk%isAvailable(), .true.) .and. checkForEvent())
        call gmsh%fltk%wait()
    end do
end if

call gmsh%finalize()

contains

subroutine createGeometryAndMesh()
    real(c_double) :: curveAngle, angle, forceParametrizablePatches
    real(c_double), allocatable :: tmp(:)
    integer(c_int), allocatable :: s(:,:)
    integer(c_int) :: l, f
    logical :: includeBoundary
    ! Clear all models and merge an STL mesh that we would like to remesh (from
    ! the parent directory):
    call gmsh%clear()
    call gmsh%merge('../t13_data.stl')

    ! We first classify ("color") the surfaces by splitting the original surface
    ! along sharp geometrical features. This will create new discrete surfaces,
    ! curves and points.

    ! Angle between two triangles above which an edge is considered as sharp,
    ! retrieved from the ONELAB database (see below):
    call gmsh%onelab%getNumber('Parameters/Angle for surface detection', tmp)
    angle = tmp(1); deallocate(tmp)

    ! For complex geometries, patches can be too complex, too elongated or too
    ! large to be parametrized; setting the following option will force the
    ! creation of patches that are amenable to reparametrization:
    call gmsh%onelab%getNumber(&
        'Parameters/Create surfaces guaranteed to be parametrizable', tmp)
    forceParametrizablePatches = tmp(1); deallocate(tmp)

    ! For open surfaces include the boundary edges in the classification
    ! process:
    includeBoundary = .true.

    ! Force curves to be split on given angle:
    curveAngle = 180

    call gmsh%model%mesh%classifySurfaces(angle * M_PI / 180., includeBoundary, &
                                          transfer(forceParametrizablePatches, .false.), &
                                          curveAngle * M_PI / 180.)

    ! Create a geometry for all the discrete curves and surfaces in the mesh, by
    ! computing a parametrization for each one
    call gmsh%model%mesh%createGeometry()

    ! Note that if a CAD model (e.g. as a STEP file, see `t20.f90') is available
    ! instead of an STL mesh, it is usually better to use that CAD model instead
    ! of the geometry created by reparametrizing the mesh% Indeed, CAD
    ! geometries will in general be more accurate, with smoother
    ! parametrizations, and will lead to more efficient and higher quality
    ! meshing. Discrete surface remeshing in Gmsh is optimized to handle dense
    ! STL meshes coming from e.g. imaging systems, where no CAD is available; it
    ! is less well suited for the poor quality STL triangulations (optimized for
    ! size, with e.g. very elongated triangles) that are usually generated by
    ! CAD tools for e.g. 3D printing.

    ! Create a volume from all the surfaces
    call gmsh%model%getEntities(s, dim=2)
    l = gmsh%model%geo%addSurfaceLoop(s(2, :))
    ret = gmsh%model%geo%addVolume([l])

    call gmsh%model%geo%synchronize()

    ! We specify element sizes imposed by a size field, just because we can :-)
    f = gmsh%model%mesh%field%add("MathEval")
    call gmsh%onelab%getNumber('Parameters/Apply funny mesh size field?', tmp)
    if (int(tmp(1)) > 0) then
        call gmsh%model%mesh%field%setString(f, "F", "2*Sin((x+y)/5) + 3")
    else
        call gmsh%model%mesh%field%setString(f, "F", "4")
    end if
    call gmsh%model%mesh%field%setAsBackgroundMesh(f)

    call gmsh%model%mesh%generate(3)
    call gmsh%write('t13.msh')

end subroutine createGeometryAndMesh


! Launch the GUI and handle the "check" event to recreate the geometry and mesh
! with new parameters if necessary:
logical function checkForEvent()
    character(len=GMSH_API_MAX_STR_LEN), allocatable :: action(:)
    checkForEvent = .false.
    call gmsh%onelab%getString("ONELAB/Action", action)
    if (size(action) > 0) then
        if (trim(action(1)) == "check") then
            call gmsh%onelab%setString("ONELAB/Action", [c_null_char])
            call createGeometryAndMesh()
            call gmsh%graphics%draw()
        end if
    end if
    checkForEvent = .true.
end function checkForEvent

end program t13