Merge pull request #96 from robertodr/restart-cavity

Fix up read-in of cavity for restart

Author: bast

CHANGELOG.md

@@ -25,6 +25,15 @@
   file, instead of overwriting the user provided file. The temporary file is
   removed after it has been parsed. Fixes #91 as noted by @ilfreddy.
 
+### Fixed
+
+- A bug in the initialization of a restart cavity from old `.npz` files.
+  Currently, the `.npz` file saves sphere center, arcs and vertices of each
+  finite element. This information is in fact needed to plot the cavity using
+  the Python script in `tools`. Older `.npz` files did not contain this
+  information and could not be read in. The additional information is read in
+  as arrays of zeros in case it is not present on file.
+
 ## [Version 1.1.10] - 2017-03-27
 
 ### Changed

src/cavity/ICavity.cpp

@@ -29,10 +29,10 @@
 
 #include <Eigen/Core>
 
-#include "utils/cnpy.hpp"
 #include "utils/MathUtils.hpp"
 #include "utils/Sphere.hpp"
 #include "utils/Symmetry.hpp"
+#include "utils/cnpy.hpp"
 
 namespace pcm {
 ICavity::ICavity() : nElements_(0), built(false) {}
@@ -102,10 +102,15 @@ void ICavity::loadCavity(const std::string & fname) {
     PCMSOLVER_ERROR("elementArea_: incoherent dimensions read in");
   elementArea_ = cnpy::custom::npy_to_eigen<double>(raw_weights);
   // 2. Get the element sphere center
-  cnpy::NpyArray raw_elSphCenter = loaded_cavity["elSphCenter"];
-  if (raw_elSphCenter.shape[1] != nElements_)
-    PCMSOLVER_ERROR("elementSphereCenter_: incoherent dimensions read in");
-  elementSphereCenter_ = cnpy::custom::npy_to_eigen<double>(raw_elSphCenter);
+  if (loaded_cavity.find("elSphCenter") == loaded_cavity.end()) {
+    // Element sphere center was not found on file, fill it up with zeros
+    elementSphereCenter_ = Eigen::Matrix3Xd::Zero(3, nElements_);
+  } else {
+    cnpy::NpyArray raw_elSphCenter = loaded_cavity["elSphCenter"];
+    if (raw_elSphCenter.shape[1] != nElements_)
+      PCMSOLVER_ERROR("elementSphereCenter_: incoherent dimensions read in");
+    elementSphereCenter_ = cnpy::custom::npy_to_eigen<double>(raw_elSphCenter);
+  }
   // 3. Get the element radius
   cnpy::NpyArray raw_elRadius = loaded_cavity["elRadius"];
   if (raw_elRadius.shape[0] != nElements_)
@@ -122,6 +127,9 @@ void ICavity::loadCavity(const std::string & fname) {
     PCMSOLVER_ERROR("elementNormal_: incoherent dimensions read in");
   elementNormal_ = cnpy::custom::npy_to_eigen<double>(raw_normals);
 
+  bool has_arcs = loaded_cavity.find("arcs_0") == loaded_cavity.end() ? false : true;
+  bool has_vertices =
+      loaded_cavity.find("vertices_0") == loaded_cavity.end() ? false : true;
   // Reconstruct the elements_ vector
   for (PCMSolverIndex i = 0; i < nElements_; ++i) {
     bool irr = false;
@@ -131,10 +139,20 @@ void ICavity::loadCavity(const std::string & fname) {
     Sphere sph(elementSphereCenter_.col(i), elementRadius_(i));
     Eigen::Matrix3Xd vertices, arcs;
     // 6. Get vertices and arcs
-    cnpy::NpyArray raw_vertices = loaded_cavity["vertices_" + pcm::to_string(i)];
-    vertices = cnpy::custom::npy_to_eigen<double>(raw_vertices);
-    cnpy::NpyArray raw_arcs = loaded_cavity["arcs_" + pcm::to_string(i)];
-    arcs = cnpy::custom::npy_to_eigen<double>(raw_arcs);
+    if (has_vertices) {
+      cnpy::NpyArray raw_vertices = loaded_cavity["vertices_" + pcm::to_string(i)];
+      vertices = cnpy::custom::npy_to_eigen<double>(raw_vertices);
+    } else {
+      // Vertices were not found on file, fill them up with zeros
+      vertices = Eigen::Matrix3Xd::Zero(3, 3);
+    }
+    if (has_arcs) {
+      cnpy::NpyArray raw_arcs = loaded_cavity["arcs_" + pcm::to_string(i)];
+      arcs = cnpy::custom::npy_to_eigen<double>(raw_arcs);
+    } else {
+      // Arcs were not found on file, fill them up with zeros
+      arcs = Eigen::Matrix3Xd::Zero(3, 3);
+    }
     if (arcs.cols() != vertices.cols())
       PCMSOLVER_ERROR("Inconsistent number of vertices read from file for element " +
                       pcm::to_string(i));

src/cavity/RestartCavity.cpp

@@ -35,7 +35,7 @@ namespace pcm {
 namespace cavity {
 std::ostream & RestartCavity::printCavity(std::ostream & os) {
   os << "Cavity type: Restart" << std::endl;
-  os << "Number of finite elements = " << nElements_;
+  os << "Number of finite elements = " << nElements_ << std::endl;
   return os;
 }