Read mesh domains from vtk file 363

Code/Source/solver/VtkData.cpp

@@ -629,6 +629,30 @@ Array<int> VtkVtpData::get_connectivity() const
   return  conn;
 }
 
+/// @brief Copy an array of cell data from an polydata mesh into the given Vector.
+//
+void VtkVtpData::copy_cell_data(const std::string& data_name, Vector<int>& mesh_data)
+{
+  auto vtk_data = vtkIntArray::SafeDownCast(impl->vtk_polydata->GetCellData()->GetArray(data_name.c_str()));
+  if (vtk_data == nullptr) { 
+    return;
+  }
+
+  int num_data = vtk_data->GetNumberOfTuples();
+  if (num_data == 0) {
+    return; 
+  }       
+
+  int num_comp = vtk_data->GetNumberOfComponents();
+
+  for (int i = 0; i < num_data; i++) {
+    auto tuple = vtk_data->GetTuple(i);
+    for (int j = 0; j < num_comp; j++) {
+      mesh_data(i) = tuple[j];
+    }
+  }
+}
+
 /// @brief Copy an array of point data from an polydata mesh into the given Array.
 //
 void VtkVtpData::copy_point_data(const std::string& data_name, Array<double>& mesh_data)
@@ -774,6 +798,19 @@ Array<double> VtkVtpData::get_points() const
   return points_array;
 }
 
+bool VtkVtpData::has_cell_data(const std::string& data_name)
+{
+  int num_arrays = impl->vtk_polydata->GetPointData()->GetNumberOfArrays();
+
+  for (int i = 0; i < num_arrays; i++) {
+    if (!strcmp(impl->vtk_polydata->GetCellData()->GetArrayName(i), data_name.c_str())) {
+      return true;
+    }
+  }
+
+  return false; 
+}
+
 bool VtkVtpData::has_point_data(const std::string& data_name)
 {
   int num_arrays = impl->vtk_polydata->GetPointData()->GetNumberOfArrays();
@@ -911,6 +948,30 @@ std::vector<std::string> VtkVtuData::get_point_data_names()
   return data_names;
 }
 
+/// @brief Copy an array of cell data from an unstructured mesh into the given Vector.
+//
+void VtkVtuData::copy_cell_data(const std::string& data_name, Vector<int>& mesh_data)
+{
+  auto vtk_data = vtkIntArray::SafeDownCast(impl->vtk_ugrid->GetCellData()->GetArray(data_name.c_str()));
+  if (vtk_data == nullptr) {
+    return;
+  }
+
+  int num_data = vtk_data->GetNumberOfTuples();
+  if (num_data == 0) {
+    return;
+  }
+
+  int num_comp = vtk_data->GetNumberOfComponents();
+
+  for (int i = 0; i < num_data; i++) {
+    auto tuple = vtk_data->GetTuple(i);
+    for (int j = 0; j < num_comp; j++) {
+      mesh_data(i) = tuple[j];
+    }
+  }
+}
+
 /// @brief Copy an array of point data from an unstructured grid into the given Array.
 //
 void VtkVtuData::copy_point_data(const std::string& data_name, Array<double>& mesh_data)
@@ -996,6 +1057,19 @@ void VtkVtuData::copy_points(Array<double>& points)
   return;
 }
 
+bool VtkVtuData::has_cell_data(const std::string& data_name)
+{
+  int num_arrays = impl->vtk_ugrid->GetPointData()->GetNumberOfArrays();
+
+  for (int i = 0; i < num_arrays; i++) {
+    if (!strcmp(impl->vtk_ugrid->GetCellData()->GetArrayName(i), data_name.c_str())) {
+      return true;
+    }
+  }
+
+  return false;
+}
+
 bool VtkVtuData::has_point_data(const std::string& data_name)
 {
   int num_arrays = impl->vtk_ugrid->GetPointData()->GetNumberOfArrays();

Code/Source/solver/VtkData.h

@@ -32,11 +32,20 @@ class VtkData {
     virtual void set_points(const Array<double>& points) = 0;
     virtual void set_connectivity(const int nsd, const Array<int>& conn, const int pid = 0) = 0;
 
+    virtual bool has_cell_data(const std::string& data_name) = 0;
     virtual bool has_point_data(const std::string& data_name) = 0;
 
     virtual void copy_points(Array<double>& points) = 0;
+
     virtual void copy_point_data(const std::string& data_name, Array<double>& mesh_data) = 0;
     virtual void copy_point_data(const std::string& data_name, Vector<double>& mesh_data) = 0;
+    virtual void copy_point_data(const std::string& data_name, Vector<int>& mesh_data) = 0;
+
+    virtual void copy_cell_data(const std::string& data_name, Vector<int>& mesh_data) = 0;
+
+    virtual Array<double> get_point_data(const std::string& data_name) = 0;
+    virtual std::vector<std::string> get_point_data_names() = 0;
+
     virtual void write() = 0;
 
     static VtkData* create_reader(const std::string& file_name);
@@ -64,13 +73,20 @@ class VtkVtpData : public VtkData {
     virtual int num_points() const override;
     virtual void read_file(const std::string& file_name) override;
 
-    void copy_points(Array<double>& points) override;
-    void copy_point_data(const std::string& data_name, Array<double>& mesh_data) override;
-    void copy_point_data(const std::string& data_name, Vector<double>& mesh_data) override;
-    void copy_point_data(const std::string& data_name, Vector<int>& mesh_data);
-    Array<double> get_point_data(const std::string& data_name);
-    std::vector<std::string> get_point_data_names();
-    bool has_point_data(const std::string& data_name) override;
+    virtual void copy_points(Array<double>& points) override;
+
+    virtual void copy_point_data(const std::string& data_name, Array<double>& mesh_data) override;
+    virtual void copy_point_data(const std::string& data_name, Vector<double>& mesh_data) override;
+    virtual void copy_point_data(const std::string& data_name, Vector<int>& mesh_data) override;
+
+    virtual void copy_cell_data(const std::string& data_name, Vector<int>& mesh_data) override;
+
+    virtual Array<double> get_point_data(const std::string& data_name) override;
+    virtual std::vector<std::string> get_point_data_names() override;
+
+    virtual bool has_cell_data(const std::string& data_name) override;
+    virtual bool has_point_data(const std::string& data_name) override;
+
     virtual void set_connectivity(const int nsd, const Array<int>& conn, const int pid = 0) override;
 
     virtual void set_element_data(const std::string& data_name, const Array<double>& data) override;
@@ -101,15 +117,22 @@ class VtkVtuData : public VtkData {
     virtual int num_points() const override;
     virtual void read_file(const std::string& file_name) override;
 
-    void copy_points(Array<double>& points) override;
-    void copy_point_data(const std::string& data_name, Array<double>& mesh_data) override;
-    void copy_point_data(const std::string& data_name, Vector<double>& mesh_data) override;
-    void copy_point_data(const std::string& data_name, Vector<int>& mesh_data);
+    virtual void copy_points(Array<double>& points) override;
+
+    virtual void copy_point_data(const std::string& data_name, Array<double>& mesh_data) override;
+    virtual void copy_point_data(const std::string& data_name, Vector<double>& mesh_data) override;
+    virtual void copy_point_data(const std::string& data_name, Vector<int>& mesh_data) override;
+
+    virtual void copy_cell_data(const std::string& data_name, Vector<int>& mesh_data) override;
+
+    virtual Array<double> get_point_data(const std::string& data_name) override;
+    virtual std::vector<std::string> get_point_data_names() override;
 
-    Array<double> get_point_data(const std::string& data_name);
-    std::vector<std::string> get_point_data_names();
     virtual Array<double> get_points() const override;
-    bool has_point_data(const std::string& data_name) override;
+
+    virtual bool has_cell_data(const std::string& data_name) override;
+    virtual bool has_point_data(const std::string& data_name) override;
+
     virtual void set_connectivity(const int nsd, const Array<int>& conn, const int pid = 0) override;
 
     virtual void set_element_data(const std::string& data_name, const Array<double>& data) override;

Code/Source/solver/read_msh.cpp

@@ -2086,16 +2086,32 @@ void set_dmn_id_ff(Simulation* simulation, mshType& lM, const std::string& file_
   }
 }
 
-/// @brief Read mesh domains from a vtu/vtp file.
+/// @brief Read mesh domains from a vtk vtu/vtp file.
 ///
-/// \todo [NOTE] Not implemented.
-//
-void set_dmn_id_vtk(Simulation* simulation, mshType& mesh, const std::string& file_name, const std::string& kwrd)
+void set_dmn_id_vtk(Simulation* simulation, mshType& lM, const std::string& file_name, const std::string& data_name)
 {
-  int btSiz = std::numeric_limits<int>::digits;
+  #define n_debug_set_dmn_id_vtk
+  #ifdef debug_set_dmn_id_vtk
+  DebugMsg dmsg(__func__, simulation->com_mod.cm.idcm());
+  dmsg.banner();
+  dmsg << "file_name: " << file_name; 
+  dmsg << "lM.gnEl: " << lM.gnEl;
+  #endif
+
+  if (lM.eId.size() == 0) {
+    lM.eId.resize(lM.gnEl);
+  }
+
+  Vector<int> element_data(lM.gnEl);
+  vtk_xml::read_element_data(lM, file_name, data_name, element_data);
+
+  for (int a = 0; a < lM.gnEl; a++) {
+    lM.eId(a) = lM.eId(a) | (1UL << element_data(a));
+  }
+
 }
 
-/// @brief This routines associates two faces with each other and sets gN.
+/// @brief Associate two faces with each other and set gN.
 ///
 /// Data set
 /// \code {.cpp}
@@ -2232,4 +2248,4 @@ void set_projector(Simulation* simulation, utils::stackType& avNds)
   }
 }
 
-};
+};

Code/Source/solver/read_msh.h

@@ -52,7 +52,7 @@ namespace read_msh_ns {
   void read_msh(Simulation* simulation);
 
   void set_dmn_id_ff(Simulation* simulation, mshType& mesh, const std::string& file_name);
-  void set_dmn_id_vtk(Simulation* simulation, mshType& mesh, const std::string& file_name, const std::string& kwrd);
+  void set_dmn_id_vtk(Simulation* simulation, mshType& lM, const std::string& file_name, const std::string& kwrd);
   void set_projector(Simulation* simulation, utils::stackType& avNds);
   void set_ris_projector(Simulation* simulation);
   void set_uris_meshes(Simulation* simulation);

Code/Source/solver/vtk_xml.cpp

@@ -454,6 +454,39 @@ void read_vtp(const std::string& file_name, faceType& face)
   }
 }
 
+/// @brief Read in element data from an integer vector array data 
+/// named `data_name` stored in VTK VTU file into 'element_data'. 
+///
+void read_element_data(const mshType& mesh, const std::string& file_name, const std::string& data_name, Vector<int>& element_data)
+{
+  if (FILE *file = fopen(file_name.c_str(), "r")) {
+    fclose(file);
+  } else {
+    throw std::runtime_error("The VTK element data file '" + file_name + "' can't be read.");
+  }
+
+  auto vtk_data = VtkData::create_reader(file_name);
+  if (vtk_data == nullptr) {
+    throw std::runtime_error("Failed to create VTK reader for file: " + file_name);
+  }
+
+  int num_elems = vtk_data->num_elems();
+  int num_points = vtk_data->num_points();
+
+  if (num_elems != mesh.gnEl) {
+    throw std::runtime_error("The number of elements (" + std::to_string(num_elems) +
+        ") in the file '" + file_name + "' is not equal to the number of elements ("
+        + std::to_string(mesh.gnEl) + ") for the mesh named '" + mesh.name + "'.");
+  }
+
+  if (!vtk_data->has_cell_data(data_name)) {
+    throw std::runtime_error("No CellData DataArray named '" + data_name +
+        "' found in the VTK file '" + file_name + "' for the mesh named '" + mesh.name + "'.");
+  }
+
+  vtk_data->copy_cell_data(data_name, element_data);
+}
+
 //----------------
 // read_vtp_pdata
 //----------------

Code/Source/solver/vtk_xml.h

@@ -22,6 +22,8 @@ void read_vtu(const std::string& file_name, mshType& mesh);
 
 void read_precomputed_solution_vtu(const std::string& file_name, const std::string& field_name, mshType& mesh);
 
+void read_element_data(const mshType& mesh, const std::string& fName, const std::string& kwrd, Vector<int>& tmpR);
+
 void read_vtu_pdata(const std::string& fName, const std::string& kwrd, const int nsd, const int m, const int idx, mshType& mesh);
 
 void read_vtus(Simulation* simulation, Array<double>& lA, Array<double>& lY, Array<double>& lD, const std::string& fName);