Changeset 093645 for src

Timestamp:

Aug 17, 2009, 10:48:03 AM (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:

0077b5, 54a746, f1cccd

Parents:

16d866

Message:

Rewrite of ConvexizeNonconvexEnvelope() with three steps: Remove concave spots, make all lines best possible, removed concave lines. Is not working yet!

• fixes to write_tecplot_file(): Checks for NULL TesselStruct, write variable extra variable U (contained in BoundaryPointSet), and use .size() instead of ...onBoundaryCount for numbers
• Find_convex_border(): while loop till break in CorrectConcaveBaselines() replaced by a single loop over all boundary lines with PickFarthestofTwoBaselines()
• ConvexizeNonconvexEnvelope(): documentation, needs molecule and filename to save tesselations, three steps consisting of single loops: RemovePointFromTesselatedSurface(), PickFarthestofTwoBaselines(), IsConvexRectangle(). And we calculate remaining concavity by CheckConvexityCriterion for all lines to a boundary point and store it as "U" in the tecplot files.

Location:

src

Files:

• boundary.cpp (modified) (7 diffs)
• boundary.hpp (modified) (1 diff)

src/boundary.cpp

@@ -229 +229 @@
 void write_tecplot_file(ofstream *out, ofstream *tecplot, class Tesselation *TesselStruct, class molecule *mol, int N)
 {
-  if (tecplot != NULL) {
+  if ((tecplot != NULL) && (TesselStruct != NULL)) {
+    // write header
     *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
-    *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\"" << endl;
-    *tecplot << "ZONE T=\"TRIANGLES" << N << "\", N=" << TesselStruct->PointsOnBoundaryCount << ", E=" << TesselStruct->TrianglesOnBoundaryCount << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
+    *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\" \"U\"" << endl;
+    *tecplot << "ZONE T=\"TRIANGLES" << N << "\", N=" << TesselStruct->PointsOnBoundary.size() << ", E=" << TesselStruct->TrianglesOnBoundary.size() << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
     int *LookupList = new int[mol->AtomCount];
     for (int i = 0; i < mol->AtomCount; i++)
@@ -244 +245 @@
       Walker = target->second->node;
       LookupList[Walker->nr] = Counter++;
-      *tecplot << Walker->node->x[0] << " " << Walker->node->x[1] << " " << Walker->node->x[2] << " " << endl;
+      *tecplot << Walker->node->x[0] << " " << Walker->node->x[1] << " " << Walker->node->x[2] << " " << target->second->value << endl;
     }
     *tecplot << endl;
@@ -476 +477 @@
   for (int axis = 0; axis < NDIM; axis++)
     for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++)
-        if (!mol->TesselStruct->AddPoint(runner->second.second, 0))
+        if (!mol->TesselStruct->AddBoundaryPoint(runner->second.second, 0))
           *out << Verbose(3) << "WARNING: Point " << *(runner->second.second) << " is already present!" << endl;
 
@@ -524 +525 @@
 
   // 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
-  if (!mol->TesselStruct->CorrectConcaveBaselines(out))
-    *out << Verbose(1) << "Correction of concave baselines failed!" << endl;
+  bool AllConvex;
+  class BoundaryLineSet *line = NULL;
+  do {
+    AllConvex = true;
+    for (LineMap::iterator LineRunner = mol->TesselStruct->LinesOnBoundary.begin(); LineRunner != mol->TesselStruct->LinesOnBoundary.end(); LineRunner++) {
+      line = LineRunner->second;
+      *out << Verbose(1) << "INFO: Current line is " << *line << "." << endl;
+      if (!line->CheckConvexityCriterion(out)) {
+        *out << Verbose(1) << "... line " << *line << " is concave, flipping it." << endl;
+
+        // flip the line
+        if (!mol->TesselStruct->PickFarthestofTwoBaselines(out, line))
+          *out << Verbose(1) << "ERROR: Correction of concave baselines failed!" << endl;
+        else
+          *out << Verbose(1) << "INFO: Correction of concave baselines worked." << endl;
+      }
+    }
+  } while (!AllConvex);
 
   // 3e. we need another correction here, for TesselPoints that are below the surface (i.e. have an odd number of concave triangles surrounding it)
@@ -562 +579 @@
 
 /** 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.
- * 
+ * -# First step, remove concave spots, i.e. singular "dents"
+ *   -# We go through all PointsOnBoundary.
+ *   -# We CheckConvexityCriterion() for all its lines.
+ *   -# If all its lines are concave, it cannot be on the convex envelope.
+ *   -# Hence, we remove it and re-create all its triangles from its getCircleOfConnectedPoints()
+ *   -# We calculate the additional volume.
+ *   -# We go over all lines until none yields a concavity anymore.
+ * -# Second step, remove concave lines, i.e. line-shape "dents"
+ *   -# We go through all LinesOnBoundary
+ *   -# We CheckConvexityCriterion()
+ *   -# If it returns concave, we flip the line in this quadrupel of points (abusing the degeneracy of the tesselation)
+ *   -# We CheckConvexityCriterion(),
+ *   -# if it's concave, we continue
+ *   -# if not, we mark an error and stop
  * 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
@@ -574 +598 @@
  * \param *out output stream for debugging
  * \param *TesselStruct non-convex envelope, is changed in return!
- * \param *configuration contains IsAngstroem()
+ * \param *mol molecule
+ * \param *filename name of file
  * \return volume difference between the non- and the created convex envelope
  */
-double ConvexizeNonconvexEnvelope(ofstream *out, class Tesselation *TesselStruct)
+double ConvexizeNonconvexEnvelope(ofstream *out, class Tesselation *TesselStruct, molecule *mol, char *filename)
 {
   double volume = 0;
   class BoundaryPointSet *point = NULL;
   class BoundaryLineSet *line = NULL;
-  class BoundaryTriangleSet *triangle = NULL;
-  Vector OldPoint, TetraederVector[3];
+  bool Concavity;
+  PointMap::iterator PointRunner, PointAdvance;
+  LineMap::iterator LineRunner, LineAdvance;
+  TriangleMap::iterator TriangleRunner, TriangleAdvance;
+
   *out << Verbose(0) << "Begin of ConvexizeNonconvexEnvelope" << endl;
 
+  // check whether there is something to work on
   if (TesselStruct == NULL) {
     *out << Verbose(1) << "ERROR: TesselStruct is empty!" << endl;
@@ -591 +620 @@
   }
 
-  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);
+  // First step: RemovePointFromTesselatedSurface
+  PointRunner = TesselStruct->PointsOnBoundary.begin();
+  PointAdvance = PointRunner; // we need an advanced point, as the PointRunner might get removed
+  while (PointRunner != TesselStruct->PointsOnBoundary.end()) {
+    PointAdvance++;
+    point = PointRunner->second;
+    *out << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
+    Concavity = true;
+    for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
+      line = LineRunner->second;
+      *out << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
+      Concavity = Concavity && (!line->CheckConvexityCriterion(out));
+    }
+    if (Concavity) {
+      *out << Verbose(1) << "... point " << *point << " cannot be on convex envelope." << endl;
+      // remove the point
+      TesselStruct->RemovePointFromTesselatedSurface(out, point);
+    }
+    PointRunner = PointAdvance;
+  }
+
+
+//  // print all lines
+//  LineRunner = TesselStruct->LinesOnBoundary.begin();
+//  LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
+//  *out << Verbose(1) << "Printing all boundary lines for debugging." << endl;
+//  while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
+//    LineAdvance++;
+//    line = LineRunner->second;
+//    *out << Verbose(2) << "INFO: Current line is " << *line << "." << endl;
+//    if (LineAdvance != TesselStruct->LinesOnBoundary.end())
+//      *out << Verbose(2) << "INFO: Next line will be " << *(LineAdvance->second) << "." << endl;
+//    LineRunner = LineAdvance;
+//  }
+//
+//  // print all triangles
+//  TriangleRunner = TesselStruct->TrianglesOnBoundary.begin();
+//  TriangleAdvance = TriangleRunner;  // we need an advanced line, as the LineRunner might get removed
+//  *out << Verbose(1) << "Printing all boundary triangles for debugging." << endl;
+//  while (TriangleRunner != TesselStruct->TrianglesOnBoundary.end()) {
+//    TriangleAdvance++;
+//    *out << Verbose(2) << "INFO: Current triangle is " << *(TriangleRunner->second) << "." << endl;
+//    if (TriangleAdvance != TesselStruct->TrianglesOnBoundary.end())
+//      *out << Verbose(2) << "INFO: Next triangle will be " << *(TriangleAdvance->second) << "." << endl;
+//    TriangleRunner = TriangleAdvance;
+//  }
+
+  // second step: PickFarthestofTwoBaselines
+  LineRunner = TesselStruct->LinesOnBoundary.begin();
+  LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
+  while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
+    LineAdvance++;
+    line = LineRunner->second;
+    *out << Verbose(1) << "INFO: Picking farthest baseline for line is " << *line << "." << endl;
+    if (LineAdvance != TesselStruct->LinesOnBoundary.end())
+    // take highest of both lines
+    TesselStruct->PickFarthestofTwoBaselines(out, line);
+    LineRunner = LineAdvance;
+  }
+
+  // calculate remaining concavity
+  for (PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
+    point = PointRunner->second;
+    *out << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
+    point->value = 0;
+    for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
+      line = LineRunner->second;
+      *out << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
+      if (!line->CheckConvexityCriterion(out))
+        point->value += 1;
+    }
+  }
+
+  // 4. Store triangles in tecplot file
+  if (filename != NULL) {
+    if (DoTecplotOutput) {
+      string OutputName(filename);
+      OutputName.append("_intermed");
+      OutputName.append(TecplotSuffix);
+      ofstream *tecplot = new ofstream(OutputName.c_str());
+      write_tecplot_file(out, tecplot, mol->TesselStruct, mol, 0);
+      tecplot->close();
+      delete(tecplot);
+    }
+    if (DoRaster3DOutput) {
+      string OutputName(filename);
+      OutputName.append("_intermed");
+      OutputName.append(Raster3DSuffix);
+      ofstream *rasterplot = new ofstream(OutputName.c_str());
+      write_raster3d_file(out, rasterplot, mol->TesselStruct, mol);
+      rasterplot->close();
+      delete(rasterplot);
+    }
+  }
+
+  // third step: IsConvexRectangle
+  LineRunner = TesselStruct->LinesOnBoundary.begin();
+  LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
+  while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
+    LineAdvance++;
+    line = LineRunner->second;
+    *out << Verbose(1) << "INFO: Current line is " << *line << "." << endl;
+    if (LineAdvance != TesselStruct->LinesOnBoundary.end())
+      *out << Verbose(1) << "INFO: Next line will be " << *(LineAdvance->second) << "." << endl;
+    if (!line->CheckConvexityCriterion(out)) {
+      *out << Verbose(1) << "... line " << *line << " is concave, flipping it." << endl;
+
+      // take highest of both lines
+      point = TesselStruct->IsConvexRectangle(out, line);
+      if (point != NULL)
+        TesselStruct->RemovePointFromTesselatedSurface(out, point);
+    }
+    LineRunner = LineAdvance;
+  }
+
+  // calculate remaining concavity
+  for (PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
+    point = PointRunner->second;
+    *out << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
+    point->value = 0;
+    for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
+      line = LineRunner->second;
+      *out << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
+      if (!line->CheckConvexityCriterion(out))
+        point->value += 1;
+    }
+  }
+
+  // 4. Store triangles in tecplot file
+  if (filename != NULL) {
+    if (DoTecplotOutput) {
+      string OutputName(filename);
+      OutputName.append(TecplotSuffix);
+      ofstream *tecplot = new ofstream(OutputName.c_str());
+      write_tecplot_file(out, tecplot, mol->TesselStruct, mol, 0);
+      tecplot->close();
+      delete(tecplot);
+    }
+    if (DoRaster3DOutput) {
+      string OutputName(filename);
+      OutputName.append(Raster3DSuffix);
+      ofstream *rasterplot = new ofstream(OutputName.c_str());
+      write_raster3d_file(out, rasterplot, mol->TesselStruct, mol);
+      rasterplot->close();
+      delete(rasterplot);
+    }
+  }
+
+  // end
   *out << Verbose(0) << "End of ConvexizeNonconvexEnvelope" << endl;
-
-  // end
   return volume;
 };

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);
+double ConvexizeNonconvexEnvelope(ofstream *out, class Tesselation *TesselStruct, molecule *mol, char *filename);
 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);