Changeset 08ef35 for src

Timestamp:

Aug 10, 2009, 4:11:47 PM (15 years ago)

Author:

Frederik Heber <heber@…>

Branches:

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

Children:

570125

Parents:

5c7bf8

Message:

New function ConvexizeNonconvexEnvelope() to calculate the volume of a non-convex envelope.

The central idea is that the volume of a convex envelope is easy to determine as a sum of pyramids with respect to a center inside the envelope. Hence, if we can "reduce" the non-convex envelope to a convex one in such a way that we know the added volume, we may determine the volume of a non-convex envelope.
The nice side effect is that we may use our Find_non_convex_border() algorithm to calculate also the convex envelope.

• We go through all BoundaryPoints and check whether one of its Baselines does not fulfill the ConvexCriterion. If so, we remove it, as it can not be a boundary point on the convex envelope, and re-construct the attached triangles. The added volume is a general tetraeder, whose formula is known.
• FIX: Find_convex_border() - We check whether AddPoint is successful or not.
• builder.cpp: case 'o' - changed to use ConvexizeNonconvexEnvelope() instead of Find_convex_border()
• Tesselation:AddPoint() - now takes second argument which is the index for BPS and always set BPS to either the newly created or the already present point. Return argument discerns between new and already present point.
• Tesselation::BPS, BLS, BTS are now public not private. We have to access them from ConvexizeNonconvexEnvelope()

Location:

src

Files:

• boundary.cpp (modified) (2 diffs)
• boundary.hpp (modified) (1 diff)
• builder.cpp (modified) (1 diff)
• tesselation.cpp (modified) (3 diffs)
• tesselation.hpp (modified) (2 diffs)

src/boundary.cpp

@@ -476 +476 @@
   for (int axis = 0; axis < NDIM; axis++)
     for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++)
-        mol->TesselStruct->AddPoint(runner->second.second);
+        if (!mol->TesselStruct->AddPoint(runner->second.second, 0))
+          *out << Verbose(3) << "WARNING: Point " << *(runner->second.second) << " is already present!" << endl;
 
   *out << Verbose(2) << "I found " << mol->TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl;
@@ -560 +561 @@
 };
 
+/** Creates a convex envelope from a given non-convex one.
+ * -# We go through all PointsOnBoundary.
+ * -# We look at each of its lines and CheckConvexityCriterion().
+ * -# If any returns concave, this Point cannot be on the final convex envelope.
+ * -# Hence, we remove it and re-create all its triangles from its getCircleOfConnectedPoints()
+ * -# We calculate the additional volume
+ * -# We go over the points until none yields a concave spot anymore.
+ * 
+ * Note: This routine - for free - calculates the difference in volume between convex and
+ * non-convex envelope, as the former is easy to calculate - VolumeOfConvexEnvelope() - it
+ * can be used to compute volumes of arbitrary shapes.
+ * \param *out output stream for debugging
+ * \param *TesselStruct non-convex envelope, is changed in return!
+ * \param *configuration contains IsAngstroem()
+ * \return volume difference between the non- and the created convex envelope
+ */
+double ConvexizeNonconvexEnvelope(ofstream *out, class Tesselation *TesselStruct)
+{
+  double volume = 0;
+  class BoundaryPointSet *point = NULL;
+  class BoundaryLineSet *line = NULL;
+  class BoundaryTriangleSet *triangle = NULL;
+  Vector OldPoint, TetraederVector[3];
+  *out << Verbose(0) << "Begin of ConvexizeNonconvexEnvelope" << endl;
+
+  if (TesselStruct == NULL) {
+    *out << Verbose(1) << "ERROR: TesselStruct is empty!" << endl;
+    return volume;
+  }
+
+  bool AllConvex = true;
+  bool Convexity = false;
+  do {
+    AllConvex = true;
+    for (PointMap::iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
+      point = PointRunner->second;
+      *out << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
+      Convexity = true;
+      for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
+        line = LineRunner->second;
+        *out << Verbose(2) << "INFO: Current line is " << *line << "." << endl;
+        if (line->CheckConvexityCriterion(out)) { 
+          // convex
+        } else { // concave
+          Convexity = false; // we have to go again through all ...
+          break;  // no need to check the others too
+        }
+      }
+      AllConvex = AllConvex && Convexity;
+      if (!Convexity) {
+        *out << Verbose(1) << "... point cannot be on convex envelope." << endl;
+        OldPoint.CopyVector(point->node->node);
+        // get list of connected points
+        BoundaryPointSet *Otherpoint = line->GetOtherEndpoint(point);
+        list<TesselPoint*> *CircleofPoints = TesselStruct->getCircleOfConnectedPoints(out, point->node, Otherpoint->node->node);
+        // remove all triangles
+        int *numbers = NULL;
+        int count = 0;
+        int i=0;
+        for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++)
+          count+=LineRunner->second->triangles.size();
+        numbers = new int[count];
+        for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
+          line = LineRunner->second;
+          for (TriangleMap::iterator TriangleRunner = line->triangles.begin(); TriangleRunner != line->triangles.end(); TriangleRunner++) {
+            triangle = TriangleRunner->second;
+            *out << Verbose(2) << "Erasing triangle " << *triangle << "." << endl;
+            numbers[i] = triangle->Nr;
+            TesselStruct->TrianglesOnBoundary.erase(triangle->Nr);
+            delete(triangle);
+            i++;
+          }
+        }
+        *out << Verbose(1) << i << " triangles were removed." << endl;
+
+        // re-create all triangles by going through connected points list
+        list<TesselPoint*>::iterator OtherCircleRunner = CircleofPoints->end();
+        OtherCircleRunner--;
+        i=0;
+        for (list<TesselPoint*>::iterator CircleRunner = CircleofPoints->begin(); CircleRunner != CircleofPoints->end(); CircleRunner++) {
+          *out << Verbose(2) << "Adding new triangle points."<< endl;
+          TesselStruct->AddPoint(point->node, 0);
+          TesselStruct->AddPoint(*CircleRunner, 1);
+          TesselStruct->AddPoint(*OtherCircleRunner, 2);
+          *out << Verbose(2) << "Adding new triangle lines."<< endl;
+          TesselStruct->AddTriangleLine(TesselStruct->BPS[0], TesselStruct->BPS[1], 0); 
+          TesselStruct->AddTriangleLine(TesselStruct->BPS[0], TesselStruct->BPS[2], 1); 
+          TesselStruct->AddTriangleLine(TesselStruct->BPS[1], TesselStruct->BPS[2], 2); 
+          TesselStruct->BTS = new class BoundaryTriangleSet(TesselStruct->BLS, numbers[i]);
+          TesselStruct->TrianglesOnBoundary.insert(TrianglePair(numbers[i], triangle));
+          *out << Verbose(2) << "Created triangle " << *TesselStruct->BTS << "." << endl;
+          // calculate volume summand as a general tetraeder
+          for (int j=0;j<3;j++) {
+            TetraederVector[j].CopyVector(TesselStruct->BPS[j]->node->node);
+            TetraederVector[j].SubtractVector(&OldPoint);
+          }
+          OldPoint.CopyVector(&TetraederVector[0]);
+          OldPoint.VectorProduct(&TetraederVector[1]);
+          volume += 1./6. * fabs(OldPoint.ScalarProduct(&TetraederVector[2]));
+          // advance other shank
+          OtherCircleRunner++;
+          if (OtherCircleRunner == CircleofPoints->end())
+            OtherCircleRunner = CircleofPoints->begin();
+          // advance number
+          i++;
+          if (i > count)
+            *out << Verbose(2) << "WARNING: Maximum of numbers reached!" << endl;
+        }
+        *out << Verbose(1) << i << " triangles were created." << endl;
+        
+        delete[](numbers);
+      }
+    }
+  } while (!AllConvex);
+  *out << Verbose(0) << "End of ConvexizeNonconvexEnvelope" << endl;
+
+  // end
+  return volume;
+};
 
 /** Determines the volume of a cluster.

src/boundary.hpp

@@ -39 +39 @@
 void Find_convex_border(ofstream *out, molecule* mol, class LinkedCell *LCList, const char *filename);
 void Find_non_convex_border(ofstream *out, molecule* mol, class LinkedCell *LC, const char *tempbasename, const double RADIUS);
+double ConvexizeNonconvexEnvelope(ofstream *out, class Tesselation *TesselStruct);
 void Find_next_suitable_point(class BoundaryTriangleSet *BaseTriangle, class BoundaryLineSet *BaseLine, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC);
 Boundaries *GetBoundaryPoints(ofstream *out, molecule *mol);

src/builder.cpp

@@ -1871 +1871 @@
                 cout << Verbose(1) << "Storing tecplot data in " << argv[argptr] << "." << endl;
                 LinkedCell LCList(mol, 10.);
-                class Tesselation *TesselStruct = NULL;
-                Find_convex_border((ofstream *)&cout, mol, &LCList, argv[argptr]);
-                double clustervolume = VolumeOfConvexEnvelope((ofstream *)&cout, TesselStruct, &configuration);
-                cout << Verbose(0) << "The tesselated surface area is " << clustervolume << "." << endl;
-                delete(TesselStruct);
+                //Find_convex_border((ofstream *)&cout, mol, &LCList, argv[argptr]);
+                Find_non_convex_border((ofstream *)&cout, mol, &LCList, argv[argptr], 10.);
+
+                double volumedifference = ConvexizeNonconvexEnvelope((ofstream *)&cout, mol->TesselStruct);
+                double clustervolume = VolumeOfConvexEnvelope((ofstream *)&cout, mol->TesselStruct, &configuration);
+                cout << Verbose(0) << "The tesselated surface area is " << clustervolume << " " << (configuration.GetIsAngstroem() ? "angstrom" : "atomiclength") << "^3." << endl;
+                cout << Verbose(0) << "The non-convex tesselated surface area is " << clustervolume-volumedifference << " " << (configuration.GetIsAngstroem() ? "angstrom" : "atomiclength") << "^3." << endl;
                 argptr+=1;
               }

src/tesselation.cpp

@@ -208 +208 @@
  * \return NULL - if endpoint not contained in BoundaryLineSet, or pointer to BoundaryPointSet otherwise
  */
-inline class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(class BoundaryPointSet *point)
+class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(class BoundaryPointSet *point)
 {
   if (endpoints[0] == point)
@@ -1044 +1044 @@
         // add Walker to boundary points
         *out << Verbose(2) << "Adding " << *Walker << " to BoundaryPoints." << endl;
-        AddPoint(Walker);
-        if (BPS[0] != NULL)
+        if (AddPoint(Walker,0))
           NewPoint = BPS[0];
         else
@@ -1100 +1099 @@
 /** Adds an point to the tesselation::PointsOnBoundary list.
  * \param *Walker point to add
- */
-void
-Tesselation::AddPoint(TesselPoint *Walker)
+ * \param n TesselStruct::BPS index to put pointer into
+ * \return true - new point was added, false - point already present
+ */
+bool
+Tesselation::AddPoint(TesselPoint *Walker, int n)
 {
   PointTestPair InsertUnique;
-  BPS[0] = new class BoundaryPointSet(Walker);
-  InsertUnique = PointsOnBoundary.insert(PointPair(Walker->nr, BPS[0]));
-  if (InsertUnique.second) // if new point was not present before, increase counter
+  BPS[n] = new class BoundaryPointSet(Walker);
+  InsertUnique = PointsOnBoundary.insert(PointPair(Walker->nr, BPS[n]));
+  if (InsertUnique.second) { // if new point was not present before, increase counter
     PointsOnBoundaryCount++;
-  else {
-    delete(BPS[0]);
-    BPS[0] = NULL;
+    return true;
+  } else {
+    delete(BPS[n]);
+    BPS[n] = InsertUnique.first->second;
+    return false;
   }
 }

src/tesselation.hpp

@@ -190 +190 @@
     virtual ~Tesselation();
 
-    void AddPoint(TesselPoint *Walker);
+    bool AddPoint(TesselPoint *Walker, int n);
     void AddTrianglePoint(TesselPoint* Candidate, int n);
     void AddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n);
@@ -236 +236 @@
     virtual bool IsEnd();
 
-  private:
-    class BoundaryPointSet *TPS[3]; //this is a Storage for pointers to triangle points, this and BPS[2] needed due to AddLine restrictions
     class BoundaryPointSet *BPS[2];
     class BoundaryLineSet *BLS[3];
     class BoundaryTriangleSet *BTS;
+
+  private:
+    class BoundaryPointSet *TPS[3]; //this is a Storage for pointers to triangle points, this and BPS[2] needed due to AddLine restrictions
 
     PointMap::iterator InternalPointer;