Commit ea9adf26 authored by ppanchal's avatar ppanchal
Browse files

new script

parent 1c713b07
function torques = compute_bem_forces(mesh,R,Psi)
S0_Gamma = fem(mesh,'P0');
Nux = @(X) (vecnorm(X,2,2)<R).* ones(size(X,1),1)*[1 0 0];
Nuy = @(X) (vecnorm(X,2,2)<R).* ones(size(X,1),1)*[0 1 0];
Nuz = @(X) (vecnorm(X,2,2)<R).* ones(size(X,1),1)*[0 0 1];
% 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);
t2matx = panel_oriented_assembly(mesh,kernelx,S0_Gamma,S0_Gamma);
forcex = 0.5 * dot(Psi,t2matx*Psi);
t2maty = panel_oriented_assembly(mesh,kernely,S0_Gamma,S0_Gamma);
forcey = 0.5 * dot(Psi,t2maty*Psi);
t2matz = panel_oriented_assembly(mesh,kernelz,S0_Gamma,S0_Gamma);
forcez = 0.5 * dot(Psi,t2matz*Psi);
torques = [forcex forcey forcez];
end
\ No newline at end of file
function torques = compute_bem_torques(mesh,R,Xcg,Psi)
S0_Gamma = fem(mesh,'P0');
Nux = @(X) (vecnorm(X,2,2)<R).* cross(ones(size(X,1),1)*[1 0 0],X-Xcg);
Nuy = @(X) (vecnorm(X,2,2)<R).* cross(ones(size(X,1),1)*[0 1 0],X-Xcg);
Nuz = @(X) (vecnorm(X,2,2)<R).* cross(ones(size(X,1),1)*[0 0 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);
t2matx = panel_oriented_assembly(mesh,kernelx,S0_Gamma,S0_Gamma);
torquex = 0.5 * dot(Psi,t2matx*Psi);
t2maty = panel_oriented_assembly(mesh,kernely,S0_Gamma,S0_Gamma);
torquey = 0.5 * dot(Psi,t2maty*Psi);
t2matz = panel_oriented_assembly(mesh,kernelz,S0_Gamma,S0_Gamma);
torquez = 0.5 * dot(Psi,t2matz*Psi);
torques = [torquex torquey torquez];
end
\ No newline at end of file
% This function computes the force and torque on object D
function [torque,force] = compute_mst_forces(mesh,Xcg,Psi)
% Creating the FEM spaces
S0_Gamma = fem(mesh,'P0');
% Creating the r X n vector
dofs = S0_Gamma.dof;
Xvec = dofs-ones(size(dofs,1),1)*Xcg;
normals = mesh.nrm;
torque = 0.5 * sum( mesh.ndv .* Psi.^2 .* cross(Xvec,normals) ,1);
force = 0.5 * sum( mesh.ndv .* Psi.^2 .* normals ,1);
end
\ No newline at end of file
addpath(genpath("../../../../"));
clear;
clc;
format long;
global X;
global W;
load('X3','X');
load('W3','W');
Nvals = 50:200:2000;
sz = size(Nvals,2);
torques_mst = zeros(sz,3);
torques_bem = zeros(sz,3);
forces_mst = zeros(sz,3);
forces_bem = zeros(sz,3);
for i = 1:sz
% Get the mesh
[mesh,mesh_in,mesh_out] = tor_tor_mesh(10,3,Nvals(i),10);
% Solve the floating potential problem on mesh
[Psif,cf] = solve_float_pt_ext(mesh,mesh_in,1,3,'gypsi');
% Definition of FEM spaces and integration rule
S1_Gamma = fem(mesh,'P1');
S0_Gamma = fem(mesh,'P0');
order = 3;
Gamma = dom(mesh,order);
Op_in = restriction(S0_Gamma,mesh_in);
% 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 = 18;
%assert(R < Rad + s);
Vin = 1;
Vout = 0;
g = @(X) (sqrt(sum(X.^2,2)) > R)* Vout + (sqrt(sum(X.^2,2)) <= R) * Vin;
% 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;
end
addpath(genpath("../../../../"));
clear;
clc;
format long;
global X;
global W;
load('X3','X');
load('W3','W');
Nvals = 50:200:2000;
sz = size(Nvals,2);
torques_mst = zeros(sz,3);
torques_bem = zeros(sz,3);
forces_mst = zeros(sz,3);
forces_bem = zeros(sz,3);
for i = 1:sz
% Get the mesh
[mesh,mesh_in,mesh_out] = tor_tor_mesh(10,3,Nvals(i),10);
% Solve the floating potential problem on mesh
[Psi,c] = solve_float_pt_ext(mesh,mesh_in,1,3,'gypsi');
% Solve the floating potential problem on mesh
[Psi1,c1] = solve_float_pt_ext(mesh,mesh_in,1,3,'ss');
err_psi = norm(Psi-Psi1)
err_c = abs(c-c1)
end
addpath(genpath("../../../../"));
clear;
clc;
format long;
global X;
global W;
load('X3','X');
load('W3','W');
Nvals = 50:200:2000;
sz = size(Nvals,2);
torques_mst = zeros(sz,3);
torques_bem = zeros(sz,3);
forces_mst = zeros(sz,3);
forces_bem = zeros(sz,3);
for i = 1:sz
% Get the mesh
[mesh,mesh_in,mesh_out] = tor_tor_mesh(10,3,Nvals(i),10);
% Solve the floating potential problem on mesh
[Psi,c] = solve_float_pt_ext(mesh,mesh_in,1,3,'gypsi');
S0_Gamma = fem(mesh,'P0');
Op_in = restriction(S0_Gamma,mesh_in);
Psi_in = Op_in * Psi;
% Computing the torques and forces using MST
[torque_mst,force_mst] = compute_mst_forces(mesh_in,[0,0,0],Psi_in);
torques_mst(i,:) = torque_mst;
forces_mst(i,:) = force_mst;
% % Computing the torques using BEM formula
% torque_bem = compute_bem_torques(mesh,18,[0,0,0],Psi);
% torques_bem(i,:) = torque_bem;
%
% % Computing the forces using BEM formula
% force_bem = compute_bem_forces(mesh,18,Psi);
% forces_bem(i,:) = force_bem;
end
save('new_script_data.mat','Nvals','forces_mst','torques_bem','torques_mst','forces_bem');
addpath(genpath("../../../../"));
clear;
clc;
format long;
global X;
global W;
load('X3','X');
load('W3','W');
Nvals = 50:200:2000;
sz = size(Nvals,2);
torques_mst = zeros(sz,3);
torques_bem = zeros(sz,3);
forces_mst = zeros(sz,3);
forces_bem = zeros(sz,3);
for i = 1:sz
% Get the mesh
[mesh,mesh_in,mesh_out] = tor_tor_mesh(10,3,Nvals(i),10);
% Solve the floating potential problem on mesh
[Psi,c] = solve_float_pt_ext(mesh,mesh_in,1,3,'gypsi');
S0_Gamma = fem(mesh,'P0');
Op_in = restriction(S0_Gamma,mesh_in);
Psi_in = Op_in * Psi;
% Computing the torques and forces using MST
[torque_mst,force_mst] = compute_mst_forces(mesh_in,[0,0,0],Psi_in);
torques_mst(i,:) = torque_mst;
forces_mst(i,:) = force_mst;
% % Computing the torques using BEM formula
% torque_bem = compute_bem_torques(mesh,18,[0,0,0],Psi);
% torques_bem(i,:) = torque_bem;
%
% % Computing the forces using BEM formula
% force_bem = compute_bem_forces(mesh,18,Psi);
% forces_bem(i,:) = force_bem;
end
save('new_script_data.mat','Nvals','forces_mst','torques_bem','torques_mst','forces_bem');
% D is the object with the floating potential / fixed charge
function [Psi,c] = solve_float_pt_ext(mesh,mesh_D,Q,order,solver)
Gamma = dom(mesh,order);
Gamma_D = dom(mesh_D,order);
S0_Gamma = fem(mesh,'P0');
switch solver
case 'gypsi'
% Assemble V
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);
case 'ss'
% Assemble V
kernel = @(x,y,z) 1./vecnorm(z,2,2)/ (4*pi);
V = panel_oriented_assembly(mesh,kernel,S0_Gamma,S0_Gamma);
end
l = integral(Gamma_D,S0_Gamma);
Block_matrix = [V l;
l' 0];
Psi_c = Block_matrix\[0*l; Q];
Psi = Psi_c(1:S0_Gamma.ndof);
c = Psi_c(S0_Gamma.ndof+1);
end
\ 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:100:1700;
sz = size(Nvals,2);
torque_mst = zeros(sz,3);
force_mst = zeros(sz,3);
torque_bem = zeros(sz,3);
% 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","torque_bem","Tn_nearest","Nvals");
\ No newline at end of file
......@@ -12,14 +12,12 @@ A = rand(3,3);
% Initializing parameters for the problem
Nvals = 100:50:1700;
Nvals = 100:100: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);
torque_bem = zeros(sz,3);
% Torus radii
r1 = 10;
......@@ -121,7 +119,7 @@ force_mst(ii,:) = 0.5 * sum( mesh_in.ndv .* Psi_in.^2 .* normals ,1)
% 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);
Nux = @(X) (sum(X.*X,2)<R*R).* cross(X(:,1)*[1 0 0],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);
......@@ -129,18 +127,23 @@ Nuz = @(X) (sum(X.*X,2)<R*R).* cross([0*X(:,1), 0*X(:,1), X(:,1)==X(:,1)],X-Xcg)
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);
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
parfor i = 1:3
disp(i)
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);
% t2maty = panel_oriented_assembly(mesh,kernely,S0_Gamma,S0_Gamma);
% torquey = 0.5 * dot(Psi,t2maty*Psi);
% t2matz = panel_oriented_assembly(mesh,kernelz,S0_Gamma,S0_Gamma);
% torquez = 0.5 * dot(Psi,t2matz*Psi);
% Visualizing the velocity fields
......@@ -153,4 +156,4 @@ 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
save('fp_tor_tor.mat',"torque_mst","force_mst","torque_bem","Nvals");
\ No newline at end of file
close all; clc; clear;
addpath(genpath("../../../../"));
format long;
global X;
global W;
load('X3','X');
load('W3','W');
% load('X','X');
% load('W','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;