ClumpInput.h
Go to the documentation of this file.
1 //Copyright (c) 2013-2023, The MercuryDPM Developers Team. All rights reserved.
2 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
3 //
4 //Redistribution and use in source and binary forms, with or without
5 //modification, are permitted provided that the following conditions are met:
6 // * Redistributions of source code must retain the above copyright
7 // notice, this list of conditions and the following disclaimer.
8 // * Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 // * Neither the name MercuryDPM nor the
12 // names of its contributors may be used to endorse or promote products
13 // derived from this software without specific prior written permission.
14 //
15 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
19 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 
26 
27 // This module loads clump configuration produced by MClump tool
28 
29 #ifndef CLUMP_INPUT_H
30 #define CLUMP_INPUT_H
31 
32 #include <algorithm>
33 #include <dirent.h>
34 #include <sys/types.h>
35 #include <iostream>
36 #include <fstream>
37 #include <vector>
38 #include <string>
39 #include <cmath>
40 #include <sstream>
41 #include "Math/Matrix.h"
42 #include <CMakeDefinitions.h>
43 
44 // Useful containers
45 typedef std::vector<double> DoubleVector;
46 typedef std::vector<DoubleVector> Double2DVector;
47 typedef std::vector<std::string> StringVector;
48 
49 // Helper function for alphabetical sorting of clump names
50 bool CompareFunction (std::string a, std::string b) {return a < b;}
51 
52 // Helper function to generate a random double in range (0, MAX)
53 double RandomDouble(double Max)
54 {return Max*((double) rand() / (RAND_MAX));}
55 
60 struct ClumpData
61 {
62 public:
63 
64  std::string path; // Path to MClump working directory
65  StringVector clumpNames; // Array of names of Clumps that will be used
66  DoubleVector mass; // clump mass
67  Double2DVector pebblesX; // Pebbles geometry (outer index goes over Clumps, inner - over pebbles)
71 
72  Double2DVector toi; // Clump tensor of inertia (I11, I12, I13, I21..I33)
73  Double2DVector pd; // Clump principal directions v1, v2, v3
74 
75 };
76 
77 
83 {
84  // Path to MCLump tool
85  a.path = getMercuryDPMSourceDir() + "/Tools/MClump/Clumps/";
86 
87  struct dirent *entry;
88  DIR *dir = opendir(a.path.c_str());
89  while ((entry = readdir(dir)) != NULL) {
90  if (entry->d_type == DT_DIR){
91  a.clumpNames.push_back(entry->d_name);
92  }
93  }
94  closedir(dir);
95 
96  // remove "." and ".." from the list of dirs (position of those in the list is OS-sensitive, hence this code)
97 
98  std::sort(a.clumpNames.begin(), a.clumpNames.end(), CompareFunction);//sort the vector
99  a.clumpNames.erase(a.clumpNames.begin());
100  a.clumpNames.erase(a.clumpNames.begin());
101 
102  // Show the names of available Clumps
103  for (int i = 0; i<a.clumpNames.size(); i++) std::cout << a.clumpNames[i] << std::endl;
104 
105 }
106 
112 {
113  logger(INFO, "Loading clump pebbles...");
114 
115  a.pebblesX.resize(a.clumpNames.size());
116  a.pebblesY.resize(a.clumpNames.size());
117  a.pebblesZ.resize(a.clumpNames.size());
118  a.pebblesR.resize(a.clumpNames.size());
119 
120  for (int i = 0; i < a.clumpNames.size(); i++ ){
121  std::ifstream infile((a.path + a.clumpNames[i] + "/Clump/Clump.txt").c_str(), std::ios::in | std::ios::binary);
122  std::string line{};
123  while (std::getline(infile, line)) {
124  std::istringstream iss(line);
125  std::string substring{};
126  StringVector val;
127  while (std::getline(iss, substring, ',')) val.push_back(substring);
128  a.pebblesX[i].push_back(std::stof(val[0]));
129  a.pebblesY[i].push_back(std::stof(val[1]));
130  a.pebblesZ[i].push_back(std::stof(val[2]));
131  a.pebblesR[i].push_back(std::stof(val[3]));
132  }
133  infile.close();
134  }
135  logger(INFO, "\t OK\n");
136 }
137 
142 {
143  logger(INFO, "Loading clump masses...");
144 
145  a.mass.resize(a.clumpNames.size());
146 
147  for (int i = 0; i < a.clumpNames.size(); i++ ){
148  std::ifstream infile((a.path + a.clumpNames[i] + "/Inertia/Mass.txt").c_str(), std::ios::in | std::ios::binary);
149  std::string mass;
150  infile >> mass;
151  a.mass[i] = std::stof(mass);
152  infile.close();
153  }
154  logger(INFO, "\t OK\n");
155 }
156 
162 {
163  logger(INFO, "Loading clump TOI..");
164  a.toi.resize(a.clumpNames.size());
165 
166  for (int i = 0; i < a.clumpNames.size(); i++ ){
167  std::ifstream infile((a.path + a.clumpNames[i] + "/Inertia/TOI.txt").c_str(), std::ios::in | std::ios::binary);
168  std::string line{};
169  while (std::getline(infile, line)) {
170  std::istringstream iss(line);
171  std::string substring{};
172  StringVector val;
173  while (std::getline(iss, substring, ',')) val.push_back(substring);
174  for (int k = 0; k < 3; k++) a.toi[i].push_back(std::stof(val[k]));
175  }
176  infile.close();
177  }
178  logger(INFO, "\t OK\n");
179 }
180 
186 {
187  logger(INFO, "Loading clump PD..");
188  a.pd.resize(a.clumpNames.size());
189 
190  for (int i = 0; i < a.clumpNames.size(); i++ ){
191  std::ifstream infile((a.path + a.clumpNames[i] + "/Inertia/PD.txt").c_str(), std::ios::in | std::ios::binary);
192  std::string line{};
193  while (std::getline(infile, line)) {
194  std::istringstream iss(line);
195  std::string substring{};
196  StringVector val;
197  while (std::getline(iss, substring, ',')) val.push_back(substring);
198  for (int k = 0; k < 3; k++) a.pd[i].push_back(std::stof(val[k]));
199  }
200  infile.close();
201  }
202  logger(INFO, "\t OK\n");
203 }
204 
205 
210 void LoadClumps(ClumpData &data, bool VERBOSE = false)
211 {
212 
213  logger(INFO, "LOAD CLUMP DATA");
214  LoadConf(data);
215  LoadPebbles(data);
216  LoadMass(data);
217  LoadTOI(data);
218  LoadPD(data);
219  if (VERBOSE) {
220  std::cout<<"LOADED CLUMPS"<<std::endl;
221  for (int i = 0; i < data.pebblesX.size(); i++) {
222  std::cout << data.clumpNames[i] << " mass:" << data.mass[i] << std::endl;
223  std::cout << data.clumpNames[i] << " list of pebbles:" << std::endl;
224  for (int j = 0; j < data.pebblesX[i].size(); j++) {
225  std::cout << "Pebble " << j << ": (" << data.pebblesX[i][j] << "," << data.pebblesY[i][j] << ","
226  << data.pebblesZ[i][j] << ")," << data.pebblesR[i][j] << std::endl;
227  }
228 
229  std::cout << data.clumpNames[i] << " TOI:" << std::endl;
230  std::cout << data.toi[i][0] << "," << data.toi[i][1] << "," << data.toi[i][2] << std::endl;
231  std::cout << data.toi[i][3] << "," << data.toi[i][4] << "," << data.toi[i][5] << std::endl;
232  std::cout << data.toi[i][6] << "," << data.toi[i][7] << "," << data.toi[i][8] << std::endl;
233 
234  std::cout << data.clumpNames[i] << " Principal directions:" << std::endl;
235  std::cout << data.pd[i][0] << "," << data.pd[i][1] << "," << data.pd[i][2] << std::endl;
236  std::cout << data.pd[i][3] << "," << data.pd[i][4] << "," << data.pd[i][5] << std::endl;
237  std::cout << data.pd[i][6] << "," << data.pd[i][7] << "," << data.pd[i][8] << std::endl;
238  }
239  }
240 
241 }
242 
247 ClumpData RotateClump(ClumpData data, int clump_index, DoubleVector new_pd)
248 {
249  ClumpData new_data = data;
250  new_data.pd[clump_index] = new_pd;
251  return new_data;
252 }
253 
254 
255 
256 
257 
258 
263 
264  Vec3D n1, n2, n3, ref;
265 
266  // basis vector n1
267  double r1 = RandomDouble(2) - 1.0;
268  double r2 = RandomDouble(1);
269 
270  double theta = acos(r1); // Note that for isotropy of n1 theta is NOT uniform!
271  double phi = 2 * M_PI * r2;
272 
273  n1 = Vec3D(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
274 
275  // basis vector n2
276  r1 = RandomDouble(1);
277  r2 = RandomDouble(1);
278 
279  theta = acos(2 * r1 - 1); // Note that for isotropy of n1 theta is NOT uniform!
280  phi = 2 * M_PI * r2;
281 
282  ref = Vec3D(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
283  n2 = Vec3D::cross(ref, n1); n2.normalise();
284  n3 = Vec3D::cross(n1, n2); n3.normalise();
285  return DoubleVector{n1.X, n1.Y, n1.Z, n2.X, n2.Y, n2.Z, n3.X, n3.Y, n3.Z};
286 
287 }
288 
289 
290 #endif // CLUMP_INPUT_H
const std::string getMercuryDPMSourceDir()
This file is used for generating definitions that give access to CMakeVariables from within a cpp fil...
Definition: CMakeDefinitions.cc:32
void LoadConf(ClumpData &a)
Definition: ClumpInput.h:82
void LoadPD(ClumpData &a)
Definition: ClumpInput.h:185
double RandomDouble(double Max)
Definition: ClumpInput.h:53
void LoadPebbles(ClumpData &a)
Definition: ClumpInput.h:111
void LoadClumps(ClumpData &data, bool VERBOSE=false)
Definition: ClumpInput.h:210
void LoadTOI(ClumpData &a)
Definition: ClumpInput.h:161
DoubleVector UniformRandomPDs()
Definition: ClumpInput.h:262
std::vector< double > DoubleVector
Definition: ClumpInput.h:45
std::vector< std::string > StringVector
Definition: ClumpInput.h:47
void LoadMass(ClumpData &a)
Definition: ClumpInput.h:141
bool CompareFunction(std::string a, std::string b)
Definition: ClumpInput.h:50
std::vector< DoubleVector > Double2DVector
Definition: ClumpInput.h:46
ClumpData RotateClump(ClumpData data, int clump_index, DoubleVector new_pd)
Definition: ClumpInput.h:247
LL< Log::VERBOSE > VERBOSE
Verbose information.
Definition: Logger.cc:57
LL< Log::INFO > INFO
Info log level.
Definition: Logger.cc:55
Logger< MERCURYDPM_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here.
Definition: Vector.h:51
Mdouble Y
Definition: Vector.h:66
Mdouble Z
Definition: Vector.h:66
static Vec3D cross(const Vec3D &a, const Vec3D &b)
Calculates the cross product of two Vec3D: .
Definition: Vector.cc:163
Mdouble X
the vector components
Definition: Vector.h:66
void normalise()
Makes this Vec3D unit length.
Definition: Vector.cc:123
const std::complex< Mdouble > i
Definition: ExtendedMath.h:51
Mdouble cos(Mdouble x)
Definition: ExtendedMath.cc:64
Mdouble sin(Mdouble x)
Definition: ExtendedMath.cc:44
Definition: ClumpInput.h:61
Double2DVector toi
Definition: ClumpInput.h:72
Double2DVector pd
Definition: ClumpInput.h:73
StringVector clumpNames
Definition: ClumpInput.h:65
Double2DVector pebblesX
Definition: ClumpInput.h:67
DoubleVector mass
Definition: ClumpInput.h:66
Double2DVector pebblesZ
Definition: ClumpInput.h:69
Double2DVector pebblesY
Definition: ClumpInput.h:68
Double2DVector pebblesR
Definition: ClumpInput.h:70
std::string path
Definition: ClumpInput.h:64