Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • src/tesselation.cpp

    raccebe r112b09  
    55 *      Author: heber
    66 */
     7
     8#include "Helpers/MemDebug.hpp"
    79
    810#include <fstream>
     
    1719#include "triangleintersectionlist.hpp"
    1820#include "vector.hpp"
     21#include "Line.hpp"
    1922#include "vector_ops.hpp"
    2023#include "verbose.hpp"
     
    441444
    442445  try {
    443     *Intersection = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(*MolCenter, *x);
    444   }
    445   catch (LinearDependenceException &excp) {
    446     Log() << Verbose(1) << excp;
    447     DoeLog(1) && (eLog() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl);
    448     return false;
    449   }
    450 
    451   DoLog(1) && (Log() << Verbose(1) << "INFO: Triangle is " << *this << "." << endl);
    452   DoLog(1) && (Log() << Verbose(1) << "INFO: Line is from " << *MolCenter << " to " << *x << "." << endl);
    453   DoLog(1) && (Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl);
    454 
    455   if (Intersection->DistanceSquared(*endpoints[0]->node->node) < MYEPSILON) {
    456     DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with first endpoint." << endl);
    457     return true;
    458   }   else if (Intersection->DistanceSquared(*endpoints[1]->node->node) < MYEPSILON) {
    459     DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with second endpoint." << endl);
    460     return true;
    461   }   else if (Intersection->DistanceSquared(*endpoints[2]->node->node) < MYEPSILON) {
    462     DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with third endpoint." << endl);
    463     return true;
    464   }
    465   // Calculate cross point between one baseline and the line from the third endpoint to intersection
    466   int i = 0;
    467   do {
    468     try {
    469       CrossPoint = GetIntersectionOfTwoLinesOnPlane(*(endpoints[i%3]->node->node),
    470                                                     *(endpoints[(i+1)%3]->node->node),
    471                                                     *(endpoints[(i+2)%3]->node->node),
    472                                                     *Intersection);
     446    Line centerLine = makeLineThrough(*MolCenter, *x);
     447    *Intersection = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(centerLine);
     448
     449    DoLog(1) && (Log() << Verbose(1) << "INFO: Triangle is " << *this << "." << endl);
     450    DoLog(1) && (Log() << Verbose(1) << "INFO: Line is from " << *MolCenter << " to " << *x << "." << endl);
     451    DoLog(1) && (Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl);
     452
     453    if (Intersection->DistanceSquared(*endpoints[0]->node->node) < MYEPSILON) {
     454      DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with first endpoint." << endl);
     455      return true;
     456    }   else if (Intersection->DistanceSquared(*endpoints[1]->node->node) < MYEPSILON) {
     457      DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with second endpoint." << endl);
     458      return true;
     459    }   else if (Intersection->DistanceSquared(*endpoints[2]->node->node) < MYEPSILON) {
     460      DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with third endpoint." << endl);
     461      return true;
     462    }
     463    // Calculate cross point between one baseline and the line from the third endpoint to intersection
     464    int i = 0;
     465    do {
     466      Line line1 = makeLineThrough(*(endpoints[i%3]->node->node),*(endpoints[(i+1)%3]->node->node));
     467      Line line2 = makeLineThrough(*(endpoints[(i+2)%3]->node->node),*Intersection);
     468      CrossPoint = line1.getIntersection(line2);
    473469      helper = (*endpoints[(i+1)%3]->node->node) - (*endpoints[i%3]->node->node);
    474470      CrossPoint -= (*endpoints[i%3]->node->node);  // cross point was returned as absolute vector
     
    477473      if ((s < -MYEPSILON) || ((s-1.) > MYEPSILON)) {
    478474        DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << "outside of triangle." << endl);
    479         i=4;
    480         break;
     475        return false;
    481476      }
    482477      i++;
    483     } catch (LinearDependenceException &excp){
    484       break;
    485     }
    486   } while (i < 3);
    487   if (i == 3) {
     478    } while (i < 3);
    488479    DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " inside of triangle." << endl);
    489480    return true;
    490   } else {
    491     DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " outside of triangle." << endl);
     481  }
     482  catch (MathException &excp) {
     483    Log() << Verbose(1) << excp;
     484    DoeLog(1) && (eLog() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl);
    492485    return false;
    493486  }
     
    516509  GetCenter(&Direction);
    517510  try {
    518     *ClosestPoint = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(*x, Direction);
    519   }
    520   catch (LinearDependenceException &excp) {
     511    Line l = makeLineThrough(*x, Direction);
     512    *ClosestPoint = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(l);
     513  }
     514  catch (MathException &excp) {
    521515    (*ClosestPoint) = (*x);
    522516  }
     
    541535    Direction = (*endpoints[(i+1)%3]->node->node) - (*endpoints[i%3]->node->node);
    542536    // calculate intersection, line can never be parallel to Direction (is the same vector as PlaneNormal);
    543     CrossPoint[i] = Plane(Direction, InPlane).GetIntersection(*(endpoints[i%3]->node->node), *(endpoints[(i+1)%3]->node->node));
     537    Line l = makeLineThrough(*(endpoints[i%3]->node->node), *(endpoints[(i+1)%3]->node->node));
     538    CrossPoint[i] = Plane(Direction, InPlane).GetIntersection(l);
    544539    CrossDirection[i] = CrossPoint[i] - InPlane;
    545540    CrossPoint[i] -= (*endpoints[i%3]->node->node);  // cross point was returned as absolute vector
Note: See TracChangeset for help on using the changeset viewer.