diff --git a/CMakeLists.txt b/CMakeLists.txt index b2a2f9f..f1355a3 100644 --- a/CMakeLists.txt +++ b/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,10 +130,8 @@ 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 @@ -140,7 +139,7 @@ set(MEPHIT_LIBS BLAS::BLAS LAPACK::LAPACK ${FFTW_LIBRARIES} - ${triangle_lib} + fortfem_triangle_compat LIBNEO::magfie LIBNEO::odeint LIBNEO::hdf5_tools diff --git a/README.md b/README.md index d397fa0..69632f8 100644 --- a/README.md +++ b/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 diff --git a/cmake/SetupCODE.cmake b/cmake/SetupCODE.cmake index 91c9e4a..d87f073 100644 --- a/cmake/SetupCODE.cmake +++ b/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= libneo branch, tag, or SHA to fetch # -DLIBNEO_PATH= local libneo source directory (empty = fetch from git) -# -DTRIANGLE_DIR= prebuilt Triangle directory (see SetupTriangle) # -DMFEM_DIR= prebuilt MFEM directory (see SetupMFEM) set(LIBNEO_REF "f3f241dea2c7dc25cf29731bb0a1328ff98f48c8" CACHE STRING "libneo branch, tag, or SHA to fetch") diff --git a/cmake/SetupFortfem.cmake b/cmake/SetupFortfem.cmake new file mode 100644 index 0000000..3f593dd --- /dev/null +++ b/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() diff --git a/cmake/SetupTriangle.cmake b/cmake/SetupTriangle.cmake deleted file mode 100644 index 25122f3..0000000 --- a/cmake/SetupTriangle.cmake +++ /dev/null @@ -1,28 +0,0 @@ -if(NOT "${TRIANGLE_DIR}" STREQUAL "") - add_custom_target(TRIANGLE ALL) -elseif(DEFINED ENV{TRIANGLE_DIR}) - add_custom_target(TRIANGLE ALL) - set(TRIANGLE_DIR $ENV{TRIANGLE_DIR}) -else()# Define Triangle repository - set(TRIANGLE_REPO "https://salsa.debian.org/science-team/triangle.git") - set(TRIANGLE_PATCHES_DIR "debian/patches") - set(TRIANGLE_PATCH_SERIES "${TRIANGLE_PATCHES_DIR}/series") - - # External Project to fetch, patch, and build Triangle - ExternalProject_Add(TRIANGLE - PREFIX ${CMAKE_BINARY_DIR}/triangle - GIT_REPOSITORY ${TRIANGLE_REPO} - GIT_TAG master - CONFIGURE_COMMAND "" # Skip configure step - BUILD_COMMAND /bin/bash ${CMAKE_SOURCE_DIR}/cmake/setup_triangle.sh - INSTALL_COMMAND "" - DOWNLOAD_EXTRACT_TIMESTAMP TRUE - BUILD_BYPRODUCTS /libtriangle.so - ) - set(TRIANGLE_DIR ${CMAKE_BINARY_DIR}/triangle/src/TRIANGLE) -endif() - -set(TRIANGLE_INCLUDE_DIR ${TRIANGLE_DIR} CACHE STRING "Path to triangle include") -set(TRIANGLE_LIB_DIR ${TRIANGLE_DIR} CACHE STRING "Path to triangle lib") - -set(triangle_lib ${TRIANGLE_LIB_DIR}/libtriangle.so CACHE STRING "TRIANGLE lib") diff --git a/cmake/setup_triangle.sh b/cmake/setup_triangle.sh deleted file mode 100644 index 1621879..0000000 --- a/cmake/setup_triangle.sh +++ /dev/null @@ -1,34 +0,0 @@ -#!/bin/bash -set -e - -cd "$1" - -CC=${CC:-gcc} - -if [ "$(uname)" == "Darwin" ]; then - CFLAGS="-I/opt/homebrew/include" - LDFLAGS="-L/opt/homebrew/lib" -fi - - -if [ -f is_patched ]; then - echo "Triangle already patched." -else - echo "Patching Triangle..." - while read -r patch; do - patch -p1 < "debian/patches/$patch" - done < debian/patches/series - touch is_patched -fi - -if [ ! -f libtriangle-1.6.so ]; then - echo "Building Triangle..." - $CC triangle.c -o triangle.o $CFLAGS -O2 -DTRILIBRARY -fPIC -DPIC -c - $CC -shared triangle.o -o libtriangle-1.6.so $LDFLAGS -lm -fi - -if [ ! -f libtriangle.so ]; then - ln -s libtriangle-1.6.so libtriangle.so -fi - -echo "Triangle built successfully." diff --git a/src/mephit_fem.c b/src/mephit_fem.c index e6829db..24728d6 100644 --- a/src/mephit_fem.c +++ b/src/mephit_fem.c @@ -8,7 +8,6 @@ #include #include #include -#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); -} diff --git a/src/mephit_fem.h b/src/mephit_fem.h index 2b7ed74..e711d36 100644 --- a/src/mephit_fem.h +++ b/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); diff --git a/src/mephit_mesh.F90 b/src/mephit_mesh.F90 index 152e31e..9192d14 100644 --- a/src/mephit_mesh.F90 +++ b/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 diff --git a/src/mephit_triangulate.f90 b/src/mephit_triangulate.f90 new file mode 100644 index 0000000..2073357 --- /dev/null +++ b/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