diff --git a/Makefile.in b/Makefile.in index c94f92fa99..d7b7714594 100755 --- a/Makefile.in +++ b/Makefile.in @@ -104,6 +104,9 @@ endif ifneq (@HAVE_PETSC33@,yes) EXCLUDE_TAGS := $(EXCLUDE_TAGS) -e petsc33 endif +ifneq (@HAVE_LIBGMSH@,yes) + EXCLUDE_TAGS := $(EXCLUDE_TAGS) -e libgmsh +endif .SUFFIXES: .f90 .F90 .c .cpp .o .a diff --git a/configure b/configure index 7b4a3d2283..bcfa8c8d66 100755 --- a/configure +++ b/configure @@ -672,6 +672,7 @@ INSTALL_SCRIPT INSTALL_PROGRAM UNROLL_LOOPS OPENMP_CXXFLAGS +HAVE_LIBGMSH LAPACK_LIBS BLAS_LIBS HAVE_CGAL @@ -9928,6 +9929,119 @@ if test "$have_udunits" = "yes"; then #define HAVE_LIBUDUNITS 1 EOF fi +{ $as_echo "$as_me:${as_lineno-$LINENO}: checking Checking for libGmsh support" >&5 +$as_echo_n "checking Checking for libGmsh support... " >&6; } +ac_ext=cpp +ac_cpp='$CXXCPP $CPPFLAGS' +ac_compile='$CXX -c $CXXFLAGS $CPPFLAGS conftest.$ac_ext >&5' +ac_link='$CXX -o conftest$ac_exeext $CXXFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' +ac_compiler_gnu=$ac_cv_cxx_compiler_gnu + +LIBS_bck="$LIBS" +LIBS="-lGmsh" +cat confdefs.h - <<_ACEOF >conftest.$ac_ext +/* end confdefs.h. */ + + #include + +#ifdef F77_DUMMY_MAIN + +# ifdef __cplusplus + extern "C" +# endif + int F77_DUMMY_MAIN() { return 1; } + +#endif +#ifdef FC_DUMMY_MAIN +#ifndef FC_DUMMY_MAIN_EQ_F77 +# ifdef __cplusplus + extern "C" +# endif + int FC_DUMMY_MAIN() { return 1; } +#endif +#endif +int +main () +{ + + GmshInitialize (); + + ; + return 0; +} +_ACEOF +if ac_fn_cxx_try_link "$LINENO"; then : + + { $as_echo "$as_me:${as_lineno-$LINENO}: result: yes" >&5 +$as_echo "yes" >&6; } + LIBS="$LIBS $LIBS_bck" + $as_echo "#define HAVE_LIBGMSH 1" >>confdefs.h + + HAVE_LIBGMSH="yes" + +else + + LIBS="-lgmsh -lgl2ps" + cat confdefs.h - <<_ACEOF >conftest.$ac_ext +/* end confdefs.h. */ + + #include +#ifdef F77_DUMMY_MAIN + +# ifdef __cplusplus + extern "C" +# endif + int F77_DUMMY_MAIN() { return 1; } + +#endif +#ifdef FC_DUMMY_MAIN +#ifndef FC_DUMMY_MAIN_EQ_F77 +# ifdef __cplusplus + extern "C" +# endif + int FC_DUMMY_MAIN() { return 1; } +#endif +#endif +int +main () +{ + + GmshInitialize (); + + ; + return 0; +} +_ACEOF +if ac_fn_cxx_try_link "$LINENO"; then : + + { $as_echo "$as_me:${as_lineno-$LINENO}: result: yes" >&5 +$as_echo "yes" >&6; } + LIBS="$LIBS $LIBS_bck" + $as_echo "#define HAVE_LIBGMSH 1" >>confdefs.h + + HAVE_LIBGMSH="yes" + +else + + HAVE_LIBGMSH="no" + { $as_echo "$as_me:${as_lineno-$LINENO}: result: no" >&5 +$as_echo "no" >&6; } + LIBS="$LIBS_bck" + +fi +rm -f core conftest.err conftest.$ac_objext \ + conftest$ac_exeext conftest.$ac_ext + +fi +rm -f core conftest.err conftest.$ac_objext \ + conftest$ac_exeext conftest.$ac_ext + +ac_ext=c +ac_cpp='$CPP $CPPFLAGS' +ac_compile='$CC -c $CFLAGS $CPPFLAGS conftest.$ac_ext >&5' +ac_link='$CC -o conftest$ac_exeext $CFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' +ac_compiler_gnu=$ac_cv_c_compiler_gnu + # Check whether --enable-openmp was given. if test "${enable_openmp+set}" = set; then : diff --git a/configure.in b/configure.in index aec487022f..d0a40fd68e 100644 --- a/configure.in +++ b/configure.in @@ -305,6 +305,37 @@ if test "$have_udunits" = "yes"; then #define HAVE_LIBUDUNITS 1 EOF fi +AC_MSG_CHECKING([Checking for libGmsh support]) +AC_LANG_PUSH([C++]) +LIBS_bck="$LIBS" +LIBS="-lGmsh" +AC_TRY_LINK([ + #include + ],[ + GmshInitialize (); + ],[ + AC_MSG_RESULT([yes]) + LIBS="$LIBS $LIBS_bck" + AC_DEFINE(HAVE_LIBGMSH, [1]) + HAVE_LIBGMSH="yes" + ],[ + LIBS="-lgmsh -lgl2ps" + AC_TRY_LINK([ + #include ],[ + GmshInitialize (); + ],[ + AC_MSG_RESULT([yes]) + LIBS="$LIBS $LIBS_bck" + AC_DEFINE(HAVE_LIBGMSH, [1]) + HAVE_LIBGMSH="yes" + ],[ + HAVE_LIBGMSH="no" + AC_MSG_RESULT([no]) + LIBS="$LIBS_bck" + ]) + ]) +AC_SUBST(HAVE_LIBGMSH) +AC_LANG_POP([C++]) AC_ARG_ENABLE(openmp, [AC_HELP_STRING([--enable-openmp], diff --git a/diagnostics/Metric_Diagnostics.F90 b/diagnostics/Metric_Diagnostics.F90 index ec027c28b3..61bcdc1f92 100644 --- a/diagnostics/Metric_Diagnostics.F90 +++ b/diagnostics/Metric_Diagnostics.F90 @@ -31,6 +31,7 @@ module metric_diagnostics use fldebug use quicksort use spud + use global_parameters, only: timestep use vector_tools use metric_tools use fields @@ -40,14 +41,18 @@ module metric_diagnostics use diagnostic_source_fields use edge_length_module use field_derivatives + use merge_tensors use form_metric_field + use metric_assemble + use simple_diagnostics implicit none private public :: calculate_scalar_edge_lengths, calculate_field_tolerance, & - & calculate_eigenvalues_symmetric + & calculate_eigenvalues_symmetric, calculate_interpolation_metric, & + calculate_minimum_metric contains @@ -130,4 +135,36 @@ subroutine calculate_eigenvalues_symmetric(state, v_field) end subroutine calculate_eigenvalues_symmetric + subroutine calculate_interpolation_metric(states, t_field) + type(state_type), dimension(:), intent(inout) :: states + type(tensor_field), intent(inout) :: t_field + + call assemble_metric(states, t_field) + + end subroutine calculate_interpolation_metric + + subroutine calculate_minimum_metric(states, t_field) + type(state_type), intent(inout), dimension(:) :: states + type(tensor_field), intent(inout) :: t_field + + type(tensor_field) :: metric + logical, save :: initialised = .false. + + if (timestep==0) then + call initialise_diagnostic_from_checkpoint(t_field) + initialised = .true. + return + end if + + if (initialised) then + call allocate(metric, t_field%mesh, "ErrorMetric") + call assemble_metric(states, metric) + call merge_tensor_fields(t_field, metric) + call deallocate(metric) + else + call assemble_metric(states, t_field) + initialised = .true. + end if + end subroutine calculate_minimum_metric + end module metric_diagnostics diff --git a/diagnostics/Simple_Diagnostics.F90 b/diagnostics/Simple_Diagnostics.F90 index f4901341e9..b51c8a9ef0 100644 --- a/diagnostics/Simple_Diagnostics.F90 +++ b/diagnostics/Simple_Diagnostics.F90 @@ -54,8 +54,10 @@ module simple_diagnostics public :: calculate_temporalmax_scalar, calculate_temporalmax_vector, calculate_temporalmin, calculate_l2norm, & calculate_time_averaged_scalar, calculate_time_averaged_vector, & + calculate_time_averaged_tensor, & calculate_time_averaged_scalar_squared, & - calculate_time_averaged_vector_times_scalar, calculate_period_averaged_scalar + calculate_time_averaged_vector_times_scalar, calculate_period_averaged_scalar, & + initialise_diagnostic_from_checkpoint ! for the period_averaged_scalar routine real, save :: last_output_time @@ -314,6 +316,39 @@ subroutine calculate_time_averaged_vector(state, v_field) end if end subroutine calculate_time_averaged_vector + subroutine calculate_time_averaged_tensor(state, t_field) + type(state_type), intent(in) :: state + type(tensor_field), intent(inout) :: t_field + + type(tensor_field), pointer :: source_field + real :: a, b, spin_up_time, current_time, dt + integer :: stat + logical :: absolute_vals + + if (timestep==0) then + call initialise_diagnostic_from_checkpoint(t_field) + return + end if + + call get_option("/timestepping/current_time", current_time) + call get_option("/timestepping/timestep", dt) + + absolute_vals=have_option(trim(t_field%option_path)//"/diagnostic/algorithm/absolute_values") + call get_option(trim(t_field%option_path)//"/diagnostic/algorithm/spin_up_time", spin_up_time, stat) + if (stat /=0) spin_up_time=0. + source_field => tensor_source_field(state, t_field) + if(absolute_vals) source_field%val = abs(source_field%val) + + if (current_time>spin_up_time) then + a = (current_time-spin_up_time-dt)/(current_time-spin_up_time); b = dt/(current_time-spin_up_time) + ! v_field = a*v_field + b*source_field + call scale(t_field, a) + call addto(t_field, source_field, b) + else + call set(t_field, source_field) + end if + end subroutine calculate_time_averaged_tensor + subroutine calculate_time_averaged_scalar_squared(state, s_field) type(state_type), intent(in) :: state type(scalar_field), intent(inout) :: s_field diff --git a/femtools/Field_Options.F90 b/femtools/Field_Options.F90 index ee1df71725..ea0a9efca7 100644 --- a/femtools/Field_Options.F90 +++ b/femtools/Field_Options.F90 @@ -98,7 +98,7 @@ module field_options & extract_pressure_mesh, extract_velocity_mesh, & & postprocess_periodic_mesh, get_diagnostic_coordinate_field, & & get_nodal_coordinate_field, extract_prognostic_pressure, & - & extract_prognostic_velocity + & extract_prognostic_velocity, get_surface_ids integer, parameter, public :: FIELD_EQUATION_UNKNOWN = 0, & FIELD_EQUATION_ADVECTIONDIFFUSION = 1, & @@ -1337,4 +1337,29 @@ function extract_prognostic_velocity(state, stat) result(vfield) end function extract_prognostic_velocity + subroutine get_surface_ids(bc_path, mesh, surface_ids) + character(len=*), intent(in) :: bc_path + type(mesh_type) :: mesh + integer, dimension(:), allocatable, intent(out) :: surface_ids + integer shape_option(2) + character(len = OPTION_PATH_LEN) :: surface_names + character(len = FIELD_NAME_LEN), dimension(:), allocatable :: name_list + + !! subroutine reads the set of integral surface ids from a + !! boundary condition path in an options file, performing + !! a lookup of named boundaries if necessary + + if (have_option(bc_path//"/surface_ids")) then + shape_option=option_shape(bc_path//"/surface_ids") + allocate(surface_ids(1:shape_option(1))) + call get_option(bc_path//"/surface_ids", surface_ids) + else if (have_option(bc_path//"/surface_names")) then + call get_option(bc_path//"/surface_names", surface_names) + call tokenize(trim(surface_names), name_list, ",") + allocate(surface_ids(size(name_list))) + call get_surface_ids_from_names(mesh, name_list, surface_ids) + end if + + end subroutine get_surface_ids + end module field_options diff --git a/femtools/Fields_Allocates.F90 b/femtools/Fields_Allocates.F90 index 1cb0a259ac..84040f4c8a 100644 --- a/femtools/Fields_Allocates.F90 +++ b/femtools/Fields_Allocates.F90 @@ -538,6 +538,11 @@ subroutine deallocate_mesh(mesh) deallocate(mesh%element_halos) end if + if (associated(mesh%surface_names)) then + deallocate(mesh%surface_names) + nullify(mesh%surface_names) + end if + call deallocate_faces(mesh) if(associated(mesh%subdomain_mesh)) then @@ -1097,6 +1102,10 @@ function make_mesh (model, shape, continuity, name) & size(mesh%faces%surface_node_list), name='Surface'//trim(mesh%name)) #endif end if + if (associated(model%surface_names)) then + allocate(mesh%surface_names(size(model%surface_names))) + mesh%surface_names = model%surface_names + end if call addref(mesh) end function make_mesh diff --git a/femtools/Fields_Base.F90 b/femtools/Fields_Base.F90 index 0b0b5f3ffb..f6c444168f 100644 --- a/femtools/Fields_Base.F90 +++ b/femtools/Fields_Base.F90 @@ -367,7 +367,8 @@ module fields_base extract_scalar_field_from_tensor_field, ele_curl_at_quad,& eval_shape, ele_jacobian_at_quad, ele_div_at_quad_tensor,& ele_2d_curl_at_quad, getsndgln, local_coords_matrix,& - local_coords_interpolation + local_coords_interpolation, set_surface_names, & + get_surface_ids_from_names contains @@ -4199,4 +4200,49 @@ subroutine write_minmax_tensor(tfield, field_expression) end subroutine write_minmax_tensor + subroutine set_surface_names(mesh, names, ids) + type(mesh_type), intent(inout) :: mesh + character(len=*), intent(in) :: names(:) + integer, intent(in) :: ids(:) + + integer :: i + + !! This function allocates memory and stores in a mesh a map from + !! an array of boundary names to an array of integral boundary ids + + assert(size(names) == size(ids)) + + allocate(mesh%surface_names(size(ids))) + + do i=1, size(ids) + mesh%surface_names(i)%id = ids(i) + mesh%surface_names(i)%name = names(i) + end do + end subroutine set_surface_names + + subroutine get_surface_ids_from_names(mesh, names, ids) + type(mesh_type), intent(in) :: mesh + character(len=*), intent(in) :: names(:) + integer, intent(out) :: ids(:) + + integer :: i, j + + !! This function performs a look-up of the integral boundary + !! stored in a mesh corresponding to a list of names + + assert(size(names) == size(ids)) + + ids = -1 + + do i=1, size(ids) + do j=1, size(mesh%surface_names) + if (trim(mesh%surface_names(j)%name) == trim(names(i))) then + ids(i) = mesh%surface_names(j)%id + exit + end if + end do + end do + + end subroutine get_surface_ids_from_names + end module fields_base diff --git a/femtools/Fields_Data_Types.F90 b/femtools/Fields_Data_Types.F90 index befa147a9e..952319afba 100644 --- a/femtools/Fields_Data_Types.F90 +++ b/femtools/Fields_Data_Types.F90 @@ -50,6 +50,11 @@ module fields_data_types integer, public, parameter :: FIELD_TYPE_NORMAL=0, FIELD_TYPE_CONSTANT=1, FIELD_TYPE_PYTHON=2, & FIELD_TYPE_DEFERRED=3 + type surface_name + integer :: id + character (len=FIELD_NAME_LEN) :: name + end type surface_name + type adjacency_cache type(csr_sparsity), pointer :: nnlist => null() type(csr_sparsity), pointer :: nelist => null() @@ -89,6 +94,8 @@ module fields_data_types !! A list of ids marking different parts of the mesh !! so that initial conditions can be associated with it. integer, dimension(:), pointer :: region_ids=>null() + !! list mapping surface names to ids + type(surface_name), dimension(:), pointer :: surface_names=>null() !! Halo information for parallel simulations. type(halo_type), dimension(:), pointer :: halos=>null() type(halo_type), dimension(:), pointer :: element_halos=>null() diff --git a/femtools/Futils.F90 b/femtools/Futils.F90 index 093ad838bc..2447adc673 100644 --- a/femtools/Futils.F90 +++ b/femtools/Futils.F90 @@ -31,6 +31,7 @@ module futils !!< Some generic fortran utility functions. use fldebug + use iso_c_binding use global_parameters, only : real_digits_10 implicit none @@ -72,13 +73,23 @@ module futils #endif end type integer_vector + interface + function strlen(s) bind(c, name='strlen') + use iso_c_binding, only: c_ptr, c_size_t + implicit none + type(c_ptr), intent(in), value :: s + integer(c_size_t) :: strlen + end function strlen + end interface + private public :: real_format_len, real_format, nullify, real_vector, real_matrix,& integer_vector, int2str, present_and_true, present_and_false, present_and_zero,& present_and_nonzero, present_and_nonempty, free_unit, nth_digit, count_chars,& multiindex, file_extension_len, file_extension, trim_file_extension_len,& - trim_file_extension, random_number_minmax, int2str_len, starts_with, tokenize + trim_file_extension, random_number_minmax, int2str_len, starts_with, tokenize,& + copy_c_string_to_fortran contains @@ -435,5 +446,40 @@ subroutine tokenize(string, tokens, delimiter) assert(start_index == len(string) + 1 + len(delimiter)) end subroutine tokenize + + subroutine copy_c_string_to_fortran(c_string, f_string) + + type(c_ptr), intent(in), target :: c_string + character(len=*), intent(out) :: f_string + + character(kind=c_char), pointer :: arr(:) + character(:,kind = c_char), pointer :: f_ptr + integer(c_size_t) :: c_len + + c_len = strlen(c_string) + +!!! associate the array with the c string. + call c_f_pointer(c_string, arr, [c_len]) +!!! make into a pointer to a Fortran string + call wrap_to_string(c_len, arr, f_ptr) + if (len(f_string) .ge. c_len) then +!!! copy happens here + f_string = f_ptr + else +!!! Warn about truncation + f_string = f_ptr(1:len(f_string)) + ewrite(-1,*) "Copy from C string truncated to "//f_string + end if + + contains + + subroutine wrap_to_string(array_len, char_array, f_string) + integer(c_size_t), intent(in) :: array_len + character(kind=c_char, len=array_len), intent(in), target :: char_array(1) + character(:,c_char), pointer :: f_string + f_string => char_array(1) + end subroutine wrap_to_string + + end subroutine copy_c_string_to_fortran end module futils diff --git a/femtools/GMSH_Common.F90 b/femtools/GMSH_Common.F90 index 64201ec32f..09b64be1cf 100644 --- a/femtools/GMSH_Common.F90 +++ b/femtools/GMSH_Common.F90 @@ -29,9 +29,17 @@ ! This module contains code and variables common to all the GMSH I/O routines +#include "fdebug.h" module gmsh_common + use iso_c_binding + use fields, only: vector_field, tensor_field, face_ele, ele_loc,& + face_loc, mesh_dim, has_discontinuous_internal_boundaries,& + node_count, element_count, getsndgln, unique_surface_element_count + + implicit none + character(len=3), parameter :: GMSHVersionStr = "2.1" integer, parameter :: asciiFormat = 0 integer, parameter :: binaryFormat = 1 @@ -58,6 +66,151 @@ module gmsh_common integer, pointer :: tags(:), nodeIDs(:) end type GMSHelement + interface + subroutine cgmsh_initialise() bind(c) + end subroutine cgmsh_initialise + end interface + + interface + subroutine cgmsh_finalise(gmodel) bind(c) + use iso_c_binding + type(c_ptr) :: gmodel + end subroutine cgmsh_finalise + end interface + + interface + subroutine cread_gmsh_file(gmodel, filename) bind(c) + use iso_c_binding + type(c_ptr) :: gmodel + character(c_char) :: filename(*) + end subroutine cread_gmsh_file + end interface + + interface + subroutine cmesh_gmsh_file(gmodel, view_name) bind(c) + use iso_c_binding + type(c_ptr) :: gmodel + character(c_char) :: view_name(*) + end subroutine cmesh_gmsh_file + end interface + + interface + subroutine cread_gmsh_sizes(gmodel, numNodes, numFaces, numElements,& + haveRegionIDs, haveBounds, haveElementOwners, & + haveColumns, dim, loc, sloc) bind(c) + use iso_c_binding + type(c_ptr) :: gmodel + integer(c_int) :: numNodes, numFaces, numElements, dim, loc, sloc + logical(c_bool) :: haveRegionIDs, haveBounds, & + haveElementOwners, haveColumns + end subroutine cread_gmsh_sizes + end interface + + interface + subroutine cread_gmsh_element_connectivity(gmodel, numElements, loc,& + ndglno, regionIDs) bind(c) + use iso_c_binding + type(c_ptr) :: gmodel + integer(c_int) :: numElements, loc + integer(c_int) :: ndglno(numElements*loc), regionIDs(numElements) + end subroutine cread_gmsh_element_connectivity + end interface + + interface + subroutine cread_gmsh_points(gmodel, dim, numNodes, val) bind(c) + use iso_c_binding + type(c_ptr) :: gmodel + integer(c_int) :: dim, numNodes + real(c_double) :: val(*) + end subroutine cread_gmsh_points + end interface + + interface + subroutine cread_gmsh_face_connectivity(gmodel, numFaces, & + sloc, sndglno, & + haveBounds, boundaryIDs, & + haveElementOwners, faceOwner) bind(c) + use iso_c_binding + type(c_ptr) :: gmodel + integer(c_int) :: numFaces, sloc + logical(c_bool) :: haveBounds, haveElementOwners + integer(c_int) :: sndglno(*), boundaryIDs(*), faceOwner(*) + end subroutine cread_gmsh_face_connectivity + end interface + + interface + function cgmsh_count_physical_names(gm, dim) bind(c) + use iso_c_binding + type(c_ptr), intent(in) :: gm + integer(c_int) :: dim + integer (c_int) :: cgmsh_count_physical_names + end function cgmsh_count_physical_names + end interface + + interface + function cget_gmsh_physical_name(gm, it, dim, idx, c_string) bind(c) + use iso_c_binding + type(c_ptr), intent(in) :: gm + type(c_ptr), intent(inout) :: it + integer(c_int) :: dim, idx + type (c_ptr), intent(out) :: c_string + logical(c_bool) :: cget_gmsh_physical_name + end function cget_gmsh_physical_name + end interface + + interface + subroutine cread_gmsh_node_data(gmodel, name, data, step) bind(c) + use iso_c_binding + type(c_ptr), intent(in) :: gmodel + character(c_char), intent(in) :: name(*) + real(c_double), intent(out) :: data(*) + integer(c_int), intent(in) :: step + end subroutine cread_gmsh_node_data + end interface + + interface + subroutine cmesh_to_gmodel(gmodel, numNodes,& + numElements, numFaces, loc, sloc, gdim, pdim, val,& + etype, eles, ftype, faces, ele_ids, face_ids, ele_owners) bind(c) + use iso_c_binding + type(c_ptr) :: gmodel + integer(c_int) :: numNodes, numElements, numFaces,& + loc, sloc, gdim, pdim, etype, ftype + real(c_double) :: val(gdim*numNodes) + integer(c_int) :: eles(numElements*loc), faces(numFaces*sloc) + type(c_ptr), value :: ele_ids, face_ids, ele_owners + end subroutine cmesh_to_gmodel + end interface + + interface + subroutine cwrite_gmsh_file(gmodel, binary, filename) bind(c) + use iso_c_binding + type(c_ptr) :: gmodel + logical(c_bool) :: binary + character(c_char) :: filename(*) + end subroutine cwrite_gmsh_file + end interface + + interface + subroutine cdata_to_pview_node_data(gm, pvdata, & + numNodes, data, & + name, numComponents) bind(c) + use iso_c_binding + type(c_ptr) :: gm, pvdata + integer(c_int) :: numNodes, numComponents + real(c_double) :: data(*) + character(c_char) :: name(*) + end subroutine cdata_to_pview_node_data + end interface + + interface + subroutine cwrite_gmsh_data_file(pvdata, binary, filename) bind(c) + use iso_c_binding + type(c_ptr) :: pvdata + logical(c_bool) :: binary + character(c_char) :: filename(*) + end subroutine cwrite_gmsh_data_file + end interface contains @@ -128,7 +281,7 @@ end subroutine binary_formatting subroutine toFluidityElementNodeOrdering( oldList, elemType ) integer, pointer :: oldList(:) integer, dimension(size(oldList)) :: nodeOrder, flNodeList - integer i, elemType + integer i, elemType, numNodes numNodes = size(oldList) @@ -162,7 +315,7 @@ end subroutine toFluidityElementNodeOrdering subroutine toGMSHElementNodeOrdering( oldList, elemType ) integer, pointer :: oldList(:) integer, dimension(size(oldList)) :: nodeOrder, gmshNodeList - integer i, elemType + integer i, elemType, numNodes numNodes = size(oldList) @@ -206,4 +359,118 @@ subroutine deallocateElementList( elements ) end subroutine deallocateElementList +#ifdef HAVE_LIBGMSH + + subroutine position_to_gmodel(positions, gmodel) + type(vector_field), intent(in) :: positions + type(c_ptr), intent(out) :: gmodel + + integer :: numNodes, numElements, numFaces, sloc, & + loc, dim, pdim, i + integer, allocatable, dimension(:) :: sndglno + integer, allocatable, dimension(:), target ::owners + logical needs_element_owners + + type(c_ptr) :: bnd_ids, reg_ids, ele_owners + + numNodes = node_count(positions) + numElements = element_count(positions) + numFaces = unique_surface_element_count(positions%mesh) + needs_element_owners = has_discontinuous_internal_boundaries(positions%mesh) + + dim = mesh_dim(positions) + pdim = size(positions%val,1) + loc = ele_loc(positions, 1) + sloc = 1 + if (numFaces>0) sloc = face_loc(positions, 1) + + allocate(sndglno(numFaces*sloc)) + if (numFaces>0) call getsndgln(positions%mesh, sndglno) + + if (associated(positions%mesh%region_ids)) then + reg_ids = c_loc(positions%mesh%region_ids(1)) + else + reg_ids = c_null_ptr + end if + bnd_ids = c_null_ptr + if (numFaces>0) then + if (associated(positions%mesh%faces%boundary_ids)) then + if (size(positions%mesh%faces%boundary_ids)>0) & + bnd_ids = c_loc(positions%mesh%faces%boundary_ids(1)) + end if + end if + if (needs_element_owners) then + allocate(owners(numFaces)) + do i=1, numFaces + owners(i) = face_ele(positions%mesh,i) + end do + ele_owners = c_loc(owners(1)) + else + ele_owners = c_null_ptr + end if + + call cmesh_to_gmodel(gmodel, numNodes,& + numElements, numFaces, loc, sloc,& + dim, pdim, positions%val, & + gmsh_type(loc, dim), positions%mesh%ndglno,& + gmsh_type(sloc, dim-1), sndglno, reg_ids,& + bnd_ids, ele_owners) + + deallocate(sndglno) + + contains + + function gmsh_type(loc, dim) + + integer, intent(in) ::loc, dim + integer gmsh_type + + + if (loc .eq. dim+1) then + select case(dim) + case(0) + gmsh_type = 15 + case(1) + gmsh_type = 1 + case(2) + gmsh_type = 2 + case(3) + gmsh_type = 4 + end select + else + select case(dim) + case(2) + gmsh_type = 3 + case(3) + gmsh_type = 5 + end select + end if + + end function gmsh_type + + end subroutine position_to_gmodel + + subroutine tensor_field_to_pview(gmodel, tfield) + type(c_ptr) :: gmodel + type(tensor_field) :: tfield + + type(c_ptr) :: pvdata + + real, dimension(3,3,node_count(tfield)) :: data + integer :: dim + + data=0 + data(1,1,:)=1.0 + data(2,2,:)=1.0 + data(3,3,:)=1.0 + dim = size(tfield%val,1) + data(1:dim,1:dim,:) = tfield%val + + call cdata_to_pview_node_data(gmodel, pvdata, & + node_count(tfield), data, trim(tfield%name)//c_null_char,9) + + end subroutine tensor_field_to_pview + +#endif + end module gmsh_common diff --git a/femtools/LibGmsh_integration.cpp b/femtools/LibGmsh_integration.cpp new file mode 100644 index 0000000000..abc9df08bb --- /dev/null +++ b/femtools/LibGmsh_integration.cpp @@ -0,0 +1,471 @@ +/* Copyright (C) 2006 Imperial College London and others. + + Please see the AUTHORS file in the main source directory for a full list + of copyright holders. + + Prof. C Pain + Applied Modelling and Computation Group + Department of Earth Science and Engineering + Imperial College London + + amcgsoftware@imperial.ac.uk + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation, + version 2.1 of the License. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 + USA +*/ +#include "fmangle.h" + +#ifdef HAVE_LIBGMSH +#include "gmsh/Context.h" +#include "gmsh/Field.h" +#include "gmsh/Gmsh.h" +#include "gmsh/GModel.h" +#include "gmsh/GEntity.h" +#include "gmsh/MElement.h" +#include "gmsh/MVertex.h" +#include "gmsh/PView.h" +#include "gmsh/PViewDataGModel.h" + +#include +#include +#include + +typedef std::vector EntityVector; +typedef std::vector VertexVector; + +// Length set to match FIELD_NAME_LEN +char buffer[102]=""; + +int FluidityNodeOrdering(int type, int node) { + // This routine maps fluidity ids to (linear) gmsh type vertex ids. + int quads[]= {0,1,3,2}; + int hexes[]= {0,1,3,2,4,5,7,6}; + switch(type) { + case TYPE_QUA: + return quads[node]; + case TYPE_HEX: + return hexes[node]; + } + return node; +} + + +int GmshNodeOrdering(int type, int node) { + // This routine maps specific gmsh type vertex ids to fluidity ids + int quads[]= {0,1,3,2}; + int hexes[]= {0,1,3,2,4,5,7,6}; + switch(type) { + case MSH_QUA_4: + return quads[node]; + case MSH_HEX_8: + return hexes[node]; + } + return node; +} +extern "C" { + + void cgmsh_initialise() { + // Routine (over) initialises Gmsh functionality + GmshInitialize(); + } + + void cgmsh_finalise(GModel *&gm) { + // Routine cleans up Gmsh functionality + delete gm; + GmshFinalize(); + } + + void cread_gmsh_file(GModel *&gm, const char* filename) { + // Fortran accessible interface to fead a Gmsh .msh file + // known to gmsh into a GModel object + gm = new GModel; + // This switch swaps how Gmsh uses the first two integer tags on elements + CTX::instance()->mesh.switchElementTags=1; + // Use generic load, rather than readXXX to support multiple formats + gm->load(filename); + // This call reads any data beyond the mesh topology + PView::readMSH(filename); + } + + void cmesh_gmsh_file(GModel *&gm, char* name) { + // Load the file "name" into the GModel pointer provided + if (strlen(name)>0) { + PView* pv=PView::getViewByName(name); + if (pv) { + // This call assigns the PView "name" as the background mesh + gm->getFields()->setBackgroundMesh(pv->getIndex()); + // use BAMG in 2d + CTX::instance()->mesh.algo2d = 7; + CTX::instance()->mesh.algo3d = 7; + } + } + // This causes the meshing to happen + gm->mesh(gm->getDim()); + CTX::instance()->mesh.switchElementTags=0; + } + + void cread_gmsh_sizes(GModel *&gm, int &numNodes, int &numFaces, + int &numElements, bool &haveRegionIDs, + bool &haveBounds, bool &haveElementOwners, + bool &haveColumns, + int &gdim, int &loc, int &sloc) { + // This utility routine provides Fortran with the sizes to allocate + // memory for. + + // Get the largest dimensionality represented + gdim = gm->getDim(); + + bool physicals = CTX::instance()->mesh.switchElementTags==1; + + numElements = 0; + numFaces = 0; + haveRegionIDs = false; + haveBounds = false; + haveElementOwners = false; + haveColumns = not (PView::getViewByName("column_ids") == NULL); + EntityVector ents; + EntityVector eles, all_eles; + EntityVector faces; + gm->getEntities(ents); + for (EntityVector::iterator it = ents.begin(); it != ents.end(); ++it) { + if ((*it)->dim() == gdim) { + if((*it)->physicals.size()>0) { + eles.push_back(*it); + } else { + all_eles.push_back(*it); + } + } + if ((*it)->dim() == gdim-1 && (physicals || (*it)->physicals.size()>0)) + faces.push_back(*it); + } // Number of mesh nodes in model + + if (eles.size()==0 && all_eles.size()==1) { + eles.push_back(all_eles[0]); + } + if (eles.size()>0) { + loc = eles[0]->getMeshElement(0)->getNumVertices(); + } else { + loc = 0; + } + if (faces.size()>0) { + sloc = faces[0]->getMeshElement(0)->getNumVertices(); + } else { + sloc = 0; + } + + for (EntityVector::iterator it = eles.begin(); it != eles.end(); ++it) { + numElements = numElements + (*it)->getNumMeshElements(); + if (physicals) { + haveRegionIDs = haveRegionIDs || (*it)->tag()>0; + } else { + haveRegionIDs = haveRegionIDs || (*it)->getPhysicalEntities().size()>0; + } + } + + if(not haveRegionIDs || physicals) { + eles[0]->addPhysicalEntity(0); + } + + for (EntityVector::iterator it = faces.begin(); it != faces.end(); ++it) { + if ((*it)->getNumMeshElements()) + haveElementOwners = haveElementOwners || + (*it)->getMeshElement(0)->getPartition()>0; + } + for (EntityVector::iterator it = faces.begin(); it != faces.end(); ++it) { + numFaces = numFaces + (*it)->getNumMeshElements(); + haveBounds = haveBounds || (*it)->tag()>0; + if (physicals) (*it)->addPhysicalEntity(0); + } + + if (CTX::instance()->mesh.switchElementTags==0) { + numNodes = gm->indexMeshVertices(false, 0, true); + } else { + numNodes = gm->indexMeshVertices(false, 0, false); + } + + + } + + bool cmp_MElement( const MElement* e1, const MElement* e2) { + return e1->getNum() < e2->getNum(); + } + + void cread_gmsh_element_connectivity(GModel *&gm, int &numElements, + int &loc, int *ndglno, int *regionIDs) { + + std::map ele_label; + std::list ele_list; + + int dim = gm->getDim(); + EntityVector ents; + gm->getEntities(ents); + for (EntityVector::iterator it = ents.begin(); it != ents.end(); ++it) { + if ((*it)->dim() == dim) { + for (unsigned int i=0; i< (*it)->getNumMeshElements(); ++i) { + MElement* e = (*it)->getMeshElement(i); + ele_list.push_back(e); + ele_label[e] = (*it)->tag(); + } + } + } + + ele_list.sort(cmp_MElement); + + int k=0; + for (std::list::iterator e = ele_list.begin(); + e != ele_list.end(); ++e) { + for (int j=0; j<(*e)->getNumVertices(); ++j) { + MVertex* v = (*e)->getVertex(FluidityNodeOrdering((*e)->getType(),j)); + ndglno[loc*k+j] = v->getIndex(); + } + if (regionIDs) { + int tag = ele_label[*e]; + if (tag>0) { + regionIDs[k] = tag; + } + } + ++k; + } + } + + void cread_gmsh_points(GModel *&gm, int &dim, int &numNodes, double *val) { + EntityVector ents; + gm->getEntities(ents); + for (EntityVector::iterator it = ents.begin(); it != ents.end(); ++it) { + for (unsigned int i=0; i<(*it)->mesh_vertices.size(); ++i) { + MVertex *v = (*it)->mesh_vertices[i]; + int idx = v->getIndex(); + if (idx>-1) { + switch(dim) { + case 3: + val[dim*(idx-1)+2] = v->z(); + case 2: + val[dim*(idx-1)+1] = v->y(); + case 1: + val[dim*(idx-1)+0] = v->x(); + } + } + } + } + } + + void cread_gmsh_face_connectivity(GModel *&gm, int &numFaces, + int &sloc, int *sndglno, + bool &haveBounds, int *boundaryIDs, + bool &haveElementOwners, int *faceOwner) { + + EntityVector ents; + EntityVector faces; + gm->getEntities(ents); + + bool physicals = CTX::instance()->mesh.switchElementTags==1; + + int dim = gm->getDim(); + + std::map face_label; + std::list face_list; + + for (EntityVector::iterator it = ents.begin(); it != ents.end(); ++it) { + if ((*it)->physicals.size()>0 && (*it)->dim() == dim-1) { + for (unsigned int i=0; i< (*it)->getNumMeshElements(); ++i) { + MElement* e = (*it)->getMeshElement(i); + face_list.push_back(e); + if (physicals) { + face_label[e] = (*it)->tag(); + } else { + face_label[e] = (*it)->getPhysicalEntities()[0]; + } + } + } + } + + face_list.sort(cmp_MElement); + + unsigned int k=0; + for (std::list::iterator e = face_list.begin(); + e != face_list.end(); ++e) { + for (int j=0; j<(*e)->getNumVertices(); ++j) { + MVertex* v = (*e)->getVertex(FluidityNodeOrdering((*e)->getType(),j)); + sndglno[sloc*k+j] = v->getIndex(); + } + if (haveBounds) { + int tag = face_label[*e]; + if (tag>0) { + boundaryIDs[k] = tag; + } else { + boundaryIDs[k] = 0; + } + } + if (haveElementOwners) { + if((*e)->getPartition()>0) { + faceOwner[k] = (*e)->getPartition(); + } + } + ++k; + } + } + + void cmesh_to_gmodel(GModel *&gm, + int &numNodes, int &numElements, int &numFaces, + int &loc, int &sloc, + int &gdim, int &pdim, double* val, + int &etype, int* eles, + int &ftype, int* faces, + int *ele_ids, int *face_ids, int *eleOwners) { + + std::map vertexMap; + std::map vertexInverseMap; + std::vector elementNum; + std::vector > vertexIndices; + std::vector elementType; + std::vector physical; + std::vector elementary; + std::vector partition; + + + double x[3] = {0,0,0}; + for (int n=0;nsetIndex(n+1); + vertexMap[n+1] = v; + vertexInverseMap[v->getNum()] = n+1; + } + + for (int i=0; i< numElements; ++i) { + elementNum.push_back(i+1); + std::vector vertices; + for (int j=0; j vertices; + for (int j=0; jgetEntities(ents); + + for (EntityVector::iterator it = ents.begin(); it != ents.end(); ++it) { + if ((*it)->dim()getDim()) continue; + for (unsigned int n=0;n<(*it)->getNumMeshElements();++n) { + MElement* e = (*it)->getMeshElement(n); + if (!e) continue; + for (int j=0; jgetNumVertices(); ++j) { + // This is necessary to get back the original node numbering + // (necessary for parallel runs). + MVertex *v = e->getVertex(j); + int idx = loc*(e->getNum()-1)+j; + v->setIndex(eles[idx]); + } + } + } + } + + void cread_gmsh_node_data(GModel *&gm, char* name, double* data, int &step) { + PView *pv = PView::getViewByName(name); + PViewDataGModel *d = dynamic_cast(pv->getData()); + int nComp = d->getNumComponents(step,0,0); + int nNodes = gm->getNumMeshVertices(); + for (int i=0; igetValueByIndex(step, i+1, 0 ,j, data[nComp*i+j]); + } + } + } + + void cdata_to_pview_node_data(GModel *&gm, PViewDataGModel *& pvd, + int &numNodes, + double *data, + char* name, + int &numComponents) { + + pvd = new PViewDataGModel(PViewDataGModel::NodeData); + std::map > vdata; + int k = 0; + for (int n=0;n lvec; + for (int i=0;iaddData( gm, vdata, 0, 0.0, 1, numComponents); + pvd->setName(name); + new PView(pvd); + } + + void cwrite_gmsh_file(GModel *&gm, bool &binary, char* filename) { + gm->writeMSH(filename, 2.2, binary); + } + + void cwrite_gmsh_data_file(PViewDataGModel *& pvd, bool &binary, + char* filename) { + pvd->writeMSH(filename, 2.2, binary, false, true); + } + + int cgmsh_count_physical_names(GModel *& gm, int &dim) { + + int count = 0; + for (GModel::piter it=gm->firstPhysicalName(); + it != gm->lastPhysicalName(); ++it) { + if (it->first.first == dim) ++count; + } + return count; + } + + bool cget_gmsh_physical_name(GModel *& gm, GModel::piter *&it, int &mdim, + int &idx, char* &c_string) { + if (it== NULL) { + it = new GModel::piter; + (*it) = gm->firstPhysicalName(); + } + while ((*it) != gm->lastPhysicalName() && (*it)->first.first !=mdim) { + ++(*it); + } + if ((*it) != gm->lastPhysicalName()) { + idx = (*it)->first.second; + buffer[(*it)->second.copy(buffer, 40)] = '\0'; + c_string = buffer; + ++(*it); + return true; + } else { + delete it; + return false; + } + } + +} + +#endif diff --git a/femtools/Makefile.in b/femtools/Makefile.in index 6cb7a2b82d..c934d7f140 100644 --- a/femtools/Makefile.in +++ b/femtools/Makefile.in @@ -107,7 +107,7 @@ OBJS = Dgtools.o Coordinates.o EventCounter.o \ GMSH_Common.o Read_GMSH.o Write_GMSH.o \ Exodusii_C_Interface.o Exodusii_F_Interface.o Exodusii_Common.o Read_Exodusii.o \ Mesh_Files.o Vertical_Extrapolation.o \ - Mesh_Quality.o Mesh_Quality_C.o + Mesh_Quality.o Mesh_Quality_C.o LibGmsh_integration.o # objects to be included in libfemtools: diff --git a/femtools/Mesh_Files.F90 b/femtools/Mesh_Files.F90 index 9ee67c88ff..e09323a47a 100644 --- a/femtools/Mesh_Files.F90 +++ b/femtools/Mesh_Files.F90 @@ -110,8 +110,10 @@ function read_mesh_simple(filename, format, quad_degree, & case("gmsh") field = read_gmsh_file(filename, quad_degree=quad_degree, quad_ngi=quad_ngi, & - quad_family=quad_family, mdim=mdim) - + quad_family=quad_family, mdim=mdim, format_string="msh") + case("geo") + field = read_gmsh_file(filename, quad_degree=quad_degree, quad_ngi=quad_ngi, & + quad_family=quad_family, mdim=mdim, format_string="geo") case("exodusii") #ifdef HAVE_LIBEXOIIV2C field = read_exodusii_file(filename, quad_degree=quad_degree, quad_ngi=quad_ngi, & diff --git a/femtools/Read_GMSH.F90 b/femtools/Read_GMSH.F90 index 9339802524..a851d41a47 100644 --- a/femtools/Read_GMSH.F90 +++ b/femtools/Read_GMSH.F90 @@ -31,9 +31,10 @@ module read_gmsh ! This module reads GMSH files and results in a vector field of ! positions. + use iso_c_binding use fldebug - use global_parameters, only : OPTION_PATH_LEN - use futils + use global_parameters, only : OPTION_PATH_LEN, FIELD_NAME_LEN + use futils, only : copy_c_string_to_fortran use quadrature use elements use spud @@ -60,12 +61,13 @@ module read_gmsh ! The main function for reading GMSH files function read_gmsh_simple( filename, quad_degree, & - quad_ngi, quad_family, mdim ) & + quad_ngi, quad_family, mdim, format_string, & + position, metric) & result (field) !!< Read a GMSH file into a coordinate field. !!< In parallel the filename must *not* include the process number. - character(len=*), intent(in) :: filename + character(len=*), intent(in) :: filename, format_string !! The degree of the quadrature. integer, intent(in), optional, target :: quad_degree !! The degree of the quadrature. @@ -74,6 +76,9 @@ function read_gmsh_simple( filename, quad_degree, & integer, intent(in), optional :: quad_family !! Dimension of mesh integer, intent(in), optional :: mdim + !! + type(vector_field), optional :: position + type(tensor_field), optional :: metric !! result: a coordinate field type(vector_field) :: field @@ -87,24 +92,153 @@ function read_gmsh_simple( filename, quad_degree, & character(len = parallel_filename_len(filename)) :: lfilename integer :: loc, sloc integer :: numNodes, numElements, numFaces - logical :: haveBounds, haveElementOwners, haveRegionIDs + logical(c_bool) :: haveBounds, haveElementOwners,& + haveRegionIDs, haveColumns integer :: dim, coordinate_dim, gdim integer :: gmshFormat integer :: n, d, e, f, nodeID +#ifdef HAVE_LIBGMSH + integer :: snames + type(c_ptr) :: gmodel, gm2, c_str, it + integer, allocatable, dimension(:) :: id_list + character(len = FIELD_NAME_LEN), dimension(:), allocatable :: name_list + real, dimension(:), allocatable :: lcolumn_ids +#endif + type(GMSHnode), pointer :: nodes(:) type(GMSHelement), pointer :: elements(:), faces(:) - + ! If running in parallel, add the process number if(isparallel()) then - lfilename = trim(parallel_filename(filename)) // ".msh" + lfilename = trim(parallel_filename(filename))& + // "." // format_string + else + lfilename = trim(filename) // "." // format_string + end if + +#if HAVE_LIBGMSH + + call cgmsh_initialise() + + call cread_gmsh_file(gmodel, trim(lfilename)//c_null_char) + + if (format_string == "geo") then + if (present(metric) ) then + call position_to_gmodel(position, gm2) + call tensor_field_to_pview(gm2, metric) + call cmesh_gmsh_file(gmodel, trim(metric%name)//c_null_char) + else + call cmesh_gmsh_file(gmodel,c_null_char) + end if + end if + + call cread_gmsh_sizes(gmodel, numNodes, numFaces, numElements,& + haveRegionIDs, haveBounds, haveElementOwners, & + haveColumns, dim, loc, sloc) + + if (present(mdim)) then + coordinate_dim = mdim + else if(have_option("/geometry/spherical_earth") ) then + ! on the n-sphere the input mesh may be 1/2d (extrusion), or 3d but + ! Coordinate is always geometry dimensional + call get_option('/geometry/dimension', gdim) + coordinate_dim = gdim + else + coordinate_dim = dim + end if + + if (present(quad_degree)) then + quad = make_quadrature(loc, dim, degree=quad_degree, family=quad_family) + else if (present(quad_ngi)) then + quad = make_quadrature(loc, dim, ngi=quad_ngi, family=quad_family) + else + FLAbort("Need to specify either quadrature degree or ngi") + end if + shape=make_element_shape(loc, dim, 1, quad) + call allocate(mesh, numNodes, numElements, shape, name="CoordinateMesh") + call allocate( field, coordinate_dim, mesh, name="Coordinate") + + snames = cgmsh_count_physical_names(gmodel, dim-1) + if (snames>0) then + allocate(name_list(snames+1), id_list(snames+1)) + n = 1 + it = c_null_ptr + do while (cget_gmsh_physical_name(gmodel,it, dim-1, id_list(n), c_str)) + call copy_c_string_to_fortran(c_str, name_list(n)) + n = n+1 + end do + call set_surface_names(field%mesh, name_list(1:snames), id_list(1:snames)) + end if + + ! deallocate our references of mesh, shape and quadrature: + call deallocate(mesh) + call deallocate(shape) + call deallocate(quad) + + if (haveRegionIDs) then + allocate( field%mesh%region_ids(numElements) ) + end if + if(haveColumns) then + allocate(field%mesh%columns(1:numNodes), lcolumn_ids(1:numNodes)) + end if + + call cread_gmsh_element_connectivity(gmodel, numElements, loc,& + field%mesh%ndglno, field%mesh%region_ids) + + call cread_gmsh_points(gmodel, coordinate_dim, numNodes, field%val) + + ! Now faces + allocate(sndglno(1:numFaces*sloc)) + sndglno=0 + if(haveBounds) then + allocate(boundaryIDs(1:numFaces)) + end if + if(haveElementOwners) then + allocate(faceOwner(1:numFaces)) + end if + + + call cread_gmsh_face_connectivity(gmodel, numFaces, sloc, & + sndglno, haveBounds, boundaryIDs, & + haveElementOwners, faceOwner) + + ! If we've got boundaries, do something + if( haveBounds ) then + if ( haveElementOwners ) then + call add_faces( field%mesh, & + sndgln = sndglno(1:numFaces*sloc), & + boundary_ids = boundaryIDs(1:numFaces), & + element_owner=faceOwner ) + else + call add_faces( field%mesh, & + sndgln = sndglno(1:numFaces*sloc), & + boundary_ids = boundaryIDs(1:numFaces) ) + end if else - lfilename = trim(filename) // ".msh" + ewrite(2,*) "WARNING: no boundaries in GMSH file "//trim(lfilename) + call add_faces( field%mesh, sndgln = sndglno(1:numFaces*sloc) ) end if + if (haveColumns) then + call cread_gmsh_node_data(gmodel,"column_ids"//c_null_char,& + lcolumn_ids, 0) + field%mesh%columns = lcolumn_ids + end if + + call cgmsh_finalise(gmodel) + + deallocate(sndglno) + if (haveBounds) deallocate(boundaryIDs) + if (haveElementOwners) deallocate(faceOwner) + +#else fd = free_unit() + if (format_string == "geo") & + FLAbort("Fluidity must be built with libgmsh support to read .geo files") + ! Open node file ewrite(2, *) "Opening "//trim(lfilename)//" for reading." open( unit=fd, file=trim(lfilename), err=43, action="read", & @@ -173,7 +307,16 @@ function read_gmsh_simple( filename, quad_degree, & haveRegionIDs = .false. end if - + + loc = size( elements(1)%nodeIDs ) + if (numFaces>0) then + sloc = size( faces(1)%nodeIDs ) + else + sloc = 0 + end if + + ! Now construct within Fluidity data structures + if (present(mdim)) then coordinate_dim = mdim else if(have_option("/geometry/spherical_earth") ) then @@ -185,15 +328,6 @@ function read_gmsh_simple( filename, quad_degree, & coordinate_dim = dim end if - loc = size( elements(1)%nodeIDs ) - if (numFaces>0) then - sloc = size( faces(1)%nodeIDs ) - else - sloc = 0 - end if - - ! Now construct within Fluidity data structures - if (present(quad_degree)) then quad = make_quadrature(loc, dim, degree=quad_degree, family=quad_family) else if (present(quad_ngi)) then @@ -280,6 +414,8 @@ function read_gmsh_simple( filename, quad_degree, & deallocate(faces) deallocate(elements) +#endif + return 43 FLExit("Unable to open "//trim(lfilename)) diff --git a/femtools/Write_GMSH.F90 b/femtools/Write_GMSH.F90 index 7c40156fc1..f030ab6547 100644 --- a/femtools/Write_GMSH.F90 +++ b/femtools/Write_GMSH.F90 @@ -30,7 +30,8 @@ module write_gmsh use fldebug - + + use iso_c_binding use global_parameters, only : OPTION_PATH_LEN use futils use elements @@ -51,7 +52,7 @@ module write_gmsh end interface ! Writes to GMSH binary format - can set to ASCII (handy for debugging) - logical, parameter :: useBinaryGMSH=.true. + logical(c_bool), parameter :: useBinaryGMSH=.true. contains @@ -82,8 +83,39 @@ subroutine write_mesh_to_gmsh(filename, state, mesh, number_of_partitions) end subroutine write_mesh_to_gmsh + ! ----------------------------------------------------------------- + +#ifdef HAVE_LIBGMSH + + subroutine write_libgmsh(filename, positions) + + character(len=*), intent(in):: filename + type(vector_field), intent(in):: positions + + character(len=longStringLen) :: meshFile + type(c_ptr) :: gmodel, pvdata + + call cgmsh_initialise() + + call position_to_gmodel(positions, gmodel) + + meshFile = trim(filename) // ".msh" + + call cwrite_gmsh_file(gmodel, useBinaryGMSH, trim(meshFile)//c_null_char) + + if (associated(positions%mesh%columns)) then + call cdata_to_pview_node_data(gmodel, pvdata, & + node_count(positions), real(positions%mesh%columns), & + "column_ids"//c_null_char,1) + call cwrite_gmsh_data_file(pvdata, useBinaryGMSH, & + trim(meshFile)//c_null_char) + end if + + call cgmsh_finalise(gmodel) + end subroutine write_libgmsh +#endif ! ----------------------------------------------------------------- @@ -105,6 +137,16 @@ subroutine write_positions_to_gmsh(filename, positions, number_of_partitions) else numParts = getnprocs() end if + +#ifdef HAVE_LIBGMSH + + if( getprocno() <= numParts ) then + + call write_libgmsh(filename, positions) + + end if + +#else ! Write out data only for those processes that contain data - SPMD requires ! that there be no early return @@ -135,7 +177,9 @@ subroutine write_positions_to_gmsh(filename, positions, number_of_partitions) ! Close GMSH file close( fileDesc ) - end if + end if + +#endif return @@ -151,7 +195,7 @@ end subroutine write_positions_to_gmsh subroutine write_gmsh_header( fd, lfilename, useBinaryGMSH ) integer :: fd character(len=*) :: lfilename - logical :: useBinaryGMSH + logical(c_bool) :: useBinaryGMSH character(len=999) :: GMSHVersionStr, GMSHFileFormat, GMSHdoubleNumBytes integer, parameter :: oneInt = 1 @@ -196,7 +240,7 @@ subroutine write_gmsh_nodes( fd, lfilename, field, useBinaryGMSH ) character(len=longStringLen) :: charBuf character(len=longStringLen) :: idStr, xStr, yStr, zStr type(vector_field), intent(in):: field - logical :: useBinaryGMSH + logical(c_bool) :: useBinaryGMSH integer numNodes, numDimen, numCoords, i, j real :: coords(3) @@ -254,7 +298,7 @@ subroutine write_gmsh_faces_and_elements( fd, lfilename, mesh, & useBinaryGMSH ) ! Writes out elements for the given mesh type(mesh_type), intent(in):: mesh - logical :: useBinaryGMSH + logical(c_bool) :: useBinaryGMSH character(len=*) :: lfilename character(len=longStringLen) :: charBuf @@ -364,17 +408,17 @@ subroutine write_gmsh_faces_and_elements( fd, lfilename, mesh, & case (2) if(useBinaryGMSH) then - write(fd, err=301) f, surface_element_id(mesh, f), 0, lnodelist + write(fd, err=301) f, surface_element_id(mesh, f), f, lnodelist else - write(fd, 6969, err=301) f, faceType, numTags, surface_element_id(mesh, f), 0, lnodelist + write(fd, 6969, err=301) f, faceType, numTags, surface_element_id(mesh, f), f, lnodelist end if case (4) if(useBinaryGMSH) then - write(fd, err=301) f, surface_element_id(mesh, f), 0, 0, face_ele(mesh, f), lnodelist + write(fd, err=301) f, surface_element_id(mesh, f), f, 0, face_ele(mesh, f), lnodelist else - write(fd, 6969, err=301) f, faceType, numTags, surface_element_id(mesh,f), 0, 0, & + write(fd, 6969, err=301) f, faceType, numTags, surface_element_id(mesh,f), f, 0, & face_ele(mesh,f), lnodelist end if @@ -401,19 +445,19 @@ subroutine write_gmsh_faces_and_elements( fd, lfilename, mesh, & if(associated(mesh%region_ids)) then if(useBinaryGMSH) then - write(fd, err=301) elemID, ele_region_id(mesh, e), 0, lnodelist + write(fd, err=301) elemID, ele_region_id(mesh, e), e, lnodelist else write(fd, 6969, err=301) elemID, elemType, 2, & - ele_region_id(mesh, e), 0, lnodelist + ele_region_id(mesh, e), e, lnodelist end if else if(useBinaryGMSH) then - write(fd, err=301) elemID, 0, 0, lnodelist + write(fd, err=301) elemID, 0, e, lnodelist else write(fd, 6969, err=301) elemID, elemType, 2, & - 0, 0, lnodelist + 0, e, lnodelist end if end if @@ -445,7 +489,7 @@ subroutine write_gmsh_node_columns( fd, meshFile, field, useBinaryGMSH ) integer :: fd character(len=*) :: meshFile type(vector_field), intent(in) :: field - logical :: useBinaryGMSH + logical(c_bool) :: useBinaryGMSH character(len=longStringLen) :: charBuf integer :: numNodes, timeStepNum, numComponents, i diff --git a/manual/configuring_fluidity.tex b/manual/configuring_fluidity.tex index 60f9bc52bf..c2eacef6b9 100644 --- a/manual/configuring_fluidity.tex +++ b/manual/configuring_fluidity.tex @@ -601,7 +601,10 @@ \subsection{Reading meshes from file} \lstinline[language=bash]+.msh+ when it runs. For information on generating meshes in Gmsh format, see chapter -\ref{chap:meshes}. +\ref{chap:meshes}. If compiled with full libgmsh support, then for serial runs +\fluidity can read the Gmsh geometry from a \lstinline[language=bash]+.geo+ file +directly by setting the \onlypdf\option{/geometry/mesh/from\_file/format}\ option +to \onlypdf\option{gmsh} \subsection{Deriving meshes from other meshes} \index{mesh!derived} diff --git a/manual/mesh_formats.tex b/manual/mesh_formats.tex index dfebaa2bea..8c5db1c1d7 100644 --- a/manual/mesh_formats.tex +++ b/manual/mesh_formats.tex @@ -50,7 +50,8 @@ \subsection{Surface IDs}\label{sec:surface_ids} faces in 1, 2 or 3 dimensions respectively) in order to set boundary conditions or other parameters. This number can then be used to specify which surface a particular boundary condition should be applied to -in \fluidity. +in \fluidity. If \fluidity is compiled with libgmsh support, then +surface names can also be used to label physical surfaces. \subsection{Region IDs}\label{sec:region_ids} \index{region ID} These are analogous to surface IDs, however they are diff --git a/manual/visualisation_and_diagnostics.tex b/manual/visualisation_and_diagnostics.tex index 31de084a58..c05696e502 100644 --- a/manual/visualisation_and_diagnostics.tex +++ b/manual/visualisation_and_diagnostics.tex @@ -443,7 +443,8 @@ \subsection{fltools} differentiate\_vtu & \ref{sec:differentiate_vtu} \\ edge\_length\_distribution & \ref{sec:edgelengthdist} \\ encode & \ref{sec:encode} \\ - fladapt & \ref{sec:fladapt} \\ + fladapt & \ref{sec:fladapt} \\ + flboundadapt & \ref{sec:flboundadapt} \\ fldecomp & \ref{sec:fldecomp} \\ fldiagnostics & \ref{sec:fldiagnostics} \\ flredecomp & \ref{sec:flredecomp} \\ @@ -703,6 +704,19 @@ \subsubsection{fladapt} \end{lstlisting} where \lstinline[language = Bash]+INPUT+ is the name of the input options file and \lstinline[language = Bash]+OUPUT+ is the name of the generated mesh. The flag \lstinline[language = Bash]+-v+ flag enables verbose output and the flag \lstinline[language = Bash]+-h+ flag displays the help message. +%%%%%%%%%%%%%%%%%%%%%%%%% FLBOUMDADAPT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsubsection{flboundadapt} +\label{sec:flboundadapt} +flbounadapt (only available if \fluidity is compiled with libgmsh support) performs a remeshing of an input geometry based on an input options file (which may be a checkpoint) and outputs the resulting mesh. It is run from the command line: + +\begin{lstlisting}[language = Bash] +flboundadapt [options] INPUT MESH OUTPUT +\end{lstlisting} +where \lstinline[language = Bash]+INPUT+ is the name of the input options file, \lstinline[language = Bash]+MESH+ is the name of the input Gmsh geometry (without the \lstinline[language = Bash]+.geo+ suffix) and \lstinline[language = Bash]+OUPUT+ is the name of the generated mesh. The flag \lstinline[language = Bash]+-v+ flag enables verbose output and the flag \lstinline[language = Bash]+-h+ flag displays the help message. + + + %%%%%%%%%%%%%%%%%%%%%%%%% FLDECOMP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsubsection{fldecomp} diff --git a/preprocessor/Boundary_Conditions_From_Options.F90 b/preprocessor/Boundary_Conditions_From_Options.F90 index 09ca3611e4..890d17a191 100644 --- a/preprocessor/Boundary_Conditions_From_Options.F90 +++ b/preprocessor/Boundary_Conditions_From_Options.F90 @@ -30,7 +30,7 @@ module boundary_conditions_from_options use fldebug use global_parameters, only: OPTION_PATH_LEN, PYTHON_FUNC_LEN, pi,& current_debug_level, FIELD_NAME_LEN - use futils, only: int2str, present_and_true + use futils, only: int2str, present_and_true, tokenize use vector_tools use quadrature use spud @@ -222,7 +222,7 @@ subroutine populate_scalar_boundary_conditions(field, bc_path, position, suppres character(len=OPTION_PATH_LEN) bc_path_i character(len=FIELD_NAME_LEN) bc_name, bc_type integer, dimension(:), allocatable:: surface_ids - integer i, nbcs, shape_option(2) + integer i, nbcs ! Get number of boundary conditions nbcs=option_count(trim(bc_path)) @@ -233,9 +233,7 @@ subroutine populate_scalar_boundary_conditions(field, bc_path, position, suppres bc_path_i=trim(bc_path)//"["//int2str(i)//"]" ! Get vector of surface ids - shape_option=option_shape(trim(bc_path_i)//"/surface_ids") - allocate(surface_ids(1:shape_option(1))) - call get_option(trim(bc_path_i)//"/surface_ids", surface_ids) + call get_surface_ids(trim(bc_path_i), field%mesh, surface_ids) ! Get name of boundary call get_option(trim(bc_path_i)//"/name", bc_name) @@ -375,15 +373,13 @@ subroutine populate_vector_boundary_conditions(state, field, bc_path, position, logical applies(3), have_sem_bc, have_smoothing, debugging_mode, prescribed(3) integer, dimension(:), allocatable:: surface_ids integer, dimension(:), pointer:: surface_element_list, surface_node_list - integer i, j, nbcs, shape_option(2) + integer i, j, nbcs nbcs=option_count(trim(bc_path)) boundary_conditions: do i=0, nbcs-1 bc_path_i=trim(bc_path)//"["//int2str(i)//"]" - shape_option=option_shape(trim(bc_path_i)//"/surface_ids") - allocate(surface_ids(1:shape_option(1))) - call get_option(trim(bc_path_i)//"/surface_ids", surface_ids) + call get_surface_ids(trim(bc_path_i), field%mesh, surface_ids) call get_option(trim(bc_path_i)//"/name", bc_name) call get_option(trim(bc_path_i)//"/type[0]/name", bc_type) diff --git a/preprocessor/Populate_State.F90 b/preprocessor/Populate_State.F90 index 1b91fddbc6..e0d3d73bca 100644 --- a/preprocessor/Populate_State.F90 +++ b/preprocessor/Populate_State.F90 @@ -268,7 +268,7 @@ subroutine insert_external_mesh(states, save_vtk_cache) if (is_active_process) then select case (mesh_file_format) - case ("triangle", "gmsh", "exodusii") + case ("triangle", "gmsh", "geo", "exodusii") ! Get mesh dimension if present call get_option(trim(mesh_path)//"/from_file/dimension", mdim, stat) ! Read mesh @@ -286,7 +286,8 @@ subroutine insert_external_mesh(states, save_vtk_cache) ! After successfully reading in an ExodusII mesh, change the option ! mesh file format to "gmsh", as the write routines for ExodusII are currently ! not implemented. Thus, checkpoints etc are dumped as gmsh mesh files - if (trim(mesh_file_format)=="exodusii") then + if (trim(mesh_file_format)=="exodusii" .or. & + trim(mesh_file_format)=="geo") then mesh_file_format = "gmsh" call set_option_attribute(trim(from_file_path)//"/format/name", trim(mesh_file_format), stat=stat) if (stat /= SPUD_NO_ERROR) then diff --git a/schemas/diagnostic_algorithms.rnc b/schemas/diagnostic_algorithms.rnc index ef24bc49e0..1d41cf7f6c 100644 --- a/schemas/diagnostic_algorithms.rnc +++ b/schemas/diagnostic_algorithms.rnc @@ -1070,6 +1070,24 @@ scalar_edge_lengths_algorithm = attribute material_phase_support { "single" } } ) +interpolation_metric_algorithm = + ( + ## give the metric which the adaptivity routine would generate + ## at the current time level + element algorithm { + attribute name { "interpolation_metric" }, + attribute material_phase_support { "multiple" } + } + ) +minimum_metric_algorithm = + ( + ## give the strictest metric which the adaptivity routine would generate + ## for the current simulation + element algorithm { + attribute name { "minimum_metric" }, + attribute material_phase_support { "multiple" } + } + ) field_tolerance_algorithm = ( ## From a field on a mesh, diagnose the anisotropic @@ -1528,6 +1546,25 @@ vector_scalar_time_averaged_algorithm = } ) +tensor_time_averaged_algorithm = + ( + ## Algorithm for time averaging tensor fields. + ## This field needs to be interpolated between adapts. + element algorithm { + attribute name { "time_averaged_tensor" }, + attribute material_phase_support { "single" }, + tensor_source_field, + ## When to start sampling? Defaults to simulation start time. + element spin_up_time { + real + }?, + ## If enabled, the average of the absolute value is taken + element absolute_values { + comment + }? + } + ) + # surface algorithm surface_horizontal_divergence_algorithm = @@ -1609,8 +1646,11 @@ vector_diagnostic_algorithms |= curl_algorithm vector_diagnostic_algorithms |= vector_python_diagnostic_algorithm tensor_diagnostic_algorithms = tensor_copy_algorithm +tensor_diagnostic_algorithms |= tensor_time_averaged_algorithm tensor_diagnostic_algorithms |= grad_vector_algorithm tensor_diagnostic_algorithms |= strain_rate_algorithm +tensor_diagnostic_algorithms |= minimum_metric_algorithm +tensor_diagnostic_algorithms |= interpolation_metric_algorithm tensor_diagnostic_algorithms |= hessian_algorithm tensor_diagnostic_algorithms |= helmholtz_smoothed_tensor_algorithm tensor_diagnostic_algorithms |= helmholtz_anisotropic_smoothed_tensor_algorithm diff --git a/schemas/diagnostic_algorithms.rng b/schemas/diagnostic_algorithms.rng index c4793dde19..0b13ade5fc 100644 --- a/schemas/diagnostic_algorithms.rng +++ b/schemas/diagnostic_algorithms.rng @@ -1430,6 +1430,30 @@ Water::Velocity + + + give the metric which the adaptivity routine would generate +at the current time level + + interpolation_metric + + + multiple + + + + + + give the strictest metric which the adaptivity routine would generate +for the current simulation + + minimum_metric + + + multiple + + + From a field on a mesh, diagnose the anisotropic @@ -2079,6 +2103,31 @@ This field needs to be interpolated between adapts. + + + Algorithm for time averaging tensor fields. +This field needs to be interpolated between adapts. + + time_averaged_tensor + + + single + + + + + When to start sampling? Defaults to simulation start time. + + + + + + If enabled, the average of the absolute value is taken + + + + + @@ -2282,12 +2331,21 @@ a horizontal vector living on the specified surface. + + + + + + + + + diff --git a/schemas/fluidity_options.rnc b/schemas/fluidity_options.rnc index 6e030a5b2a..027cf43ffd 100644 --- a/schemas/fluidity_options.rnc +++ b/schemas/fluidity_options.rnc @@ -2354,10 +2354,19 @@ prognostic_density_field = ## Boundary conditions element boundary_conditions { attribute name { string }, - ## Surface id: - element surface_ids { - integer_vector - }, + ( + ## Surface id: + element surface_ids { + integer_vector + } + | + ## Surface names + ## + ## Comma delimited list of Gmsh physical names + element surface_names { + anystring + } + ), ## Type ( element type { @@ -6824,10 +6833,19 @@ prognostic_photosynthetic_radiation = ## Boundary conditions element boundary_conditions { attribute name { string }, - ## Surface id: - element surface_ids { - integer_vector - }, + ( + ## Surface id: + element surface_ids { + integer_vector + } + | + ## Surface names + ## + ## Comma delimited list of Gmsh physical names + element surface_names { + anystring + } + ), ## Type ( element type { diff --git a/schemas/fluidity_options.rng b/schemas/fluidity_options.rng index 12cc2860cd..52c839a026 100644 --- a/schemas/fluidity_options.rng +++ b/schemas/fluidity_options.rng @@ -3181,10 +3181,18 @@ requires a distinct name for the options dictionary. - - Surface id: - - + + + Surface id: + + + + Surface names + +Comma delimited list of Gmsh physical names + + + Type @@ -8576,10 +8584,18 @@ radiation for water and phytoplankton. - - Surface id: - - + + + Surface id: + + + + Surface names + +Comma delimited list of Gmsh physical names + + + Type diff --git a/schemas/mesh_options.rnc b/schemas/mesh_options.rnc index 84437e9165..3c56c186d2 100644 --- a/schemas/mesh_options.rnc +++ b/schemas/mesh_options.rnc @@ -8,6 +8,11 @@ mesh_info = attribute name { "gmsh" }, comment }| + ## GMSH geometry format [requires libgmsh support] + element format { + attribute name { "geo" }, + comment + }| ## Read the mesh from a vtu. Note that the mesh will have no ## surface or region IDs. element format { diff --git a/schemas/mesh_options.rng b/schemas/mesh_options.rng index 195638ee82..6f7a45befc 100644 --- a/schemas/mesh_options.rng +++ b/schemas/mesh_options.rng @@ -12,6 +12,13 @@ + + GMSH geometry format [requires libgmsh support] + + geo + + + Read the mesh from a vtu. Note that the mesh will have no surface or region IDs. diff --git a/schemas/prognostic_field_options.rnc b/schemas/prognostic_field_options.rnc index 5a7b4bf8bb..7475f87b35 100644 --- a/schemas/prognostic_field_options.rnc +++ b/schemas/prognostic_field_options.rnc @@ -3,10 +3,19 @@ scalar_boundary_conditions = ## Boundary conditions element boundary_conditions { attribute name { string }, - ## Surface id: - element surface_ids { - integer_vector - }, + ( + ## Surface id: + element surface_ids { + integer_vector + } + | + ## Surface names + ## + ## Comma delimited list of Gmsh physical names + element surface_names { + anystring + } + ), ## Type ( element type { @@ -1186,10 +1195,19 @@ prognostic_velocity_field = ## Boundary conditions element boundary_conditions { attribute name { string }, - ## Surface id: - element surface_ids { - integer_vector - }, + ( + ## Surface id: + element surface_ids { + integer_vector + } + | + ## Surface names + ## + ## Comma delimited list of Gmsh physical names + element surface_names { + anystring + } + ), velocity_boundary_conditions }*, ## For a Newtonian fluid this is the shear viscosity. @@ -1566,10 +1584,19 @@ prognostic_pressure_field = )*, element boundary_conditions { attribute name { string }, - ## Surface id: - element surface_ids { - integer_vector - }, + ( + ## Surface id: + element surface_ids { + integer_vector + } + | + ## Surface names + ## + ## Comma delimited list of Gmsh physical names + element surface_names { + anystring + } + ), ## Type element type { attribute name { "dirichlet" }, @@ -1670,10 +1697,19 @@ prognostic_geostrophic_pressure_field = ## this only makes sense when excluding coriolis. element boundary_conditions { attribute name { string }, - ## Surface id: - element surface_ids { - integer_vector - }, + ( + ## Surface id: + element surface_ids { + integer_vector + } + | + ## Surface names + ## + ## Comma delimited list of Gmsh physical names + element surface_names { + anystring + } + ), ## Type element type { attribute name { "dirichlet" }, @@ -1702,10 +1738,19 @@ prognostic_vertical_balance_pressure_field = ## This is normally be a homogeneous bc on the top surface. element boundary_conditions { attribute name { string }, - ## Surface id: - element surface_ids { - integer_vector - }, + ( + ## Surface id: + element surface_ids { + integer_vector + } + | + ## Surface names + ## + ## Comma delimited list of Gmsh physical names + element surface_names { + anystring + } + ), ## Type element type { attribute name { "dirichlet" }, @@ -1832,10 +1877,19 @@ prognostic_foam_velocity_potential_field = element boundary_conditions { attribute replaces { "boundary, TTPER1 TTPER2 TTPERI" }, attribute name { string }, - ## Surface id: - element surface_ids { - integer_vector - }, + ( + ## Surface id: + element surface_ids { + integer_vector + } + | + ## Surface names + ## + ## Comma delimited list of Gmsh physical names + element surface_names { + anystring + } + ), ## Type ( element type { @@ -1916,10 +1970,19 @@ prognostic_multipath_stream_function_field = ## Boundary conditions element boundary_conditions { attribute name { string }, - ## Surface id: - element surface_ids { - integer_vector - }, + ( + ## Surface id: + element surface_ids { + integer_vector + } + | + ## Surface names + ## + ## Comma delimited list of Gmsh physical names + element surface_names { + anystring + } + ), ## The streamfunction will be zero on the primary boundary. There must be exactly one primary boundary. element primary_boundary{ empty @@ -1937,10 +2000,19 @@ prognostic_multipath_stream_function_field = ## Boundary conditions element boundary_conditions { attribute name { string }, - ## Surface id: - element surface_ids { - integer_vector - }, + ( + ## Surface id: + element surface_ids { + integer_vector + } + | + ## Surface names + ## + ## Comma delimited list of Gmsh physical names + element surface_names { + anystring + } + ), ## Secondary boundaries have a value given by the flux between this boundary and the primary boundary. element secondary_boundary{ ## A point on or behind the primary boundary *from* which the flux line should extend. diff --git a/schemas/prognostic_field_options.rng b/schemas/prognostic_field_options.rng index 012c30792f..ea489ff647 100644 --- a/schemas/prognostic_field_options.rng +++ b/schemas/prognostic_field_options.rng @@ -6,10 +6,18 @@ - - Surface id: - - + + + Surface id: + + + + Surface names + +Comma delimited list of Gmsh physical names + + + Type @@ -1447,10 +1455,18 @@ requires a distinct name for the options dictionary. - - Surface id: - - + + + Surface id: + + + + Surface names + +Comma delimited list of Gmsh physical names + + + @@ -1943,10 +1959,18 @@ requires a distinct name for the options dictionary. - - Surface id: - - + + + Surface id: + + + + Surface names + +Comma delimited list of Gmsh physical names + + + Type @@ -2077,10 +2101,18 @@ this only makes sense when excluding coriolis. - - Surface id: - - + + + Surface id: + + + + Surface names + +Comma delimited list of Gmsh physical names + + + Type @@ -2117,10 +2149,18 @@ This is normally be a homogeneous bc on the top surface. - - Surface id: - - + + + Surface id: + + + + Surface names + +Comma delimited list of Gmsh physical names + + + Type @@ -2266,10 +2306,18 @@ requires a distinct name for the options dictionary. - - Surface id: - - + + + Surface id: + + + + Surface names + +Comma delimited list of Gmsh physical names + + + Type @@ -2368,10 +2416,18 @@ Only applicable to control volume spatial discretisations. - - Surface id: - - + + + Surface id: + + + + Surface names + +Comma delimited list of Gmsh physical names + + + The streamfunction will be zero on the primary boundary. There must be exactly one primary boundary. @@ -2394,10 +2450,18 @@ Only applicable to control volume spatial discretisations. - - Surface id: - - + + + Surface id: + + + + Surface names + +Comma delimited list of Gmsh physical names + + + Secondary boundaries have a value given by the flux between this boundary and the primary boundary. diff --git a/tests/flboundadapt/circle.geo b/tests/flboundadapt/circle.geo new file mode 100644 index 0000000000..508357e766 --- /dev/null +++ b/tests/flboundadapt/circle.geo @@ -0,0 +1,40 @@ +Point(1) = {-1,-1,0}; +Point(2) = { 1,-1,0}; +Point(3) = { 1, 1,0}; +Point(4) = {-1, 1,0}; + +Point(5) = {0,0,0}; + +r=0.3; + +Point(6) = {-r, 0,0}; +Point(7) = { 0, r,0}; +Point(8) = { r, 0,0}; +Point(9) = { 0,-r,0}; + +Line(1) = {1,2}; +Line(2) = {2,3}; +Line(3) = {3,4}; +Line(4) = {4,1}; + +Circle(5) = {6,5,7}; +Circle(6) = {7,5,8}; +Circle(7) = {8,5,9}; +Circle(8) = {9,5,6}; + +Line Loop(1) = {1,2,3,4}; +Line Loop(2) = {5,6,7,8}; + +Plane Surface(1) = {1,2}; + +Physical Surface(1) = {1}; + +Physical Line(1) = {4}; +Physical Line(2) = {1,3}; +Physical Line(4) = {2}; +Physical Line(5) = {5,6,7,8}; + +Field[1]=MathEval; +Field[1].F="0.2"; + +Background Field =2; \ No newline at end of file diff --git a/tests/flboundadapt/flboundadapt.flml b/tests/flboundadapt/flboundadapt.flml new file mode 100644 index 0000000000..00a6edce64 --- /dev/null +++ b/tests/flboundadapt/flboundadapt.flml @@ -0,0 +1,133 @@ + + + + fladapt + + + fluids + + + + 2 + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2 + + + + + + vtk + + + + 1 + + + + + + + + 0.0 + + + 1.0 + + + 1.0 + + + + + + + + + + def val(x, t): + r = 1.0/(x[0]**2+x[1]**2) + return [r, 0.0] + + + + + + + + + + + + + + + + 1.0e-1 0.0 + + + + + + + + + + + + + + + + + + + + 1 + + + 10000 + + + + + + 1.0e-2 0.0 0.0 1.0e-2 + + + + + + + 1.0e0 0.0 0.0 1.0e0 + + + + + + diff --git a/tests/flboundadapt/flboundadapt.xml b/tests/flboundadapt/flboundadapt.xml new file mode 100644 index 0000000000..04a1f722a5 --- /dev/null +++ b/tests/flboundadapt/flboundadapt.xml @@ -0,0 +1,24 @@ + + + Offline geometry adaptivity + + flml 2dadapt libgmsh + + make clean-run-debug; flboundadapt -v flboundadapt circle output > flboundadapt.log + + + +import fluidity.diagnostics.gmshtools as gmshtools + +output = gmshtools.ReadMsh("output.msh") + + + + +import fluidity_tools +fluidity_tools.compare_variable(float(output.NodeCoordsCount()), 800.0, 0.05) + + + + + diff --git a/tests/named_boundary/named_boundary.flml b/tests/named_boundary/named_boundary.flml new file mode 100644 index 0000000000..0749ae2252 --- /dev/null +++ b/tests/named_boundary/named_boundary.flml @@ -0,0 +1,240 @@ + + + + named_boundary + + + fluids + + + + 2 + + + + + + + + + + + + + + discontinuous + + + + + + + + + + + + 2 + + + + + + + + + + 5 + + + + + + vtk + + + + 0 + + + + + + + + + 0.0 + + + 1.0 + + + 1.0 + + + 3 + + + + + + + + + + + + only first timestep + + + + + + + + 1.0e-7 + + + 1000 + + + + + + + + + right + + + + + 0.0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1.0 + + + + + 1.0 + + + 1.0 + + + + + + + 1.0e-7 + + + 1000 + + + + + + + + + 0.0 0.0 + + + + + left,right + + + + + + + 1.0 + + + + + 0.0 + + + + + + + + top,bottom + + + + + + + 0.0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/named_boundary/named_boundary.xml b/tests/named_boundary/named_boundary.xml new file mode 100644 index 0000000000..ea163d8033 --- /dev/null +++ b/tests/named_boundary/named_boundary.xml @@ -0,0 +1,51 @@ + + + + + named_boundary + + flml libgmsh + + + fluidity -v2 -l named_boundary.flml + + + + + +from fluidity_tools import stat_parser +s = stat_parser("named_boundary.stat") +time=s["ElapsedTime"]["value"] + + + +import vtktools + +file = vtktools.vtu('named_boundary_0.vtu') +file.GetFieldNames() +velocity = file.GetVectorField('Velocity') + + + + + + +import fluidity_tools +fluidity_tools.compare_variable(abs(time[-1]), 1.0, 1.0e-14) + + + +import fluidity_tools + +nnodes = len(velocity) +print 'max X error:', max(abs(velocity[:,0]- 1.0)) +print 'max Y error:', max(abs(velocity[:,1]- 0.0)) +assert(all(abs(velocity[:,0]- 1.0)< 1.0e-6)) +assert(all(abs(velocity[:,1])<1.0e-6)) + + + + + + + diff --git a/tests/named_boundary/src/square.geo b/tests/named_boundary/src/square.geo new file mode 100644 index 0000000000..143f9bde67 --- /dev/null +++ b/tests/named_boundary/src/square.geo @@ -0,0 +1,20 @@ +Point(1) = {0,0,0,0.1}; +Point(2) = {0,1,0,0.1}; +Point(3) = {1,1,0,0.1}; +Point(4) = {1,0,0,0.1}; + +Line(1) = {1,2}; +Line(2) = {2,3}; +Line(3) = {3,4}; +Line(4) = {4,1}; + +Line Loop(1) = {1,2,3,4}; + +Plane Surface(1) = 1; + +Physical Line("left") = 1; +Physical Line("top") = 2; +Physical Line("right") = 3; +Physical Line("bottom") = 4; + +Physical Surface("square") = 1; diff --git a/tools/Flboundadapt.F90 b/tools/Flboundadapt.F90 new file mode 100644 index 0000000000..d8dca2720a --- /dev/null +++ b/tools/Flboundadapt.F90 @@ -0,0 +1,186 @@ +! Copyright (C) 2006 Imperial College London and others. +! +! Please see the AUTHORS file in the main source directory for a full list +! of copyright holders. +! +! Prof. C Pain +! Applied Modelling and Computation Group +! Department of Earth Science and Engineering +! Imperial College London +! +! amcgsoftware@imperial.ac.uk +! +! This library is free software; you can redistribute it and/or +! modify it under the terms of the GNU Lesser General Public +! License as published by the Free Software Foundation, +! version 2.1 of the License. +! +! This library is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +! Lesser General Public License for more details. +! +! You should have received a copy of the GNU Lesser General Public +! License along with this library; if not, write to the Free Software +! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +! USA + +#include "fdebug.h" + +subroutine flboundadapt(input_flmlname_c, & + input_geometryname_c, & + output_meshname_c, metric_name_c) bind(c) + !!< Peforms a mesh adapt based on the supplied input options file. + !!< Outputs the resulting mesh. + + use iso_c_binding + use futils + use global_parameters, only: FIELD_NAME_LEN, OPTION_PATH_LEN,& + topology_mesh_name + use fldebug + use parallel_tools, only: isparallel + use reference_counting + use spud + use fields + use edge_length_module + use state_module + use field_options + use vtk_interfaces + use diagnostic_fields_wrapper, only: calculate_diagnostic_variables + use limit_metric_module + use mesh_files + use Read_gmsh + use populate_state_module + use metric_assemble + use adapt_state_module + use diagnostic_fields_new, only : & + & calculate_diagnostic_variables_new => calculate_diagnostic_variables + + implicit none + + interface + subroutine check_options() + end subroutine check_options + +#ifdef HAVE_PYTHON + subroutine python_init() + end subroutine python_init +#endif + end interface + + type(c_ptr), value :: input_flmlname_c, input_geometryname_c,& + output_meshname_c, metric_name_c + + character(len=FIELD_NAME_LEN):: input_flmlname + character(len=FIELD_NAME_LEN):: input_geometryname + character(len=FIELD_NAME_LEN):: output_meshname + character(len=FIELD_NAME_LEN):: metric_name + integer :: i + type(mesh_type), pointer :: old_mesh + type(state_type), dimension(:), pointer :: states + type(vector_field) :: new_mesh_field + type(vector_field), pointer :: position + type(tensor_field), pointer :: metric + type(tensor_field) :: t_edge_lengths + logical :: have_metric_name + + ! now turn into proper fortran string + call copy_c_string_to_fortran(input_flmlname_c, input_flmlname) + call copy_c_string_to_fortran(input_geometryname_c, input_geometryname) + call copy_c_string_to_fortran(output_meshname_c, output_meshname) + call copy_c_string_to_fortran(metric_name_c, metric_name) + + ewrite(1, *) "In flboundadapt" + +#ifdef HAVE_PYTHON + call python_init() +#endif + + have_metric_name = len(trim(metric_name)) > 0 + + ewrite(2, *) "Input file name: " // trim(input_flmlname) + ewrite(2, *) "Input geometry name: " // trim(input_geometryname) + ewrite(2, *) "Output mesh name: " // trim(output_meshname) + if (have_metric_name) then + ewrite(2, *) "Metric name: " // trim(output_meshname) + end if + + if (have_metric_name) then + allocate(states(1)) + call vtk_read_state(trim(input_flmlname)//".vtu", states(1)) + else + ! Load the options tree + call load_options(trim(input_flmlname) // ".flml") + if(.not. have_option("/simulation_name")) then + FLExit("Failed to find simulation name after loading options file") + end if + if(debug_level() >= 1) then + ewrite(1, *) "Options tree:" + call print_options() + end if +#ifdef DDEBUG + ewrite(1, *) "Performing options sanity check" + call check_options() + ewrite(1, *) "Options sanity check successful" +#endif + +! Populate the system state + call populate_state(states) + + ! Calculate diagnostic fields + call calculate_diagnostic_variables(states) + call calculate_diagnostic_variables_new(states) + end if + ! Find the external mesh field + + ! Note that Gmsh expects metric data on a mesh covering the geometry + ! so this may not work for periodic data. + + position => extract_vector_field(states(1), "Coordinate") + + if (have_metric_name) then + metric => extract_tensor_field(states(1), metric_name) + old_mesh => metric%mesh + else + ! Assemble the error metric + allocate(metric) + old_mesh => extract_mesh(states(1), topology_mesh_name) + call allocate(metric, old_mesh, "ErrorMetric") + call assemble_metric(states, metric) + end if + + ewrite(0, *) "Expected nodes = ", expected_nodes(position, metric) + + call allocate(t_edge_lengths, metric%mesh, "TensorEdgeLengths") + call get_edge_lengths(metric, t_edge_lengths) + call vtk_write_fields(trim(output_meshname) // "EdgeLengths", position = position, model = position%mesh, & + & tfields = (/t_edge_lengths, metric/)) + call deallocate(t_edge_lengths) + + ! Generate and the mesh + + new_mesh_field = read_gmsh_file(trim(input_geometryname),& + format_string="geo", & + quad_ngi=ele_ngi(old_mesh,1),& + position = position, metric = metric) + + ! Write the output mesh + call write_mesh_files(output_meshname, "gmsh", new_mesh_field) + + ! Deallocate + do i = 1, size(states) + call deallocate(states(i)) + end do + if (have_metric_name) then + deallocate(states) + else + call deallocate(metric) + deallocate(metric) + end if + call deallocate(new_mesh_field) + + call print_references(0) + + ewrite(1, *) "Exiting flboundadapt" + +end subroutine flboundadapt diff --git a/tools/Flboundadapt_main.cpp b/tools/Flboundadapt_main.cpp new file mode 100644 index 0000000000..c0f915e0f7 --- /dev/null +++ b/tools/Flboundadapt_main.cpp @@ -0,0 +1,155 @@ +/* Copyright (C) 2006 Imperial College London and others. + + Please see the AUTHORS file in the main source directory for a full list + of copyright holders. + + Prof. C Pain + Applied Modelling and Computation Group + Department of Earth Science and Engineering + Imperial College London + + amcgsoftware@imperial.ac.uk + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation, + version 2.1 of the License. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 + USA +*/ + +#include +#include +#include + +#include "confdefs.h" + +#ifdef HAVE_MPI +#include +#endif + +#include "Usage.h" +#include "c++debug.h" + +using namespace std; + +extern "C"{ + void flboundadapt(const char *, const char *, const char *, const char *); +#ifdef HAVE_PYTHON +#include "python_statec.h" +#endif +} + +void Usage(){ + cerr << "Usage: flboundadapt INPUT_FLML INPUT_GEO_FILE OUTPUT_MESH\n" + << "or\n" + << "flboundadapt -i INPUT_VTU METRIC_NAME INPUT_GEO_FILE OUTPUT_MESH\n" + << "\n" + << "Performs mesh generation based on the input options file (which may be a\n" + << "checkpoint). Outputs the resulting mesh.\n" + << "\n" + << "Options:\n" + << "\n" + << "-h\t\tDisplay this help\n" + << "-i\t\t read from VTU file\n" + << "-v\t\tVerbose mode" << endl; +} + +int main(int argc, char** argv){ +#ifdef HAVE_MPI + MPI::Init(argc, argv); + // Undo some MPI init shenanigans + chdir(getenv("PWD")); +#endif + PetscInit(argc, argv); +#ifdef HAVE_PYTHON + // Initialize the Python Interpreter + python_init_(); +#endif + + // Modified version of flredecomp argument parsing + // Get any command line arguments + // Reset optarg so we can detect changes + optarg = NULL; + char c; + map args; + bool vtk_mode = false; + while((c = getopt(argc, argv, "hiv")) != -1){ + if (c != '?'){ + if (c == 'i') vtk_mode=true; + if(optarg == NULL){ + args[c] = "true"; + }else{ + args[c] = optarg; + } + }else{ + if(isprint(optopt)){ + cerr << "Unknown option " << optopt << endl; + }else{ + cerr << "Unknown option " << hex << optopt << endl; + } + Usage(); + exit(-1); + } + } + + // Help + if(args.count('h')){ + Usage(); + exit(0); + } + + // Verbosity + int verbosity = 0; + if(args.count('v') > 0){ + verbosity = 3; + } + set_global_debug_level_fc(&verbosity); + + // Input and output base names + int I = vtk_mode?1:0; + string input_name, input_geometryname, output_meshname, metric_name=""; + if(argc > optind + 3 +I){ + input_name = argv[optind + 1]; + if (vtk_mode) metric_name = argv[optind + 2]; + input_geometryname = argv[optind + 2 + I]; + output_meshname = argv[optind + 3 + I]; + }else if(argc == optind + 3+(vtk_mode?1:0) ){ + input_name = argv[optind]; + if (vtk_mode) metric_name = argv[optind + 1]; + input_geometryname = argv[optind + 1 +I]; + output_meshname = argv[optind + 2 + I]; + }else{ + Usage(); + exit(-1); + } + + + + flboundadapt(input_name.c_str(), + input_geometryname.c_str(), + output_meshname.c_str(), + metric_name.c_str()); + +#ifdef HAVE_PYTHON + // Finalize the Python Interpreter + python_end_(); +#endif + +#ifdef HAVE_PETSC + PetscFinalize(); +#endif + +#ifdef HAVE_MPI + MPI::Finalize(); +#endif + return 0; +} diff --git a/tools/Makefile.in b/tools/Makefile.in index 2e4491cea5..ff57f8a23c 100644 --- a/tools/Makefile.in +++ b/tools/Makefile.in @@ -75,6 +75,7 @@ SUPERMESH_DIFFERENCE=../bin/supermesh_difference DIFFERENTIATE_VTU=../bin/differentiate_vtu VTU_BINS=../bin/vtu_bins FLADAPT=../bin/fladapt +FLBOUNDADAPT=../bin/flboundadapt VTKPROJECTION=../bin/vtk_projection PERIODISE=../bin/periodise @@ -94,6 +95,9 @@ ifeq (@GFORTRAN_4_5_OR_NEWER@,yes) BINARIES += $(VISUALISE_ELEMENTS) endif +ifeq (@HAVE_LIBGMSH@,yes) +BINARIES += $(FLBOUNDADAPT) +endif .SUFFIXES: .f90 .F90 .c .cpp .o .a @@ -208,6 +212,9 @@ $(VTU_BINS): Vtu_Bins.o Vtu_Bins_main.o lib/ $(FLADAPT): Fladapt.o Fladapt_main.o lib/ $(LINKER) -o $@ $(filter %.o,$^) -l$(FLUIDITY) $(LIBS) +$(FLBOUNDADAPT): Flboundadapt.o Flboundadapt_main.o lib/ + $(LINKER) -o $@ $(filter %.o,$^) -l$(FLUIDITY) $(LIBS) + $(UNIFIEDMESH): unifiedmesh.o unifiedmesh_main.o lib/ $(LINKER) -o $@ $(filter %.o,$^) -l$(FLUIDITY) $(LIBS) diff --git a/tools/Vertical_Integration.F90 b/tools/Vertical_Integration.F90 index 5d4b175f62..895543b470 100644 --- a/tools/Vertical_Integration.F90 +++ b/tools/Vertical_Integration.F90 @@ -78,7 +78,7 @@ subroutine vertical_integration(target_basename_, target_basename_len, & ! Step 1: Read in the data write(*,*) target_basename, output_basename, integrated_filename - positions_b_surf = read_gmsh_file(target_basename, quad_degree = quad_degree) + positions_b_surf = read_gmsh_file(target_basename, quad_degree = quad_degree,format_string="msh") dim = positions_b_surf%dim + 1 call vtk_read_state(integrated_filename, state_a, quad_degree = quad_degree) diff --git a/tools/gmsh2vtu.F90 b/tools/gmsh2vtu.F90 index 163bb17a5a..e0100f2818 100644 --- a/tools/gmsh2vtu.F90 +++ b/tools/gmsh2vtu.F90 @@ -47,7 +47,7 @@ subroutine gmsh2vtu(filename_, filename_len) bind(c) filename(i:i)=filename_(i) end do - positions=read_gmsh_file(filename, quad_degree=3) + positions=read_gmsh_file(filename, quad_degree=3, format_string="msh") if (associated(positions%mesh%region_ids)) then regions=piecewise_constant_field(positions%mesh, name="Regions")