local refinement for 2d curves

libsrc/geom2d/genmesh2d.cpp

@@ -9,12 +9,14 @@ namespace netgen
 
 
 
-  void CalcPartition (double l, double h, double h1, double h2,
+  void CalcPartition (const SplineSegExt & spline, double l, 
-		      double hcurve, double elto0, Array<double> & points);
+		      MeshingParameters & mp, double h, double h1, double h2,
+		      double hcurve, double elto0, 
+		      Array<double> & points);
 
   // partitionizes spline curve
   void Partition (const SplineSegExt & spline,
-		  double h, double elto0,
+		  MeshingParameters & mp, double h, double elto0,
 		  Mesh & mesh, Point3dTree & searchtree, int segnr) 
   {
     enum { D = 2 };
@@ -31,7 +33,7 @@ namespace netgen
     double h2 = min (spline.EndPI().hmax, h/spline.EndPI().refatpoint);
     double hcurve = min (spline.hmax, h/spline.reffak);
 
-    CalcPartition (l, h, h1, h2, hcurve, elto0, curvepoints);
+    CalcPartition (spline, l, mp, h, h1, h2, hcurve, elto0, curvepoints);
     //  cout << "curvepoints = " << curvepoints << endl;
 
     dt = 1.0 / n;
@@ -130,7 +132,7 @@ namespace netgen
 
 
 
-  void SplineGeometry2d :: PartitionBoundary (double h, Mesh & mesh2d)
+  void SplineGeometry2d :: PartitionBoundary (MeshingParameters & mp, double h, Mesh & mesh2d)
   {
     enum { D = 2 };
     Box<D> bbox;
@@ -166,13 +168,13 @@ namespace netgen
 	  maximum = min2 ( maximum, h);
 	  if ( minimum > 0 )
 	    // GetSpline(i).Partition(minimum, elto0, mesh2d, searchtree, i+1);
-	    Partition(GetSpline(i), minimum, elto0, mesh2d, searchtree, i+1);	    
+	    Partition(GetSpline(i), mp, minimum, elto0, mesh2d, searchtree, i+1);	    
 	  else if ( maximum > 0 )
 	    // GetSpline(i).Partition(maximum, elto0, mesh2d, searchtree, i+1);
-	    Partition(GetSpline(i), maximum, elto0, mesh2d, searchtree, i+1);
+	    Partition(GetSpline(i), mp, maximum, elto0, mesh2d, searchtree, i+1);
 	  else
 	    // GetSpline(i).Partition(h, elto0, mesh2d, searchtree, i+1);
-	    Partition(GetSpline(i), h, elto0, mesh2d, searchtree, i+1);
+	    Partition(GetSpline(i), mp, h, elto0, mesh2d, searchtree, i+1);
 	}
       else
 	{
@@ -296,7 +298,7 @@ namespace netgen
     mesh = new Mesh;
     mesh->SetDimension (2);
 
-    geometry.PartitionBoundary (h, *mesh);
+    geometry.PartitionBoundary (mp, h, *mesh);
 
 
     // marks mesh points for hp-refinement

libsrc/geom2d/geometry2d.cpp

@@ -827,47 +827,52 @@ namespace netgen
 
 
 
-
+  /*
-
+  void CalcPartition (const SplineSegExt & spline,
-
+		      double l, double h, double h1, double h2,
-  void CalcPartition (double l, double h, double h1, double h2,
 		      double hcurve, double elto0, Array<double> & points)
   {
-    // cout << "calcpart, h = " << h << ", h1 = " << h1 << ", h2 = " << h2 << ", hcurve = " << hcurve << endl;
+    double fperel, oldf, f;
-    int i, j, n, nel;
-    double sum, t, dt, fun, fperel, oldf, f;
 
-    n = 1000;
+    int n = 1000;
 
     points.SetSize (0);
 
-    sum = 0;
+    double dt = l / n;
-    dt = l / n;
+
-    t = 0.5 * dt;
+    double sum = 0;
-    for (i = 1; i <= n; i++)
+    for (int i = 1; i <= n; i++)
       {
-	fun = min3 (hcurve, t/elto0 + h1, (l-t)/elto0 + h2);
+	double t = (i-0.5)*dt;
+	double fun = min3 (hcurve, t/elto0 + h1, (l-t)/elto0 + h2);
+	double curv = spline.CalcCurvature (t/l);
+	cout << "curv = " << curv << endl;
+	if (curv < 1e-10) curv = 1e-10;
+	fun = min2 (fun, 0.1/curv);
 	sum += dt / fun;
-	t += dt;
       }
 
-    nel = int (sum+1);
+    int nel = int (sum+1);
     fperel = sum / nel;
 
     points.Append (0);
 
-    i = 1;
+    int i = 1;
     oldf = 0;
-    t = 0.5 * dt;
+    // t = 0.5 * dt;
-    for (j = 1; j <= n && i < nel; j++)
+    for (int j = 1; j <= n && i < nel; j++)
       {
-	fun = min3 (hcurve, t/elto0 + h1, (l-t)/elto0 + h2);
+	double t = (j-0.5)*dt;
+	double fun = min3 (hcurve, t/elto0 + h1, (l-t)/elto0 + h2);
+	double curv = spline.CalcCurvature (t/l);
+	if (curv < 1e-10) curv = 1e-10;
+	fun = min2 (fun, 0.1/curv);
 
 	f = oldf + dt / fun;
 
-	while (f > i * fperel && i < nel)
+	while (i * fperel < f && i < nel)
 	  {
-	    points.Append ( (l/n) * (j-1 +  (i * fperel - oldf) / (f - oldf)) );
+	    points.Append ( dt * (j-1) +  (i * fperel - oldf) * fun);
 	    i++;
 	  }
 	oldf = f;
@@ -875,6 +880,104 @@ namespace netgen
       }
     points.Append (l);
   }
+  */
+
+
+
+ void CalcPartition (const SplineSegExt & spline, double l, 
+		     MeshingParameters & mp, double h, double h1, double h2,
+		     double hcurve, double elto0, Array<double> & points)
+  {
+    double fperel, oldf, f;
+
+    int n = 10000;
+    
+    Array<Point<2> > xi(n);
+    Array<double> hi(n);
+    
+    hi = h; 
+    for (int i = 0; i < n; i++)
+      {
+	xi[i] = spline.GetPoint ( (i+0.5) / n );
+	double hc = 1.0/mp.curvaturesafety / (1e-20+spline.CalcCurvature ( (i+0.5) / n ));
+	// *testout << "t = " << (i+0.5)/n << ", xi = " << xi[i] << ", hi = " << hc << endl;
+	hi[i] = min2(hi[i], hc);
+      }
+    hi[0] = min (hi[0], h1);
+    hi[n-1] = min (hi[n-1], h2); 
+    
+    if ( (spline.GetPoint(0)-spline.GetPoint(1)).Length() < 1e-12 )
+      hi[0] = hi[n-1] = min2 (hi[0], hi[n-1]);
+
+    // limit slope
+    double gradh = 1/elto0;
+    for (int i = 0; i < n-1; i++)
+      {
+	double hnext = hi[i] + gradh * (xi[i+1]-xi[i]).Length();
+	hi[i+1] = min(hi[i+1], hnext);
+      } 
+    for (int i = n-1; i > 1; i--)
+      {
+	double hnext = hi[i] + gradh * (xi[i-1]-xi[i]).Length();
+	hi[i-1] = min(hi[i-1], hnext);
+      }
+
+
+    points.SetSize (0);
+
+    double dt = l / n;
+
+    double sum = 0;
+    for (int i = 1; i <= n; i++)
+      {
+	double t = (i-0.5)*dt;
+	/*
+	double fun = min3 (hcurve, t/elto0 + h1, (l-t)/elto0 + h2);
+	double curv = spline.CalcCurvature (t/l);
+	cout << "curv = " << curv << endl;
+	if (curv < 1e-10) curv = 1e-10;
+	fun = min2 (fun, 0.1/curv);
+	*/
+	double fun = hi[i-1];
+	sum += dt / fun;
+      }
+
+    int nel = int (sum+1);
+    fperel = sum / nel;
+
+    points.Append (0);
+
+    int i = 1;
+    oldf = 0;
+    // t = 0.5 * dt;
+    for (int j = 1; j <= n && i < nel; j++)
+      {
+	double t = (j-0.5)*dt;
+	double fun = hi[j-1];
+	/*
+	double fun = min3 (hcurve, t/elto0 + h1, (l-t)/elto0 + h2);
+	double curv = spline.CalcCurvature (t/l);
+	if (curv < 1e-10) curv = 1e-10;
+	fun = min2 (fun, 0.1/curv);
+	*/
+	f = oldf + dt / fun;
+
+	while (i * fperel < f && i < nel)
+	  {
+	    points.Append ( dt * (j-1) +  (i * fperel - oldf) * fun);
+	    i++;
+	  }
+	oldf = f;
+	t += dt;
+      }
+    points.Append (l);
+  }
+
+
+
+
+
+
 
 
 

libsrc/geom2d/geometry2d.hpp

@@ -97,6 +97,15 @@ namespace netgen
       return seg.GetType();
     }
 
+    virtual double CalcCurvature (double t) const
+    {
+      Point<2> point;
+      Vec<2> first, second;
+      GetDerivatives (t, point, first, second);
+      double curv = fabs(first(0)*second(1)-first(1)*second(0)) / pow(first.Length(), 3);
+      return curv;
+    }
+    
   };
 
 
@@ -141,7 +150,7 @@ namespace netgen
     DLL_HEADER virtual int GenerateMesh (Mesh*& mesh, MeshingParameters & mparam,
 			      int perfstepsstart, int perfstepsend);
     
-    void PartitionBoundary (double h, Mesh & mesh2d);
+    void PartitionBoundary (MeshingParameters & mp, double h, Mesh & mesh2d);
 
     void CopyEdgeMesh (int from, int to, Mesh & mesh2d, Point3dTree & searchtree);
 

libsrc/gprim/spline.hpp

@@ -54,11 +54,23 @@ namespace netgen
     virtual Point<D> GetPoint (double t) const = 0;
     /// returns a (not necessarily unit-length) tangent vector for 0 <= t <= 1
     virtual Vec<D> GetTangent (const double t) const
-    { cerr << "GetTangent not implemented for spline base-class"  << endl; Vec<D> dummy; return dummy;}
+    { 
+      cerr << "GetTangent not implemented for spline base-class"  << endl; 
+      Vec<D> dummy; return dummy;
+    }
+
     virtual void GetDerivatives (const double t, 
 				 Point<D> & point,
 				 Vec<D> & first,
-				 Vec<D> & second) const {;}
+				 Vec<D> & second) const 
+    {
+      double eps = 1e-6;
+      point = GetPoint (t);
+      Point<D> pl = GetPoint (t-eps);
+      Point<D> pr = GetPoint (t+eps);
+      first = 1.0/(2*eps) * (pr-pl);
+      second = 1.0/sqr(eps) * ( (pr-point)+(pl-point));
+    }
 
 
     /// returns initial point on curve