From f0430e75a1207cb6c163f1fe569f9de3511f708b Mon Sep 17 00:00:00 2001 From: ahbarnett Date: Thu, 11 Apr 2024 22:26:28 -0400 Subject: [PATCH] add traction extraction demo to stfmm3dPerfTest.m --- matlab/Contents.m | 2 +- matlab/stfmm3dPerfTest.m | 43 ++++++++++++++++++++++++++++++++++++---- 2 files changed, 40 insertions(+), 5 deletions(-) diff --git a/matlab/Contents.m b/matlab/Contents.m index febb3f9..5d39faa 100644 --- a/matlab/Contents.m +++ b/matlab/Contents.m @@ -16,7 +16,7 @@ % hfmm3dTest % emfmm3dTest % stfmm3dTest -% stfmm3dPerfTest +% stfmm3dPerfTest - demo and math test for stokeslet & stresslet, traction % % For examples of use see: % lfmm3d_example diff --git a/matlab/stfmm3dPerfTest.m b/matlab/stfmm3dPerfTest.m index d33008c..d14bcb8 100644 --- a/matlab/stfmm3dPerfTest.m +++ b/matlab/stfmm3dPerfTest.m @@ -1,11 +1,13 @@ disp("Perf and math test matlab driver for FMM3D/Stokes") -disp("Source to target, stokeslet+stresslet (1-vec), vel pot only..."); -% Barnett 8/4/23 +% Barnett 8/4/23. Added traction output test, 4/12/24. clear ns = 100000; nt = 100000; eps = 1e-3; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +disp("Source to target, stokeslet+stresslet (1-vec), vel pot only..."); ifstrs = 1; % include stresslets? ifppreg = 0; % no eval at sources ifppregtarg = 1; % just vel out @@ -24,7 +26,7 @@ t = toc; fprintf("%d to %d points, eps=%.3g, done in %.3g s (%.3g tot pts/sec)\n",ns,nt,eps,t,(ns+nt)/t) u = U.pottarg; -i = randi(nt); % which targ +i = randi(nt); % which targ to check fprintf("velpot targ(i=%d) =\n",i); disp(u(:,i)) @@ -42,6 +44,39 @@ end end ui = 0.5 * ui; % FMM3D is 1/4pi off from true 1/8pi prefactor -fprintf("rel err vs direct at ith targ: %.3g\n",norm(u(:,i)-ui)/norm(ui)) +fprintf("rel err vs direct at ith targ: %.3g\n\n",norm(u(:,i)-ui)/norm(ui)) + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +disp("Stokeslet sources to targets: compute vel pot, traction at oriented targs..."); +% just give the changes from the above setup... +ifppregtarg = 3; % request everything +srcinfo = rmfield(srcinfo,{'strslet','strsvec'}); +targnor = rand(3,nt); % normal vectors for targets (eg, surface normals) + +tic; +U = stfmm3d(eps,srcinfo,ifppreg,targ,ifppregtarg); +t = toc; +fprintf("%d to %d points with u,p,gradu: done in %.3g s (%.3g tot pts/sec)\n",ns,nt,t,(ns+nt)/t) +p = U.pretarg; % pressure (1*ntarg) +gradu = U.gradtarg; % grad vel (3*3*ntarg) +shearstress = gradu + permute(gradu,[2 1 3]); % gradu + gradu^T +% stress sigma = -pI + mu*(gradu+gradu^T), then T = sigma.n We assume mu=1. +T = -(ones(3,1)*p).*targnor + squeeze(sum(shearstress .* reshape(kron(ones(3,1),targnor),3,3,nt),1)); % compute all tractions (painful tensor contraction in matlab) +i = randi(nt); % which targ to check +%fprintf("targ(i=%d): p=%.15g, gradu = \n",i,p(i)); disp(gradu(:,:,i)) +fprintf("targ(i=%d) traction T = \n",i); +disp(T(:,i)) +% check T (target tractions) vs direct formula +Ti = zeros(3,1); +nori = targnor(:,i); +for j=1:ns + R = targ(:,i) - srcinfo.sources(:,j); % x-y with std sign (targ-src). + r = sqrt(sum(R.^2)); % dist + f = srcinfo.stoklet(:,j); % strength + Ti = Ti - (3/(4*pi))*(1/r^5)*R*dot(nori,R)*dot(f,R); % true T_{ijk} n_j f_k +end +Ti = Ti * (4*pi); % FMM3D has 1/4pi missing in prefactor +fprintf("rel err vs direct at ith targ: %.3g\n",norm(T(:,i)-Ti)/norm(Ti))