Replace Shewchuk Triangle with fortfem's MIT-licensed Triangle-compatible mesher

CMakeLists.txt

@@ -58,7 +58,7 @@ find_package(FFTW REQUIRED COMPONENTS DOUBLE_LIB)
 
 include(ExternalProject)
 include(SetupCODE)
-include(SetupTriangle)
+include(SetupFortfem)
 if(WITH_MFEM)
   include(SetupMFEM)
   add_definitions(-DUSE_MFEM)
@@ -94,6 +94,7 @@ set_source_files_properties(${COMMON_FORTRAN_SRC_FILES}
 
 set(MEPHIT_FORTRAN_SRC_FILES
   src/mephit_conf.f90
+  src/mephit_triangulate.f90
   src/mephit_util.f90
   src/mephit_mesh.F90
   src/mephit_pert.f90
@@ -118,7 +119,7 @@ if(WITH_MFEM)
   )
 endif()
 set_source_files_properties(src/mephit_run.c ${MEPHIT_C_SRC_FILES} ${MEPHIT_CPP_SRC_FILES}
-	PROPERTIES COMPILE_FLAGS "${CMAKE_C_FLAGS} -DREAL=double -I${TRIANGLE_INCLUDE_DIR} -I${SUITESPARSE_INCLUDE_DIRS} -L${TRIANGLE_LIB_DIR}")
+	PROPERTIES COMPILE_FLAGS "${CMAKE_C_FLAGS} -I${SUITESPARSE_INCLUDE_DIRS}")
 
 
 ### Define library
@@ -129,18 +130,16 @@ add_library(mephit SHARED
   ${MEPHIT_C_SRC_FILES}
   ${MEPHIT_CPP_SRC_FILES}
 )
-add_dependencies(mephit TRIANGLE)
 target_include_directories(mephit PUBLIC
   ${magfie_include_dir}
-  ${TRIANGLE_INCLUDE_DIR}
   ${FFTW_INCLUDE_DIRS}
 )
 set(MEPHIT_LIBS
   GSL::gsl
   BLAS::BLAS
   LAPACK::LAPACK
   ${FFTW_LIBRARIES}
-  ${triangle_lib}
+  fortfem_triangle_compat
   LIBNEO::magfie
   LIBNEO::odeint
   LIBNEO::hdf5_tools

README.md

@@ -12,7 +12,7 @@ Prerequisites from external sources for running MEPHIT are as follows.
 - [SuiteSparse](https://github.com/DrTimothyAldenDavis/SuiteSparse)
 - [GSL](https://www.gnu.org/software/gsl/)
 - [FFTW3](http://fftw.org/)
-- [Triangle](https://www.cs.cmu.edu/~quake/triangle.html)
+- [fortfem](https://github.com/lazy-fortran/fortfem), fetched automatically; its MIT-licensed mesher replaces Shewchuk's Triangle
 - [Boost](https://www.boost.org/)
 - [FreeFem++](https://github.com/FreeFem/FreeFem-sources)
 - [MFEM](https://mfem.org/) is optional

cmake/SetupCODE.cmake

@@ -1,9 +1,8 @@
-# Dependency overrides for a hermetic build. With nothing set, libneo, Triangle,
+# Dependency overrides for a hermetic build. With nothing set, libneo, fortfem,
 # and (optionally) MFEM are fetched and built from source.
 #
 #   -DLIBNEO_REF=<ref>   libneo branch, tag, or SHA to fetch
 #   -DLIBNEO_PATH=<dir>  local libneo source directory (empty = fetch from git)
-#   -DTRIANGLE_DIR=<dir>  prebuilt Triangle directory (see SetupTriangle)
 #   -DMFEM_DIR=<dir>      prebuilt MFEM directory (see SetupMFEM)
 set(LIBNEO_REF "f3f241dea2c7dc25cf29731bb0a1328ff98f48c8" CACHE STRING
     "libneo branch, tag, or SHA to fetch")

cmake/SetupFortfem.cmake

@@ -0,0 +1,20 @@
+# Fetch fortfem for its Triangle-compatible mesher (target
+# fortfem_triangle_compat, MIT-licensed). This is the current fortfem main
+# commit as of 2026-08-17, pinned so the MEPHIT branch remains reproducible.
+include(FetchContent)
+
+set(FORTFEM_REF "605dc7f056aa9b15cf8e0673eaa13b13ce12a273" CACHE STRING
+    "fortfem commit to fetch")
+set(FORTFEM_PATH "" CACHE PATH
+    "local fortfem source directory (empty = fetch from git)")
+set(FORTFEM_TRIANGLE_COMPAT_ONLY ON CACHE BOOL "" FORCE)
+
+if(FORTFEM_PATH)
+    add_subdirectory(${FORTFEM_PATH} ${CMAKE_CURRENT_BINARY_DIR}/fortfem EXCLUDE_FROM_ALL)
+else()
+    FetchContent_Declare(fortfem
+        GIT_REPOSITORY https://github.com/lazy-fortran/fortfem.git
+        GIT_TAG ${FORTFEM_REF}
+    )
+    FetchContent_MakeAvailable(fortfem)
+endif()

src/mephit_fem.c

@@ -8,7 +8,6 @@
 #include <sys/stat.h>
 #include <gsl/gsl_errno.h>
 #include <gsl/gsl_integration.h>
-#include "triangle.h"
 #include "mephit_util.h"
 #include "mephit_fem.h"
 
@@ -243,86 +242,3 @@ void gauss_legendre_unit_interval(int order, double *points, double *weights)
   gsl_integration_glfixed_table_free(table);
 }
 
-void FEM_triangulate_external(const int npt_inner,
-                              const int npt_outer,
-                              const double *bdry_R,
-                              const double *bdry_Z,
-                              const double R_mid,
-                              const double Z_mid,
-                              const char *fname)
-{
-  int k;
-  FILE *fid;
-  struct triangulateio in, out, vorout;
-
-  // initialize all fields to zero or NULL
-  memset(&in, 0, sizeof(struct triangulateio));
-  memset(&out, 0, sizeof(struct triangulateio));
-  memset(&vorout, 0, sizeof(struct triangulateio));
-
-  in.numberofpoints = npt_inner + npt_outer;
-  in.numberofsegments = npt_inner + npt_outer;
-  in.numberofholes = 1;
-  in.pointlist = (REAL *) calloc(2 * (size_t) in.numberofpoints, sizeof(REAL));
-  in.pointmarkerlist = (int *) calloc((size_t) in.numberofpoints, sizeof(int));
-  in.segmentlist = (int *) calloc(2 * (size_t) in.numberofsegments, sizeof(int));
-  in.holelist = (REAL *) calloc(2 * (size_t) in.numberofholes, sizeof(REAL));
-  for (k = 0; k < in.numberofpoints; ++k) {
-    in.pointlist[2 * k] = bdry_R[k];
-    in.pointlist[2 * k + 1] = bdry_Z[k];
-  }
-  for (k = 0; k < npt_inner; ++k) {
-    in.segmentlist[2 * k] = k;
-    in.segmentlist[2 * k + 1] = (k + 1) % npt_inner;
-  }
-  for (k = 0; k < npt_outer; ++k) {
-    in.segmentlist[2 * (npt_inner + k)] = npt_inner + k;
-    in.segmentlist[2 * (npt_inner + k) + 1] = npt_inner + (k + 1) % npt_outer;
-  }
-  in.holelist[0] = R_mid;
-  in.holelist[1] = Z_mid;
-
-  // triangulate options:
-  // B - omit boundary markers in output
-  // e - generate edge list (to be used later)
-  // j - clean point list
-  // n - generate neighbor list (to be used later)
-  // p - triangulate from given polygon boundary
-  // q - minimum angle of 20 degrees
-  // Y - don't modify boundary edges
-  // z - use zero indexing
-  triangulate("BejnpqYz", &in, &out, &vorout);
-
-  fid = fopen(fname, "w");
-  fprintf(fid, "%i %i %i\n",
-          out.numberofpoints,
-          out.numberoftriangles,
-          out.numberofsegments);
-  for (k = 0; k < out.numberofpoints; ++k) {
-    fprintf(fid, "%.16e %.16e 0\n",
-            out.pointlist[2 * k],
-            out.pointlist[2 * k + 1]);
-  }
-  for (k = 0; k < out.numberoftriangles; ++k) {
-    fprintf(fid, "%i %i %i 1\n",
-            out.trianglelist[3 * k] + 1,
-            out.trianglelist[3 * k + 1] + 1,
-            out.trianglelist[3 * k + 2] + 1);
-  }
-  for (k = 0; k < out.numberofsegments; ++k) {
-    fprintf(fid, "%i %i 2\n",
-            out.segmentlist[2 * k] + 1,
-            out.segmentlist[2 * k + 1] + 1);
-  }
-  fclose(fid);
-
-  free(in.pointlist);
-  free(in.pointmarkerlist);
-  free(in.segmentlist);
-  free(in.holelist);  // same as out.holelist
-  trifree(out.pointlist);
-  trifree(out.trianglelist);
-  trifree(out.neighborlist);
-  trifree(out.segmentlist);
-  trifree(out.edgelist);
-}

src/mephit_fem.h

@@ -34,13 +34,6 @@ void FEM_deinit(void);
 
 void gauss_legendre_unit_interval(int order, double *points, double *weights);
 
-void FEM_triangulate_external(const int npt_inner,
-                              const int npt_outer,
-                              const double *bdry_R,
-                              const double *bdry_Z,
-                              const double R_mid,
-                              const double Z_mid,
-                              const char *fname);
 void Rtree_init(int ntri, double *tri_bb);
 void Rtree_query(double R, double Z, int *result_size, int **results);
 

src/mephit_mesh.F90

@@ -259,15 +259,6 @@ subroutine gauss_legendre_unit_interval(order, points, weights) &
       real(c_double), intent(out), dimension(1:order) :: points, weights
     end subroutine gauss_legendre_unit_interval
 
-    subroutine FEM_triangulate_external(npt_inner, npt_outer, node_R, node_Z, R_O, Z_O, fname) &
-      bind(C, name = 'FEM_triangulate_external')
-      use iso_c_binding, only: c_char, c_int, c_double
-      integer(c_int), intent(in), value :: npt_inner, npt_outer
-      real(c_double), intent(in), dimension(1:npt_inner + npt_outer) :: node_R, node_Z
-      real(c_double), intent(in), value :: R_O, Z_O
-      character(c_char), intent(in) :: fname(*)
-    end subroutine FEM_triangulate_external
-
     subroutine Rtree_init(ntri, tri_bb) bind(C, name = 'Rtree_init')
       use iso_c_binding, only: c_int, c_double
       integer(c_int), intent(in), value :: ntri
@@ -3052,7 +3043,7 @@ end subroutine mesh_write_MFEM
 #endif
 
   subroutine write_FreeFem_mesh
-    use iso_c_binding, only: c_null_char
+    use mephit_triangulate, only: FEM_triangulate_external
     use mephit_util, only: pi, linspace
     use mephit_conf, only: basename_suffix, decorate_filename
     integer :: fid, kpoi, ktri, kp, kedge, npt_inner, npt_outer
@@ -3109,7 +3100,7 @@ subroutine write_FreeFem_mesh
     bdry_R(npt_inner+1:) = R_mid + R_rad * cos(theta)
     bdry_Z(npt_inner+1:) = Z_mid + Z_rad * sin(theta)
     call FEM_triangulate_external(npt_inner, npt_outer, bdry_R, bdry_Z, R_mid, Z_mid, &
-      decorate_filename('outer.msh', '', basename_suffix) // c_null_char)
+      decorate_filename('outer.msh', '', basename_suffix))
     deallocate(bdry_R, bdry_Z, theta)
   end subroutine write_FreeFem_mesh
 

src/mephit_triangulate.f90

@@ -0,0 +1,69 @@
+module mephit_triangulate
+  use iso_fortran_env, only: dp => real64
+  implicit none
+  private
+  public :: FEM_triangulate_external
+
+contains
+
+  !> Triangulate the region between two closed boundary loops and write the
+  !> result in FreeFem msh format. Replaces Shewchuk's Triangle with flags
+  !> "BejnpqYz" by fortfem's clean-room reimplementation.
+  subroutine FEM_triangulate_external(npt_inner, npt_outer, bdry_R, bdry_Z, &
+    R_hole, Z_hole, fname)
+    use triangle_compat, only: tc_result_t, triangulate_compat
+    integer, intent(in) :: npt_inner, npt_outer
+    real(dp), intent(in) :: bdry_R(npt_inner + npt_outer), bdry_Z(npt_inner + npt_outer)
+    real(dp), intent(in) :: R_hole, Z_hole
+    character(len = *), intent(in) :: fname
+    type(tc_result_t) :: res
+    real(dp) :: points(2, npt_inner + npt_outer), holes(2, 1)
+    integer :: segments(2, npt_inner + npt_outer)
+    integer :: fid, k, stat
+
+    points(1, :) = bdry_R
+    points(2, :) = bdry_Z
+    do k = 1, npt_inner
+      segments(:, k) = [k, mod(k, npt_inner) + 1]
+    end do
+    do k = 1, npt_outer
+      segments(:, npt_inner + k) = [npt_inner + k, npt_inner + mod(k, npt_outer) + 1]
+    end do
+    holes(:, 1) = [R_hole, Z_hole]
+    call triangulate_compat(points, segments, holes, res, stat, &
+      min_angle = 20d0, quality = .true., nobisect = 1)
+    if (stat /= 0) then
+      error stop 'FEM_triangulate_external: triangulate_compat failed'
+    end if
+    open(newunit = fid, file = fname, status = 'replace', form = 'formatted', &
+      action = 'write')
+    write (fid, '(i0, 2(1x, i0))') res%npoints, res%ntriangles, res%nsegments
+    do k = 1, res%npoints
+      write (fid, '(a, 1x, a, " 0")') &
+        c_double_string(res%points(1, k)), c_double_string(res%points(2, k))
+    end do
+    do k = 1, res%ntriangles
+      write (fid, '(i0, 2(1x, i0), " 1")') res%triangles(:, k)
+    end do
+    do k = 1, res%nsegments
+      write (fid, '(i0, 1x, i0, " 2")') res%segments(:, k)
+    end do
+    close(fid)
+  end subroutine FEM_triangulate_external
+
+  !> Render a double like C's printf("%.16e", ...) to keep the msh file
+  !> format byte-identical to the former C writer.
+  function c_double_string(val) result(str)
+    real(dp), intent(in) :: val
+    character(len = :), allocatable :: str
+    character(len = 24) :: buf
+    integer :: k
+
+    write (buf, '(es23.16e2)') val
+    str = trim(adjustl(buf))
+    do k = 1, len(str)
+      if (str(k:k) == 'E') str(k:k) = 'e'
+    end do
+  end function c_double_string
+
+end module mephit_triangulate