From a93459225814e54ceb903b7f23c4f97e781b598a Mon Sep 17 00:00:00 2001 From: James Percival Date: Mon, 2 Oct 2017 19:47:38 +0100 Subject: [PATCH 1/9] Read/write .msh files with libgmsh. --- configure | 78 +++++++++++++++--- configure.in | 17 ++++ femtools/Makefile.in | 2 +- femtools/Read_GMSH.F90 | 170 +++++++++++++++++++++++++++++++++++++--- femtools/Write_GMSH.F90 | 141 ++++++++++++++++++++++++++++++--- 5 files changed, 372 insertions(+), 36 deletions(-) diff --git a/configure b/configure index 7b4a3d2283..9f21c797f1 100755 --- a/configure +++ b/configure @@ -734,7 +734,6 @@ infodir docdir oldincludedir includedir -runstatedir localstatedir sharedstatedir sysconfdir @@ -850,7 +849,6 @@ datadir='${datarootdir}' sysconfdir='${prefix}/etc' sharedstatedir='${prefix}/com' localstatedir='${prefix}/var' -runstatedir='${localstatedir}/run' includedir='${prefix}/include' oldincludedir='/usr/include' docdir='${datarootdir}/doc/${PACKAGE_TARNAME}' @@ -1103,15 +1101,6 @@ do | -silent | --silent | --silen | --sile | --sil) silent=yes ;; - -runstatedir | --runstatedir | --runstatedi | --runstated \ - | --runstate | --runstat | --runsta | --runst | --runs \ - | --run | --ru | --r) - ac_prev=runstatedir ;; - -runstatedir=* | --runstatedir=* | --runstatedi=* | --runstated=* \ - | --runstate=* | --runstat=* | --runsta=* | --runst=* | --runs=* \ - | --run=* | --ru=* | --r=*) - runstatedir=$ac_optarg ;; - -sbindir | --sbindir | --sbindi | --sbind | --sbin | --sbi | --sb) ac_prev=sbindir ;; -sbindir=* | --sbindir=* | --sbindi=* | --sbind=* | --sbin=* \ @@ -1249,7 +1238,7 @@ fi for ac_var in exec_prefix prefix bindir sbindir libexecdir datarootdir \ datadir sysconfdir sharedstatedir localstatedir includedir \ oldincludedir docdir infodir htmldir dvidir pdfdir psdir \ - libdir localedir mandir runstatedir + libdir localedir mandir do eval ac_val=\$$ac_var # Remove trailing slashes. @@ -1402,7 +1391,6 @@ Fine tuning of the installation directories: --sysconfdir=DIR read-only single-machine data [PREFIX/etc] --sharedstatedir=DIR modifiable architecture-independent data [PREFIX/com] --localstatedir=DIR modifiable single-machine data [PREFIX/var] - --runstatedir=DIR modifiable per-process data [LOCALSTATEDIR/run] --libdir=DIR object code libraries [EPREFIX/lib] --includedir=DIR C header files [PREFIX/include] --oldincludedir=DIR C header files for non-gcc [/usr/include] @@ -9928,6 +9916,70 @@ 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 + + +else + + { $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 +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..f7765656c7 100644 --- a/configure.in +++ b/configure.in @@ -305,6 +305,23 @@ 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]) + ],[ + AC_MSG_RESULT([no]) + LIBS="$LIBS_bck" + ]) +AC_LANG_POP([C++]) AC_ARG_ENABLE(openmp, [AC_HELP_STRING([--enable-openmp], 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/Read_GMSH.F90 b/femtools/Read_GMSH.F90 index 9339802524..80732403bd 100644 --- a/femtools/Read_GMSH.F90 +++ b/femtools/Read_GMSH.F90 @@ -31,6 +31,7 @@ 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 @@ -87,15 +88,78 @@ 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 integer :: dim, coordinate_dim, gdim integer :: gmshFormat integer :: n, d, e, f, nodeID +#ifdef HAVE_LIBGMSH + type(c_ptr) :: gmodel +#endif + type(GMSHnode), pointer :: nodes(:) type(GMSHelement), pointer :: elements(:), faces(:) - + 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 cread_gmsh_sizes(gmodel, numNodes, numFaces, numElements,& + haveRegionIDs, haveBounds, haveElementOwners, 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 + 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, & + boundaryIDs, faceOwner) bind(c) + use iso_c_binding + type(c_ptr) :: gmodel + integer(c_int) :: numFaces, sloc + integer(c_int) :: sndglno(*), boundaryIDs(*),& + faceOwner(numFaces) + end subroutine cread_gmsh_face_connectivity + end interface + ! If running in parallel, add the process number if(isparallel()) then lfilename = trim(parallel_filename(filename)) // ".msh" @@ -103,6 +167,86 @@ function read_gmsh_simple( filename, quad_degree, & lfilename = trim(filename) // ".msh" end if +#if HAVE_LIBGMSH + + call cread_gmsh_file(gmodel, trim(lfilename)//c_null_char) + call cread_gmsh_sizes(gmodel, numNodes, numFaces, numElements,& + haveRegionIDs, haveBounds, haveElementOwners, 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") + + ! 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 + + 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, boundaryIDs, 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 + ewrite(2,*) "WARNING: no boundaries in GMSH file "//trim(lfilename) + call add_faces( field%mesh, sndgln = sndglno(1:numFaces*sloc) ) + end if + + call cgmsh_finalise(gmodel) + + deallocate(sndglno) + if (haveBounds) deallocate(boundaryIDs) + if (haveElementOwners) deallocate(faceOwner) + +#else fd = free_unit() ! Open node file @@ -173,7 +317,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 +338,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 +424,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..4fbec831ac 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 @@ -82,9 +83,117 @@ 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 + + integer :: numNodes, numElements, numFaces, sloc, loc, dim + integer, allocatable, dimension(:) :: sndglno + character(len=longStringLen) :: meshFile + + type(c_ptr) :: gmodel + + 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 cmesh_to_gmodel(gmodel, binary, numNodes,& + numElements, numFaces, loc, sloc, gdim, val,& + etype, eles, ftype, faces, ele_ids, face_ids) bind(c) + use iso_c_binding + type(c_ptr) :: gmodel + integer(c_int) :: binary, numNodes, numElements, numFaces,& + loc, sloc, gdim, 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 + end subroutine cmesh_to_gmodel + end interface + + interface + subroutine cwrite_gmsh_file(gmodel, filename) bind(c) + use iso_c_binding + type(c_ptr) :: gmodel + character(c_char) :: filename(*) + end subroutine cwrite_gmsh_file + end interface + + numNodes = node_count(positions) + numElements = element_count(positions) + numFaces = unique_surface_element_count(positions%mesh) + + dim = mesh_dim(positions) + loc = ele_loc(positions, 1) + sloc = 1 + if (numFaces>0) sloc = face_loc(positions, 1) + + allocate(sndglno(numFaces*sloc)) + call getsndgln(positions%mesh, sndglno) + + call cmesh_to_gmodel(gmodel, merge(1,0,useBinaryGmsh), numNodes,& + numElements, numFaces, loc, sloc,& + dim, positions%val, & + gmsh_type(loc, dim), positions%mesh%ndglno,& + gmsh_type(sloc, dim-1), sndglno, c_loc(positions%mesh%region_ids),& + c_loc(positions%mesh%faces%boundary_ids)) + + deallocate(sndglno) + + meshFile = trim(filename) // ".msh" + + call cwrite_gmsh_file(gmodel, trim(meshFile)//c_null_char) + + call cgmsh_finalise(gmodel) + + 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 write_libgmsh + +#endif + ! ----------------------------------------------------------------- @@ -105,6 +214,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 +254,9 @@ subroutine write_positions_to_gmsh(filename, positions, number_of_partitions) ! Close GMSH file close( fileDesc ) - end if + end if + +#endif return @@ -364,17 +485,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 +522,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 From 490349c9dbf4ffdc59ba382e94d2ce2b859c49ff Mon Sep 17 00:00:00 2001 From: James Percival Date: Wed, 4 Oct 2017 13:57:33 +0100 Subject: [PATCH 2/9] Missing file. --- femtools/LibGmsh_integration.cpp | 335 +++++++++++++++++++++++++++++++ 1 file changed, 335 insertions(+) create mode 100644 femtools/LibGmsh_integration.cpp diff --git a/femtools/LibGmsh_integration.cpp b/femtools/LibGmsh_integration.cpp new file mode 100644 index 0000000000..408293b4c1 --- /dev/null +++ b/femtools/LibGmsh_integration.cpp @@ -0,0 +1,335 @@ +/* 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/Gmsh.h" +#include "gmsh/GModel.h" +#include "gmsh/GEntity.h" +#include "gmsh/MVertex.h" +#include "gmsh/Context.h" +#include "gmsh/PView.h" +#include "gmsh/discreteVertex.h" +#include "gmsh/discreteEdge.h" +#include "gmsh/discreteFace.h" +#include "gmsh/discreteRegion.h" + +#include +#include +#include + +typedef std::vector EntityVector; +typedef std::vector VertexVector; + +int dim; + +int FluidityNodeOrdering(int type, int node) { + 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; +} + +extern "C" { + + void cgmsh_initialise() { + GmshInitialize(); + } + + void cgmsh_finalise(GModel *&gm) { + delete gm; + GmshFinalize(); + } + + void cread_gmsh_file(GModel *&gm, const char* filename) { + GmshInitialize(); + gm = new GModel; + CTX::instance()->mesh.switchElementTags=1; + gm->readMSH(filename); + PView::readMSH(filename); + } + + void cread_gmsh_sizes(GModel *&gm, int &numNodes, int &numFaces, + int &numElements, bool &haveRegionIDs, + bool &haveBounds, bool &haveElementOwners, + int &gdim, int &loc, int &sloc) { + + numNodes = gm->getNumMeshVertices(); + + dim = 0; + if (gm->getNumRegions()>0) { + dim = 3; + } else if (gm->getNumFaces()>0) { + dim = 2; + } else if (gm->getNumEdges()>0) { + dim = 1; + } + gdim = dim; + + numElements = 0; + numFaces = 0; + haveRegionIDs = false; + haveBounds = false; + haveElementOwners = false; + std::vector eles; + std::vector faces; + gm->getEntities(eles, dim); + gm->getEntities(faces, dim-1); + 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(); + haveRegionIDs = haveRegionIDs || (*it)->tag()>0; + } + for (EntityVector::iterator it = faces.begin(); it != faces.end(); ++it) { + numFaces = numFaces + (*it)->getNumMeshElements(); + haveBounds = haveBounds || (*it)->tag()>0; + haveElementOwners = haveElementOwners || (*it)->getPhysicalEntities().size()==3; + } + } + + void cread_gmsh_element_connectivity(void *&gmv, int &numElements, + int &loc, int *ndglno, int *regionIDs) { + + GModel *gm = (GModel*) gmv; + int k=0; + std::vector ents; + gm->getEntities(ents, dim); + for (EntityVector::iterator ent = ents.begin(); ent != ents.end(); ++ent) { + for (unsigned int i=0; i< (*ent)->getNumMeshElements(); ++i) { + MElement* e = (*ent)->getMeshElement(i); + for (int j=0; jgetNumVertices(); ++j) { + MVertex* v = e->getVertex(FluidityNodeOrdering(e->getType(),j)); + ndglno[loc*k+j] = v->getIndex(); + } + if ((*ent)->tag()>0) regionIDs[k] = (*ent)->tag(); + ++k; + } + } + } + + void cread_gmsh_points(GModel *&gm, int &dim, int &numNodes, double *val) { + for (int i=0; igetMeshVertexByTag(i+1); + switch(dim) { + case 3: + val[dim*i+2] = v->z(); + case 2: + val[dim*i+1] = v->y(); + case 1: + val[dim*i+0] = v->x(); + } + } + } + + void cread_gmsh_face_connectivity(GModel *&gm, int *numFaces, + int &sloc, int *sndglno, + int *boundaryIDs, int *faceOwner) { + int k=0; + std::vector ents; + gm->getEntities(ents, dim-1); + for (EntityVector::iterator ent = ents.begin(); ent != ents.end(); ++ent) { + for (unsigned int i=0; i< (*ent)->getNumMeshElements(); ++i) { + MElement* e = (*ent)->getMeshElement(i); + for (int j=0; jgetNumVertices(); ++j) { + MVertex* v = e->getVertex(FluidityNodeOrdering(e->getType(),j)); + sndglno[sloc*k+j] = v->getIndex(); + } + if ((*ent)->tag()>0) boundaryIDs[k] = (*ent)->tag(); + if ((*ent)->getPhysicalEntities().size()==3) + faceOwner[k] = (*ent)->getPhysicalEntities()[2]; + ++k; + } + } + } + + void cmesh_to_gmodel(GModel *&gm, int &binary, + int &numNodes, int &numElements, int &numFaces, + int &loc, int &sloc, + int &gdim, double* val, + int &etype, int* eles, + int &ftype, int* faces, + int *ele_ids, int *face_ids) { + + GmshInitialize(); + gm = new GModel; + + GEntity *e; + std::map point; + std::map edge; + std::map face; + std::map region; + + std::set uele_ids; + std::set uface_ids; + + uele_ids.insert(0); + if (ele_ids) { + for (int i=0;i::iterator it=uele_ids.begin(); + it != uele_ids.end(); ++it) { + edge[*it] = new discreteEdge(gm,0,NULL,NULL); + edge[*it]->addPhysicalEntity(*it); + edge[*it]->setTag(*it); + gm->add(edge[*it]); + } + for (std::set::iterator it=uface_ids.begin(); + it != uface_ids.end(); ++it) { + point[*it] = new discreteVertex(gm,0); + point[*it]->addPhysicalEntity(*it); + point[*it]->setTag(*it); + gm->add(point[*it]); + } + e = edge[0]; + break; + case 2: + for (std::set::iterator it=uele_ids.begin(); + it != uele_ids.end(); ++it) { + face[*it] = new discreteFace(gm,0); + face[*it]->addPhysicalEntity(*it); + face[*it]->setTag(*it); + gm->add(face[*it]); + } + for (std::set::iterator it=uface_ids.begin(); + it != uface_ids.end(); ++it) { + edge[*it] = new discreteEdge(gm,0,NULL,NULL); + edge[*it]->addPhysicalEntity(*it); + edge[*it]->setTag(*it); + gm->add(edge[*it]); + } + e=face[0]; + break; + case 3: + for (std::set::iterator it=uele_ids.begin(); + it != uele_ids.end(); ++it) { + region[*it] = new discreteRegion(gm,0); + region[*it]->addPhysicalEntity(*it); + region[*it]->setTag(*it); + gm->add(region[*it]); + } + for (std::set::iterator it=uface_ids.begin(); + it != uface_ids.end(); ++it) { + face[*it] = new discreteFace(gm,0); + face[*it]->addPhysicalEntity(*it); + face[*it]->setTag(*it); + gm->add(face[*it]); + } + e=region[0]; + break; + } + + MElementFactory f; + + double x[3] = {0,0,0}; + for (int n=0;nsetIndex(n+1); + e->addMeshVertex(v); + } + + for (int n=0;ngetMeshVertex(eles[loc*n+i]-1)); + } + MElement *ele = f.create(etype, vertices); + switch(gdim) { + case 1: + edge[id]->addLine((MLine*) ele); + break; + case 2: + face[id]->addTriangle((MTriangle*) ele); + break; + case 3: + region[id]->addTetrahedron((MTetrahedron*) ele); + break; + } + } + + for (int n=0;ngetMeshVertex(faces[sloc*n+i]-1)); + } + MElement *ele = f.create(ftype, vertices); + switch(gdim) { + case 1: + point[id]->addPoint((MPoint*) ele); + break; + case 2: + edge[id]->addLine((MLine*) ele); + break; + case 3: + face[id]->addTriangle((MTriangle*) ele); + break; + } + } + + } + + void cwrite_gmsh_file(GModel *&gm, char* filename) { + gm->save(filename); + } + +} + +#endif From 07240189e8639bebd0752787c87b062828d54192 Mon Sep 17 00:00:00 2001 From: James Percival Date: Mon, 9 Oct 2017 15:56:52 +0100 Subject: [PATCH 3/9] Fully support column ids and faceOwners with libgmsh. --- Makefile.in | 3 + configure | 58 +++++++- configure.in | 24 +++- femtools/LibGmsh_integration.cpp | 218 +++++++++++++++++++++++++------ femtools/Read_GMSH.F90 | 45 +++++-- femtools/Write_GMSH.F90 | 93 ++++++++++--- 6 files changed, 364 insertions(+), 77 deletions(-) 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 9f21c797f1..314add3329 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 @@ -9929,7 +9930,7 @@ LIBS="-lGmsh" cat confdefs.h - <<_ACEOF >conftest.$ac_ext /* end confdefs.h. */ - #include + #include #ifdef F77_DUMMY_MAIN @@ -9964,16 +9965,65 @@ $as_echo "yes" >&6; } LIBS="$LIBS $LIBS_bck" $as_echo "#define HAVE_LIBGMSH 1" >>confdefs.h + HAVE_LIBGMSH = yes else - { $as_echo "$as_me:${as_lineno-$LINENO}: result: no" >&5 + 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 + + HAVE_LIBGMSH = no + { $as_echo "$as_me:${as_lineno-$LINENO}: result: no" >&5 $as_echo "no" >&6; } - LIBS="$LIBS_bck" + 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' @@ -13929,7 +13979,7 @@ fi ##### # -LIBS="$LIBS -L./lib" +LIB="$LIBS -L./lib" ############################################################## # Stream I/O diff --git a/configure.in b/configure.in index f7765656c7..ef8fe58ac9 100644 --- a/configure.in +++ b/configure.in @@ -310,17 +310,31 @@ AC_LANG_PUSH([C++]) LIBS_bck="$LIBS" LIBS="-lGmsh" AC_TRY_LINK([ - #include + #include ],[ GmshInitialize (); ],[ AC_MSG_RESULT([yes]) LIBS="$LIBS $LIBS_bck" - AC_DEFINE(HAVE_LIBGMSH, [1]) + AC_DEFINE(HAVE_LIBGMSH, [1]) + HAVE_LIBGMSH = yes ],[ - AC_MSG_RESULT([no]) - LIBS="$LIBS_bck" + LIBS="-lgmsh" + 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, @@ -1084,7 +1098,7 @@ fi ##### # -LIBS="$LIBS -L./lib" +LIB="$LIBS -L./lib" ############################################################## # Stream I/O diff --git a/femtools/LibGmsh_integration.cpp b/femtools/LibGmsh_integration.cpp index 408293b4c1..e4a4ca4493 100644 --- a/femtools/LibGmsh_integration.cpp +++ b/femtools/LibGmsh_integration.cpp @@ -30,10 +30,16 @@ #ifdef HAVE_LIBGMSH #include "gmsh/Gmsh.h" #include "gmsh/GModel.h" +#include "gmsh/JacobianBasis.h" #include "gmsh/GEntity.h" +#include "gmsh/MElement.h" +#include "gmsh/MTetrahedron.h" +#include "gmsh/MLine.h" +#include "gmsh/MTriangle.h" #include "gmsh/MVertex.h" #include "gmsh/Context.h" #include "gmsh/PView.h" +#include "gmsh/PViewDataGModel.h" #include "gmsh/discreteVertex.h" #include "gmsh/discreteEdge.h" #include "gmsh/discreteFace.h" @@ -48,6 +54,9 @@ typedef std::vector VertexVector; int dim; +// Length set to match FIELD_NAME_LEN +char buffer[102]=""; + int FluidityNodeOrdering(int type, int node) { int quads[]= {0,1,3,2}; int hexes[]= {0,1,3,2,4,5,7,6}; @@ -60,20 +69,37 @@ int FluidityNodeOrdering(int type, int node) { return node; } + +int GmshNodeOrdering(int type, int node) { + 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 GmshInitialize(); gm = new GModel; + // This swith swaps how Gmsh uses the first two integer tags on elements CTX::instance()->mesh.switchElementTags=1; gm->readMSH(filename); PView::readMSH(filename); @@ -82,10 +108,13 @@ extern "C" { 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) { + // Number of mesh nodes in model numNodes = gm->getNumMeshVertices(); + // Get the largest dimensionality represented dim = 0; if (gm->getNumRegions()>0) { dim = 3; @@ -101,10 +130,15 @@ extern "C" { haveRegionIDs = false; haveBounds = false; haveElementOwners = false; + haveColumns = not (PView::getViewByName("column_ids") == NULL); + std::vector ents; std::vector eles; std::vector faces; - gm->getEntities(eles, dim); - gm->getEntities(faces, dim-1); + gm->getEntities(ents); + for (EntityVector::iterator it = ents.begin(); it != ents.end(); ++it) { + if ((*it)->dim() == dim) eles.push_back(*it); + if ((*it)->dim() == dim-1) faces.push_back(*it); + } if (eles.size()>0) { loc = eles[0]->getMeshElement(0)->getNumVertices(); } else { @@ -119,10 +153,15 @@ extern "C" { numElements = numElements + (*it)->getNumMeshElements(); haveRegionIDs = haveRegionIDs || (*it)->tag()>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; - haveElementOwners = haveElementOwners || (*it)->getPhysicalEntities().size()==3; } } @@ -132,8 +171,12 @@ extern "C" { GModel *gm = (GModel*) gmv; int k=0; std::vector ents; - gm->getEntities(ents, dim); - for (EntityVector::iterator ent = ents.begin(); ent != ents.end(); ++ent) { + std::vector eles; + gm->getEntities(ents); + for (EntityVector::iterator it = ents.begin(); it != ents.end(); ++it) { + if ((*it)->dim() == dim) eles.push_back(*it); + } + for (EntityVector::iterator ent = eles.begin(); ent != eles.end(); ++ent) { for (unsigned int i=0; i< (*ent)->getNumMeshElements(); ++i) { MElement* e = (*ent)->getMeshElement(i); for (int j=0; jgetNumVertices(); ++j) { @@ -160,34 +203,48 @@ extern "C" { } } - void cread_gmsh_face_connectivity(GModel *&gm, int *numFaces, + void cread_gmsh_face_connectivity(GModel *&gm, int &numFaces, int &sloc, int *sndglno, - int *boundaryIDs, int *faceOwner) { + bool &haveBounds, int *boundaryIDs, + bool &haveElementOwners, int *faceOwner) { int k=0; std::vector ents; - gm->getEntities(ents, dim-1); - for (EntityVector::iterator ent = ents.begin(); ent != ents.end(); ++ent) { + std::vector faces; + gm->getEntities(ents); + for (EntityVector::iterator it = ents.begin(); it != ents.end(); ++it) { + if ((*it)->dim() == dim-1) faces.push_back(*it); + } + for (EntityVector::iterator ent = faces.begin(); ent != faces.end(); ++ent) { for (unsigned int i=0; i< (*ent)->getNumMeshElements(); ++i) { MElement* e = (*ent)->getMeshElement(i); for (int j=0; jgetNumVertices(); ++j) { MVertex* v = e->getVertex(FluidityNodeOrdering(e->getType(),j)); sndglno[sloc*k+j] = v->getIndex(); } - if ((*ent)->tag()>0) boundaryIDs[k] = (*ent)->tag(); - if ((*ent)->getPhysicalEntities().size()==3) - faceOwner[k] = (*ent)->getPhysicalEntities()[2]; + if (haveBounds) { + if ((*ent)->tag()>0) { + boundaryIDs[k] = (*ent)->tag(); + } else { + boundaryIDs[k] = 0; + } + } + if (haveElementOwners) { + if(e->getPartition()>0) { + faceOwner[k] = e->getPartition(); + } + } ++k; } } } - void cmesh_to_gmodel(GModel *&gm, int &binary, + void cmesh_to_gmodel(GModel *&gm, int &numNodes, int &numElements, int &numFaces, int &loc, int &sloc, - int &gdim, double* val, + int &gdim, int &pdim, double* val, int &etype, int* eles, int &ftype, int* faces, - int *ele_ids, int *face_ids) { + int *ele_ids, int *face_ids, int *eleOwners) { GmshInitialize(); gm = new GModel; @@ -206,7 +263,7 @@ extern "C" { for (int i=0;i::iterator it=uele_ids.begin(); it != uele_ids.end(); ++it) { - edge[*it] = new discreteEdge(gm,0,NULL,NULL); + edge[*it] = new discreteEdge(gm,*it,NULL,NULL); edge[*it]->addPhysicalEntity(*it); edge[*it]->setTag(*it); gm->add(edge[*it]); } + e = edge[0]; for (std::set::iterator it=uface_ids.begin(); it != uface_ids.end(); ++it) { - point[*it] = new discreteVertex(gm,0); + point[*it] = new discreteVertex(gm,*it); point[*it]->addPhysicalEntity(*it); point[*it]->setTag(*it); gm->add(point[*it]); } - e = edge[0]; break; case 2: for (std::set::iterator it=uele_ids.begin(); it != uele_ids.end(); ++it) { - face[*it] = new discreteFace(gm,0); + face[*it] = new discreteFace(gm,*it); face[*it]->addPhysicalEntity(*it); face[*it]->setTag(*it); gm->add(face[*it]); } + e=face[0]; for (std::set::iterator it=uface_ids.begin(); it != uface_ids.end(); ++it) { - edge[*it] = new discreteEdge(gm,0,NULL,NULL); + edge[*it] = new discreteEdge(gm,*it,NULL,NULL); edge[*it]->addPhysicalEntity(*it); edge[*it]->setTag(*it); gm->add(edge[*it]); } - e=face[0]; break; case 3: for (std::set::iterator it=uele_ids.begin(); it != uele_ids.end(); ++it) { - region[*it] = new discreteRegion(gm,0); + region[*it] = new discreteRegion(gm,*it); region[*it]->addPhysicalEntity(*it); region[*it]->setTag(*it); gm->add(region[*it]); } + e=region[0]; for (std::set::iterator it=uface_ids.begin(); it != uface_ids.end(); ++it) { - face[*it] = new discreteFace(gm,0); + face[*it] = new discreteFace(gm,*it); face[*it]->addPhysicalEntity(*it); face[*it]->setTag(*it); gm->add(face[*it]); } - e=region[0]; break; } @@ -272,8 +329,8 @@ extern "C" { double x[3] = {0,0,0}; for (int n=0;nsetIndex(n+1); @@ -286,19 +343,26 @@ extern "C" { VertexVector vertices; for (int i=0;igetMeshVertex(eles[loc*n+i]-1)); + vertices.push_back(e->getMeshVertex(eles[loc*n+ + GmshNodeOrdering(etype,i)]-1)); } MElement *ele = f.create(etype, vertices); - switch(gdim) { - case 1: + switch(etype) { + case MSH_LIN_2: edge[id]->addLine((MLine*) ele); break; - case 2: + case MSH_TRI_3: face[id]->addTriangle((MTriangle*) ele); break; - case 3: + case MSH_QUA_4: + face[id]->addQuadrangle((MQuadrangle*) ele); + break; + case MSH_TET_4: region[id]->addTetrahedron((MTetrahedron*) ele); break; + case MSH_HEX_8: + region[id]->addHexahedron((MHexahedron*) ele); + break; } } @@ -308,26 +372,100 @@ extern "C" { VertexVector vertices; for (int i=0;igetMeshVertex(faces[sloc*n+i]-1)); + vertices.push_back(e->getMeshVertex(faces[sloc*n + +GmshNodeOrdering(ftype,i)]-1)); } MElement *ele = f.create(ftype, vertices); - switch(gdim) { - case 1: + if (eleOwners) ele->setPartition(eleOwners[n]); + switch(ftype) { + case MSH_PNT: point[id]->addPoint((MPoint*) ele); break; - case 2: + case MSH_LIN_2: edge[id]->addLine((MLine*) ele); break; - case 3: + case MSH_TRI_3: face[id]->addTriangle((MTriangle*) ele); break; + case MSH_QUA_4: + face[id]->addQuadrangle((MQuadrangle*) ele); + break; } - } + } + } + + 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; } - void cwrite_gmsh_file(GModel *&gm, char* filename) { - gm->save(filename); + 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) { + std::cout << ((*it)->first.second) << std::endl; + ++(*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; + } } } diff --git a/femtools/Read_GMSH.F90 b/femtools/Read_GMSH.F90 index 80732403bd..dc6d013634 100644 --- a/femtools/Read_GMSH.F90 +++ b/femtools/Read_GMSH.F90 @@ -88,13 +88,15 @@ function read_gmsh_simple( filename, quad_degree, & character(len = parallel_filename_len(filename)) :: lfilename integer :: loc, sloc integer :: numNodes, numElements, numFaces - logical(c_bool) :: 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 type(c_ptr) :: gmodel + real, dimension(:), allocatable :: lcolumn_ids #endif type(GMSHnode), pointer :: nodes(:) @@ -122,11 +124,13 @@ end subroutine cread_gmsh_file interface subroutine cread_gmsh_sizes(gmodel, numNodes, numFaces, numElements,& - haveRegionIDs, haveBounds, haveElementOwners, dim, loc, sloc) bind(c) + 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 + logical(c_bool) :: haveRegionIDs, haveBounds, & + haveElementOwners, haveColumns end subroutine cread_gmsh_sizes end interface @@ -150,15 +154,27 @@ end subroutine cread_gmsh_points end interface interface - subroutine cread_gmsh_face_connectivity(gmodel, numFaces, sloc, sndglno, & - boundaryIDs, faceOwner) bind(c) + 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 - integer(c_int) :: sndglno(*), boundaryIDs(*),& - faceOwner(numFaces) + logical(c_bool) :: haveBounds, haveElementOwners + integer(c_int) :: sndglno(*), boundaryIDs(*), faceOwner(*) end subroutine cread_gmsh_face_connectivity 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 ! If running in parallel, add the process number if(isparallel()) then @@ -171,7 +187,8 @@ end subroutine cread_gmsh_face_connectivity call cread_gmsh_file(gmodel, trim(lfilename)//c_null_char) call cread_gmsh_sizes(gmodel, numNodes, numFaces, numElements,& - haveRegionIDs, haveBounds, haveElementOwners, dim, loc, sloc) + haveRegionIDs, haveBounds, haveElementOwners, & + haveColumns, dim, loc, sloc) if (present(mdim)) then coordinate_dim = mdim @@ -203,6 +220,9 @@ end subroutine cread_gmsh_face_connectivity 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) @@ -221,7 +241,8 @@ end subroutine cread_gmsh_face_connectivity call cread_gmsh_face_connectivity(gmodel, numFaces, sloc, & - sndglno, boundaryIDs, faceOwner) + sndglno, haveBounds, boundaryIDs, & + haveElementOwners, faceOwner) ! If we've got boundaries, do something if( haveBounds ) then @@ -240,6 +261,12 @@ end subroutine cread_gmsh_face_connectivity 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) diff --git a/femtools/Write_GMSH.F90 b/femtools/Write_GMSH.F90 index 4fbec831ac..c6b724c926 100644 --- a/femtools/Write_GMSH.F90 +++ b/femtools/Write_GMSH.F90 @@ -52,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=.false. contains @@ -92,11 +92,14 @@ subroutine write_libgmsh(filename, positions) character(len=*), intent(in):: filename type(vector_field), intent(in):: positions - integer :: numNodes, numElements, numFaces, sloc, loc, dim + integer :: numNodes, numElements, numFaces, sloc, & + loc, dim, pdim, i integer, allocatable, dimension(:) :: sndglno + integer, allocatable, dimension(:), target ::owners character(len=longStringLen) :: meshFile + logical needs_element_owners - type(c_ptr) :: gmodel + type(c_ptr) :: gmodel, bnd_ids, reg_ids, pvdata, ele_owners interface subroutine cgmsh_initialise() bind(c) @@ -111,32 +114,56 @@ end subroutine cgmsh_finalise end interface interface - subroutine cmesh_to_gmodel(gmodel, binary, numNodes,& - numElements, numFaces, loc, sloc, gdim, val,& - etype, eles, ftype, faces, ele_ids, face_ids) bind(c) + 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) :: binary, numNodes, numElements, numFaces,& - loc, sloc, gdim, etype, ftype + 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 + type(c_ptr), value :: ele_ids, face_ids, ele_owners end subroutine cmesh_to_gmodel end interface interface - subroutine cwrite_gmsh_file(gmodel, filename) bind(c) + 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 + 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) @@ -144,18 +171,46 @@ end subroutine cwrite_gmsh_file allocate(sndglno(numFaces*sloc)) call getsndgln(positions%mesh, sndglno) - call cmesh_to_gmodel(gmodel, merge(1,0,useBinaryGmsh), numNodes,& + if (associated(positions%mesh%region_ids)) then + reg_ids = c_loc(positions%mesh%region_ids(1)) + else + reg_ids = c_null_ptr + end if + 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)) + else + bnd_ids = c_null_ptr + 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, positions%val, & + dim, pdim, positions%val, & gmsh_type(loc, dim), positions%mesh%ndglno,& - gmsh_type(sloc, dim-1), sndglno, c_loc(positions%mesh%region_ids),& - c_loc(positions%mesh%faces%boundary_ids)) + gmsh_type(sloc, dim-1), sndglno, reg_ids, bnd_ids, ele_owners) deallocate(sndglno) meshFile = trim(filename) // ".msh" - call cwrite_gmsh_file(gmodel, trim(meshFile)//c_null_char) + 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, & + numNodes, 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) @@ -272,7 +327,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 @@ -317,7 +372,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) @@ -375,7 +430,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 @@ -566,7 +621,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 From 875dd34ee20e31a4db951ae720f912f6e64f49aa Mon Sep 17 00:00:00 2001 From: James Percival Date: Mon, 9 Oct 2017 16:39:46 +0100 Subject: [PATCH 4/9] Add named boundaries. --- femtools/Field_Options.F90 | 23 +- femtools/Fields_Allocates.F90 | 9 + femtools/Fields_Base.F90 | 38 ++- femtools/Fields_Data_Types.F90 | 7 + femtools/Futils.F90 | 39 ++- femtools/LibGmsh_integration.cpp | 1 - femtools/Read_GMSH.F90 | 39 ++- .../Boundary_Conditions_From_Options.F90 | 14 +- schemas/fluidity_options.rnc | 34 ++- schemas/fluidity_options.rng | 32 ++- schemas/prognostic_field_options.rnc | 136 +++++++--- schemas/prognostic_field_options.rng | 128 +++++++--- tests/named_boundary/named_boundary.flml | 240 ++++++++++++++++++ tests/named_boundary/named_boundary.xml | 51 ++++ tests/named_boundary/src/square.geo | 20 ++ 15 files changed, 716 insertions(+), 95 deletions(-) create mode 100644 tests/named_boundary/named_boundary.flml create mode 100644 tests/named_boundary/named_boundary.xml create mode 100644 tests/named_boundary/src/square.geo diff --git a/femtools/Field_Options.F90 b/femtools/Field_Options.F90 index ee1df71725..d1ddd16d2e 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,25 @@ 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 + + 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..bac378c117 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..ff13c8d436 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,39 @@ 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 + + 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 + + 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..0cb6fe3d07 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,31 @@ 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=*) :: f_string + + character(kind=c_char), pointer :: arr(:) + character(:,kind = c_char), pointer :: f_ptr + +!!! associate the array with the c string. + call c_f_pointer(c_string, arr, [strlen(c_string)]) +!!! make into a pointer to a Fortran string + call wrap_to_string(strlen(c_string), arr, f_ptr) +!!! copy happens here + f_string = f_ptr + + 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/LibGmsh_integration.cpp b/femtools/LibGmsh_integration.cpp index e4a4ca4493..3cd68cb4e2 100644 --- a/femtools/LibGmsh_integration.cpp +++ b/femtools/LibGmsh_integration.cpp @@ -453,7 +453,6 @@ extern "C" { (*it) = gm->firstPhysicalName(); } while ((*it) != gm->lastPhysicalName() && (*it)->first.first !=mdim) { - std::cout << ((*it)->first.second) << std::endl; ++(*it); } if ((*it) != gm->lastPhysicalName()) { diff --git a/femtools/Read_GMSH.F90 b/femtools/Read_GMSH.F90 index dc6d013634..38596da3e1 100644 --- a/femtools/Read_GMSH.F90 +++ b/femtools/Read_GMSH.F90 @@ -33,7 +33,7 @@ module read_gmsh use iso_c_binding use fldebug - use global_parameters, only : OPTION_PATH_LEN + use global_parameters, only : OPTION_PATH_LEN, FIELD_NAME_LEN use futils use quadrature use elements @@ -95,7 +95,10 @@ function read_gmsh_simple( filename, quad_degree, & integer :: n, d, e, f, nodeID #ifdef HAVE_LIBGMSH - type(c_ptr) :: gmodel + integer :: snames + type(c_ptr) :: gmodel, c_str, it + integer, allocatable, dimension(:) :: id_list + character(len = FIELD_NAME_LEN), dimension(:), allocatable :: name_list real, dimension(:), allocatable :: lcolumn_ids #endif @@ -166,6 +169,26 @@ subroutine cread_gmsh_face_connectivity(gmodel, numFaces, & 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 @@ -212,6 +235,18 @@ end subroutine cread_gmsh_node_data 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) 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/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/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/named_boundary/named_boundary.flml b/tests/named_boundary/named_boundary.flml new file mode 100644 index 0000000000..944a142891 --- /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; From ef94e5e4e4b000896751348ef4d84a16aa558b12 Mon Sep 17 00:00:00 2001 From: James Percival Date: Mon, 9 Oct 2017 20:58:19 +0100 Subject: [PATCH 5/9] Basic integration for libgmsh meshing inside fluidty executable. --- femtools/LibGmsh_integration.cpp | 39 +++++++++++++++++++++--- femtools/Mesh_Files.F90 | 6 ++-- femtools/Read_GMSH.F90 | 24 ++++++++++++--- preprocessor/Populate_State.F90 | 5 +-- schemas/mesh_options.rnc | 5 +++ schemas/mesh_options.rng | 7 +++++ tests/named_boundary/named_boundary.flml | 4 +-- 7 files changed, 75 insertions(+), 15 deletions(-) diff --git a/femtools/LibGmsh_integration.cpp b/femtools/LibGmsh_integration.cpp index 3cd68cb4e2..6bdd3354a1 100644 --- a/femtools/LibGmsh_integration.cpp +++ b/femtools/LibGmsh_integration.cpp @@ -101,10 +101,25 @@ extern "C" { gm = new GModel; // This swith swaps how Gmsh uses the first two integer tags on elements CTX::instance()->mesh.switchElementTags=1; - gm->readMSH(filename); + gm->load(filename); PView::readMSH(filename); } + void cmesh_gmsh_file(GModel *&gm) { + // Get the largest dimensionality represented + dim = 0; + if (gm->getNumRegions()>0) { + dim = 3; + } else if (gm->getNumFaces()>0) { + dim = 2; + } else if (gm->getNumEdges()>0) { + dim = 1; + } + gm->mesh(dim); + gm->indexMeshVertices(false, 0, true); + CTX::instance()->mesh.switchElementTags=0; + } + void cread_gmsh_sizes(GModel *&gm, int &numNodes, int &numFaces, int &numElements, bool &haveRegionIDs, bool &haveBounds, bool &haveElementOwners, @@ -149,9 +164,16 @@ extern "C" { } else { sloc = 0; } + + bool physicals = CTX::instance()->mesh.switchElementTags==1; + for (EntityVector::iterator it = eles.begin(); it != eles.end(); ++it) { numElements = numElements + (*it)->getNumMeshElements(); - haveRegionIDs = haveRegionIDs || (*it)->tag()>0; + if (physicals) { + haveRegionIDs = haveRegionIDs || (*it)->tag()>0; + } else { + haveRegionIDs = haveRegionIDs || (*it)->getPhysicalEntities().size()>0; + } } for (EntityVector::iterator it = faces.begin(); it != faces.end(); ++it) { @@ -211,6 +233,9 @@ extern "C" { std::vector ents; std::vector faces; gm->getEntities(ents); + + bool physicals = CTX::instance()->mesh.switchElementTags==1; + for (EntityVector::iterator it = ents.begin(); it != ents.end(); ++it) { if ((*it)->dim() == dim-1) faces.push_back(*it); } @@ -222,10 +247,14 @@ extern "C" { sndglno[sloc*k+j] = v->getIndex(); } if (haveBounds) { - if ((*ent)->tag()>0) { - boundaryIDs[k] = (*ent)->tag(); + if (physicals) { + if ((*ent)->tag()>0) { + boundaryIDs[k] = (*ent)->tag(); + } else { + boundaryIDs[k] = 0; + } } else { - boundaryIDs[k] = 0; + boundaryIDs[k] = (*ent)->getPhysicalEntities()[0]; } } if (haveElementOwners) { 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 38596da3e1..ab84375017 100644 --- a/femtools/Read_GMSH.F90 +++ b/femtools/Read_GMSH.F90 @@ -61,12 +61,12 @@ 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 ) & 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. @@ -124,6 +124,13 @@ subroutine cread_gmsh_file(gmodel, filename) bind(c) character(c_char) :: filename(*) end subroutine cread_gmsh_file end interface + + interface + subroutine cmesh_gmsh_file(gmodel) bind(c) + use iso_c_binding + type(c_ptr) :: gmodel + end subroutine cmesh_gmsh_file + end interface interface subroutine cread_gmsh_sizes(gmodel, numNodes, numFaces, numElements,& @@ -201,14 +208,20 @@ end subroutine cread_gmsh_node_data ! 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) // ".msh" + lfilename = trim(filename) // "." // format_string end if #if HAVE_LIBGMSH call cread_gmsh_file(gmodel, trim(lfilename)//c_null_char) + + if (format_string == "geo") then + call cmesh_gmsh_file(gmodel) + end if + call cread_gmsh_sizes(gmodel, numNodes, numFaces, numElements,& haveRegionIDs, haveBounds, haveElementOwners, & haveColumns, dim, loc, sloc) @@ -311,6 +324,9 @@ end subroutine cread_gmsh_node_data #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", & 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/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/tests/named_boundary/named_boundary.flml b/tests/named_boundary/named_boundary.flml index 944a142891..0749ae2252 100644 --- a/tests/named_boundary/named_boundary.flml +++ b/tests/named_boundary/named_boundary.flml @@ -11,8 +11,8 @@ 2 - - + + From 8cb5b0f44a67d1c5e1cd873c39b5ac004187471f Mon Sep 17 00:00:00 2001 From: James Percival Date: Thu, 12 Oct 2017 15:44:35 +0100 Subject: [PATCH 6/9] Add flboundadapt tool to generate GMSH from a mesh metric set up in a .flml file. --- configure | 6 +- configure.in | 6 +- femtools/GMSH_Common.F90 | 268 ++++++++++++++++++++++++++- femtools/LibGmsh_integration.cpp | 230 ++++++++++++++--------- femtools/Read_GMSH.F90 | 120 ++---------- femtools/Write_GMSH.F90 | 142 +------------- tests/flboundadapt/circle.geo | 40 ++++ tests/flboundadapt/flboundadapt.flml | 133 +++++++++++++ tests/flboundadapt/flboundadapt.xml | 24 +++ tools/Flboundadapt.F90 | 187 +++++++++++++++++++ tools/Flboundadapt_main.cpp | 146 +++++++++++++++ tools/Makefile.in | 4 + 12 files changed, 970 insertions(+), 336 deletions(-) create mode 100644 tests/flboundadapt/circle.geo create mode 100644 tests/flboundadapt/flboundadapt.flml create mode 100644 tests/flboundadapt/flboundadapt.xml create mode 100644 tools/Flboundadapt.F90 create mode 100644 tools/Flboundadapt_main.cpp diff --git a/configure b/configure index 314add3329..77b00fe6b8 100755 --- a/configure +++ b/configure @@ -9965,7 +9965,7 @@ $as_echo "yes" >&6; } LIBS="$LIBS $LIBS_bck" $as_echo "#define HAVE_LIBGMSH 1" >>confdefs.h - HAVE_LIBGMSH = yes + HAVE_LIBGMSH="yes" else @@ -10007,11 +10007,11 @@ $as_echo "yes" >&6; } LIBS="$LIBS $LIBS_bck" $as_echo "#define HAVE_LIBGMSH 1" >>confdefs.h - HAVE_LIBGMSH = yes + HAVE_LIBGMSH="yes" else - HAVE_LIBGMSH = no + HAVE_LIBGMSH="no" { $as_echo "$as_me:${as_lineno-$LINENO}: result: no" >&5 $as_echo "no" >&6; } LIBS="$LIBS_bck" diff --git a/configure.in b/configure.in index ef8fe58ac9..c86660027d 100644 --- a/configure.in +++ b/configure.in @@ -317,7 +317,7 @@ AC_TRY_LINK([ AC_MSG_RESULT([yes]) LIBS="$LIBS $LIBS_bck" AC_DEFINE(HAVE_LIBGMSH, [1]) - HAVE_LIBGMSH = yes + HAVE_LIBGMSH="yes" ],[ LIBS="-lgmsh" AC_TRY_LINK([ @@ -327,9 +327,9 @@ AC_TRY_LINK([ AC_MSG_RESULT([yes]) LIBS="$LIBS $LIBS_bck" AC_DEFINE(HAVE_LIBGMSH, [1]) - HAVE_LIBGMSH = yes + HAVE_LIBGMSH="yes" ],[ - HAVE_LIBGMSH = no + HAVE_LIBGMSH="no" AC_MSG_RESULT([no]) LIBS="$LIBS_bck" ]) diff --git a/femtools/GMSH_Common.F90 b/femtools/GMSH_Common.F90 index 64201ec32f..4a870ce792 100644 --- a/femtools/GMSH_Common.F90 +++ b/femtools/GMSH_Common.F90 @@ -29,9 +29,15 @@ ! 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 + + implicit none + character(len=3), parameter :: GMSHVersionStr = "2.1" integer, parameter :: asciiFormat = 0 integer, parameter :: binaryFormat = 1 @@ -58,6 +64,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 +279,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 +313,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 +357,117 @@ 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)) + 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 + 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)) + else + bnd_ids = c_null_ptr + 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 index 6bdd3354a1..6233231855 100644 --- a/femtools/LibGmsh_integration.cpp +++ b/femtools/LibGmsh_integration.cpp @@ -30,7 +30,7 @@ #ifdef HAVE_LIBGMSH #include "gmsh/Gmsh.h" #include "gmsh/GModel.h" -#include "gmsh/JacobianBasis.h" +#include "gmsh/Field.h" #include "gmsh/GEntity.h" #include "gmsh/MElement.h" #include "gmsh/MTetrahedron.h" @@ -52,12 +52,11 @@ typedef std::vector EntityVector; typedef std::vector VertexVector; -int dim; - // 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) { @@ -71,6 +70,7 @@ int FluidityNodeOrdering(int type, int 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) { @@ -97,26 +97,29 @@ extern "C" { 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 - GmshInitialize(); gm = new GModel; - // This swith swaps how Gmsh uses the first two integer tags on elements + // 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) { - // Get the largest dimensionality represented - dim = 0; - if (gm->getNumRegions()>0) { - dim = 3; - } else if (gm->getNumFaces()>0) { - dim = 2; - } else if (gm->getNumEdges()>0) { - dim = 1; + 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; + } } - gm->mesh(dim); - gm->indexMeshVertices(false, 0, true); + // This causes the meshing to happen + gm->mesh(gm->getDim()); CTX::instance()->mesh.switchElementTags=0; } @@ -125,20 +128,13 @@ extern "C" { bool &haveBounds, bool &haveElementOwners, bool &haveColumns, int &gdim, int &loc, int &sloc) { - - // Number of mesh nodes in model - numNodes = gm->getNumMeshVertices(); + // This utility routine provides Fortran with the sizes to allocate + // memory for. // Get the largest dimensionality represented - dim = 0; - if (gm->getNumRegions()>0) { - dim = 3; - } else if (gm->getNumFaces()>0) { - dim = 2; - } else if (gm->getNumEdges()>0) { - dim = 1; - } - gdim = dim; + gdim = gm->getDim(); + + bool physicals = CTX::instance()->mesh.switchElementTags==1; numElements = 0; numFaces = 0; @@ -146,13 +142,24 @@ extern "C" { haveBounds = false; haveElementOwners = false; haveColumns = not (PView::getViewByName("column_ids") == NULL); - std::vector ents; - std::vector eles; - std::vector faces; + EntityVector ents; + EntityVector eles, all_eles; + EntityVector faces; gm->getEntities(ents); for (EntityVector::iterator it = ents.begin(); it != ents.end(); ++it) { - if ((*it)->dim() == dim) eles.push_back(*it); - if ((*it)->dim() == dim-1) faces.push_back(*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(); @@ -165,8 +172,6 @@ extern "C" { sloc = 0; } - bool physicals = CTX::instance()->mesh.switchElementTags==1; - for (EntityVector::iterator it = eles.begin(); it != eles.end(); ++it) { numElements = numElements + (*it)->getNumMeshElements(); if (physicals) { @@ -176,6 +181,10 @@ extern "C" { } } + if(not haveRegionIDs || physicals) { + eles[0]->addPhysicalEntity(0); + } + for (EntityVector::iterator it = faces.begin(); it != faces.end(); ++it) { if ((*it)->getNumMeshElements()) haveElementOwners = haveElementOwners || @@ -184,43 +193,77 @@ extern "C" { 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); } + + } - void cread_gmsh_element_connectivity(void *&gmv, int &numElements, + 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) { - GModel *gm = (GModel*) gmv; - int k=0; - std::vector ents; - std::vector eles; + 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) eles.push_back(*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(); + } + } } - for (EntityVector::iterator ent = eles.begin(); ent != eles.end(); ++ent) { - for (unsigned int i=0; i< (*ent)->getNumMeshElements(); ++i) { - MElement* e = (*ent)->getMeshElement(i); - for (int j=0; jgetNumVertices(); ++j) { - MVertex* v = e->getVertex(FluidityNodeOrdering(e->getType(),j)); - ndglno[loc*k+j] = v->getIndex(); + + 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; } - if ((*ent)->tag()>0) regionIDs[k] = (*ent)->tag(); - ++k; } + ++k; } } void cread_gmsh_points(GModel *&gm, int &dim, int &numNodes, double *val) { - for (int i=0; igetMeshVertexByTag(i+1); - switch(dim) { - case 3: - val[dim*i+2] = v->z(); - case 2: - val[dim*i+1] = v->y(); - case 1: - val[dim*i+0] = v->x(); + 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(); + } + } } } } @@ -229,42 +272,56 @@ extern "C" { int &sloc, int *sndglno, bool &haveBounds, int *boundaryIDs, bool &haveElementOwners, int *faceOwner) { - int k=0; - std::vector ents; - std::vector faces; + + 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)->dim() == dim-1) faces.push_back(*it); - } - for (EntityVector::iterator ent = faces.begin(); ent != faces.end(); ++ent) { - for (unsigned int i=0; i< (*ent)->getNumMeshElements(); ++i) { - MElement* e = (*ent)->getMeshElement(i); - for (int j=0; jgetNumVertices(); ++j) { - MVertex* v = e->getVertex(FluidityNodeOrdering(e->getType(),j)); - sndglno[sloc*k+j] = v->getIndex(); - } - if (haveBounds) { + 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) { - if ((*ent)->tag()>0) { - boundaryIDs[k] = (*ent)->tag(); - } else { - boundaryIDs[k] = 0; - } + face_label[e] = (*it)->tag(); } else { - boundaryIDs[k] = (*ent)->getPhysicalEntities()[0]; + face_label[e] = (*it)->getPhysicalEntities()[0]; } } - if (haveElementOwners) { - if(e->getPartition()>0) { - faceOwner[k] = e->getPartition(); - } - } - ++k; } } + + 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, @@ -275,7 +332,6 @@ extern "C" { int &ftype, int* faces, int *ele_ids, int *face_ids, int *eleOwners) { - GmshInitialize(); gm = new GModel; GEntity *e; @@ -305,7 +361,7 @@ extern "C" { for (std::set::iterator it=uele_ids.begin(); it != uele_ids.end(); ++it) { edge[*it] = new discreteEdge(gm,*it,NULL,NULL); - edge[*it]->addPhysicalEntity(*it); + edge[*it]->addPhysicalEntity(std::max(1,*it)); edge[*it]->setTag(*it); gm->add(edge[*it]); } @@ -322,7 +378,7 @@ extern "C" { for (std::set::iterator it=uele_ids.begin(); it != uele_ids.end(); ++it) { face[*it] = new discreteFace(gm,*it); - face[*it]->addPhysicalEntity(*it); + face[*it]->addPhysicalEntity(std::max(1,*it)); face[*it]->setTag(*it); gm->add(face[*it]); } @@ -339,7 +395,7 @@ extern "C" { for (std::set::iterator it=uele_ids.begin(); it != uele_ids.end(); ++it) { region[*it] = new discreteRegion(gm,*it); - region[*it]->addPhysicalEntity(*it); + region[*it]->addPhysicalEntity(std::max(1,*it)); region[*it]->setTag(*it); gm->add(region[*it]); } diff --git a/femtools/Read_GMSH.F90 b/femtools/Read_GMSH.F90 index ab84375017..6bd63d5cc6 100644 --- a/femtools/Read_GMSH.F90 +++ b/femtools/Read_GMSH.F90 @@ -61,7 +61,8 @@ module read_gmsh ! The main function for reading GMSH files function read_gmsh_simple( filename, quad_degree, & - quad_ngi, quad_family, mdim, format_string ) & + 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. @@ -75,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 @@ -96,7 +100,7 @@ function read_gmsh_simple( filename, quad_degree, & #ifdef HAVE_LIBGMSH integer :: snames - type(c_ptr) :: gmodel, c_str, it + 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 @@ -105,106 +109,6 @@ function read_gmsh_simple( filename, quad_degree, & type(GMSHnode), pointer :: nodes(:) type(GMSHelement), pointer :: elements(:), faces(:) - 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) bind(c) - use iso_c_binding - type(c_ptr) :: gmodel - 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 ! If running in parallel, add the process number if(isparallel()) then @@ -216,10 +120,18 @@ end subroutine cread_gmsh_node_data #if HAVE_LIBGMSH + call cgmsh_initialise() + call cread_gmsh_file(gmodel, trim(lfilename)//c_null_char) if (format_string == "geo") then - call cmesh_gmsh_file(gmodel) + 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,& @@ -325,7 +237,7 @@ end subroutine cread_gmsh_node_data fd = free_unit() if (format_string == "geo") & - FlAbort(" Fluidity must be built with libgmsh support to read .geo files") + FLAbort("Fluidity must be built with libgmsh support to read .geo files") ! Open node file ewrite(2, *) "Opening "//trim(lfilename)//" for reading." diff --git a/femtools/Write_GMSH.F90 b/femtools/Write_GMSH.F90 index c6b724c926..f030ab6547 100644 --- a/femtools/Write_GMSH.F90 +++ b/femtools/Write_GMSH.F90 @@ -52,7 +52,7 @@ module write_gmsh end interface ! Writes to GMSH binary format - can set to ASCII (handy for debugging) - logical(c_bool), parameter :: useBinaryGMSH=.false. + logical(c_bool), parameter :: useBinaryGMSH=.true. contains @@ -92,113 +92,12 @@ subroutine write_libgmsh(filename, positions) character(len=*), intent(in):: filename type(vector_field), intent(in):: positions - integer :: numNodes, numElements, numFaces, sloc, & - loc, dim, pdim, i - integer, allocatable, dimension(:) :: sndglno - integer, allocatable, dimension(:), target ::owners character(len=longStringLen) :: meshFile - logical needs_element_owners - - type(c_ptr) :: gmodel, bnd_ids, reg_ids, pvdata, ele_owners - - 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 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 - - 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)) - 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 - 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)) - else - bnd_ids = c_null_ptr - 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 + type(c_ptr) :: gmodel, pvdata - 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) + call cgmsh_initialise() - deallocate(sndglno) + call position_to_gmodel(positions, gmodel) meshFile = trim(filename) // ".msh" @@ -206,7 +105,7 @@ end subroutine cwrite_gmsh_data_file if (associated(positions%mesh%columns)) then call cdata_to_pview_node_data(gmodel, pvdata, & - numNodes, real(positions%mesh%columns), & + 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) @@ -214,37 +113,6 @@ end subroutine cwrite_gmsh_data_file call cgmsh_finalise(gmodel) - 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 write_libgmsh #endif 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/tools/Flboundadapt.F90 b/tools/Flboundadapt.F90 new file mode 100644 index 0000000000..e3662b0ba0 --- /dev/null +++ b/tools/Flboundadapt.F90 @@ -0,0 +1,187 @@ +! 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) 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 + 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 + + character(len=FIELD_NAME_LEN):: input_flmlname + character(len=FIELD_NAME_LEN):: input_geometryname + character(len=FIELD_NAME_LEN):: output_meshname + integer :: i, stat + type(mesh_type), pointer :: old_mesh + type(state_type), dimension(:), pointer :: states + type(vector_field) :: new_mesh_field + type(vector_field), pointer :: new_mesh_field_ptr, old_mesh_field + type(tensor_field) :: metric, t_edge_lengths + character(len=FIELD_NAME_LEN) :: mesh_format + character(len=OPTION_PATH_LEN) :: mesh_file_name + real, dimension(:,:,:), allocatable :: metric_data + + ! 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) + + ewrite(1, *) "In flboundadapt" + +#ifdef HAVE_PYTHON + call python_init() +#endif + + ewrite(2, *) "Input flml name: " // trim(input_flmlname) + ewrite(2, *) "Input geometry name: " // trim(input_geometryname) + ewrite(2, *) "Output mesh name: " // trim(output_meshname) + + ! 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) + + ! Find the external mesh field + !call find_mesh_field_to_adapt(states(1), old_mesh_field) + old_mesh_field => extract_vector_field(states(1), "Coordinate") + old_mesh => old_mesh_field%mesh + + ! Assemble the error metric + call allocate(metric, old_mesh, "ErrorMetric") + call assemble_metric(states, metric) + + ewrite(0, *) "Expected nodes = ", expected_nodes(old_mesh_field, 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 = old_mesh_field, model = metric%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 = old_mesh_field, 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 + call deallocate(metric) + call deallocate(new_mesh_field) + + call print_references(0) + + ewrite(1, *) "Exiting fladapt" + +contains + + subroutine find_mesh_field_to_adapt(state, mesh_field) + !!< Find the external mesh field to be used by adaptivity + + type(state_type), intent(in) :: state + type(vector_field), pointer :: mesh_field + + character(len = FIELD_NAME_LEN) :: mesh_field_name + type(mesh_type), pointer :: mesh + + call find_mesh_to_adapt(state, mesh) + if(trim(mesh%name) == "CoordinateMesh") then + mesh_field_name = "Coordinate" + else + mesh_field_name = trim(mesh%name) // "Coordinate" + end if + + if(.not. has_vector_field(state, mesh_field_name)) then + FLAbort("External mesh field " // trim(mesh_field_name) // " not found in the system state") + end if + + mesh_field => extract_vector_field(state, mesh_field_name) + + end subroutine find_mesh_field_to_adapt + +end subroutine flboundadapt diff --git a/tools/Flboundadapt_main.cpp b/tools/Flboundadapt_main.cpp new file mode 100644 index 0000000000..dafad8ca1d --- /dev/null +++ b/tools/Flboundadapt_main.cpp @@ -0,0 +1,146 @@ +/* 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 *); +#ifdef HAVE_PYTHON +#include "python_statec.h" +#endif +} + +void Usage(){ + cerr << "Usage: flboundadapt INPUT_FLML 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" + << "-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; + while((c = getopt(argc, argv, "hv")) != -1){ + if (c != '?'){ + 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 + string input_flmlname, input_geometryname, output_meshname; + if(argc > optind + 3){ + input_flmlname = argv[optind + 1]; + input_geometryname = argv[optind + 2]; + output_meshname = argv[optind + 3]; + }else if(argc == optind + 3){ + input_flmlname = argv[optind]; + input_geometryname = argv[optind+1]; + output_meshname = argv[optind + 2]; + }else{ + Usage(); + exit(-1); + } + + + + flboundadapt(input_flmlname.c_str(), + input_geometryname.c_str(), + output_meshname.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..9cae1fe898 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 @@ -208,6 +209,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) From d1e4739c73ddefa44506558990b6ef4e033c26d7 Mon Sep 17 00:00:00 2001 From: James Percival Date: Mon, 16 Oct 2017 18:54:07 +0100 Subject: [PATCH 7/9] Fixes for some of @stephankramer's and @drhodrid's notes, also including better support for standard libgmsh-dev on 16.04. --- configure | 18 ++++++++++--- configure.in | 4 +-- femtools/Field_Options.F90 | 4 +++ femtools/Fields_Allocates.F90 | 2 +- femtools/Fields_Base.F90 | 10 +++++++ femtools/Futils.F90 | 17 +++++++++--- femtools/GMSH_Common.F90 | 4 ++- femtools/Read_GMSH.F90 | 2 +- manual/configuring_fluidity.tex | 5 +++- manual/mesh_formats.tex | 3 ++- manual/visualisation_and_diagnostics.tex | 16 +++++++++++- tools/Flboundadapt.F90 | 33 +++--------------------- tools/Makefile.in | 3 +++ tools/Vertical_Integration.F90 | 2 +- tools/gmsh2vtu.F90 | 2 +- 15 files changed, 78 insertions(+), 47 deletions(-) diff --git a/configure b/configure index 77b00fe6b8..bcfa8c8d66 100755 --- a/configure +++ b/configure @@ -735,6 +735,7 @@ infodir docdir oldincludedir includedir +runstatedir localstatedir sharedstatedir sysconfdir @@ -850,6 +851,7 @@ datadir='${datarootdir}' sysconfdir='${prefix}/etc' sharedstatedir='${prefix}/com' localstatedir='${prefix}/var' +runstatedir='${localstatedir}/run' includedir='${prefix}/include' oldincludedir='/usr/include' docdir='${datarootdir}/doc/${PACKAGE_TARNAME}' @@ -1102,6 +1104,15 @@ do | -silent | --silent | --silen | --sile | --sil) silent=yes ;; + -runstatedir | --runstatedir | --runstatedi | --runstated \ + | --runstate | --runstat | --runsta | --runst | --runs \ + | --run | --ru | --r) + ac_prev=runstatedir ;; + -runstatedir=* | --runstatedir=* | --runstatedi=* | --runstated=* \ + | --runstate=* | --runstat=* | --runsta=* | --runst=* | --runs=* \ + | --run=* | --ru=* | --r=*) + runstatedir=$ac_optarg ;; + -sbindir | --sbindir | --sbindi | --sbind | --sbin | --sbi | --sb) ac_prev=sbindir ;; -sbindir=* | --sbindir=* | --sbindi=* | --sbind=* | --sbin=* \ @@ -1239,7 +1250,7 @@ fi for ac_var in exec_prefix prefix bindir sbindir libexecdir datarootdir \ datadir sysconfdir sharedstatedir localstatedir includedir \ oldincludedir docdir infodir htmldir dvidir pdfdir psdir \ - libdir localedir mandir + libdir localedir mandir runstatedir do eval ac_val=\$$ac_var # Remove trailing slashes. @@ -1392,6 +1403,7 @@ Fine tuning of the installation directories: --sysconfdir=DIR read-only single-machine data [PREFIX/etc] --sharedstatedir=DIR modifiable architecture-independent data [PREFIX/com] --localstatedir=DIR modifiable single-machine data [PREFIX/var] + --runstatedir=DIR modifiable per-process data [LOCALSTATEDIR/run] --libdir=DIR object code libraries [EPREFIX/lib] --includedir=DIR C header files [PREFIX/include] --oldincludedir=DIR C header files for non-gcc [/usr/include] @@ -9969,7 +9981,7 @@ $as_echo "yes" >&6; } else - LIBS="-lgmsh" + LIBS="-lgmsh -lgl2ps" cat confdefs.h - <<_ACEOF >conftest.$ac_ext /* end confdefs.h. */ @@ -13979,7 +13991,7 @@ fi ##### # -LIB="$LIBS -L./lib" +LIBS="$LIBS -L./lib" ############################################################## # Stream I/O diff --git a/configure.in b/configure.in index c86660027d..d0a40fd68e 100644 --- a/configure.in +++ b/configure.in @@ -319,7 +319,7 @@ AC_TRY_LINK([ AC_DEFINE(HAVE_LIBGMSH, [1]) HAVE_LIBGMSH="yes" ],[ - LIBS="-lgmsh" + LIBS="-lgmsh -lgl2ps" AC_TRY_LINK([ #include ],[ GmshInitialize (); @@ -1098,7 +1098,7 @@ fi ##### # -LIB="$LIBS -L./lib" +LIBS="$LIBS -L./lib" ############################################################## # Stream I/O diff --git a/femtools/Field_Options.F90 b/femtools/Field_Options.F90 index d1ddd16d2e..ea0a9efca7 100644 --- a/femtools/Field_Options.F90 +++ b/femtools/Field_Options.F90 @@ -1345,6 +1345,10 @@ subroutine get_surface_ids(bc_path, mesh, surface_ids) 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))) diff --git a/femtools/Fields_Allocates.F90 b/femtools/Fields_Allocates.F90 index bac378c117..84040f4c8a 100644 --- a/femtools/Fields_Allocates.F90 +++ b/femtools/Fields_Allocates.F90 @@ -1104,7 +1104,7 @@ function make_mesh (model, shape, continuity, name) & end if if (associated(model%surface_names)) then allocate(mesh%surface_names(size(model%surface_names))) - mesh%surface_names(:) = model%surface_names(:) + mesh%surface_names = model%surface_names end if call addref(mesh) diff --git a/femtools/Fields_Base.F90 b/femtools/Fields_Base.F90 index ff13c8d436..f6c444168f 100644 --- a/femtools/Fields_Base.F90 +++ b/femtools/Fields_Base.F90 @@ -4207,6 +4207,11 @@ subroutine set_surface_names(mesh, names, 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) @@ -4221,6 +4226,11 @@ subroutine get_surface_ids_from_names(mesh, names, ids) 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 diff --git a/femtools/Futils.F90 b/femtools/Futils.F90 index 0cb6fe3d07..2447adc673 100644 --- a/femtools/Futils.F90 +++ b/femtools/Futils.F90 @@ -450,17 +450,26 @@ end subroutine tokenize subroutine copy_c_string_to_fortran(c_string, f_string) type(c_ptr), intent(in), target :: c_string - character(len=*) :: f_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, [strlen(c_string)]) + call c_f_pointer(c_string, arr, [c_len]) !!! make into a pointer to a Fortran string - call wrap_to_string(strlen(c_string), arr, f_ptr) + call wrap_to_string(c_len, arr, f_ptr) + if (len(f_string) .ge. c_len) then !!! copy happens here - f_string = f_ptr + 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 diff --git a/femtools/GMSH_Common.F90 b/femtools/GMSH_Common.F90 index 4a870ce792..9918ec0b83 100644 --- a/femtools/GMSH_Common.F90 +++ b/femtools/GMSH_Common.F90 @@ -34,7 +34,9 @@ module gmsh_common use iso_c_binding - use fields + 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 diff --git a/femtools/Read_GMSH.F90 b/femtools/Read_GMSH.F90 index 6bd63d5cc6..a851d41a47 100644 --- a/femtools/Read_GMSH.F90 +++ b/femtools/Read_GMSH.F90 @@ -34,7 +34,7 @@ module read_gmsh use iso_c_binding use fldebug use global_parameters, only : OPTION_PATH_LEN, FIELD_NAME_LEN - use futils + use futils, only : copy_c_string_to_fortran use quadrature use elements use spud 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/tools/Flboundadapt.F90 b/tools/Flboundadapt.F90 index e3662b0ba0..28a404c436 100644 --- a/tools/Flboundadapt.F90 +++ b/tools/Flboundadapt.F90 @@ -121,8 +121,7 @@ end subroutine python_init call calculate_diagnostic_variables_new(states) ! Find the external mesh field - !call find_mesh_field_to_adapt(states(1), old_mesh_field) - old_mesh_field => extract_vector_field(states(1), "Coordinate") + old_mesh_field => extract_vector_field(states(1), topology_mesh_name) old_mesh => old_mesh_field%mesh ! Assemble the error metric @@ -156,32 +155,6 @@ end subroutine python_init call print_references(0) - ewrite(1, *) "Exiting fladapt" - -contains - - subroutine find_mesh_field_to_adapt(state, mesh_field) - !!< Find the external mesh field to be used by adaptivity - - type(state_type), intent(in) :: state - type(vector_field), pointer :: mesh_field - - character(len = FIELD_NAME_LEN) :: mesh_field_name - type(mesh_type), pointer :: mesh - - call find_mesh_to_adapt(state, mesh) - if(trim(mesh%name) == "CoordinateMesh") then - mesh_field_name = "Coordinate" - else - mesh_field_name = trim(mesh%name) // "Coordinate" - end if - - if(.not. has_vector_field(state, mesh_field_name)) then - FLAbort("External mesh field " // trim(mesh_field_name) // " not found in the system state") - end if - - mesh_field => extract_vector_field(state, mesh_field_name) - - end subroutine find_mesh_field_to_adapt - + ewrite(1, *) "Exiting flboundadapt" + end subroutine flboundadapt diff --git a/tools/Makefile.in b/tools/Makefile.in index 9cae1fe898..ff57f8a23c 100644 --- a/tools/Makefile.in +++ b/tools/Makefile.in @@ -95,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 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") From 3a4200667b33d11f6a856b160aad4af23550af98 Mon Sep 17 00:00:00 2001 From: James Percival Date: Wed, 18 Oct 2017 18:00:07 +0100 Subject: [PATCH 8/9] Plus bugfix. Squash me. --- tools/Flboundadapt.F90 | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/tools/Flboundadapt.F90 b/tools/Flboundadapt.F90 index 28a404c436..a7c84600d7 100644 --- a/tools/Flboundadapt.F90 +++ b/tools/Flboundadapt.F90 @@ -35,7 +35,8 @@ subroutine flboundadapt(input_flmlname_c, & use iso_c_binding use futils - use global_parameters, only: FIELD_NAME_LEN, OPTION_PATH_LEN + use global_parameters, only: FIELD_NAME_LEN, OPTION_PATH_LEN,& + topology_mesh_name use fldebug use parallel_tools, only: isparallel use reference_counting @@ -77,7 +78,7 @@ end subroutine python_init type(mesh_type), pointer :: old_mesh type(state_type), dimension(:), pointer :: states type(vector_field) :: new_mesh_field - type(vector_field), pointer :: new_mesh_field_ptr, old_mesh_field + type(vector_field), pointer :: new_mesh_field_ptr, position type(tensor_field) :: metric, t_edge_lengths character(len=FIELD_NAME_LEN) :: mesh_format character(len=OPTION_PATH_LEN) :: mesh_file_name @@ -121,18 +122,22 @@ end subroutine python_init call calculate_diagnostic_variables_new(states) ! Find the external mesh field - old_mesh_field => extract_vector_field(states(1), topology_mesh_name) - old_mesh => old_mesh_field%mesh + + ! 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") + old_mesh => extract_mesh(states(1), topology_mesh_name) ! Assemble the error metric call allocate(metric, old_mesh, "ErrorMetric") call assemble_metric(states, metric) - ewrite(0, *) "Expected nodes = ", expected_nodes(old_mesh_field, metric) + 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 = old_mesh_field, model = metric%mesh, & + call vtk_write_fields(trim(output_meshname) // "EdgeLengths", position = position, model = position%mesh, & & tfields = (/t_edge_lengths, metric/)) call deallocate(t_edge_lengths) @@ -141,7 +146,7 @@ end subroutine python_init new_mesh_field = read_gmsh_file(trim(input_geometryname),& format_string="geo", & quad_ngi=ele_ngi(old_mesh,1),& - position = old_mesh_field, metric = metric) + position = position, metric = metric) ! Write the output mesh call write_mesh_files(output_meshname, "gmsh", new_mesh_field) From 8e425e12b1e5b8cad3f2f468b455bceea1c40b22 Mon Sep 17 00:00:00 2001 From: James Percival Date: Wed, 25 Oct 2017 19:21:21 +0100 Subject: [PATCH 9/9] Run flboundadapt from a vtu with a metric in. --- diagnostics/Metric_Diagnostics.F90 | 39 +++++- diagnostics/Simple_Diagnostics.F90 | 37 +++++- femtools/GMSH_Common.F90 | 13 +- femtools/LibGmsh_integration.cpp | 192 ++++++++--------------------- schemas/diagnostic_algorithms.rnc | 40 ++++++ schemas/diagnostic_algorithms.rng | 58 +++++++++ tools/Flboundadapt.F90 | 83 ++++++++----- tools/Flboundadapt_main.cpp | 35 ++++-- 8 files changed, 306 insertions(+), 191 deletions(-) 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/GMSH_Common.F90 b/femtools/GMSH_Common.F90 index 9918ec0b83..09b64be1cf 100644 --- a/femtools/GMSH_Common.F90 +++ b/femtools/GMSH_Common.F90 @@ -385,18 +385,19 @@ subroutine position_to_gmodel(positions, gmodel) if (numFaces>0) sloc = face_loc(positions, 1) allocate(sndglno(numFaces*sloc)) - call getsndgln(positions%mesh, sndglno) + 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 - 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)) - else - bnd_ids = c_null_ptr + 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)) diff --git a/femtools/LibGmsh_integration.cpp b/femtools/LibGmsh_integration.cpp index 6233231855..abc9df08bb 100644 --- a/femtools/LibGmsh_integration.cpp +++ b/femtools/LibGmsh_integration.cpp @@ -28,22 +28,15 @@ #include "fmangle.h" #ifdef HAVE_LIBGMSH +#include "gmsh/Context.h" +#include "gmsh/Field.h" #include "gmsh/Gmsh.h" #include "gmsh/GModel.h" -#include "gmsh/Field.h" #include "gmsh/GEntity.h" #include "gmsh/MElement.h" -#include "gmsh/MTetrahedron.h" -#include "gmsh/MLine.h" -#include "gmsh/MTriangle.h" #include "gmsh/MVertex.h" -#include "gmsh/Context.h" #include "gmsh/PView.h" #include "gmsh/PViewDataGModel.h" -#include "gmsh/discreteVertex.h" -#include "gmsh/discreteEdge.h" -#include "gmsh/discreteFace.h" -#include "gmsh/discreteRegion.h" #include #include @@ -331,86 +324,16 @@ extern "C" { int &etype, int* eles, int &ftype, int* faces, int *ele_ids, int *face_ids, int *eleOwners) { - - gm = new GModel; - GEntity *e; - std::map point; - std::map edge; - std::map face; - std::map region; - - std::set uele_ids; - std::set uface_ids; - - uele_ids.insert(0); - if (ele_ids) { - for (int i=0;i::iterator it=uele_ids.begin(); - it != uele_ids.end(); ++it) { - edge[*it] = new discreteEdge(gm,*it,NULL,NULL); - edge[*it]->addPhysicalEntity(std::max(1,*it)); - edge[*it]->setTag(*it); - gm->add(edge[*it]); - } - e = edge[0]; - for (std::set::iterator it=uface_ids.begin(); - it != uface_ids.end(); ++it) { - point[*it] = new discreteVertex(gm,*it); - point[*it]->addPhysicalEntity(*it); - point[*it]->setTag(*it); - gm->add(point[*it]); - } - break; - case 2: - for (std::set::iterator it=uele_ids.begin(); - it != uele_ids.end(); ++it) { - face[*it] = new discreteFace(gm,*it); - face[*it]->addPhysicalEntity(std::max(1,*it)); - face[*it]->setTag(*it); - gm->add(face[*it]); - } - e=face[0]; - for (std::set::iterator it=uface_ids.begin(); - it != uface_ids.end(); ++it) { - edge[*it] = new discreteEdge(gm,*it,NULL,NULL); - edge[*it]->addPhysicalEntity(*it); - edge[*it]->setTag(*it); - gm->add(edge[*it]); - } - break; - case 3: - for (std::set::iterator it=uele_ids.begin(); - it != uele_ids.end(); ++it) { - region[*it] = new discreteRegion(gm,*it); - region[*it]->addPhysicalEntity(std::max(1,*it)); - region[*it]->setTag(*it); - gm->add(region[*it]); - } - e=region[0]; - for (std::set::iterator it=uface_ids.begin(); - it != uface_ids.end(); ++it) { - face[*it] = new discreteFace(gm,*it); - face[*it]->addPhysicalEntity(*it); - face[*it]->setTag(*it); - gm->add(face[*it]); - } - break; - } + std::map vertexMap; + std::map vertexInverseMap; + std::vector elementNum; + std::vector > vertexIndices; + std::vector elementType; + std::vector physical; + std::vector elementary; + std::vector partition; - MElementFactory f; double x[3] = {0,0,0}; for (int n=0;nsetIndex(n+1); - e->addMeshVertex(v); + vertexMap[n+1] = v; + vertexInverseMap[v->getNum()] = n+1; } - for (int n=0;ngetMeshVertex(eles[loc*n+ - GmshNodeOrdering(etype,i)]-1)); - } - MElement *ele = f.create(etype, vertices); - switch(etype) { - case MSH_LIN_2: - edge[id]->addLine((MLine*) ele); - break; - case MSH_TRI_3: - face[id]->addTriangle((MTriangle*) ele); - break; - case MSH_QUA_4: - face[id]->addQuadrangle((MQuadrangle*) ele); - break; - case MSH_TET_4: - region[id]->addTetrahedron((MTetrahedron*) ele); - break; - case MSH_HEX_8: - region[id]->addHexahedron((MHexahedron*) ele); - break; + for (int i=0; i< numElements; ++i) { + elementNum.push_back(i+1); + std::vector vertices; + for (int j=0; jgetMeshVertex(faces[sloc*n - +GmshNodeOrdering(ftype,i)]-1)); + for (int i=0; i< numFaces; ++i) { + elementNum.push_back(numElements+i+1); + std::vector vertices; + for (int j=0; jsetPartition(eleOwners[n]); - switch(ftype) { - case MSH_PNT: - point[id]->addPoint((MPoint*) ele); - break; - case MSH_LIN_2: - edge[id]->addLine((MLine*) ele); - break; - case MSH_TRI_3: - face[id]->addTriangle((MTriangle*) ele); - break; - case MSH_QUA_4: - face[id]->addQuadrangle((MQuadrangle*) ele); - break; + elementType.push_back(ftype); + vertexIndices.push_back(vertices); + physical.push_back( face_ids ? face_ids[i] : 1 ); + elementary.push_back( face_ids ? face_ids[i] : i+1); + partition.push_back( eleOwners ? eleOwners[i] : 0); + } + + gm = GModel::createGModel(vertexMap, elementNum, vertexIndices, + elementType, physical, elementary, partition); + EntityVector ents; + gm->getEntities(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); 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/tools/Flboundadapt.F90 b/tools/Flboundadapt.F90 index a7c84600d7..d8dca2720a 100644 --- a/tools/Flboundadapt.F90 +++ b/tools/Flboundadapt.F90 @@ -29,7 +29,7 @@ subroutine flboundadapt(input_flmlname_c, & input_geometryname_c, & - output_meshname_c) bind(c) + output_meshname_c, metric_name_c) bind(c) !!< Peforms a mesh adapt based on the supplied input options file. !!< Outputs the resulting mesh. @@ -69,25 +69,26 @@ end subroutine python_init end interface type(c_ptr), value :: input_flmlname_c, input_geometryname_c,& - output_meshname_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 - integer :: i, stat + 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 :: new_mesh_field_ptr, position - type(tensor_field) :: metric, t_edge_lengths - character(len=FIELD_NAME_LEN) :: mesh_format - character(len=OPTION_PATH_LEN) :: mesh_file_name - real, dimension(:,:,:), allocatable :: metric_data + 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" @@ -95,43 +96,58 @@ end subroutine python_init call python_init() #endif - ewrite(2, *) "Input flml name: " // trim(input_flmlname) + 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) - - ! 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() + 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" + ewrite(1, *) "Performing options sanity check" + call check_options() + ewrite(1, *) "Options sanity check successful" #endif - ! Populate the system state - call populate_state(states) +! Populate the system state + call populate_state(states) ! Calculate diagnostic fields - call calculate_diagnostic_variables(states) - call calculate_diagnostic_variables_new(states) - + 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") - old_mesh => extract_mesh(states(1), topology_mesh_name) + if (have_metric_name) then + metric => extract_tensor_field(states(1), metric_name) + old_mesh => metric%mesh + else ! Assemble the error metric - call allocate(metric, old_mesh, "ErrorMetric") - call assemble_metric(states, 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) @@ -153,9 +169,14 @@ end subroutine python_init ! Deallocate do i = 1, size(states) - call deallocate(states(i)) + call deallocate(states(i)) end do - call deallocate(metric) + if (have_metric_name) then + deallocate(states) + else + call deallocate(metric) + deallocate(metric) + end if call deallocate(new_mesh_field) call print_references(0) diff --git a/tools/Flboundadapt_main.cpp b/tools/Flboundadapt_main.cpp index dafad8ca1d..c0f915e0f7 100644 --- a/tools/Flboundadapt_main.cpp +++ b/tools/Flboundadapt_main.cpp @@ -42,7 +42,7 @@ using namespace std; extern "C"{ - void flboundadapt(const char *, const char *, const char *); + void flboundadapt(const char *, const char *, const char *, const char *); #ifdef HAVE_PYTHON #include "python_statec.h" #endif @@ -50,6 +50,8 @@ extern "C"{ 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" @@ -57,6 +59,7 @@ void Usage(){ << "Options:\n" << "\n" << "-h\t\tDisplay this help\n" + << "-i\t\t read from VTU file\n" << "-v\t\tVerbose mode" << endl; } @@ -78,8 +81,10 @@ int main(int argc, char** argv){ optarg = NULL; char c; map args; - while((c = getopt(argc, argv, "hv")) != -1){ + 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{ @@ -110,15 +115,18 @@ int main(int argc, char** argv){ set_global_debug_level_fc(&verbosity); // Input and output base names - string input_flmlname, input_geometryname, output_meshname; - if(argc > optind + 3){ - input_flmlname = argv[optind + 1]; - input_geometryname = argv[optind + 2]; - output_meshname = argv[optind + 3]; - }else if(argc == optind + 3){ - input_flmlname = argv[optind]; - input_geometryname = argv[optind+1]; - output_meshname = argv[optind + 2]; + 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); @@ -126,9 +134,10 @@ int main(int argc, char** argv){ - flboundadapt(input_flmlname.c_str(), + flboundadapt(input_name.c_str(), input_geometryname.c_str(), - output_meshname.c_str()); + output_meshname.c_str(), + metric_name.c_str()); #ifdef HAVE_PYTHON // Finalize the Python Interpreter