Commit 1c713b07 authored by ppanchal's avatar ppanchal
Browse files

Sauter schwab based computation

parent 892b22ff
close all; clc; clear;
addpath(genpath("../../../"));
format long;
global X;
global W;
load('X3','X');
load('W3','W');
rng(10);
A = rand(3,3);
[Q,R] = qr(A);
% Initializing parameters for the problem
Nvals = 100:50:1700;
sz = size(Nvals,2);
l2errs = zeros(sz,1);
averrs = zeros(sz,1);
linferrs = zeros(sz,1);
Tn_nearest = zeros(sz,1);
Tn_plane = zeros(sz,1);
% Torus radii
r1 = 10;
r2 = 3;
N = 60;
for ii = 1:sz
N = Nvals(ii)
% distance between the centers
dist = 40;
% Mesh for the geometry
mesh_in = mshTorus(N,r1,r2);
mesh_out = mesh_in;
% Rotate the outer torus
mesh_out.vtx = mesh_out.vtx * Q;
% Translate the outer torus
N_vtcs = size(mesh_out.vtx,1);
trans = ones(N_vtcs,1) * [dist 0 0];
mesh_out.vtx = mesh_out.vtx + trans;
% Join to create the final mesh
mesh = union(mesh_in,mesh_out);
figure;
plot(mesh);
title('Mesh and normals');
hold on;
% Checking the normal direction for the mesh
ctrs = mesh.ctr;
nrms = mesh.nrm;
quiver3(ctrs(:,1),ctrs(:,2),ctrs(:,3),nrms(:,1),nrms(:,2),nrms(:,3));
scatter3(r1+r2,0,0);
%%
% Definition of FEM spaces and integration rule
S1_Gamma = fem(mesh,'P1');
S0_Gamma = fem(mesh,'P0');
S0_Gamma_in = fem(mesh_in,'P0');
order = 3;
Gamma = dom(mesh,order);
Op_in = restriction(S0_Gamma,mesh_in);
Op_out = restriction(S0_Gamma,mesh_out);
% Solving a Direct first kind BVP to get the representation of the state
% V Psi = (0.5*M+K) g_N (Interior problem)
% V Psi = (-0.5*M+K) g_N (Exterior problem)
% Getting the Single Layer matrix V using Gypsilab implementation
Gxy = @(X,Y)femGreenKernel(X,Y,'[1/r]',0); % 0 wave number
V = 1/(4*pi)*integral(Gamma,Gamma,S0_Gamma,Gxy,S0_Gamma);
V = V + 1/(4*pi)*regularize(Gamma,Gamma,S0_Gamma,'[1/r]',S0_Gamma);
% Getting the Double Layer matrix K using Gypsilab implementation
GradG = cell(3,1);
GradG{1} = @(X,Y)femGreenKernel(X,Y,'grady[1/r]1',0);
GradG{2} = @(X,Y)femGreenKernel(X,Y,'grady[1/r]2',0);
GradG{3} = @(X,Y)femGreenKernel(X,Y,'grady[1/r]3',0);
K = 1/(4*pi)*integral(Gamma,Gamma,S0_Gamma,GradG,ntimes(S0_Gamma));
K = K +1/(4*pi)*regularize(Gamma,Gamma,S0_Gamma,'grady[1/r]',ntimes(S0_Gamma));
% Defining the mass matrix M
M = integral(Gamma,S0_Gamma,S0_Gamma);
% Defining the Dirichlet boundary condition
% Cutoff radius for the BC
%R = 1;
R = dist/2;
%assert(R < Rad + s);
Vin = 10;
Vout = 20;
% Defining the Dirichlet boundary condition
g = @(x) 1/4/pi./vecnorm(x,2,2); % Point charge inside first torus
% Checking the boundary condition
figure;
plot(mesh);
hold on;
plot(mesh,g(S0_Gamma.dof));
title('Dirichlet Boundary Condition');
colorbar;
% Constructing the RHS
g_N = integral(Gamma,S0_Gamma,g);
% Exterior problem
Psi = V\((-0.5 * M + K)* (M\g_N));
Psi_in = Op_in * Psi;
Psi_out = Op_out * Psi;
% Solving the adjoint problem to get the adjoint solution
Rho = V\(-0.5 * g_N);
% Visualizing the Neumann trace
dofs = S0_Gamma.dof;
normals = mesh.nrm;
sPsi = -2*Psi;
figure;
plot(mesh);
hold on;
quiver3(dofs(:,1),dofs(:,2),dofs(:,3),normals(:,1).*sPsi,normals(:,2).*sPsi,normals(:,3).*sPsi,0);
title('Field on the surface');
%quiver3(dofs(:,1),dofs(:,2),dofs(:,3),normals(:,1),normals(:,2),normals(:,3));
% Visualizing the charge density on the surface of the objects
figure;
plot(mesh);
hold on;
plot(mesh,Psi);
title("Surface charge density");
colorbar;
%% Checking the convergence of the Neumann trace at a point
testpt = [r1+r2 0 0];
dofs = S0_Gamma.dof;
Ndofs = size(dofs,1);
distances = vecnorm((dofs - ones(Ndofs,1) * testpt),2,2);
[min_d,min_ind] = min(distances);
Tn_nearest(ii) = Psi(min_ind)
min_d
dofs(min_ind,:)
%% Using the plane method
elts = mesh.elt;
vtcs = mesh.vtx;
Nelts = size(elts,1);
dist_elts = zeros(Nelts,1);
for i = 1:Nelts
dist_elts(i) = dist_plane(testpt,vtcs(elts(i,1),:),vtcs(elts(i,2),:),vtcs(elts(i,3),:));
end
[mind_plane,minind_plane] = min(dist_elts)
Tn_plane(ii) = Psi(minind_plane)
%% Exact solution
err_vec = Psi;
% L2 Error
l2errs(ii) = err_vec' * M * err_vec %norm(err_vec)
averrs(ii) = err_vec' * V * err_vec
linferrs(ii) = max(abs(err_vec))
close all;
end
save('dirichlet_sol_tor_tor.mat',"averrs","l2errs","Tn_plane","Tn_nearest","Nvals","linferrs");
\ No newline at end of file
...@@ -8,7 +8,7 @@ load('W3','W'); ...@@ -8,7 +8,7 @@ load('W3','W');
rng(10); rng(10);
A = rand(3,3); A = rand(3,3);
[Q,R] = qr(A); [Qrot,R] = qr(A);
% Initializing parameters for the problem % Initializing parameters for the problem
...@@ -36,25 +36,26 @@ dist = 40; ...@@ -36,25 +36,26 @@ dist = 40;
mesh_in = mshTorus(N,r1,r2); mesh_in = mshTorus(N,r1,r2);
mesh_out = mesh_in; mesh_out = mesh_in;
% Rotate the outer torus
mesh_out.vtx = mesh_out.vtx * Q;
% Translate the outer torus % Translate the outer torus
N_vtcs = size(mesh_out.vtx,1); N_vtcs = size(mesh_out.vtx,1);
trans = ones(N_vtcs,1) * [dist 0 0]; trans = ones(N_vtcs,1) * [dist 0 0];
mesh_out.vtx = mesh_out.vtx + trans; mesh_out.vtx = mesh_out.vtx + trans;
% Rotate the inner torus
mesh_in.vtx = mesh_in.vtx * Qrot;
% Join to create the final mesh % Join to create the final mesh
mesh = union(mesh_in,mesh_out); mesh = union(mesh_in,mesh_out);
figure; figure;
plot(mesh); plot(mesh);
title('Mesh and normals'); title('Mesh');
hold on; hold on;
% Checking the normal direction for the mesh % Checking the normal direction for the mesh
ctrs = mesh.ctr; ctrs = mesh.ctr;
nrms = mesh.nrm; nrms = mesh.nrm;
quiver3(ctrs(:,1),ctrs(:,2),ctrs(:,3),nrms(:,1),nrms(:,2),nrms(:,3)); quiver3(ctrs(:,1),ctrs(:,2),ctrs(:,3),nrms(:,1),nrms(:,2),nrms(:,3));
scatter3(r1+r2,0,0);
%% %%
...@@ -93,8 +94,8 @@ M = integral(Gamma,S0_Gamma,S0_Gamma); ...@@ -93,8 +94,8 @@ M = integral(Gamma,S0_Gamma,S0_Gamma);
%R = 1; %R = 1;
R = dist/2; R = dist/2;
%assert(R < Rad + s); %assert(R < Rad + s);
Vin = 10; Vin = -0.007487063249394;
Vout = 20; Vout = 0;
g = @(X) (sqrt(sum(X.^2,2)) > R)* Vout + (sqrt(sum(X.^2,2)) <= R) * Vin; g = @(X) (sqrt(sum(X.^2,2)) > R)* Vout + (sqrt(sum(X.^2,2)) <= R) * Vin;
% Checking the boundary condition % Checking the boundary condition
...@@ -144,6 +145,7 @@ distances = vecnorm((dofs - ones(Ndofs,1) * testpt),2,2); ...@@ -144,6 +145,7 @@ distances = vecnorm((dofs - ones(Ndofs,1) * testpt),2,2);
Tn_nearest(ii) = Psi(min_ind) Tn_nearest(ii) = Psi(min_ind)
min_d min_d
dofs(min_ind,:) dofs(min_ind,:)
charge_in = sum( mesh_in.ndv .* Psi_in, 1)
%% Using the plane method %% Using the plane method
elts = mesh.elt; elts = mesh.elt;
......
loglog(Nvals,abs(torque_mst(:,1)))
hold on
loglog(Nvals,abs(torque_mst(:,2)))
loglog(Nvals,abs(torque_mst(:,3)))
figure
loglog(Nvals,abs(force_mst(:,1)))
hold on
loglog(Nvals,abs(force_mst(:,2)))
loglog(Nvals,abs(force_mst(:,3)))
\ No newline at end of file
close all; clc; clear;
addpath(genpath("../../../../"));
format long;
global X;
global W;
load('X3','X');
load('W3','W');
rng(10);
A = rand(3,3);
[Qrot,R] = qr(A);
% Initializing parameters for the problem
Nvals = 100:50:1700;
sz = size(Nvals,2);
torque_mst = zeros(sz,3);
force_mst = zeros(sz,3);
linferrs = zeros(sz,1);
Tn_nearest = zeros(sz,1);
Tn_plane = zeros(sz,1);
% Torus radii
r1 = 10;
r2 = 3;
N = 60;
for ii = 1:sz
N = Nvals(ii)
% distance between the centers
dist = 40;
% Mesh for the geometry
mesh_in = mshTorus(N,r1,r2);
mesh_out = mesh_in;
% Translate the outer torus
N_vtcs = size(mesh_out.vtx,1);
trans = ones(N_vtcs,1) * [dist 0 0];
mesh_out.vtx = mesh_out.vtx + trans;
% Rotate the inner torus
mesh_in.vtx = mesh_in.vtx * Qrot;
% Join to create the final mesh
mesh = union(mesh_in,mesh_out);
figure;
plot(mesh);
title('Mesh');
%%
% Definition of FEM spaces and integration rule
S1_Gamma = fem(mesh,'P1');
S0_Gamma = fem(mesh,'P0');
S0_Gamma_in = fem(mesh_in,'P0');
order = 3;
Gamma = dom(mesh,order);
Gamma_D = dom(mesh_in,order);
Op_in = restriction(S0_Gamma,mesh_in);
Op_out = restriction(S0_Gamma,mesh_out);
% Getting the Single Layer matrix V using Sauter Schwab
kernel = @(x,y,z) 1./vecnorm(z,2,2)/ (4*pi);
V = panel_oriented_assembly(mesh,kernel,S0_Gamma,S0_Gamma);
% l vector for fixed charge formulation
l = integral(Gamma_D, S0_Gamma);
% Block matrix
Q = 10;
Vblock = [V l; l' 0];
temp = Vblock\[0*l; Q];
Psi = temp(1:S0_Gamma.ndof);
c = temp(S0_Gamma.ndof+1);
% Psi in
Psi_in = Op_in * Psi;
% Visualizing the Neumann trace
dofs = S0_Gamma.dof;
normals = mesh.nrm;
sPsi = 2*Psi;
figure;
plot(mesh);
hold on;
quiver3(dofs(:,1),dofs(:,2),dofs(:,3),normals(:,1).*sPsi,normals(:,2).*sPsi,normals(:,3).*sPsi,0);
title('Field on the surface');
%quiver3(dofs(:,1),dofs(:,2),dofs(:,3),normals(:,1),normals(:,2),normals(:,3));
% Visualizing the charge density on the surface of the objects
figure;
plot(mesh);
hold on;
plot(mesh,Psi);
title("Surface charge density");
colorbar;
%% Torque evaluation using maxwell stress tensor on the inner torus
Xcg = [0 0 0];
dofs = S0_Gamma_in.dof;
Xvec = dofs-(dofs(:,1)==dofs(:,1))*Xcg;
normals = mesh_in.nrm;
torque_mst(ii,:) = 0.5 * sum( mesh_in.ndv .* Psi_in.^2 .* cross(Xvec,normals) ,1)
force_mst(ii,:) = 0.5 * sum( mesh_in.ndv .* Psi_in.^2 .* normals ,1)
%%
% Choice of velocity field for force computation
% Input = N X 3, output = N X 3
% Velocity fields for the far away sphere
%Nux = @(X) (sum(X.*X,2)>R*R).*[X(:,1)==X(:,1), 0 * X(:,1) , 0 * X(:,1)];
%Nuy = @(X) (sum(X.*X,2)>R*R).*[0*X(:,1), X(:,2)==X(:,2) , 0 * X(:,1)];
%Nuz = @(X) (sum(X.*X,2)>R*R).*[0* X(:,1), 0 * X(:,1) , X(:,3)==X(:,3)];
% Velocity fields for the near sphere
% Cutoff radius
R = 15;
Nux = @(X) (sum(X.*X,2)<R*R).* cross([X(:,1)==X(:,1), 0*X(:,1), 0*X(:,1)],X-Xcg);
Nuy = @(X) (sum(X.*X,2)<R*R).* cross([0*X(:,1), X(:,1)==X(:,1), 0*X(:,1)],X-Xcg);
Nuz = @(X) (sum(X.*X,2)<R*R).* cross([0*X(:,1), 0*X(:,1), X(:,1)==X(:,1)],X-Xcg);
% Kernels
kernelx = @(x,y,z) sum(z.*(Nux(x) - Nux(y)), 2)./(vecnorm(z,2,2).^3)/ (4*pi);
kernely = @(x,y,z) sum(z.*(Nuy(x) - Nuy(y)), 2)./(vecnorm(z,2,2).^3)/ (4*pi);
kernelz = @(x,y,z) sum(z.*(Nuz(x) - Nuz(y)), 2)./(vecnorm(z,2,2).^3)/ (4*pi);
% kernels = cell(3,1);
% kernels = {kernelx,kernely,kernelz};
% torques = cell(3,1);
% t2mats = cell(3,1);
%
% parfor i = 1:3
% t2mats(i) = panel_oriented_assembly(mesh,kernels(i),S0_Gamma,S0_Gamma);
% torques(i) = dot(Psi,t2mats(i)*Psi);
% end
% t2matx = panel_oriented_assembly(mesh,kernelx,S0_Gamma,S0_Gamma);
% torquex = 0.5 * dot(Psi,t2matx*Psi);
% Visualizing the velocity fields
figure;
plot(mesh);
hold on;
vtcs = S0_Gamma.dof;
vels = Nux(vtcs);
quiver3(vtcs(:,1),vtcs(:,2),vtcs(:,3),vels(:,1),vels(:,2),vels(:,3));
title('Perturbation field');
end
save('fp_tor_tor.mat',"torque_mst","force_mst","Tn_plane","Tn_nearest","Nvals","linferrs");
\ No newline at end of file
close all; clc; clear;
addpath(genpath("../../../../"));
format long;
global X;
global W;
load('X3','X');
load('W3','W');
rng(10);
A = rand(3,3);
[Q,R] = qr(A);
% Initializing parameters for the problem
Nvals = 100:50:1700;
sz = size(Nvals,2);
torque_mst = zeros(sz,1);
force_mst = zeros(sz,1);
linferrs = zeros(sz,1);
Tn_nearest = zeros(sz,1);
Tn_plane = zeros(sz,1);
% Torus radii
r1 = 10;
r2 = 3;
N = 60;
for ii = 1:sz
N = Nvals(ii)
% distance between the centers
dist = 40;
% Mesh for the geometry
mesh_in = mshTorus(N,r1,r2);
mesh_out = mesh_in;
% Translate the outer torus
N_vtcs = size(mesh_out.vtx,1);
trans = ones(N_vtcs,1) * [dist 0 0];
mesh_out.vtx = mesh_out.vtx + trans;
% Rotate the inner torus
mesh_in.vtx = mesh_in.vtx * Q;
% Join to create the final mesh
mesh = union(mesh_in,mesh_out);
figure;
plot(mesh);
title('Mesh');
%%
% Definition of FEM spaces and integration rule
S1_Gamma = fem(mesh,'P1');
S0_Gamma = fem(mesh,'P0');
S0_Gamma_in = fem(mesh_in,'P0');
order = 3;
Gamma = dom(mesh,order);
Gamma_D = dom(mesh_in,order);
Op_in = restriction(S0_Gamma,mesh_in);
Op_out = restriction(S0_Gamma,mesh_out);
% Getting the Single Layer matrix V using Gypsilab implementation
Gxy = @(X,Y)femGreenKernel(X,Y,'[1/r]',0); % 0 wave number
V = 1/(4*pi)*integral(Gamma,Gamma,S0_Gamma,Gxy,S0_Gamma);
V = V + 1/(4*pi)*regularize(Gamma,Gamma,S0_Gamma,'[1/r]',S0_Gamma);
% l vector for fixed charge formulation
l = integral(Gamma_D, S0_Gamma);
% Block matrix
Q = 1;
Vblock = [V l; l' 0];
temp = Vblock\[0*l; Q];
Psi = temp(1:S0_Gamma.ndof);
c = temp(S0_Gamma.ndof+1);
% Psi in
Psi_in = Op_in * Psi;
% Visualizing the Neumann trace
dofs = S0_Gamma.dof;
normals = mesh.nrm;
sPsi = 2*Psi;
figure;
plot(mesh);
hold on;
quiver3(dofs(:,1),dofs(:,2),dofs(:,3),normals(:,1).*sPsi,normals(:,2).*sPsi,normals(:,3).*sPsi,0);
title('Field on the surface');
%quiver3(dofs(:,1),dofs(:,2),dofs(:,3),normals(:,1),normals(:,2),normals(:,3));
% Visualizing the charge density on the surface of the objects
figure;
plot(mesh);
hold on;
plot(mesh,Psi);
title("Surface charge density");
colorbar;
%% Torque evaluation using maxwell stress tensor on the inner torus
Xcg = [0 0 0];
dofs = S0_Gamma_in.dof;
Xvec = dofs-(dofs(:,1)==dofs(:,1))*Xcg;
normals = mesh_in.nrm;
torque_mst(ii) = 0.5 * sum( mesh_in.ndv .* Psi_in.^2 .* cross(Xvec,normals) ,1)
force_mst(ii) = 0.5 * sum( mesh_in.ndv .* Psi_in.^2 .* normals ,1)
%%
% Choice of velocity field for force computation
% Input = N X 3, output = N X 3
% Velocity fields for the far away sphere
%Nux = @(X) (sum(X.*X,2)>R*R).*[X(:,1)==X(:,1), 0 * X(:,1) , 0 * X(:,1)];
%Nuy = @(X) (sum(X.*X,2)>R*R).*[0*X(:,1), X(:,2)==X(:,2) , 0 * X(:,1)];
%Nuz = @(X) (sum(X.*X,2)>R*R).*[0* X(:,1), 0 * X(:,1) , X(:,3)==X(:,3)];
% Velocity fields for the near sphere
% Cutoff radius
R = 10;
Nux = @(X) (sum(X.*X,2)<R*R).* cross([X(:,1)==X(:,1), 0*X(:,1), 0*X(:,1)],X-Xcg);
Nuy = @(X) (sum(X.*X,2)<R*R).* cross([0*X(:,1), X(:,1)==X(:,1), 0*X(:,1)],X-Xcg);
Nuz = @(X) (sum(X.*X,2)<R*R).* cross([0*X(:,1), 0*X(:,1), X(:,1)==X(:,1)],X-Xcg);
% Kernels
kernelx = @(x,y,z) sum(z.*(Nux(x) - Nux(y)), 2)./(vecnorm(z,2,2).^3)/ (4*pi);
kernely = @(x,y,z) sum(z.*(Nuy(x) - Nuy(y)), 2)./(vecnorm(z,2,2).^3)/ (4*pi);
kernelz = @(x,y,z) sum(z.*(Nuz(x) - Nuz(y)), 2)./(vecnorm(z,2,2).^3)/ (4*pi);