From db5ad09b790f4522c887f437dc01d9e5e1e1f734 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Wed, 11 Mar 2020 11:48:05 +0100 Subject: [PATCH] Initial CGNS format read support --- CMakeLists.txt | 10 + cmake/SuperBuild.cmake | 1 + libsrc/interface/CMakeLists.txt | 4 +- libsrc/interface/readuser.cpp | 7 + libsrc/interface/rw_cgns.cpp | 466 ++++++++++++++++++++++++++++++++ libsrc/interface/writeuser.hpp | 9 + libsrc/meshing/python_mesh.cpp | 1 + ng/menustat.tcl | 1 + ng/onetcl.cpp | 1 + tests/build_debug.sh | 8 +- tests/dockerfile | 2 +- 11 files changed, 506 insertions(+), 4 deletions(-) create mode 100644 libsrc/interface/rw_cgns.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 009660aa..8f11ccca 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -17,6 +17,7 @@ option( USE_MPI "enable mpi parallelization" OFF ) option( USE_OCC "(not supported) compile with OpenCascade geometry kernel" OFF) option( USE_JPEG "enable snapshots using library libjpeg" OFF ) option( USE_MPEG "enable video recording with FFmpeg, uses libavcodec" OFF ) +option( USE_CGNS "enable CGNS file read/write support" OFF ) option( INTEL_MIC "cross compile for intel xeon phi") option( INSTALL_PROFILES "install environment variable settings to /etc/profile.d" OFF ) option( USE_CCACHE "use ccache") @@ -388,6 +389,15 @@ if(ENABLE_CPP_CORE_GUIDELINES_CHECK) endif() endif(ENABLE_CPP_CORE_GUIDELINES_CHECK) +add_library(netgen_cgns INTERFACE) +if(USE_CGNS) + find_library( CGNS_LIBRARY cgns ) + find_path( CGNS_INCLUDE_DIR cgnslib.h ) + target_compile_definitions(netgen_cgns INTERFACE NG_CGNS) + target_include_directories(netgen_cgns INTERFACE ${CGNS_INCLUDE_DIR}) + target_link_libraries(netgen_cgns INTERFACE ${CGNS_LIBRARY}) +endif(USE_CGNS) + add_subdirectory(libsrc) add_subdirectory(ng) add_subdirectory(tutorials) diff --git a/cmake/SuperBuild.cmake b/cmake/SuperBuild.cmake index a04e4164..83dfe269 100644 --- a/cmake/SuperBuild.cmake +++ b/cmake/SuperBuild.cmake @@ -128,6 +128,7 @@ set_vars( NETGEN_CMAKE_ARGS USE_OCC USE_MPEG USE_JPEG + USE_CGNS USE_INTERNAL_TCL INSTALL_PROFILES INTEL_MIC diff --git a/libsrc/interface/CMakeLists.txt b/libsrc/interface/CMakeLists.txt index bc1bc2f3..47661a0c 100644 --- a/libsrc/interface/CMakeLists.txt +++ b/libsrc/interface/CMakeLists.txt @@ -4,10 +4,10 @@ add_library(interface ${NG_LIB_TYPE} read_fnf_mesh.cpp readtetmesh.cpp readuser.cpp writeabaqus.cpp writediffpack.cpp writedolfin.cpp writeelmer.cpp writefeap.cpp writefluent.cpp writegmsh.cpp writejcm.cpp writepermas.cpp writetecplot.cpp writetet.cpp writetochnog.cpp writeuser.cpp - wuchemnitz.cpp writegmsh2.cpp writeOpenFOAM15x.cpp + wuchemnitz.cpp writegmsh2.cpp writeOpenFOAM15x.cpp rw_cgns.cpp ) -target_link_libraries(interface mesh csg geom2d stl visual) +target_link_libraries(interface PUBLIC mesh csg geom2d stl visual PRIVATE netgen_cgns) if(NOT WIN32) install( TARGETS interface ${NG_INSTALL_DIR}) diff --git a/libsrc/interface/readuser.cpp b/libsrc/interface/readuser.cpp index 472592ac..d3d50cf2 100644 --- a/libsrc/interface/readuser.cpp +++ b/libsrc/interface/readuser.cpp @@ -649,6 +649,13 @@ namespace netgen ReadFNFFormat (mesh, filename); } + // .cgns file - CFD General Notation System + if ( (strlen (filename) > 5) && + strcmp (&filename[strlen (filename)-5], ".cgns") == 0 ) + { + ReadCGNSMesh (mesh, filename); + } + if ( ( (strlen (filename) > 4) && strcmp (&filename[strlen (filename)-4], ".stl") == 0 ) || ( (strlen (filename) > 5) && strcmp (&filename[strlen (filename)-5], ".stlb") == 0 ) ) { diff --git a/libsrc/interface/rw_cgns.cpp b/libsrc/interface/rw_cgns.cpp new file mode 100644 index 00000000..7de2aefb --- /dev/null +++ b/libsrc/interface/rw_cgns.cpp @@ -0,0 +1,466 @@ +#include +#include "writeuser.hpp" + +#ifdef NG_CGNS +#include + +#include + +namespace netgen::cg +{ + int getDim(ElementType_t type) + { + switch(type) + { + case BAR_2: + case BAR_3: + return 1; + case TRI_3: + case TRI_6: + case QUAD_4: + case QUAD_8: + return 2; + case TETRA_4: + case TETRA_10: + case PYRA_5: + case PYRA_13: + case HEXA_8: + case HEXA_20: + case PENTA_6: + case PENTA_15: + return 3; + default: + throw Exception("Read CGNS: unknown element type " + string(ElementTypeName[type])); + } + } + + Segment ReadCGNSElement1D( ElementType_t type, FlatArray verts, int vert_offset=0 ) + { + int np; + cg_npe(type, &np); + + Segment s; + for (auto i : Range(np)) + s[i] = vert_offset+verts[i]; + return s; + } + + Element2d ReadCGNSElement2D( ElementType_t type, FlatArray verts, int vert_offset=0 ) + { + static constexpr int map_tri3[] = {0,2,1}; + static constexpr int map_tri6[] = {0,2,1,3,5,4}; // untested + static constexpr int map_quad4[] = {0,3,2,1}; + static constexpr int map_quad8[] = {0,3,2,1,4,7,6,5}; // untested + + const int * map = nullptr; + switch(type) + { + case TRI_3: + map = map_tri3; + break; + case QUAD_4: + map = map_quad4; + break; + case TRI_6: + map = map_tri6; + break; + case QUAD_8: + map = map_quad8; + break; + default: + throw Exception("Read CGNS: unknown element type " + string(ElementTypeName[type])); + } + + int np; + cg_npe(type, &np); + + Element2d el(np); + for (auto i : Range(np)) + el[i] = vert_offset+verts[map[i]]; + return el; + } + + Element ReadCGNSElement3D( ElementType_t type, FlatArray verts, int vert_offset=0 ) + { + static constexpr int map_tet4[] = {0,2,1,3}; + static constexpr int map_prism6[] = {0,2,1,3,5,4}; + static constexpr int map_pyra5[] = {0,3,2,1,4}; + static constexpr int map_hexa8[] = {0,3,2,1,4,7,6,5}; + int np; + cg_npe(type, &np); + + const int * map = nullptr; + switch(type) + { + case TETRA_4: + map = map_tet4; break; + case PYRA_5: + map = map_pyra5; break; + case PENTA_6: + map = map_prism6; break; + case HEXA_8: + map = map_hexa8; break; + // TODO: Second order elements + case TETRA_10: + case PYRA_13: + case HEXA_20: + case PENTA_15: + default: + throw Exception("Read CGNS: unknown element type " + string(ElementTypeName[type])); + } + + Element el(np); + for (auto i : Range(np)) + el[i] = vert_offset+verts[map[i]]; + return el; + } + + // maps cgns node type to ngsolve node type + // enum NODE_TYPE { NT_VERTEX = 0, NT_EDGE = 1, NT_FACE = 2, NT_CELL = 3, NT_ELEMENT = 4, NT_FACET = 5 }; + int getNodeType( GridLocation_t location ) + { + switch(location) + { + case Vertex: + return 0; + case CellCenter: + return 3; + case FaceCenter: + return 2; + case EdgeCenter: + return 1; + default: + throw Exception("Read CGNS: unknown grid location " + string(GridLocationName[location])); + } + } + + + struct Solution + { + int fn, base, zone, solution; + string name; + GridLocation_t location; // solution is defined on either cells, faces, edges or vertices + PointSetType_t point_type; + cgsize_t n_points; + + Array field_names; + Array field_datatypes; + + Solution() = default; + + Solution(int fn_, int base_, int zone_, int solution_) + : fn(fn_), base(base_), zone(zone_), solution(solution_) + { + char solname[100]; + cg_sol_info(fn, base, zone, solution, solname, &location); + name = solname; + cg_sol_ptset_info(fn, base, zone, solution, &point_type, &n_points); + + int n_fields = 0; + cg_nfields(fn, base, zone, solution, &n_fields); + + field_names.SetSize(n_fields); + field_datatypes.SetSize(n_fields); + for(auto fi : Range(n_fields)) + { + char buf[100]; + cg_field_info(fn, base, zone, solution, fi+1, &field_datatypes[fi], buf); + field_names[fi] = buf; + } + } + }; + + struct Zone + { + ZoneType_t zone_type; + int fn, base, zone; + int nv, ne, first_vertex, first_mat, first_bc; + Array materials; + Array boundaries; + string name; + cgsize_t size[3]; + + Array solutions; + + Zone(int fn_, int base_, int zone_) + : fn(fn_), base(base_), zone(zone_) + { + cg_zone_type(fn, base, zone, &zone_type); + char zone_name[100]; + cg_zone_read(fn,base,zone, zone_name, size); + nv = size[0]; + + int n_solutions; + cg_nsols(fn, base, zone, &n_solutions); + + solutions.SetSize(n_solutions); + for(auto si : Range(n_solutions)) + solutions[si] = Solution{fn, base, zone, si+1}; + } + + void ReadSolutions( std::vector & sol_names, std::vector> & sol_values, std::vector & sol_locations ) + { + static Timer tall("CGNS::ReadSolutions"); RegionTimer rtall(tall); + for (auto & sol : solutions) + { + for (auto fi : Range(sol.field_names.Size())) + { + cgsize_t size = sol.n_points; + if(size==0) + { + switch(sol.location) + { + case Vertex: + size = nv; + break; + case CellCenter: + size = ne; + break; + case FaceCenter: + case IFaceCenter: + case JFaceCenter: + case KFaceCenter: + case EdgeCenter: + default: + throw Exception("Read CGNS: unknown grid location " + string(GridLocationName[sol.location])); + } + } + + size = size==0 ? nv : size; + auto values = Array(size); + + cgsize_t imin = 1UL; + cg_field_read(fn, base, zone, sol.solution, sol.field_names[fi].c_str(), RealDouble, &imin, &size, &values[0]); + sol_names.push_back(sol.field_names[fi]); + sol_values.emplace_back(std::move(values)); + sol_locations.push_back(sol.location); + } + } + } + + void ReadMesh( Mesh & mesh ) + { + static Timer tall("CGNS::ReadMesh-Zone"); RegionTimer rtall(tall); + static Timer tsection("CGNS::ReadMesh-Section"); + first_vertex = mesh.GetNP(); + first_mat = mesh.GetRegionNamesCD(0).Size(); + first_bc = mesh.GetRegionNamesCD(1).Size(); + ne = 0; + + Array x(nv), y(nv), z(nv); + cgsize_t imin=1; + cg_coord_read(fn,base,zone, "CoordinateX", RealSingle, &imin, &nv, &x[0]); + cg_coord_read(fn,base,zone, "CoordinateY", RealSingle, &imin, &nv, &y[0]); + cg_coord_read(fn,base,zone, "CoordinateZ", RealSingle, &imin, &nv, &z[0]); + + for(auto i : Range(nv)) + mesh.AddPoint( {x[i], y[i], z[i]} ); + + int nsections; + cg_nsections(fn, base, zone, &nsections); + + int bc = first_bc; + int material = first_mat; + for (auto section : Range(1,nsections+1)) + { + RegionTimer rtsection(tsection); + char name[100]; + ElementType_t type; + cgsize_t start, end; + int nbndry, parent_flag; + + cg_section_read(fn, base, zone, section, name, &type, &start, &end, &nbndry, &parent_flag); + PrintMessage(4, "Read section ", section, " with name ", name, " and element type ", cg_ElementTypeName(type)); + + string ngname{name}; + + for (char & c : ngname) + if(c==' ') + c = '_'; + + + if(type==MIXED) + { + bc++; + material++; + mesh.AddFaceDescriptor(FaceDescriptor(bc, 1, 0, 1)); + mesh.SetBCName(bc-1, ngname); + mesh.SetMaterial(material, ngname); + + cgsize_t nv; + cg_ElementDataSize(fn, base, zone, section, &nv); + + Array vertices(nv); + cg_poly_elements_read(fn, base, zone, section, &vertices[0], nullptr, nullptr); + + size_t vi = 0; + while(vi(vertices[vi++]); + int dim = getDim(type); + + if(dim==1) + { + auto el = ReadCGNSElement1D(type, vertices.Range(vi, vertices.Size()), first_vertex); + mesh.AddSegment(el); + vi += el.GetNP(); + } + + if(dim==2) + { + auto el = ReadCGNSElement2D(type, vertices.Range(vi, vertices.Size()), first_vertex); + el.SetIndex(bc); + mesh.AddSurfaceElement(el); + vi += el.GetNP(); + } + + if(dim==3) + { + auto el = ReadCGNSElement3D(type, vertices.Range(vi, vertices.Size()), first_vertex); + el.SetIndex(material); + mesh.AddVolumeElement(el); + vi += el.GetNP(); + ne++; + } + } + } + else + { + int dim = getDim(type); + + cgsize_t nv; + cg_ElementDataSize(fn, base, zone, section, &nv); + int np=0; + cg_npe(type, &np); + + Array vertices(nv); + cg_elements_read(fn, base, zone, section, &vertices[0], nullptr); + int ne_section = nv/np; + + if(dim==1) + { + for(auto i : Range(ne_section)) + { + auto el = ReadCGNSElement1D(type, vertices.Range(np*i, np*(i+1)), first_vertex); + mesh.AddSegment(el); + } + } + + if(dim==2) + { + bc++; + mesh.AddFaceDescriptor(FaceDescriptor(bc, 1, 0, 1)); + for(auto i : Range(ne_section)) + { + auto el = ReadCGNSElement2D(type, vertices.Range(np*i, np*(i+1)), first_vertex); + el.SetIndex(bc); + mesh.AddSurfaceElement(el); + } + mesh.SetBCName(bc-1, ngname); + } + + if(dim==3) + { + material++; + for(auto i : Range(ne_section)) + { + auto el = ReadCGNSElement3D(type, vertices.Range(np*i, np*(i+1)), first_vertex); + el.SetIndex(material); + mesh.AddVolumeElement(el); + } + mesh.SetMaterial(material, ngname); + ne += ne_section; + } + } + } + } + }; +} + +namespace netgen +{ + void ReadCGNSMesh (Mesh & mesh, const string & filename) + { + static Timer tall("CGNS::ReadMesh"); RegionTimer rtall(tall); + int fn; + cg_open(filename.c_str(),CG_MODE_READ,&fn); + + int base = 1; + int nzones; + cg_nzones(fn, base, &nzones); + + int bc = 0; + int material = 0; + + for (auto zi : Range(1, nzones+1)) + { + ZoneType_t zone_type; + cg_zone_type(fn, base, zi, &zone_type); + if(zone_type != Unstructured ) + { + PrintMessage(2, "skipping zone with type ", cg_ZoneTypeName(zone_type) ); + continue; + } + cg::Zone zone(fn, base, zi); + zone.ReadMesh( mesh ); + } + } + + // Reads mesh and solutions of .csns file + tuple, vector, vector>, vector> ReadCGNSFile(string filename, int base) + { + static Timer tall("CGNS::ReadFile"); RegionTimer rtall(tall); + int fn; + cg_open(filename.c_str(),CG_MODE_READ,&fn); + + int nbases; + cg_nbases(fn, &nbases); + + int nzones; + cg_nzones(fn, base, &nzones); + + auto mesh = make_shared(); + + int bc = 0; + int material = 0; + + + std::vector names; + std::vector> values; + std::vector locations; + + for (auto zi : Range(1, nzones+1)) + { + ZoneType_t zone_type; + cg_zone_type(fn, base, zi, &zone_type); + if(zone_type != Unstructured ) + { + clog << "skipping zone with type " << cg_ZoneTypeName(zone_type) << endl; + continue; + } + cg::Zone zone(fn, base, zi); + zone.ReadMesh( *mesh ); + zone.ReadSolutions( names, values, locations ); + } + + cg_close(fn); + return std::make_tuple(mesh, names, values, locations); + } +} + +#else // NG_CGNS + +namespace netgen +{ + void ReadCGNSMesh (Mesh & mesh, const string & filename) + { + PrintMessage(1, "Could not import CGNS mesh: Netgen was built without CGNS support"); + } + + tuple, vector, vector>, vector> ReadCGNSFile(string filename, int base) + { + throw Exception("Netgen was built without CGNS support"); + } +} + +#endif // NG_CGNS diff --git a/libsrc/interface/writeuser.hpp b/libsrc/interface/writeuser.hpp index 3e98c433..8c58ea5f 100644 --- a/libsrc/interface/writeuser.hpp +++ b/libsrc/interface/writeuser.hpp @@ -150,6 +150,15 @@ extern void ReadFNFFormat (Mesh & mesh, +extern void DLL_HEADER ReadCGNSMesh (Mesh & mesh, + const string & filename); + +// Read Mesh and solutions from CGNS file +extern tuple, vector, vector>, vector> +DLL_HEADER ReadCGNSFile(string filename, int base); + + + void WriteDolfinFormat (const Mesh & mesh, const string & filename); diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 45adc9fd..23cb5d3e 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -1069,6 +1069,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) })); + m.def("ReadCGNSFile", &ReadCGNSFile, py::arg("filename"), py::arg("base")=1, "Read mesh and solution vectors from CGNS file"); } PYBIND11_MODULE(libmesh, m) { diff --git a/ng/menustat.tcl b/ng/menustat.tcl index 6f870a81..329d59bd 100644 --- a/ng/menustat.tcl +++ b/ng/menustat.tcl @@ -233,6 +233,7 @@ loadmeshinifile; {"TET format" {.tet} } {"STL format" {.stl .stlb} } {"Pro/ENGINEER neutral format" {.fnf} } + {"CFD General Notation System" {.cgns} } } set file [tk_getOpenFile -filetypes $types ] if {$file != ""} { diff --git a/ng/onetcl.cpp b/ng/onetcl.cpp index a0a31da4..bb98f97f 100644 --- a/ng/onetcl.cpp +++ b/ng/onetcl.cpp @@ -868,6 +868,7 @@ DLL_HEADER const char * ngscript[] = {"" ,"{\"TET format\" {.tet} }\n" ,"{\"STL format\" {.stl .stlb} }\n" ,"{\"Pro/ENGINEER neutral format\" {.fnf} }\n" +,"{\"CFD General Notation System\" {.cgns} }\n" ,"}\n" ,"set file [tk_getOpenFile -filetypes $types ]\n" ,"if {$file != \"\"} {\n" diff --git a/tests/build_debug.sh b/tests/build_debug.sh index fe650b8e..349722d8 100755 --- a/tests/build_debug.sh +++ b/tests/build_debug.sh @@ -1,6 +1,12 @@ cd mkdir -p build/netgen cd build/netgen -cmake ../../src/netgen -DUSE_CCACHE=ON -DBUILD_TYPE=DEBUG -DENABLE_UNIT_TESTS=ON -DUSE_OCC=ON +cmake \ + -DUSE_CCACHE=ON \ + -DBUILD_TYPE=DEBUG \ + -DENABLE_UNIT_TESTS=ON \ + -DUSE_OCC=ON \ + -DUSE_CGNS=ON \ + ../../src/netgen make -j12 make install diff --git a/tests/dockerfile b/tests/dockerfile index 3415481f..e14f64ff 100644 --- a/tests/dockerfile +++ b/tests/dockerfile @@ -1,5 +1,5 @@ FROM ubuntu:19.10 ENV DEBIAN_FRONTEND=noninteractive MAINTAINER Matthias Hochsteger -RUN apt-get update && apt-get -y install python3 libpython3-dev libxmu-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev ccache python3-pytest python3-numpy python3-tk clang-tidy python3-distutils clang libocct-data-exchange-dev +RUN apt-get update && apt-get -y install python3 libpython3-dev libxmu-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev ccache python3-pytest python3-numpy python3-tk clang-tidy python3-distutils clang libocct-data-exchange-dev libcgns-dev ADD . /root/src/netgen