| 56 |
|
const std::string& filename, |
| 57 |
|
const std::string& sele1, |
| 58 |
|
const std::string& sele2, |
| 59 |
< |
RealType rCut, RealType voxelSize, |
| 59 |
> |
RealType rCut, |
| 60 |
> |
RealType voxelSize, |
| 61 |
|
RealType gaussWidth) |
| 62 |
< |
: StaticAnalyser(info, filename), selectionScript1_(sele1), |
| 63 |
< |
evaluator1_(info), seleMan1_(info), selectionScript2_(sele2), |
| 64 |
< |
evaluator2_(info), seleMan2_(info), rCut_(rCut), voxelSize_(voxelSize), |
| 65 |
< |
gaussWidth_(gaussWidth) { |
| 62 |
> |
: StaticAnalyser(info, filename), |
| 63 |
> |
selectionScript1_(sele1), selectionScript2_(sele2), |
| 64 |
> |
seleMan1_(info), seleMan2_(info), evaluator1_(info), evaluator2_(info), |
| 65 |
> |
rCut_(rCut), voxelSize_(voxelSize), gaussWidth_(gaussWidth) { |
| 66 |
|
|
| 67 |
|
evaluator1_.loadScriptString(sele1); |
| 68 |
|
if (!evaluator1_.isDynamic()) { |
| 115 |
|
RealType cospsi; |
| 116 |
|
RealType Qk; |
| 117 |
|
std::vector<std::pair<RealType,StuntDouble*> > myNeighbors; |
| 118 |
< |
std::vector<std::pair<Vector3d, RealType> > qvals; |
| 119 |
< |
std::vector<std::pair<Vector3d, RealType> >::iterator qiter; |
| 118 |
> |
//std::vector<std::pair<Vector3d, RealType> > qvals; |
| 119 |
> |
//std::vector<std::pair<Vector3d, RealType> >::iterator qiter; |
| 120 |
|
int isd1; |
| 121 |
|
int isd2; |
| 122 |
|
|
| 123 |
+ |
|
| 124 |
+ |
int kMax = int(5.0 * gaussWidth_ / voxelSize_); |
| 125 |
+ |
int kSqLim = kMax*kMax; |
| 126 |
+ |
cerr << "gw = " << gaussWidth_ << " vS = " << voxelSize_ << " kMax = " |
| 127 |
+ |
<< kMax << " kSqLim = " << kSqLim << "\n"; |
| 128 |
+ |
|
| 129 |
|
DumpReader reader(info_, dumpFilename_); |
| 130 |
|
int nFrames = reader.getNFrames(); |
| 131 |
|
|
| 132 |
|
for (int istep = 0; istep < nFrames; istep += step_) { |
| 133 |
|
reader.readFrame(istep); |
| 134 |
|
|
| 135 |
+ |
currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot(); |
| 136 |
+ |
Mat3x3d hmat = currentSnapshot_->getHmat(); |
| 137 |
+ |
Vector3d halfBox = Vector3d(hmat(0,0), hmat(1,1), hmat(2,2)) / 2.0; |
| 138 |
+ |
|
| 139 |
|
if (evaluator1_.isDynamic()) { |
| 140 |
|
seleMan1_.setSelectionSet(evaluator1_.evaluate()); |
| 141 |
|
} |
| 153 |
|
} |
| 154 |
|
} |
| 155 |
|
|
| 156 |
< |
qvals.clear(); |
| 156 |
> |
//qvals.clear(); |
| 157 |
|
|
| 158 |
|
// outer loop is over the selected StuntDoubles: |
| 159 |
|
for (sd = seleMan1_.beginSelected(isd1); sd != NULL; |
| 224 |
|
if (nang > 0) { |
| 225 |
|
if (usePeriodicBoundaryConditions_) |
| 226 |
|
currentSnapshot_->wrapVector(rk); |
| 227 |
< |
qvals.push_back(std::make_pair(rk, Qk)); |
| 228 |
< |
} |
| 229 |
< |
} |
| 227 |
> |
//qvals.push_back(std::make_pair(rk, Qk)); |
| 228 |
> |
|
| 229 |
> |
Vector3d pos = rk + halfBox; |
| 230 |
|
|
| 231 |
< |
for (int i = 0; i < nBins_(0); ++i) { |
| 232 |
< |
for(int j = 0; j < nBins_(1); ++j) { |
| 233 |
< |
for(int k = 0; k < nBins_(2); ++k) { |
| 234 |
< |
Vector3d pos = Vector3d(i, j, k) * voxelSize_; |
| 235 |
< |
for(qiter = qvals.begin(); qiter != qvals.end(); ++qiter) { |
| 236 |
< |
Vector3d d = pos - (*qiter).first; |
| 237 |
< |
RealType denom = pow(2.0 * sqrt(M_PI) * gaussWidth_, 3); |
| 238 |
< |
RealType exponent = -dot(d,d) / pow(2.0*gaussWidth_, 2); |
| 239 |
< |
RealType weight = exp(exponent) / denom; |
| 240 |
< |
count_[i][j][k] += weight; |
| 241 |
< |
hist_[i][j][k] += weight * (*qiter).second; |
| 231 |
> |
|
| 232 |
> |
Vector3i whichVoxel(int(pos[0] / voxelSize_), |
| 233 |
> |
int(pos[1] / voxelSize_), |
| 234 |
> |
int(pos[2] / voxelSize_)); |
| 235 |
> |
|
| 236 |
> |
for (int l = -kMax; l <= kMax; l++) { |
| 237 |
> |
for (int m = -kMax; m <= kMax; m++) { |
| 238 |
> |
for (int n = -kMax; n <= kMax; n++) { |
| 239 |
> |
int kk = l*l + m*m + n*n; |
| 240 |
> |
if(kk <= kSqLim) { |
| 241 |
> |
|
| 242 |
> |
int ll = (whichVoxel[0] + l) % nBins_(0); |
| 243 |
> |
ll = ll < 0 ? nBins_(0) + ll : ll; |
| 244 |
> |
int mm = (whichVoxel[1] + m) % nBins_(1); |
| 245 |
> |
mm = mm < 0 ? nBins_(1) + mm : mm; |
| 246 |
> |
int nn = (whichVoxel[2] + n) % nBins_(2); |
| 247 |
> |
nn = nn < 0 ? nBins_(2) + nn : nn; |
| 248 |
> |
|
| 249 |
> |
Vector3d bPos = Vector3d(ll,mm,nn) * voxelSize_ - halfBox; |
| 250 |
> |
Vector3d d = bPos - rk; |
| 251 |
> |
currentSnapshot_->wrapVector(d); |
| 252 |
> |
RealType denom = pow(2.0 * sqrt(M_PI) * gaussWidth_, 3); |
| 253 |
> |
RealType exponent = -dot(d,d) / pow(2.0*gaussWidth_, 2); |
| 254 |
> |
RealType weight = exp(exponent) / denom; |
| 255 |
> |
count_[ll][mm][nn] += weight; |
| 256 |
> |
hist_[ll][mm][nn] += weight * Qk; |
| 257 |
> |
} |
| 258 |
> |
} |
| 259 |
|
} |
| 260 |
|
} |
| 261 |
|
} |
| 262 |
|
} |
| 263 |
+ |
|
| 264 |
+ |
// for (int i = 0; i < nBins_(0); ++i) { |
| 265 |
+ |
// for(int j = 0; j < nBins_(1); ++j) { |
| 266 |
+ |
// for(int k = 0; k < nBins_(2); ++k) { |
| 267 |
+ |
// Vector3d pos = Vector3d(i, j, k) * voxelSize_ - halfBox; |
| 268 |
+ |
// for(qiter = qvals.begin(); qiter != qvals.end(); ++qiter) { |
| 269 |
+ |
// Vector3d d = pos - (*qiter).first; |
| 270 |
+ |
// currentSnapshot_->wrapVector(d); |
| 271 |
+ |
// RealType denom = pow(2.0 * sqrt(M_PI) * gaussWidth_, 3); |
| 272 |
+ |
// RealType exponent = -dot(d,d) / pow(2.0*gaussWidth_, 2); |
| 273 |
+ |
// RealType weight = exp(exponent) / denom; |
| 274 |
+ |
// count_[i][j][k] += weight; |
| 275 |
+ |
// hist_[i][j][k] += weight * (*qiter).second; |
| 276 |
+ |
// } |
| 277 |
+ |
// } |
| 278 |
+ |
// } |
| 279 |
+ |
// } |
| 280 |
|
} |
| 281 |
|
writeQxyz(); |
| 282 |
|
} |
| 283 |
|
|
| 284 |
|
void TetrahedralityParamXYZ::writeQxyz() { |
| 285 |
+ |
|
| 286 |
+ |
Mat3x3d hmat = info_->getSnapshotManager()->getCurrentSnapshot()->getHmat(); |
| 287 |
+ |
|
| 288 |
|
// normalize by total weight in voxel: |
| 289 |
|
for (unsigned int i = 0; i < hist_.size(); ++i) { |
| 290 |
|
for(unsigned int j = 0; j < hist_[i].size(); ++j) { |
| 291 |
|
for(unsigned int k = 0;k < hist_[i][j].size(); ++k) { |
| 292 |
< |
hist_[i][j][k] /= count_[i][j][k]; |
| 292 |
> |
hist_[i][j][k] = hist_[i][j][k] / count_[i][j][k]; |
| 293 |
|
} |
| 294 |
|
} |
| 295 |
|
} |
| 296 |
|
|
| 297 |
|
std::ofstream qXYZstream(outputFilename_.c_str()); |
| 298 |
|
if (qXYZstream.is_open()) { |
| 299 |
+ |
qXYZstream << "# AmiraMesh ASCII 1.0\n\n"; |
| 300 |
+ |
qXYZstream << "# Dimensions in x-, y-, and z-direction\n"; |
| 301 |
+ |
qXYZstream << " define Lattice " << hist_.size() << " " |
| 302 |
+ |
<< hist_[0].size() << " " << hist_[0][0].size() << "\n"; |
| 303 |
+ |
|
| 304 |
+ |
qXYZstream << "Parameters {\n"; |
| 305 |
+ |
qXYZstream << " CoordType \"uniform\",\n"; |
| 306 |
+ |
qXYZstream << " # BoundingBox is xmin xmax ymin ymax zmin zmax\n"; |
| 307 |
+ |
qXYZstream << " BoundingBox 0.0 " << hmat(0,0) << |
| 308 |
+ |
" 0.0 " << hmat(1,1) << |
| 309 |
+ |
" 0.0 " << hmat(2,2) << "\n"; |
| 310 |
+ |
qXYZstream << "}\n"; |
| 311 |
+ |
|
| 312 |
+ |
qXYZstream << "Lattice { double ScalarField } = @1\n"; |
| 313 |
+ |
|
| 314 |
+ |
qXYZstream << "@1\n"; |
| 315 |
|
|
| 316 |
< |
for (unsigned int i = 0; i < hist_.size(); ++i) { |
| 317 |
< |
for(unsigned int j = 0; j < hist_[i].size(); ++j) { |
| 318 |
< |
for(unsigned int k = 0;k < hist_[i][j].size(); ++k) { |
| 319 |
< |
qXYZstream.write(reinterpret_cast<char *>( &hist_[i][j][k] ), |
| 320 |
< |
sizeof( hist_[i][j][k] )); |
| 316 |
> |
for (std::size_t k = 0; k < hist_[0][0].size(); ++k) { |
| 317 |
> |
for(std::size_t j = 0; j < hist_[0].size(); ++j) { |
| 318 |
> |
for(std::size_t i = 0; i < hist_.size(); ++i) { |
| 319 |
> |
qXYZstream << hist_[i][j][k] << " "; |
| 320 |
> |
|
| 321 |
> |
//qXYZstream.write(reinterpret_cast<char *>( &hist_[i][j][k] ), |
| 322 |
> |
// sizeof( hist_[i][j][k] )); |
| 323 |
|
} |
| 324 |
|
} |
| 325 |
|
} |
