Diff [3b179825d4e6f17bd22fbd7f880f59cc7fd8924e:0516fd9f440d63cdfbb63b946c9933f28f304b68] for / – MoleCuilder

LinearAlgebra/src/LinearAlgebra/Line.cpp

-              r3b1798
+              r0516fd
  */
 Vector Line::getIntersection(const Line& otherLine) const{
-  Info FunctionInfo(__func__);
   pointset line1Points = getPointsOnLine();
 …
   //}
   if (fabs(M->Determinant()) > LINALG_MYEPSILON()) {
     Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl;
+    eLog() << Verbose(2) << "Determinant of coefficient matrix is NOT zero." << endl;
     throw SkewException() << LinearAlgebraDeterminant(M->Determinant()) << LinearAlgebraMatrixContent(&(*M));
+  }
   Log() << Verbose(1) << "INFO: Line1a = " << Line1a << ", Line1b = " << Line1b << ", Line2a = " << Line2a << ", Line2b = " << Line2b << "." << endl;
+  Log() << Verbose(6) << "INFO: Line1a = " << Line1a << ", Line1b = " << Line1b << ", Line2a = " << Line2a << ", Line2b = " << Line2b << "." << endl;
 …
   Vector c = Line2a - Line1a;
   Vector d = Line2b - Line1b;
   Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl;
+  Log() << Verbose(7) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl;
   if ((a.NormSquared() <= LINALG_MYEPSILON()) || (b.NormSquared() <= LINALG_MYEPSILON())) {
    res.Zero();
    Log() << Verbose(1) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl;
+   eLog() << Verbose(2) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl;
    throw LinearDependenceException() << LinearAlgebraVectorPair( make_pair(&a, &b) );
+  }
 …
     if ((factor > -LINALG_MYEPSILON()) && (factor - 1. <= LINALG_MYEPSILON())) {
       res = Line2a;
       Log() << Verbose(1) << "Lines conincide." << endl;
+      Log() << Verbose(5) << "Lines conincide." << endl;
       return res;
     } else {
 …
       if ((factor > -LINALG_MYEPSILON()) && (factor - 1. <= LINALG_MYEPSILON())) {
         res = Line2b;
         Log() << Verbose(1) << "Lines conincide." << endl;
+        Log() << Verbose(5) << "Lines conincide." << endl;
         return res;
+      }
+    }
     Log() << Verbose(1) << "Lines are parallel." << endl;
+    eLog() << Verbose(2) << "Lines are parallel." << endl;
     res.Zero();
     throw LinearDependenceException() << LinearAlgebraVectorPair( make_pair(&a, &b) );
 …
   temp2 = a;
   temp2.VectorProduct(b);
   Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl;
+  Log() << Verbose(7) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl;
   if (fabs(temp2.NormSquared()) > LINALG_MYEPSILON())
     s = temp1.ScalarProduct(temp2)/temp2.NormSquared();
   else
     s = 0.;
   Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl;
+  Log() << Verbose(6) << "Factor s is " << temp1.ScalarProduct(temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl;
   // construct intersection
 …
   res.Scale(s);
   res += Line1a;
   Log() << Verbose(1) << "Intersection is at " << res << "." << endl;
+  Log() << Verbose(5) << "Intersection is at " << res << "." << endl;
   return res;

configure.ac

r3b1798	r0516fd
207	207	BOOST_ITERATOR
208	208	BOOST_MULTIARRAY
	209	BOOST_MULTIINDEXCONTAINER
209	210	BOOST_MPL
210	211	BOOST_PREPROCESSOR

doc/userguide/userguide.xml

-              r3b1798
+              r0516fd
           Non-convex envelope</title>
+          <para>This will create a non-convex envelope for a molecule.</para>
+          <para>This will create a non-convex envelope for a molecule and store
+          it to a file for viewing with external programs.</para>
           <programlisting>
 …
           <para>This tesselation file can be conveniently viewed with
           <productname>TecPlot</productname> or with one of the Tcl script
+          in the util folder with <productname>VMD</productname>.</para>
+          in the util folder with <productname>VMD</productname>. Also,
+          still pictures can be produced with <productname>Raster3D
+          </productname>.
+          <note>The required file header.r3d can be found in a subfolder of
+          the util folder.</note>
+          </para>
         </section>
 …
           envelope</title>
+          <para>This will create a convex envelope for a molecule.</para>
+          <para>This will create a convex envelope for a molecule and give the
+          volumes of both the non-convex and the convex envelope. This is good
+          for measuring the space a molecule takes up, e.g. when filling a
+          domain and taking care of correct densities.</para>
           <programlisting>
 …
           </programlisting>
           <para>This tesselation file can be conveniently viewed with
+          <para>This tesselation file can be likewise viewed with
           <productname>TecPlot</productname> or with one of the Tcl script
           in the util folder with <productname>VMD</productname>.</para>

m4/boost.m4

-              r3b1798
+              r0516fd
 [BOOST_FIND_HEADER([boost/multi_array.hpp])])
+# BOOST_MULTIINDEXCCONTAINER()
+# ------------------
+# Look for Boost.MultiIndexContainer
+BOOST_DEFUN([MultiIndexContainer],
+[BOOST_FIND_HEADER([boost/multi_index_container.hpp])])
 # BOOST_NUMERIC_UBLAS()

src/Actions/ActionQueue.cpp

-              r3b1798
+              r0516fd
 void ActionQueue::queueAction(const Action * const _action, enum Action::QueryOptions state)
+{
-  OBSERVE;
-  NOTIFY(ActionQueued);
   Action *newaction = _action->clone(state);
   newaction->prepare(state);
 …
     std::cerr << "ActionQueue cleared." << std::endl;
+  }
+  if (lastActionOk) {
+    OBSERVE;
+    NOTIFY(ActionQueued);
+    _lastchangedaction = newaction;
+  }
 #else
+  {
 …
   mtx_queue.unlock();
 #endif
-  _lastchangedaction = newaction;
+}
 …
         CurrentAction = (size_t)-1;
+      }
+      if (lastActionOk) {
+        OBSERVE;
+        NOTIFY(ActionQueued);
+        _lastchangedaction = actionqueue[CurrentAction];
+      }
       // access actionqueue, hence using mutex
       mtx_queue.lock();

src/Actions/TesselationAction/ConvexEnvelopeAction.cpp

-              r3b1798
+              r0516fd
     LOG(1, "Storing tecplot non-convex data in " << params.filenameNonConvex.get() << ".");
     PointCloudAdaptor<molecule> cloud(mol, mol->name);
     LCList = new LinkedCell_deprecated(cloud, 100.);
+    LCList = new LinkedCell_deprecated(cloud, 2.*params.SphereRadius.get());
     //Boundaries *BoundaryPoints = NULL;
     //FindConvexBorder(mol, BoundaryPoints, TesselStruct, LCList, argv[argptr]);
     // TODO: Beide Funktionen sollten streams anstelle des Filenamen benutzen, besser fuer unit tests
     Success &= FindNonConvexBorder(mol, TesselStruct, LCList, 50., params.filenameNonConvex.get().string().c_str());
+    Success &= FindNonConvexBorder(mol, TesselStruct, LCList, params.SphereRadius.get(), params.filenameNonConvex.get().string().c_str());
     //RemoveAllBoundaryPoints(TesselStruct, mol, argv[argptr]);
+    const double volumedifference = ConvexizeNonconvexEnvelope(TesselStruct, mol, params.filenameConvex.get().string().c_str());
+    const double volumedifference =
+        ConvexizeNonconvexEnvelope(TesselStruct, mol, params.filenameConvex.get().string().c_str(),
+            params.DoOutputEveryStep.get());
+    // check whether tesselated structure is truly convex
+    if (!TesselStruct->isConvex()) {
+      ELOG(1, "Tesselated surface has not been properly convexized.");
+        Success = false;
+    } else {
+      LOG(2, "DEBUG: We do have a convex surface tesselation now.");
+    }
+    // we check whether all molecule's atoms are still inside
+    std::vector<std::string> outside_atoms;
+    for(molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter)
+      if (!TesselStruct->IsInnerPoint((*iter)->getPosition(), LCList))
+        outside_atoms.push_back((*iter)->getName());
+    if (outside_atoms.empty())
+      LOG(2, "DEBUG: All molecule's atoms are inside the tesselated, convex surface.");
+    else {
+      ELOG(1, "The following atoms are not inside the tesselated, convex surface:"
+          << outside_atoms);
+      Success = false;
+    }
     const double clustervolume = TesselStruct->getVolumeOfConvexEnvelope(configuration->GetIsAngstroem());
     LOG(0, "The tesselated volume area is " << clustervolume << " " << (configuration->GetIsAngstroem() ? "angstrom" : "atomiclength") << "^3.");
 …
     return Action::success;
   else {
+    STATUS("Failed to find the non convex border.");
+    STATUS(std::string("Failed to find the non convex border")
+        +std::string("containing all atoms")
+        +std::string("."));
     return Action::failure;
+  }

src/Actions/TesselationAction/ConvexEnvelopeAction.def

-              r3b1798
+              r0516fd
 class TesselationListClass;
+#include "Parameters/Validators/DummyValidator.hpp"
 #include "Parameters/Validators/Ops_Validator.hpp"
+#include "Parameters/Validators/Specific/BoxLengthValidator.hpp"
 #include "Parameters/Validators/Specific/FilePresentValidator.hpp"
 …
 // ValueStorage by the token "Z" -> first column: int, Z, "Z"
 // "undefine" if no parameters are required, use (NOPARAM_DEFAULT) for each (undefined) default value
 #define paramtypes (boost::filesystem::path)(boost::filesystem::path)
 #define paramtokens ("convex-file")("nonconvex-file")
 #define paramdescriptions ("filename of the convex envelope")("filename of the non-convex envelope")
 #undef paramdefaults
 #define paramreferences (filenameConvex)(filenameNonConvex)
+#define paramtypes (double)(boost::filesystem::path)(boost::filesystem::path)(bool)
+#define paramtokens ("convex-envelope")("convex-file")("nonconvex-file")("DoOutputEveryStep")
+#define paramdescriptions ("radius of the rolling sphere")("filename of the convex envelope")("filename of the non-convex envelope")("whether to write an extra file for debugging of each convexization step")
+#define paramdefaults (NOPARAM_DEFAULT)(NOPARAM_DEFAULT)(NOPARAM_DEFAULT)(PARAM_DEFAULT(0))
+#define paramreferences (SphereRadius)(filenameConvex)(filenameNonConvex)(DoOutputEveryStep)
 #define paramvalids \
+(BoxLengthValidator()) \
 (!FilePresentValidator()) \
+(!FilePresentValidator())
+(!FilePresentValidator()) \
+(DummyValidator<bool>())
 #undef statetypes

src/Potentials/Specifics/ConstantPotential.cpp

r3b1798	r0516fd
63	63	;
64	64	const std::string ConstantPotential::potential_token("constant");
65		Coordinator::ptr ConstantPotential::coordinator(~~new OneBody_Constant(~~));
	65	Coordinator::ptr ConstantPotential::coordinator(Memory::ignore(new OneBody_Constant()));
66	66
67	67	ConstantPotential::ConstantPotential() :

src/Potentials/Specifics/FourBodyPotential_Improper.cpp

r3b1798	r0516fd
49	49	;
50	50	const std::string FourBodyPotential_Improper::improper_token("improper");
51		Coordinator::ptr FourBodyPotential_Improper::coordinator(~~new FourBody_ImproperAngle(~~));
	51	Coordinator::ptr FourBodyPotential_Improper::coordinator(Memory::ignore(new FourBody_ImproperAngle()));
52	52
53	53	FourBodyPotential_Improper::FourBodyPotential_Improper() :

src/Potentials/Specifics/FourBodyPotential_Torsion.cpp

r3b1798	r0516fd
64	64	;
65	65	const std::string FourBodyPotential_Torsion::potential_token("torsion");
66		Coordinator::ptr FourBodyPotential_Torsion::coordinator(~~new FourBody_TorsionAngle(~~));
	66	Coordinator::ptr FourBodyPotential_Torsion::coordinator(Memory::ignore(new FourBody_TorsionAngle()));
67	67
68	68	FourBodyPotential_Torsion::FourBodyPotential_Torsion() :

src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp

r3b1798	r0516fd
80	80	;
81	81	const std::string ManyBodyPotential_Tersoff::potential_token("tersoff");
82		Coordinator::ptr ManyBodyPotential_Tersoff::coordinator(~~new OneBody_Constant(~~));
	82	Coordinator::ptr ManyBodyPotential_Tersoff::coordinator(Memory::ignore(new OneBody_Constant()));
83	83
84	84	ManyBodyPotential_Tersoff::ManyBodyPotential_Tersoff() :

src/Potentials/Specifics/PairPotential_Harmonic.cpp

r3b1798	r0516fd
65	65	;
66	66	const std::string PairPotential_Harmonic::potential_token("harmonic_bond");
67		Coordinator::ptr PairPotential_Harmonic::coordinator(~~new TwoBody_Length(~~));
	67	Coordinator::ptr PairPotential_Harmonic::coordinator(Memory::ignore(new TwoBody_Length()));
68	68
69	69	PairPotential_Harmonic::PairPotential_Harmonic() :

src/Potentials/Specifics/PairPotential_LennardJones.cpp

r3b1798	r0516fd
65	65	;
66	66	const std::string PairPotential_LennardJones::potential_token("lennardjones");
67		Coordinator::ptr PairPotential_LennardJones::coordinator(~~new TwoBody_Length(~~));
	67	Coordinator::ptr PairPotential_LennardJones::coordinator(Memory::ignore(new TwoBody_Length()));
68	68
69	69	void PairPotential_LennardJones::setDefaultParameters()

src/Potentials/Specifics/PairPotential_Morse.cpp

r3b1798	r0516fd
67	67	;
68	68	const std::string PairPotential_Morse::potential_token("morse");
69		Coordinator::ptr PairPotential_Morse::coordinator(~~new TwoBody_Length(~~));
	69	Coordinator::ptr PairPotential_Morse::coordinator(Memory::ignore(new TwoBody_Length()));
70	70
71	71	PairPotential_Morse::PairPotential_Morse() :

src/Potentials/Specifics/ThreeBodyPotential_Angle.cpp

r3b1798	r0516fd
65	65	;
66	66	const std::string ThreeBodyPotential_Angle::potential_token("harmonic_angle");
67		Coordinator::ptr ThreeBodyPotential_Angle::coordinator(~~new ThreeBody_Angle(~~));
	67	Coordinator::ptr ThreeBodyPotential_Angle::coordinator(Memory::ignore(new ThreeBody_Angle()));
68	68
69	69	ThreeBodyPotential_Angle::ThreeBodyPotential_Angle() :

src/Tesselation/BoundaryLineSet.cpp

r3b1798	r0516fd
309	309	ostream & operator <<(ostream &ost, const BoundaryLineSet &a)
310	310	{
311		ost << "[" << a.Nr << "\|" << a.endpoints[0]->node->getName() << " at " << a.endpoints[0]->node->getPosition() << "," << a.endpoints[1]->node->getName() << " at " << a.endpoints[1]->node->getPosition() << "]";
	311	ost << "[" << a.Nr << "\|" << a.endpoints[0]->node->getName() << "," << a.endpoints[1]->node->getName() << "]";
	312	//ost << "[" << a.Nr << "\|" << a.endpoints[0]->node->getName() << " at " << a.endpoints[0]->node->getPosition() << "," << a.endpoints[1]->node->getName() << " at " << a.endpoints[1]->node->getPosition() << "]";
312	313	return ost;
313	314	}

src/Tesselation/BoundaryPointSet.cpp

r3b1798	r0516fd
133	133	ostream & operator <<(ostream &ost, const BoundaryPointSet &a)
134	134	{
135		ost << "[" << a.Nr << "\|" << a.getName() << " at " << a.getPosition() << "]";
	135	ost << "[" << a.Nr << "\|" << a.getName() << "]"; // " at " << a.getPosition() << "]";
136	136	return ost;
137	137	}

src/Tesselation/BoundaryTriangleSet.cpp

-              r3b1798
+              r0516fd
   // set endpoints
   int Counter = 0;
   LOG(4, "DEBUG: New triangle " << Nr << " with end points: ");
+  LOG(4, "DEBUG: New triangle " << Nr << " with end points ... and lines:");
   for (PointMap::iterator runner = OrderMap.begin(); runner != OrderMap.end(); runner++) {
     endpoints[Counter] = runner->second;
     LOG(4, "DEBUG:    " << *endpoints[Counter]);
+    LOG(4, "DEBUG:    " << *endpoints[Counter] << "\t\t" << *lines[Counter]);
     Counter++;
+  }
 …
   // 1. get intersection with plane
   LOG(3, "DEBUG: Looking for closest point of triangle " << *this << " to " << x << ".");
   LOG(3, "DEBUG: endpoints are " << endpoints[0]->node->getPosition() << ","
+  LOG(4, "DEBUG: Looking for closest point of triangle " << *this << " to " << x << ".");
+  LOG(5, "DEBUG: endpoints are " << endpoints[0]->node->getPosition() << ","
       << endpoints[1]->node->getPosition() << ", and " << endpoints[2]->node->getPosition() << ".");
   try {
 …
+  }
   Vector InPlane(ClosestPoint); // points from plane intersection to straight-down point
   LOG(4, "DEBUG: Closest point on triangle plane is " << ClosestPoint << ".");
+  LOG(5, "DEBUG: Closest point on triangle plane is " << ClosestPoint << ".");
   // 2. Calculate in plane part of line (x, intersection)
 …
     CrossPoint[i] = l.getClosestPoint(InPlane);
     // NOTE: direction of line is normalized, hence s must not necessarily be in [0,1] for the baseline
     LOG(4, "DEBUG: Closest point on line from " << (endpoints[(i+1)%3]->node->getPosition())
+    LOG(5, "DEBUG: Closest point on line from " << (endpoints[(i+1)%3]->node->getPosition())
         << " to " << (endpoints[i%3]->node->getPosition()) << " is " << CrossPoint[i] << ".");
     CrossPoint[i] -= (endpoints[(i+1)%3]->node->getPosition());  // cross point was returned as absolute vector
     const double s = CrossPoint[i].ScalarProduct(Direction)/Direction.NormSquared();
     LOG(5, "DEBUG: Factor s is " << s << ".");
+    LOG(6, "DEBUG: Factor s is " << s << ".");
     if ((s >= -MYEPSILON) && ((s-1.) <= MYEPSILON)) {
       CrossPoint[i] += (endpoints[(i+1)%3]->node->getPosition());  // make cross point absolute again
       LOG(5, "DEBUG: Crosspoint is " << CrossPoint[i] << ", intersecting BoundaryLine between "
+      LOG(6, "DEBUG: Crosspoint is " << CrossPoint[i] << ", intersecting BoundaryLine between "
           << endpoints[i % 3]->node->getPosition() << " and "
           << endpoints[(i + 1) % 3]->node->getPosition() << ".");
 …
       else
         CrossPoint[i] = (endpoints[i%3]->node->getPosition());
       LOG(5, "DEBUG: Crosspoint is " << CrossPoint[i] << ", intersecting outside of BoundaryLine between "
+      LOG(6, "DEBUG: Crosspoint is " << CrossPoint[i] << ", intersecting outside of BoundaryLine between "
           << endpoints[i % 3]->node->getPosition() << " and "
           << endpoints[(i + 1) % 3]->node->getPosition() << ".");
 …
     ShortestDistance = InPlane.DistanceSquared(x);
+  }
   LOG(3, "DEBUG: Closest Point is " << ClosestPoint << " with shortest squared distance is " << ShortestDistance << ".");
+  LOG(4, "DEBUG: Closest Point is " << ClosestPoint << " with shortest squared distance is " << ShortestDistance << ".");
   return ShortestDistance;
 …
+{
   //Info FunctionInfo(__func__);
+  LOG(1, "INFO: Checking " << Points[0] << "," << Points[1] << "," << Points[2] << " against " << endpoints[0] << "," << endpoints[1] << "," << endpoints[2] << ".");
+  return (((endpoints[0] == Points[0]) || (endpoints[0] == Points[1]) || (endpoints[0] == Points[2])) && ((endpoints[1] == Points[0]) || (endpoints[1] == Points[1]) || (endpoints[1] == Points[2])) && ((endpoints[2] == Points[0]) || (endpoints[2] == Points[1]) || (endpoints[2] == Points[2])
+  LOG(5, "DEBUG: Checking " << *Points[0] << "," << *Points[1] << "," << *Points[2]
+      << " against " << *this); //*endpoints[0] << "," << *endpoints[1] << "," << *endpoints[2] << ".");
+  return (((endpoints[0] == Points[0]) || (endpoints[0] == Points[1]) || (endpoints[0] == Points[2]))
+      && ((endpoints[1] == Points[0]) || (endpoints[1] == Points[1]) || (endpoints[1] == Points[2]))
+      && ((endpoints[2] == Points[0]) || (endpoints[2] == Points[1]) || (endpoints[2] == Points[2])
   ));

src/Tesselation/CandidateForTesselation.cpp

-              r3b1798
+              r0516fd
   if (!pointlist.empty())
     LOG(3, "DEBUG: Checking whether sphere contains candidate list and baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ...");
+    LOG(3, "DEBUG: Checking validity whether sphere contains candidate list and baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ...");
   else
     LOG(3, "DEBUG: Checking whether sphere with no candidates contains baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ...");
+    LOG(3, "DEBUG: Checking validity whether sphere with no candidates contains baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ...");
   // check baseline for OptCenter and OtherOptCenter being on sphere's surface
   for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) {
 …
         return false;
       } else {
         LOG(3, "DEBUG: Candidate " << *Walker << " is inside by " << distance << ".");
+        LOG(4, "DEBUG: Candidate " << *Walker << " is inside by " << distance << ".");
+      }
+    }
+  }
   LOG(2, "DEBUG: Checking whether sphere contains no others points ...");
+  LOG(3, "DEBUG: Checking validity whether sphere contains no others points ...");
   bool flag = true;
   for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) {
 …
+    {
       LOG(3, "DEBUG: The following atoms are inside sphere at " << (*VRunner) << ":");
+      LOG(4, "DEBUG: The following atoms are inside sphere at " << (*VRunner) << ":");
       for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
         LOG(3, "DEBUG:  " << *(*Runner) << " with distance " << (*Runner)->distance(*(*VRunner)) << ".");
+        LOG(4, "DEBUG:  " << *(*Runner) << " with distance " << (*Runner)->distance(*(*VRunner)) << ".");
+    }
     // remove baseline's endpoints and candidates
     for (int i = 0; i < 2; i++) {
       LOG(3, "DEBUG: removing baseline tesselpoint " << *BaseLine->endpoints[i]->node << ".");
+      LOG(5, "DEBUG: removing baseline tesselpoint " << *BaseLine->endpoints[i]->node << ".");
       ListofPoints->remove(BaseLine->endpoints[i]->node);
+    }
     for (TesselPointList::const_iterator Runner = pointlist.begin(); Runner != pointlist.end(); ++Runner) {
       LOG(3, "DEBUG: removing candidate tesselpoint " << *(*Runner) << ".");
+      LOG(5, "DEBUG: removing candidate tesselpoint " << *(*Runner) << ".");
       ListofPoints->remove(*Runner);
+    }

src/Tesselation/boundary.cpp

-              r3b1798
+              r0516fd
  * \return volume difference between the non- and the created convex envelope
  */
+double ConvexizeNonconvexEnvelope(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename)
+double ConvexizeNonconvexEnvelope(
+    Tesselation *&TesselStruct,
+    const molecule * const mol,
+    const char * const filename,
+    bool DebugOutputEveryStep)
+{
         //Info FunctionInfo(__func__);
 …
+  }
+  // First step: RemovePointFromTesselatedSurface
+  LOG(1, "INFO: Making tesselated surface with " << TesselStruct->TrianglesOnBoundaryCount
+      << " convex ...");
+  // first purge all degenerate triangles
+  TesselStruct->RemoveDegeneratedTriangles();
   do {
     Concavity = false;
+    sprintf(dummy, "-first-%d", run);
+    //CalculateConcavityPerBoundaryPoint(TesselStruct);
+    StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);
+    if (DebugOutputEveryStep) {
+      sprintf(dummy, "-%d", run++);
+      //CalculateConcavityPerBoundaryPoint(TesselStruct);
+      LOG(1, "INFO: Writing " << run << "th tesselation file.");
+      StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);
+    }
+    // first step: remove all full-concave point
     PointRunner = TesselStruct->PointsOnBoundary.begin();
     PointAdvance = PointRunner; // we need an advanced point, as the PointRunner might get removed
 …
       PointAdvance++;
       point = PointRunner->second;
+      LOG(1, "INFO: Current point is " << *point << ".");
+      for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
+        line = LineRunner->second;
+        LOG(1, "INFO: Current line of point " << *point << " is " << *line << ".");
+        if (!line->CheckConvexityCriterion()) {
+          // remove the point if needed
+          LOG(1, "... point " << *point << " cannot be on convex envelope.");
+          volume += TesselStruct->RemovePointFromTesselatedSurface(point);
+          sprintf(dummy, "-first-%d", ++run);
+          StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);
+      LOG(2, "INFO: Current point is " << *point << ".");
+      // check that at least a single line is concave
+      LineMap::iterator LineRunner = point->lines.begin();
+      for (; LineRunner != point->lines.end(); LineRunner++) {
+        const BoundaryLineSet * line = LineRunner->second;
+        LOG(3, "INFO: Current line of point " << *point << " is " << *line << ".");
+        if (!line->CheckConvexityCriterion())
+          break;
+      }
+      // remove the point if needed
+      if (LineRunner != point->lines.end()) {
+        const double tmp = TesselStruct->RemoveFullConcavePointFromTesselatedSurface(point);
+        if (tmp > 0.) {
+          volume += tmp;
           Concavity = true;
-          break;
+        }
+      }
 …
+    }
+    sprintf(dummy, "-second-%d", run);
+    //CalculateConcavityPerBoundaryPoint(TesselStruct);
+    StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);
+    // second step: PickFarthestofTwoBaselines
+    if (DebugOutputEveryStep) {
+      sprintf(dummy, "-%d", run++);
+      //CalculateConcavityPerBoundaryPoint(TesselStruct);
+      LOG(1, "INFO: Writing " << run << "th tesselation file.");
+      StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);
+    }
+    // second step: flip baselines, i.e. add general tetraeder at concave lines
+    // when the tetraeder does not intersect with other already present triangles
     LineRunner = TesselStruct->LinesOnBoundary.begin();
     LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
+    std::map<double, std::pair<BoundaryLineSet *, double> > GainMap;
     while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
       LineAdvance++;
       line = LineRunner->second;
+      LOG(1, "INFO: Picking farthest baseline for line is " << *line << ".");
+      // take highest of both lines
+      if (TesselStruct->IsConvexRectangle(line) == NULL) {
+        const double tmp = TesselStruct->PickFarthestofTwoBaselines(line);
+        volume += tmp;
+        if (tmp != 0.) {
+          TesselStruct->FlipBaseline(line);
+          Concavity = true;
+      if (!line->CheckConvexityCriterion()) {
+        LOG(2, "INFO: concave line is " << *line << ".");
+        // gather the other points
+        BoundaryPointSet *BPS[4];
+        int m = 0;
+        {
+          for (TriangleMap::iterator runner = line->triangles.begin(); runner != line->triangles.end(); runner++)
+            for (int j = 0; j < 3; j++) // all of their endpoints and baselines
+              if (!line->ContainsBoundaryPoint(runner->second->endpoints[j])) // and neither of its endpoints
+                BPS[m++] = runner->second->endpoints[j];
+        }
+        BPS[2] = line->endpoints[0];
+        BPS[3] = line->endpoints[1];
+        LOG(3, "DEBUG: other line would consist of " << *BPS[0] << " and "
+            << *BPS[1] << ".");
+        // check for already present (third) side of the tetraeder as we then
+        // would create a degenerate triangle
+        bool TetraederSidePresent = false;
+        {
+          class TesselPoint *TriangleCandidates[3];
+          TriangleCandidates[0] = BPS[0]->node;
+          TriangleCandidates[1] = BPS[1]->node;
+          TriangleCandidates[2] = BPS[2]->node;
+          if ((TesselStruct->GetPresentTriangle(TriangleCandidates) != NULL)) {
+            LOG(2, "REJECT: Triangle side " << *TriangleCandidates[0] << ","
+                << *TriangleCandidates[1] << "," << *TriangleCandidates[2] << " present.");
+            TetraederSidePresent = true;
+          }
+          TriangleCandidates[2] = BPS[3]->node;
+          if ((TesselStruct->GetPresentTriangle(TriangleCandidates) != NULL)) {
+            LOG(2, "REJECT: Triangle side " << *TriangleCandidates[0] << ","
+                << *TriangleCandidates[1] << "," << *TriangleCandidates[2] << " present.");
+            TetraederSidePresent = true;
+          }
+        }
+        if ((BPS[0] != BPS[1]) && (m == 2) && (!TetraederSidePresent)) {
+          // check whether all adjacent triangles do not intersect with new line
+          bool no_line_intersects = true;
+          Vector Intersection;
+          TriangleSet triangles;
+          TriangleSet *firsttriangles = TesselStruct->GetAllTriangles(line->endpoints[0]);
+          TriangleSet *secondtriangles = TesselStruct->GetAllTriangles(line->endpoints[1]);
+          triangles.insert( firsttriangles->begin(), firsttriangles->end() );
+          triangles.insert( secondtriangles->begin(), secondtriangles->end() );
+          delete firsttriangles;
+          delete secondtriangles;
+          for (TriangleSet::const_iterator triangleiter = triangles.begin();
+              triangleiter != triangles.end(); ++triangleiter) {
+            const BoundaryTriangleSet * triangle = *triangleiter;
+            bool line_intersects = triangle->GetIntersectionInsideTriangle(
+                BPS[0]->node->getPosition(),
+                BPS[1]->node->getPosition(),
+                Intersection);
+            // switch result when coinciding with endpoint
+            bool concave_adjacent_line = false;
+            bool intersection_is_endnode = false;
+            for (int j=0;j<2;++j) {
+              if (Intersection.DistanceSquared(BPS[j]->node->getPosition()) < MYEPSILON) {
+                intersection_is_endnode = true;
+                // check whether its an adjacent triangle and if it's concavely connected
+                // only then are we in danger of cutting through it and need to check
+                // sign of normal vector and intersecting line
+                for (int i=0;i<2;++i)
+                  for (int lineindex=0;lineindex < 3;++lineindex)
+                    if ((triangle->lines[lineindex]->ContainsBoundaryPoint(line->endpoints[i]))
+                        && (triangle->lines[lineindex]->ContainsBoundaryPoint(BPS[j]))) {
+                      concave_adjacent_line = !triangle->lines[lineindex]->CheckConvexityCriterion();
+                  }
+                if (concave_adjacent_line) {
+                  const Vector intersector =
+                      BPS[(j+1)%2]->node->getPosition() - Intersection;
+                  if (triangle->NormalVector.ScalarProduct(intersector) >= -MYEPSILON) {
+                    LOG(4, "ACCEPT: Intersection coincides with first endpoint "
+                        << *BPS[j] << ".");
+                    line_intersects = false;
+                  } else {
+                    LOG(4, "REJECT: Intersection ends on wrong side of triangle.");
+                  }
+                } else {
+                  LOG(4, "ACCEPT: Intersection coincides with first endpoint "
+                      << *BPS[j] << ".");
+                  line_intersects = false;
+                }
+              }
+            }
+            // if we have an intersection, check that it is within either
+            // endpoint, i.e. check that scalar product between vectors going
+            // from intersction to either endpoint has negative sign (both
+            // vectors point in opposite directions)
+            if (!intersection_is_endnode && line_intersects) {
+              const Vector firstvector = BPS[0]->node->getPosition() - Intersection;
+              const Vector secondvector = BPS[1]->node->getPosition() - Intersection;
+              if (firstvector.ScalarProduct(secondvector) >= 0)
+                line_intersects = false;
+            }
+            no_line_intersects &= !line_intersects;
+          }
+          if (no_line_intersects) {
+            // calculate the volume
+            const double tmp = line->CalculateConvexity();
+            const double gain =
+              CalculateVolumeofGeneralTetraeder(
+                  BPS[0]->node->getPosition(),
+                  BPS[1]->node->getPosition(),
+                  BPS[2]->node->getPosition(),
+                  BPS[3]->node->getPosition());
+            GainMap.insert(std::make_pair(tmp, std::make_pair(line,gain) ));
+            LOG(2, "DEBUG: Adding concave line " << *line << " with gain of "
+                << gain << ".");
+          } else {
+            // if 2 or 3 don't
+            LOG(2, "DEBUG: We don't added concave line " << *line
+                << " as other line intersects with adjacent triangles.");
+          }
+        }
+      }
       LineRunner = LineAdvance;
+    }
+    run++;
+  } while (Concavity);
+  //CalculateConcavityPerBoundaryPoint(TesselStruct);
+  //StoreTrianglesinFile(mol, filename, "-third");
+  // 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;
+//    LOG(1, "INFO: Current line is " << *line << ".");
+//    //if (LineAdvance != TesselStruct->LinesOnBoundary.end())
+//      //LOG(1, "INFO: Next line will be " << *(LineAdvance->second) << ".");
+//    if (!line->CheckConvexityCriterion(out)) {
+//      LOG(1, "INFO: ... line " << *line << " is concave, flipping it.");
+//
+//      // take highest of both lines
+//      point = TesselStruct->IsConvexRectangle(line);
+//      if (point != NULL)
+//        volume += TesselStruct->RemovePointFromTesselatedSurface(point);
+//    }
+//    LineRunner = LineAdvance;
+//  }
+    // flip line with most gain
+    if (!GainMap.empty()) {
+      line = GainMap.begin()->second.first;
+      const double tmp = GainMap.begin()->second.second;
+      volume += tmp;
+//      GainMap.clear();
+      // and flip the line
+      LOG(1, "INFO: Flipping current most concave line " << *line << " with gain of "
+          << tmp << ".");
+      TesselStruct->FlipBaseline(line);
+      Concavity = true;
+    }
+  } while ((Concavity)); // && (run < 100)
   CalculateConcavityPerBoundaryPoint(TesselStruct);
 …
   // end
   LOG(0, "Volume is " << volume << ".");
+  LOG(0, "RESULT: Added volume in convexization is " << volume << ".");
   return volume;
 };
 …
 //  status = CheckListOfBaselines(TesselStruct);
   cout << "before correction" << endl;
+//  cout << "before correction" << endl;
   // store before correction
 …
   // write final envelope
   CalculateConcavityPerBoundaryPoint(TesselStruct);
   cout << "after correction" << endl;
+//  cout << "after correction" << endl;
   StoreTrianglesinFile(mol, TesselStruct, filename, "");

src/Tesselation/boundary.hpp

r3b1798	r0516fd
39	39	/******************************************** declarations *****************************/
40	40
41		double ConvexizeNonconvexEnvelope(class Tesselation &TesselStruct, const molecule const mol, const char * const filename);
	41	double ConvexizeNonconvexEnvelope(class Tesselation &TesselStruct, const molecule const mol, const char * const filename, bool DebugOutputEveryStep = false);
42	42	void FindConvexBorder(const molecule* const mol, Boundaries BoundaryPts, Tesselation &TesselStruct, const LinkedCell_deprecated LCList, const char filename);
43	43	Vector* FindEmbeddingHole(MoleculeListClass mols, molecule srcmol);

src/Tesselation/tesselation.cpp

-              r3b1798
+              r0516fd
 #include "CodePatterns/MemDebug.hpp"
+#include <algorithm>
+#include <boost/foreach.hpp>
 #include <fstream>
 #include <iomanip>
+#include <iterator>
 #include <sstream>
 …
 class molecule;
 const char *TecplotSuffix=".dat";
 const char *Raster3DSuffix=".r3d";
 const char *VRMLSUffix=".wrl";
 const double ParallelEpsilon=1e-3;
+const char *TecplotSuffix = ".dat";
+const char *Raster3DSuffix = ".r3d";
+const char *VRMLSUffix = ".wrl";
+const double ParallelEpsilon = 1e-3;
 const double Tesselation::HULLEPSILON = 1e-9;
 …
  */
 Tesselation::Tesselation() :
+  PointsOnBoundaryCount(0),
+  LinesOnBoundaryCount(0),
+  TrianglesOnBoundaryCount(0),
+  LastTriangle(NULL),
+  TriangleFilesWritten(0),
+  InternalPointer(PointsOnBoundary.begin())
+    PointsOnBoundaryCount(0), LinesOnBoundaryCount(0), TrianglesOnBoundaryCount(0), LastTriangle(NULL), TriangleFilesWritten(0), InternalPointer(PointsOnBoundary.begin())
+{
   //Info FunctionInfo(__func__);
 …
+{
   // create linkedcell
         LinkedCell_deprecated *LinkedList = new LinkedCell_deprecated(cloud, 2.*SPHERERADIUS);
         // check for at least three points
+        {
           bool ThreePointsFound = true;
           cloud.GoToFirst();
     for (size_t i=0;i<3;++i, cloud.GoToNext())
+  LinkedCell_deprecated *LinkedList = new LinkedCell_deprecated(cloud, 2. * SPHERERADIUS);
+  // check for at least three points
+  {
+    bool ThreePointsFound = true;
+    cloud.GoToFirst();
+    for (size_t i = 0; i < 3; ++i, cloud.GoToNext())
       ThreePointsFound &= (!cloud.IsEnd());
     cloud.GoToFirst();
 …
       return;
+    }
+        }
         // find a starting triangle
+  }
+  // find a starting triangle
   FindStartingTriangle(SPHERERADIUS, LinkedList);
 …
         //the line is there, so there is a triangle, but only one.
         const bool TesselationFailFlag = FindNextSuitableTriangle(*baseline, *T, SPHERERADIUS, LinkedList);
+        ASSERT( TesselationFailFlag,
+            "Tesselation::operator() - no suitable candidate triangle found.");
+        ASSERT(TesselationFailFlag, "Tesselation::operator() - no suitable candidate triangle found.");
+      }
+    }
     // 2b. search for smallest ShortestAngle among all candidates
     double ShortestAngle = 4.*M_PI;
+    double ShortestAngle = 4. * M_PI;
     for (CandidateMap::iterator Runner = OpenLines.begin(); Runner != OpenLines.end(); Runner++) {
       if (Runner->second->ShortestAngle < ShortestAngle) {
 …
+      }
+    }
     if ((ShortestAngle == 4.*M_PI) || (baseline->pointlist.empty()))
+    if ((ShortestAngle == 4. * M_PI) || (baseline->pointlist.empty()))
       OneLoopWithoutSuccessFlag = false;
     else {
 …
 double Tesselation::getVolumeOfConvexEnvelope(const bool IsAngstroem) const
+{
+  // calculate center of gravity
+  Vector center;
+  if (!PointsOnBoundary.empty()) {
+    for (PointMap::const_iterator iter = PointsOnBoundary.begin();
+        iter != PointsOnBoundary.end(); ++iter)
+      center += iter->second->node->getPosition();
+    center *= 1./(double)PointsOnBoundary.size();
+  }
+  // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
   double volume = 0.;
+  Vector x;
+  Vector y;
+  // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
+  for (TriangleMap::const_iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++)
+    { // go through every triangle, calculate volume of its pyramid with CoG as peak
+      x = runner->second->getEndpoint(0) - runner->second->getEndpoint(1);
+      const double G = runner->second->getArea();
+      x = runner->second->getPlane().getNormal();
+      x.Scale(runner->second->getEndpoint(1).ScalarProduct(x));
+      const double h = x.Norm(); // distance of CoG to triangle
+      const double PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
+      LOG(1, "INFO: Area of triangle is " << setprecision(10) << G << " "
+          << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is "
+          << h << " and the volume is " << PyramidVolume << " "
+          << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
+      volume += PyramidVolume;
+    }
+  LOG(0, "RESULT: The summed volume is " << setprecision(6)
+      << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
+  for (TriangleMap::const_iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) { // go through every triangle, calculate volume of its pyramid with CoG as peak
+    const double TetrahedronVolume = CalculateVolumeofGeneralTetraeder(
+        runner->second->endpoints[0]->getPosition(),
+        runner->second->endpoints[1]->getPosition(),
+        runner->second->endpoints[2]->getPosition(),
+        center);
+    LOG(1, "INFO: volume of tetrahedron is " << setprecision(10) << TetrahedronVolume
+        << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
+    volume += TetrahedronVolume;
+  }
+  LOG(0, "RESULT: The summed volume is " << setprecision(6) << volume
+      << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
   return volume;
 …
   // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
+  for (TriangleMap::const_iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++)
+    { // go through every triangle, calculate volume of its pyramid with CoG as peak
+                const double area = runner->second->getArea();
+                LOG(1, "INFO: Area of triangle is " << setprecision(10) << area << " "
+          << (IsAngstroem ? "angstrom" : "atomiclength") << "^2.");
+      surfacearea += area;
+    }
+  LOG(0, "RESULT: The summed surface area is " << setprecision(6)
+      << surfacearea << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
+  for (TriangleMap::const_iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) { // go through every triangle, calculate volume of its pyramid with CoG as peak
+    const double area = runner->second->getArea();
+    LOG(1, "INFO: Area of triangle is " << setprecision(10) << area << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^2.");
+    surfacearea += area;
+  }
+  LOG(0, "RESULT: The summed surface area is " << setprecision(6) << surfacearea << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
   return surfacearea;
+}
 /** Gueses first starting triangle of the convex envelope.
 …
         tmp = B->second->node->DistanceSquared(C->second->node->getPosition());
         distance += tmp * tmp;
         DistanceMMap.insert(DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator> (B, C)));
+        DistanceMMap.insert(DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator>(B, C)));
+      }
+    }
 …
     // 2. next, we have to check whether all points reside on only one side of the triangle
     // 3. construct plane vector
+    PlaneVector = Plane(A->second->node->getPosition(),
+                        baseline->second.first->second->node->getPosition(),
+                        baseline->second.second->second->node->getPosition()).getNormal();
+    PlaneVector = Plane(A->second->node->getPosition(), baseline->second.first->second->node->getPosition(), baseline->second.second->second->node->getPosition()).getNormal();
     LOG(2, "Plane vector of candidate triangle is " << PlaneVector);
     // 4. loop over all points
 …
         // prepare some auxiliary vectors
         Vector BaseLineCenter, BaseLine;
+        BaseLineCenter = 0.5 * ((baseline->second->endpoints[0]->node->getPosition()) +
+                                (baseline->second->endpoints[1]->node->getPosition()));
+        BaseLineCenter = 0.5 * ((baseline->second->endpoints[0]->node->getPosition()) + (baseline->second->endpoints[1]->node->getPosition()));
         BaseLine = (baseline->second->endpoints[0]->node->getPosition()) - (baseline->second->endpoints[1]->node->getPosition());
 …
         // vector in propagation direction (out of triangle)
         // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
         PropagationVector = Plane(BaseLine, NormalVector,0).getNormal();
+        PropagationVector = Plane(BaseLine, NormalVector, 0).getNormal();
         TempVector = CenterVector - (baseline->second->endpoints[0]->node->getPosition()); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
         //LOG(0, "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << ".");
 …
             // in case NOT both were found, create virtually this triangle, get its normal vector, calculate angle
             flag = true;
+            VirtualNormalVector = Plane((baseline->second->endpoints[0]->node->getPosition()),
+                                        (baseline->second->endpoints[1]->node->getPosition()),
+                                        (target->second->node->getPosition())).getNormal();
+            TempVector = (1./3.) * ((baseline->second->endpoints[0]->node->getPosition()) +
+                                    (baseline->second->endpoints[1]->node->getPosition()) +
+                                    (target->second->node->getPosition()));
+            VirtualNormalVector = Plane((baseline->second->endpoints[0]->node->getPosition()), (baseline->second->endpoints[1]->node->getPosition()), (target->second->node->getPosition())).getNormal();
+            TempVector = (1. / 3.) * ((baseline->second->endpoints[0]->node->getPosition()) + (baseline->second->endpoints[1]->node->getPosition()) + (target->second->node->getPosition()));
             TempVector -= (*Center);
             // make it always point outward
 …
               winner = target;
               LOG(5, "DEBUG: New winner " << *winner->second->node << " due to smaller angle between normal vectors.");
+            } else if (fabs(SmallestAngle - TempAngle) < MYEPSILON) { // check the angle to propagation, both possible targets are in one plane! (their normals have same angle)
+              // hence, check the angles to some normal direction from our base line but in this common plane of both targets...
+              helper = (target->second->node->getPosition()) - BaseLineCenter;
+              helper.ProjectOntoPlane(BaseLine);
+              // ...the one with the smaller angle is the better candidate
+              TempVector = (target->second->node->getPosition()) - BaseLineCenter;
+              TempVector.ProjectOntoPlane(VirtualNormalVector);
+              TempAngle = TempVector.Angle(helper);
+              TempVector = (winner->second->node->getPosition()) - BaseLineCenter;
+              TempVector.ProjectOntoPlane(VirtualNormalVector);
+              if (TempAngle < TempVector.Angle(helper)) {
+                TempAngle = NormalVector.Angle(VirtualNormalVector);
+                SmallestAngle = TempAngle;
+                winner = target;
+                LOG(5, "DEBUG: New winner " << *winner->second->node << " due to smaller angle " << TempAngle << " to propagation direction.");
+            } else
+              if (fabs(SmallestAngle - TempAngle) < MYEPSILON) { // check the angle to propagation, both possible targets are in one plane! (their normals have same angle)
+                // hence, check the angles to some normal direction from our base line but in this common plane of both targets...
+                helper = (target->second->node->getPosition()) - BaseLineCenter;
+                helper.ProjectOntoPlane(BaseLine);
+                // ...the one with the smaller angle is the better candidate
+                TempVector = (target->second->node->getPosition()) - BaseLineCenter;
+                TempVector.ProjectOntoPlane(VirtualNormalVector);
+                TempAngle = TempVector.Angle(helper);
+                TempVector = (winner->second->node->getPosition()) - BaseLineCenter;
+                TempVector.ProjectOntoPlane(VirtualNormalVector);
+                if (TempAngle < TempVector.Angle(helper)) {
+                  TempAngle = NormalVector.Angle(VirtualNormalVector);
+                  SmallestAngle = TempAngle;
+                  winner = target;
+                  LOG(5, "DEBUG: New winner " << *winner->second->node << " due to smaller angle " << TempAngle << " to propagation direction.");
+                } else
+                  LOG(5, "DEBUG: Keeping old winner " << *winner->second->node << " due to smaller angle to propagation direction.");
               } else
+                LOG(5, "DEBUG: Keeping old winner " << *winner->second->node << " due to smaller angle to propagation direction.");
+            } else
+              LOG(5, "DEBUG: Keeping old winner " << *winner->second->node << " due to smaller angle between normal vectors.");
+                LOG(5, "DEBUG: Keeping old winner " << *winner->second->node << " due to smaller angle between normal vectors.");
+          }
         } // end of loop over all boundary points
 …
   cloud.GoToFirst();
   PointCloudAdaptor< Tesselation, MapValueIterator<Tesselation::iterator> > newcloud(this, cloud.GetName());
+  PointCloudAdaptor<Tesselation, MapValueIterator<Tesselation::iterator> > newcloud(this, cloud.GetName());
   BoundaryPoints = new LinkedCell_deprecated(newcloud, 5.);
   while (!cloud.IsEnd()) { // we only have to go once through all points, as boundary can become only bigger
 …
  * @param *b second endpoint
  * @param n index of Tesselation::BLS giving the line with both endpoints
+ */
+void Tesselation::AddTesselationLine(const Vector * const OptCenter, const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n)
+ * @return true - inserted a new line, false - present line used
+ */
+bool Tesselation::AddTesselationLine(const Vector * const OptCenter, const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n)
+{
   bool insertNewLine = true;
 …
       if (FindLine->second->triangles.size() == 1) {
         CandidateMap::iterator Finder = OpenLines.find(FindLine->second);
+        if (!Finder->second->pointlist.empty())
+          LOG(4, "DEBUG: line " << *(FindLine->second) << " is open with candidate " << **(Finder->second->pointlist.begin()) << ".");
+        else
+          LOG(4, "DEBUG: line " << *(FindLine->second) << " is open with no candidate.");
+        // get open line
+        for (TesselPointList::const_iterator CandidateChecker = Finder->second->pointlist.begin(); CandidateChecker != Finder->second->pointlist.end(); ++CandidateChecker) {
+          if ((*(CandidateChecker) == candidate->node) && (OptCenter == NULL || OptCenter->DistanceSquared(Finder->second->OptCenter) < MYEPSILON )) { // stop searching if candidate matches
+            LOG(4, "ACCEPT: Candidate " << *(*CandidateChecker) << " has the right center " << Finder->second->OptCenter << ".");
+        ASSERT(Finder != OpenLines.end(), "Tesselation::AddTesselationLine() - " + toString(*FindLine->second) + " is not a new line and not in OpenLines.");
+        if (Finder->second != NULL) {
+          if (!Finder->second->pointlist.empty())
+            LOG(4, "DEBUG: line " << *(FindLine->second) << " is open with candidate " << **(Finder->second->pointlist.begin()) << ".");
+          else {
+            LOG(4, "ACCEPT: line " << *(FindLine->second) << " is open with no candidate.");
             insertNewLine = false;
             WinningLine = FindLine->second;
-            break;
-          } else {
-            LOG(5, "REJECT: Candidate " << *(*CandidateChecker) << "'s center " << Finder->second->OptCenter << " does not match desired on " << *OptCenter << ".");
+          }
+          // get open line
+          for (TesselPointList::const_iterator CandidateChecker = Finder->second->pointlist.begin(); CandidateChecker != Finder->second->pointlist.end(); ++CandidateChecker) {
+            if ((*(CandidateChecker) == candidate->node) && (OptCenter == NULL || OptCenter->DistanceSquared(Finder->second->OptCenter) < MYEPSILON)) { // stop searching if candidate matches
+              LOG(4, "ACCEPT: Candidate " << *(*CandidateChecker) << " has the right center " << Finder->second->OptCenter << ".");
+              insertNewLine = false;
+              WinningLine = FindLine->second;
+              break;
+            } else {
+              LOG(5, "REJECT: Candidate " << *(*CandidateChecker) << "'s center " << Finder->second->OptCenter << " does not match desired on " << *OptCenter << ".");
+            }
+          }
+        } else {
+          LOG(4, "ACCEPT: There are no candidates for present open line, use what we have.");
+          insertNewLine = false;
+          WinningLine = FindLine->second;
+        }
+      }
 …
     AddExistingTesselationTriangleLine(WinningLine, n);
+  }
+  return insertNewLine;
+}
+;
 …
   // also add to open lines
   CandidateForTesselation *CFT = new CandidateForTesselation(BLS[n]);
   OpenLines.insert(pair<BoundaryLineSet *, CandidateForTesselation *> (BLS[n], CFT));
+  OpenLines.insert(pair<BoundaryLineSet *, CandidateForTesselation *>(BLS[n], CFT));
+}
+;
 …
       } else {
         output << "still attached to another triangle: ";
-        OpenLines.insert(pair<BoundaryLineSet *, CandidateForTesselation *> (triangle->lines[i], NULL));
         for (TriangleMap::iterator TriangleRunner = triangle->lines[i]->triangles.begin(); TriangleRunner != triangle->lines[i]->triangles.end(); TriangleRunner++)
           output << "\t[" << (TriangleRunner->second)->Nr << "|" << *((TriangleRunner->second)->endpoints[0]) << ", " << *((TriangleRunner->second)->endpoints[1]) << ", " << *((TriangleRunner->second)->endpoints[2]) << "] \t";
+        OpenLines.insert(pair<BoundaryLineSet *, CandidateForTesselation *>(triangle->lines[i], NULL));
+      }
       LOG(3, "DEBUG: " << output.str());
 …
   else
     Numbers[1] = -1;
+  // erase from OpenLines if present
+  {
+    CandidateMap::iterator openlineiter = OpenLines.find(line);
+    if (openlineiter != OpenLines.end())
+      OpenLines.erase(openlineiter);
+  }
   for (int i = 0; i < 2; i++) {
 …
         LOG(4, "DEBUG: " << *line->endpoints[i] << " has no more lines it's attached to, erasing.");
         RemoveTesselationPoint(line->endpoints[i]);
+      } else if (DoLog(0)) {
+        std::stringstream output;
+        output << "DEBUG: " << *line->endpoints[i] << " has still lines it's attached to: ";
+        for (LineMap::iterator LineRunner = line->endpoints[i]->lines.begin(); LineRunner != line->endpoints[i]->lines.end(); LineRunner++)
+          output << "[" << *(LineRunner->second) << "] \t";
+        LOG(4, output.str());
+      }
+      } else
+        if (DoLog(0)) {
+          std::stringstream output;
+          output << "DEBUG: " << *line->endpoints[i] << " has still lines it's attached to: ";
+          for (LineMap::iterator LineRunner = line->endpoints[i]->lines.begin(); LineRunner != line->endpoints[i]->lines.end(); LineRunner++)
+            output << "[" << *(LineRunner->second) << "] \t";
+          LOG(4, output.str());
+        }
       line->endpoints[i] = NULL; // free'd or not: disconnect
     } else
 …
   //Info FunctionInfo(__func__);
   LOG(3, "DEBUG: Checking whether sphere contains no others points ...");
+  LOG(3, "DEBUG: Checking degeneracy by whether sphere contains no others points ...");
   bool flag = true;
 …
   TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, &CandidateLine.OtherOptCenter);
+  LOG(3, "DEBUG: The following atoms are inside sphere at " << CandidateLine.OtherOptCenter << ":");
+  LOG(3, "DEBUG: CheckDegeneracy: There are " << ListofPoints->size() << " points inside the sphere.");
+  LOG(4, "DEBUG: The following atoms are inside sphere at " << CandidateLine.OtherOptCenter << ":");
   for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
     LOG(3, "DEBUG:   " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << ".");
+    LOG(4, "DEBUG:   " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << ".");
   // remove triangles's endpoints
 …
     LOG(3, "DEBUG: CheckDegeneracy: There are still " << ListofPoints->size() << " points inside the sphere.");
     flag = false;
     LOG(3, "DEBUG: External atoms inside of sphere at " << CandidateLine.OtherOptCenter << ":");
+    LOG(4, "DEBUG: External atoms inside of sphere at " << CandidateLine.OtherOptCenter << ":");
     for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
       LOG(3, "DEBUG:   " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << ".");
+      LOG(4, "DEBUG:   " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << ".");
+  }
   delete ListofPoints;
 …
   // 1. searching topmost point with respect to each axis
+  LOG(2, "INFO: Searching for topmost pointer from each axis.");
   for (int i = 0; i < NDIM; i++) { // each axis
     LC->n[i] = LC->N[i] - 1; // current axis is topmost cell
 …
+  }
   if (DoLog(1)) {
+  if (DoLog(3)) {
     std::stringstream output;
     output << "Found maximum coordinates: ";
 …
     BaseLine = new BoundaryLineSet();
     BaseLine->endpoints[0] = new BoundaryPointSet(MaxPoint[k]);
     LOG(2, "DEBUG: Coordinates of start node at " << *BaseLine->endpoints[0]->node << ".");
+    LOG(2, "INFO: Coordinates of start node at " << *BaseLine->endpoints[0]->node << ".");
     double ShortestAngle;
 …
+    }
     BaseLine->endpoints[1] = new BoundaryPointSet(Temporary);
     LOG(1, "INFO: Second node is at " << *Temporary << ".");
+    LOG(2, "INFO: Second node is at " << *Temporary << ".");
     // construct center of circle
     CircleCenter = 0.5 * ((BaseLine->endpoints[0]->node->getPosition()) + (BaseLine->endpoints[1]->node->getPosition()));
     LOG(1, "INFO: CircleCenter is at " << CircleCenter << ".");
+    LOG(4, "DEBUG: CircleCenter is at " << CircleCenter << ".");
     // construct normal vector of circle
     CirclePlaneNormal = (BaseLine->endpoints[0]->node->getPosition()) - (BaseLine->endpoints[1]->node->getPosition());
     LOG(1, "INFO: CirclePlaneNormal is at " << CirclePlaneNormal << ".");
+    LOG(4, "DEBUG: CirclePlaneNormal is at " << CirclePlaneNormal << ".");
     double radius = CirclePlaneNormal.NormSquared();
 …
     NormalVector.ProjectOntoPlane(CirclePlaneNormal);
     NormalVector.Normalize();
     LOG(1, "INFO: NormalVector is at " << NormalVector << ".");
+    LOG(4, "DEBUG: NormalVector is at " << NormalVector << ".");
     ShortestAngle = 2. * M_PI; // This will indicate the quadrant.
     SphereCenter = (CircleRadius * NormalVector) +  CircleCenter;
+    SphereCenter = (CircleRadius * NormalVector) + CircleCenter;
     // Now, NormalVector and SphereCenter are two orthonormalized vectors in the plane defined by CirclePlaneNormal (not normalized)
     // look in one direction of baseline for initial candidate
     try {
+      SearchDirection = Plane(CirclePlaneNormal, NormalVector,0).getNormal();  // whether we look "left" first or "right" first is not important ...
+    } catch(LinearAlgebraException) {
+      ELOG(1, "Vectors are linear dependent: "
+          << CirclePlaneNormal << ", " << NormalVector << ".");
+      SearchDirection = Plane(CirclePlaneNormal, NormalVector, 0).getNormal(); // whether we look "left" first or "right" first is not important ...
+    } catch (LinearAlgebraException &e) {
+      ELOG(1, "Vectors are linear dependent: " << CirclePlaneNormal << ", " << NormalVector << ".");
       delete BaseLine;
       continue;
 …
     // adding point 1 and point 2 and add the line between them
     LOG(2, "DEBUG: Found second point is at " << *BaseLine->endpoints[1]->node << ".");
+    LOG(3, "DEBUG: Found second point is at " << *BaseLine->endpoints[1]->node << ".");
     //LOG(1, "INFO: OldSphereCenter is at " << helper << ".");
 …
       for (TesselPointList::iterator it = OptCandidates.pointlist.begin(); it != OptCandidates.pointlist.end(); it++)
         output << *(*it);
       LOG(2, "DEBUG: List of third Points is: " << output.str());
+      LOG(3, "DEBUG: List of third Points is: " << output.str());
+    }
     if (!OptCandidates.pointlist.empty()) {
 …
 //  return result;
 //};
 /** This function finds a triangle to a line, adjacent to an existing one.
  * @param out output stream for debugging
 …
   // construct center of circle
+  CircleCenter = 0.5 * ((CandidateLine.BaseLine->endpoints[0]->node->getPosition()) +
+                        (CandidateLine.BaseLine->endpoints[1]->node->getPosition()));
+  CircleCenter = 0.5 * ((CandidateLine.BaseLine->endpoints[0]->node->getPosition()) + (CandidateLine.BaseLine->endpoints[1]->node->getPosition()));
   // construct normal vector of circle
+  CirclePlaneNormal = (CandidateLine.BaseLine->endpoints[0]->node->getPosition()) -
+                      (CandidateLine.BaseLine->endpoints[1]->node->getPosition());
+  CirclePlaneNormal = (CandidateLine.BaseLine->endpoints[0]->node->getPosition()) - (CandidateLine.BaseLine->endpoints[1]->node->getPosition());
   // calculate squared radius of circle
 …
     CircleRadius = RADIUS * RADIUS - radius / 4.;
     CirclePlaneNormal.Normalize();
     LOG(3, "DEBUG: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << ".");
     LOG(3, "DEBUG: OldSphereCenter is at " << T.SphereCenter << ".");
+    LOG(4, "DEBUG: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << ".");
+    LOG(4, "DEBUG: OldSphereCenter is at " << T.SphereCenter << ".");
     // construct SearchDirection and an "outward pointer"
     SearchDirection = Plane(RelativeSphereCenter, CirclePlaneNormal,0).getNormal();
+    SearchDirection = Plane(RelativeSphereCenter, CirclePlaneNormal, 0).getNormal();
     helper = CircleCenter - (ThirdPoint->node->getPosition());
     if (helper.ScalarProduct(SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
+    if (helper.ScalarProduct(SearchDirection) < -HULLEPSILON) // ohoh, SearchDirection points inwards!
       SearchDirection.Scale(-1.);
     LOG(3, "DEBUG: SearchDirection is " << SearchDirection << ".");
+    LOG(4, "DEBUG: SearchDirection is " << SearchDirection << ".");
     if (fabs(RelativeSphereCenter.ScalarProduct(SearchDirection)) > HULLEPSILON) {
       // rotated the wrong way!
 …
   } else {
     LOG(3, "DEBUG: Circumcircle for base line " << *CandidateLine.BaseLine << " and base triangle " << T << " is too big!");
+    LOG(4, "DEBUG: Circumcircle for base line " << *CandidateLine.BaseLine << " and base triangle " << T << " is too big!");
+  }
 …
     baseline = Runner->second;
     if (baseline->pointlist.empty()) {
       ASSERT((baseline->BaseLine->triangles.size() == 1),"Open line without exactly one attached triangle");
+      ASSERT((baseline->BaseLine->triangles.size() == 1), "Open line without exactly one attached triangle");
       T = (((baseline->BaseLine->triangles.begin()))->second);
       LOG(4, "DEBUG: Finding best candidate for open line " << *baseline->BaseLine << " of triangle " << *T);
 …
   TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, CandidateLine.BaseLine->endpoints[1]->node->getPosition());
+  {
+  if (DoLog(3)) {
     std::stringstream output;
     for (TesselPointList::iterator TesselRunner = connectedClosestPoints->begin(); TesselRunner != connectedClosestPoints->end(); ++TesselRunner)
       output << **TesselRunner;
     LOG(3, "DEBUG: List of Candidates for Turning Point " << *TurningPoint << ":");
+    LOG(3, "DEBUG: List of Candidates for Turning Point " << *TurningPoint << ":" << output.str());
+  }
 …
+  }
   delete (connectedClosestPoints);
+};
+}
+;
 /** for polygons (multiple candidates for a baseline) sets internal edges to the correct next candidate.
 …
       else {
         LOG(4, "DEBUG: line " << *(FindLine->second) << " is open with no candidate, setting to next Sprinter" << (*Sprinter));
         Finder->second->T = BTS;  // is last triangle
+        Finder->second->T = BTS; // is last triangle
         Finder->second->pointlist.push_back(Sprinter);
         Finder->second->ShortestAngle = 0.;
 …
+    }
+  }
+};
+}
+;
 /** If a given \a *triangle is degenerated, this adds both sides.
 …
   // give some verbose output about the whole procedure
   if (CandidateLine.T != NULL)
     LOG(2, "DEBUG: --> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << ".");
+    LOG(2, "INFO: --> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << ".");
   else
     LOG(2, "DEBUG: --> New starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle.");
+    LOG(2, "INFO: --> New starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle.");
   triangle = BTS;
 …
+;
+bool Tesselation::isConvex() const
+{
+  bool status = true;
+  // go through all lines on boundary
+  for (LineMap::const_iterator lineiter = LinesOnBoundary.begin();
+      lineiter != LinesOnBoundary.end(); ++lineiter) {
+    if (!lineiter->second->CheckConvexityCriterion()) {
+      LOG(2, "DEBUG: Line " << *lineiter->second << " is not convex.");
+      status = false;
+    }
+  }
+  return status;
+}
 /** Checks whether the quadragon of the two triangles connect to \a *Base is convex.
  * We look whether the closest point on \a *Base with respect to the other baseline is outside
 …
+    }
     for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
     LOG(4, "DEBUG: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << ".");
+      LOG(4, "DEBUG: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << ".");
       BaseLineNormal += (runner->second->NormalVector);
+    }
 …
   class BoundaryLineSet *OldLines[4], *NewLine;
   class BoundaryPointSet *OldPoints[2];
+  Vector BaseLineNormal;
+  Vector BaseLineNormal[2];
+  Vector OtherEndpoint[2]; // fourth point to either triangle
   int OldTriangleNrs[2], OldBaseLineNr;
   int i, m;
   // calculate NormalVector for later use
-  BaseLineNormal.Zero();
   if (Base->triangles.size() < 2) {
     ELOG(1, "Less than two triangles are attached to this baseline!");
     return NULL;
+  }
+  for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
+    LOG(1, "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << ".");
+    BaseLineNormal += (runner->second->NormalVector);
+  }
+  BaseLineNormal.Scale(-1. / 2.); // has to point inside for BoundaryTriangleSet::GetNormalVector()
+  {
+    int i=0;
+    for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
+      LOG(5, "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << ".");
+      OtherEndpoint[(i+1)%2] = runner->second->GetThirdEndpoint(Base)->node->getPosition();
+      BaseLineNormal[i++] = (runner->second->NormalVector);
+      ASSERT( i <= 2,
+          "Tesselation::FlipBaseline() - not exactly two triangles found.");
+    }
+  }
   // get the two triangles
 …
   // index OldLines and OldPoints
+  // note that oldlines[0,1] belong to first triangle and hence first normal
+  // vector and oldlines[2,3] belong to second triangle and its normal vector
   for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
     for (int j = 0; j < 3; j++) // all of their endpoints and baselines
 …
         OldPoints[m++] = runner->second->endpoints[j];
+  Vector BasePoints[2]; // endpoints of Base
+  if (OldLines[0]->ContainsBoundaryPoint(Base->endpoints[0])) {
+    BasePoints[0] = Base->endpoints[0]->node->getPosition();
+    BasePoints[1] = Base->endpoints[1]->node->getPosition();
+  } else {
+    BasePoints[0] = Base->endpoints[1]->node->getPosition();
+    BasePoints[1] = Base->endpoints[0]->node->getPosition();
+  }
   // check whether everything is in place to create new lines and triangles
   if (i < 4) {
 …
       return NULL;
+    }
+  // construct plane of first triangle for calculating normal vectors later
+  const Plane triangleplane = Base->triangles.begin()->second->getPlane();
+  // get fourth point projected down onto this plane
+  const Vector Intersection = triangleplane.getClosestPoint(OtherEndpoint[0]);
   // remove triangles and baseline removes itself
 …
   m = 0;
   // first obtain all triangle to delete ... (otherwise we pull the carpet (Base) from under the for-loop's feet)
   list <BoundaryTriangleSet *> TrianglesOfBase;
+  list<BoundaryTriangleSet *> TrianglesOfBase;
   for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); ++runner)
     TrianglesOfBase.push_back(runner->second);
   // .. then delete each triangle (which deletes the line as well)
   for (list <BoundaryTriangleSet *>::iterator runner = TrianglesOfBase.begin(); !TrianglesOfBase.empty(); runner = TrianglesOfBase.begin()) {
+  for (list<BoundaryTriangleSet *>::iterator runner = TrianglesOfBase.begin(); !TrianglesOfBase.empty(); runner = TrianglesOfBase.begin()) {
     LOG(3, "DEBUG: Deleting triangle " << *(*runner) << ".");
     OldTriangleNrs[m++] = (*runner)->Nr;
 …
   LinesOnBoundary.insert(LinePair(OldBaseLineNr, NewLine)); // no need for check for unique insertion as NewLine is definitely a new one
   LOG(3, "DEBUG: Created new baseline " << *NewLine << ".");
+  // Explanation for normal vector choice:
+  // Decisive for the normal vector of the new triangle is whether the fourth
+  // endpoint is with respect to the joining boundary line on one side or on
+  // the other with respect to the endpoint of the plane triangle that is not
+  // contained in the joining boundary line.
   // construct new triangles with flipped baseline
 …
     BLS[2] = NewLine;
     BTS = new class BoundaryTriangleSet(BLS, OldTriangleNrs[0]);
+    BTS->GetNormalVector(BaseLineNormal);
+    {
+      Line joiningline = makeLineThrough(
+          OldLines[0]->endpoints[0]->node->getPosition(),
+          OldLines[0]->endpoints[1]->node->getPosition());
+      // BasePoints[1] is not contained in OldLines[0], hence is third endpoint
+      // of plane triangle. BasePoints[0] is in the joining OldLines[0] and
+      // we check whether Intersection is on same side as BasePoints[1] or not.
+      const Vector pointinginside =
+          joiningline.getClosestPoint(BasePoints[1]) - BasePoints[1];
+      const Vector pointingtointersection =
+          joiningline.getClosestPoint(Intersection) - Intersection;
+      const double sign_of_intersection =
+          pointingtointersection.ScalarProduct(pointinginside);
+      LOG(4, "DEBUG: Sign of SKP between both lines w.r.t " << *OldLines[0]
+          << " is " << sign_of_intersection << ".");
+      const double sign_of_normal = (sign_of_intersection >= 0) ? -1. : +1.;
+      LOG(4, "DEBUG: Opposite normal direction is "
+          << sign_of_normal * BaseLineNormal[0] << ".");
+      BTS->GetNormalVector(sign_of_normal * BaseLineNormal[0]);
+    }
     AddTesselationTriangle(OldTriangleNrs[0]);
     LOG(3, "DEBUG: Created new triangle " << *BTS << ".");
 …
     BLS[2] = NewLine;
     BTS = new class BoundaryTriangleSet(BLS, OldTriangleNrs[1]);
+    BTS->GetNormalVector(BaseLineNormal);
+    {
+      Line joiningline = makeLineThrough(
+          OldLines[1]->endpoints[0]->node->getPosition(),
+          OldLines[1]->endpoints[1]->node->getPosition());
+      // BasePoints[0] is not contained in OldLines[1], hence is third endpoint
+      // of plane triangle. BasePoints[1] is in th1e joining OldLines[1] and
+      // we check whether Intersection is on same side as BasePoints[0] or not.
+      const Vector pointinginside =
+          joiningline.getClosestPoint(BasePoints[0]) - BasePoints[0];
+      const Vector pointingtointersection =
+          joiningline.getClosestPoint(Intersection) - Intersection;
+      const double sign_of_intersection =
+          pointingtointersection.ScalarProduct(pointinginside);
+      LOG(4, "DEBUG: Sign of SKP between both lines w.r.t " << *OldLines[1]
+          << " is " << sign_of_intersection << ".");
+      const double sign_of_normal = (sign_of_intersection >= 0) ? -1. : +1.;
+      LOG(4, "DEBUG: Opposite normal direction is "
+          << sign_of_normal * BaseLineNormal[0] << ".");
+      BTS->GetNormalVector(sign_of_normal * BaseLineNormal[0]);
+    }
     AddTesselationTriangle(OldTriangleNrs[1]);
     LOG(3, "DEBUG: Created new triangle " << *BTS << ".");
 …
   TesselPoint *Candidate = NULL;
   LOG(3, "DEBUG: NormalVector of BaseTriangle is " << NormalVector << ".");
+  LOG(4, "DEBUG: NormalVector of BaseTriangle is " << NormalVector << ".");
   // copy old center
 …
   // construct center of circle
+  CircleCenter = 0.5 * ((CandidateLine.BaseLine->endpoints[0]->node->getPosition()) +
+                        (CandidateLine.BaseLine->endpoints[1]->node->getPosition()));
+  CircleCenter = 0.5 * ((CandidateLine.BaseLine->endpoints[0]->node->getPosition()) + (CandidateLine.BaseLine->endpoints[1]->node->getPosition()));
   // construct normal vector of circle
+  CirclePlaneNormal = (CandidateLine.BaseLine->endpoints[0]->node->getPosition()) -
+                      (CandidateLine.BaseLine->endpoints[1]->node->getPosition());
+  CirclePlaneNormal = (CandidateLine.BaseLine->endpoints[0]->node->getPosition()) - (CandidateLine.BaseLine->endpoints[1]->node->getPosition());
   RelativeOldSphereCenter = OldSphereCenter - CircleCenter;
 …
     CircleRadius = RADIUS * RADIUS - radius;
     CirclePlaneNormal.Normalize();
     LOG(3, "DEBUG: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << ".");
+    LOG(4, "DEBUG: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << ".");
     // test whether old center is on the band's plane
 …
     radius = RelativeOldSphereCenter.NormSquared();
     if (fabs(radius - CircleRadius) < HULLEPSILON) {
       LOG(3, "DEBUG: RelativeOldSphereCenter is at " << RelativeOldSphereCenter << ".");
+      LOG(4, "DEBUG: RelativeOldSphereCenter is at " << RelativeOldSphereCenter << ".");
       // check SearchDirection
       LOG(3, "DEBUG: SearchDirection is " << SearchDirection << ".");
+      LOG(4, "DEBUG: SearchDirection is " << SearchDirection << ".");
       if (fabs(RelativeOldSphereCenter.ScalarProduct(SearchDirection)) > HULLEPSILON) { // rotated the wrong way!
         ELOG(1, "SearchDirection and RelativeOldSphereCenter are not orthogonal!");
 …
                   // find center on the plane
                   GetCenterofCircumcircle(NewPlaneCenter, CandidateLine.BaseLine->endpoints[0]->node->getPosition(), CandidateLine.BaseLine->endpoints[1]->node->getPosition(), Candidate->getPosition());
                   LOG(3, "DEBUG: NewPlaneCenter is " << NewPlaneCenter << ".");
+                  LOG(5, "DEBUG: NewPlaneCenter is " << NewPlaneCenter << ".");
                   try {
+                    NewNormalVector = Plane((CandidateLine.BaseLine->endpoints[0]->node->getPosition()),
+                                            (CandidateLine.BaseLine->endpoints[1]->node->getPosition()),
+                                            (Candidate->getPosition())).getNormal();
+                    LOG(3, "DEBUG: NewNormalVector is " << NewNormalVector << ".");
+                    NewNormalVector = Plane((CandidateLine.BaseLine->endpoints[0]->node->getPosition()), (CandidateLine.BaseLine->endpoints[1]->node->getPosition()), (Candidate->getPosition())).getNormal();
+                    LOG(5, "DEBUG: NewNormalVector is " << NewNormalVector << ".");
                     radius = CandidateLine.BaseLine->endpoints[0]->node->DistanceSquared(NewPlaneCenter);
                     LOG(3, "DEBUG: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << ".");
                     LOG(3, "DEBUG: SearchDirection is " << SearchDirection << ".");
                     LOG(3, "DEBUG: Radius of CircumCenterCircle is " << radius << ".");
+                    LOG(5, "DEBUG: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << ".");
+                    LOG(5, "DEBUG: SearchDirection is " << SearchDirection << ".");
+                    LOG(5, "DEBUG: Radius of CircumCenterCircle is " << radius << ".");
                     if (radius < RADIUS * RADIUS) {
                       otherradius = CandidateLine.BaseLine->endpoints[1]->node->DistanceSquared(NewPlaneCenter);
 …
                         // construct both new centers
                         NewSphereCenter = NewPlaneCenter;
                         OtherNewSphereCenter= NewPlaneCenter;
+                        OtherNewSphereCenter = NewPlaneCenter;
                         helper = NewNormalVector;
                         helper.Scale(sqrt(RADIUS * RADIUS - radius));
                         LOG(4, "DEBUG: Distance of NewPlaneCenter " << NewPlaneCenter << " to either NewSphereCenter is " << helper.Norm() << " of vector " << helper << " with sphere radius " << RADIUS << ".");
+                        LOG(5, "DEBUG: Distance of NewPlaneCenter " << NewPlaneCenter << " to either NewSphereCenter is " << helper.Norm() << " of vector " << helper << " with sphere radius " << RADIUS << ".");
                         NewSphereCenter += helper;
                         LOG(4, "DEBUG: NewSphereCenter is at " << NewSphereCenter << ".");
+                        LOG(5, "DEBUG: NewSphereCenter is at " << NewSphereCenter << ".");
                         // OtherNewSphereCenter is created by the same vector just in the other direction
                         helper.Scale(-1.);
                         OtherNewSphereCenter += helper;
                         LOG(4, "DEBUG: OtherNewSphereCenter is at " << OtherNewSphereCenter << ".");
+                        LOG(5, "DEBUG: OtherNewSphereCenter is at " << OtherNewSphereCenter << ".");
                         alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection, HULLEPSILON);
                         Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection, HULLEPSILON);
 …
                           if ((CandidateLine.ShortestAngle - HULLEPSILON) < alpha) {
                             CandidateLine.pointlist.push_back(Candidate);
                             LOG(2, "ACCEPT: We have found an equally good candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << ".");
+                            LOG(3, "ACCEPT: We have found an equally good candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << ".");
                           } else {
                             // remove all candidates from the list and then the list itself
                             CandidateLine.pointlist.clear();
                             CandidateLine.pointlist.push_back(Candidate);
                             LOG(2, "ACCEPT: We have found a better candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << ".");
+                            LOG(3, "ACCEPT: We have found a better candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << ".");
+                          }
                           CandidateLine.ShortestAngle = alpha;
                           LOG(2, "DEBUG: There are " << CandidateLine.pointlist.size() << " candidates in the list now.");
+                          LOG(4, "DEBUG: There are " << CandidateLine.pointlist.size() << " candidates in the list now.");
                         } else {
                           if ((Candidate != NULL) && (CandidateLine.pointlist.begin() != CandidateLine.pointlist.end())) {
                             LOG(3, "REJECT: Old candidate " << *(*CandidateLine.pointlist.begin()) << " with " << CandidateLine.ShortestAngle << " is better than new one " << *Candidate << " with " << alpha << " .");
+                            LOG(4, "REJECT: Old candidate " << *(*CandidateLine.pointlist.begin()) << " with " << CandidateLine.ShortestAngle << " is better than new one " << *Candidate << " with " << alpha << " .");
                           } else {
                             LOG(3, "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected.");
+                            LOG(4, "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected.");
+                          }
+                        }
                       } else {
                         ELOG(0, "REJECT: Distance to center of circumcircle is not the same from each corner of the triangle: " << fabs(radius - otherradius));
+                        ASSERT(0, std::string("FindThirdPointForTesselation() - ") + std::string("REJECT: Distance to center of circumcircle is not the same from each corner of the triangle: ") + toString(fabs(radius - otherradius)));
+                      }
                     } else {
                       LOG(3, "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << ".");
+                      LOG(4, "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << ".");
+                    }
+                  }
+                  catch (LinearDependenceException &excp){
+                    LOG(3, boost::diagnostic_information(excp));
+                    LOG(3, "REJECT: Three points from " << *CandidateLine.BaseLine << " and Candidate " << *Candidate << " are linear-dependent.");
+                  } catch (LinearDependenceException &excp) {
+                    LOG(4, boost::diagnostic_information(excp));
+                    LOG(4, "REJECT: Three points from " << *CandidateLine.BaseLine << " and Candidate " << *Candidate << " are linear-dependent.");
+                  }
                 } else {
                   if (ThirdPoint != NULL) {
                     LOG(3, "REJECT: Base triangle " << *CandidateLine.BaseLine << " and " << *ThirdPoint << " contains Candidate " << *Candidate << ".");
+                    LOG(4, "REJECT: Base triangle " << *CandidateLine.BaseLine << " and " << *ThirdPoint << " contains Candidate " << *Candidate << ".");
                   } else {
                     LOG(3, "REJECT: Base triangle " << *CandidateLine.BaseLine << " contains Candidate " << *Candidate << ".");
+                    LOG(4, "REJECT: Base triangle " << *CandidateLine.BaseLine << " contains Candidate " << *Candidate << ".");
+                  }
+                }
 …
   } else {
     if (ThirdPoint != NULL)
       LOG(3, "Circumcircle for base line " << *CandidateLine.BaseLine << " and third node " << *ThirdPoint << " is too big!");
+      LOG(4, "DEBUG: Circumcircle for base line " << *CandidateLine.BaseLine << " and third node " << *ThirdPoint << " is too big!");
     else
       LOG(3, "Circumcircle for base line " << *CandidateLine.BaseLine << " is too big!");
+  }
   LOG(2, "DEBUG: Sorting candidate list ...");
+      LOG(4, "DEBUG: Circumcircle for base line " << *CandidateLine.BaseLine << " is too big!");
+  }
+  LOG(5, "DEBUG: Sorting candidate list ...");
   if (CandidateLine.pointlist.size() > 1) {
     CandidateLine.pointlist.unique();
 …
+  }
+  if ((!CandidateLine.pointlist.empty()) && (!CandidateLine.CheckValidity(RADIUS, LC))) {
+    ELOG(0, "There were other points contained in the rolling sphere as well!");
+    performCriticalExit();
+  }
+  ASSERT(CandidateLine.pointlist.empty() || (CandidateLine.CheckValidity(RADIUS, LC)), std::string("Tesselation::FindThirdPointForTesselation()") + std::string("There were other points contained in the rolling sphere as well!"));
+}
+;
 …
+{
   //Info FunctionInfo(__func__);
   const BoundaryLineSet * lines[2] = { line1, line2 };
+  const BoundaryLineSet * lines[2] = {line1, line2};
   class BoundaryPointSet *node = NULL;
   PointMap OrderMap;
 …
     // for both lines
     for (int j = 0; j < 2; j++) { // for both endpoints
       OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *> (lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j]));
+      OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *>(lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j]));
       if (!OrderTest.second) { // if insertion fails, we have common endpoint
         node = OrderTest.first->second;
 …
               // we should make sure that both triangles end up in the same entry
               // in the distance multimap. Hence, we round to 6 digit precision.
+              const double distance =
+e-6*floor(FindPoint->second->node->DistanceSquared(x)*1e+6);
+              const double distance = 1e-6 * floor(FindPoint->second->node->DistanceSquared(x) * 1e+6);
               points->insert(DistanceToPointPair(distance, FindPoint->second));
               LOG(3, "DEBUG: Putting " << *FindPoint->second << " into the list.");
 …
     for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
       // calculate closest point on line to desired point
+      helper = 0.5 * (((LineRunner->second)->endpoints[0]->node->getPosition()) +
+                      ((LineRunner->second)->endpoints[1]->node->getPosition()));
+      helper = 0.5 * (((LineRunner->second)->endpoints[0]->node->getPosition()) + ((LineRunner->second)->endpoints[1]->node->getPosition()));
       Center = (x) - helper;
+      BaseLine = ((LineRunner->second)->endpoints[0]->node->getPosition()) -
+                 ((LineRunner->second)->endpoints[1]->node->getPosition());
+      BaseLine = ((LineRunner->second)->endpoints[0]->node->getPosition()) - ((LineRunner->second)->endpoints[1]->node->getPosition());
       Center.ProjectOntoPlane(BaseLine);
       const double distance = Center.NormSquared();
 …
     for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
+      BaseLine = ((LineRunner->second)->endpoints[0]->node->getPosition()) -
+                 ((LineRunner->second)->endpoints[1]->node->getPosition());
+      BaseLine = ((LineRunner->second)->endpoints[0]->node->getPosition()) - ((LineRunner->second)->endpoints[1]->node->getPosition());
       const double lengthBase = BaseLine.NormSquared();
 …
           MinDistance = lengthEnd;
           LOG(1, "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[0]->node << " is closer with " << lengthEnd << ".");
+        } else if (fabs(lengthEnd - MinDistance) < MYEPSILON) { // additional best candidate
+          ClosestLines.insert(LineRunner->second);
+          LOG(1, "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[1]->node << " is equally good with " << lengthEnd << ".");
+        } else { // line is worse
+          LOG(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 << ".");
+        }
+        } else
+          if (fabs(lengthEnd - MinDistance) < MYEPSILON) { // additional best candidate
+            ClosestLines.insert(LineRunner->second);
+            LOG(1, "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[1]->node << " is equally good with " << lengthEnd << ".");
+          } else { // line is worse
+            LOG(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 << ".");
+          }
       } else { // intersection is closer, calculate
         // calculate closest point on line to desired point
 …
   double distance = 0.;
   if (triangle == NULL) {// is boundary point or only point in point cloud?
+  if (triangle == NULL) { // is boundary point or only point in point cloud?
     LOG(1, "No triangle given!");
     return -1.;
 …
       takePoint = true;
       current = findLines->second->endpoints[1]->node;
+    } else if (findLines->second->endpoints[1]->Nr == Point->getNr()) {
+      takePoint = true;
+      current = findLines->second->endpoints[0]->node;
+    }
+    } else
+      if (findLines->second->endpoints[1]->Nr == Point->getNr()) {
+        takePoint = true;
+        current = findLines->second->endpoints[0]->node;
+      }
     if (takePoint) {
 …
   Vector OrthogonalVector;
   Vector helper;
   const TesselPoint * const TrianglePoints[3] = { Point, NULL, NULL };
+  const TesselPoint * const TrianglePoints[3] = {Point, NULL, NULL};
   TriangleList *triangles = NULL;
 …
   // calculate central point
   triangles = FindTriangles(TrianglePoints);
+  if ((triangles != NULL) && (!triangles->empty())) {
+    for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
+      PlaneNormal += (*Runner)->NormalVector;
+  } else {
+    ELOG(0, "Could not find any triangles for point " << *Point << ".");
+    performCriticalExit();
+  }
+  ASSERT((triangles == NULL) || (triangles->empty()), std::string("Tesselation::GetCircleOfConnectedTriangles()") + std::string("Could not find any triangles for point " + toString(*Point) + "."));
+  for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
+    PlaneNormal += (*Runner)->NormalVector;
   PlaneNormal.Scale(1.0 / triangles->size());
   LOG(4, "DEBUG: Calculated PlaneNormal of all circle points is " << PlaneNormal << ".");
 …
     AngleZero = ((*SetOfNeighbours->begin())->getPosition()) - (Point->getPosition());
     AngleZero.ProjectOntoPlane(PlaneNormal);
+    if (AngleZero.NormSquared() < MYEPSILON) {
+      ELOG(0, "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!");
+      performCriticalExit();
+    }
+    ASSERT(AngleZero.NormSquared() > MYEPSILON, std::string("Tesselation::GetCircleOfConnectedTriangles() - ") + std::string("AngleZero is 0 even with alternative reference.") + std::string("The algorithm has to be changed here!"));
+  }
   LOG(4, "DEBUG: Reference vector on this plane representing angle 0 is " << AngleZero << ".");
   if (AngleZero.NormSquared() > MYEPSILON)
     OrthogonalVector = Plane(PlaneNormal, AngleZero,0).getNormal();
+    OrthogonalVector = Plane(PlaneNormal, AngleZero, 0).getNormal();
   else
     OrthogonalVector.MakeNormalTo(PlaneNormal);
 …
     double angle = GetAngle(helper, AngleZero, OrthogonalVector);
     LOG(4, "DEBUG" << angle << " for point " << **listRunner << ".");
     anglesOfPoints.insert(pair<double, TesselPoint*> (angle, (*listRunner)));
+    anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
+  }
 …
   int counter = 0;
   while (TesselC != SetOfNeighbours->end()) {
+    helper = Plane(((*TesselA)->getPosition()),
+                   ((*TesselB)->getPosition()),
+                   ((*TesselC)->getPosition())).getNormal();
+    helper = Plane(((*TesselA)->getPosition()), ((*TesselB)->getPosition()), ((*TesselC)->getPosition())).getNormal();
     LOG(5, "DEBUG: Making normal vector out of " << *(*TesselA) << ", " << *(*TesselB) << " and " << *(*TesselC) << ":" << helper);
     counter++;
 …
+  }
   //LOG(0, "Summed vectors " << center << "; number of points " << connectedPoints.size() << "; scale factor " << counter);
   PlaneNormal.Scale(1.0 / (double) counter);
+  PlaneNormal.Scale(1.0 / (double)counter);
   //  LOG(1, "INFO: Calculated center of all circle points is " << center << ".");
   //
 …
     AngleZero.ProjectOntoPlane(PlaneNormal);
+  }
   if ((Reference.IsZero()) || (AngleZero.NormSquared() < MYEPSILON )) {
+  if ((Reference.IsZero()) || (AngleZero.NormSquared() < MYEPSILON)) {
     LOG(4, "DEBUG: Using alternatively " << (*SetOfNeighbours->begin())->getPosition() << " as angle 0 referencer.");
     AngleZero = ((*SetOfNeighbours->begin())->getPosition()) - (Point->getPosition());
     AngleZero.ProjectOntoPlane(PlaneNormal);
+    if (AngleZero.NormSquared() < MYEPSILON) {
+      ELOG(0, "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!");
+      performCriticalExit();
+    }
+    ASSERT(AngleZero.NormSquared() > MYEPSILON, std::string("Tesselation::GetCircleOfSetOfPoints() - ") + std::string("AngleZero is 0 even with alternative reference.") + std::string("The algorithm has to be changed here!"));
+  }
   LOG(4, "DEBUG: Reference vector on this plane representing angle 0 is " << AngleZero << ".");
   if (AngleZero.NormSquared() > MYEPSILON)
     OrthogonalVector = Plane(PlaneNormal, AngleZero,0).getNormal();
+    OrthogonalVector = Plane(PlaneNormal, AngleZero, 0).getNormal();
   else
     OrthogonalVector.MakeNormalTo(PlaneNormal);
 …
       angle = 2. * M_PI - angle;
     LOG(4, "DEBUG: Calculated angle between " << helper << " and " << AngleZero << " is " << angle << " for point " << **listRunner << ".");
+    InserterTest = anglesOfPoints.insert(pair<double, TesselPoint*> (angle, (*listRunner)));
+    if (!InserterTest.second) {
+      ELOG(0, "GetCircleOfSetOfPoints() got two atoms with same angle: " << *((InserterTest.first)->second) << " and " << (*listRunner));
+      performCriticalExit();
+    }
+    InserterTest = anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
+    ASSERT(InserterTest.second, std::string("Tesselation::GetCircleOfSetOfPoints() - ") + std::string("got two atoms with same angle " + toString(*((InserterTest.first)->second))) + std::string(" and " + toString((*listRunner))));
+  }
 …
   //Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
   list<TesselPointList *> *ListOfPaths = new list<TesselPointList *> ;
+  list<TesselPointList *> *ListOfPaths = new list<TesselPointList *>;
   TesselPointList *connectedPath = NULL;
   Vector center;
 …
   map<class BoundaryTriangleSet *, bool>::iterator TriangleRunner;
   for (LineMap::iterator Runner = ReferencePoint->lines.begin(); Runner != ReferencePoint->lines.end(); Runner++) {
+    TouchedLine.insert(pair<class BoundaryLineSet *, bool> (Runner->second, false));
+    for (TriangleMap::iterator Sprinter = Runner->second->triangles.begin(); Sprinter != Runner->second->triangles.end(); Sprinter++)
+      TouchedTriangle.insert(pair<class BoundaryTriangleSet *, bool> (Sprinter->second, false));
+    LOG(4, "DEBUG: Adding " << *Runner->second << " to TouchedLine map.");
+    TouchedLine.insert(pair<class BoundaryLineSet *, bool>(Runner->second, false));
+    for (TriangleMap::iterator Sprinter = Runner->second->triangles.begin(); Sprinter != Runner->second->triangles.end(); Sprinter++) {
+      LOG(4, "DEBUG: Adding " << *Sprinter->second << " to TouchedTriangle map.");
+      TouchedTriangle.insert(pair<class BoundaryTriangleSet *, bool>(Sprinter->second, false));
+    }
+  }
   if (!ReferencePoint->lines.empty()) {
 …
       if (LineRunner == TouchedLine.end()) {
         ELOG(1, "I could not find " << *runner->second << " in the touched list.");
+      } else if (!LineRunner->second) {
+        LineRunner->second = true;
+        connectedPath = new TesselPointList;
+        triangle = NULL;
+        CurrentLine = runner->second;
+        StartLine = CurrentLine;
+        CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint);
+        LOG(1, "INFO: Beginning path retrieval at " << *CurrentPoint << " of line " << *CurrentLine << ".");
+        do {
+          // push current one
+          LOG(1, "INFO: Putting " << *CurrentPoint << " at end of path.");
+          connectedPath->push_back(CurrentPoint->node);
+          // find next triangle
+          for (TriangleMap::iterator Runner = CurrentLine->triangles.begin(); Runner != CurrentLine->triangles.end(); Runner++) {
+            LOG(1, "INFO: Inspecting triangle " << *Runner->second << ".");
+            if ((Runner->second != triangle)) { // look for first triangle not equal to old one
+              triangle = Runner->second;
+              TriangleRunner = TouchedTriangle.find(triangle);
+              if (TriangleRunner != TouchedTriangle.end()) {
+                if (!TriangleRunner->second) {
+                  TriangleRunner->second = true;
+                  LOG(1, "INFO: Connecting triangle is " << *triangle << ".");
+                  break;
+      } else
+        if (!LineRunner->second) {
+          LineRunner->second = true;
+          connectedPath = new TesselPointList;
+          triangle = NULL;
+          CurrentLine = runner->second;
+          StartLine = CurrentLine;
+          CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint);
+          LOG(3, "INFO: Beginning path retrieval at " << *CurrentPoint << " of line " << *CurrentLine << ".");
+          do {
+            // push current one
+            LOG(3, "INFO: Putting " << *CurrentPoint << " at end of path.");
+            connectedPath->push_back(CurrentPoint->node);
+            // find next triangle
+            for (TriangleMap::iterator Runner = CurrentLine->triangles.begin(); Runner != CurrentLine->triangles.end(); Runner++) {
+              LOG(4, "DEBUG: Inspecting triangle " << *Runner->second << ".");
+              if ((Runner->second != triangle)) { // look for first triangle not equal to old one
+                triangle = Runner->second;
+                TriangleRunner = TouchedTriangle.find(triangle);
+                if (TriangleRunner != TouchedTriangle.end()) {
+                  if (!TriangleRunner->second) {
+                    TriangleRunner->second = true;
+                    LOG(4, "DEBUG: Connecting triangle is " << *triangle << ".");
+                    break;
+                  } else {
+                    LOG(4, "DEBUG: Skipping " << *triangle << ", as we have already visited it.");
+                    triangle = NULL;
+                  }
                 } else {
                   LOG(1, "INFO: Skipping " << *triangle << ", as we have already visited it.");
+                  ELOG(1, "I could not find " << *triangle << " in the touched list.");
                   triangle = NULL;
+                }
               } else {
                 ELOG(1, "I could not find " << *triangle << " in the touched list.");
+                // as we have stumbled upon the same triangle, we don't need the check anymore
                 triangle = NULL;
+              }
+            }
+            if (triangle == NULL)
+              break;
+            // find next line
+            for (int i = 0; i < 3; i++) {
+              if ((triangle->lines[i] != CurrentLine) && (triangle->lines[i]->ContainsBoundaryPoint(ReferencePoint))) { // not the current line and still containing Point
+                CurrentLine = triangle->lines[i];
+                LOG(3, "INFO: Connecting line is " << *CurrentLine << ".");
+                break;
+              }
+            }
+            LineRunner = TouchedLine.find(CurrentLine);
+            if (LineRunner == TouchedLine.end())
+              ELOG(1, "I could not find " << *CurrentLine << " in the touched list.");
+            else
+              LineRunner->second = true;
+            // find next point
+            CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint);
+          } while (CurrentLine != StartLine);
+          // last point is missing, as it's on start line
+          if (StartLine->GetOtherEndpoint(ReferencePoint)->node != connectedPath->back()) {
+            LOG(3, "INFO: Putting " << *CurrentPoint << " at end of path to close it.");
+            connectedPath->push_back(StartLine->GetOtherEndpoint(ReferencePoint)->node);
+          }
+          if (triangle == NULL)
+            break;
+          // find next line
+          for (int i = 0; i < 3; i++) {
+            if ((triangle->lines[i] != CurrentLine) && (triangle->lines[i]->ContainsBoundaryPoint(ReferencePoint))) { // not the current line and still containing Point
+              CurrentLine = triangle->lines[i];
+              LOG(1, "INFO: Connecting line is " << *CurrentLine << ".");
+              break;
+            }
+          }
+          LineRunner = TouchedLine.find(CurrentLine);
+          if (LineRunner == TouchedLine.end())
+            ELOG(1, "I could not find " << *CurrentLine << " in the touched list.");
+          else
+            LineRunner->second = true;
+          // find next point
+          CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint);
+        } while (CurrentLine != StartLine);
+        // last point is missing, as it's on start line
+        LOG(1, "INFO: Putting " << *CurrentPoint << " at end of path.");
+        if (StartLine->GetOtherEndpoint(ReferencePoint)->node != connectedPath->back())
+          connectedPath->push_back(StartLine->GetOtherEndpoint(ReferencePoint)->node);
+        ListOfPaths->push_back(connectedPath);
+      } else {
+        LOG(1, "INFO: Skipping " << *runner->second << ", as we have already visited it.");
+      }
+          ListOfPaths->push_back(connectedPath);
+        } else {
+          LOG(3, "DEBUG: Skipping " << *runner->second << ", as we have already visited it.");
+        }
+    }
   } else {
 …
   //Info FunctionInfo(__func__);
   list<TesselPointList *> *ListofPaths = GetPathsOfConnectedPoints(Point);
   list<TesselPointList *> *ListofClosedPaths = new list<TesselPointList *> ;
+  list<TesselPointList *> *ListofClosedPaths = new list<TesselPointList *>;
   TesselPointList *connectedPath = NULL;
   TesselPointList *newPath = NULL;
 …
     connectedPath = *ListRunner;
+    LOG(1, "INFO: Current path is " << connectedPath << ".");
+    if (DoLog(2)) {
+      std::stringstream output;
+      output << "INFO: Current path is ";
+      BOOST_FOREACH(const TesselPoint * const item, *connectedPath) {
+        output << *item << " ";
+      }
+      LOG(1, output.str());
+    }
     // go through list, look for reappearance of starting Point and count
 …
       if ((*CircleRunner == *CircleStart) && (CircleRunner != CircleStart)) { // is not the very first point
         // we have a closed circle from Marker to new Marker
         if (DoLog(1)) {
+        if (DoLog(3)) {
           std::stringstream output;
+          output << count + 1 << ". closed path consists of: ";
+          for (TesselPointList::iterator CircleSprinter = Marker;
+              CircleSprinter != CircleRunner;
+              CircleSprinter++)
+          output << "DEBUG: " << count + 1 << ". closed path consists of: ";
+          for (TesselPointList::iterator CircleSprinter = Marker; CircleSprinter != CircleRunner; CircleSprinter++)
             output << (**CircleSprinter) << " <-> ";
           LOG(1, output.str());
 …
+    }
+  }
   LOG(1, "INFO: " << count << " closed additional path(s) have been created.");
+  LOG(2, "DEBUG: " << count << " closed additional path(s) have been created.");
   // delete list of paths
 …
+}
+;
+struct CloserToPiHalf
+{
+  bool operator()(double angle, double smallestangle)
+  {
+    return (fabs(angle - M_PI / 2.) < fabs(smallestangle - M_PI / 2.));
+  }
+};
+bool Tesselation::IsPointBelowSurroundingPolygon(const BoundaryPointSet *_point) const
+{
+  // check for NULL
+  if (_point == NULL) {
+    return false;
+  }
+  // get list of connected points
+  if (_point->lines.empty()) {
+    LOG(1, "INFO: point " << *_point << " is not connected to any lines.");
+    return false;
+  }
+  bool PointIsBelow = true;
+  // create Orthogonal vector as reference for angle (pointing into [pi,2pi) interval)
+  Vector OrthogonalVector;
+  for (LineMap::const_iterator lineiter = _point->lines.begin();
+      lineiter != _point->lines.end(); ++lineiter)
+    for (TriangleMap::const_iterator triangleiter = lineiter->second->triangles.begin();
+        triangleiter != lineiter->second->triangles.end();
+        ++triangleiter)
+        OrthogonalVector += triangleiter->second->NormalVector;
+  OrthogonalVector.Normalize();
+  // go through all closed paths for this point
+  typedef list<TesselPointList *> TPL_list_t;
+  const TPL_list_t *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(_point->node);
+  for (TPL_list_t::const_iterator closedpathsiter = ListOfClosedPaths->begin();
+      (closedpathsiter != ListOfClosedPaths->end()) && PointIsBelow;
+      ++closedpathsiter) {
+    const TesselPointList *connectedPath = *closedpathsiter;
+    TesselPointList::const_iterator ListAdvance = connectedPath->begin(); // gives angle direction
+    TesselPointList::const_iterator ListRunner = ListAdvance++;
+    for (; (ListAdvance != connectedPath->end()) && PointIsBelow;
+        ListRunner = ListAdvance++) { // go through all closed paths
+      LOG(2, "DEBUG: Current reference node is " << **ListRunner
+          << ", advanced node is " << **ListAdvance);
+      // reference vector to point to check for this point of connected path
+      const Vector Reference =
+          ((*ListAdvance)->getPosition()) - ((*ListRunner)->getPosition());
+      // go through all other points in this connected path
+      for (TesselPointList::const_iterator OtherRunner = ListAdvance;
+          (OtherRunner != ListRunner) && PointIsBelow;
+          ++OtherRunner == connectedPath->end() ?
+              OtherRunner = connectedPath->begin() :
+              OtherRunner) {
+        if (OtherRunner == ListAdvance)
+          continue;
+        LOG(3, "DEBUG: Current other node is " << **OtherRunner);
+        // build the plane with normal vector
+        const Vector onebeam =
+            ((*OtherRunner)->getPosition()) - ((*ListAdvance)->getPosition());
+        Vector NormalVector = Reference;
+        NormalVector.VectorProduct(onebeam);
+        // needs to point in same general direction as average NormalVector of all triangles
+        if (NormalVector.ScalarProduct(OrthogonalVector) < 0)
+          NormalVector *= -1.;
+        try {
+          Plane plane(NormalVector, ((*ListRunner)->getPosition()));
+          // check whether point is below
+          if (plane.SignedDistance(_point->node->getPosition()) > 0) {
+            LOG(2, "DEBUG: For plane " << plane << " point " << *_point
+                << " would remain above.");
+            PointIsBelow = false;
+          }
+        } catch (ZeroVectorException &e) {
+          ELOG(3, "Vectors are linear dependent, skipping.");
+        }
+      }
+    }
+  }
+  delete ListOfClosedPaths;
+  return PointIsBelow;
+}
+double Tesselation::RemovePointSurroundedByPolygon(
+    TesselPointList *connectedPath,
+    BoundaryPointSet *point)
+{
+  double volume = 0.;
+  const Vector OldPoint = point->node->getPosition();
+  TesselPoint *oldNode = point->node;
+  // remove present triangles for this connectedPath
+  unsigned int count = 0;
+  for (TesselPointList::const_iterator iter = connectedPath->begin();
+      iter != connectedPath->end(); ++iter)
+    LOG(4, "DEBUG: Node in connectedPath is " << **iter);
+  {
+    TesselPointList::iterator FirstNode, SecondNode;
+    SecondNode = connectedPath->begin();
+    FirstNode = SecondNode++;
+    for (;FirstNode != connectedPath->end(); ++SecondNode, ++FirstNode) {
+      LOG(3, "DEBUG: MiddleNode is " << **FirstNode << ".");
+      if (SecondNode == connectedPath->end())
+        SecondNode = connectedPath->begin();
+      TesselPoint *TriangleCandidates[3];
+      TriangleCandidates[0] = *FirstNode;
+      TriangleCandidates[1] = *SecondNode;
+      TriangleCandidates[2] = oldNode;
+      BoundaryTriangleSet *triangle = GetPresentTriangle(TriangleCandidates);
+      ASSERT( triangle != NULL,
+          "Tesselation::RemovePointSurroundedByPolygon() - triangle to points "
+          +toString((**SecondNode))+", "+toString((**FirstNode))+" and "
+          +toString(*oldNode)+" does not exist.");
+      LOG(3, "DEBUG: Removing triangle " << *triangle << ".");
+      RemoveTesselationTriangle(triangle);
+      ++count;
+    }
+    LOG(2, "INFO: " << count << " triangles were removed.");
+  }
+  // re-create all triangles by going through connected points list
+  LineList NewLines;
+  typedef std::vector<double> angles_t;
+  angles_t angles;
+  count = 0;
+  for (; !connectedPath->empty();) {
+    // search middle node with widest angle to next neighbours
+    TesselPointList::iterator StartNode, MiddleNode, EndNode;
+    for (MiddleNode = connectedPath->begin(); MiddleNode != connectedPath->end(); MiddleNode++) {
+      LOG(3, "INFO: MiddleNode is " << **MiddleNode << ".");
+      // construct vectors to next and previous neighbour
+      StartNode = MiddleNode;
+      if (StartNode == connectedPath->begin())
+        StartNode = connectedPath->end();
+      StartNode--;
+      //LOG(3, "INFO: StartNode is " << **StartNode << ".");
+      const Vector Point = ((*StartNode)->getPosition()) - ((*MiddleNode)->getPosition());
+      EndNode = MiddleNode;
+      EndNode++;
+      if (EndNode == connectedPath->end())
+        EndNode = connectedPath->begin();
+      //LOG(3, "INFO: EndNode is " << **StartNode << ".");
+      const Vector Reference = ((*EndNode)->getPosition()) - ((*MiddleNode)->getPosition());
+      Vector OrthogonalVector = ((*MiddleNode)->getPosition()) - OldPoint;
+      OrthogonalVector.MakeNormalTo(Reference);
+      angles.push_back(GetAngle(Point, Reference, OrthogonalVector));
+    }
+    ASSERT( !angles.empty(),
+        "Tesselation::RemovePointSurroundedByPolygon() - angles empty");
+    const angles_t::const_iterator maxiter = std::max_element(angles.begin(), angles.end());
+    angles_t::const_iterator miniter = angles.begin();
+    // distinguish between convex and nonconvex polygon
+    if (*maxiter > M_PI) {
+      // connectedPath is not convex: The idea is to fill any kinks pointing
+      // inside into the connectedPath close to this concave spot first, making
+      // it eventually become convex.
+      // Hence, use adjacent (and convex) fill-in point
+      miniter = (maxiter == angles.begin()) ? angles.end()-1 : maxiter-1;
+      if (*miniter > M_PI) {
+        miniter = (maxiter+1 == angles.end()) ? angles.begin() : maxiter+1;
+        if (*miniter > M_PI) {
+          miniter = std::min_element(angles.begin(), angles.end());
+        }
+      }
+    } else {
+      // is convex
+      miniter = std::min_element(angles.begin(), angles.end(), CloserToPiHalf());
+    }
+    MiddleNode = connectedPath->begin();
+    std::advance(MiddleNode, std::distance(const_cast<const angles_t &>(angles).begin(), miniter));
+    ASSERT(MiddleNode != connectedPath->end(),
+        "Tesselation::RemovePointSurroundedByPolygon() - Could not find a smallest angle!");
+    angles.clear();
+    StartNode = MiddleNode;
+    EndNode = MiddleNode;
+    if (StartNode == connectedPath->begin())
+      StartNode = connectedPath->end();
+    StartNode--;
+    EndNode++;
+    if (EndNode == connectedPath->end())
+      EndNode = connectedPath->begin();
+    LOG(2, "INFO: StartNode is " << **StartNode << ".");
+    LOG(2, "INFO: MiddleNode is " << **MiddleNode << ".");
+    LOG(2, "INFO: EndNode is " << **EndNode << ".");
+    LOG(1, "INFO: Attempting to create triangle " << (*StartNode)->getName() << ", " << (*MiddleNode)->getName() << " and " << (*EndNode)->getName() << ".");
+    TesselPoint *TriangleCandidates[3];
+    TriangleCandidates[0] = *StartNode;
+    TriangleCandidates[1] = *MiddleNode;
+    TriangleCandidates[2] = *EndNode;
+    BoundaryTriangleSet *triangle = GetPresentTriangle(TriangleCandidates);
+    if (triangle != NULL) {
+      const Vector center = 1./3. * ((*StartNode)->getPosition()
+              + (*MiddleNode)->getPosition()
+              + (*EndNode)->getPosition());
+      const Vector NormalVector = OldPoint - center;
+      // check orientation of normal vector (that points inside)
+      ASSERT( triangle->NormalVector.ScalarProduct(NormalVector) > std::numeric_limits<double>::epsilon()*1e2,
+          "Tesselation::RemovePointSurroundedByPolygon() - New triangle with same orientation already present as "
+          +toString(*triangle)+"!");
+    }
+    LOG(3, "DEBUG: Adding new triangle points.");
+    AddTesselationPoint(*StartNode, 0);
+    AddTesselationPoint(*MiddleNode, 1);
+    AddTesselationPoint(*EndNode, 2);
+    LOG(3, "DEBUG: Adding new triangle lines.");
+    AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
+    // line between start and end must be new (except for very last triangle)
+    if (AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1))
+      NewLines.push_back(BLS[1]);
+    AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
+    BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+    const Vector center = 1./3. * ((*StartNode)->getPosition()
+            + (*MiddleNode)->getPosition()
+            + (*EndNode)->getPosition());
+    const Vector NormalVector = OldPoint - center;
+    BTS->GetNormalVector(NormalVector);
+    AddTesselationTriangle();
+    // calculate volume summand as a general tetraeder
+    volume += CalculateVolumeofGeneralTetraeder(
+        TPS[0]->node->getPosition(),
+        TPS[1]->node->getPosition(),
+        TPS[2]->node->getPosition(),
+        OldPoint);
+    // advance number
+    ++count;
+    // prepare nodes for next triangle
+    LOG(3, "DEBUG: Removing " << **MiddleNode << " from closed path.");
+    connectedPath->remove(*MiddleNode); // remove the middle node (it is surrounded by triangles)
+    LOG(3, "DEBUG: Remaining points: " << connectedPath->size() << ".");
+    ASSERT(connectedPath->size() >= 2,
+        "Tesselation::RemovePointSurroundedByPolygon() - There are only two endpoints left!");
+    if (connectedPath->size() == 2) { // we are done
+      connectedPath->remove(*StartNode); // remove the start node
+      connectedPath->remove(*EndNode); // remove the end node
+    }
+  }
+  LOG(1, "INFO: " << count << " triangles were created.");
+  return volume;
+}
 /** Removes a boundary point from the envelope while keeping it closed.
 …
 double Tesselation::RemovePointFromTesselatedSurface(class BoundaryPointSet *point)
+{
-  class BoundaryLineSet *line = NULL;
-  class BoundaryTriangleSet *triangle = NULL;
-  Vector OldPoint, NormalVector;
   double volume = 0;
-  int count = 0;
   if (point == NULL) {
 …
     LOG(4, "DEBUG: Removing point " << *point << " from tesselated boundary ...");
-  // copy old location for the volume
-  OldPoint = (point->node->getPosition());
   // get list of connected points
   if (point->lines.empty()) {
 …
   list<TesselPointList *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
+  TesselPointList *connectedPath = NULL;
+  // gather all triangles
+  for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++)
+    count += LineRunner->second->triangles.size();
+  TriangleMap Candidates;
+  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;
+      Candidates.insert(TrianglePair(triangle->Nr, triangle));
+    }
+  }
+  // remove all triangles
+  count = 0;
+  NormalVector.Zero();
+  for (TriangleMap::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
+    LOG(1, "INFO: Removing triangle " << *(Runner->second) << ".");
+    NormalVector -= Runner->second->NormalVector; // has to point inward
+    RemoveTesselationTriangle(Runner->second);
+    count++;
+  }
+  LOG(1, count << " triangles were removed.");
+  list<TesselPointList *>::iterator ListAdvance = ListOfClosedPaths->begin();
+  list<TesselPointList *>::iterator ListRunner = ListAdvance;
+//  TriangleMap::iterator NumberRunner = Candidates.begin();
+  Vector Point, Reference, OrthogonalVector;
+  for (; ListRunner != ListOfClosedPaths->end(); ListRunner = ListAdvance) { // go through all closed paths
+    if (ListAdvance != ListOfClosedPaths->end())
+      ListAdvance++;
+    TesselPointList *connectedPath = *ListRunner;
+    volume += RemovePointSurroundedByPolygon(connectedPath, point);
+    ListOfClosedPaths->remove(connectedPath);
+    delete (connectedPath);
+  }
+  delete (ListOfClosedPaths);
+  LOG(1, "INFO: Removed volume is " << volume << ".");
+  return volume;
+}
+;
+bool Tesselation::CheckAllConcaveInPolygon(
+    const TesselPointList *connectedPath,
+    const BoundaryPointSet *point
+    )
+{
+  // check whether lines in this path to point to remove are all concave
+  bool all_lines_concave = true;
+  if (connectedPath->size() >= 2) {
+    TesselPointList::const_iterator StartNode, MiddleNode, EndNode;
+    // have nearest neighbors to a middle node to know adjacent triangles
+    for (MiddleNode = connectedPath->begin();
+        all_lines_concave && (MiddleNode != connectedPath->end());
+        MiddleNode++) {
+      LOG(3, "INFO: MiddleNode is " << **MiddleNode << ".");
+      EndNode = MiddleNode;
+      if (EndNode == connectedPath->begin())
+        EndNode = connectedPath->end();
+      --EndNode;
+      StartNode = MiddleNode;
+      ++StartNode;
+      if (StartNode == connectedPath->end())
+        StartNode = connectedPath->begin();
+      AddTesselationPoint(point->node, 0);
+      AddTesselationPoint(*MiddleNode, 1);
+      ASSERT( point->lines.find((*MiddleNode)->getNr()) != point->lines.end(),
+          "Tesselation::CheckAllConcaveInPolygon() - MiddleNode "
+          +toString(**MiddleNode)+" not present in "
+          +toString(*point)+"'s lines.");
+      // get line between Node and point and check
+      std::pair<LineMap::const_iterator, LineMap::const_iterator> FindPair =
+          point->lines.equal_range((*MiddleNode)->getNr());
+      LineMap::const_iterator FindLine = FindPair.first;
+      for (; FindLine != FindPair.second; ++FindLine) {
+        // line  has got two triangles, check whether they resemble those
+        // with start and endnode
+        const BoundaryLineSet *currentline = FindLine->second;
+        unsigned int matching_triangles = 0;
+        for (TriangleMap::const_iterator triangleiter = currentline->triangles.begin();
+            triangleiter != currentline->triangles.end(); ++triangleiter) {
+          const BoundaryTriangleSet *triangle = triangleiter->second;
+          AddTesselationPoint(*StartNode, 2);
+          if (triangle->IsPresentTupel(TPS))
+            ++matching_triangles;
+          AddTesselationPoint(*EndNode, 2);
+          if (triangle->IsPresentTupel(TPS))
+            ++matching_triangles;
+        }
+        if (matching_triangles == 2)
+          break;
+      }
+      if (FindLine != FindPair.second) {
+        LOG(3, "INFO: Current line of point " << *point << " is " << *FindLine->second << ".");
+        all_lines_concave &= !FindLine->second->CheckConvexityCriterion();
+      }
+    }
+  } else {
+    // check the single line
+    if (connectedPath->empty())
+      return false;
+    LineMap::const_iterator FindLine = point->lines.find((*connectedPath->begin())->getNr());
+    ASSERT( FindLine != point->lines.end(),
+        "Tesselation::CheckAllConcaveInPolygon() - point "
+        +toString((*connectedPath->begin())->getNr())+" not present in "
+        +toString(*point)+"'s lines.");
+    return !FindLine->second->CheckConvexityCriterion();
+  }
+  return all_lines_concave;
+}
+/** Removes a boundary point from the envelope while keeping it closed.
+ * We remove the old triangles connected to the point and re-create new triangles to close the surface following this ansatz:
+ *  -# a closed path(s) of boundary points surrounding the point to be removed is constructed
+ *  -# on each closed path, we pick three adjacent points, create a triangle with them and subtract the middle point from the path
+ *  -# we advance two points (i.e. the next triangle will start at the ending point of the last triangle) and continue as before
+ *  -# the surface is closed, when the path is empty
+ * Thereby, we (hopefully) make sure that the removed points remains beneath the surface (this is checked via IsInnerPoint eventually).
+ * \param *out output stream for debugging
+ * \param *point point to be removed
+ * \return volume added to the volume inside the tesselated surface by the removal
+ */
+double Tesselation::RemoveFullConcavePointFromTesselatedSurface(class BoundaryPointSet *point)
+{
+  double volume = 0;
+  if (point == NULL) {
+    ELOG(1, "Cannot remove the point " << point << ", it's NULL!");
+    return 0.;
+  } else
+    LOG(4, "DEBUG: Removing point " << *point << " from tesselated boundary ...");
+  // get list of connected points
+  if (point->lines.empty()) {
+    ELOG(1, "Cannot remove the point " << *point << ", it's connected to no lines!");
+    return 0.;
+  }
+  list<TesselPointList *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
   list<TesselPointList *>::iterator ListAdvance = ListOfClosedPaths->begin();
   list<TesselPointList *>::iterator ListRunner = ListAdvance;
 //  TriangleMap::iterator NumberRunner = Candidates.begin();
   TesselPointList::iterator StartNode, MiddleNode, EndNode;
-  double angle;
-  double smallestangle;
   Vector Point, Reference, OrthogonalVector;
+  if (count > 2) { // less than three triangles, then nothing will be created
+    class TesselPoint *TriangleCandidates[3];
+    count = 0;
+    for (; ListRunner != ListOfClosedPaths->end(); ListRunner = ListAdvance) { // go through all closed paths
+      if (ListAdvance != ListOfClosedPaths->end())
+        ListAdvance++;
+      connectedPath = *ListRunner;
+      // re-create all triangles by going through connected points list
+      LineList NewLines;
+      for (; !connectedPath->empty();) {
+        // search middle node with widest angle to next neighbours
+        EndNode = connectedPath->end();
+        smallestangle = 0.;
+        for (MiddleNode = connectedPath->begin(); MiddleNode != connectedPath->end(); MiddleNode++) {
+          LOG(1, "INFO: MiddleNode is " << **MiddleNode << ".");
+          // construct vectors to next and previous neighbour
+          StartNode = MiddleNode;
+          if (StartNode == connectedPath->begin())
+            StartNode = connectedPath->end();
+          StartNode--;
+          //LOG(3, "INFO: StartNode is " << **StartNode << ".");
+          Point = ((*StartNode)->getPosition()) - ((*MiddleNode)->getPosition());
+          StartNode = MiddleNode;
+          StartNode++;
+          if (StartNode == connectedPath->end())
+            StartNode = connectedPath->begin();
+          //LOG(3, "INFO: EndNode is " << **StartNode << ".");
+          Reference = ((*StartNode)->getPosition()) - ((*MiddleNode)->getPosition());
+          OrthogonalVector = ((*MiddleNode)->getPosition()) - OldPoint;
+          OrthogonalVector.MakeNormalTo(Reference);
+          angle = GetAngle(Point, Reference, OrthogonalVector);
+          //if (angle < M_PI)  // no wrong-sided triangles, please?
+          if (fabs(angle - M_PI) < fabs(smallestangle - M_PI)) { // get straightest angle (i.e. construct those triangles with smallest area first)
+            smallestangle = angle;
+            EndNode = MiddleNode;
+          }
+        }
+        MiddleNode = EndNode;
+        if (MiddleNode == connectedPath->end()) {
+          ELOG(0, "CRITICAL: Could not find a smallest angle!");
+          performCriticalExit();
+        }
+        StartNode = MiddleNode;
+        if (StartNode == connectedPath->begin())
+          StartNode = connectedPath->end();
+        StartNode--;
+        EndNode++;
+        if (EndNode == connectedPath->end())
+          EndNode = connectedPath->begin();
+        LOG(2, "INFO: StartNode is " << **StartNode << ".");
+        LOG(2, "INFO: MiddleNode is " << **MiddleNode << ".");
+        LOG(2, "INFO: EndNode is " << **EndNode << ".");
+        LOG(1, "INFO: Attempting to create triangle " << (*StartNode)->getName() << ", " << (*MiddleNode)->getName() << " and " << (*EndNode)->getName() << ".");
+        TriangleCandidates[0] = *StartNode;
+        TriangleCandidates[1] = *MiddleNode;
+        TriangleCandidates[2] = *EndNode;
+        triangle = GetPresentTriangle(TriangleCandidates);
+        if (triangle != NULL) {
+          ELOG(0, "New triangle already present, skipping!");
+          StartNode++;
+          MiddleNode++;
+          EndNode++;
+          if (StartNode == connectedPath->end())
+            StartNode = connectedPath->begin();
+          if (MiddleNode == connectedPath->end())
+            MiddleNode = connectedPath->begin();
+          if (EndNode == connectedPath->end())
+            EndNode = connectedPath->begin();
+          continue;
+        }
+        LOG(3, "Adding new triangle points.");
+        AddTesselationPoint(*StartNode, 0);
+        AddTesselationPoint(*MiddleNode, 1);
+        AddTesselationPoint(*EndNode, 2);
+        LOG(3, "Adding new triangle lines.");
+        AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
+        AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
+        NewLines.push_back(BLS[1]);
+        AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
+        BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+        BTS->GetNormalVector(NormalVector);
+        AddTesselationTriangle();
+        // calculate volume summand as a general tetraeder
+        volume += CalculateVolumeofGeneralTetraeder(TPS[0]->node->getPosition(), TPS[1]->node->getPosition(), TPS[2]->node->getPosition(), OldPoint);
+        // advance number
+        count++;
+        // prepare nodes for next triangle
+        StartNode = EndNode;
+        LOG(2, "Removing " << **MiddleNode << " from closed path, remaining points: " << connectedPath->size() << ".");
+        connectedPath->remove(*MiddleNode); // remove the middle node (it is surrounded by triangles)
+        if (connectedPath->size() == 2) { // we are done
+          connectedPath->remove(*StartNode); // remove the start node
+          connectedPath->remove(*EndNode); // remove the end node
+          break;
+        } else if (connectedPath->size() < 2) { // something's gone wrong!
+          ELOG(0, "CRITICAL: There are only two endpoints left!");
+          performCriticalExit();
+        } else {
+          MiddleNode = StartNode;
+          MiddleNode++;
+          if (MiddleNode == connectedPath->end())
+            MiddleNode = connectedPath->begin();
+          EndNode = MiddleNode;
+          EndNode++;
+          if (EndNode == connectedPath->end())
+            EndNode = connectedPath->begin();
+        }
+      }
+      // maximize the inner lines (we preferentially created lines with a huge angle, which is for the tesselation not wanted though useful for the closing)
+      if (NewLines.size() > 1) {
+        LineList::iterator Candidate;
+        class BoundaryLineSet *OtherBase = NULL;
+        double tmp, maxgain;
+        do {
+          maxgain = 0;
+          for (LineList::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
+            tmp = PickFarthestofTwoBaselines(*Runner);
+            if (maxgain < tmp) {
+              maxgain = tmp;
+              Candidate = Runner;
+            }
+          }
+          if (maxgain != 0) {
+            volume += maxgain;
+            LOG(1, "Flipping baseline with highest volume" << **Candidate << ".");
+            OtherBase = FlipBaseline(*Candidate);
+            NewLines.erase(Candidate);
+            NewLines.push_back(OtherBase);
+          }
+        } while (maxgain != 0.);
+      }
+      ListOfClosedPaths->remove(connectedPath);
+      delete (connectedPath);
+    }
+    LOG(1, "INFO: " << count << " triangles were created.");
+  } else {
+    while (!ListOfClosedPaths->empty()) {
+      ListRunner = ListOfClosedPaths->begin();
+      connectedPath = *ListRunner;
+      ListOfClosedPaths->remove(connectedPath);
+      delete (connectedPath);
+    }
+    LOG(3, "DEBUG: No need to create any triangles.");
+  for (; ListRunner != ListOfClosedPaths->end(); ListRunner = ListAdvance) { // go through all closed paths
+    if (ListAdvance != ListOfClosedPaths->end())
+      ListAdvance++;
+    TesselPointList *connectedPath = *ListRunner;
+    if (CheckAllConcaveInPolygon(connectedPath, point)) {
+      LOG(1, "INFO: ... point " << *point << " cannot be on convex envelope, all lines concave.");
+      volume += RemovePointSurroundedByPolygon(connectedPath, point);
+    }
+    ListOfClosedPaths->remove(connectedPath);
+    delete (connectedPath);
+  }
   delete (ListOfClosedPaths);
+  LOG(1, "INFO: Removed volume is " << volume << ".");
+  if (volume > 0.)
+    LOG(1, "INFO: Removed volume is " << volume << ".");
   return volume;
 …
             if (TrianglePoints[j] != NULL) {
               for (FindLine = TrianglePoints[i]->lines.find(TrianglePoints[j]->node->getNr()); // is a multimap!
               (FindLine != TrianglePoints[i]->lines.end()) && (FindLine->first == TrianglePoints[j]->node->getNr()); FindLine++) {
+                  (FindLine != TrianglePoints[i]->lines.end()) && (FindLine->first == TrianglePoints[j]->node->getNr()); FindLine++) {
                 for (FindTriangle = FindLine->second->triangles.begin(); FindTriangle != FindLine->second->triangles.end(); FindTriangle++) {
                   if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
 …
           break;
       for (FindLine = TrianglePoints[(i + 1) % 3]->lines.find(TrianglePoints[(i + 2) % 3]->node->getNr()); // is a multimap!
       (FindLine != TrianglePoints[(i + 1) % 3]->lines.end()) && (FindLine->first == TrianglePoints[(i + 2) % 3]->node->getNr()); FindLine++) {
+          (FindLine != TrianglePoints[(i + 1) % 3]->lines.end()) && (FindLine->first == TrianglePoints[(i + 2) % 3]->node->getNr()); FindLine++) {
         for (FindTriangle = FindLine->second->triangles.begin(); FindTriangle != FindLine->second->triangles.end(); FindTriangle++) {
           if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
 …
+    }
     default:
+      ELOG(0, "Number of wildcards is greater than 3, cannot happen!");
+      performCriticalExit();
+      ASSERT(0, "Tesselation::FindTriangles() - Number of wildcards is greater than 3, cannot happen!");
       break;
+  }
 …
     if (a->endpoints[lowerNra] < b->endpoints[lowerNrb])
       return true;
+    else if (a->endpoints[lowerNra] > b->endpoints[lowerNrb])
+      return false;
+    else { // both lower-numbered endpoints are the same ...
+      if (a->endpoints[(lowerNra + 1) % 2] < b->endpoints[(lowerNrb + 1) % 2])
+        return true;
+      else if (a->endpoints[(lowerNra + 1) % 2] > b->endpoints[(lowerNrb + 1) % 2])
+    else
+      if (a->endpoints[lowerNra] > b->endpoints[lowerNrb])
         return false;
+    }
+      else { // both lower-numbered endpoints are the same ...
+        if (a->endpoints[(lowerNra + 1) % 2] < b->endpoints[(lowerNrb + 1) % 2])
+          return true;
+        else
+          if (a->endpoints[(lowerNra + 1) % 2] > b->endpoints[(lowerNrb + 1) % 2])
+            return false;
+      }
     return false;
+  }
 …
   // sanity check
   if (LinesOnBoundary.empty()) {
     ELOG(2, "FindAllDegeneratedTriangles() was called without any tesselation structure.");
     return DegeneratedLines;
+  }
   LineMap::iterator LineRunner1;
   pair<UniqueLines::iterator, bool> tester;
+      if (LinesOnBoundary.empty()) {
+        ELOG(2, "FindAllDegeneratedTriangles() was called without any tesselation structure.");
+        return DegeneratedLines;
+      }
+      LineMap::iterator LineRunner1;
+      pair<UniqueLines::iterator, bool> tester;
   for (LineRunner1 = LinesOnBoundary.begin(); LineRunner1 != LinesOnBoundary.end(); ++LineRunner1) {
     tester = AllLines.insert(LineRunner1->second);
 …
     const LineMap::const_iterator Line2 = LinesOnBoundary.find((*it).second);
     if (Line1 != LinesOnBoundary.end() && Line2 != LinesOnBoundary.end())
       LOG(3, "DEBUG: " << *Line1->second << " => " << *Line2->second);
+    LOG(3, "DEBUG: " << *Line1->second << " => " << *Line2->second);
     else
       ELOG(1, "Either " << (*it).first << " or " << (*it).second << " are not in LinesOnBoundary!");
+    ELOG(1, "Either " << (*it).first << " or " << (*it).second << " are not in LinesOnBoundary!");
+  }
 …
       for (TriangleRunner2 = line2->triangles.begin(); TriangleRunner2 != line2->triangles.end(); ++TriangleRunner2) {
         if ((TriangleRunner1->second != TriangleRunner2->second) && (TriangleRunner1->second->IsPresentTupel(TriangleRunner2->second))) {
           DegeneratedTriangles->insert(pair<int, int> (TriangleRunner1->second->Nr, TriangleRunner2->second->Nr));
           DegeneratedTriangles->insert(pair<int, int> (TriangleRunner2->second->Nr, TriangleRunner1->second->Nr));
+          DegeneratedTriangles->insert(pair<int, int>(TriangleRunner1->second->Nr, TriangleRunner2->second->Nr));
+          DegeneratedTriangles->insert(pair<int, int>(TriangleRunner2->second->Nr, TriangleRunner1->second->Nr));
+        }
+      }
 …
   // iterate over all degenerated triangles
+  for (IndexToIndex::iterator TriangleKeyRunner = DegeneratedTriangles->begin(); !DegeneratedTriangles->empty(); TriangleKeyRunner = DegeneratedTriangles->begin()) {
+  for (IndexToIndex::iterator TriangleKeyRunner = DegeneratedTriangles->begin();
+      !DegeneratedTriangles->empty();
+      TriangleKeyRunner = DegeneratedTriangles->begin()) {
     LOG(3, "DEBUG: Checking presence of triangles " << TriangleKeyRunner->first << " and " << TriangleKeyRunner->second << ".");
     // both ways are stored in the map, only use one
 …
 //                OtherpartnerTriangle = TriangleRunner->second;
 //              }
             /// interchanges their lines so that triangle->lines[i] == partnerTriangle->lines[j]
             // the line of triangle receives the degenerated ones
             triangle->lines[i]->triangles.erase(Othertriangle->Nr);
+              /// interchanges their lines so that triangle->lines[i] == partnerTriangle->lines[j]
+              // the line of triangle receives the degenerated ones
+              triangle->lines[i]->triangles.erase(Othertriangle->Nr);
             triangle->lines[i]->triangles.insert(TrianglePair(partnerTriangle->Nr, partnerTriangle));
             for (int k = 0; k < 3; k++)
 …
       // erase the pair
       count += (int) DegeneratedTriangles->erase(triangle->Nr);
+      count += (int)DegeneratedTriangles->erase(triangle->Nr);
       LOG(4, "DEBUG: RemoveDegeneratedTriangles() removes triangle " << *triangle << ".");
       RemoveTesselationTriangle(triangle);
       count += (int) DegeneratedTriangles->erase(partnerTriangle->Nr);
+      count += (int)DegeneratedTriangles->erase(partnerTriangle->Nr);
       LOG(4, "DEBUG: RemoveDegeneratedTriangles() removes triangle " << *partnerTriangle << ".");
       RemoveTesselationTriangle(partnerTriangle);
     } else {
       LOG(4, "DEBUG: RemoveDegeneratedTriangles() does not remove triangle " << *triangle << " and its partner " << *partnerTriangle << " because it is essential for at" << " least one of the endpoints to be kept in the tesselation structure.");
+      // we need to remove them from the list nonetheless
+      DegeneratedTriangles->erase(triangle->Nr);
+      DegeneratedTriangles->erase(partnerTriangle->Nr);
+    }
+  }
 …
   class BoundaryLineSet *BestLine = NULL;
   for (LineMap::iterator Runner = NearestBoundaryPoint->lines.begin(); Runner != NearestBoundaryPoint->lines.end(); Runner++) {
+    BaseLine = (Runner->second->endpoints[0]->node->getPosition()) -
+               (Runner->second->endpoints[1]->node->getPosition());
+    CenterToPoint = 0.5 * ((Runner->second->endpoints[0]->node->getPosition()) +
+                           (Runner->second->endpoints[1]->node->getPosition()));
+    BaseLine = (Runner->second->endpoints[0]->node->getPosition()) - (Runner->second->endpoints[1]->node->getPosition());
+    CenterToPoint = 0.5 * ((Runner->second->endpoints[0]->node->getPosition()) + (Runner->second->endpoints[1]->node->getPosition()));
     CenterToPoint -= (point->getPosition());
     angle = CenterToPoint.Angle(BaseLine);
     if (fabs(angle - M_PI/2.) < fabs(BestAngle - M_PI/2.)) {
+    if (fabs(angle - M_PI / 2.) < fabs(BestAngle - M_PI / 2.)) {
       BestAngle = angle;
       BestLine = Runner->second;
 …
   for (int i = 0; i < 3; i++) { // look for the same line as BestLine (only it's its degenerated companion)
     if ((BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[0])) && (BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[1]))) {
+      if (BestLine == BTS->lines[i]) {
+        ELOG(0, "BestLine is same as found line, something's wrong here!");
+        performCriticalExit();
+      }
+      BTS->lines[i]->triangles.insert(pair<int, class BoundaryTriangleSet *> (TempTriangle->Nr, TempTriangle));
+      ASSERT(BestLine != BTS->lines[i], std::string("Tesselation::AddBoundaryPointByDegeneratedTriangle() - ") + std::string("BestLine is same as found line, something's wrong here!"));
+      BTS->lines[i]->triangles.insert(pair<int, class BoundaryTriangleSet *>(TempTriangle->Nr, TempTriangle));
       TempTriangle->lines[nr] = BTS->lines[i];
       break;
 …
   if (LastTriangle != NULL) {
     stringstream sstr;
     sstr << "-"<< TrianglesOnBoundary.size() << "-" << LastTriangle->getEndpointName(0) << "_" << LastTriangle->getEndpointName(1) << "_" << LastTriangle->getEndpointName(2);
+    sstr << "-" << TrianglesOnBoundary.size() << "-" << LastTriangle->getEndpointName(0) << "_" << LastTriangle->getEndpointName(1) << "_" << LastTriangle->getEndpointName(2);
     NumberName = sstr.str();
     if (DoTecplotOutput) {
 …
     if (s1->endpoints.size() < s2->endpoints.size())
       return true;
+    else if (s1->endpoints.size() > s2->endpoints.size())
+      return false;
+    else { // equality of number of endpoints
+      PointSet::const_iterator Walker1 = s1->endpoints.begin();
+      PointSet::const_iterator Walker2 = s2->endpoints.begin();
+      while ((Walker1 != s1->endpoints.end()) || (Walker2 != s2->endpoints.end())) {
+        if ((*Walker1)->Nr < (*Walker2)->Nr)
+          return true;
+        else if ((*Walker1)->Nr > (*Walker2)->Nr)
+          return false;
+        Walker1++;
+        Walker2++;
+    else
+      if (s1->endpoints.size() > s2->endpoints.size())
+        return false;
+      else { // equality of number of endpoints
+        PointSet::const_iterator Walker1 = s1->endpoints.begin();
+        PointSet::const_iterator Walker2 = s2->endpoints.begin();
+        while ((Walker1 != s1->endpoints.end()) || (Walker2 != s2->endpoints.end())) {
+          if ((*Walker1)->Nr < (*Walker2)->Nr)
+            return true;
+          else
+            if ((*Walker1)->Nr > (*Walker2)->Nr)
+              return false;
+          Walker1++;
+          Walker2++;
+        }
+        return false;
+      }
-      return false;
+    }
+  }
 };
 …
       for (TriangleMap::const_iterator TriangleRunner = (LineRunner->second)->triangles.begin(); TriangleRunner != (LineRunner->second)->triangles.end(); TriangleRunner++) {
         if (DegeneratedTriangles->find(TriangleRunner->second->Nr) == DegeneratedTriangles->end()) {
           TriangleInsertionTester = TriangleVectors.insert(pair<int, Vector *> ((TriangleRunner->second)->Nr, &((TriangleRunner->second)->NormalVector)));
+          TriangleInsertionTester = TriangleVectors.insert(pair<int, Vector *>((TriangleRunner->second)->Nr, &((TriangleRunner->second)->NormalVector)));
           if (TriangleInsertionTester.second)
             LOG(5, "DEBUG:  Adding triangle " << *(TriangleRunner->second) << " to triangles to check-list.");
 …
         if (VectorWalker != VectorRunner) { // skip equals
           const double SCP = VectorWalker->second->ScalarProduct(*VectorRunner->second); // ScalarProduct should result in -1. for degenerated triangles
           LOG(4, "DEBUG: Checking " << *VectorWalker->second << " against " << *VectorRunner->second << ": " << SCP);
+          LOG(4, "DEBUG: Checking " << *(VectorWalker->second) << " against " << *(VectorRunner->second) << ": " << SCP);
           if (fabs(SCP + 1.) < ParallelEpsilon) {
             InsertionTester = EndpointCandidateList.insert((Runner->second));
 …
     // The Problem is probably due to two degenerated polygons being connected by a bridging, non-degenerated polygon, as somehow one node has
     // connections to either polygon ...
+    if (T->size() % 2 != 0) {
+      ELOG(0, " degenerated polygon contains an odd number of triangles, probably contains bridging non-degenerated ones, too!");
+      performCriticalExit();
+    }
+    ASSERT(T->size() % 2 == 0, std::string("Tesselation::CorrectAllDegeneratedPolygons() - ") + std::string(" degenerated polygon contains an odd number of triangles,") + std::string(" probably contains bridging non-degenerated ones, too!"));
     TriangleSet::iterator TriangleWalker = T->begin(); // is the inner iterator
     /// 4a. Get NormalVector for one side (this is "front")

src/Tesselation/tesselation.hpp

-              r3b1798
+              r0516fd
     void AddTesselationPoint(TesselPoint* Candidate, const int n);
     void SetTesselationPoint(TesselPoint* Candidate, const int n) const;
     void AddTesselationLine(const Vector * OptCenter, const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n);
+    bool AddTesselationLine(const Vector * OptCenter, const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n);
     void AddNewTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n);
     void AddExistingTesselationTriangleLine(class BoundaryLineSet *FindLine, int n);
 …
     void RemoveTesselationPoint(class BoundaryPointSet *point);
     bool CheckDegeneracy(CandidateForTesselation &CandidateLine, const double RADIUS, const LinkedCell_deprecated *LC) const;
+    bool isConvex() const;
     // concave envelope
 …
     void GuessStartingTriangle();
     bool InsertStraddlingPoints(IPointCloud & cloud, const LinkedCell_deprecated *LC);
+    double RemovePointSurroundedByPolygon(
+        TesselPointList *connectedPath,
+        BoundaryPointSet *point);
+    bool CheckAllConcaveInPolygon(
+        const TesselPointList *connectedPath,
+        const BoundaryPointSet *point
+        );
+    double RemoveFullConcavePointFromTesselatedSurface(class BoundaryPointSet *point);
     double RemovePointFromTesselatedSurface(class BoundaryPointSet *point);
     class BoundaryLineSet * FlipBaseline(class BoundaryLineSet *Base);
 …
     DistanceToPointMap * FindClosestBoundaryPointsToVector(const Vector &x, const LinkedCell_deprecated* LC) const;
     BoundaryLineSet * FindClosestBoundaryLineToVector(const Vector &x, const LinkedCell_deprecated* LC) const;
+    bool IsPointBelowSurroundingPolygon(const BoundaryPointSet *_point) const;
     // print for debugging

src/Tesselation/tesselationhelpers.cpp

-              r3b1798
+              r0516fd
   if (helper.ScalarProduct(SearchDirection) < -HULLEPSILON)  // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
     alpha = 2.*M_PI - alpha;
   LOG(3, "DEBUG: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << RelativeOldSphereCenter << " and resulting angle is " << alpha << ".");
+  LOG(5, "DEBUG: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << RelativeOldSphereCenter << " and resulting angle is " << alpha << ".");
   radius = helper.distance(RelativeOldSphereCenter);
   helper.ProjectOntoPlane(NormalVector);
   // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
   if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {
     LOG(4, "DEBUG: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << ".");
+    LOG(6, "DEBUG: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << ".");
     return alpha;
   } else {
     LOG(3, "DEBUG: NewSphereCenter " << RelativeNewSphereCenter << " is too close to RelativeOldSphereCenter" << RelativeOldSphereCenter << ".");
+    LOG(5, "DEBUG: NewSphereCenter " << RelativeNewSphereCenter << " is too close to RelativeOldSphereCenter" << RelativeOldSphereCenter << ".");
     return 2.*M_PI;
+  }
 …
+  }
   LOG(1, "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << ".");
+  LOG(4, "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << ".");
   return phi;
 …
  * \param *c third vector
  * \param *d fourth vector
  * \return \f$ \frac{1}{6} \cdot ((a-d) \times (a-c) \cdot  (a-b)) \f$
+ * \return \f$ \frac{1}{6} | (a-d) \cdot  \bigl ( (b-d) \times (c-d) \bigr ) | \f$
  */
 double CalculateVolumeofGeneralTetraeder(const Vector &a, const Vector &b, const Vector &c, const Vector &d)
 …
   for (int j=0;j<3;j++)
     TetraederVector[j].SubtractVector(d);
   Point = TetraederVector[0];
   Point.VectorProduct(TetraederVector[1]);
   volume = 1./6. * fabs(Point.ScalarProduct(TetraederVector[2]));
+  Point = TetraederVector[1];
+  Point.VectorProduct(TetraederVector[2]);
+  volume = 1./6. * fabs(Point.ScalarProduct(TetraederVector[0]));
   return volume;
 };
 …
   Normal.VectorProduct(OtherBaseline);
   Normal.Normalize();
   LOG(1, "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << ".");
+  LOG(3, "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << ".");
   // project one offset point of OtherBase onto this plane (and add plane offset vector)
 …
   *Intersection = line1.getIntersection(line2);
   Normal = (*Intersection) - (Base->endpoints[0]->node->getPosition());
   LOG(1, "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(Baseline)/Baseline.NormSquared()) << ".");
+  LOG(3, "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(Baseline)/Baseline.NormSquared()) << ".");
   return Intersection;
 …
     *tecplot << endl;
     // print connectivity
     LOG(1, "The following triangles were created:");
+    LOG(4, "DEBUG: The following triangles were created:");
     for (TriangleMap::const_iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) {
       LOG(1, " " << runner->second->endpoints[0]->node->getName() << "<->" << runner->second->endpoints[1]->node->getName() << "<->" << runner->second->endpoints[2]->node->getName());
+      LOG(4, " " << runner->second->endpoints[0]->node->getName() << "<->" << runner->second->endpoints[1]->node->getName() << "<->" << runner->second->endpoints[2]->node->getName());
       *tecplot << LookupList[runner->second->endpoints[0]->node->getNr()] << " " << LookupList[runner->second->endpoints[1]->node->getNr()] << " " << LookupList[runner->second->endpoints[2]->node->getNr()] << endl;
+    }
 …
   for (PointMap::const_iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
     point = PointRunner->second;
     LOG(1, "INFO: Current point is " << *point << ".");
+    LOG(2, "INFO: Current point is " << *point << ".");
     // calculate mean concavity over all connected line

src/Tesselation/triangleintersectionlist.cpp

-              r3b1798
+              r0516fd
     // if we have found one, check Scalarproduct between the vector
     Vector TestVector = (Point) - (*(*runner).second);
     LOG(1, "INFO: Distance vector between point " << Point
+    LOG(4, "DEBUG: Distance vector between point " << Point
         << " and triangle intersection at " << (*(*runner).second) << " is "
         << TestVector << ".");
     if (fabs(TestVector.NormSquared()) < MYEPSILON)  {//
       LOG(1, "ACCEPT: Point is on the intersected triangle.");
+      LOG(3, "ACCEPT: Point is on the intersected triangle.");
       return true;
+    }
     const double sign = (*runner).first->NormalVector.ScalarProduct(TestVector);
     LOG(1, "INFO: Checking sign of SKP between Distance vector and triangle's NormalVector "
+    LOG(4, "DEBUG: Checking sign of SKP between Distance vector and triangle's NormalVector "
         << (*runner).first->NormalVector << ": " << sign << ".");
     if (sign < 0) {
       LOG(1, "ACCEPT: Point lies on inner side of intersected triangle.");
+      LOG(3, "ACCEPT: Point lies on inner side of intersected triangle.");
       return true;
     } else {
       LOG(1, "REJECT: Point lies on outer side of intersected triangle.");
+      LOG(3, "REJECT: Point lies on outer side of intersected triangle.");
       return false;
+    }
 …
       for (DistanceTriangleMap::const_iterator iter = DistanceRange.first;
           iter != DistanceRange.second; ++iter) {
         LOG(3, "DEBUG: " << *(iter->second) << " is in the list with " << iter->first << ".");
+        LOG(4, "DEBUG: " << *(iter->second) << " is in the list with " << iter->first << ".");
         ++count;
+      }
       LOG(1, "INFO: There are " << count << " possible triangles at the smallest distance.");
+      LOG(3, "INFO: There are " << count << " possible triangles at the smallest distance.");
+    }
     // if there is more than one, check all of their normal vectors
 …
     // would be truely inside, all of the closest triangles would have to show
     // their inner side
     LOG(1, "INFO: Looking for first SKP greater than zero ...");
+    LOG(4, "DEBUG: Looking for first SKP greater than zero ...");
     for (DistanceTriangleMap::const_iterator iter = DistanceRange.first;
         iter != DistanceRange.second; ++iter) {
 …
       // construct SKP
       const double sign = (*iter).second->NormalVector.ScalarProduct(TestVector);
       LOG(1, "INFO: Checking SKP of closest triangle " << *(iter->second) << " = " << sign << ".");
+      LOG(4, "DEBUG: Checking SKP of closest triangle " << *(iter->second) << " = " << sign << ".");
       // if positive return
       if (sign > 0)

src/UIElements/Views/Qt4/Plotting/QSeisPlotPage.cpp

r3b1798	r0516fd
39	39	QSeisPlotPage::~QSeisPlotPage()
40	40	{
	41	if (layoutInitialized)
41	42	deinitializeLayout();
42	43	}

src/documentation/constructs/tesselation.dox

-              r3b1798
+              r0516fd
  * to such a surface.
+ *
  * \section tesselation-procedure
+ * \section tesselation-procedure Tesselation procedure
+ *
  * In the tesselation all atoms act as possible hindrance to a rolling sphere
 …
  * other (although the code for that is not yet fully implemented).
+ *
+ * \section tesselation-extension
+ * \section tesselation-convexization Making a surface convex
+ *
+ * A closed surface created by the aforementioned procedure can be made convex
+ * by a combination of the following two ways:
+ * -# Removing a point from the surface that is connected to other points only
+ *    by concave lines. This might also be imagined as removing bumps or
+ *    craters in the surface.
+ * -# flipping a base line or rather adding a general tetrahedron to remove a
+ *    concave line on the surface.
+ *
+ * With the first way one has to pay attention to possible degenerated lines
+ * and triangles. That's why prior to the this convexization procedure all
+ * possible degenerated triangles are removed. Furthermore, when looking at
+ * a removal candidate and its connected points, all these points are split
+ * up into so-called connected paths. The crater to be removed or filled-up
+ * has a low point -- the point to be removed -- and a rim, defined by all
+ * points connected by concave lines to the low point. However, when a point
+ * has degenerated lines attached to it (i.e. two lines with the same
+ * endpoints), it may have multiple rims (imagine two craters on either
+ * side of the surface and the volume being so small/slim that they reach
+ * through to the same low point). We have to discern between these multiple
+ * rims, therefore the connected points are placed into different closed
+ * rings, so-called polygons. The point is the removed and the polygon re-
+ * tesselated which essentially fills the crater.
+ *
+ * With the second way, we have to pay attention that the filled-in tetrahedron
+ * does not intersect already present triangles. The baseline defines two
+ * points and as the tesselated surface is closed, it must be connected to two
+ * triangles. These together define a set of four points that make up the
+ * tetrahedron. Naturally, two sides of the tetrahedron are always present
+ * already (and will become removed in place of the other two, effectively
+ * adding more volume to the tesselation).
+ * Now first, we only flip base lines that are concave. Second, none of the two
+ * other sides of the tetrahedron must be present. And lastly, we must check
+ * for all surrounding triangles that the new baseline (formed by point 3
+ * and 4) does not intersect these. Essentially, if we imagine a plane
+ * containing this new baseline, then each possibly intersecting triangle shows
+ * up as a brief line segment. We have to make sure that all of these segments
+ * remain below the new baseline in this plane. Also, things are complicated
+ * as the first and last line segment will always intersect with the baseline
+ * at the endpoint. There, we basically have to make sure that the line segment
+ * goes off in the right direction, namely outward.
+ *
+ * \section tesselation-volume Measuring the volume contained
+ *
+ * There is no straight-forward way to measure the volume contained in a
+ * non-convex tesselated surface. However, there is for a convex surface
+ * because convexity means that a line between any inner point and a point on
+ * the boundary will not intersect the surface anywhere else. Hence, we may
+ * use the center of gravity of all boundary points (by the same argument it
+ * must be an inner point) and calculate the volume of the general tetrahedron
+ * by looking at each of the tesselation's triangles in turn and summing up.
+ *
+ * We can calculate the volume of the original non-convex tesselation because
+ * the two ways mentioned above -- removing points and flipping baselines --
+ * both involve just addding general tetrahedron whose volume we may easily
+ * calculate. By bookkeeping of how much volume is added and calculating the
+ * total convex volume in the end, we also get the volume contained in the
+ * prior non-convex surface.
+ *
+ * \section tesselation-extension Issues whebn extended a tesselated surface
+ *
  * The main problem for extending the mesh to match with the normal sense is
 …
  * from a combination of spheres with van der Waals radii.
+ *
  * \date 2014-03-10
+ * \date 2014-10-09
+ *
  */

tests/Python/AllActions/options.dat

-              r3b1798
+              r0516fd
 change-element  "1"
 change-molname  "water"
+convex-envelope "50."
 convex-file     "convexfile"
 copy-molecule   "0"
 …
 DoCyclesFull    "0"
 DoLongrange     "0"
+DoOutputEveryStep       "0"
 DoPrintDebug    "0"
 DoRotate        "0"

tests/Tesselations/Makefile.am

-              r3b1798
+              r0516fd
         molecuilder.in \
         package.m4 \
+        $(srcdir)/1_2-dimethoxyethane \
+        $(srcdir)/1_2-dimethylbenzene \
+        $(srcdir)/2-methylcyclohexanone \
+        $(srcdir)/benzene \
+        $(srcdir)/cholesterol \
+        $(srcdir)/cycloheptane \
+        $(srcdir)/dimethyl_bromomalonate \
+        $(srcdir)/glucose \
+        $(srcdir)/heptan \
+        $(srcdir)/isoleucine \
+        $(srcdir)/neohexane \
+        $(srcdir)/N_N-dimethylacetamide \
+        $(srcdir)/proline \
+        $(srcdir)/putrescine \
+        $(srcdir)/tartaric_acid
+        $(srcdir)/Convex/1_2-dimethoxyethane \
+        $(srcdir)/Convex/1_2-dimethylbenzene \
+        $(srcdir)/Convex/2-methylcyclohexanone \
+        $(srcdir)/Convex/amylose \
+        $(srcdir)/Convex/benzene \
+        $(srcdir)/Convex/cholesterol \
+        $(srcdir)/Convex/cycloheptane \
+        $(srcdir)/Convex/dimethyl_bromomalonate \
+        $(srcdir)/Convex/glucose \
+        $(srcdir)/Convex/heptan \
+        $(srcdir)/Convex/isoleucine \
+        $(srcdir)/Convex/neohexane \
+        $(srcdir)/Convex/N_N-dimethylacetamide \
+        $(srcdir)/Convex/proline \
+        $(srcdir)/Convex/putrescine \
+        $(srcdir)/Convex/tartaric_acid \
+        $(srcdir)/NonConvex/1_2-dimethoxyethane \
+        $(srcdir)/NonConvex/1_2-dimethylbenzene \
+        $(srcdir)/NonConvex/2-methylcyclohexanone \
+        $(srcdir)/NonConvex/benzene \
+        $(srcdir)/NonConvex/cholesterol \
+        $(srcdir)/NonConvex/cycloheptane \
+        $(srcdir)/NonConvex/dimethyl_bromomalonate \
+        $(srcdir)/NonConvex/glucose \
+        $(srcdir)/NonConvex/heptan \
+        $(srcdir)/NonConvex/isoleucine \
+        $(srcdir)/NonConvex/neohexane \
+        $(srcdir)/NonConvex/N_N-dimethylacetamide \
+        $(srcdir)/NonConvex/proline \
+        $(srcdir)/NonConvex/putrescine \
+        $(srcdir)/NonConvex/tartaric_acid
 TESTSUITE = $(srcdir)/testsuite
 …
 TESTSCRIPTS = \
+        benzene/testsuite-benzene.at \
+_2-dimethoxyethane/testsuite-1_2-dimethoxyethane.at \
+_2-dimethylbenzene/testsuite-1_2-dimethylbenzene.at \
+-methylcyclohexanone/testsuite-2-methylcyclohexanone.at \
+        cholesterol/testsuite-cholesterol.at \
+        cycloheptane/testsuite-cycloheptane.at \
+        dimethyl_bromomalonate/testsuite-dimethyl_bromomalonate.at \
+        glucose/testsuite-glucose.at \
+        isoleucine/testsuite-isoleucine.at \
+        neohexane/testsuite-neohexane.at \
+        N_N-dimethylacetamide/testsuite-N_N-dimethylacetamide.at \
+        proline/testsuite-proline.at \
+        putrescine/testsuite-putrescine.at \
+        tartaric_acid/testsuite-tartaric_acid.at
+        Convex/1_2-dimethoxyethane/testsuite-1_2-dimethoxyethane.at \
+        Convex/1_2-dimethylbenzene/testsuite-1_2-dimethylbenzene.at \
+        Convex/2-methylcyclohexanone/testsuite-2-methylcyclohexanone.at \
+        Convex/benzene/testsuite-benzene.at \
+        Convex/amylose/testsuite-amylose.at \
+        Convex/cholesterol/testsuite-cholesterol.at \
+        Convex/cycloheptane/testsuite-cycloheptane.at \
+        Convex/dimethyl_bromomalonate/testsuite-dimethyl_bromomalonate.at \
+        Convex/glucose/testsuite-glucose.at \
+        Convex/isoleucine/testsuite-isoleucine.at \
+        Convex/neohexane/testsuite-neohexane.at \
+        Convex/N_N-dimethylacetamide/testsuite-N_N-dimethylacetamide.at \
+        Convex/proline/testsuite-proline.at \
+        Convex/putrescine/testsuite-putrescine.at \
+        Convex/tartaric_acid/testsuite-tartaric_acid.at \
+        NonConvex/1_2-dimethoxyethane/testsuite-1_2-dimethoxyethane.at \
+        NonConvex/1_2-dimethylbenzene/testsuite-1_2-dimethylbenzene.at \
+        NonConvex/2-methylcyclohexanone/testsuite-2-methylcyclohexanone.at \
+        NonConvex/benzene/testsuite-benzene.at \
+        NonConvex/cholesterol/testsuite-cholesterol.at \
+        NonConvex/cycloheptane/testsuite-cycloheptane.at \
+        NonConvex/dimethyl_bromomalonate/testsuite-dimethyl_bromomalonate.at \
+        NonConvex/glucose/testsuite-glucose.at \
+        NonConvex/isoleucine/testsuite-isoleucine.at \
+        NonConvex/neohexane/testsuite-neohexane.at \
+        NonConvex/N_N-dimethylacetamide/testsuite-N_N-dimethylacetamide.at \
+        NonConvex/proline/testsuite-proline.at \
+        NonConvex/putrescine/testsuite-putrescine.at \
+        NonConvex/tartaric_acid/testsuite-tartaric_acid.at
 max_jobs = 4

tests/Tesselations/testsuite.at

-              r3b1798
+              r0516fd
 m4_ifdef([AT_COLOR_TESTS], [AT_COLOR_TESTS])
+### NONCONVEX TESSELATION
 # tesselation of 1_2-dimethoxyethane
 m4_include(1_2-dimethoxyethane/testsuite-1_2-dimethoxyethane.at)
+m4_include(NonConvex/1_2-dimethoxyethane/testsuite-1_2-dimethoxyethane.at)
 # tesselation of 1_2-dimethylbenzene
 m4_include(1_2-dimethylbenzene/testsuite-1_2-dimethylbenzene.at)
+m4_include(NonConvex/1_2-dimethylbenzene/testsuite-1_2-dimethylbenzene.at)
 # tesselation of 2-methylcyclohexanone
 m4_include(2-methylcyclohexanone/testsuite-2-methylcyclohexanone.at)
+m4_include(NonConvex/2-methylcyclohexanone/testsuite-2-methylcyclohexanone.at)
 # tesselation of benzene
 m4_include(benzene/testsuite-benzene.at)
+m4_include(NonConvex/benzene/testsuite-benzene.at)
 # tesselation of cholesterol
 m4_include(cholesterol/testsuite-cholesterol.at)
+m4_include(NonConvex/cholesterol/testsuite-cholesterol.at)
 # tesselation of cycloheptane
 m4_include(cycloheptane/testsuite-cycloheptane.at)
+m4_include(NonConvex/cycloheptane/testsuite-cycloheptane.at)
 # tesselation of dimethyl_bromomalonate
 m4_include(dimethyl_bromomalonate/testsuite-dimethyl_bromomalonate.at)
+m4_include(NonConvex/dimethyl_bromomalonate/testsuite-dimethyl_bromomalonate.at)
 # tesselation of glucose
 m4_include(glucose/testsuite-glucose.at)
+m4_include(NonConvex/glucose/testsuite-glucose.at)
 # tesselation of heptan
 m4_include(heptan/testsuite-heptan.at)
+m4_include(NonConvex/heptan/testsuite-heptan.at)
 # tesselation of isoleucine
 m4_include(isoleucine/testsuite-isoleucine.at)
+m4_include(NonConvex/isoleucine/testsuite-isoleucine.at)
 # tesselation of neohexane
 m4_include(neohexane/testsuite-neohexane.at)
+m4_include(NonConvex/neohexane/testsuite-neohexane.at)
 # tesselation of N_N-dimethylacetamide
 m4_include(N_N-dimethylacetamide/testsuite-N_N-dimethylacetamide.at)
+m4_include(NonConvex/N_N-dimethylacetamide/testsuite-N_N-dimethylacetamide.at)
 # tesselation of proline
 m4_include(proline/testsuite-proline.at)
+m4_include(NonConvex/proline/testsuite-proline.at)
 # tesselation of putrescine
 m4_include(putrescine/testsuite-putrescine.at)
+m4_include(NonConvex/putrescine/testsuite-putrescine.at)
 # tesselation of tartaric_acid
+m4_include(tartaric_acid/testsuite-tartaric_acid.at)
+m4_include(NonConvex/tartaric_acid/testsuite-tartaric_acid.at)
+### NONCONVEX TESSELATION
+# tesselation of 1_2-dimethoxyethane
+m4_include(Convex/1_2-dimethoxyethane/testsuite-1_2-dimethoxyethane.at)
+# tesselation of 1_2-dimethylbenzene
+m4_include(Convex/1_2-dimethylbenzene/testsuite-1_2-dimethylbenzene.at)
+# tesselation of 2-methylcyclohexanone
+m4_include(Convex/2-methylcyclohexanone/testsuite-2-methylcyclohexanone.at)
+# tesselation of amylose
+m4_include(Convex/amylose/testsuite-amylose.at)
+# tesselation of benzene
+m4_include(Convex/benzene/testsuite-benzene.at)
+# tesselation of cholesterol
+m4_include(Convex/cholesterol/testsuite-cholesterol.at)
+# tesselation of cycloheptane
+m4_include(Convex/cycloheptane/testsuite-cycloheptane.at)
+# tesselation of dimethyl_bromomalonate
+m4_include(Convex/dimethyl_bromomalonate/testsuite-dimethyl_bromomalonate.at)
+# tesselation of glucose
+m4_include(Convex/glucose/testsuite-glucose.at)
+# tesselation of heptan
+m4_include(Convex/heptan/testsuite-heptan.at)
+# tesselation of isoleucine
+m4_include(Convex/isoleucine/testsuite-isoleucine.at)
+# tesselation of neohexane
+m4_include(Convex/neohexane/testsuite-neohexane.at)
+# tesselation of N_N-dimethylacetamide
+m4_include(Convex/N_N-dimethylacetamide/testsuite-N_N-dimethylacetamide.at)
+# tesselation of proline
+m4_include(Convex/proline/testsuite-proline.at)
+# tesselation of putrescine
+m4_include(Convex/putrescine/testsuite-putrescine.at)
+# tesselation of tartaric_acid
+m4_include(Convex/tartaric_acid/testsuite-tartaric_acid.at)

tests/regression/Tesselation/BigConvex/testsuite-tesselation-big-convex-envelope.at

-              r3b1798
+              r0516fd
 AT_SETUP([Tesselation - big convex Envelope])
 AT_KEYWORDS([Tesselation convex convex-envelope])
-AT_XFAIL_IF([/bin/true])
 file=test.conf
 AT_CHECK([/bin/cp -f ${abs_top_srcdir}/tests/regression/Tesselation/BigConvex/pre/test./conf $file], 0)
+AT_CHECK([/bin/cp -f ${abs_top_srcdir}/tests/regression/Tesselation/BigConvex/pre/test.conf $file], 0)
 convexfile=ConvexEnvelope
 AT_CHECK([../../molecuilder -i $file  --select-molecule-by-id 0 --convex-envelope --convex-file $convexfile --nonconvex-file NonConvexEnvelope], 0, [stdout], [stderr])
+AT_CHECK([../../molecuilder -i $file  --select-molecule-by-id 0 --convex-envelope 5. --convex-file $convexfile --nonconvex-file NonConvexEnvelope], 0, [stdout], [stderr])
 AT_CHECK([diff ${convexfile}.dat ${abs_top_srcdir}/tests/regression/Tesselation/BigConvex/post/ConvexEnvelope.dat], 0, [ignore], [ignore])
 #AT_CHECK([diff ${convexfile}.r3d ${abs_top_srcdir}/tests/regression/Tesselation/BigConvex/post/ConvexEnvelope.r3d], 0, [ignore], [ignore])
 AT_CHECK([fgrep "tesselated volume area is 16.4016 angstrom^3" stdout], 0, [ignore], [ignore])
+AT_CHECK([fgrep "tesselated volume area is 223.577 angstrom^3" stdout], 0, [ignore], [ignore])
 AT_CLEANUP

tests/regression/Tesselation/Convex/testsuite-tesselation-convex-envelope.at

-              r3b1798
+              r0516fd
 AT_CHECK([/bin/cp -f ${abs_top_srcdir}/tests/regression/Tesselation/Convex/pre/test.conf $file], 0)
 convexfile=ConvexEnvelope
 AT_CHECK([../../molecuilder -i $file  --select-molecule-by-id 0 --convex-envelope --convex-file $convexfile --nonconvex-file NonConvexEnvelope], 0, [stdout], [stderr])
+AT_CHECK([../../molecuilder -i $file  --select-molecule-by-id 0 --convex-envelope 50. --convex-file $convexfile --nonconvex-file NonConvexEnvelope], 0, [stdout], [stderr])
 AT_CHECK([diff ${convexfile}.dat ${abs_top_srcdir}/tests/regression/Tesselation/Convex/post/ConvexEnvelope.dat], 0, [ignore], [ignore])
 #AT_CHECK([diff ${convexfile}.r3d ${abs_top_srcdir}/tests/regression/Tesselation/Convex/post/ConvexEnvelope.r3d], 0, [ignore], [ignore])
 AT_CHECK([fgrep "tesselated volume area is 16.4016 angstrom^3" stdout], 0, [ignore], [ignore])
+AT_CHECK([fgrep "tesselated volume area is 8.89109 angstrom^3" stdout], 0, [ignore], [ignore])
 AT_CHECK([diff ConvexEnvelope.dat NonConvexEnvelope.dat], 0, [ignore], [ignore])

Context Navigation

Changes in / [3b1798:0516fd]

Legend:

LinearAlgebra/src/LinearAlgebra/Line.cpp

configure.ac

doc/userguide/userguide.xml

m4/boost.m4

src/Actions/ActionQueue.cpp

src/Actions/TesselationAction/ConvexEnvelopeAction.cpp

src/Actions/TesselationAction/ConvexEnvelopeAction.def

src/Potentials/Specifics/ConstantPotential.cpp

src/Potentials/Specifics/FourBodyPotential_Improper.cpp

src/Potentials/Specifics/FourBodyPotential_Torsion.cpp

src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp

src/Potentials/Specifics/PairPotential_Harmonic.cpp

src/Potentials/Specifics/PairPotential_LennardJones.cpp

src/Potentials/Specifics/PairPotential_Morse.cpp

src/Potentials/Specifics/ThreeBodyPotential_Angle.cpp

src/Tesselation/BoundaryLineSet.cpp

src/Tesselation/BoundaryPointSet.cpp

src/Tesselation/BoundaryTriangleSet.cpp

src/Tesselation/CandidateForTesselation.cpp

src/Tesselation/boundary.cpp

src/Tesselation/boundary.hpp

src/Tesselation/tesselation.cpp

src/Tesselation/tesselation.hpp

src/Tesselation/tesselationhelpers.cpp

src/Tesselation/triangleintersectionlist.cpp

src/UIElements/Views/Qt4/Plotting/QSeisPlotPage.cpp

src/documentation/constructs/tesselation.dox

tests/Python/AllActions/options.dat

tests/Tesselations/Makefile.am

tests/Tesselations/testsuite.at

tests/regression/Tesselation/BigConvex/testsuite-tesselation-big-convex-envelope.at

tests/regression/Tesselation/Convex/testsuite-tesselation-convex-envelope.at

Download in other formats: