Initial CGNS format read support

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)

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

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})

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 ) )
       {

libsrc/interface/rw_cgns.cpp

@@ -0,0 +1,466 @@
+#include <meshing.hpp>
+#include "writeuser.hpp"
+
+#ifdef NG_CGNS
+#include <variant>
+
+#include <cgnslib.h>
+
+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<cgsize_t> 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<cgsize_t> 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<cgsize_t> 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<string> field_names;
+      Array<DataType_t> 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<int> materials;
+      Array<int> boundaries;
+      string name;
+      cgsize_t size[3];
+
+      Array<Solution> 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<string> & sol_names, std::vector<Array<double>> & sol_values, std::vector<int> & 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<double>(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<float> 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<cgsize_t> vertices(nv);
+                  cg_poly_elements_read(fn, base, zone, section, &vertices[0], nullptr, nullptr);
+
+                  size_t vi = 0;
+                  while(vi<nv)
+                    {
+                      auto type = static_cast<ElementType_t>(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<cgsize_t> 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<shared_ptr<Mesh>, vector<string>, vector<Array<double>>, vector<int>> 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<Mesh>();
+
+      int bc = 0;
+      int material = 0;
+
+
+      std::vector<string> names;
+      std::vector<Array<double>> values;
+      std::vector<int> 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<shared_ptr<Mesh>, vector<string>, vector<Array<double>>, vector<int>> ReadCGNSFile(string filename, int base)
+    {
+      throw Exception("Netgen was built without CGNS support");
+    }
+}
+
+#endif // NG_CGNS

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<shared_ptr<Mesh>, vector<string>, vector<Array<double>>, vector<int>>
+DLL_HEADER ReadCGNSFile(string filename, int base);
+
+
+
 void WriteDolfinFormat (const Mesh & mesh,
                         const string & filename);
 

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) {

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 != ""} {

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"

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

tests/dockerfile

@@ -1,5 +1,5 @@
 FROM ubuntu:19.10
 ENV DEBIAN_FRONTEND=noninteractive
 MAINTAINER Matthias Hochsteger <matthias.hochsteger@tuwien.ac.at>
-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