reading via meshio

libsrc/meshing/python_mesh.cpp

@@ -913,6 +913,89 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
           {
             return self.AddFaceDescriptor (fd);
           })
+
+    .def ("AddPoints", [](Mesh & self, py::buffer b)
+          {
+            py::buffer_info info = b.request();
+            if (info.ndim != 2)
+              throw std::runtime_error("AddPoints needs buffer of dimension 2");
+            if (info.format != py::format_descriptor<double>::format())
+              throw std::runtime_error("AddPoints needs buffer of type double");
+            if (info.strides[0] != sizeof(double)*info.shape[1])
+              throw std::runtime_error("AddPoints needs packed array");              
+            double * ptr = static_cast<double*> (info.ptr);
+            if (info.shape[1]==2)
+              for (auto i : Range(info.shape[0]))
+                {
+                  self.AddPoint (Point<3>(ptr[0], ptr[1], 0));
+                  ptr += 2;
+                }
+            if (info.shape[1]==3)
+              for (auto i : Range(info.shape[0]))
+                {
+                  self.AddPoint (Point<3>(ptr[0], ptr[1], ptr[2]));
+                  ptr += 3;
+                }
+          })
+    .def ("AddElements", [](Mesh & self, int dim, int index, py::buffer b, int base)
+          {
+            py::buffer_info info = b.request();
+            if (info.ndim != 2)
+              throw std::runtime_error("AddElements needs buffer of dimension 2");
+            if (info.format != py::format_descriptor<int>::format())
+              throw std::runtime_error("AddPoints needs buffer of type int");
+
+            int * ptr = static_cast<int*> (info.ptr);
+            if (dim == 2)
+              {
+                ELEMENT_TYPE type;
+                int np = info.shape[1];
+                switch (np)
+                  {
+                  case 3: type = TRIG; break;
+                  case 4: type = QUAD; break;
+                  case 6: type = TRIG6; break;
+                  case 8: type = QUAD8; break;
+                  default:
+                    throw Exception("unsupported 2D element with "+ToString(np)+" points");
+                  }
+                for (auto i : Range(info.shape[0]))
+                  {
+                    Element2d el(type);
+                    for (int j = 0; j < np;j ++)
+                      el[j] = ptr[j]+PointIndex::BASE-base;
+                    el.SetIndex(index);
+                    self.AddSurfaceElement (el);
+                    ptr += info.strides[0]/sizeof(int);
+                  }
+              }
+            if (dim == 3)
+              {
+                ELEMENT_TYPE type;
+                int np = info.shape[1];
+                switch (np)
+                  {
+                  case 4: type = TET; break;
+                    /* // have to check ordering of points
+                       case 10: type = TET10; break;
+                       case 8: type = HEX; break;
+                       case 6: type = PRISM; break;
+                    */
+                  default:
+                    throw Exception("unsupported 3D element with "+ToString(np)+" points");
+                  }
+                for (auto i : Range(info.shape[0]))
+                  {
+                    Element el(type);
+                    for (int j = 0; j < np;j ++)
+                      el[j] = ptr[j]+PointIndex::BASE-base;
+                    el.SetIndex(index);
+                    self.AddVolumeElement (el);
+                    ptr += info.strides[0]/sizeof(int);
+                  }
+              }
+            
+          }, py::arg("dim"), py::arg("index"), py::arg("data"), py::arg("base")=0)
     
     .def ("DeleteSurfaceElement",
           [](Mesh & self, SurfaceElementIndex i)

libsrc/meshing/topology.cpp

@@ -1559,7 +1559,7 @@ namespace netgen
 	      }
 
 	    if (cnt_err && ntasks == 1)
-	      cout << cnt_err << " elements are not matching !!!" << endl;
+	      cout << IM(5) << cnt_err << " elements are not matching !!!" << endl;
 	  }
         // NgProfiler::StopTimer (timer2c);
 

python/CMakeLists.txt

@@ -11,7 +11,8 @@ install(FILES
         ${CMAKE_CURRENT_BINARY_DIR}/config.py
         ${CMAKE_CURRENT_BINARY_DIR}/version.py
         __main__.py __init__.py
-        meshing.py csg.py geom2d.py stl.py gui.py NgOCC.py occ.py read_gmsh.py
+        meshing.py csg.py geom2d.py stl.py gui.py NgOCC.py occ.py
+        read_gmsh.py read_meshio.py
         webgui.py
         DESTINATION ${NG_INSTALL_DIR_PYTHON}/${NG_INSTALL_SUFFIX}
         COMPONENT netgen

python/read_meshio.py

@@ -0,0 +1,20 @@
+from netgen.meshing import *
+
+def ReadViaMeshIO(filename):
+    import meshio
+    import numpy as np
+    
+    # print ("reading via meshio:", filename)
+
+    m = meshio.read(filename)
+    pts = m.points
+    
+    mesh = Mesh(dim=pts.shape[1])
+    mesh.AddPoints ( np.asarray(pts, dtype=np.float64) )  # needed for correct little/big endian
+
+    fd = mesh.Add (FaceDescriptor(bc=1))
+    for cb in m.cells:
+        mesh.AddElements(dim=cb.dim, index=1, data=np.asarray(cb.data, dtype=np.int32), base=0)
+            
+    return mesh
+