Skip to content

Commit 520ad9e

Browse files
author
Roberto Di Remigio
committed
Add some custom overloads to cnpy functions
1 parent 6923a5d commit 520ad9e

25 files changed

+323
-400
lines changed

include/Config.hpp.in

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,8 @@
4141
#define attribute(x)
4242
#endif
4343

44+
#include <boost/current_function.hpp>
45+
4446
#include "Cxx11Workarounds.hpp"
4547
#include "ErrorHandling.hpp"
4648
#include "LoggerInterface.hpp"

include/ErrorHandling.hpp

Lines changed: 15 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -68,35 +68,30 @@
6868
* See also here: http://www.boost.org/doc/libs/1_59_0/doc/html/boost_staticassert.html
6969
*/
7070

71-
/*! \brief Kills execution and prints out error message */
72-
inline void pcmsolver_die(const std::string & message)
73-
{
74-
pcmsolver_die(message.c_str());
75-
}
76-
77-
/*! \brief Kills execution and prints out error code and error message */
78-
inline void pcmsolver_die(const std::string & message, int code)
79-
{
80-
pcmsolver_die(message.c_str(), code);
81-
}
82-
83-
84-
/*! \brief Kills execution and prints out error message */
85-
inline void pcmsolver_die(const char * message)
71+
/*! \brief Kills execution and prints out error message to stderr
72+
* \param message Error message
73+
* \param function Name of the function killing execution
74+
* \param code Error code. Defaults to EXIT_FAILURE
75+
*/
76+
inline void pcmsolver_die(const std::string & message, const std::string & function, int code = EXIT_FAILURE)
8677
{
87-
std::fprintf(stderr, "PCMSolver fatal error: %s\n", message);
88-
std::exit(EXIT_FAILURE);
78+
pcmsolver_die(message.c_str(), function.c_str(), code);
8979
}
9080

91-
/*! \brief Kills execution and prints out error code and error message */
92-
inline void pcmsolver_die(const char * message, int code)
81+
/*! \brief Kills execution and prints out error message to stderr
82+
* \param message Error message
83+
* \param function Name of the function killing execution
84+
* \param code Error code. Defaults to EXIT_FAILURE
85+
*/
86+
inline void pcmsolver_die(const char * message, const char * function, int code = EXIT_FAILURE)
9387
{
88+
std::fprintf(stderr, "In function: %s\n", function);
9489
std::fprintf(stderr, "PCMSolver fatal error %i: %s\n", code, message);
9590
std::exit(EXIT_FAILURE);
9691
}
9792

9893
/// Macro to be used to signal error conditions
99-
#define PCMSOLVER_ERROR(arg) pcmsolver_die(arg)
94+
#define PCMSOLVER_ERROR(arg, func) pcmsolver_die(arg, func)
10095

10196
/// Macro to be used for assertions
10297
#define PCMSOLVER_ASSERT(arg) assert(arg)

src/bi_operators/CollocationIntegrator.hpp

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -112,25 +112,25 @@ struct CollocationIntegrator
112112
/**@{ Single and double layer potentials for a IonicLiquid Green's function by collocation */
113113
template <typename DerivativeTraits>
114114
Eigen::MatrixXd singleLayer(const IonicLiquid<DerivativeTraits, CollocationIntegrator> & /* gf */, const std::vector<Element> & e) const {
115-
PCMSOLVER_ERROR("CollocationIntegrator::singleLayer not implemented yet for IonicLiquid");
115+
PCMSOLVER_ERROR("CollocationIntegrator::singleLayer not implemented yet for IonicLiquid", BOOST_CURRENT_FUNCTION);
116116
return Eigen::MatrixXd::Zero(e.size(), e.size());
117117
}
118118
template <typename DerivativeTraits>
119119
Eigen::MatrixXd doubleLayer(const IonicLiquid<DerivativeTraits, CollocationIntegrator> & /* gf */, const std::vector<Element> & e) const {
120-
PCMSOLVER_ERROR("CollocationIntegrator::doubleLayer not implemented yet for IonicLiquid");
120+
PCMSOLVER_ERROR("CollocationIntegrator::doubleLayer not implemented yet for IonicLiquid", BOOST_CURRENT_FUNCTION);
121121
return Eigen::MatrixXd::Zero(e.size(), e.size());
122122
}
123123
/**@}*/
124124

125125
/**@{ Single and double layer potentials for an AnisotropicLiquid Green's function by collocation */
126126
template <typename DerivativeTraits>
127127
Eigen::MatrixXd singleLayer(const AnisotropicLiquid<DerivativeTraits, CollocationIntegrator> & /* gf */, const std::vector<Element> & e) const {
128-
PCMSOLVER_ERROR("CollocationIntegrator::singleLayer not implemented yet for AnisotropicLiquid");
128+
PCMSOLVER_ERROR("CollocationIntegrator::singleLayer not implemented yet for AnisotropicLiquid", BOOST_CURRENT_FUNCTION);
129129
return Eigen::MatrixXd::Zero(e.size(), e.size());
130130
}
131131
template <typename DerivativeTraits>
132132
Eigen::MatrixXd doubleLayer(const AnisotropicLiquid<DerivativeTraits, CollocationIntegrator> & /* gf */, const std::vector<Element> & e) const {
133-
PCMSOLVER_ERROR("CollocationIntegrator::doubleLayer not implemented yet for AnisotropicLiquid");
133+
PCMSOLVER_ERROR("CollocationIntegrator::doubleLayer not implemented yet for AnisotropicLiquid", BOOST_CURRENT_FUNCTION);
134134
return Eigen::MatrixXd::Zero(e.size(), e.size());
135135
}
136136
/**@}*/

src/bi_operators/PurisimaIntegrator.hpp

Lines changed: 12 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -121,29 +121,29 @@ struct PurisimaIntegrator
121121

122122
/**@{ Single and double layer potentials for a IonicLiquid Green's function by collocation */
123123
template <typename DerivativeTraits>
124-
Eigen::MatrixXd singleLayer(const IonicLiquid<DerivativeTraits, PurisimaIntegrator> & /* gf */, const std::vector<Element> & e) const {
125-
PCMSOLVER_ERROR("PurisimaIntegrator::singleLayer not implemented yet for IonicLiquid");
124+
Eigen::MatrixXd singleLayer(const IonicLiquid<DerivativeTraits, PurisimaIntegrator> & /* gf */, const std::vector<Element> & e) const {
125+
PCMSOLVER_ERROR("PurisimaIntegrator::singleLayer not implemented yet for IonicLiquid", BOOST_CURRENT_FUNCTION);
126126
return Eigen::MatrixXd::Zero(e.size(), e.size());
127-
}
127+
}
128128
template <typename DerivativeTraits>
129-
Eigen::MatrixXd doubleLayer(const IonicLiquid<DerivativeTraits, PurisimaIntegrator> & /* gf */, const std::vector<Element> & e) const {
130-
PCMSOLVER_ERROR("PurisimaIntegrator::doubleLayer not implemented yet for IonicLiquid");
129+
Eigen::MatrixXd doubleLayer(const IonicLiquid<DerivativeTraits, PurisimaIntegrator> & /* gf */, const std::vector<Element> & e) const {
130+
PCMSOLVER_ERROR("PurisimaIntegrator::doubleLayer not implemented yet for IonicLiquid", BOOST_CURRENT_FUNCTION);
131131
return Eigen::MatrixXd::Zero(e.size(), e.size());
132-
}
132+
}
133133
/**@}*/
134134

135135
/**@{ Single and double layer potentials for an AnisotropicLiquid Green's function by collocation */
136136
template <typename DerivativeTraits>
137-
Eigen::MatrixXd singleLayer(const AnisotropicLiquid<DerivativeTraits, PurisimaIntegrator> & /* gf */, const std::vector<Element> & e) const {
138-
PCMSOLVER_ERROR("PurisimaIntegrator::singleLayer not implemented yet for AnisotropicLiquid");
137+
Eigen::MatrixXd singleLayer(const AnisotropicLiquid<DerivativeTraits, PurisimaIntegrator> & /* gf */, const std::vector<Element> & e) const {
138+
PCMSOLVER_ERROR("PurisimaIntegrator::singleLayer not implemented yet for AnisotropicLiquid", BOOST_CURRENT_FUNCTION);
139139
return Eigen::MatrixXd::Zero(e.size(), e.size());
140140

141-
}
141+
}
142142
template <typename DerivativeTraits>
143-
Eigen::MatrixXd doubleLayer(const AnisotropicLiquid<DerivativeTraits, PurisimaIntegrator> & /* gf */, const std::vector<Element> & e) const {
144-
PCMSOLVER_ERROR("PurisimaIntegrator::doubleLayer not implemented yet for AnisotropicLiquid");
143+
Eigen::MatrixXd doubleLayer(const AnisotropicLiquid<DerivativeTraits, PurisimaIntegrator> & /* gf */, const std::vector<Element> & e) const {
144+
PCMSOLVER_ERROR("PurisimaIntegrator::doubleLayer not implemented yet for AnisotropicLiquid", BOOST_CURRENT_FUNCTION);
145145
return Eigen::MatrixXd::Zero(e.size(), e.size());
146-
}
146+
}
147147
/**@}*/
148148

149149
/**@{ Single and double layer potentials for a SphericalDiffuse Green's function by collocation */

src/bin/run_pcm.cpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -65,7 +65,7 @@ int main(int argc, char * argv[])
6565
std::ofstream out_stream;
6666
out_stream.open("pcmsolver.out");
6767

68-
if (argc > 2) PCMSOLVER_ERROR("Too many arguments supplied to run_pcm");
68+
if (argc > 2) PCMSOLVER_ERROR("Too many arguments supplied", "run_pcm");
6969
Input parsed(argv[1]);
7070
parsed.initMolecule();
7171

src/cavity/Cavity.cpp

Lines changed: 81 additions & 112 deletions
Original file line numberDiff line numberDiff line change
@@ -40,123 +40,92 @@
4040

4141
void Cavity::saveCavity(const std::string & fname)
4242
{
43-
/*
44-
std::ofstream weights("weights.txt", std::ios_base::out);
45-
// First line in the weights file is the number of elements.
46-
// This is for a sanity-check on the save/load operations.
47-
weights << std::setprecision(std::numeric_limits<double>::digits10) << nElements_ << std::endl;
48-
weights << std::setprecision(std::numeric_limits<double>::digits10) << elementArea_ << std::endl;
49-
std::ofstream elRadius("element_radius.txt", std::ios_base::out);
50-
elRadius << std::setprecision(std::numeric_limits<double>::digits10) << elementRadius_ << std::endl;
51-
std::ofstream centers("centers.txt", std::ios_base::out);
52-
centers << std::setprecision(std::numeric_limits<double>::digits10) << elementCenter_ << std::endl;
53-
std::ofstream normals("normals.txt", std::ios_base::out);
54-
normals << std::setprecision(std::numeric_limits<double>::digits10) << elementNormal_ << std::endl;
55-
*/
43+
/*
44+
std::ofstream weights("weights.txt", std::ios_base::out);
45+
// First line in the weights file is the number of elements.
46+
// This is for a sanity-check on the save/load operations.
47+
weights << std::setprecision(std::numeric_limits<double>::digits10) << nElements_ << std::endl;
48+
weights << std::setprecision(std::numeric_limits<double>::digits10) << elementArea_ << std::endl;
49+
std::ofstream elRadius("element_radius.txt", std::ios_base::out);
50+
elRadius << std::setprecision(std::numeric_limits<double>::digits10) << elementRadius_ << std::endl;
51+
std::ofstream centers("centers.txt", std::ios_base::out);
52+
centers << std::setprecision(std::numeric_limits<double>::digits10) << elementCenter_ << std::endl;
53+
std::ofstream normals("normals.txt", std::ios_base::out);
54+
normals << std::setprecision(std::numeric_limits<double>::digits10) << elementNormal_ << std::endl;
55+
*/
5656

57-
// Write everything in a single .npz binary file
58-
unsigned int dim = static_cast<unsigned int>(nElements_);
59-
// Write the number of elements, it will be used to check sanity of the save/load operations.
60-
const unsigned int shape[] = {1};
61-
cnpy::npz_save(fname, "elements", &dim, shape, 1, "w", false);
62-
// Write weights
63-
const unsigned int weights_shape[] = {dim};
64-
cnpy::npz_save(fname, "weights", elementArea_.data(), weights_shape, 1, "a", true);
65-
// Write element sphere center
66-
const unsigned int elSphCenter_shape[] = {3, dim};
67-
cnpy::npz_save(fname, "elSphCenter", elementSphereCenter_.data(), elSphCenter_shape, 2, "a",
68-
true);
69-
// Write element radius
70-
const unsigned int elRadius_shape[] = {dim};
71-
cnpy::npz_save(fname, "elRadius", elementRadius_.data(), elRadius_shape, 1, "a",
72-
true);
73-
// Write centers
74-
const unsigned int centers_shape[] = {3, dim};
75-
cnpy::npz_save(fname, "centers", elementCenter_.data(), centers_shape, 2, "a", true);
76-
// Write normals
77-
const unsigned int normals_shape[] = {3, dim};
78-
cnpy::npz_save(fname, "normals", elementNormal_.data(), normals_shape, 2, "a", true);
79-
// Write vertices TODO!!!!
80-
//const unsigned int vertices_shape[] = {dim, 10, 3};
81-
// Write arcs TODO!!!!
82-
//const unsigned int arcs_shape[] = {dim, 10, 3};
57+
// Write everything in a single .npz binary file
58+
unsigned int dim = static_cast<unsigned int>(nElements_);
59+
// Write the number of elements, it will be used to check sanity of the save/load operations.
60+
const unsigned int shape[] = {1};
61+
cnpy::npz_save(fname, "elements", &dim, shape, 1, "w", false);
62+
// Write weights
63+
cnpy::custom::npz_save(fname, "weights", elementArea_);
64+
// Write element sphere center
65+
cnpy::custom::npz_save(fname, "elSphCenter", elementSphereCenter_);
66+
// Write element radius
67+
cnpy::custom::npz_save(fname, "elRadius", elementRadius_);
68+
// Write centers
69+
cnpy::custom::npz_save(fname, "centers", elementCenter_);
70+
// Write normals
71+
cnpy::custom::npz_save(fname, "normals", elementNormal_);
72+
// Write vertices TODO!!!!
73+
//const unsigned int vertices_shape[] = {dim, 10, 3};
74+
// Write arcs TODO!!!!
75+
//const unsigned int arcs_shape[] = {dim, 10, 3};
8376
}
8477

8578
void Cavity::loadCavity(const std::string & fname)
8679
{
87-
// We initialize molecule_ to a dummy Molecule
88-
molecule_ = Molecule();
89-
// Load the .npz binary file and then traverse it to get the data needed to rebuild the cavity.
90-
cnpy::npz_t loaded_cavity = cnpy::npz_load(fname);
91-
// 0. Get the number of elements
92-
cnpy::NpyArray raw_ele = loaded_cavity["elements"];
93-
int * ne = reinterpret_cast<int*>(raw_ele.data);
94-
nElements_ = *ne;
95-
// Set the size of the irreducible portion of the cavity
96-
// it will be the same as the total size, since a restarted cavity is always C1
97-
nIrrElements_ = nElements_;
80+
// We initialize molecule_ to a dummy Molecule
81+
molecule_ = Molecule();
82+
// Load the .npz binary file and then traverse it to get the data needed to rebuild the cavity.
83+
cnpy::npz_t loaded_cavity = cnpy::npz_load(fname);
84+
// 0. Get the number of elements
85+
cnpy::NpyArray raw_ele = loaded_cavity["elements"];
86+
int * ne = reinterpret_cast<int*>(raw_ele.data);
87+
nElements_ = *ne;
88+
// Set the size of the irreducible portion of the cavity
89+
// it will be the same as the total size, since a restarted cavity is always C1
90+
nIrrElements_ = nElements_;
9891

99-
// 1. Get the weights
100-
cnpy::NpyArray raw_weights = loaded_cavity["weights"];
101-
int dim = raw_weights.shape[0];
102-
if (dim != nElements_) {
103-
PCMSOLVER_ERROR("A problem occurred while loading the cavity. Inconsistent dimension of weights vector!");
104-
} else {
105-
elementArea_ = getFromRawBuffer<double>(dim, 1, raw_weights.data);
106-
}
92+
// 1. Get the weights
93+
cnpy::NpyArray raw_weights = loaded_cavity["weights"];
94+
if (raw_weights.shape[0] != nElements_) PCMSOLVER_ERROR("elementArea_: incoherent dimensions read in", BOOST_CURRENT_FUNCTION);
95+
elementArea_ = cnpy::custom::npy_to_eigen(raw_weights);
96+
// 2. Get the element sphere center
97+
cnpy::NpyArray raw_elSphCenter = loaded_cavity["elSphCenter"];
98+
if (raw_elSphCenter.shape[1] != nElements_) PCMSOLVER_ERROR("elementSphereCenter_: incoherent dimensions read in", BOOST_CURRENT_FUNCTION);
99+
elementSphereCenter_ = cnpy::custom::npy_to_eigen(raw_elSphCenter);
100+
// 3. Get the element radius
101+
cnpy::NpyArray raw_elRadius = loaded_cavity["elRadius"];
102+
if (raw_elRadius.shape[0] != nElements_) PCMSOLVER_ERROR("elementRadius_: incoherent dimensions read in", BOOST_CURRENT_FUNCTION);
103+
elementRadius_ = cnpy::custom::npy_to_eigen(raw_elRadius);
104+
// 4. Get the centers
105+
cnpy::NpyArray raw_centers = loaded_cavity["centers"];
106+
if (raw_centers.shape[1] != nElements_) PCMSOLVER_ERROR("elementCenter_: incoherent dimensions read in", BOOST_CURRENT_FUNCTION);
107+
elementCenter_ = cnpy::custom::npy_to_eigen(raw_centers);
108+
// 5. Get the normal vectors
109+
cnpy::NpyArray raw_normals = loaded_cavity["normals"];
110+
if (raw_normals.shape[1] != nElements_) PCMSOLVER_ERROR("elementNormal_: incoherent dimensions read in", BOOST_CURRENT_FUNCTION);
111+
elementNormal_ = cnpy::custom::npy_to_eigen(raw_normals);
107112

108-
// 2. Get the element sphere center
109-
cnpy::NpyArray raw_elSphCenter = loaded_cavity["elSphCenter"];
110-
dim = raw_elSphCenter.shape[1];
111-
if (dim != nElements_) {
112-
PCMSOLVER_ERROR("A problem occurred while loading the cavity. Inconsistent dimension of element sphere radius matrix!");
113-
} else {
114-
elementSphereCenter_ = getFromRawBuffer<double>(3, dim, raw_elSphCenter.data);
115-
}
116-
117-
// 3. Get the element radius
118-
cnpy::NpyArray raw_elRadius = loaded_cavity["elRadius"];
119-
dim = raw_elRadius.shape[0];
120-
if (dim != nElements_) {
121-
PCMSOLVER_ERROR("A problem occurred while loading the cavity. Inconsistent dimension of element radius vector!");
122-
} else {
123-
elementRadius_ = getFromRawBuffer<double>(dim, 1, raw_elRadius.data);
124-
}
125-
126-
// 4. Get the centers
127-
cnpy::NpyArray raw_centers = loaded_cavity["centers"];
128-
dim = raw_centers.shape[1];
129-
if (dim != nElements_) {
130-
PCMSOLVER_ERROR("A problem occurred while loading the cavity. Inconsistent dimension of centers matrix!");
131-
} else {
132-
elementCenter_ = getFromRawBuffer<double>(3, dim, raw_centers.data);
133-
}
134-
135-
// 5. Get the normal vectors
136-
cnpy::NpyArray raw_normals = loaded_cavity["normals"];
137-
dim = raw_normals.shape[1];
138-
if (dim != nElements_) {
139-
PCMSOLVER_ERROR("A problem occurred while loading the cavity. Inconsistent dimension of normals matrix!");
140-
} else {
141-
elementNormal_ = getFromRawBuffer<double>(3, dim, raw_normals.data);
142-
}
143-
144-
// Reconstruct the elements_ vector
145-
for (int i = 0; i < nElements_; ++i) {
146-
bool irr = false;
147-
// PEDRA puts the irreducible tesserae first
148-
if (i < nIrrElements_) irr = true;
149-
Sphere sph(elementSphereCenter_.col(i), elementRadius_(i));
150-
int nv = 3; // BOGUS!!!
151-
Eigen::Matrix3Xd vertices, arcs;
152-
vertices.resize(Eigen::NoChange, nv); // BOGUS!!!
153-
arcs.resize(Eigen::NoChange, nv); // BOGUS!!
154-
// Populate vertices and arcs
155-
elements_.push_back(Element(nv, 0,
156-
elementArea_(i),
157-
elementCenter_.col(i),
158-
elementNormal_.col(i),
159-
irr, sph,
160-
vertices, arcs));
161-
}
113+
// Reconstruct the elements_ vector
114+
for (int i = 0; i < nElements_; ++i) {
115+
bool irr = false;
116+
// PEDRA puts the irreducible tesserae first
117+
if (i < nIrrElements_) irr = true;
118+
Sphere sph(elementSphereCenter_.col(i), elementRadius_(i));
119+
int nv = 3; // BOGUS!!!
120+
Eigen::Matrix3Xd vertices, arcs;
121+
vertices.resize(Eigen::NoChange, nv); // BOGUS!!!
122+
arcs.resize(Eigen::NoChange, nv); // BOGUS!!
123+
// Populate vertices and arcs
124+
elements_.push_back(Element(nv, 0,
125+
elementArea_(i),
126+
elementCenter_.col(i),
127+
elementNormal_.col(i),
128+
irr, sph,
129+
vertices, arcs));
130+
}
162131
}

0 commit comments

Comments
 (0)