Update with accumulated changes

compositemodel/include/GoTools/compositemodel/CurveModel.h

@@ -53,6 +53,7 @@ namespace Go
 {
 
  class CompositeCurve;
+  class Vertex;
 
 //===========================================================================
 /** A curve model including topological information
@@ -107,6 +108,12 @@ class CurveModel : public CompositeModel
   /// \return Index to curve
   int getIndex(ParamCurve* curve) const;
 
+    /// Fetch all edges
+    std::vector<shared_ptr<ftEdge> > allEdges()
+    {
+      return edges_;
+    }
+
   /// Evaluate position
   /// \param idx Index of curve
   /// \param par[] Parameter value
@@ -225,6 +232,11 @@ class CurveModel : public CompositeModel
   /// \return Vector of pointers to the composite curves
   std::vector<shared_ptr<CompositeCurve> > fetchCompositeCurves() const;
 
+  /// Return all vertices associated with this surface model
+  /// \retval vertices Vector of pointers to all vertices.
+  void getVertices(std::vector<shared_ptr<Vertex> >& vertices) const;
+
+
 private:
   std::vector<shared_ptr<ftEdge> > edges_;
 

compositemodel/include/GoTools/compositemodel/SurfaceModel.h

@@ -285,6 +285,21 @@ class GO_API SurfaceModel : public CompositeModel
   			std::vector<Point>& der) const;  // Result
 
 
+  /// Closest point between a given point and the outer boundary/boundaries
+  /// of this surface model
+  /// Returns one point
+  /// \param pnt Input point
+  /// \param clo_pnt Found closest point
+  /// \param idx Index of surface where the closest point is found
+  /// \param clo_par[] Parameter value corresponding to the closest point
+  /// \param dist Distance between input point and found closest point
+    void
+      closestBoundaryPoint(Point& pnt,     // Input point
+			   Point& clo_pnt, // Found closest point
+			   int& idx,          // Index of surface where the closest point is found
+			   double clo_par[],  // Parameter value corresponding to the closest point
+			   double& dist);     // Distance between input point and found closest point
+
   /// Closest point between a given point and this surface model
   /// Returns one point
   /// \param pnt Input point
@@ -304,7 +319,11 @@ class GO_API SurfaceModel : public CompositeModel
   /// \return Closest point
   ftPoint closestPoint(const Point& point);
 
-  /// Closest point between a given point and this surface model
+ void closestPoint(Point& point, int seed_ix, double seed[],
+		   Point& clo_pt, int& idx, double clo_par[],
+		   double& dist);
+
+ /// Closest point between a given point and this surface model
   /// \param point Input point
   /// \return Closest point
   ftPoint closestPoint(const ftPoint& point) { return closestPoint(point.position()); }
@@ -913,7 +932,7 @@ class GO_API SurfaceModel : public CompositeModel
 		      std::vector<std::pair<double,double> >& crv_bound,
 		      bool compute_curves=true) const;
 
-  ftPoint closestPointLocal(const ftPoint& point) const;
+  ftPoint closestPointLocal(const ftPoint& point, bool use_seed=false) const;
 
   void localExtreme(ftSurface *face, Point& dir, 
 		    Point& ext_pnt, int& ext_id,

compositemodel/include/GoTools/compositemodel/ftPointSet.h

@@ -118,6 +118,15 @@ class ftSamplePoint
     /// Am I on sub surface boundary?
     bool isOnSubSurfaceBoundary() const
     { return ((at_boundary_ == 2) || (at_boundary_ == 1)); }
+    /// Am I on an inner boundary?
+  bool isOnInnerBoundary() const
+  { return at_boundary_ == 2; }
+
+  virtual bool isCorner() const
+  {
+    return false;
+  }
+
     /// Get number of neighbours.
     int getNmbNeighbour() const
     { return (int)next_.size();}

compositemodel/src/AdaptSurface.C

@@ -1010,7 +1010,7 @@ namespace Go
     // Fetch constant parameter curves in the selected curve direction
     double u1 = (cv_dir == 0) ? dom.umin() : dom.vmin();
     double u2 = (cv_dir == 0) ? dom.umax() : dom.vmax();
-    double udel = (u2 - u1)/(double)(nmb_cv+1);
+    double udel = (u2 - u1)/(double)(nmb_cv+2);
     double tol1 = std::max(1.0e-7, 1.0e-5*(u2-u1));
     double par[2];
     //int ki, kj;
@@ -1071,7 +1071,7 @@ namespace Go
 	    nmb = std::max(nmb, min_samples);
 	    double v1 = crvs[kr]->startparam();
 	    double v2 = crvs[kr]->endparam();
-	    double vdel = (v2 - v1)/(double)(nmb+1);
+	    double vdel = (v2 - v1)/(double)(nmb+2);
 	    double tol2 = std::max(1.0e-7, 1.0e-5*(v2-v1));
 	    if (consider_joint)
 	      par[pt_dir] = std::min(v1+vdel, surf->nextSegmentVal(pt_dir, v1, true, tol2));

compositemodel/src/CurveModel.C

@@ -47,8 +47,8 @@
 
 using std::vector;
 
-namespace Go
-{
+using namespace Go;
+
   //===========================================================================
   CurveModel::CurveModel(double gap,   // Gap between adjacent curves
 			 double neighbour,  // Threshold for whether curves are adjacent
@@ -391,20 +391,40 @@ void CurveModel::tesselate(vector<shared_ptr<GeneralMesh> >& meshes) const
     vector<shared_ptr<Vertex> > vertices(all_vertices.begin(), 
 					 all_vertices.end());
     for (ki=0; ki<vertices.size(); ++ki)
-      for (kj=ki+1; kj<vertices.size(); ++kj)
-	{
-	  double dist = 
-	    vertices[ki]->getVertexPoint().dist(vertices[kj]->getVertexPoint());
-	  if (dist < toptol_.neighbour)
-	    {
-	      vector<ftEdge*> connected_edges = vertices[kj]->allEdges();
-	      vertices[ki]->joinVertex(vertices[kj]);
-	      for (size_t kr=0; kr<connected_edges.size(); ++kr)
-		connected_edges[kr]->replaceVertex(vertices[kj], vertices[ki]);
-	      vertices.erase(vertices.begin() + kj);
-	      break;
-	    }
+      {
+	for (kj=ki+1; kj<vertices.size(); )
+	  {
+	    double dist = 
+	      vertices[ki]->getVertexPoint().dist(vertices[kj]->getVertexPoint());
+	    if (dist < toptol_.neighbour)
+	      {
+		vector<ftEdge*> connected_edges = vertices[kj]->allEdges();
+		vertices[ki]->joinVertex(vertices[kj]);
+		for (size_t kr=0; kr<connected_edges.size(); ++kr)
+		  connected_edges[kr]->replaceVertex(vertices[kj], vertices[ki]);
+		vertices.erase(vertices.begin() + kj);
+	      }
+	    else
+	      ++kj;
+	  }
 	}
 
   }
-} // namespace Go
+
+  //===========================================================================
+  void CurveModel::getVertices(vector<shared_ptr<Vertex> >& vertices) const
+  //===========================================================================
+  {
+    std::set<shared_ptr<Vertex> > all_vertices;  // All vertices in the model represented once
+    for (size_t ki=0; ki<edges_.size(); ++ki)
+      {
+	shared_ptr<Vertex> vert = edges_[ki]->getVertex(true);
+	all_vertices.insert(vert);
+	vert = edges_[ki]->getVertex(false);
+	all_vertices.insert(vert);
+      }
+    vertices.clear();
+    vertices.insert(vertices.end(), all_vertices.begin(), all_vertices.end());
+  }
+
+

compositemodel/src/SurfaceModel.C

@@ -746,6 +746,53 @@ namespace Go
 #endif
   }
 
+  //===========================================================================
+  void SurfaceModel::closestBoundaryPoint(Point& pnt,     // Input point
+					  Point& clo_pnt, // Found closest point
+					  int& idx,           // Index of surface where the closest point is found
+					  double clo_par[],   // Parameter value corresponding to the closest point
+					  double& dist)       // Distance between input point and found closest point
+  //===========================================================================
+  {
+    idx = -1;
+    int loop_ix1 = -1, loop_ix2 = -1;
+    size_t cv_ix;
+    dist = std::numeric_limits<double>::max();
+    double cv_par;
+    for (size_t ki=0; ki<boundary_curves_.size(); ++ki)
+      {
+	for (size_t kj=0; kj<boundary_curves_[ki].size(); ++kj)
+	  {
+	    int cv_ind;
+	    double tpar, tdist;
+	    Point cv_close;
+	    boundary_curves_[ki][kj]->closestPoint(pnt, cv_ind, tpar,
+						   cv_close, tdist);
+	    if (tdist < dist)
+	      {
+		cv_ix = cv_ind;
+		loop_ix1 = (int)ki;
+		loop_ix2 = (int)kj;
+		cv_par = tpar;
+		dist = tdist;
+		clo_pnt = cv_close;
+	      }
+	  }
+      }
+
+    if (loop_ix1 >= 0)
+      {
+	shared_ptr<ftEdgeBase> edg0 =
+	  boundary_curves_[loop_ix1][loop_ix2]->getEdge(cv_ix);
+	ftEdge* edg = edg0->geomEdge();
+	ftFaceBase *face0 = edg->face();
+	ftSurface *face = face0->asFtSurface();
+	idx = getIndex(face);
+	Point fpar = edg->faceParameter(cv_par);
+	clo_par[0] = fpar[0];
+	clo_par[1] = fpar[1];
+      }
+  }
 
 
 
@@ -870,7 +917,34 @@ namespace Go
     return ftPoint(bestcp, bestface, bestu, bestv);
   }
 
-
+  //===========================================================================
+  void SurfaceModel::closestPoint(Point& point, int seed_ix, double seed[],
+				  Point& clo_pt, int& idx, double clo_par[],
+				  double& dist)
+  //===========================================================================
+  {
+    fill(face_checked_.begin(), face_checked_.end(), false);
+   ftPoint inpt;
+    bool use_seed = (seed_ix >= 0);
+    int ix = (seed_ix >= 0) ? seed_ix : 0;
+    ftSurface *first = static_pointer_cast<ftSurface>(faces_[ix]).get();
+    if (seed_ix >= 0)
+      inpt = ftPoint(point, first, seed[0], seed[1]);
+    else
+      inpt = ftPoint(point, first);
+    ftPoint close = closestPointLocal(inpt, use_seed);
+    if (close.face() != 0)
+      {
+	clo_pt = close.position();
+	idx = getIndex(close.face());
+	clo_par[0] = close.u();
+	clo_par[1] = close.v();
+	dist = point.dist(clo_pt);
+      }
+    else
+      closestPoint(point, clo_pt, idx, clo_par, dist);
+  }
+
 
   //===========================================================================
   int SurfaceModel::nmbEntities() const
@@ -1397,9 +1471,10 @@ void SurfaceModel::swapFaces(int idx1, int idx2)
   }
 
   //===========================================================================
-  ftPoint SurfaceModel::closestPointLocal(const ftPoint& point) const
+  ftPoint SurfaceModel::closestPointLocal(const ftPoint& point, bool use_seed) const
   //===========================================================================
   {
+    double tol = toptol_.gap;  // Should be a parameter tolerance
     const Point& pt = point.position();
     int id = getIndex(point.face());
     ftSurface* curface = 0;
@@ -1413,13 +1488,33 @@ void SurfaceModel::swapFaces(int idx1, int idx2)
     double closestpt_epsilon = toptol_.neighbour; // Maybe gap instead?
     ftSurface* bestface = 0;
     int nmb_checked = 0;
+    double seedval[2];
+    double *seed = use_seed ? seedval : 0;
+    seedval[0] = point.u();
+    seedval[1] = point.v();
     while (!finished) {
       if (!face_checked_[id]) {
 	curface = dynamic_cast<ftSurface*>(faces_[id].get());
 	face_checked_[id] = true;
 	//  	    cout << "Face: " << id << endl;
 	ASSERT(curface != 0);
-	curface->closestPoint(pt, u, v, cp, dist, closestpt_epsilon);
+	const Domain& dom = curface->surface()->parameterDomain();
+	if (seed)
+	  {
+	    bool in_domain = false;
+	    try {
+	      Vector2D seed2(seed[0], seed[1]);
+	      in_domain = dom.isInDomain(seed2, tol);
+             }
+	    catch (...)
+	      {
+		in_domain = false;
+	      }
+	    if (!in_domain)
+	      seed = 0;
+	  }
+	curface->closestPoint(pt, u, v, cp, dist, closestpt_epsilon, NULL,
+			      seed);
 	nmb_checked++;
 	if (dist < bestdist) {
 	  bestdist = dist;
@@ -1429,7 +1524,6 @@ void SurfaceModel::swapFaces(int idx1, int idx2)
 	  bestface = curface;
 	}
 	// Check if the point was on the boundary
-	const Domain& dom = curface->surface()->parameterDomain();
 	bool on_boundary = dom.isOnBoundary(Vector2D(u, v),
 					    toptol_.neighbour);
 	if (on_boundary) {
@@ -1438,11 +1532,59 @@ void SurfaceModel::swapFaces(int idx1, int idx2)
 	    = curface->edgeClosestToPoint(u, v);
 	  ftEdgeBase* twin = boundary_edge->twin();
 	  if (!twin) // That is, there is no neighbour
-	    finished = true;
-	  else {
+	    {
+	      bool corner = dom.isOnCorner(Vector2D(u, v),
+					   toptol_.neighbour);
+	      if (corner)
+		{
+		  shared_ptr<Vertex> vx1, vx2;
+		  boundary_edge->geomEdge()->getVertices(vx1, vx2);
+		  Point pt1 = vx1->getVertexPoint();
+		  Point pt2 = vx2->getVertexPoint();
+		  shared_ptr<Vertex> vx = (pt1.dist(pt) <= pt2.dist(pt)) ? vx1 : vx2;
+		  vector<ftEdge*> edgs = vx->getEdges(curface);
+		  for (size_t kr=0; kr<edgs.size(); ++kr)
+		    {
+		      if (edgs[kr] != boundary_edge)
+			{
+			  twin = edgs[kr]->twin();
+			  if (twin)
+			    break;
+			}
+		    }
+		}
+	      if (!twin)
+		finished = true;
+	    }
+	  if (twin)
+	    {
 	    //  		    cout << "We're crossing a boundary!" << endl;
-	    id = getIndex(twin->face()->asFtSurface());
-	  }
+	    int id2 = getIndex(twin->face()->asFtSurface());
+	    if (id2 == id)
+	      {
+		shared_ptr<ParamSurface> surf = getSurface(id);
+		SurfaceTools::surface_seedfind(pt, *surf, 0,
+					       seedval[0], seedval[1]);
+	      }
+	    else
+	      {
+		id = id2;
+		if (use_seed)
+		  {
+		    ftEdge *twin2 = twin->geomEdge();
+		    if (twin2)
+		      {
+			double tpar, tdist;
+			Point close;
+			twin2->closestPoint(bestcp, tpar, close, tdist);
+			Point fpar = twin2->faceParameter(tpar);
+			seedval[0] = fpar[0];
+			seedval[1] = fpar[1];
+			seed = seedval;
+		      }
+		  }
+	      }
+	    }
 	} else // point was in the interior
 	  finished = true;
       } else { // if face_checked_[id]
@@ -2149,7 +2291,8 @@ void SurfaceModel::swapFaces(int idx1, int idx2)
       }
   }
 
-//===========================================================================
+
+///===========================================================================
 vector<shared_ptr<ftEdge> > SurfaceModel::getBoundaryEdges() const
 //===========================================================================
 {