Changeset 5de9da

Timestamp:

Jul 14, 2010, 2:52:07 PM (15 years ago)

Author:

Tillmann Crueger <crueger@…>

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:

ef84ca

Parents:

06300d

Message:

Added some important Methods to Shape class hierachy

• Shape::isOnSurface()
• Shape::getNormal()

Location:

src/Shapes

Files:

• BaseShapes.cpp (modified) (2 diffs)
• BaseShapes_impl.hpp (modified) (1 diff)
• Shape.cpp (modified) (5 diffs)
• Shape.hpp (modified) (2 diffs)
• ShapeOps.cpp (modified) (7 diffs)
• ShapeOps_impl.hpp (modified) (1 diff)
• Shape_impl.hpp (modified) (7 diffs)

src/Shapes/BaseShapes.cpp

@@ -10 +10 @@
 
 #include "vector.hpp"
+#include "Helpers/Assert.hpp"
+
+#include <cmath>
 
 bool Sphere_impl::isInside(const Vector &point){
   return point.NormSquared()<=1;
+}
+
+bool Sphere_impl::isOnSurface(const Vector &point){
+  return fabs(point.NormSquared()-1)<MYEPSILON;
+}
+
+Vector Sphere_impl::getNormal(const Vector &point) throw(NotOnSurfaceException){
+  if(!isOnSurface(point)){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  return point;
 }
 
@@ -21 +35 @@
 
 bool Cuboid_impl::isInside(const Vector &point){
-  return point[0]<=1 && point[1]<=1 && point[2]<=1;
+  return fabs(point[0])<=1 && fabs(point[1])<=1 && fabs(point[2])<=1;
+}
+
+bool Cuboid_impl::isOnSurface(const Vector &point){
+  bool retVal = isInside(point);
+  // test all borders of the cuboid
+  // double fabs
+  retVal = retVal &&
+           ((fabs(fabs(point[0])-1)  < MYEPSILON) ||
+            (fabs(fabs(point[1])-1)  < MYEPSILON) ||
+            (fabs(fabs(point[2])-1)  < MYEPSILON));
+  return retVal;
+}
+
+Vector Cuboid_impl::getNormal(const Vector &point) throw(NotOnSurfaceException){
+  if(!isOnSurface(point)){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  Vector res;
+  // figure out on which sides the Vector lies (maximum 3, when it is in a corner)
+  for(int i=NDIM;i--;){
+    if(fabs(fabs(point[i])-1)<MYEPSILON){
+      // add the scaled (-1/+1) Vector to the set of surface vectors
+      res[i] = point[i];
+    }
+  }
+  ASSERT(res.NormSquared()>=1 && res.NormSquared()<=3,"To many or to few sides found for this Vector");
+
+  res.Normalize();
+  return res;
 }
 
 Shape Cuboid(){
-  Shape::impl_ptr impl = Shape::impl_ptr(new Sphere_impl());
+  Shape::impl_ptr impl = Shape::impl_ptr(new Cuboid_impl());
   return Shape(impl);
 }

src/Shapes/BaseShapes_impl.hpp

@@ -13 +13 @@
 class Sphere_impl : public Shape_impl {
   virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException);
 };
 
 class Cuboid_impl : public Shape_impl {
   virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException);
 };
 

src/Shapes/Shape.cpp

@@ -9 +9 @@
 #include "Shape_impl.hpp"
 
+#include "Helpers/Assert.hpp"
+
 Shape::Shape(const Shape& src) :
   impl(src.getImpl())
@@ -17 +19 @@
 bool Shape::isInside(const Vector &point) const{
   return impl->isInside(point);
+}
+
+bool Shape::isOnSurface(const Vector &point) const{
+  return impl->isOnSurface(point);
+}
+
+Vector Shape::getNormal(const Vector &point) const throw (NotOnSurfaceException){
+  return impl->getNormal(point);
 }
 
@@ -65 +75 @@
 }
 
+bool AndShape_impl::isOnSurface(const Vector &point){
+  // check the number of surfaces that this point is on
+  int surfaces =0;
+  surfaces += lhs->isOnSurface(point);
+  surfaces += rhs->isOnSurface(point);
+
+  switch(surfaces){
+    case 0:
+      return false;
+      // no break necessary
+    case 1:
+      // if it is inside for the object where it does not lie on
+      // the surface the whole point lies inside
+      return (lhs->isOnSurface(point) && !rhs->isInside(point)) ||
+             (rhs->isOnSurface(point) && !lhs->isInside(point));
+      // no break necessary
+    case 2:
+      {
+        // it lies on both Shapes... could be an edge or an inner point
+        // test the direction of the normals
+        Vector direction=lhs->getNormal(point)+rhs->getNormal(point);
+        // if the directions are opposite we lie on the inside
+        return !direction.IsZero();
+      }
+      // no break necessary
+    default:
+      // if this happens there is something wrong
+      ASSERT(0,"Default case should have never been used");
+  }
+  return false; // never reached
+}
+
+Vector AndShape_impl::getNormal(const Vector &point) throw (NotOnSurfaceException){
+  Vector res;
+  if(!isOnSurface(point)){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  res += lhs->isOnSurface(point)?lhs->getNormal(point):zeroVec;
+  res += rhs->isOnSurface(point)?rhs->getNormal(point):zeroVec;
+  res.Normalize();
+  return res;
+}
+
 Shape operator&&(const Shape &lhs,const Shape &rhs){
   Shape::impl_ptr newImpl = Shape::impl_ptr(new AndShape_impl(getShapeImpl(lhs),getShapeImpl(rhs)));
@@ -82 +135 @@
 }
 
+bool OrShape_impl::isOnSurface(const Vector &point){
+  // check the number of surfaces that this point is on
+  int surfaces =0;
+  surfaces += lhs->isOnSurface(point);
+  surfaces += rhs->isOnSurface(point);
+
+  switch(surfaces){
+    case 0:
+      return false;
+      // no break necessary
+    case 1:
+      // if it is inside for the object where it does not lie on
+      // the surface the whole point lies inside
+      return (lhs->isOnSurface(point) && rhs->isInside(point)) ||
+             (rhs->isOnSurface(point) && lhs->isInside(point));
+      // no break necessary
+    case 2:
+      {
+        // it lies on both Shapes... could be an edge or an inner point
+        // test the direction of the normals
+        Vector direction=lhs->getNormal(point)+rhs->getNormal(point);
+        // if the directions are opposite we lie on the inside
+        return !direction.IsZero();
+      }
+      // no break necessary
+    default:
+      // if this happens there is something wrong
+      ASSERT(0,"Default case should have never been used");
+  }
+  return false; // never reached
+}
+
+Vector OrShape_impl::getNormal(const Vector &point) throw (NotOnSurfaceException){
+  Vector res;
+  if(!isOnSurface(point)){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  res += lhs->isOnSurface(point)?lhs->getNormal(point):zeroVec;
+  res += rhs->isOnSurface(point)?rhs->getNormal(point):zeroVec;
+  res.Normalize();
+  return res;
+}
+
 Shape operator||(const Shape &lhs,const Shape &rhs){
   Shape::impl_ptr newImpl = Shape::impl_ptr(new OrShape_impl(getShapeImpl(lhs),getShapeImpl(rhs)));
@@ -99 +195 @@
 }
 
+bool NotShape_impl::isOnSurface(const Vector &point){
+  return arg->isOnSurface(point);
+}
+
+Vector NotShape_impl::getNormal(const Vector &point) throw(NotOnSurfaceException){
+  return -1*arg->getNormal(point);
+}
+
 Shape operator!(const Shape &arg){
   Shape::impl_ptr newImpl = Shape::impl_ptr(new NotShape_impl(getShapeImpl(arg)));

src/Shapes/Shape.hpp

@@ -10 +10 @@
 
 #include <boost/shared_ptr.hpp>
+
+#include "Exceptions/NotOnSurfaceException.hpp"
 
 class Vector;
@@ -25 +27 @@
 
   bool isInside(const Vector &point) const;
+  bool isOnSurface(const Vector &point) const;
+  Vector getNormal(const Vector &point) const throw(NotOnSurfaceException);
 
   Shape &operator=(const Shape& rhs);

src/Shapes/ShapeOps.cpp

@@ -11 +11 @@
 #include "Helpers/Assert.hpp"
 
+/*************** Base case ***********************/
+
+ShapeOpsBase_impl::ShapeOpsBase_impl(const Shape::impl_ptr &_arg) :
+  arg(_arg){}
+
+ShapeOpsBase_impl::~ShapeOpsBase_impl(){}
+
+bool ShapeOpsBase_impl::isInside(const Vector &point){
+  return arg->isInside(translateIn(point));
+}
+
+bool ShapeOpsBase_impl::isOnSurface(const Vector &point){
+  return arg->isOnSurface(translateIn(point));
+}
+
+Vector ShapeOpsBase_impl::getNormal(const Vector &point) throw (NotOnSurfaceException){
+  Vector helper = translateIn(point);
+  if(!arg->isOnSurface(helper)){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  return translateOutNormal(arg->getNormal(helper));
+}
+
 /********************* Resize ********************/
 
 Resize_impl::Resize_impl(const Shape::impl_ptr &_arg,double _size) :
-  arg(_arg), size(_size)
+  ShapeOpsBase_impl(_arg), size(_size)
 {
   ASSERT(size>0,"Cannot resize a Shape to size zero or below");
@@ -21 +44 @@
 Resize_impl::~Resize_impl(){}
 
-bool Resize_impl::isInside(const Vector& point){
-  return arg->isInside((1/size) * point);
+Vector Resize_impl::translateIn(const Vector& point){
+  return (1/size) * point;
+}
+
+Vector Resize_impl::translateOutPos(const Vector& point){
+  return size * point;
+}
+
+Vector Resize_impl::translateOutNormal(const Vector& point){
+  return point;
 }
 
@@ -33 +64 @@
 
 Translate_impl::Translate_impl(const Shape::impl_ptr &_arg, const Vector &_offset) :
-  arg(_arg),offset(_offset)
+  ShapeOpsBase_impl(_arg),offset(_offset)
 {}
 
 Translate_impl::~Translate_impl(){}
 
-bool Translate_impl::isInside(const Vector& point){
-  return arg->isInside(point-offset);
+Vector Translate_impl::translateIn(const Vector& point){
+  return point-offset;
+}
+
+Vector Translate_impl::translateOutPos(const Vector& point){
+  return point+offset;
+}
+
+Vector Translate_impl::translateOutNormal(const Vector& point){
+  return point;
 }
 
@@ -50 +89 @@
 
 Stretch_impl::Stretch_impl(const Shape::impl_ptr &_arg, const Vector &_factors) :
-  arg(_arg),factors(_factors)
+  ShapeOpsBase_impl(_arg),factors(_factors)
 {
   ASSERT(factors[0]>0,"cannot stretch a shape by a negative amount");
@@ -62 +101 @@
 Stretch_impl::~Stretch_impl(){}
 
-bool Stretch_impl::isInside(const Vector& point){
+Vector Stretch_impl::translateIn(const Vector& point){
   Vector helper=point;
   helper.ScaleAll(reciFactors);
-  return arg->isInside(helper);
+  return helper;
+}
+
+Vector Stretch_impl::translateOutPos(const Vector& point){
+  Vector helper=point;
+  helper.ScaleAll(factors);
+  return helper;
+}
+
+Vector Stretch_impl::translateOutNormal(const Vector& point){
+  Vector helper=point;
+  // the normalFactors are derived from appearances of the factors
+  // with in the vectorproduct
+  Vector normalFactors;
+  normalFactors[0]=factors[1]*factors[2];
+  normalFactors[1]=factors[0]*factors[2];
+  normalFactors[2]=factors[0]*factors[1];
+  helper.ScaleAll(normalFactors);
+  return helper;
 }
 
@@ -76 +133 @@
 
 Transform_impl::Transform_impl(const Shape::impl_ptr &_arg, const Matrix &_transformation) :
-  arg(_arg),transformation(_transformation)
+  ShapeOpsBase_impl(_arg),transformation(_transformation)
 {
   transformationInv = transformation.invert();
@@ -83 +140 @@
 Transform_impl::~Transform_impl(){}
 
-bool Transform_impl::isInside(const Vector& point){
-  return arg->isInside(transformationInv * point);
+Vector Transform_impl::translateIn(const Vector& point){
+  return transformationInv * point;
+}
+
+Vector Transform_impl::translateOutPos(const Vector& point){
+  return transformation * point;
+}
+
+Vector Transform_impl::translateOutNormal(const Vector& point){
+  Matrix mat = transformation.determinant() * transformation.invert().transpose();
+  return mat * point;
 }
 

src/Shapes/ShapeOps_impl.hpp

@@ -13 +13 @@
 #include "Matrix.hpp"
 
-class Resize_impl :  public Shape_impl
+class ShapeOpsBase_impl : public Shape_impl{
+public:
+  ShapeOpsBase_impl(const Shape::impl_ptr&);
+  virtual ~ShapeOpsBase_impl();
+  virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw (NotOnSurfaceException);
+protected:
+  virtual Vector translateIn(const Vector &point)=0;
+  virtual Vector translateOutPos(const Vector &point)=0;
+  virtual Vector translateOutNormal(const Vector &point)=0;
+private:
+  Shape::impl_ptr arg;
+};
+
+class Resize_impl :  public ShapeOpsBase_impl
 {
 public:
   Resize_impl(const Shape::impl_ptr&,double);
   virtual ~Resize_impl();
-  virtual bool isInside(const Vector& point);
+protected:
+  virtual Vector translateIn(const Vector &point);
+  virtual Vector translateOutPos(const Vector &point);
+  virtual Vector translateOutNormal(const Vector &point);
 private:
-  Shape::impl_ptr arg;
   double size;
 };
 
-class Translate_impl :  public Shape_impl
+class Translate_impl :  public ShapeOpsBase_impl
 {
 public:
   Translate_impl(const Shape::impl_ptr&, const Vector&);
   virtual ~Translate_impl();
-  virtual bool isInside(const Vector& point);
+protected:
+  virtual Vector translateIn(const Vector &point);
+  virtual Vector translateOutPos(const Vector &point);
+  virtual Vector translateOutNormal(const Vector &point);
 private:
-  Shape::impl_ptr arg;
   Vector offset;
 };
 
-class Stretch_impl : public Shape_impl
+class Stretch_impl : public ShapeOpsBase_impl
 {
 public:
   Stretch_impl(const Shape::impl_ptr&, const Vector&);
   virtual ~Stretch_impl();
-  virtual bool isInside(const Vector& point);
+protected:
+  virtual Vector translateIn(const Vector &point);
+  virtual Vector translateOutPos(const Vector &point);
+  virtual Vector translateOutNormal(const Vector &point);
 private:
-  Shape::impl_ptr arg;
   Vector factors;
   Vector reciFactors;
 };
 
-class Transform_impl : public Shape_impl
+class Transform_impl : public ShapeOpsBase_impl
 {
 public:
   Transform_impl(const Shape::impl_ptr&, const Matrix&);
   virtual ~Transform_impl();
-  virtual bool isInside(const Vector& point);
+protected:
+  virtual Vector translateIn(const Vector &point);
+  virtual Vector translateOutPos(const Vector &point);
+  virtual Vector translateOutNormal(const Vector &point);
 private:
-  Shape::impl_ptr arg;
   Matrix transformation;
   Matrix transformationInv;

src/Shapes/Shape_impl.hpp

@@ -10 +10 @@
 
 #include "Shapes/Shape.hpp"
+#include "vector.hpp"
 
 class Shape_impl {
@@ -16 +17 @@
   virtual ~Shape_impl(){};
   virtual bool isInside(const Vector &point)=0;
+  virtual bool isOnSurface(const Vector &point)=0;
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException)=0;
 };
 
@@ -23 +26 @@
     return true;
   }
+  virtual bool isOnSurface(const Vector &point){
+    return false;
+  }
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
 };
 
@@ -28 +37 @@
   virtual bool isInside(const Vector &point){
     return false;
+  }
+  virtual bool isOnSurface(const Vector &point){
+    return false;
+  }
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
   }
 };
@@ -36 +51 @@
   virtual ~AndShape_impl();
   virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException);
 private:
   Shape::impl_ptr lhs;
@@ -46 +63 @@
   virtual ~OrShape_impl();
   virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException);
 private:
   Shape::impl_ptr lhs;
@@ -56 +75 @@
   virtual ~NotShape_impl();
   virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException);
 private:
   Shape::impl_ptr arg;