Merge pull request #323 from VibekeSkytt/lrref_affected

Modified search for affected B-splines in lr refinement

Author: sbriseid

lrsplines2D/include/GoTools/lrsplines2D/LRSplineUtils.h

@@ -135,14 +135,16 @@ namespace Go
 			 LRSplineSurface::BSplineMap& bmap,
 			 double domain[], 
 			 std::vector<std::unique_ptr<BSplineUniLR> >& bspline_vec1,
-			 std::vector<std::unique_ptr<BSplineUniLR> >& bspline_vec2);
+			 std::vector<std::unique_ptr<BSplineUniLR> >& bspline_vec2,
+			 bool support);
 
     std::tuple<int, int, int, int>
       refine_mesh(Direction2D d, double fixed_val, double start, double end, 
 		  int mult, bool absolute,
 		  int spline_degree, double knot_tol,
 		  Mesh2D& mesh, 
-		  std::vector<std::unique_ptr<BSplineUniLR> >& bsplines);
+		  std::vector<std::unique_ptr<BSplineUniLR> >& bsplines,
+		  bool& refined);
 
     bool support_equal(const LRBSpline2D* b1, const LRBSpline2D* b2);
 
@@ -184,6 +186,12 @@ namespace Go
 			   bool u_at_end, bool v_at_end, 
 			   std::vector<Point>& result);
 
+    void
+      get_affected_bsplines(const std::vector<LRSplineSurface::Refinement2D>& refs, 
+			    const LRSplineSurface::ElementMap& emap,
+			    double knot_tol, const Mesh2D& mesh,
+			    std::vector<LRBSpline2D*>& affected);
+    
     //==============================================================================
     struct support_compare
     //==============================================================================

lrsplines2D/include/GoTools/lrsplines2D/Mesh2D.h

@@ -57,8 +57,8 @@ namespace Go
 struct GPos { 
   // due to Visual Studio 2010 not supporting initializer lists, we have
   // to make an explicit constructor here.
-  GPos(int i, int m) : ix(i), mult(m) {}
-  GPos() : ix(-1), mult(-1) {}
+GPos(int i, int m) : ix(i), mult(m) {}
+GPos() : ix(-1), mult(-1) {}
 
   int ix; 
   int mult;
@@ -87,12 +87,11 @@ class Mesh2D : public MeshLR
   // values.
   template<typename Array>
   Mesh2D(const Array& xknots, const Array& yknots);
-  
 
   Mesh2D(const std::vector<double>& xknots, const std::vector<double>& yknots,
 	 const std::vector<std::vector<int> >& mrvecx,
 	 const std::vector<std::vector<int> >& mrvecy);
-   
+  
    // Read the mesh from a stream
   virtual void read(std::istream& is);        
 
@@ -240,7 +239,7 @@ class Mesh2D : public MeshLR
   // ix - the index of the row/column of the meshrectangles
   // start - the index to the start of the first consecutive meshrectangle along the line
   // end   - the index to the one-past-end of the last consecutive meshrectangle along the line
-  void setMult(Direction2D d, int ix, int start, int end, int mult);
+  bool setMult(Direction2D d, int ix, int start, int end, int mult);
 
   // increment multiplicity of a consecutive set of meshrectangles by 'mult'
   // The consecutive set is specified by:
@@ -398,6 +397,7 @@ std::vector<int> Mesh2D::compactify_ixvec_(Iterator kvec_start, Iterator kvec_en
   mult.clear();  
   for (auto i = result.begin(); i != result.end(); ++i) 
       mult.push_back(std::count(kvec_start, kvec_end, *i));
+
   return result;
 }
 
@@ -517,9 +517,11 @@ inline Direction2D flip(Direction2D d)
 } 
 
 // defining streaming operators
-inline std::ostream& operator<<(std::ostream& os, const GPos& g)  { return os << g.ix << " " << g.mult << " ";}
-inline std::istream& operator>>(std::istream& is, GPos& g)        { return is >> g.ix >> g.mult;}
-inline std::ostream& operator<<(std::ostream& os, const Mesh2D& m){ m.write(os); return os;}
+inline std::ostream& operator<<(std::ostream& os, const GPos& g)  
+{ return os << g.ix << " " << g.mult << " ";}
+inline std::istream& operator>>(std::istream& is, GPos& g)
+{ return is >> g.ix >> g.mult;}
+inline std::ostream& operator<<(std::ostream& os, const Mesh2D& m) { m.write(os); return os;}
 inline std::istream& operator>>(std::istream& is, Mesh2D& m)      { m.read(is); return is;}
 
 }; // end namespace Go

lrsplines2D/src/LRSplineSurface.C

@@ -44,6 +44,7 @@
 #include <iostream> // @@ debug
 #include <fstream>
 #include <iterator> // @@ debug - remove
+#include <string.h>
 //#include <chrono>   // @@ debug
 #include <set>
 #include <tuple>
@@ -716,10 +717,12 @@ void LRSplineSurface::refine(Direction2D d, double fixed_val, double start,
   // Make a copy of the initial mesh
   Mesh2D mesh2 = mesh_;
 
+  bool refined;
   const auto indices = // tuple<int, int, int, int>
-  LRSplineUtils::refine_mesh(d, fixed_val, start, end, mult, absolute, 
-			     degree(d), knot_tol_, mesh_,
-			     (d == XFIXED) ?  bsplinesuni1_ : bsplinesuni2_);
+    LRSplineUtils::refine_mesh(d, fixed_val, start, end, mult, 
+			       absolute, degree(d), knot_tol_, mesh_,
+			       (d == XFIXED) ?  bsplinesuni1_ : bsplinesuni2_,
+			       refined);
 
 #ifdef DEBUG
   std::ofstream of2("mesh1.eps");
@@ -788,6 +791,22 @@ void LRSplineSurface::refine(Direction2D d, double fixed_val, double start,
 #ifdef DEBUG
 	elems_affected.push_back(curr_el);
 #endif
+	// if (d == XFIXED)
+	//   {
+	//     double u1 = curr_el->umin();
+	//     double u2 = curr_el->umax();
+	//     if ((u2 - fixed_val > 0.55*(u2-u1) || fixed_val - u1 > 0.55*(u2-u1)) &&
+	// 	u2-fixed_val > 0.0001 && fixed_val-u1 > 0.0001)
+	//       std::cout << "Unbalanced element, 1. par: " << fixed_val << " in [" << u1 << "," << u2 << "]" << std::endl;
+	//   }
+	// else
+	//   {
+	//     double v1 = curr_el->vmin();
+	//     double v2 = curr_el->vmax();
+	//     if ((v2 - fixed_val > 0.55*(v2-v1) || fixed_val - v1 > 0.55*(v2-v1)) &&
+	// 	v2-fixed_val > 0.0001 && fixed_val-v1 > 0.0001)
+	//       std::cout << "Unbalanced element, 2. par: " << fixed_val << " in [" << v1 << "," << v2 << "]" << std::endl;
+	//   }
       }
   }
   vector<LRBSpline2D*> bsplines_affected(all_bsplines.begin(), all_bsplines.end());
@@ -891,7 +910,7 @@ void LRSplineSurface::refine(Direction2D d, double fixed_val, double start,
   // combination of objects and pointers will not work.
     LRSplineUtils::iteratively_split2(bsplines_affected, mesh_, 
 				      bsplines_, domain, 
-				      bsplinesuni1_, bsplinesuni2_);
+				      bsplinesuni1_, bsplinesuni2_, true);
 				      // bsplinesuni1_, iu1, iu2,
 				      // bsplinesuni2_, iv1, iv2);
 
@@ -999,8 +1018,7 @@ void LRSplineSurface::refine(Direction2D d, double fixed_val, double start,
 
 #ifdef DEBUG
 	    if (it2 == emap_.end())
-	      int stop_break = 1;
-	    //std::cout << "LRSplineSurface::refine : Element not found" << std::endl;
+	      std::cout << "LRSplineSurface::refine : Element not found" << std::endl;
 #endif
 	  }
 
@@ -1031,14 +1049,9 @@ void LRSplineSurface::refine(Direction2D d, double fixed_val, double start,
 	    it2->second->getOutsidePoints(data_points, d, sort_in_u);
 	    it2->second->getOutsideSignificantPoints(significant_points, 
 						     d, sort_in_u_significant);
- 	    it2->second->getOutsideGhostPoints(ghost_points, d, 
+	    it2->second->getOutsideGhostPoints(ghost_points, d, 
 					       sort_in_u_ghost);
 	    pt_del = it2->second->getNmbValPrPoint();
-	    // it2->second->getAccuracyInfo(averr, maxerr, nmbout);
-	    // accerr = it2->second->getAccumulatedError();
-
-	    // Update accuracy statistices in element
-	    it2->second->updateAccuracyInfo();
 
 	    // Update supported LRBsplines
 	    for (size_t kb=0; kb<bsplines_affected.size(); ++kb)
@@ -1079,17 +1092,14 @@ void LRSplineSurface::refine(Direction2D d, double fixed_val, double start,
 	    if (data_points.size() > 0)
 	      elem->addDataPoints(data_points.begin(), data_points.end(),
 				  sort_in_u, pt_del);
-           if (significant_points.size() > 0)
-             elem->addSignificantPoints(significant_points.begin(), 
-                                        significant_points.end(),
-                                        sort_in_u_significant, pt_del);
+	    if (significant_points.size() > 0)
+	      elem->addSignificantPoints(significant_points.begin(), 
+					 significant_points.end(),
+					 sort_in_u_significant, pt_del);
+	      
 	    if (ghost_points.size() > 0)
 	      elem->addGhostPoints(ghost_points.begin(), ghost_points.end(),
 				   sort_in_u_ghost, pt_del);
-	    //elem->setAccuracyInfo(accerr, averr, maxerr, nmbout);  // Not exact info as the
-	    // element has been split
-	    elem->updateAccuracyInfo();  // Accuracy statistic in element
-
 	    emap_.insert(std::make_pair(key, std::move(elem)));
 	    //auto it3 = emap_.find(key);
 
@@ -1107,6 +1117,22 @@ void LRSplineSurface::refine(Direction2D d, double fixed_val, double start,
 #endif
 }
 
+  struct refval
+  {
+    refval(int val, int m)
+    {
+      kval = val;
+      mult = m;
+    }
+    int kval, mult;
+  };
+
+  int compare_refs_sf(refval r1, refval r2)
+  {
+    return (r1.kval < r2.kval);
+  }
+
+
 //==============================================================================
   void LRSplineSurface::refine(const vector<Refinement2D>& refs, 
 			     bool absolute)
@@ -1123,7 +1149,14 @@ void LRSplineSurface::refine(Direction2D d, double fixed_val, double start,
     int stop_break = 1;
   }
 #endif
+  // Collect B-splines affected by the refinements
+  vector<LRBSpline2D*> affected;
+  LRSplineUtils::get_affected_bsplines(refs, emap_, knot_tol_, mesh_, 
+  				       affected);
+
+  vector<refval> u_refs, v_refs;
   for (size_t i = 0; i != refs.size(); ++i) {
+    bool refined;
     const Refinement2D& r = refs[i];
     const auto indices = // tuple<int, int, int, int>
       LRSplineUtils::refine_mesh(r.d, 
@@ -1135,51 +1168,109 @@ void LRSplineSurface::refine(Direction2D d, double fixed_val, double start,
 				 degree(r.d), 
 				 knot_tol_, 
 				 mesh_, 
-				 (r.d == XFIXED) ?  bsplinesuni1_ : bsplinesuni2_);
+				 (r.d == XFIXED) ?  bsplinesuni1_ : bsplinesuni2_,
+				 refined);
 
-    // Not efficient
-    int fixed_ix = get<1>(indices); // Index of fixed_val in global knot vector.
-    int last_ix = 
-      BSplineUniUtils::last_overlapping_bsplineuni(fixed_ix,
-						   (r.d == XFIXED) ? bsplinesuni1_ : bsplinesuni2_);
-    
+    if (!refined)
+      continue;
+    refval curr(r.kval, r.multiplicity);
+    size_t kr;
     if (r.d == XFIXED)
+      {
+	for (kr=0; kr<u_refs.size(); ++kr)
+	  if (u_refs[kr].kval == curr.kval && 
+	      u_refs[kr].mult == curr.mult)
+	    break;
+	if (kr == u_refs.size())
+	  u_refs.push_back(curr);
+      }
+    else if (r.d == YFIXED)
+      {
+	for (kr=0; kr<v_refs.size(); ++kr)
+	  if (v_refs[kr].kval == curr.kval && 
+	      v_refs[kr].mult == curr.mult)
+	    break;
+	if (kr == v_refs.size())
+	  v_refs.push_back(curr);
+      }
+  }
+
+  if (u_refs.size() + v_refs.size() == 0)
+    return; // Mesh rectangles already existing
+  
+  //std::cout << "Post refine mesh" << std::endl;
+  std::sort(u_refs.begin(), u_refs.end(), compare_refs_sf);
+  std::sort(v_refs.begin(), v_refs.end(), compare_refs_sf);
+
+  for (int ki=(int)u_refs.size()-1; ki>=0; --ki)
+    {
+      int kval = u_refs[ki].kval;
+      int fixed_ix = Mesh2DUtils::last_nonlarger_knotvalue_ix(mesh_, XFIXED, kval);
+
+      int last_ix = 
+	BSplineUniUtils::last_overlapping_bsplineuni(fixed_ix, bsplinesuni1_);
+    
       LRSplineUtils::split_univariate(bsplinesuni1_, last_ix, fixed_ix, 
-				      (absolute) ? r.multiplicity : 1);
-    else
+				      absolute? u_refs[ki].mult : 1);
+    }
+
+   for (int ki=(int)v_refs.size()-1; ki>=0; --ki)
+    {
+      int kval = v_refs[ki].kval;
+      int fixed_ix = Mesh2DUtils::last_nonlarger_knotvalue_ix(mesh_, YFIXED, kval);
+
+      int last_ix = 
+	BSplineUniUtils::last_overlapping_bsplineuni(fixed_ix, bsplinesuni2_);
+    
       LRSplineUtils::split_univariate(bsplinesuni2_, last_ix, fixed_ix, 
-				      (absolute) ? r.multiplicity : 1);
-  }
+				      absolute? v_refs[ki].mult : 1); 
+    }
+ 
+  //   // Not efficient
+  //   int fixed_ix = get<1>(indices); // Index of fixed_val in global knot vector.
+  //   int last_ix = 
+  //     BSplineUniUtils::last_overlapping_bsplineuni(fixed_ix,
+  // 						   (r.d == XFIXED) ? bsplinesuni1_ : bsplinesuni2_);
+    
+  //   if (r.d == XFIXED)
+  //     LRSplineUtils::split_univariate(bsplinesuni1_, last_ix, fixed_ix, 
+  // 				      (absolute) ? r.multiplicity : 1);
+  //   else
+  //     LRSplineUtils::split_univariate(bsplinesuni2_, last_ix, fixed_ix, 
+  // 				      (absolute) ? r.multiplicity : 1);
+  // }
 
 
   //std::wcout << "Preparing for iterative splitting." << std::endl;
-  vector<unique_ptr<LRBSpline2D> > affected;
-  affected.reserve(bsplines_.size());
-//  for_each(bsplines_.begin(), bsplines_.end(), [&](const BSplineMap::value_type& b) {
-  for (auto it = bsplines_.begin(); it!= bsplines_.end(); ++it)
-    {
-      // @@@ VSK. This is maybe the place to remove element information from the bsplines?
-      unique_ptr<LRBSpline2D> ptr = std::move(it->second);
-      affected.emplace_back(std::move(ptr));//b.second);
-    };
+  //vector<unique_ptr<LRBSpline2D> > affected;
+//   vector<LRBSpline2D*> affected;
+//   affected.reserve(bsplines_.size());
+// //  for_each(bsplines_.begin(), bsplines_.end(), [&](const BSplineMap::value_type& b) {
+//   for (auto it = bsplines_.begin(); it!= bsplines_.end(); ++it)
+//     {
+//       affected.push_back((*it).second.get());
+//       // @@@ VSK. This is maybe the place to remove element information from the bsplines?
+//       // unique_ptr<LRBSpline2D> ptr = std::move(it->second);
+//       // affected.emplace_back(std::move(ptr));//b.second);
+//     };
 
   // @@@ VSK. In this case, we should not bother about splitting elements. They will
   // be regenerated later. Thus, the bsplines should NOT be updated with elements during
   // splitting
   // The bsplines should not have any pointers to elements. They will be set later
   //std::wcout << "Iteratively splitting." << std::endl;
 
-  LRSplineUtils::iteratively_split(affected, mesh_, 
-				   bsplinesuni1_, bsplinesuni2_);
-  bsplines_.clear();
+  LRSplineUtils::iteratively_split2(affected, mesh_, bsplines_, NULL, 
+				    bsplinesuni1_, bsplinesuni2_, false);
+  //bsplines_.clear();
 
   //std::wcout << "Splitting finished, now inserting resulting functions" << std::endl;
   // The bsplines are checked for duplicates and inserted in the global bspline map
 //  for_each(affected.begin(), affected.end(), [&](unique_ptr<LRBSpline2D> b) {
-for (auto it = affected.begin(); it != affected.end(); ++it)
-  {
-    LRSplineUtils::insert_basis_function(*it, mesh_, bsplines_);
-  };
+// for (auto it = affected.begin(); it != affected.end(); ++it)
+//   {
+//     LRSplineUtils::insert_basis_function(*it, mesh_, bsplines_);
+//   };
 
 #if 0//ndef NDEBUG
   {
@@ -1218,7 +1309,6 @@ for (auto it = affected.begin(); it != affected.end(); ++it)
 
 }
 
-
 //==============================================================================
   void LRSplineSurface::addSurface(const LRSplineSurface& other_sf, double fac)
 //==============================================================================
@@ -3573,6 +3663,7 @@ LRSplineSurface::collect_basis(int from_u, int to_u,
   return b_splines;
 }
 
+
 //===========================================================================
 void
 LRSplineSurface::checkSupport(LRBSpline2D* basis) const