Add Libgmsh support/ flboundadapt

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
 

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 <gmsh/Gmsh.h>
+
+#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 <gmsh/Gmsh.h>
+#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 :

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 <gmsh/Gmsh.h>
+  ],[
+  GmshInitialize ();
+  ],[
+  AC_MSG_RESULT([yes])
+  LIBS="$LIBS $LIBS_bck"
+  AC_DEFINE(HAVE_LIBGMSH, [1])
+  HAVE_LIBGMSH="yes"	
+  ],[
+  LIBS="-lgmsh -lgl2ps"
+  AC_TRY_LINK([
+	#include <gmsh/Gmsh.h>],[
+  	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],

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

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

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

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