source: src/boundary.cpp@ ba9f5b

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 ba9f5b was bdc91e, checked in by Frederik Heber <heber@…>, 14 years ago

MEMFIXES: Tesselation routines were leaking memory.

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

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