Changeset 2319ed for src/vector.cpp

Timestamp:

Aug 3, 2009, 8:10:09 AM (16 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, Candidate_v1.7.0, 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:

0dbddc, 357fba

Parents:

e1589e

Message:

We are one step further in fixing the convex hull: There are two functions of Saskia Metzler missing, but then we may proceed with testing whether the simple correction scheme of the convex envelope works, but one thing: Right now we cannot associate a Tesselation to its molecule as the snake bites it's one tail. Hence, the next commit will consist of fixing this bad-OOP issue.

• Makefile.am: Just some alphabetical resorting.
• atom::atom() new copy constructor
• builder.cpp: some output for cluster volume, molecule::AddCopyAtom() uses new copy constructor
• FillBoxWithMolecule() - new function to fill the remainder of the simulation box with some given filler molecules. Makes explicit use of the tesselated surfaces
• find_convex_border() - InsertStraddlingPoints() and CorrectConcaveBaselines() is called to correct for atoms outside the envelope and caused-by concave points
• Tesselation::InsertStraddlingPoints() enlarges the envelope for all atoms found outside, Tesselation::CorrectConcaveBaselines() corrects all found baselines if the adjacent triangles are concave by flipping.
• boundary.cpp: Lots of helper routines for stuff further below:
  • BoundaryLineSet::IsConnectedTo() checks whether lines have common endpoint
  • BoundaryLineSet::CheckConvexityCriterion() checks whether adjacent triangles are convex or concave
  • BoundaryLineSet::ContainsBoundaryPoint() checks whether line contains a given boundary point
  • BoundaryTriangleSet::ContainsBoundaryLine() and BoundaryTriangleSet::ContainsBoundaryPoint() check whether triangle contains boundary line or point
  • BoundaryTriangleSet::GetIntersectionInsideTriangle() calculates intersection point between line and triangle,
• The following routines are needed to check whether point is in- or outside:
  • Vector::GetIntersectionWithPlane() calculates intersection point of a plane and a line.
  • Vector::GetIntersectionOfTwoLinesOnPlane() intersection of two lines on the same plane.
• FIX: Tesselation::AddPoint() - newly created BoundaryPoint is removed if already present.

Problem: We want to associate a Tesselation class with each molecule class. However, so far we have to know about atoms and bond and molecules inside the Tesselation. We have to remove this dependency and create some intermediate class which enables/encapsulates the access to Vectors, e.g. hidden inside the atom class. This is also good OOP! The Tesselation also only needs a set of Vectors, not more!

File:

• src/vector.cpp (modified) (2 diffs)

src/vector.cpp

@@ -196 +196 @@
 };
 
+/** Calculates the intersection point between a line defined by \a *LineVector and \a *LineVector2 and a plane defined by \a *Normal and \a *PlaneOffset.
+ * According to [Bronstein] the vectorial plane equation is:
+ *   -# \f$\stackrel{r}{\rightarrow} \cdot \stackrel{N}{\rightarrow} + D = 0\f$,
+ * where \f$\stackrel{r}{\rightarrow}\f$ is the vector to be testet, \f$\stackrel{N}{\rightarrow}\f$ is the plane's normal vector and
+ * \f$D = - \stackrel{a}{\rightarrow} \stackrel{N}{\rightarrow}\f$, the offset with respect to origin, if \f$\stackrel{a}{\rightarrow}\f$,
+ * is an offset vector onto the plane. The line is parametrized by \f$\stackrel{x}{\rightarrow} + k \stackrel{t}{\rightarrow}\f$, where
+ * \f$\stackrel{x}{\rightarrow}\f$ is the offset and \f$\stackrel{t}{\rightarrow}\f$ the directional vector (NOTE: No need to normalize
+ * the latter). Inserting the parametrized form into the plane equation and solving for \f$k\f$, which we insert then into the parametrization
+ * of the line yields the intersection point on the plane.
+ * \param *out output stream for debugging
+ * \param *PlaneNormal Plane's normal vector
+ * \param *PlaneOffset Plane's offset vector
+ * \param *LineVector first vector of line
+ * \param *LineVector2 second vector of line
+ * \return true -  \a this contains intersection point on return, false - line is parallel to plane
+ */
+bool Vector::GetIntersectionWithPlane(ofstream *out, Vector *PlaneNormal, Vector *PlaneOffset, Vector *LineVector, Vector *LineVector2)
+{
+  double factor;
+  Vector Direction;
+
+  // find intersection of a line defined by Offset and Direction with a  plane defined by triangle
+  Direction.CopyVector(LineVector2);
+  Direction.SubtractVector(LineVector);
+  if (Direction.ScalarProduct(PlaneNormal) < MYEPSILON)
+    return false;
+  factor = LineVector->ScalarProduct(PlaneNormal)*(-PlaneOffset->ScalarProduct(PlaneNormal))/(Direction.ScalarProduct(PlaneNormal));
+  Direction.Scale(factor);
+  CopyVector(LineVector);
+  SubtractVector(&Direction);
+
+  // test whether resulting vector really is on plane
+  Direction.CopyVector(this);
+  Direction.SubtractVector(PlaneOffset);
+  if (Direction.ScalarProduct(PlaneNormal) > MYEPSILON)
+    return true;
+  else
+    return false;
+};
+
+/** Calculates the intersection of the two lines that are both on the same plane.
+ * Note that we do not check whether they are on the same plane.
+ * \param *out output stream for debugging
+ * \param *Line1a first vector of first line
+ * \param *Line1b second vector of first line
+ * \param *Line2a first vector of second line
+ * \param *Line2b second vector of second line
+ * \return true - \a this will contain the intersection on return, false - lines are parallel
+ */
+bool Vector::GetIntersectionOfTwoLinesOnPlane(ofstream *out, Vector *Line1a, Vector *Line1b, Vector *Line2a, Vector *Line2b)
+{
+  double k1,a1,h1,b1,k2,a2,h2,b2;
+  // equation for line 1
+  if (fabs(Line1a->x[0] - Line2a->x[0]) < MYEPSILON) {
+    k1 = 0;
+    h1 = 0;
+  } else {
+    k1 = (Line1a->x[1] - Line2a->x[1])/(Line1a->x[0] - Line2a->x[0]);
+    h1 = (Line1a->x[2] - Line2a->x[2])/(Line1a->x[0] - Line2a->x[0]);
+  }
+  a1 = 0.5*((Line1a->x[1] + Line2a->x[1]) - k1*(Line1a->x[0] + Line2a->x[0]));
+  b1 = 0.5*((Line1a->x[2] + Line2a->x[2]) - h1*(Line1a->x[0] + Line2a->x[0]));
+
+  // equation for line 2
+  if (fabs(Line2a->x[0] - Line2a->x[0]) < MYEPSILON) {
+    k1 = 0;
+    h1 = 0;
+  } else {
+    k1 = (Line2a->x[1] - Line2a->x[1])/(Line2a->x[0] - Line2a->x[0]);
+    h1 = (Line2a->x[2] - Line2a->x[2])/(Line2a->x[0] - Line2a->x[0]);
+  }
+  a1 = 0.5*((Line2a->x[1] + Line2a->x[1]) - k1*(Line2a->x[0] + Line2a->x[0]));
+  b1 = 0.5*((Line2a->x[2] + Line2a->x[2]) - h1*(Line2a->x[0] + Line2a->x[0]));
+
+  // calculate cross point
+  if (fabs((a1-a2)*(h1-h2) - (b1-b2)*(k1-k2)) < MYEPSILON) {
+    x[0] = (a2-a1)/(k1-k2);
+    x[1] = (k1*a2-k2*a1)/(k1-k2);
+    x[2] = (h1*b2-h2*b1)/(h1-h2);
+    *out << Verbose(4) << "Lines do intersect at " << this << "." << endl;
+    return true;
+  } else {
+    *out << Verbose(4) << "Lines do not intersect." << endl;
+    return false;
+  }
+};
+
 /** Calculates the projection of a vector onto another \a *y.
  * \param *y array to second vector
@@ -453 +540 @@
 };
 
-/** Do a matrix multiplication with \a *matrix' inverse.
+/** Do a matrix multiplication with the \a *A' inverse.
  * \param *matrix NDIM_NDIM array
  */