Mercurial > repos > public > sbplib

changeset 99:ecf77a50d4fe feature/arclen-param

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
The arclength parameter function was slow because it called fzero. Now it is constructed once and for all with splines. Much better performance.
author Martin Almquist <martin.almquist@it.uu.se>
date Sun, 06 Dec 2015 15:41:42 +0100
parents b30f3d8845f4
children 2fd0fd3e75eb bc5db54f9efd
files +grid/Curve.m testArcLength.m
diffstat 2 files changed, 3 insertions(+), 57 deletions(-) [+]
line wrap: on
line diff
diff -r b30f3d8845f4 -r ecf77a50d4fe +grid/Curve.m
--- a/+grid/Curve.m	Sun Dec 06 13:47:12 2015 +0100
+++ b/+grid/Curve.m	Sun Dec 06 15:41:42 2015 +0100
@@ -108,8 +108,9 @@
             arcVec = grid.Curve.arcLength(obj.gp,0,tvec);
             arcLength = grid.Curve.spline(tvec,arcVec);
    
-            % Stretch the parameter
-            arcPar = @(t) util.fzero_vec(@(s)arcLength(s) - t*arcLength(1),[0-10*eps,1+10*eps]);
+            % Stretch the parameter, construct function with splines
+            arcParVec = util.fzero_vec(@(s)arcLength(s) - tvec*arcLength(1),[0-10*eps,1+10*eps]);
+            arcPar = grid.Curve.spline(tvec,arcParVec);
             
             % New function and derivative
             g_new = @(t)obj.g(arcPar(t));
diff -r b30f3d8845f4 -r ecf77a50d4fe testArcLength.m
--- a/testArcLength.m	Sun Dec 06 13:47:12 2015 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,55 +0,0 @@
-m = 201; L = 1; order = 4;
-close all;
-tic
-
-g1 = @(t) 2*pi*t.*sin(2*pi*t);
-g2 = @(t) 2*pi*t.*cos(2*pi*t);
-g = @(t) [g1(t);g2(t)];
-
-t = linspace(0,L,m)'; dt = L/(m-1);
-ops = sbp.Ordinary(m,dt,order);
-D = ops.derivatives.D1;
-
-C = grid.Curve(g,[],[],D);
-
-% Function
-figure;
-C.plot([],'bo');
-hold on
-plot(g1(t),g2(t),'r-');
-drawnow
-
-% Derivative
-figure
-C.plot_derivative([],'bo');
-hold on;
-plot(2*pi*sin(2*pi*t) + (2*pi)^2*t.*cos(2*pi*t),...
-    2*pi*cos(2*pi*t) - (2*pi)^2*t.*sin(2*pi*t),'r-')
-drawnow
-
-% Arc length
-L = C.arc_length_fun(t);
-figure;
-plot(t,L)
-drawnow
-
-% Stretch curve
-C2 = C.stretch_parameter();
-z = linspace(0,1,m);
-gnew = C2.g(z);
-gpnew = C2.gp(z);
-
-% Compare stretched and unstretched curves.
-figure
-plot(g1(t),g2(t),'b*',gnew(1,:),gnew(2,:),'ro');
-
-% Compare stretched and unstretched derivatives. 
-figure
-theta = linspace(0,2*pi,100);
-plot(cos(theta),sin(theta),'-b',gpnew(1,:),gpnew(2,:),'rx');
-
-toc
-
-
-
-