~ubuntu-branches/ubuntu/vivid/travis/vivid-proposed

Committer: Package Import Robot
Author(s): Daniel Leidert, Daniel Leidert, Michael Banck
Date: 2014-09-22 10:00:17 UTC
mfrom: (1.1.9)
Revision ID: package-import@ubuntu.com-20140922100017-voentkpft033vxa6

Tags: 140902-1

* New upstream release.

[ Daniel Leidert ]
* debian/copyright: Updated.
* debian/upstream: Renamed to debian/upstream/metadata.

[ Michael Banck ]
* debian/travis.1: Add LAMMPS and DLPOLY to trajectory filetypes and
document -stream option.
* debian/control (Description): List supported fileteypes and
available analysis features.

files added:
debian/upstream/metadata

src/chiral.cpp

src/chiral.h

src/ir.cpp

src/ir.h

src/linalg.cpp

src/linalg.h

src/normalcoordinate.cpp

src/normalcoordinate.h

src/qr_fact.cpp

src/qr_fact.h

src/sdfmap.cpp

src/sdfmap.h

src/xdvectorn.cpp

src/xdvectorn.h

files removed:
src/cluster.cpp

src/cluster.h

src/fesa.cpp

src/fesa.h

src/ffgen.cpp

src/ffgen.h

src/ffparm.cpp

src/ltspec.cpp

src/ltspec.h

src/microhet.cpp

src/microhet.h

src/react.cpp

src/react.h

src/resp.cpp

src/resp.h

src/sdengine.cpp

src/sdengine.h

src/structurefactor.cpp

src/structurefactor.h

src/thermo.cpp

src/thermo.h

src/v_base.cpp

src/v_base.h

src/v_base_wl.cpp

src/v_c_loops.h

src/v_cell.cpp

src/v_cell.h

src/v_common.cpp

src/v_common.h

src/v_compute.cpp

src/v_compute.h

src/v_config.h

src/v_container.cpp

src/v_container.h

src/v_container_prd.cpp

src/v_container_prd.h

src/v_pre_container.cpp

src/v_pre_container.h

src/v_rad_option.h

src/v_unitcell.cpp

src/v_unitcell.h

src/v_wall.cpp

src/v_wall.h

src/v_worklist.h

src/voro++.h

src/voroanalysis.cpp

src/voroanalysis.h

src/vorowrapper.cpp

src/vorowrapper.h

files modified:
Makefile

debian/changelog

debian/control

debian/copyright

debian/travis.1

debian/upstream

src/2df.cpp

src/2df.h

src/3df.cpp

src/3df.h

src/acf.cpp

src/config.h

src/dacf.cpp

src/df.cpp

src/df.h

src/elementdata.cpp

src/gather.cpp

src/globalvar.cpp

src/globalvar.h

src/interface.cpp

src/interface.h

src/lmmin.cpp

src/maintools.cpp

src/maintools.h

src/moltools.cpp

src/moltools.h

src/plproj.cpp

src/plproj.h

src/raman.cpp

src/raman.h

src/reordyn.cpp

src/timestep.cpp

src/timestep.h

src/tools.cpp

src/tools.h

src/travis.cpp

src/travis.h

src/xdf.cpp

src/xdf.h

src/xdvector3.cpp

src/xmatrix3.cpp

src/xmatrix3.h

src/xobarray.cpp

src/xobarray.h

src/xvector3.h

Show diffs side-by-side

added added

removed removed

src/voroanalysis.cpp

/*****************************************************************************

TRAVIS - Trajectory Analyzer and Visualizer

http://www.travis-analyzer.de/

2012-2014 Martin Thomas

This file written by Martin Brehm.

This program is free software: you can redistribute it and/or modify

it under the terms of the GNU General Public License as published by

the Free Software Foundation, either version 3 of the License, or

(at your option) any later version.

This program is distributed in the hope that it will be useful,

but WITHOUT ANY WARRANTY; without even the implied warranty of

MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

GNU General Public License for more details.

You should have received a copy of the GNU General Public License

along with this program. If not, see <http://www.gnu.org/licenses/>.

*****************************************************************************/

#include "voroanalysis.h"

#include "globalvar.h"

#include "maintools.h"

CVoroPoint::CVoroPoint()

{

m_iaVoroCells.SetName("CVoroPoint::m_iaVoroCells");

}

CVoroPoint::~CVoroPoint()

{

}

/*CVoroFace::CVoroFace()

{

}

CVoroFace::~CVoroFace()

{

}*/

CDelaunayPoint::CDelaunayPoint()

{

}

CDelaunayPoint::~CDelaunayPoint()

{

}

CDelaunayTetrahedron::CDelaunayTetrahedron()

{

}

CDelaunayTetrahedron::~CDelaunayTetrahedron()

{

}

CVoroCell::CVoroCell()

{

m_iaVoroPoints.SetName("CVoroCell::m_iaVoroPoints");

m_iaNeighborCells.SetName("CVoroCell::m_iaNeighborCells");

}

CVoroCell::~CVoroCell()

{

}

CVoroAnalysis::CVoroAnalysis()

{

m_oaVoroCells.SetName("CVoroAnalysis::m_oaVoroCells");

m_oaDelaunayTetrahedra.SetName("CVoroAnalysis::m_oaDelaunayTetrahedra");

m_oaDelaunayPoints.SetName("CVoroAnalysis::m_oaDelaunayPoints");

m_oaVoroPoints.SetName("CVoroAnalysis::m_oaVoroPoints");

}

CVoroAnalysis::~CVoroAnalysis()

{

}

void CVoroAnalysis::Build(CTimeStep *ts, bool sanity, bool verbose)

{

voronoicell_neighbor c;

container_periodic_poly *con;

100

CVoroCell *vc;

101

int ijk, q, z, z2, z3, faces, fc, id, ti, i, cc;

102

double *pp;

103

vector<int> nb;

104

vector<int> fv;

105

CVoroPoint *vp;

106

CDelaunayTetrahedron *dt;

107

CDelaunayPoint *dp;

108

double px, py, pz, pk[4];

109

CxVector3 vec;

110

double tfa[8], tf;

111

112

#define EPS 0.000001

113

114

try { con = new container_periodic_poly(g_fBoxX/1000.0,0,g_fBoxY/1000.0,0,0,g_fBoxZ/1000.0,g_pVoroWrapper->m_iBlocksX,g_pVoroWrapper->m_iBlocksY,g_pVoroWrapper->m_iBlocksZ,g_iVoroMemory); } catch(...) { con = NULL; }

115

if (con == NULL) NewException((double)sizeof(container_periodic_poly),__FILE__,__LINE__,__PRETTY_FUNCTION__);

116

117

m_oaDelaunayPoints.SetMaxSize(g_iGesAtomCount);

118

m_oaVoroCells.SetSize(g_iGesAtomCount);

119

120

for (z=0;z<g_iGesAtomCount;z++)

121

{

122

try { dp = new CDelaunayPoint(); } catch(...) { dp = NULL; }

123

if (dp == NULL) NewException((double)sizeof(CDelaunayPoint),__FILE__,__LINE__,__PRETTY_FUNCTION__);

124

125

dp->m_vPos = ts->m_vaCoords[z];

126

m_oaDelaunayPoints.Add(dp);

127

con->put(z,ts->m_vaCoords[z][0]/1000.0,ts->m_vaCoords[z][1]/1000.0,ts->m_vaCoords[z][2]/1000.0,0.0005);

128

}

129

130

c_loop_all_periodic vl(*con);

131

132

if (sanity)

133

{

134

mprintf(" Voronoi decomposition:\n");

135

mprintf(WHITE," [");

136

}

137

138

cc = 0;

139

if (vl.start())

140

{

141

142

{

143

if (sanity)

144

if (fmod(cc,g_iGesAtomCount/60.0) < 1.0)

145

mprintf(WHITE,"#");

146

147

if (con->compute_cell(c,vl))

148

{

149

ijk=vl.ijk;

150

q=vl.q;

151

pp=con->p[ijk]+con->ps*q;

152

153

id = con->id[ijk][q];

154

155

try { vc = new CVoroCell(); } catch(...) { vc = NULL; }

156

if (vc == NULL) NewException((double)sizeof(CVoroCell),__FILE__,__LINE__,__PRETTY_FUNCTION__);

157

158

m_oaVoroCells[id] = vc;

159

vc->m_vPos[0] = ts->m_vaCoords[id][0]/1000.0;

160

vc->m_vPos[1] = ts->m_vaCoords[id][1]/1000.0;

161

vc->m_vPos[2] = ts->m_vaCoords[id][2]/1000.0;

162

163

c.neighbors(nb);

164

faces = c.number_of_faces();

165

c.face_vertices(fv);

166

167

ti = 0;

168

169

if (verbose)

170

{

171

mprintf("### Zelle %d\n",id+1);

172

mprintf(" %d Nachbarn: ",nb.size());

173

}

174

175

vc->m_iaNeighborCells.SetSize(nb.size());

176

for (z=0;z<(int)nb.size();z++)

177

{

178

if (verbose)

179

mprintf("%d, ",nb[z]+1);

180

vc->m_iaNeighborCells[z] = nb[z];

181

}

182

183

if (verbose)

184

{

185

mprintf("\n");

186

mprintf(" %d Flaechen, %d Vertices.\n",faces,fv.size());

187

}

188

189

i = 0;

190

for (z=0;z<faces;z++)

191

{

192

fc = fv[ti];

193

194

if (verbose)

195

mprintf(" - Flaeche %d: %d Punkte - ",z+1,fc);

196

197

for (z2=0;z2<fc;z2++)

198

{

199

if (verbose)

200

mprintf("%d, ",fv[ti+z2+1]+1);

201

if (i < fv[ti+z2+1])

202

i = fv[ti+z2+1];

203

}

204

if (verbose)

205

mprintf("\n");

206

207

ti += fv[ti]+1;

208

}

209

210

if (verbose)

211

{

212

mprintf(" %d Voronoi-Punkte.\n",i+1);

213

mprintf(" Base Point ist %.6f | %.6f | %.6f .\n",*pp,pp[1],pp[2]);

214

}

215

216

for (z=0;z<=i;z++)

217

{

218

px = *pp+c.pts[z*3]*0.5;

219

py = pp[1]+c.pts[z*3+1]*0.5;

220

pz = pp[2]+c.pts[z*3+2]*0.5;

221

222

while (px < 0) px += g_fBoxX/1000.0;

223

while (px >= g_fBoxX/1000.0) px -= g_fBoxX/1000.0;

224

while (py < 0) py += g_fBoxY/1000.0;

225

while (py >= g_fBoxY/1000.0) py -= g_fBoxY/1000.0;

226

while (pz < 0) pz += g_fBoxZ/1000.0;

227

while (pz >= g_fBoxZ/1000.0) pz -= g_fBoxZ/1000.0;

228

229

/* while (px < -0.5) px += 1.0;

230

while (px >= 0.5) px -= 1.0;

231

while (py < -0.5) py += 1.0;

232

while (py >= 0.5) py -= 1.0;

233

while (pz < -0.5) pz += 1.0;

234

while (pz >= 0.5) pz -= 1.0;*/

235

236

if (verbose)

237

mprintf(" %2d: %.6f | %.6f | %.6f\n",z+1,px,py,pz);

238

239

for (z2=0;z2<m_oaVoroPoints.GetSize();z2++)

240

{

241

vp = (CVoroPoint*)m_oaVoroPoints[z2];

242

if (fabs(vp->m_vPos[0]-px) > EPS)

243

continue;

244

if (fabs(vp->m_vPos[1]-py) > EPS)

245

continue;

246

if (fabs(vp->m_vPos[2]-pz) > EPS)

247

continue;

248

if (verbose)

249

mprintf(" Punkt bereits vorhanden: %d.\n",z2+1);

250

vc->m_iaVoroPoints.Add(z2);

251

goto _done;

252

}

253

vc->m_iaVoroPoints.Add(m_oaVoroPoints.GetSize());

254

255

try { vp = new CVoroPoint(); } catch(...) { vp = NULL; }

256

if (vp == NULL) NewException((double)sizeof(CVoroPoint),__FILE__,__LINE__,__PRETTY_FUNCTION__);

257

258

vp->m_vPos[0] = px;

259

vp->m_vPos[1] = py;

260

vp->m_vPos[2] = pz;

261

m_oaVoroPoints.Add(vp);

262

if (verbose)

263

mprintf(" * Fuege Punkt als %d hinzu.\n",m_oaVoroPoints.GetSize());

264

_done:;

265

}

266

}

267

cc++;

268

} while (vl.inc());

269

270

if (sanity)

271

m_iMaxVoroPoints = m_oaVoroPoints.GetSize() * 2;

272

273

if (verbose)

274

mprintf("Insgesamt %d verschiedene Voronoi-Punkte.\n\n",m_oaVoroPoints.GetSize());

275

276

if (sanity)

277

{

278

mprintf(WHITE,"]\n");

279

mprintf("\n Found %d different voronoi points.\n\n",m_oaVoroPoints.GetSize());

280

}

281

}

282

283

if (sanity)

284

{

285

mprintf(" Delaunay tetrahedron generation:\n");

286

mprintf(WHITE," [");

287

}

288

289

for (z=0;z<m_oaVoroPoints.GetSize();z++)

290

{

291

if (sanity)

292

if (fmod(z,m_oaVoroPoints.GetSize()/60.0) < 1.0)

293

mprintf(WHITE,"#");

294

295

vp = (CVoroPoint*)m_oaVoroPoints[z];

296

297

if (verbose)

298

mprintf("### Voronoi-Punkt %d: %.6f | %.6f | %.6f\n",z+1,vp->m_vPos[0],vp->m_vPos[1],vp->m_vPos[2]);

299

300

try { dt = new CDelaunayTetrahedron(); } catch(...) { dt = NULL; }

301

if (dt == NULL) NewException((double)sizeof(CDelaunayTetrahedron),__FILE__,__LINE__,__PRETTY_FUNCTION__);

302

303

ti = 0;

304

305

for (z2=0;z2<m_oaVoroCells.GetSize();z2++)

306

{

307

vc = (CVoroCell*)m_oaVoroCells[z2];

308

for (z3=0;z3<vc->m_iaVoroPoints.GetSize();z3++)

309

{

310

if (vc->m_iaVoroPoints[z3] == z)

311

{

312

if (verbose)

313

{

314

vec = vp->m_vPos - vc->m_vPos;

315

316

while (vec[0] < -g_fBoxX/2000.0) vec[0] += g_fBoxX/1000.0;

317

while (vec[0] >= g_fBoxX/2000.0) vec[0] -= g_fBoxX/1000.0;

318

while (vec[1] < -g_fBoxY/2000.0) vec[1] += g_fBoxY/1000.0;

319

while (vec[1] >= g_fBoxY/2000.0) vec[1] -= g_fBoxY/1000.0;

320

while (vec[2] < -g_fBoxZ/2000.0) vec[2] += g_fBoxZ/1000.0;

321

while (vec[2] >= g_fBoxZ/2000.0) vec[2] -= g_fBoxZ/1000.0;

322

323

/* while (vec[0] < -0.5) vec[0] += 1.0;

324

while (vec[0] >= 0.5) vec[0] -= 1.0;

325

while (vec[1] < -0.5) vec[1] += 1.0;

326

while (vec[1] >= 0.5) vec[1] -= 1.0;

327

while (vec[2] < -0.5) vec[2] += 1.0;

328

while (vec[2] >= 0.5) vec[2] -= 1.0;*/

329

330

px = vec.GetLength();

331

332

mprintf(" - In Zelle %d als Punkt %d, Abstand %.8f. ",z2+1,z3+1,px);

333

vc->m_vPos.Dump();

334

mprintf("\n");

335

}

336

337

if (ti < 4)

338

dt->m_iPointID[ti] = z2;

339

340

ti++;

341

342

}

343

}

344

}

345

346

if (ti != 4)

347

{

348

eprintf("\nError: %d instead of 4 points in tetrahedron %d.\n",ti,z+1);

349

for (z2=0;z2<m_oaVoroCells.GetSize();z2++)

350

{

351

vc = (CVoroCell*)m_oaVoroCells[z2];

352

for (z3=0;z3<vc->m_iaVoroPoints.GetSize();z3++)

353

{

354

if (vc->m_iaVoroPoints[z3] == z)

355

{

356

vec = vp->m_vPos - vc->m_vPos;

357

358

while (vec[0] < -g_fBoxX/2000.0) vec[0] += g_fBoxX/1000.0;

359

while (vec[0] >= g_fBoxX/2000.0) vec[0] -= g_fBoxX/1000.0;

360

while (vec[1] < -g_fBoxY/2000.0) vec[1] += g_fBoxY/1000.0;

361

while (vec[1] >= g_fBoxY/2000.0) vec[1] -= g_fBoxY/1000.0;

362

while (vec[2] < -g_fBoxZ/2000.0) vec[2] += g_fBoxZ/1000.0;

363

while (vec[2] >= g_fBoxZ/2000.0) vec[2] -= g_fBoxZ/1000.0;

364

365

px = vec.GetLength();

366

367

mprintf(" - In cell %d as point %d, distance %.8f. ",z2+1,z3+1,px);

368

vc->m_vPos.Dump();

369

mprintf("\n");

370

}

371

}

372

}

373

delete dt;

374

continue;

375

}

376

377

m_oaDelaunayTetrahedra.Add(dt);

378

379

if (verbose || sanity)

380

{

381

pk[0] = 1E10;

382

pk[1] = 1E10;

383

pk[2] = 1E10;

384

pk[3] = 1E10;

385

386

for (z2=0;z2<m_oaVoroCells.GetSize();z2++)

387

{

388

vc = (CVoroCell*)m_oaVoroCells[z2];

389

390

vec = vp->m_vPos - vc->m_vPos;

391

392

while (vec[0] < -g_fBoxX/2000.0) vec[0] += g_fBoxX/1000.0;

393

while (vec[0] >= g_fBoxX/2000.0) vec[0] -= g_fBoxX/1000.0;

394

while (vec[1] < -g_fBoxY/2000.0) vec[1] += g_fBoxY/1000.0;

395

while (vec[1] >= g_fBoxY/2000.0) vec[1] -= g_fBoxY/1000.0;

396

while (vec[2] < -g_fBoxZ/2000.0) vec[2] += g_fBoxZ/1000.0;

397

while (vec[2] >= g_fBoxZ/2000.0) vec[2] -= g_fBoxZ/1000.0;

398

399

/* while (vec[0] < -0.5) vec[0] += 1.0;

400

while (vec[0] >= 0.5) vec[0] -= 1.0;

401

while (vec[1] < -0.5) vec[1] += 1.0;

402

while (vec[1] >= 0.5) vec[1] -= 1.0;

403

while (vec[2] < -0.5) vec[2] += 1.0;

404

while (vec[2] >= 0.5) vec[2] -= 1.0;*/

405

406

px = vec.GetLength();

407

408

// px = FoldVector1(vp->m_vPos - vc->m_vPos).GetLength();

409

410

if (px < pk[3])

411

pk[3] = px;

412

if (pk[3] < pk[2])

413

{

414

px = pk[2];

415

pk[2] = pk[3];

416

pk[3] = px;

417

}

418

if (pk[2] < pk[1])

419

{

420

px = pk[1];

421

pk[1] = pk[2];

422

pk[2] = px;

423

}

424

if (pk[1] < pk[0])

425

{

426

px = pk[0];

427

pk[0] = pk[1];

428

pk[1] = px;

429

}

430

}

431

432

if (verbose)

433

mprintf(" Kleinste 4 Abstaende:\n %.8f\n %.8f\n %.8f\n %.8f\n",pk[0],pk[1],pk[2],pk[3]);

434

435

if (sanity || verbose)

436

{

437

for (z2=0;z2<4;z2++)

438

tfa[z2] = pk[z2];

439

440

for (z2=0;z2<4;z2++)

441

{

442

vc = (CVoroCell*)m_oaVoroCells[dt->m_iPointID[z2]];

443

for (z3=0;z3<vc->m_iaVoroPoints.GetSize();z3++)

444

{

445

if (vc->m_iaVoroPoints[z3] == z)

446

{

447

vec = vp->m_vPos - vc->m_vPos;

448

449

while (vec[0] < -g_fBoxX/2000.0) vec[0] += g_fBoxX/1000.0;

450

while (vec[0] >= g_fBoxX/2000.0) vec[0] -= g_fBoxX/1000.0;

451

while (vec[1] < -g_fBoxY/2000.0) vec[1] += g_fBoxY/1000.0;

452

while (vec[1] >= g_fBoxY/2000.0) vec[1] -= g_fBoxY/1000.0;

453

while (vec[2] < -g_fBoxZ/2000.0) vec[2] += g_fBoxZ/1000.0;

454

while (vec[2] >= g_fBoxZ/2000.0) vec[2] -= g_fBoxZ/1000.0;

455

456

/* while (vec[0] < -0.5) vec[0] += 1.0;

457

while (vec[0] >= 0.5) vec[0] -= 1.0;

458

while (vec[1] < -0.5) vec[1] += 1.0;

459

while (vec[1] >= 0.5) vec[1] -= 1.0;

460

while (vec[2] < -0.5) vec[2] += 1.0;

461

while (vec[2] >= 0.5) vec[2] -= 1.0;*/

462

463

tfa[z2+4] = vec.GetLength();

464

}

465

}

466

}

467

468

tf = MaxDiff_DoubleArray(tfa,8);

469

if (tf > EPS)

470

{

471

eprintf("\nSanity check failed: Largest Deviation is %.8f > EPS = %.8f.\n%.8f\n%.8f\n%.8f\n%.8f\n%.8f\n%.8f\n%.8f\n%.8f\n\n",tf,EPS,tfa[0],tfa[1],tfa[2],tfa[3],tfa[4],tfa[5],tfa[6],tfa[7]);

472

}

473

}

474

}

475

476

}

477

478

if (sanity)

479

{

480

mprintf(WHITE,"]\n");

481

mprintf("\n Test passed.\n\n");

482

}

483

484

for (z=0;z<m_oaDelaunayTetrahedra.GetSize();z++)

485

{

486

dt = (CDelaunayTetrahedron*)m_oaDelaunayTetrahedra[z];

487

488

dt->m_vCenter[0] = ((CVoroPoint*)m_oaVoroPoints[z])->m_vPos[0] * 1000.0;

489

dt->m_vCenter[1] = ((CVoroPoint*)m_oaVoroPoints[z])->m_vPos[1] * 1000.0;

490

dt->m_vCenter[2] = ((CVoroPoint*)m_oaVoroPoints[z])->m_vPos[2] * 1000.0;

491

492

// vec = dt->m_vCenter - ((CDelaunayPoint*)m_oaDelaunayPoints[dt->m_iPointID[0]])->m_vPos;

493

494

dt->m_fRadius = FoldVector(dt->m_vCenter - ((CDelaunayPoint*)m_oaDelaunayPoints[dt->m_iPointID[0]])->m_vPos).GetLength();

495

}

496

497

delete con;

498

}

499

500

501

void CVoroAnalysis::Parse()

502

{

503

CTimeStep *t;

504

CMolecule *m;

505

CSingleMolecule *sm;

506

int z;

507

508

mprintf(WHITE,">>> Void Analysis >>>\n\n");

509

510

mprintf(" Initializing...\n");

511

g_pVoroWrapper->Init();

512

mprintf("\n");

513

514

mprintf("*** Voro: Box density is %f particles / Angstrom^3.\n",g_pVoroWrapper->m_fBoxDens);

515

mprintf("*** Voro: Using %d x %d x %d blocks.\n",g_pVoroWrapper->m_iBlocksX,g_pVoroWrapper->m_iBlocksY,g_pVoroWrapper->m_iBlocksZ);

516

517

try { t = new CTimeStep(); } catch(...) { t = NULL; }

518

if (t == NULL) NewException((double)sizeof(CTimeStep),__FILE__,__LINE__,__PRETTY_FUNCTION__);

519

520

t->CopyFrom(&g_TimeStep);

521

522

m = (CMolecule*)g_oaMolecules[0];

523

sm = (CSingleMolecule*)g_oaSingleMolecules[m->m_laSingleMolIndex[0]];

524

525

t->CenterPos(t->m_vaCoords[((CxIntArray*)sm->m_oaAtomOffset[sm->m_baAtomIndex.GetSize()-1])->GetAt(1)]);

526

t->CenterPos(CxVector3(-g_fBoxX/2.0,-g_fBoxY/2.0,-g_fBoxZ/2.0));

527

t->FoldAtomsPositive();

528

g_pVoroWrapper->Dump("voro.txt",t);

529

mprintf("\n");

530

mprintf(" Voro++: Using cell memory for %d particles.\n\n",g_iVoroMemory);

531

532

mprintf(" Performing sanity check...\n\n");

533

Build(t,true,false);

534

535

// g_pVoroWrapper->WritePOV(t,"voro.pov",0,0);

536

delete t;

537

538

// AskYesNo(" This is hard-coded for Ulrikes Bmim-Br ^^ [accept] ",false);

539

540

m_bSphereHole = AskYesNo(" Do you want to create a spherical hole distribution function (y/n)? [yes] ",true);

541

542

if (m_bSphereHole)

543

{

544

try { m_pSphereHoleDF = new CDF(); } catch(...) { m_pSphereHoleDF = NULL; }

545

if (m_pSphereHoleDF == NULL) NewException((double)sizeof(CDF),__FILE__,__LINE__,__PRETTY_FUNCTION__);

546

547

m_pSphereHoleDF->m_fMinVal = 0;

548

m_pSphereHoleDF->m_fMaxVal = AskFloat(" Enter the maximal hole diameter (in pm): [1500.0] ",1500.0f);

549

m_pSphereHoleDF->m_iResolution = AskUnsignedInteger(" Enter the binning resolution: [%d] ",int(m_pSphereHoleDF->m_fMaxVal),int(m_pSphereHoleDF->m_fMaxVal));

550

m_pSphereHoleDF->Create();

551

m_pSphereHoleDF->SetLabelX("Hole diameter");

552

m_pSphereHoleDF->SetLabelY("Occurrence");

553

mprintf("\n");

554

}

555

556

m_bSphereHoleRDF = AskYesNo(" Do you want to create a spherical hole RDF (y/n)? [yes] ",true);

557

558

if (m_bSphereHoleRDF)

559

{

560

try { m_pSphereHoleRDF = new CDF(); } catch(...) { m_pSphereHoleRDF = NULL; }

561

if (m_pSphereHoleRDF == NULL) NewException((double)sizeof(CDF),__FILE__,__LINE__,__PRETTY_FUNCTION__);

562

563

m_pSphereHoleRDF->m_fMinVal = 0;

564

m_pSphereHoleRDF->m_fMaxVal = AskUnsignedInteger(" Enter the RDF range (in pm): [%d] ",HalfBox(),HalfBox());

565

m_pSphereHoleRDF->m_iResolution = AskUnsignedInteger(" Enter the RDF resolution: [300] ",300);

566

m_pSphereHoleRDF->Create();

567

mprintf("\n");

568

}

569

570

m_bSphereHole2DF = AskYesNo(" Do you want to create a spherical hole 2D DF (y/n)? [no] ",false);

571

572

if (m_bSphereHole2DF)

573

{

574

try { m_pSphereHole2DF = new C2DF(); } catch(...) { m_pSphereHole2DF = NULL; }

575

if (m_pSphereHole2DF == NULL) NewException((double)sizeof(C2DF),__FILE__,__LINE__,__PRETTY_FUNCTION__);

576

577

m_pSphereHole2DF->m_fMinVal[0] = 0;

578

m_pSphereHole2DF->m_fMaxVal[0] = AskUnsignedInteger(" Enter the RDF channel range (in pm): [%d] ",HalfBox(),HalfBox());

579

m_pSphereHole2DF->m_iRes[0] = AskUnsignedInteger(" Enter the RDF channel resolution: [100] ",100);

580

m_pSphereHole2DF->m_fMinVal[1] = 0;

581

m_pSphereHole2DF->m_fMaxVal[1] = AskUnsignedInteger(" Enter the hole radius range (in pm): [300] ",300);

582

m_pSphereHole2DF->m_iRes[1] = AskUnsignedInteger(" Enter the hole radius channel resolution: [50] ",50);

583

m_pSphereHole2DF->SetLabelX("Void distance [pm]");

584

m_pSphereHole2DF->SetLabelY("Minimal void radius [pm]");

585

m_pSphereHole2DF->SetLabelZ("Occurence");

586

m_pSphereHole2DF->Create();

587

588

try { m_pSphereHole2DFCounter = new unsigned long[m_pSphereHole2DF->m_iRes[1]]; } catch(...) { m_pSphereHole2DFCounter = NULL; }

589

if (m_pSphereHole2DFCounter == NULL) NewException((double)m_pSphereHole2DF->m_iRes[1]*sizeof(unsigned long),__FILE__,__LINE__,__PRETTY_FUNCTION__);

590

591

for (z=0;z<m_pSphereHole2DF->m_iRes[1];z++)

592

m_pSphereHole2DFCounter[z] = 0;

593

594

mprintf("\n");

595

}

596

597

m_bAllAtomRDF = AskYesNo(" Do you want to create an all-atom RDF (y/n)? [no] ",false);

598

599

if (m_bAllAtomRDF)

600

{

601

try { m_pAllAtomRDF = new CDF(); } catch(...) { m_pAllAtomRDF = NULL; }

602

if (m_pAllAtomRDF == NULL) NewException((double)sizeof(CDF),__FILE__,__LINE__,__PRETTY_FUNCTION__);

603

604

m_pAllAtomRDF->m_fMinVal = 0;

605

m_pAllAtomRDF->m_fMaxVal = AskUnsignedInteger(" Enter the RDF range (in pm): [%d] ",HalfBox(),HalfBox());

606

m_pAllAtomRDF->m_iResolution = AskUnsignedInteger(" Enter the RDF resolution: [300] ",300);

607

m_pAllAtomRDF->Create();

608

mprintf("\n");

609

}

610

611

612

m_bHoleDump = AskYesNo(" Do you want to write out a XYZ hole trajectory (y/n)? [no] ",false);

613

614

if (m_bHoleDump)

615

{

616

m_fDumpMinHole = AskFloat(" Enter the minimal hole diameter to output (in pm): [400] ",400) / 2.0;

617

m_fHoleDump = OpenFileWrite("holedump.xyz",true);

618

}

619

620

m_bSphereHoleSDF = AskYesNo(" Create a spherical hole SDF (y/n)? [no] ",false);

621

622

if (m_bSphereHoleSDF)

623

{

624

g_bVoidSDF = true;

625

626

try { m_pSphereHoleSDF = new CSDF(); } catch(...) { m_pSphereHoleSDF = NULL; }

627

if (m_pSphereHoleSDF == NULL) NewException((double)sizeof(CSDF),__FILE__,__LINE__,__PRETTY_FUNCTION__);

628

629

m_pSphereHoleSDF->Parse(true);

630

631

m_pSphereHoleSDF->m_pSDF->m_fMinVal[0] = -m_pSphereHoleSDF->m_fRadius;

632

m_pSphereHoleSDF->m_pSDF->m_fMaxVal[0] = m_pSphereHoleSDF->m_fRadius;

633

m_pSphereHoleSDF->m_pSDF->m_fMinVal[1] = -m_pSphereHoleSDF->m_fRadius;

634

m_pSphereHoleSDF->m_pSDF->m_fMaxVal[1] = m_pSphereHoleSDF->m_fRadius;

635

m_pSphereHoleSDF->m_pSDF->m_fMinVal[2] = -m_pSphereHoleSDF->m_fRadius;

636

m_pSphereHoleSDF->m_pSDF->m_fMaxVal[2] = m_pSphereHoleSDF->m_fRadius;

637

m_pSphereHoleSDF->m_pSDF->m_iRes[0] = m_pSphereHoleSDF->m_iResolution;

638

m_pSphereHoleSDF->m_pSDF->m_iRes[1] = m_pSphereHoleSDF->m_iResolution;

639

m_pSphereHoleSDF->m_pSDF->m_iRes[2] = m_pSphereHoleSDF->m_iResolution;

640

m_pSphereHoleSDF->m_pSDF->m_iHistogramRes = m_pSphereHoleSDF->m_iHistogramRes;

641

m_pSphereHoleSDF->m_pSDF->Create();

642

}

643

644

m_bEmptySpaceSDF = AskYesNo(" Create an empty space SDF (y/n)? [no] ",false);

645

646

if (m_bEmptySpaceSDF)

647

{

648

g_bVoidSDF = true;

649

650

m_bNewEmptyMode = AskYesNo(" Use new method (y/n)? [no] ",false);

651

652

if (m_bNewEmptyMode)

653

{

654

m_iTempRes = AskUnsignedInteger(" Enter resoultion of the temporary SDF: [300] ",300);

655

656

try { m_pTempSDF = new C3DF(); } catch(...) { m_pTempSDF = NULL; }

657

if (m_pTempSDF == NULL) NewException((double)sizeof(C3DF),__FILE__,__LINE__,__PRETTY_FUNCTION__);

658

659

m_pTempSDF->m_fMinVal[0] = 0;

660

m_pTempSDF->m_fMaxVal[0] = g_fBoxX;

661

m_pTempSDF->m_fMinVal[1] = 0;

662

m_pTempSDF->m_fMaxVal[1] = g_fBoxY;

663

m_pTempSDF->m_fMinVal[2] = 0;

664

m_pTempSDF->m_fMaxVal[2] = g_fBoxZ;

665

m_pTempSDF->m_iRes[0] = m_iTempRes;

666

m_pTempSDF->m_iRes[1] = m_iTempRes;

667

m_pTempSDF->m_iRes[2] = m_iTempRes;

668

m_pTempSDF->m_iHistogramRes = m_pEmptySpaceSDF->m_iHistogramRes;

669

} else

670

{

671

m_bEmptySpaceSDF_RefMol = AskYesNo(" Include the atoms of the reference molecule (y/n)? [yes] ",true);

672

m_iEmptySpaceBorder = AskUnsignedInteger(" Enter the size of the border in bins: [10] ",10);

673

674

try { m_pTempSDF = new C3DF(); } catch(...) { m_pTempSDF = NULL; }

675

if (m_pTempSDF == NULL) NewException((double)sizeof(C3DF),__FILE__,__LINE__,__PRETTY_FUNCTION__);

676

677

m_pTempSDF->m_fMinVal[0] = -m_pEmptySpaceSDF->m_fRadius;

678

m_pTempSDF->m_fMaxVal[0] = m_pEmptySpaceSDF->m_fRadius;

679

m_pTempSDF->m_fMinVal[1] = -m_pEmptySpaceSDF->m_fRadius;

680

m_pTempSDF->m_fMaxVal[1] = m_pEmptySpaceSDF->m_fRadius;

681

m_pTempSDF->m_fMinVal[2] = -m_pEmptySpaceSDF->m_fRadius;

682

m_pTempSDF->m_fMaxVal[2] = m_pEmptySpaceSDF->m_fRadius;

683

m_pTempSDF->m_iRes[0] = m_pEmptySpaceSDF->m_iResolution;

684

m_pTempSDF->m_iRes[1] = m_pEmptySpaceSDF->m_iResolution;

685

m_pTempSDF->m_iRes[2] = m_pEmptySpaceSDF->m_iResolution;

686

m_pTempSDF->m_iHistogramRes = m_pEmptySpaceSDF->m_iHistogramRes;

687

}

688

m_pTempSDF->Create();

689

690

try { m_pEmptySpaceSDF = new CSDF(); } catch(...) { m_pEmptySpaceSDF = NULL; }

691

if (m_pEmptySpaceSDF == NULL) NewException((double)sizeof(CSDF),__FILE__,__LINE__,__PRETTY_FUNCTION__);

692

693

m_pEmptySpaceSDF->Parse(true);

694

695

m_pEmptySpaceSDF->m_pSDF->m_fMinVal[0] = -m_pEmptySpaceSDF->m_fRadius;

696

m_pEmptySpaceSDF->m_pSDF->m_fMaxVal[0] = m_pEmptySpaceSDF->m_fRadius;

697

m_pEmptySpaceSDF->m_pSDF->m_fMinVal[1] = -m_pEmptySpaceSDF->m_fRadius;

698

m_pEmptySpaceSDF->m_pSDF->m_fMaxVal[1] = m_pEmptySpaceSDF->m_fRadius;

699

m_pEmptySpaceSDF->m_pSDF->m_fMinVal[2] = -m_pEmptySpaceSDF->m_fRadius;

700

m_pEmptySpaceSDF->m_pSDF->m_fMaxVal[2] = m_pEmptySpaceSDF->m_fRadius;

701

m_pEmptySpaceSDF->m_pSDF->m_iRes[0] = m_pEmptySpaceSDF->m_iResolution;

702

m_pEmptySpaceSDF->m_pSDF->m_iRes[1] = m_pEmptySpaceSDF->m_iResolution;

703

m_pEmptySpaceSDF->m_pSDF->m_iRes[2] = m_pEmptySpaceSDF->m_iResolution;

704

m_pEmptySpaceSDF->m_pSDF->m_iHistogramRes = m_pEmptySpaceSDF->m_iHistogramRes;

705

m_pEmptySpaceSDF->m_pSDF->Create();

706

}

707

708

mprintf(WHITE,"\n<<< End of Void Analysis <<<\n\n");

709

}

710

711

712

void CVoroAnalysis::Step(CTimeStep *ts)

713

{

714

int z, z2, z3, z4, ti;

715

CxVec3Array tva;

716

CDelaunayTetrahedron *dt, *dt2;

717

CMolecule *m;

718

CSingleMolecule *sm;

719

double tf, tf2, tfx;

720

721

722

if (m_bSphereHole || m_bSphereHoleRDF || m_bSphereHole2DF || m_bHoleDump || m_bSphereHoleRDF)

723

{

724

Clean();

725

726

Build(ts,false,false);

727

}

728

729

if (m_bSphereHole)

730

{

731

for (z=0;z<m_oaDelaunayTetrahedra.GetSize();z++)

732

{

733

dt = (CDelaunayTetrahedron*)m_oaDelaunayTetrahedra[z];

734

m_pSphereHoleDF->AddToBin(dt->m_fRadius*2.0);

735

736

if (dt->m_fRadius > 250.0)

737

tva.Add(dt->m_vCenter);

738

}

739

}

740

741

if (m_bSphereHoleRDF)

742

{

743

for (z=0;z<m_oaDelaunayTetrahedra.GetSize();z++)

744

{

745

dt = (CDelaunayTetrahedron*)m_oaDelaunayTetrahedra[z];

746

747

for (z2=z+1;z2<m_oaDelaunayTetrahedra.GetSize();z2++)

748

{

749

dt2 = (CDelaunayTetrahedron*)m_oaDelaunayTetrahedra[z2];

750

751

tf = FoldVector(dt->m_vCenter - dt2->m_vCenter).GetLength();

752

753

m_pSphereHoleRDF->AddToBin(tf);

754

}

755

}

756

}

757

758

if (m_bSphereHole2DF)

759

{

760

tfx = 1.0 / m_pSphereHole2DF->m_iRes[1] * (m_pSphereHole2DF->m_fMaxVal[1]-m_pSphereHole2DF->m_fMinVal[1]);

761

762

for (z2=0;z2<m_oaDelaunayTetrahedra.GetSize();z2++)

763

{

764

dt = (CDelaunayTetrahedron*)m_oaDelaunayTetrahedra[z2];

765

766

for (z3=z2+1;z3<m_oaDelaunayTetrahedra.GetSize();z3++)

767

{

768

dt2 = (CDelaunayTetrahedron*)m_oaDelaunayTetrahedra[z3];

769

770

tf2 = FoldVector(dt->m_vCenter - dt2->m_vCenter).GetLength();

771

772

tf = m_pSphereHole2DF->m_fMinVal[1];

773

774

for (z=0;z<m_pSphereHole2DF->m_iRes[1];z++)

775

{

776

tf += tfx;

777

778

if (dt->m_fRadius < tf)

779

continue;

780

781

if (dt2->m_fRadius < tf)

782

continue;

783

784

m_pSphereHole2DFCounter[z]++;

785

786

m_pSphereHole2DF->AddToBin(tf2,tf);

787

}

788

}

789

}

790

}

791

792

if (m_bAllAtomRDF)

793

{

794

for (z=0;z<g_iGesAtomCount;z++)

795

{

796

for (z2=z+1;z2<g_iGesAtomCount;z2++)

797

{

798

tf = FoldVector(ts->m_vaCoords[z] - ts->m_vaCoords[z2]).GetLength();

799

m_pAllAtomRDF->AddToBin(tf);

800

}

801

}

802

}

803

804

if (m_bHoleDump)

805

{

806

mfprintf(m_fHoleDump," %d\n",g_iGesAtomCount+m_iMaxVoroPoints);

807

mfprintf(m_fHoleDump,"\n");

808

809

for (z=0;z<g_oaMolecules.GetSize();z++)

810

{

811

m = (CMolecule*)g_oaMolecules[z];

812

if (m->m_bPseudo)

813

continue;

814

for (z2=0;z2<m->m_laSingleMolIndex.GetSize();z2++)

815

{

816

sm = (CSingleMolecule*)g_oaSingleMolecules[m->m_laSingleMolIndex[z2]];

817

for (z3=0;z3<m->m_baAtomIndex.GetSize();z3++)

818

{

819

if ((!g_bSaveVirtAtoms) && (m->m_baAtomIndex[z3] == g_iVirtAtomType))

820

continue;

821

for (z4=0;z4<((CxIntArray*)sm->m_oaAtomOffset[z3])->GetSize();z4++)

822

mfprintf(m_fHoleDump," %s %8.5f %8.5f %8.5f\n",((CAtom*)g_oaAtoms[m->m_baAtomIndex[z3]])->m_sName,ts->m_vaCoords[((CxIntArray*)sm->m_oaAtomOffset[z3])->GetAt(z4)][0]/100.0f,ts->m_vaCoords[((CxIntArray*)sm->m_oaAtomOffset[z3])->GetAt(z4)][1]/100.0f,ts->m_vaCoords[((CxIntArray*)sm->m_oaAtomOffset[z3])->GetAt(z4)][2]/100.0f);

823

}

824

}

825

}

826

827

ti = 0;

828

for (z=0;z<m_oaDelaunayTetrahedra.GetSize();z++)

829

{

830

dt = (CDelaunayTetrahedron*)m_oaDelaunayTetrahedra[z];

831

if (dt->m_fRadius > m_fDumpMinHole)

832

{

833

ti++;

834

mfprintf(m_fHoleDump," He %8.5f %8.5f %8.5f\n",dt->m_vCenter[0]/100.0,dt->m_vCenter[1]/100.0,dt->m_vCenter[2]/100.0);

835

}

836

}

837

for (z=ti;z<m_iMaxVoroPoints;z++)

838

mfprintf(m_fHoleDump," He %8.5f %8.5f %8.5f\n",ts->m_vaCoords[0][0]/100.0f,ts->m_vaCoords[0][1]/100.0f,ts->m_vaCoords[0][2]/100.0f);

839

}

840

}

841

842

843

void CVoroAnalysis::Finish()

844

{

845

int z, z2;

846

C3DF *tempSDF;

847

char buf[256];

848

849

if (m_bSphereHole)

850

{

851

m_pSphereHoleDF->NormBinIntegral(100000.0);

852

m_pSphereHoleDF->Write("","sphere_hole_df.csv","",false);

853

}

854

855

if (m_bSphereHoleRDF)

856

{

857

m_pSphereHoleRDF->CorrectRadialDist();

858

m_pSphereHoleRDF->MultiplyBin(g_fBoxX*g_fBoxY*g_fBoxZ / (4.0/3.0*Pi) / (m_pSphereHoleRDF->m_fBinEntries+m_pSphereHoleRDF->m_fSkipEntries));

859

m_pSphereHoleRDF->Write("","sphere_hole_rdf.csv","",false);

860

}

861

862

if (m_bSphereHole2DF)

863

{

864

m_pSphereHole2DF->CorrectRadialDist(0);

865

for (z=0;z<m_pSphereHole2DF->m_iRes[1];z++)

866

{

867

if (m_pSphereHole2DFCounter[z] != 0)

868

{

869

for (z2=0;z2<m_pSphereHole2DF->m_iRes[0];z2++)

870

m_pSphereHole2DF->m_pBin[z2+m_pSphereHole2DF->m_iRes[0]*z] *= 10000.0 / m_pSphereHole2DFCounter[z];

871

}

872

}

873

m_pSphereHole2DF->Log();

874

// m_pSphereHoleRDF->MultiplyBin(g_fBoxX*g_fBoxY*g_fBoxZ / (4.0/3.0*Pi) / (m_pSphereHoleRDF->m_fBinEntries+m_pSphereHoleRDF->m_fSkipEntries));

875

m_pSphereHole2DF->WriteMathematicaNb("","sphere_hole_2df.nb","",false);

876

m_pSphereHole2DF->WriteGnuplotInput("","sphere_hole_2df","",false);

877

}

878

879

if (m_bAllAtomRDF)

880

{

881

m_pAllAtomRDF->CorrectRadialDist();

882

m_pAllAtomRDF->MultiplyBin(g_fBoxX*g_fBoxY*g_fBoxZ / (4.0/3.0*Pi) / g_iSteps / g_iGesAtomCount / g_iGesAtomCount * 2.0 * g_iStride);

883

m_pAllAtomRDF->Write("","all_atom_rdf.csv","",false);

884

}

885

886

if (m_bHoleDump)

887

{

888

fclose(m_fHoleDump);

889

}

890

891

if (m_bEmptySpaceSDF)

892

{

893

mprintf(WHITE,"*** Empty Space SDF ***\n");

894

m_pEmptySpaceSDF->m_pSDF->MultiplyBin(pow(m_pEmptySpaceSDF->m_iResolution/m_pEmptySpaceSDF->m_fRadius*1000.0,3) / (double)g_iSteps / ((CMolecule*)g_oaMolecules[g_iFixMol])->m_laSingleMolIndex.GetSize() * (g_bDoubleBox?g_iDoubleBoxFactor:1));

895

for (z2=0;z2<=g_iSDFSmoothGrade;z2++)

896

{

897

try { tempSDF = new C3DF(); } catch(...) { tempSDF = NULL; }

898

if (tempSDF == NULL) NewException((double)sizeof(C3DF),__FILE__,__LINE__,__PRETTY_FUNCTION__);

899

900

tempSDF->CopyFrom(m_pEmptySpaceSDF->m_pSDF);

901

if (z2 != 0)

902

{

903

tempSDF->Smooth(z2);

904

sprintf(buf,".s%d.plt",z2);

905

} else sprintf(buf,".plt");

906

mprintf(" Saving SDF as \"sdf_emptyspace%s\"...\n",buf);

907

tempSDF->WritePLT("sdf_","emptyspace",buf,true);

908

909

if (z2 != 0)

910

sprintf(buf,".s%d.cube",z2);

911

else sprintf(buf,".cube");

912

mprintf(" Saving SDF as \"sdf_emptyspace%s\"...\n",buf);

913

tempSDF->WriteCube("sdf_","emptyspace",buf,true);

914

}

915

}

916

917

if (m_bSphereHoleSDF)

918

{

919

mprintf(WHITE,"*** Spherical Hole SDF ***\n");

920

m_pSphereHoleSDF->m_pSDF->MultiplyBin(pow(m_pSphereHoleSDF->m_iResolution/m_pSphereHoleSDF->m_fRadius*1000.0,3) / (double)g_iSteps / ((CMolecule*)g_oaMolecules[g_iFixMol])->m_laSingleMolIndex.GetSize() * (g_bDoubleBox?g_iDoubleBoxFactor:1));

921

for (z2=0;z2<=g_iSDFSmoothGrade;z2++)

922

{

923

try { tempSDF = new C3DF(); } catch(...) { tempSDF = NULL; }

924

if (tempSDF == NULL) NewException((double)sizeof(C3DF),__FILE__,__LINE__,__PRETTY_FUNCTION__);

925

926

tempSDF->CopyFrom(m_pSphereHoleSDF->m_pSDF);

927

if (z2 != 0)

928

{

929

tempSDF->Smooth(z2);

930

sprintf(buf,".s%d.plt",z2);

931

} else sprintf(buf,".plt");

932

mprintf(" Saving SDF as \"sdf_spherehole%s\"...\n",buf);

933

tempSDF->WritePLT("sdf_","spherehole",buf,true);

934

935

if (z2 != 0)

936

sprintf(buf,".s%d.cube",z2);

937

else sprintf(buf,".cube");

938

mprintf(" Saving SDF as \"sdf_spherehole%s\"...\n",buf);

939

tempSDF->WriteCube("sdf_","spherehole",buf,true);

940

}

941

}

942

}

943

944

945

void CVoroAnalysis::Clean()

946

{

947

int z;

948

949

for (z=0;z<m_oaDelaunayPoints.GetSize();z++)

950

delete (CDelaunayPoint*)m_oaDelaunayPoints[z];

951

m_oaDelaunayPoints.RemoveAll_KeepSize();

952

953

for (z=0;z<m_oaDelaunayTetrahedra.GetSize();z++)

954

delete (CDelaunayTetrahedron*)m_oaDelaunayTetrahedra[z];

955

m_oaDelaunayTetrahedra.RemoveAll_KeepSize();

956

957

for (z=0;z<m_oaVoroCells.GetSize();z++)

958

delete (CVoroCell*)m_oaVoroCells[z];

959

m_oaVoroCells.RemoveAll_KeepSize();

960

961

for (z=0;z<m_oaVoroPoints.GetSize();z++)

962

delete (CVoroPoint*)m_oaVoroPoints[z];

963

m_oaVoroPoints.RemoveAll_KeepSize();

964

}

965

966

967

void CVoroAnalysis::BinEmptySDF()

968

{

969

int x, y, z, ti, ti2, ti3;

970

971

for (z=m_iEmptySpaceBorder;z<m_pTempSDF->m_iRes[2]-m_iEmptySpaceBorder;z++)

972

{

973

ti = z*m_pTempSDF->m_iResXY;

974

for (y=m_iEmptySpaceBorder;y<m_pTempSDF->m_iRes[1]-m_iEmptySpaceBorder;y++)

975

{

976

ti2 = ti + y * m_pTempSDF->m_iRes[0];

977

for (x=m_iEmptySpaceBorder;x<m_pTempSDF->m_iRes[0]-m_iEmptySpaceBorder;x++)

978

{

979

ti3 = ti2 + x;

980

if (m_pTempSDF->m_pBin[ti3] == 0)

981

m_pEmptySpaceSDF->m_pSDF->m_pBin[ti3]++;

982

}

983

}

984

}

985

}

986

987

988

void CVoroAnalysis::BinEmptySDF_New(CxMatrix3 *m, CxVector3 *c)

989

{

990

int x, y, z, ti, ti2, ti3;

991

CxVector3 v, v2;

992

993

for (z=0;z<m_pEmptySpaceSDF->m_pSDF->m_iRes[2];z++)

994

{

995

ti = z*m_pEmptySpaceSDF->m_pSDF->m_iResXY;

996

v[2] = ((double)z)/m_pEmptySpaceSDF->m_pSDF->m_iRes[2]*(m_pEmptySpaceSDF->m_pSDF->m_fMaxVal[2]-m_pEmptySpaceSDF->m_pSDF->m_fMinVal[2])+m_pEmptySpaceSDF->m_pSDF->m_fMinVal[2] + (*c)[2];

997

for (y=0;y<m_pEmptySpaceSDF->m_pSDF->m_iRes[1];y++)

998

{

999

ti2 = ti + y * m_pEmptySpaceSDF->m_pSDF->m_iRes[0];

1000

v[1] = ((double)y)/m_pEmptySpaceSDF->m_pSDF->m_iRes[1]*(m_pEmptySpaceSDF->m_pSDF->m_fMaxVal[1]-m_pEmptySpaceSDF->m_pSDF->m_fMinVal[1])+m_pEmptySpaceSDF->m_pSDF->m_fMinVal[1] + (*c)[1];

1001

for (x=0;x<m_pEmptySpaceSDF->m_pSDF->m_iRes[0];x++)

1002

{

1003

ti3 = ti2 + x;

1004

v[0] = ((double)x)/m_pEmptySpaceSDF->m_pSDF->m_iRes[0]*(m_pEmptySpaceSDF->m_pSDF->m_fMaxVal[0]-m_pEmptySpaceSDF->m_pSDF->m_fMinVal[0])+m_pEmptySpaceSDF->m_pSDF->m_fMinVal[0] + (*c)[0];

1005

1006

v2 = *m * v;

1007

1008

// mprintf(" %f %f %f --> %f %f %f\n",v[0],v[1],v[2],v2[0],v2[1],v2[2]);

1009

1010

while (v2[0] >= g_fBoxX)

1011

v2[0] -= g_fBoxX;

1012

while (v2[0] < 0)

1013

v2[0] += g_fBoxX;

1014

while (v2[1] >= g_fBoxY)

1015

v2[1] -= g_fBoxY;

1016

while (v2[1] < 0)

1017

v2[1] += g_fBoxY;

1018

while (v2[2] >= g_fBoxZ)

1019

v2[2] -= g_fBoxZ;

1020

while (v2[2] < 0)

1021

v2[2] += g_fBoxZ;

1022

1023

// mprintf(" %f %f %f\n",v2[0],v2[1],v2[2]);

1024

1025

if (m_pTempSDF->IsZero(v2[0],v2[1],v2[2]))

1026

m_pEmptySpaceSDF->m_pSDF->m_pBin[ti3]++;

1027

}

1028

}

1029

}

1030

// m_pEmptySpaceSDF->m_pSDF->WritePLT("","bla2.plt","",true);

1031

// abort();

1032

}

1033

1034

1035

void CVoroAnalysis::BinSphereHoleSDF(CxMatrix3 *m, CxVector3 *c)

1036

{

1037

int z;

1038

CDelaunayTetrahedron *dt;

1039

CxVector3 vec;

1040

1041

m_pTempVA->RemoveAll_KeepSize();

1042

1043

for (z=0;z<m_oaDelaunayTetrahedra.GetSize();z++)

1044

{

1045

dt = (CDelaunayTetrahedron*)m_oaDelaunayTetrahedra[z];

1046

vec = dt->m_vCenter - *c;

1047

1048

while (vec[0] < -g_fBoxX/2.0) vec[0] += g_fBoxX;

1049

while (vec[0] >= g_fBoxX/2.0) vec[0] -= g_fBoxX;

1050

while (vec[1] < -g_fBoxY/2.0) vec[1] += g_fBoxY;

1051

while (vec[1] >= g_fBoxY/2.0) vec[1] -= g_fBoxY;

1052

while (vec[2] < -g_fBoxZ/2.0) vec[2] += g_fBoxZ;

1053

while (vec[2] >= g_fBoxZ/2.0) vec[2] -= g_fBoxZ;

1054

1055

m_pTempVA->Add(vec);

1056

}

1057

1058

for (z=0;z<m_oaDelaunayTetrahedra.GetSize();z++)

1059

{

1060

vec = *m * (*m_pTempVA)[z];

1061

1062

if (m_pSphereHoleSDF->m_bVdWSpheres)

1063

m_pSphereHoleSDF->m_pSDF->AddToBin_Sphere(vec,((CDelaunayTetrahedron*)m_oaDelaunayTetrahedra[z])->m_fRadius);

1064

else m_pSphereHoleSDF->m_pSDF->AddToBin(vec);

1065

}

1066

1067

}

Older »