source: src/boundary.cpp@ af374d

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults AddingActions_SaveParseParticleParameters Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests AtomFragments Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_BoundInBox_CenterInBox_MoleculeActions Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_PopActions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes Gui_displays_atomic_force_velocity ImplicitCharges IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix MoreRobust_FragmentAutomation ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable
Last change on this file since af374d was af374d, checked in by Frederik Heber <heber@…>, 15 years ago

Correction of degenerated polygons is working.

Signed-off-by: Frederik Heber <heber@…>

  • Property mode set to 100755
File size: 49.2 KB
Line 
1/** \file boundary.cpp
2 *
3 * Implementations and super-function for envelopes
4 */
5
6#include "atom.hpp"
7#include "bond.hpp"
8#include "boundary.hpp"
9#include "config.hpp"
10#include "element.hpp"
11#include "helpers.hpp"
12#include "info.hpp"
13#include "linkedcell.hpp"
14#include "log.hpp"
15#include "memoryallocator.hpp"
16#include "molecule.hpp"
17#include "tesselation.hpp"
18#include "tesselationhelpers.hpp"
19
20#include<gsl/gsl_poly.h>
21
22// ========================================== F U N C T I O N S =================================
23
24
25/** Determines greatest diameters of a cluster defined by its convex envelope.
26 * Looks at lines parallel to one axis and where they intersect on the projected planes
27 * \param *out output stream for debugging
28 * \param *BoundaryPoints NDIM set of boundary points defining the convex envelope on each projected plane
29 * \param *mol molecule structure representing the cluster
30 * \param *&TesselStruct Tesselation structure with triangles
31 * \param IsAngstroem whether we have angstroem or atomic units
32 * \return NDIM array of the diameters
33 */
34double *GetDiametersOfCluster(const Boundaries *BoundaryPtr, const molecule *mol, Tesselation *&TesselStruct, const bool IsAngstroem)
35{
36 Info FunctionInfo(__func__);
37 // get points on boundary of NULL was given as parameter
38 bool BoundaryFreeFlag = false;
39 double OldComponent = 0.;
40 double tmp = 0.;
41 double w1 = 0.;
42 double w2 = 0.;
43 Vector DistanceVector;
44 Vector OtherVector;
45 int component = 0;
46 int Othercomponent = 0;
47 Boundaries::const_iterator Neighbour;
48 Boundaries::const_iterator OtherNeighbour;
49 double *GreatestDiameter = new double[NDIM];
50
51 const Boundaries *BoundaryPoints;
52 if (BoundaryPtr == NULL) {
53 BoundaryFreeFlag = true;
54 BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
55 } else {
56 BoundaryPoints = BoundaryPtr;
57 Log() << Verbose(0) << "Using given boundary points set." << endl;
58 }
59 // determine biggest "diameter" of cluster for each axis
60 for (int i = 0; i < NDIM; i++)
61 GreatestDiameter[i] = 0.;
62 for (int axis = 0; axis < NDIM; axis++)
63 { // regard each projected plane
64 //Log() << Verbose(1) << "Current axis is " << axis << "." << endl;
65 for (int j = 0; j < 2; j++)
66 { // and for both axis on the current plane
67 component = (axis + j + 1) % NDIM;
68 Othercomponent = (axis + 1 + ((j + 1) & 1)) % NDIM;
69 //Log() << Verbose(1) << "Current component is " << component << ", Othercomponent is " << Othercomponent << "." << endl;
70 for (Boundaries::const_iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
71 //Log() << Verbose(1) << "Current runner is " << *(runner->second.second) << "." << endl;
72 // seek for the neighbours pair where the Othercomponent sign flips
73 Neighbour = runner;
74 Neighbour++;
75 if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
76 Neighbour = BoundaryPoints[axis].begin();
77 DistanceVector.CopyVector(&runner->second.second->x);
78 DistanceVector.SubtractVector(&Neighbour->second.second->x);
79 do { // seek for neighbour pair where it flips
80 OldComponent = DistanceVector.x[Othercomponent];
81 Neighbour++;
82 if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
83 Neighbour = BoundaryPoints[axis].begin();
84 DistanceVector.CopyVector(&runner->second.second->x);
85 DistanceVector.SubtractVector(&Neighbour->second.second->x);
86 //Log() << Verbose(2) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;
87 } while ((runner != Neighbour) && (fabs(OldComponent / fabs(
88 OldComponent) - DistanceVector.x[Othercomponent] / fabs(
89 DistanceVector.x[Othercomponent])) < MYEPSILON)); // as long as sign does not flip
90 if (runner != Neighbour) {
91 OtherNeighbour = Neighbour;
92 if (OtherNeighbour == BoundaryPoints[axis].begin()) // make it wrap around
93 OtherNeighbour = BoundaryPoints[axis].end();
94 OtherNeighbour--;
95 //Log() << Verbose(1) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl;
96 // now we have found the pair: Neighbour and OtherNeighbour
97 OtherVector.CopyVector(&runner->second.second->x);
98 OtherVector.SubtractVector(&OtherNeighbour->second.second->x);
99 //Log() << Verbose(1) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl;
100 //Log() << Verbose(1) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl;
101 // do linear interpolation between points (is exact) to extract exact intersection between Neighbour and OtherNeighbour
102 w1 = fabs(OtherVector.x[Othercomponent]);
103 w2 = fabs(DistanceVector.x[Othercomponent]);
104 tmp = fabs((w1 * DistanceVector.x[component] + w2
105 * OtherVector.x[component]) / (w1 + w2));
106 // mark if it has greater diameter
107 //Log() << Verbose(1) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;
108 GreatestDiameter[component] = (GreatestDiameter[component]
109 > tmp) ? GreatestDiameter[component] : tmp;
110 } //else
111 //Log() << Verbose(1) << "Saw no sign flip, probably top or bottom node." << endl;
112 }
113 }
114 }
115 Log() << Verbose(0) << "RESULT: The biggest diameters are "
116 << GreatestDiameter[0] << " and " << GreatestDiameter[1] << " and "
117 << GreatestDiameter[2] << " " << (IsAngstroem ? "angstrom"
118 : "atomiclength") << "." << endl;
119
120 // free reference lists
121 if (BoundaryFreeFlag)
122 delete[] (BoundaryPoints);
123
124 return GreatestDiameter;
125}
126;
127
128
129/** Determines the boundary points of a cluster.
130 * Does a projection per axis onto the orthogonal plane, transforms into spherical coordinates, sorts them by the angle
131 * and looks at triples: if the middle has less a distance than the allowed maximum height of the triangle formed by the plane's
132 * center and first and last point in the triple, it is thrown out.
133 * \param *out output stream for debugging
134 * \param *mol molecule structure representing the cluster
135 * \param *&TesselStruct pointer to Tesselation structure
136 */
137Boundaries *GetBoundaryPoints(const molecule *mol, Tesselation *&TesselStruct)
138{
139 Info FunctionInfo(__func__);
140 atom *Walker = NULL;
141 PointMap PointsOnBoundary;
142 LineMap LinesOnBoundary;
143 TriangleMap TrianglesOnBoundary;
144 Vector *MolCenter = mol->DetermineCenterOfAll();
145 Vector helper;
146 BoundariesTestPair BoundaryTestPair;
147 Vector AxisVector;
148 Vector AngleReferenceVector;
149 Vector AngleReferenceNormalVector;
150 Vector ProjectedVector;
151 Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr)
152 double angle = 0.;
153
154 // 3a. Go through every axis
155 for (int axis = 0; axis < NDIM; axis++) {
156 AxisVector.Zero();
157 AngleReferenceVector.Zero();
158 AngleReferenceNormalVector.Zero();
159 AxisVector.x[axis] = 1.;
160 AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;
161 AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
162
163 Log() << Verbose(1) << "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << endl;
164
165 // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
166 Walker = mol->start;
167 while (Walker->next != mol->end) {
168 Walker = Walker->next;
169 ProjectedVector.CopyVector(&Walker->x);
170 ProjectedVector.SubtractVector(MolCenter);
171 ProjectedVector.ProjectOntoPlane(&AxisVector);
172
173 // correct for negative side
174 const double radius = ProjectedVector.NormSquared();
175 if (fabs(radius) > MYEPSILON)
176 angle = ProjectedVector.Angle(&AngleReferenceVector);
177 else
178 angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
179
180 //Log() << Verbose(1) << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
181 if (ProjectedVector.ScalarProduct(&AngleReferenceNormalVector) > 0) {
182 angle = 2. * M_PI - angle;
183 }
184 Log() << Verbose(1) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl;
185 BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, Walker)));
186 if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity
187 Log() << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl;
188 Log() << Verbose(2) << "Present vector: " << *BoundaryTestPair.first->second.second << endl;
189 Log() << Verbose(2) << "New vector: " << *Walker << endl;
190 const double ProjectedVectorNorm = ProjectedVector.NormSquared();
191 if ((ProjectedVectorNorm - BoundaryTestPair.first->second.first) > MYEPSILON) {
192 BoundaryTestPair.first->second.first = ProjectedVectorNorm;
193 BoundaryTestPair.first->second.second = Walker;
194 Log() << Verbose(2) << "Keeping new vector due to larger projected distance " << ProjectedVectorNorm << "." << endl;
195 } else if (fabs(ProjectedVectorNorm - BoundaryTestPair.first->second.first) < MYEPSILON) {
196 helper.CopyVector(&Walker->x);
197 helper.SubtractVector(MolCenter);
198 const double oldhelperNorm = helper.NormSquared();
199 helper.CopyVector(&BoundaryTestPair.first->second.second->x);
200 helper.SubtractVector(MolCenter);
201 if (helper.NormSquared() < oldhelperNorm) {
202 BoundaryTestPair.first->second.second = Walker;
203 Log() << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl;
204 } else {
205 Log() << Verbose(2) << "Keeping present vector due to larger distance to molecule center " << oldhelperNorm << "." << endl;
206 }
207 } else {
208 Log() << Verbose(2) << "Keeping present vector due to larger projected distance " << ProjectedVectorNorm << "." << endl;
209 }
210 }
211 }
212 // printing all inserted for debugging
213 // {
214 // Log() << Verbose(1) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
215 // int i=0;
216 // for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
217 // if (runner != BoundaryPoints[axis].begin())
218 // Log() << Verbose(0) << ", " << i << ": " << *runner->second.second;
219 // else
220 // Log() << Verbose(0) << i << ": " << *runner->second.second;
221 // i++;
222 // }
223 // Log() << Verbose(0) << endl;
224 // }
225 // 3c. throw out points whose distance is less than the mean of left and right neighbours
226 bool flag = false;
227 Log() << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl;
228 do { // do as long as we still throw one out per round
229 flag = false;
230 Boundaries::iterator left = BoundaryPoints[axis].end();
231 Boundaries::iterator right = BoundaryPoints[axis].end();
232 for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
233 // set neighbours correctly
234 if (runner == BoundaryPoints[axis].begin()) {
235 left = BoundaryPoints[axis].end();
236 } else {
237 left = runner;
238 }
239 left--;
240 right = runner;
241 right++;
242 if (right == BoundaryPoints[axis].end()) {
243 right = BoundaryPoints[axis].begin();
244 }
245 // check distance
246
247 // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
248 {
249 Vector SideA, SideB, SideC, SideH;
250 SideA.CopyVector(&left->second.second->x);
251 SideA.SubtractVector(MolCenter);
252 SideA.ProjectOntoPlane(&AxisVector);
253 // Log() << Verbose(1) << "SideA: " << SideA << endl;
254
255 SideB.CopyVector(&right->second.second->x);
256 SideB.SubtractVector(MolCenter);
257 SideB.ProjectOntoPlane(&AxisVector);
258 // Log() << Verbose(1) << "SideB: " << SideB << endl;
259
260 SideC.CopyVector(&left->second.second->x);
261 SideC.SubtractVector(&right->second.second->x);
262 SideC.ProjectOntoPlane(&AxisVector);
263 // Log() << Verbose(1) << "SideC: " << SideC << endl;
264
265 SideH.CopyVector(&runner->second.second->x);
266 SideH.SubtractVector(MolCenter);
267 SideH.ProjectOntoPlane(&AxisVector);
268 // Log() << Verbose(1) << "SideH: " << SideH << endl;
269
270 // calculate each length
271 const double a = SideA.Norm();
272 //const double b = SideB.Norm();
273 //const double c = SideC.Norm();
274 const double h = SideH.Norm();
275 // calculate the angles
276 const double alpha = SideA.Angle(&SideH);
277 const double beta = SideA.Angle(&SideC);
278 const double gamma = SideB.Angle(&SideH);
279 const double delta = SideC.Angle(&SideH);
280 const double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
281 //Log() << Verbose(1) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
282 Log() << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
283 if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance)) < MYEPSILON) && ((h - MinDistance)) < -MYEPSILON) {
284 // throw out point
285 Log() << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
286 BoundaryPoints[axis].erase(runner);
287 flag = true;
288 }
289 }
290 }
291 } while (flag);
292 }
293 delete(MolCenter);
294 return BoundaryPoints;
295};
296
297/** Tesselates the convex boundary by finding all boundary points.
298 * \param *out output stream for debugging
299 * \param *mol molecule structure with Atom's and Bond's.
300 * \param *TesselStruct Tesselation filled with points, lines and triangles on boundary on return
301 * \param *LCList atoms in LinkedCell list
302 * \param *filename filename prefix for output of vertex data
303 * \return *TesselStruct is filled with convex boundary and tesselation is stored under \a *filename.
304 */
305void FindConvexBorder(const molecule* mol, Tesselation *&TesselStruct, const LinkedCell *LCList, const char *filename)
306{
307 Info FunctionInfo(__func__);
308 bool BoundaryFreeFlag = false;
309 Boundaries *BoundaryPoints = NULL;
310
311 if (TesselStruct != NULL) // free if allocated
312 delete(TesselStruct);
313 TesselStruct = new class Tesselation;
314
315 // 1. Find all points on the boundary
316 if (BoundaryPoints == NULL) {
317 BoundaryFreeFlag = true;
318 BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
319 } else {
320 Log() << Verbose(0) << "Using given boundary points set." << endl;
321 }
322
323// printing all inserted for debugging
324 for (int axis=0; axis < NDIM; axis++)
325 {
326 Log() << Verbose(1) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
327 int i=0;
328 for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
329 if (runner != BoundaryPoints[axis].begin())
330 Log() << Verbose(0) << ", " << i << ": " << *runner->second.second;
331 else
332 Log() << Verbose(0) << i << ": " << *runner->second.second;
333 i++;
334 }
335 Log() << Verbose(0) << endl;
336 }
337
338 // 2. fill the boundary point list
339 for (int axis = 0; axis < NDIM; axis++)
340 for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++)
341 if (!TesselStruct->AddBoundaryPoint(runner->second.second, 0))
342 eLog() << Verbose(2) << "Point " << *(runner->second.second) << " is already present!" << endl;
343
344 Log() << Verbose(0) << "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl;
345 // now we have the whole set of edge points in the BoundaryList
346
347 // listing for debugging
348 // Log() << Verbose(1) << "Listing PointsOnBoundary:";
349 // for(PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
350 // Log() << Verbose(0) << " " << *runner->second;
351 // }
352 // Log() << Verbose(0) << endl;
353
354 // 3a. guess starting triangle
355 TesselStruct->GuessStartingTriangle();
356
357 // 3b. go through all lines, that are not yet part of two triangles (only of one so far)
358 TesselStruct->TesselateOnBoundary(mol);
359
360 // 3c. check whether all atoms lay inside the boundary, if not, add to boundary points, segment triangle into three with the new point
361 if (!TesselStruct->InsertStraddlingPoints(mol, LCList))
362 eLog() << Verbose(1) << "Insertion of straddling points failed!" << endl;
363
364 Log() << Verbose(0) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " intermediate triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;
365
366 // 4. Store triangles in tecplot file
367 if (filename != NULL) {
368 if (DoTecplotOutput) {
369 string OutputName(filename);
370 OutputName.append("_intermed");
371 OutputName.append(TecplotSuffix);
372 ofstream *tecplot = new ofstream(OutputName.c_str());
373 WriteTecplotFile(tecplot, TesselStruct, mol, 0);
374 tecplot->close();
375 delete(tecplot);
376 }
377 if (DoRaster3DOutput) {
378 string OutputName(filename);
379 OutputName.append("_intermed");
380 OutputName.append(Raster3DSuffix);
381 ofstream *rasterplot = new ofstream(OutputName.c_str());
382 WriteRaster3dFile(rasterplot, TesselStruct, mol);
383 rasterplot->close();
384 delete(rasterplot);
385 }
386 }
387
388 // 3d. check all baselines whether the peaks of the two adjacent triangles with respect to center of baseline are convex, if not, make the baseline between the two peaks and baseline endpoints become the new peaks
389 bool AllConvex = true;
390 class BoundaryLineSet *line = NULL;
391 do {
392 AllConvex = true;
393 for (LineMap::iterator LineRunner = TesselStruct->LinesOnBoundary.begin(); LineRunner != TesselStruct->LinesOnBoundary.end(); LineRunner++) {
394 line = LineRunner->second;
395 Log() << Verbose(1) << "INFO: Current line is " << *line << "." << endl;
396 if (!line->CheckConvexityCriterion()) {
397 Log() << Verbose(1) << "... line " << *line << " is concave, flipping it." << endl;
398
399 // flip the line
400 if (TesselStruct->PickFarthestofTwoBaselines(line) == 0.)
401 eLog() << Verbose(1) << "Correction of concave baselines failed!" << endl;
402 else {
403 TesselStruct->FlipBaseline(line);
404 Log() << Verbose(1) << "INFO: Correction of concave baselines worked." << endl;
405 }
406 }
407 }
408 } while (!AllConvex);
409
410 // 3e. we need another correction here, for TesselPoints that are below the surface (i.e. have an odd number of concave triangles surrounding it)
411// if (!TesselStruct->CorrectConcaveTesselPoints(out))
412// Log() << Verbose(1) << "Correction of concave tesselpoints failed!" << endl;
413
414 Log() << Verbose(0) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;
415
416 // 4. Store triangles in tecplot file
417 if (filename != NULL) {
418 if (DoTecplotOutput) {
419 string OutputName(filename);
420 OutputName.append(TecplotSuffix);
421 ofstream *tecplot = new ofstream(OutputName.c_str());
422 WriteTecplotFile(tecplot, TesselStruct, mol, 0);
423 tecplot->close();
424 delete(tecplot);
425 }
426 if (DoRaster3DOutput) {
427 string OutputName(filename);
428 OutputName.append(Raster3DSuffix);
429 ofstream *rasterplot = new ofstream(OutputName.c_str());
430 WriteRaster3dFile(rasterplot, TesselStruct, mol);
431 rasterplot->close();
432 delete(rasterplot);
433 }
434 }
435
436
437 // free reference lists
438 if (BoundaryFreeFlag)
439 delete[] (BoundaryPoints);
440};
441
442/** For testing removes one boundary point after another to check for leaks.
443 * \param *out output stream for debugging
444 * \param *TesselStruct Tesselation containing envelope with boundary points
445 * \param *mol molecule
446 * \param *filename name of file
447 * \return true - all removed, false - something went wrong
448 */
449bool RemoveAllBoundaryPoints(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename)
450{
451 Info FunctionInfo(__func__);
452 int i=0;
453 char number[MAXSTRINGSIZE];
454
455 if ((TesselStruct == NULL) || (TesselStruct->PointsOnBoundary.empty())) {
456 eLog() << Verbose(1) << "TesselStruct is empty." << endl;
457 return false;
458 }
459
460 PointMap::iterator PointRunner;
461 while (!TesselStruct->PointsOnBoundary.empty()) {
462 Log() << Verbose(1) << "Remaining points are: ";
463 for (PointMap::iterator PointSprinter = TesselStruct->PointsOnBoundary.begin(); PointSprinter != TesselStruct->PointsOnBoundary.end(); PointSprinter++)
464 Log() << Verbose(0) << *(PointSprinter->second) << "\t";
465 Log() << Verbose(0) << endl;
466
467 PointRunner = TesselStruct->PointsOnBoundary.begin();
468 // remove point
469 TesselStruct->RemovePointFromTesselatedSurface(PointRunner->second);
470
471 // store envelope
472 sprintf(number, "-%04d", i++);
473 StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, number);
474 }
475
476 return true;
477};
478
479/** Creates a convex envelope from a given non-convex one.
480 * -# First step, remove concave spots, i.e. singular "dents"
481 * -# We go through all PointsOnBoundary.
482 * -# We CheckConvexityCriterion() for all its lines.
483 * -# If all its lines are concave, it cannot be on the convex envelope.
484 * -# Hence, we remove it and re-create all its triangles from its getCircleOfConnectedPoints()
485 * -# We calculate the additional volume.
486 * -# We go over all lines until none yields a concavity anymore.
487 * -# Second step, remove concave lines, i.e. line-shape "dents"
488 * -# We go through all LinesOnBoundary
489 * -# We CheckConvexityCriterion()
490 * -# If it returns concave, we flip the line in this quadrupel of points (abusing the degeneracy of the tesselation)
491 * -# We CheckConvexityCriterion(),
492 * -# if it's concave, we continue
493 * -# if not, we mark an error and stop
494 * Note: This routine - for free - calculates the difference in volume between convex and
495 * non-convex envelope, as the former is easy to calculate - VolumeOfConvexEnvelope() - it
496 * can be used to compute volumes of arbitrary shapes.
497 * \param *out output stream for debugging
498 * \param *TesselStruct non-convex envelope, is changed in return!
499 * \param *mol molecule
500 * \param *filename name of file
501 * \return volume difference between the non- and the created convex envelope
502 */
503double ConvexizeNonconvexEnvelope(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename)
504{
505 Info FunctionInfo(__func__);
506 double volume = 0;
507 class BoundaryPointSet *point = NULL;
508 class BoundaryLineSet *line = NULL;
509 bool Concavity = false;
510 char dummy[MAXSTRINGSIZE];
511 PointMap::iterator PointRunner;
512 PointMap::iterator PointAdvance;
513 LineMap::iterator LineRunner;
514 LineMap::iterator LineAdvance;
515 TriangleMap::iterator TriangleRunner;
516 TriangleMap::iterator TriangleAdvance;
517 int run = 0;
518
519 // check whether there is something to work on
520 if (TesselStruct == NULL) {
521 eLog() << Verbose(1) << "TesselStruct is empty!" << endl;
522 return volume;
523 }
524
525 // First step: RemovePointFromTesselatedSurface
526 do {
527 Concavity = false;
528 sprintf(dummy, "-first-%d", run);
529 //CalculateConcavityPerBoundaryPoint(TesselStruct);
530 StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);
531
532 PointRunner = TesselStruct->PointsOnBoundary.begin();
533 PointAdvance = PointRunner; // we need an advanced point, as the PointRunner might get removed
534 while (PointRunner != TesselStruct->PointsOnBoundary.end()) {
535 PointAdvance++;
536 point = PointRunner->second;
537 Log() << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
538 for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
539 line = LineRunner->second;
540 Log() << Verbose(1) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
541 if (!line->CheckConvexityCriterion()) {
542 // remove the point if needed
543 Log() << Verbose(1) << "... point " << *point << " cannot be on convex envelope." << endl;
544 volume += TesselStruct->RemovePointFromTesselatedSurface(point);
545 sprintf(dummy, "-first-%d", ++run);
546 StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);
547 Concavity = true;
548 break;
549 }
550 }
551 PointRunner = PointAdvance;
552 }
553
554 sprintf(dummy, "-second-%d", run);
555 //CalculateConcavityPerBoundaryPoint(TesselStruct);
556 StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);
557
558 // second step: PickFarthestofTwoBaselines
559 LineRunner = TesselStruct->LinesOnBoundary.begin();
560 LineAdvance = LineRunner; // we need an advanced line, as the LineRunner might get removed
561 while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
562 LineAdvance++;
563 line = LineRunner->second;
564 Log() << Verbose(1) << "INFO: Picking farthest baseline for line is " << *line << "." << endl;
565 // take highest of both lines
566 if (TesselStruct->IsConvexRectangle(line) == NULL) {
567 const double tmp = TesselStruct->PickFarthestofTwoBaselines(line);
568 volume += tmp;
569 if (tmp != 0.) {
570 TesselStruct->FlipBaseline(line);
571 Concavity = true;
572 }
573 }
574 LineRunner = LineAdvance;
575 }
576 run++;
577 } while (Concavity);
578 //CalculateConcavityPerBoundaryPoint(TesselStruct);
579 //StoreTrianglesinFile(mol, filename, "-third");
580
581 // third step: IsConvexRectangle
582// LineRunner = TesselStruct->LinesOnBoundary.begin();
583// LineAdvance = LineRunner; // we need an advanced line, as the LineRunner might get removed
584// while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
585// LineAdvance++;
586// line = LineRunner->second;
587// Log() << Verbose(1) << "INFO: Current line is " << *line << "." << endl;
588// //if (LineAdvance != TesselStruct->LinesOnBoundary.end())
589// //Log() << Verbose(1) << "INFO: Next line will be " << *(LineAdvance->second) << "." << endl;
590// if (!line->CheckConvexityCriterion(out)) {
591// Log() << Verbose(1) << "... line " << *line << " is concave, flipping it." << endl;
592//
593// // take highest of both lines
594// point = TesselStruct->IsConvexRectangle(line);
595// if (point != NULL)
596// volume += TesselStruct->RemovePointFromTesselatedSurface(point);
597// }
598// LineRunner = LineAdvance;
599// }
600
601 CalculateConcavityPerBoundaryPoint(TesselStruct);
602 StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, "");
603
604 // end
605 Log() << Verbose(0) << "Volume is " << volume << "." << endl;
606 return volume;
607};
608
609
610/** Determines the volume of a cluster.
611 * Determines first the convex envelope, then tesselates it and calculates its volume.
612 * \param *out output stream for debugging
613 * \param *TesselStruct Tesselation filled with points, lines and triangles on boundary on return
614 * \param *configuration needed for path to store convex envelope file
615 * \return determined volume of the cluster in cubed config:GetIsAngstroem()
616 */
617double VolumeOfConvexEnvelope(class Tesselation *TesselStruct, class config *configuration)
618{
619 Info FunctionInfo(__func__);
620 bool IsAngstroem = configuration->GetIsAngstroem();
621 double volume = 0.;
622 Vector x;
623 Vector y;
624
625 // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
626 for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++)
627 { // go through every triangle, calculate volume of its pyramid with CoG as peak
628 x.CopyVector(runner->second->endpoints[0]->node->node);
629 x.SubtractVector(runner->second->endpoints[1]->node->node);
630 y.CopyVector(runner->second->endpoints[0]->node->node);
631 y.SubtractVector(runner->second->endpoints[2]->node->node);
632 const double a = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(runner->second->endpoints[1]->node->node));
633 const double b = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(runner->second->endpoints[2]->node->node));
634 const double c = sqrt(runner->second->endpoints[2]->node->node->DistanceSquared(runner->second->endpoints[1]->node->node));
635 const double G = sqrt(((a + b + c) * (a + b + c) - 2 * (a * a + b * b + c * c)) / 16.); // area of tesselated triangle
636 x.MakeNormalVector(runner->second->endpoints[0]->node->node, runner->second->endpoints[1]->node->node, runner->second->endpoints[2]->node->node);
637 x.Scale(runner->second->endpoints[1]->node->node->ScalarProduct(&x));
638 const double h = x.Norm(); // distance of CoG to triangle
639 const double PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
640 Log() << Verbose(1) << "Area of triangle is " << setprecision(10) << G << " "
641 << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is "
642 << h << " and the volume is " << PyramidVolume << " "
643 << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
644 volume += PyramidVolume;
645 }
646 Log() << Verbose(0) << "RESULT: The summed volume is " << setprecision(6)
647 << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3."
648 << endl;
649
650 return volume;
651};
652
653/** Stores triangles to file.
654 * \param *out output stream for debugging
655 * \param *mol molecule with atoms and bonds
656 * \param *&TesselStruct Tesselation with boundary triangles
657 * \param *filename prefix of filename
658 * \param *extraSuffix intermediate suffix
659 */
660void StoreTrianglesinFile(const molecule * const mol, const Tesselation *&TesselStruct, const char *filename, const char *extraSuffix)
661{
662 Info FunctionInfo(__func__);
663 // 4. Store triangles in tecplot file
664 if (filename != NULL) {
665 if (DoTecplotOutput) {
666 string OutputName(filename);
667 OutputName.append(extraSuffix);
668 OutputName.append(TecplotSuffix);
669 ofstream *tecplot = new ofstream(OutputName.c_str());
670 WriteTecplotFile(tecplot, TesselStruct, mol, -1);
671 tecplot->close();
672 delete(tecplot);
673 }
674 if (DoRaster3DOutput) {
675 string OutputName(filename);
676 OutputName.append(extraSuffix);
677 OutputName.append(Raster3DSuffix);
678 ofstream *rasterplot = new ofstream(OutputName.c_str());
679 WriteRaster3dFile(rasterplot, TesselStruct, mol);
680 rasterplot->close();
681 delete(rasterplot);
682 }
683 }
684};
685
686/** Creates multiples of the by \a *mol given cluster and suspends them in water with a given final density.
687 * We get cluster volume by VolumeOfConvexEnvelope() and its diameters by GetDiametersOfCluster()
688 * \param *out output stream for debugging
689 * \param *configuration needed for path to store convex envelope file
690 * \param *mol molecule structure representing the cluster
691 * \param *&TesselStruct Tesselation structure with triangles on return
692 * \param ClusterVolume guesstimated cluster volume, if equal 0 we used VolumeOfConvexEnvelope() instead.
693 * \param celldensity desired average density in final cell
694 */
695void PrepareClustersinWater(config *configuration, molecule *mol, double ClusterVolume, double celldensity)
696{
697 Info FunctionInfo(__func__);
698 bool IsAngstroem = true;
699 double *GreatestDiameter = NULL;
700 Boundaries *BoundaryPoints = NULL;
701 class Tesselation *TesselStruct = NULL;
702 Vector BoxLengths;
703 int repetition[NDIM] = { 1, 1, 1 };
704 int TotalNoClusters = 1;
705 atom *Walker = NULL;
706 double totalmass = 0.;
707 double clustervolume = 0.;
708 double cellvolume = 0.;
709
710 // transform to PAS
711 mol->PrincipalAxisSystem(true);
712
713 IsAngstroem = configuration->GetIsAngstroem();
714 GreatestDiameter = GetDiametersOfCluster(BoundaryPoints, mol, TesselStruct, IsAngstroem);
715 BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
716 LinkedCell LCList(mol, 10.);
717 FindConvexBorder(mol, TesselStruct, &LCList, NULL);
718
719 // some preparations beforehand
720 if (ClusterVolume == 0)
721 clustervolume = VolumeOfConvexEnvelope(TesselStruct, configuration);
722 else
723 clustervolume = ClusterVolume;
724
725 for (int i = 0; i < NDIM; i++)
726 TotalNoClusters *= repetition[i];
727
728 // sum up the atomic masses
729 Walker = mol->start;
730 while (Walker->next != mol->end) {
731 Walker = Walker->next;
732 totalmass += Walker->type->mass;
733 }
734 Log() << Verbose(0) << "RESULT: The summed mass is " << setprecision(10) << totalmass << " atomicmassunit." << endl;
735 Log() << Verbose(0) << "RESULT: The average density is " << setprecision(10) << totalmass / clustervolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
736
737 // solve cubic polynomial
738 Log() << Verbose(1) << "Solving equidistant suspension in water problem ..." << endl;
739 if (IsAngstroem)
740 cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_A - (totalmass / clustervolume)) / (celldensity - 1);
741 else
742 cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_a0 - (totalmass / clustervolume)) / (celldensity - 1);
743 Log() << Verbose(1) << "Cellvolume needed for a density of " << celldensity << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
744
745 double minimumvolume = TotalNoClusters * (GreatestDiameter[0] * GreatestDiameter[1] * GreatestDiameter[2]);
746 Log() << Verbose(1) << "Minimum volume of the convex envelope contained in a rectangular box is " << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
747 if (minimumvolume > cellvolume) {
748 eLog() << Verbose(1) << "the containing box already has a greater volume than the envisaged cell volume!" << endl;
749 Log() << Verbose(0) << "Setting Box dimensions to minimum possible, the greatest diameters." << endl;
750 for (int i = 0; i < NDIM; i++)
751 BoxLengths.x[i] = GreatestDiameter[i];
752 mol->CenterEdge(&BoxLengths);
753 } else {
754 BoxLengths.x[0] = (repetition[0] * GreatestDiameter[0] + repetition[1] * GreatestDiameter[1] + repetition[2] * GreatestDiameter[2]);
755 BoxLengths.x[1] = (repetition[0] * repetition[1] * GreatestDiameter[0] * GreatestDiameter[1] + repetition[0] * repetition[2] * GreatestDiameter[0] * GreatestDiameter[2] + repetition[1] * repetition[2] * GreatestDiameter[1] * GreatestDiameter[2]);
756 BoxLengths.x[2] = minimumvolume - cellvolume;
757 double x0 = 0.;
758 double x1 = 0.;
759 double x2 = 0.;
760 if (gsl_poly_solve_cubic(BoxLengths.x[0], BoxLengths.x[1], BoxLengths.x[2], &x0, &x1, &x2) == 1) // either 1 or 3 on return
761 Log() << Verbose(0) << "RESULT: The resulting spacing is: " << x0 << " ." << endl;
762 else {
763 Log() << Verbose(0) << "RESULT: The resulting spacings are: " << x0 << " and " << x1 << " and " << x2 << " ." << endl;
764 x0 = x2; // sorted in ascending order
765 }
766
767 cellvolume = 1.;
768 for (int i = 0; i < NDIM; i++) {
769 BoxLengths.x[i] = repetition[i] * (x0 + GreatestDiameter[i]);
770 cellvolume *= BoxLengths.x[i];
771 }
772
773 // set new box dimensions
774 Log() << Verbose(0) << "Translating to box with these boundaries." << endl;
775 mol->SetBoxDimension(&BoxLengths);
776 mol->CenterInBox();
777 }
778 // update Box of atoms by boundary
779 mol->SetBoxDimension(&BoxLengths);
780 Log() << Verbose(0) << "RESULT: The resulting cell dimensions are: " << BoxLengths.x[0] << " and " << BoxLengths.x[1] << " and " << BoxLengths.x[2] << " with total volume of " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
781};
782
783
784/** Fills the empty space of the simulation box with water/
785 * \param *out output stream for debugging
786 * \param *List list of molecules already present in box
787 * \param *TesselStruct contains tesselated surface
788 * \param *filler molecule which the box is to be filled with
789 * \param configuration contains box dimensions
790 * \param distance[NDIM] distance between filling molecules in each direction
791 * \param RandAtomDisplacement maximum distance for random displacement per atom
792 * \param RandMolDisplacement maximum distance for random displacement per filler molecule
793 * \param DoRandomRotation true - do random rotiations, false - don't
794 * \return *mol pointer to new molecule with filled atoms
795 */
796molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, double distance[NDIM], double RandomAtomDisplacement, double RandomMolDisplacement, bool DoRandomRotation)
797{
798 Info FunctionInfo(__func__);
799 molecule *Filling = new molecule(filler->elemente);
800 Vector CurrentPosition;
801 int N[NDIM];
802 int n[NDIM];
803 double *M = ReturnFullMatrixforSymmetric(filler->cell_size);
804 double Rotations[NDIM*NDIM];
805 Vector AtomTranslations;
806 Vector FillerTranslations;
807 Vector FillerDistance;
808 double FillIt = false;
809 atom *Walker = NULL;
810 bond *Binder = NULL;
811 int i = 0;
812 LinkedCell *LCList[List->ListOfMolecules.size()];
813 double phi[NDIM];
814 class Tesselation *TesselStruct[List->ListOfMolecules.size()];
815
816 i=0;
817 for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
818 Log() << Verbose(1) << "Pre-creating linked cell lists for molecule " << *ListRunner << "." << endl;
819 LCList[i] = new LinkedCell((*ListRunner), 5.); // get linked cell list
820 if (TesselStruct[i] == NULL) {
821 Log() << Verbose(1) << "Pre-creating tesselation for molecule " << *ListRunner << "." << endl;
822 FindNonConvexBorder((*ListRunner), TesselStruct[i], (const LinkedCell *&)LCList[i], 5., NULL);
823 }
824 i++;
825 }
826
827 // Center filler at origin
828 filler->CenterOrigin();
829 filler->Center.Zero();
830
831 filler->CountAtoms();
832 atom * CopyAtoms[filler->AtomCount];
833
834 // calculate filler grid in [0,1]^3
835 FillerDistance.Init(distance[0], distance[1], distance[2]);
836 FillerDistance.InverseMatrixMultiplication(M);
837 Log() << Verbose(1) << "INFO: Grid steps are ";
838 for(int i=0;i<NDIM;i++) {
839 N[i] = (int) ceil(1./FillerDistance.x[i]);
840 Log() << Verbose(1) << N[i];
841 if (i != NDIM-1)
842 Log() << Verbose(1)<< ", ";
843 else
844 Log() << Verbose(1) << "." << endl;
845 }
846
847 // go over [0,1]^3 filler grid
848 for (n[0] = 0; n[0] < N[0]; n[0]++)
849 for (n[1] = 0; n[1] < N[1]; n[1]++)
850 for (n[2] = 0; n[2] < N[2]; n[2]++) {
851 // calculate position of current grid vector in untransformed box
852 CurrentPosition.Init((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
853 CurrentPosition.MatrixMultiplication(M);
854 Log() << Verbose(2) << "INFO: Current Position is " << CurrentPosition << "." << endl;
855 // Check whether point is in- or outside
856 FillIt = true;
857 i=0;
858 for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
859 // get linked cell list
860 if (TesselStruct[i] == NULL) {
861 eLog() << Verbose(1) << "TesselStruct of " << (*ListRunner) << " is NULL. Didn't we pre-create it?" << endl;
862 FillIt = false;
863 } else {
864 FillIt = FillIt && (!TesselStruct[i]->IsInnerPoint(CurrentPosition, LCList[i]));
865 i++;
866 }
867 }
868
869 if (FillIt) {
870 // fill in Filler
871 Log() << Verbose(2) << "Space at " << CurrentPosition << " is unoccupied by any molecule, filling in." << endl;
872
873 // create molecule random translation vector ...
874 for (int i=0;i<NDIM;i++)
875 FillerTranslations.x[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
876 Log() << Verbose(2) << "INFO: Translating this filler by " << FillerTranslations << "." << endl;
877
878 // go through all atoms
879 Walker = filler->start;
880 while (Walker->next != filler->end) {
881 Walker = Walker->next;
882 // copy atom ...
883 CopyAtoms[Walker->nr] = new atom(Walker);
884
885 // create atomic random translation vector ...
886 for (int i=0;i<NDIM;i++)
887 AtomTranslations.x[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.);
888
889 // ... and rotation matrix
890 if (DoRandomRotation) {
891 for (int i=0;i<NDIM;i++) {
892 phi[i] = rand()/(RAND_MAX/(2.*M_PI));
893 }
894
895 Rotations[0] = cos(phi[0]) *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]));
896 Rotations[3] = sin(phi[0]) *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]));
897 Rotations[6] = cos(phi[1])*sin(phi[2]) ;
898 Rotations[1] = - sin(phi[0])*cos(phi[1]) ;
899 Rotations[4] = cos(phi[0])*cos(phi[1]) ;
900 Rotations[7] = sin(phi[1]) ;
901 Rotations[3] = - cos(phi[0]) *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]));
902 Rotations[5] = - sin(phi[0]) *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]));
903 Rotations[8] = cos(phi[1])*cos(phi[2]) ;
904 }
905
906 // ... and put at new position
907 if (DoRandomRotation)
908 CopyAtoms[Walker->nr]->x.MatrixMultiplication(Rotations);
909 CopyAtoms[Walker->nr]->x.AddVector(&AtomTranslations);
910 CopyAtoms[Walker->nr]->x.AddVector(&FillerTranslations);
911 CopyAtoms[Walker->nr]->x.AddVector(&CurrentPosition);
912
913 // insert into Filling
914
915 // FIXME: gives completely different results if CopyAtoms[..] used instead of Walker, why???
916 Log() << Verbose(4) << "Filling atom " << *Walker << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[Walker->nr]->x) << "." << endl;
917 Filling->AddAtom(CopyAtoms[Walker->nr]);
918 }
919
920 // go through all bonds and add as well
921 Binder = filler->first;
922 while(Binder->next != filler->last) {
923 Binder = Binder->next;
924 Log() << Verbose(3) << "Adding Bond between " << *CopyAtoms[Binder->leftatom->nr] << " and " << *CopyAtoms[Binder->rightatom->nr]<< "." << endl;
925 Filling->AddBond(CopyAtoms[Binder->leftatom->nr], CopyAtoms[Binder->rightatom->nr], Binder->BondDegree);
926 }
927 } else {
928 // leave empty
929 Log() << Verbose(2) << "Space at " << CurrentPosition << " is occupied." << endl;
930 }
931 }
932 Free(&M);
933 for (size_t i=0;i<List->ListOfMolecules.size();i++) {
934 delete(LCList[i]);
935 delete(TesselStruct[i]);
936 }
937 return Filling;
938};
939
940
941/** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule.
942 * \param *out output stream for debugging
943 * \param *mol molecule structure with Atom's and Bond's
944 * \param *&TesselStruct Tesselation filled with points, lines and triangles on boundary on return
945 * \param *&LCList atoms in LinkedCell list
946 * \param RADIUS radius of the virtual sphere
947 * \param *filename filename prefix for output of vertex data
948 * \return true - tesselation successful, false - tesselation failed
949 */
950bool FindNonConvexBorder(const molecule* const mol, Tesselation *&TesselStruct, const LinkedCell *&LCList, const double RADIUS, const char *filename = NULL)
951{
952 Info FunctionInfo(__func__);
953 bool freeLC = false;
954 bool status = false;
955 CandidateForTesselation *baseline;
956 LineMap::iterator testline;
957 bool OneLoopWithoutSuccessFlag = true; // marks whether we went once through all baselines without finding any without two triangles
958 bool TesselationFailFlag = false;
959 BoundaryTriangleSet *T = NULL;
960
961 if (TesselStruct == NULL) {
962 Log() << Verbose(1) << "Allocating Tesselation struct ..." << endl;
963 TesselStruct= new Tesselation;
964 } else {
965 delete(TesselStruct);
966 Log() << Verbose(1) << "Re-Allocating Tesselation struct ..." << endl;
967 TesselStruct = new Tesselation;
968 }
969
970 // initialise Linked Cell
971 if (LCList == NULL) {
972 LCList = new LinkedCell(mol, 2.*RADIUS);
973 freeLC = true;
974 }
975
976 // 1. get starting triangle
977 TesselStruct->FindStartingTriangle(RADIUS, LCList);
978
979 // 2. expand from there
980 while ((!TesselStruct->OpenLines.empty()) && (OneLoopWithoutSuccessFlag)) {
981 // 2a. fill all new OpenLines
982 Log() << Verbose(1) << "There are " << TesselStruct->OpenLines.size() << " open lines to scan for candidates:" << endl;
983 for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++)
984 Log() << Verbose(2) << *(Runner->second) << endl;
985
986 for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
987 baseline = Runner->second;
988 if (baseline->pointlist.empty()) {
989 T = (((baseline->BaseLine->triangles.begin()))->second);
990 Log() << Verbose(1) << "Finding best candidate for open line " << *baseline->BaseLine << " of triangle " << *T << endl;
991 TesselationFailFlag = TesselStruct->FindNextSuitableTriangle(*baseline, *T, RADIUS, LCList); //the line is there, so there is a triangle, but only one.
992 }
993 }
994
995 // 2b. search for smallest ShortestAngle among all candidates
996 double ShortestAngle = 4.*M_PI;
997 Log() << Verbose(1) << "There are " << TesselStruct->OpenLines.size() << " open lines to scan for the best candidates:" << endl;
998 for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++)
999 Log() << Verbose(2) << *(Runner->second) << endl;
1000
1001 for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
1002 if (Runner->second->ShortestAngle < ShortestAngle) {
1003 baseline = Runner->second;
1004 ShortestAngle = baseline->ShortestAngle;
1005 //Log() << Verbose(1) << "New best candidate is " << *baseline->BaseLine << " with point " << *baseline->point << " and angle " << baseline->ShortestAngle << endl;
1006 }
1007 }
1008 if ((ShortestAngle == 4.*M_PI) || (baseline->pointlist.empty()))
1009 OneLoopWithoutSuccessFlag = false;
1010 else {
1011 TesselStruct->AddCandidateTriangle(*baseline);
1012 }
1013
1014 // write temporary envelope
1015 if (filename != NULL) {
1016 if ((DoSingleStepOutput && ((TesselStruct->TrianglesOnBoundary.size() % SingleStepWidth == 0)))) { // if we have a new triangle and want to output each new triangle configuration
1017 TesselStruct->Output(filename, mol);
1018 }
1019 }
1020 }
1021// // check envelope for consistency
1022// status = CheckListOfBaselines(TesselStruct);
1023//
1024// // look whether all points are inside of the convex envelope, otherwise add them via degenerated triangles
1025// //->InsertStraddlingPoints(mol, LCList);
1026// mol->GoToFirst();
1027// class TesselPoint *Runner = NULL;
1028// while (!mol->IsEnd()) {
1029// Runner = mol->GetPoint();
1030// Log() << Verbose(1) << "Checking on " << Runner->Name << " ... " << endl;
1031// if (!->IsInnerPoint(Runner, LCList)) {
1032// Log() << Verbose(2) << Runner->Name << " is outside of envelope, adding via degenerated triangles." << endl;
1033// ->AddBoundaryPointByDegeneratedTriangle(Runner, LCList);
1034// } else {
1035// Log() << Verbose(2) << Runner->Name << " is inside of or on envelope." << endl;
1036// }
1037// mol->GoToNext();
1038// }
1039
1040// // Purges surplus triangles.
1041// TesselStruct->RemoveDegeneratedTriangles();
1042
1043 // correct degenerated polygons
1044 TesselStruct->CorrectAllDegeneratedPolygons();
1045
1046 // check envelope for consistency
1047 status = CheckListOfBaselines(TesselStruct);
1048
1049 // write final envelope
1050 CalculateConcavityPerBoundaryPoint(TesselStruct);
1051 StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, "");
1052
1053 if (freeLC)
1054 delete(LCList);
1055
1056 return status;
1057};
1058
1059
1060/** Finds a hole of sufficient size in \a *mols to embed \a *srcmol into it.
1061 * \param *out output stream for debugging
1062 * \param *mols molecules in the domain to embed in between
1063 * \param *srcmol embedding molecule
1064 * \return *Vector new center of \a *srcmol for embedding relative to \a this
1065 */
1066Vector* FindEmbeddingHole(MoleculeListClass *mols, molecule *srcmol)
1067{
1068 Info FunctionInfo(__func__);
1069 Vector *Center = new Vector;
1070 Center->Zero();
1071 // calculate volume/shape of \a *srcmol
1072
1073 // find embedding holes
1074
1075 // if more than one, let user choose
1076
1077 // return embedding center
1078 return Center;
1079};
1080
Note: See TracBrowser for help on using the repository browser.