Changeset 97498a for src


Ignore:
Timestamp:
Jan 8, 2010, 2:04:22 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:
9fb860
Parents:
c15ca2
Message:

Attempt to fix the tesselation::IsInnerPoint().

We try the IsInnerPoint() as follows:

  1. Find nearest BoundaryPoints - working
  2. Find Closest BoundaryLine's - working
  3. Find closest Triangle that is well aligned (wrt to NormalVector and Distance) - unsure whether correctly working
  4. Check whether alignment is on boundary or inside/outside - working
  5. If on boundary, we check whether it's inside of triangle by intersecting with boundary lines - not working

Hence, we code a wrapper for GSL routines, to - finally - allow for solution of linear system of equations.

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

Location:
src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • src/tesselation.cpp

    rc15ca2 r97498a  
    388388};
    389389
    390 /** Finds the point on the triangle \a *BTS the line defined by \a *MolCenter and \a *x crosses through.
     390/** Finds the point on the triangle \a *BTS through which the line defined by \a *MolCenter and \a *x crosses.
    391391 * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
    392392 * This we test if it's really on the plane and whether it's inside the triangle on the plane or not.
    … …  
    411411  }
    412412
     413  Log() << Verbose(1) << "INFO: Triangle is " << *this << "." << endl;
     414  Log() << Verbose(1) << "INFO: Line is from " << *MolCenter << " to " << *x << "." << endl;
     415  Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl;
     416
    413417  // Calculate cross point between one baseline and the line from the third endpoint to intersection
    414418  int i=0;
    415419  do {
    416420    if (CrossPoint.GetIntersectionOfTwoLinesOnPlane(endpoints[i%3]->node->node, endpoints[(i+1)%3]->node->node, endpoints[(i+2)%3]->node->node, Intersection, &NormalVector)) {
     421      CrossPoint.SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
    417422      helper.CopyVector(endpoints[(i+1)%3]->node->node);
    418423      helper.SubtractVector(endpoints[i%3]->node->node);
     424      break;
    419425    } else
    420426      i++;
    421     if (i>2)
    422       break;
    423   } while (CrossPoint.NormSquared() < MYEPSILON);
     427  } while (i<3);
    424428  if (i==3) {
    425429    eLog() << Verbose(0) << "Could not find any cross points, something's utterly wrong here!" << endl;
    426430  }
    427   CrossPoint.SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
     431  Log() << Verbose(1) << "INFO: Crosspoint is " << CrossPoint << "." << endl;
    428432
    429433  // check whether intersection is inside or not by comparing length of intersection and length of cross point
    … …  
    29872991 * \return map of BoundaryPointSet of closest points sorted by squared distance or NULL.
    29882992 */
    2989 DistanceMap * Tesselation::FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const
     2993DistanceToPointMap * Tesselation::FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const
    29902994{
    29912995  Info FunctionInfo(__func__);
    … …  
    30043008  Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
    30053009
    3006   DistanceMap * points = new DistanceMap;
     3010  DistanceToPointMap * points = new DistanceToPointMap;
    30073011  LC->GetNeighbourBounds(Nlower, Nupper);
    30083012  //Log() << Verbose(1) << endl;
    … …  
    30153019          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
    30163020            FindPoint = PointsOnBoundary.find((*Runner)->nr);
    3017             if (FindPoint != PointsOnBoundary.end())
    3018               points->insert(DistancePair (FindPoint->second->node->node->DistanceSquared(x), FindPoint->second) );
     3021            if (FindPoint != PointsOnBoundary.end()) {
     3022              points->insert(DistanceToPointPair (FindPoint->second->node->node->DistanceSquared(x), FindPoint->second) );
     3023              Log() << Verbose(1) << "INFO: Putting " << *FindPoint->second << " into the list." << endl;
     3024            }
    30193025          }
    30203026        } else {
    … …  
    30423048
    30433049  // get closest points
    3044   DistanceMap * points = FindClosestBoundaryPointsToVector(x,LC);
     3050  DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x,LC);
    30453051  if (points == NULL) {
    30463052    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
    … …  
    30553061  Vector Center;
    30563062  Vector BaseLine;
    3057   for (DistanceMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
     3063  for (DistanceToPointMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
    30583064    for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
    30593065      // calculate closest point on line to desired point
    … …  
    31093115
    31103116        // get closest points
    3111         DistanceMap * points = FindClosestBoundaryPointsToVector(x,LC);
     3117        DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x,LC);
    31123118  if (points == NULL) {
    31133119    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
    … …  
    31183124  Log() << Verbose(1) << "Finding closest BoundaryTriangle to " << *x << " ... " << endl;
    31193125  BoundaryLineSet *ClosestLine = NULL;
    3120   double MinDistance = -1.;
    3121   Vector helper;
     3126  double MinDistance = 1e+16;
     3127  Vector BaseLineIntersection;
    31223128  Vector Center;
    31233129  Vector BaseLine;
    3124   for (DistanceMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
     3130  Vector BaseLineCenter;
     3131  for (DistanceToPointMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
    31253132    for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
    3126       for (TriangleMap::iterator TriangleRunner = LineRunner->second->triangles.begin(); TriangleRunner != LineRunner->second->triangles.end(); TriangleRunner++) {
     3133
     3134      BaseLine.CopyVector((LineRunner->second)->endpoints[0]->node->node);
     3135      BaseLine.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
     3136      const double lengthBase = BaseLine.NormSquared();
     3137
     3138      BaseLineIntersection.CopyVector(x);
     3139      BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[0]->node->node);
     3140      const double lengthEndA = BaseLineIntersection.NormSquared();
     3141
     3142      BaseLineIntersection.CopyVector(x);
     3143      BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
     3144      const double lengthEndB = BaseLineIntersection.NormSquared();
     3145
     3146      if ((lengthEndA > lengthBase) || (lengthEndB > lengthBase) || ((lengthEndA < MYEPSILON) || (lengthEndB < MYEPSILON))) {  // intersection would be outside, take closer endpoint
     3147        if ((lengthEndA < MinDistance) && (lengthEndA < lengthEndB)) {
     3148          ClosestLine = LineRunner->second;
     3149          MinDistance = lengthEndA;
     3150          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[0]->node << " is closer with " << lengthEndA << "." << endl;
     3151        } else if (lengthEndB < MinDistance) {
     3152          ClosestLine = LineRunner->second;
     3153          MinDistance = lengthEndB;
     3154          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[1]->node << " is closer with " << lengthEndB << "." << endl;
     3155        } else {
     3156          Log() << Verbose(1) << "REJECT: Line " << *LineRunner->second << " to either endpoints is further away than present closest line candidate: " << lengthEndA << ", " << lengthEndB << ", and distance is longer than baseline:" << lengthBase << "." << endl;
     3157        }
     3158      } else { // intersection is closer, calculate
    31273159        // calculate closest point on line to desired point
    3128         helper.CopyVector((LineRunner->second)->endpoints[0]->node->node);
    3129         helper.AddVector((LineRunner->second)->endpoints[1]->node->node);
    3130         helper.Scale(0.5);
    3131         Center.CopyVector(x);
    3132         Center.SubtractVector(&helper);
    3133         BaseLine.CopyVector((LineRunner->second)->endpoints[0]->node->node);
    3134         BaseLine.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
     3160        BaseLineIntersection.CopyVector(x);
     3161        BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
     3162        Center.CopyVector(&BaseLineIntersection);
    31353163        Center.ProjectOntoPlane(&BaseLine);
    3136         const double distance = Center.NormSquared();
     3164        BaseLineIntersection.SubtractVector(&Center);
     3165        const double distance = BaseLineIntersection.NormSquared();
     3166        if (Center.NormSquared() > BaseLine.NormSquared()) {
     3167          eLog() << Verbose(0) << "Algorithmic error: In second case we have intersection outside of baseline!" << endl;
     3168        }
    31373169        if ((ClosestLine == NULL) || (distance < MinDistance)) {
    3138           // additionally calculate intersection on line (whether it's on the line section or not)
    3139           helper.CopyVector(x);
    3140           helper.SubtractVector((LineRunner->second)->endpoints[0]->node->node);
    3141           helper.SubtractVector(&Center);
    3142           const double lengthA = helper.ScalarProduct(&BaseLine);
    3143           helper.CopyVector(x);
    3144           helper.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
    3145           helper.SubtractVector(&Center);
    3146           const double lengthB = helper.ScalarProduct(&BaseLine);
    3147           if (lengthB*lengthA < 0) {  // if have different sign
    3148             ClosestLine = LineRunner->second;
    3149             MinDistance = distance;
    3150             Log() << Verbose(1) << "ACCEPT: New closest line is " << *ClosestLine << " with projected distance " << MinDistance << "." << endl;
    3151           } else {
    3152             Log() << Verbose(1) << "REJECT: Intersection is outside of the line section: " << lengthA << " and " << lengthB << "." << endl;
    3153           }
     3170          ClosestLine = LineRunner->second;
     3171          MinDistance = distance;
     3172          Log() << Verbose(1) << "ACCEPT: Intersection in between endpoints, new closest line " << *LineRunner->second << " is " << *ClosestLine << " with projected distance " << MinDistance << "." << endl;
    31543173        } else {
    3155           Log() << Verbose(1) << "REJECT: Point is too further away than present line: " << distance << " >> " << MinDistance << "." << endl;
     3174          Log() << Verbose(2) << "REJECT: Point is further away from line " << *LineRunner->second << " than present closest line: " << distance << " >> " << MinDistance << "." << endl;
    31563175        }
    31573176      }
    … …  
    31903209    return NULL;
    31913210
    3192   // reject all which are not closest
    3193   double MinDistance = -1.;
     3211  // go through all and pick the one with the best alignment to x
     3212  double MinAlignment = 2.*M_PI;
    31943213  for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++) {
    31953214    (*Runner)->GetCenter(&Center);
    31963215    helper.CopyVector(x);
    31973216    helper.SubtractVector(&Center);
    3198     helper.ProjectOntoPlane(&(*Runner)->NormalVector);
    3199     const double distance = helper.NormSquared();
    3200     if (fabs(distance - MinDistance) < MYEPSILON) {  // degenerate case, keep both
    3201       candidates.push_back(*Runner);
    3202     } else if (distance < MinDistance) {
    3203       candidates.clear();
    3204       candidates.push_back(*Runner);
    3205       MinDistance = distance;
     3217    const double Alignment = helper.Angle(&(*Runner)->NormalVector);
     3218    if (Alignment < MinAlignment) {
     3219      result = *Runner;
     3220      MinAlignment = Alignment;
     3221      Log() << Verbose(1) << "ACCEPT: Triangle " << *result << " is better aligned with " << MinAlignment << "." << endl;
     3222    } else {
     3223      Log() << Verbose(1) << "REJECT: Triangle " << *result << " is worse aligned with " << MinAlignment << "." << endl;
    32063224    }
    32073225  }
    32083226  delete(triangles);
    3209   if (!candidates.empty()) {
    3210     if (candidates.size() == 1) { // there is no degenerate case
    3211       result = candidates.front();
    3212       Log() << Verbose(1) << "Normal Vector of this triangle is " << result->NormalVector << "." << endl;
    3213     } else {
    3214       result = candidates.front();
    3215       result->GetCenter(&Center);
    3216       Center.SubtractVector(x);
    3217       Log() << Verbose(1) << "Normal Vector of this front side is " << result->NormalVector << "." << endl;
    3218       if (Center.ScalarProduct(&result->NormalVector) < 0) {
    3219         result = candidates.back();
    3220         Log() << Verbose(1) << "Normal Vector of this back side is " << result->NormalVector << "." << endl;
    3221         if (Center.ScalarProduct(&result->NormalVector) < 0) {
    3222           eLog() << Verbose(1) << "Front and back side yield NormalVector in wrong direction!" << endl;
    3223         }
    3224       }
    3225     }
    3226   } else {
    3227     result = NULL;
    3228     Log() << Verbose(0) << "No closest triangle found." << endl;
    3229   }
     3227
    32303228  return result;
    32313229};
    32323230
    32333231/** Checks whether the provided Vector is within the tesselation structure.
     3232 * Calls FindClosestTriangleToVector() and checks whether the resulting triangle's BoundaryTriangleSet#NormalVector points
     3233 * towards or away from the given \a &Point.
    32343234 *
    32353235 * @param point of which to check the position
    … …  
    32453245  Vector Center;
    32463246  Vector helper;
     3247  Vector DistanceToCenter;
     3248  Vector Intersection;
    32473249
    32483250  if (result == NULL) {// is boundary point or only point in point cloud?
    … …  
    32553257  result->GetCenter(&Center);
    32563258  Log() << Verbose(2) << "INFO: Central point of the triangle is " << Center << "." << endl;
    3257   Center.SubtractVector(&Point);
     3259  DistanceToCenter.CopyVector(&Center);
     3260  DistanceToCenter.SubtractVector(&Point);
    32583261  Log() << Verbose(2) << "INFO: Vector from point to test to center is " << Center << "." << endl;
    3259   if (Center.ScalarProduct(&result->NormalVector) > -MYEPSILON) {
     3262
     3263  // check whether we are on boundary
     3264  if (fabs(DistanceToCenter.ScalarProduct(&result->NormalVector)) < MYEPSILON) {
     3265    // calculate whether inside of triangle
     3266    DistanceToCenter.ProjectOntoPlane(&result->NormalVector);
     3267    DistanceToCenter.AddVector(&Center);
     3268    Center.CopyVector(&DistanceToCenter);
     3269    Center.SubtractVector(&result->NormalVector); // points towards MolCenter
     3270    DistanceToCenter.AddVector(&result->NormalVector); // points outside
     3271    Log() << Verbose(1) << "INFO: Calling Intersection with " << Center << " and " << DistanceToCenter << "." << endl;
     3272    if (result->GetIntersectionInsideTriangle(&Center, &DistanceToCenter, &Intersection)) {
     3273      Log() << Verbose(1) << Point << " is inner point: sufficiently close to boundary, " << Intersection << "." << endl;
     3274      return true;
     3275    } else {
     3276      Log() << Verbose(1) << Point << " is NOT an inner point: on triangle plane but outside of triangle bounds." << endl;
     3277      return false;
     3278    }
     3279  }
     3280
     3281  // then check direction to boundary
     3282  if (DistanceToCenter.ScalarProduct(&result->NormalVector) > MYEPSILON) {
    32603283    Log() << Verbose(1) << Point << " is an inner point." << endl;
    32613284    return true;
    32623285  } else {
    3263     for (int i=0;i<3;i++) {
    3264       helper.CopyVector(result->endpoints[i]->node->node);
    3265       helper.SubtractVector(&Point);
    3266       if (Center.NormSquared() > helper.NormSquared())
    3267         Center.CopyVector(&helper);
    3268     }
    3269     if (Center.NormSquared() < epsilon*epsilon) {
    3270       Log() << Verbose(1) << Point << " is inner point, being sufficiently close (wrt " << epsilon << ") to boundary." << endl;
    3271       return true;
    3272     }
    32733286    Log() << Verbose(1) << Point << " is NOT an inner point." << endl;
    32743287    return false;
    32753288  }
    3276 }
    3277 
    3278 /** Checks whether the provided TesselPoint is within the tesselation structure.
    3279  *
    3280  * @param *Point of which to check the position
    3281  * @param *LC Linked Cell structure
    3282  * @param epsilon Distance of \a &Point to Center of triangle is tested greater against this value (standard: -MYEPSILON)
    3283  *
    3284  * @return true if the point is inside the tesselation structure, false otherwise
    3285  */
    3286 bool Tesselation::IsInnerPoint(const TesselPoint * const Point, const LinkedCell* const LC, const double epsilon) const
    3287 {
    3288         Info FunctionInfo(__func__);
    3289   return IsInnerPoint(*(Point->node), LC, epsilon);
    32903289}
    32913290
    … …  
    34603459  }
    34613460
     3461  // check whether there's something to do
     3462  if (SetOfNeighbours->size() < 3) {
     3463    for (TesselPointSet::iterator TesselRunner = SetOfNeighbours->begin(); TesselRunner != SetOfNeighbours->end(); TesselRunner++)
     3464      connectedCircle->push_back(*TesselRunner);
     3465    return connectedCircle;
     3466  }
     3467
     3468  Log() << Verbose(1) << "INFO: Point is " << *Point << " and Reference is " << *Reference << "." << endl;
    34623469  // calculate central point
    3463   for (TesselPointSet::const_iterator TesselRunner = SetOfNeighbours->begin(); TesselRunner != SetOfNeighbours->end(); TesselRunner++)
    3464     center.AddVector((*TesselRunner)->node);
     3470
     3471  TesselPointSet::const_iterator TesselA = SetOfNeighbours->begin();
     3472  TesselPointSet::const_iterator TesselB = SetOfNeighbours->begin();
     3473  TesselPointSet::const_iterator TesselC = SetOfNeighbours->begin();
     3474  TesselB++;
     3475  TesselC++;
     3476  TesselC++;
     3477  int counter = 0;
     3478  while (TesselC != SetOfNeighbours->end()) {
     3479    helper.MakeNormalVector((*TesselA)->node, (*TesselB)->node, (*TesselC)->node);
     3480    Log() << Verbose(0) << "Making normal vector out of " << *(*TesselA) << ", " << *(*TesselB) << " and " << *(*TesselC) << ":" << helper << endl;
     3481    counter++;
     3482    TesselA++;
     3483    TesselB++;
     3484    TesselC++;
     3485    PlaneNormal.AddVector(&helper);
     3486  }
    34653487  //Log() << Verbose(0) << "Summed vectors " << center << "; number of points " << connectedPoints.size()
    3466   //  << "; scale factor " << 1.0/connectedPoints.size();
    3467   center.Scale(1.0/SetOfNeighbours->size());
    3468   Log() << Verbose(1) << "INFO: Calculated center of all circle points is " << center << "." << endl;
    3469 
    3470   // projection plane of the circle is at the closes Point and normal is pointing away from center of all circle points
    3471   PlaneNormal.CopyVector(Point->node);
    3472   PlaneNormal.SubtractVector(&center);
    3473   PlaneNormal.Normalize();
     3488  //  << "; scale factor " << counter;
     3489  PlaneNormal.Scale(1.0/(double)counter);
     3490//  Log() << Verbose(1) << "INFO: Calculated center of all circle points is " << center << "." << endl;
     3491//
     3492//  // projection plane of the circle is at the closes Point and normal is pointing away from center of all circle points
     3493//  PlaneNormal.CopyVector(Point->node);
     3494//  PlaneNormal.SubtractVector(&center);
     3495//  PlaneNormal.Normalize();
    34743496  Log() << Verbose(1) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl;
    34753497
    … …  
    35043526    helper.ProjectOntoPlane(&PlaneNormal);
    35053527    double angle = GetAngle(helper, AngleZero, OrthogonalVector);
     3528    if (angle > M_PI) // the correction is of no use here (and not desired)
     3529      angle = 2.*M_PI - angle;
    35063530    Log() << Verbose(0) << "INFO: Calculated angle between " << helper << " and " << AngleZero << " is " << angle << " for point " << **listRunner << "." << endl;
    35073531    InserterTest = anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));

  • src/tesselation.hpp

    rc15ca2 r97498a  
    7979#define PolygonList list < class BoundaryPolygonSet * >
    8080
    81 #define DistanceMap multimap <double, class BoundaryPointSet * >
    82 #define DistancePair pair <double, class BoundaryPointSet * >
     81#define DistanceToPointMap multimap <double, class BoundaryPointSet * >
     82#define DistanceToPointPair pair <double, class BoundaryPointSet * >
    8383
    8484#define DistanceMultiMap multimap <double, pair < PointMap::iterator, PointMap::iterator> >
    … …  
    313313    class BoundaryTriangleSet * FindClosestTriangleToPoint(const Vector *x, const LinkedCell* LC) const;
    314314    bool IsInnerPoint(const Vector &Point, const LinkedCell* const LC, const double epsilon = -MYEPSILON) const;
    315     bool IsInnerPoint(const TesselPoint * const Point, const LinkedCell* const LC, const double epsilon = -MYEPSILON) const;
    316315    bool AddBoundaryPoint(TesselPoint * Walker, const int n);
    317     DistanceMap * FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const;
     316    DistanceToPointMap * FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const;
    318317    BoundaryLineSet * FindClosestBoundaryLineToVector(const Vector *x, const LinkedCell* LC) const;
    319318    TriangleList * FindClosestTrianglesToVector(const Vector *x, const LinkedCell* LC) const;

  • src/unittests/tesselation_insideoutsideunittest.cpp

    rc15ca2 r97498a  
    2626void TesselationInOutsideTest::setUp()
    2727{
    28   setVerbosity(3);
     28  setVerbosity(5);
    2929
    3030  // create corners
    … …  
    3232  Walker = new TesselPoint;
    3333  Walker->node = new Vector(0., 0., 0.);
    34   Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
     34  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
    3535  strcpy(Walker->Name, "1");
    3636  Walker->nr = 1;
    … …  
    3838  Walker = new TesselPoint;
    3939  Walker->node = new Vector(0., 1., 0.);
    40   Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
     40  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
    4141  strcpy(Walker->Name, "2");
    4242  Walker->nr = 2;
    … …  
    4444  Walker = new TesselPoint;
    4545  Walker->node = new Vector(1., 0., 0.);
    46   Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
     46  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
    4747  strcpy(Walker->Name, "3");
    4848  Walker->nr = 3;
    … …  
    5050  Walker = new TesselPoint;
    5151  Walker->node = new Vector(1., 1., 0.);
    52   Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
     52  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
    5353  strcpy(Walker->Name, "4");
    5454  Walker->nr = 4;
    … …  
    5656  Walker = new TesselPoint;
    5757  Walker->node = new Vector(0., 0., 1.);
    58   Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
     58  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
    5959  strcpy(Walker->Name, "5");
    6060  Walker->nr = 5;
    … …  
    6262  Walker = new TesselPoint;
    6363  Walker->node = new Vector(0., 1., 1.);
    64   Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
     64  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
    6565  strcpy(Walker->Name, "6");
    6666  Walker->nr = 6;
    … …  
    6868  Walker = new TesselPoint;
    6969  Walker->node = new Vector(1., 0., 1.);
    70   Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
     70  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
    7171  strcpy(Walker->Name, "7");
    7272  Walker->nr = 7;
    … …  
    7474  Walker = new TesselPoint;
    7575  Walker->node = new Vector(1., 1., 1.);
    76   Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
     76  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
    7777  strcpy(Walker->Name, "8");
    7878  Walker->nr = 8;
    … …  
    146146{
    147147  double n[3];
    148   const double boundary = 5.;
     148  const double boundary = 2.;
    149149  const double step = 1.;
    150150
    … …  
    153153    for (n[1] = -boundary; n[1] <= boundary; n[1]+=step)
    154154      for (n[2] = -boundary; n[2] <= boundary; n[2]+=step) {
    155         if ( ((n[0] > 0) && (n[0] > 0) && (n[0] > 0)) && (n[0]+n[1]+n[2] <= 1.))
     155        if ( ((n[0] >= 0.) && (n[1] >= 0.) && (n[2] >= 0.)) && (fabs(n[0])+fabs(n[1])+fabs(n[2]) <= 1.))
    156156          CPPUNIT_ASSERT_EQUAL( true , TesselStruct->IsInnerPoint(Vector(n[0], n[1], n[2]), LinkedList) );
    157157        else

  • src/vector.cpp

    rc15ca2 r97498a  
    88#include "defs.hpp"
    99#include "helpers.hpp"
    10 #include "memoryallocator.hpp"
     10#include "info.hpp"
    1111#include "leastsquaremin.hpp"
    1212#include "log.hpp"
     13#include "memoryallocator.hpp"
    1314#include "vector.hpp"
    1415#include "verbose.hpp"
     16
     17#include <gsl/gsl_linalg.h>
     18#include <gsl/gsl_matrix.h>
     19#include <gsl/gsl_permutation.h>
     20#include <gsl/gsl_vector.h>
    1521
    1622/************************************ Functions for class vector ************************************/
    … …  
    219225bool Vector::GetIntersectionWithPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset, const Vector * const Origin, const Vector * const LineVector)
    220226{
     227  Info FunctionInfo(__func__);
    221228  double factor;
    222229  Vector Direction, helper;
    … …  
    226233  Direction.SubtractVector(Origin);
    227234  Direction.Normalize();
    228   //Log() << Verbose(4) << "INFO: Direction is " << Direction << "." << endl;
     235  Log() << Verbose(1) << "INFO: Direction is " << Direction << "." << endl;
    229236  factor = Direction.ScalarProduct(PlaneNormal);
    230237  if (factor < MYEPSILON) { // Uniqueness: line parallel to plane?
    … …  
    236243  factor = helper.ScalarProduct(PlaneNormal)/factor;
    237244  if (factor < MYEPSILON) { // Origin is in-plane
    238     //Log() << Verbose(2) << "Origin of line is in-plane, simple." << endl;
     245    Log() << Verbose(1) << "Origin of line is in-plane, simple." << endl;
    239246    CopyVector(Origin);
    240247    return true;
    … …  
    243250  Direction.Scale(factor);
    244251  CopyVector(Origin);
    245   //Log() << Verbose(4) << "INFO: Scaled direction is " << Direction << "." << endl;
     252  Log() << Verbose(1) << "INFO: Scaled direction is " << Direction << "." << endl;
    246253  AddVector(&Direction);
    247254
    … …  
    250257  helper.SubtractVector(PlaneOffset);
    251258  if (helper.ScalarProduct(PlaneNormal) < MYEPSILON) {
    252     //Log() << Verbose(2) << "INFO: Intersection at " << *this << " is good." << endl;
     259    Log() << Verbose(1) << "INFO: Intersection at " << *this << " is good." << endl;
    253260    return true;
    254261  } else {
    … …  
    299306bool Vector::GetIntersectionOfTwoLinesOnPlane(const Vector * const Line1a, const Vector * const Line1b, const Vector * const Line2a, const Vector * const Line2b, const Vector *PlaneNormal)
    300307{
    301   bool result = true;
     308  Info FunctionInfo(__func__);
    302309  Vector Direction, OtherDirection;
    303   Vector AuxiliaryNormal;
    304   Vector Distance;
    305   const Vector *Normal = NULL;
    306   Vector *ConstructedNormal = NULL;
    307   bool FreeNormal = false;
     310  gsl_matrix *A = gsl_matrix_alloc(NDIM,NDIM);
     311  gsl_vector *b = gsl_vector_alloc(NDIM);
     312  gsl_vector *x = gsl_vector_alloc(NDIM);
     313  gsl_permutation *perm = NULL;
     314  int signum;
    308315
    309316  // construct both direction vectors
    310   Zero();
    311317  Direction.CopyVector(Line1b);
    312318  Direction.SubtractVector(Line1a);
    … …  
    318324    return false;
    319325
    320   Direction.Normalize();
    321   OtherDirection.Normalize();
    322 
    323   //Log() << Verbose(4) << "INFO: Normalized Direction " << Direction << " and OtherDirection " << OtherDirection << "." << endl;
    324 
    325   if (fabs(OtherDirection.ScalarProduct(&Direction) - 1.) < MYEPSILON) { // lines are parallel
    326     if ((Line1a == Line2a) || (Line1a == Line2b))
    327       CopyVector(Line1a);
    328     else if ((Line1b == Line2b) || (Line1b == Line2b))
    329         CopyVector(Line1b);
    330     else
    331       return false;
    332     Log() << Verbose(4) << "INFO: Intersection is " << *this << "." << endl;
     326  // set vector
     327  for (int i=0;i<NDIM;i++)
     328    gsl_vector_set(b, i, Line1a->x[i]-Line2a->x[i]);
     329  Log() << Verbose(1) << "b = " << endl;
     330  gsl_vector_fprintf(stdout, b, "%g");
     331
     332  // set matrix
     333  for (int i=0;i<NDIM;i++)
     334    gsl_matrix_set(A, 0, i, -Direction.x[i]);
     335  for (int i=0;i<NDIM;i++)
     336    gsl_matrix_set(A, 1, i, OtherDirection.x[i]);
     337  for (int i=0;i<NDIM;i++)
     338    gsl_matrix_set(A, 2, i, 1.);
     339  Log() << Verbose(1) << "A = " << endl;
     340  gsl_matrix_fprintf(stdout, A, "%g");
     341
     342  // Solve A x=b for x
     343  perm = gsl_permutation_alloc(NDIM);
     344  gsl_linalg_LU_decomp(A, perm, &signum);
     345  gsl_linalg_LU_solve(A, perm, b, x);
     346  gsl_permutation_free(perm);
     347  gsl_vector_free(b);
     348  gsl_matrix_free(A);
     349
     350  Log() << Verbose(1) << "Solution is " << gsl_vector_get(x,0) << ", " << gsl_vector_get(x,1) << "." << endl;
     351
     352  CopyVector(&Direction);
     353  Scale(gsl_vector_get(x,0));
     354  AddVector(Line1a);
     355  Log() << Verbose(1) << "INFO: First intersection is " << *this << "." << endl;
     356
     357  Vector test;
     358  test.CopyVector(&OtherDirection);
     359  test.Scale(gsl_vector_get(x,1));
     360  test.AddVector(Line2a);
     361  Log() << Verbose(1) << "INFO: Second intersection is " << test << "." << endl;
     362  test.SubtractVector(this);
     363
     364  gsl_vector_free(x);
     365
     366  if (test.IsZero())
    333367    return true;
    334   } else {
    335     // check whether we have a plane normal vector
    336     if (PlaneNormal == NULL) {
    337       ConstructedNormal = new Vector;
    338       ConstructedNormal->MakeNormalVector(&Direction, &OtherDirection);
    339       Normal = ConstructedNormal;
    340       FreeNormal = true;
    341     } else
    342       Normal = PlaneNormal;
    343 
    344     AuxiliaryNormal.MakeNormalVector(&OtherDirection, Normal);
    345     //Log() << Verbose(4) << "INFO: PlaneNormal is " << *Normal << " and AuxiliaryNormal " << AuxiliaryNormal << "." << endl;
    346 
    347     Distance.CopyVector(Line2a);
    348     Distance.SubtractVector(Line1a);
    349     //Log() << Verbose(4) << "INFO: Distance is " << Distance << "." << endl;
    350     if (Distance.IsZero()) {
    351       // offsets are equal, match found
    352       CopyVector(Line1a);
    353       result = true;
    354     } else {
    355       CopyVector(Distance.Projection(&AuxiliaryNormal));
    356       //Log() << Verbose(4) << "INFO: Projected Distance is " << *this << "." << endl;
    357       double factor = Direction.ScalarProduct(&AuxiliaryNormal);
    358       //Log() << Verbose(4) << "INFO: Scaling factor is " << factor << "." << endl;
    359       Scale(1./(factor*factor));
    360       //Log() << Verbose(4) << "INFO: Scaled Distance is " << *this << "." << endl;
    361       CopyVector(Projection(&Direction));
    362       //Log() << Verbose(4) << "INFO: Distance, projected into Direction, is " << *this << "." << endl;
    363       if (this->IsZero())
    364         result = false;
    365       else
    366         result = true;
    367       AddVector(Line1a);
    368     }
    369 
    370     if (FreeNormal)
    371       delete(ConstructedNormal);
    372   }
    373   if (result)
    374     Log() << Verbose(4) << "INFO: Intersection is " << *this << "." << endl;
    375 
    376   return result;
     368  else
     369    return false;
    377370};
    378371
Note: See TracChangeset for help on using the changeset viewer.