Commit 92300c5c authored by ppanchal's avatar ppanchal
Browse files

Cube Cube computations on ada

parent 05b529b9
function out = G(n,alpha)
out = exp(-(n+0.5).*alpha)./sinh((n+0.5)*alpha);
end
\ No newline at end of file
function out = Gsum(n,alpha)
% nidx = 0:n;
% nidx2 = 2 * nidx;
% Gvec = G(nidx', alpha * (nidx'==nidx'));
%
% out = sum(Gvec.*((2*nidx'+1).*Gvec - )
out = 0;
for i = 0:n
out = out + G(i,alpha) * ( (2*i+1)*G(i,alpha) - (2*i+2)*G(i+1,alpha) );
end
end
\ No newline at end of file
function out = correction(n)
for i = 0:n
for j = 0:n
end
end
end
\ No newline at end of file
% Testing the force formula
% Testing convergence for some alpha values
alpha1 = acosh(10);
alpha2 = acosh(20);
force1 = [];
force2 = [];
for i = 1 : 500
force1 = [force1 Gsum(i,alpha1)];
force2 = [force2 Gsum(i,alpha2)];
end
figure;
plot(1:500,force1);
hold on;
plot(1:500,force2);
%alpha = acosh(1+x);
% Next plot
xvals = 0.001:0.001:0.01;
Gvals = [];
for i = 1 : size(xvals,2)
Gvals = [Gvals Gsum(100,acosh(1+xvals(i)))];
end
figure;
plot(xvals,Gvals);
\ No newline at end of file
% Testing the force formula
% Testing convergence for some alpha values
alpha1 = acosh(10);
alpha2 = acosh(20);
force1 = [];
force2 = [];
for i = 1 : 500
force1 = [force1 Gsum(i,alpha1)];
force2 = [force2 Gsum(i,alpha2)];
end
figure;
plot(1:500,force1);
hold on;
plot(1:500,force2);
%alpha = acosh(1+x);
% Next plot
xvals = 0.001:0.001:0.01;
%function [] = force_cube_cube_ext(N) %function [] = force_cube_cube_ext(N)
close all; clc; clear; close all; clc; clear;
%distances = [];
%forcevals = [];
addpath(genpath("../../")); addpath(genpath("../../"));
% Initializing parameters for the problem % Initializing parameters for the problem
N = 20; N = 50;
%N = getenv('TESTVAR');
%disp(N);
% Dimensions of the cube % Dimensions of the cube
L = [1,1,1]; L = [1,1,1];
%distances = [distances; offset]; %distances = [distances; offset];
...@@ -16,21 +15,26 @@ L = [1,1,1]; ...@@ -16,21 +15,26 @@ L = [1,1,1];
cube_vol_mesh = mshCube(N,L); cube_vol_mesh = mshCube(N,L);
mesh_in = cube_vol_mesh.bnd; mesh_in = cube_vol_mesh.bnd;
cube_vol_mesh = mshCube(N/2.,L/2.); cube_vol_mesh = mshCube(N,L);
mesh_out = cube_vol_mesh.bnd; mesh_out = cube_vol_mesh.bnd;
% Translate the outer cube % Translate the outer cube
tx = 1;
ty = 3;
tz = 2;
dist = norm([tx ty tz]);
N_vtcs = size(mesh_out.vtx,1); N_vtcs = size(mesh_out.vtx,1);
trans = [ones(N_vtcs,1), 3 * ones(N_vtcs,1), 2 * ones(N_vtcs,1)]; trans = [tx * ones(N_vtcs,1), ty * ones(N_vtcs,1), tz * ones(N_vtcs,1)];
mesh_out.vtx = mesh_out.vtx + trans; mesh_out.vtx = mesh_out.vtx + trans;
% 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); hold on; plot(mesh);
title('Mesh');
%% %% BEM
% Definition of FEM spaces and integration rule % Definition of FEM spaces and integration rule
S1_Gamma = fem(mesh,'P1'); S1_Gamma = fem(mesh,'P1');
...@@ -62,68 +66,47 @@ K = K +1/(4*pi)*regularize(Gamma,Gamma,S0_Gamma,'grady[1/r]',ntimes(S0_Gamma)); ...@@ -62,68 +66,47 @@ K = K +1/(4*pi)*regularize(Gamma,Gamma,S0_Gamma,'grady[1/r]',ntimes(S0_Gamma));
M = integral(Gamma,S0_Gamma,S0_Gamma); M = integral(Gamma,S0_Gamma,S0_Gamma);
% Defining the Dirichlet boundary condition % Defining the Dirichlet boundary condition
R = 1; R = 1.1; % Cutoff radius, also used for the velocity field
g = @(X) (sqrt(sum(X.^2,2)) > R*R)*(1) + (sqrt(sum(X.^2,2)) <= R*R)*3; g = @(X) (sqrt(sum(X.^2,2)) > R)*(1) + (sqrt(sum(X.^2,2)) <= R)*3;
%g = @(X) 0 * sum(X,2); figure;
%g = @(X) X(:,1)+X(:,2)+X(:,3); plot(mesh);
% g = @(X) X(:,1); hold on;
%g = @(X) 1./sqrt(sum(X.^2,2)); plot(mesh,g(S0_Gamma.dof));
title('Dirichlet Boundary Condition');
colorbar;
% Constructing the RHS % Constructing the RHS
% g_N = g(S0_Gamma.dof);
g_N = integral(Gamma,S0_Gamma,g); g_N = integral(Gamma,S0_Gamma,g);
% PI_g = M\g_N;
% Checking the boundary condition
%PI_g = g(S0_Gamma.dof);
%plot(mesh,PI_g)
% Solving for state solution
% Interior problem
%Psi = V\((0.5*M + K)*(M\g_N));
gradU{1} = @(X)(-X(:,1)./(X(:,1).^2 + X(:,2).^2 +X(:,3).^2 ).^(3/2));
gradU{2} = @(X)(-X(:,2)./(X(:,1).^2 + X(:,2).^2 +X(:,3).^2 ).^(3/2));
gradU{3} = @(X)(-X(:,3)./(X(:,1).^2 + X(:,2).^2 +X(:,3).^2 ).^(3/2));
Psi_exact = M\integral(Gamma,ntimes(S0_Gamma),gradU);
% Exterior problem % Exterior problem
% More exact RHS (Martin's way)
Psi = V\((-0.5 * M + K)* (M\g_N)); Psi = V\((-0.5 * M + K)* (M\g_N));
% Approximating the RHS
%Psi = V\((-0.5 * M + K)* g_N);
% Getting Neumann trace on the two distinct surfaces
Psi_in = Op_in * Psi; Psi_in = Op_in * Psi;
Psi_out = Op_out * Psi; Psi_out = Op_out * Psi;
% Solving the adjoint problem to get the adjoint solution % Solving the adjoint problem to get the adjoint solution
% More exact way for RHS (Martin)
Rho = V\(-0.5 * g_N); Rho = V\(-0.5 * g_N);
%Rho = V\(-0.5 * M * g_N);
% Visualizing the Neumann trace % Visualizing the Neumann trace
dofs = S0_Gamma.dof; dofs = S0_Gamma.dof;
normals = mesh.nrm; normals = mesh.nrm;
sPsi = -2*Psi; 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); 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)); %quiver3(dofs(:,1),dofs(:,2),dofs(:,3),normals(:,1),normals(:,2),normals(:,3));
% Checking the solutions % Visualizing the charge density on the surface of the objects
%figure;
%plot(Psi_in);
%hold on;
%plot(Psi_out);
figure; figure;
plot(Psi); plot(mesh);
hold on; hold on;
% centers = mesh.ctr; plot(mesh,Psi);
%exact_sol = - sum(centers .* normals,2)./(sqrt(sum(centers.^2,2))).^3; title("Surface charge density");
%plot(exact_solution); colorbar;
legend(['Psi','Exact solution']);
% plot(sum(normals,2));
plot(Psi_exact);
norm(Psi-Psi_exact)
%% %%
% Classical formula for force evaluation % Classical formula for force evaluation
normals_in = mesh_in.nrm; normals_in = mesh_in.nrm;
...@@ -133,18 +116,41 @@ normals_out = mesh_out.nrm; ...@@ -133,18 +116,41 @@ normals_out = mesh_out.nrm;
classical_force_out = 0.5 * sum( mesh_out.ndv .* Psi_out.^2.* normals_out, 1) classical_force_out = 0.5 * sum( mesh_out.ndv .* Psi_out.^2.* normals_out, 1)
charge_in = sum( mesh_in.ndv .* Psi_in, 1)
charge_out = sum( mesh_out.ndv .* Psi_out, 1)
coulomb_force = charge_in * charge_out / 4 / pi / dist^2
%% %%
% Choice of velocity field for force computation % Choice of velocity field for force computation
% Input = N X 3, output = N X 3 % Input = N X 3, output = N X 3
%Nu = @(X) (sum(X.*X,2)<=R*R).*[X(:,1)==X(:,1), 0 * X(:,1) , 0 * X(:,1)]; %Nu = @(X) (sum(X.*X,2)<=R*R).*[X(:,1)==X(:,1), 0 * X(:,1) , 0 * X(:,1)];
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)]; % Velocity fields for the far away cube
Nuz = @(X) (sum(X.*X,2)>R*R).*[0* X(:,1), 0 * X(:,1) , X(:,3)==X(:,3)]; %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 cube
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)];
%Nu = @(X) (sum(X.*X,2)>R*R).*[X(:,1)==X(:,1), 0 * X(:,1) , 0 * X(:,1)]; %Nu = @(X) (sum(X.*X,2)>R*R).*[X(:,1)==X(:,1), 0 * X(:,1) , 0 * X(:,1)];
%Nu = @(X) (sum(X.*X,2)<=R*R).*[0 * X(:,1), X(:,1)==X(:,1) , 0 * X(:,1)]; %Nu = @(X) (sum(X.*X,2)<=R*R).*[0 * X(:,1), X(:,1)==X(:,1) , 0 * X(:,1)];
%Nu = @(X) (sum(X.*X,2)<=R*R).*[X(:,1)==X(:,1), 0 * X(:,1) , 0 * X(:,1)]; %Nu = @(X) (sum(X.*X,2)<=R*R).*[X(:,1)==X(:,1), 0 * X(:,1) , 0 * X(:,1)];
%Nu = @(X) (sum(X.*X,2)<=R*R).*[X(:,1)==X(:,1), 0 * X(:,1) , 0 * X(:,1)]; %Nu = @(X) (sum(X.*X,2)<=R*R).*[X(:,1)==X(:,1), 0 * X(:,1) , 0 * X(:,1)];
% 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');
% Definition of the kernel for T2 % Definition of the kernel for T2
kernelx = @(x,y,z) sum(z.*(Nux(x) - Nux(y)), 2)./(sqrt(sum(z.^2,2)).^3)/ (4*pi); kernelx = @(x,y,z) sum(z.*(Nux(x) - Nux(y)), 2)./(sqrt(sum(z.^2,2)).^3)/ (4*pi);
kernely = @(x,y,z) sum(z.*(Nuy(x) - Nuy(y)), 2)./(sqrt(sum(z.^2,2)).^3)/ (4*pi); kernely = @(x,y,z) sum(z.*(Nuy(x) - Nuy(y)), 2)./(sqrt(sum(z.^2,2)).^3)/ (4*pi);
...@@ -161,4 +167,8 @@ forcey = dot(Psi,t2maty * Rho) ...@@ -161,4 +167,8 @@ forcey = dot(Psi,t2maty * Rho)
forcez = dot(Psi,t2matz * Rho) forcez = dot(Psi,t2matz * Rho)
%forcevals = [forcevals; force]; %forcevals = [forcevals; force];
str1 = "Cube_Cube_";
fname = append(str1,int2str(N));
save(fname);
exit;
%end %end
\ No newline at end of file
%function [] = force_cube_cube_ext(N) %function [] = force_cube_cube_ext(N)
close all; clc; clear; close all; clc; clear;
%distances = [];
%forcevals = [];
addpath(genpath("../../")); addpath(genpath("../../"));
% Initializing parameters for the problem % Initializing parameters for the problem
N = 20; %N = 20;
N = getenv('TESTVAR');
% Dimensions of the cube % Dimensions of the cube
L = [1,1,1]; L = [1,1,1];
%distances = [distances; offset]; %distances = [distances; offset];
...@@ -20,17 +18,22 @@ cube_vol_mesh = mshCube(N/2.,L/2.); ...@@ -20,17 +18,22 @@ cube_vol_mesh = mshCube(N/2.,L/2.);
mesh_out = cube_vol_mesh.bnd; mesh_out = cube_vol_mesh.bnd;
% Translate the outer cube % Translate the outer cube
tx = 1;
ty = 3;
tz = 2;
dist = norm([tx ty tz]);
N_vtcs = size(mesh_out.vtx,1); N_vtcs = size(mesh_out.vtx,1);
trans = [ones(N_vtcs,1), 3 * ones(N_vtcs,1), 2 * ones(N_vtcs,1)]; trans = [tx * ones(N_vtcs,1), ty * ones(N_vtcs,1), tz * ones(N_vtcs,1)];
mesh_out.vtx = mesh_out.vtx + trans; mesh_out.vtx = mesh_out.vtx + trans;
% 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); hold on; plot(mesh);
title('Mesh');
%% %% BEM
% Definition of FEM spaces and integration rule % Definition of FEM spaces and integration rule
S1_Gamma = fem(mesh,'P1'); S1_Gamma = fem(mesh,'P1');
...@@ -62,68 +65,47 @@ K = K +1/(4*pi)*regularize(Gamma,Gamma,S0_Gamma,'grady[1/r]',ntimes(S0_Gamma)); ...@@ -62,68 +65,47 @@ K = K +1/(4*pi)*regularize(Gamma,Gamma,S0_Gamma,'grady[1/r]',ntimes(S0_Gamma));
M = integral(Gamma,S0_Gamma,S0_Gamma); M = integral(Gamma,S0_Gamma,S0_Gamma);
% Defining the Dirichlet boundary condition % Defining the Dirichlet boundary condition
R = 1; R = 1.1; % Cutoff radius, also used for the velocity field
g = @(X) (sqrt(sum(X.^2,2)) > R*R)*(1) + (sqrt(sum(X.^2,2)) <= R*R)*3; g = @(X) (sqrt(sum(X.^2,2)) > R)*(1) + (sqrt(sum(X.^2,2)) <= R)*3;
%g = @(X) 0 * sum(X,2); figure;
%g = @(X) X(:,1)+X(:,2)+X(:,3); plot(mesh);
% g = @(X) X(:,1); hold on;
%g = @(X) 1./sqrt(sum(X.^2,2)); plot(mesh,g(S0_Gamma.dof));
title('Dirichlet Boundary Condition');
colorbar;
% Constructing the RHS % Constructing the RHS
% g_N = g(S0_Gamma.dof);
g_N = integral(Gamma,S0_Gamma,g); g_N = integral(Gamma,S0_Gamma,g);
% PI_g = M\g_N;
% Checking the boundary condition
%PI_g = g(S0_Gamma.dof);
%plot(mesh,PI_g)
% Solving for state solution
% Interior problem
%Psi = V\((0.5*M + K)*(M\g_N));
gradU{1} = @(X)(-X(:,1)./(X(:,1).^2 + X(:,2).^2 +X(:,3).^2 ).^(3/2));
gradU{2} = @(X)(-X(:,2)./(X(:,1).^2 + X(:,2).^2 +X(:,3).^2 ).^(3/2));
gradU{3} = @(X)(-X(:,3)./(X(:,1).^2 + X(:,2).^2 +X(:,3).^2 ).^(3/2));
Psi_exact = M\integral(Gamma,ntimes(S0_Gamma),gradU);
% Exterior problem % Exterior problem
% More exact RHS (Martin's way)
Psi = V\((-0.5 * M + K)* (M\g_N)); Psi = V\((-0.5 * M + K)* (M\g_N));
% Approximating the RHS
%Psi = V\((-0.5 * M + K)* g_N);
% Getting Neumann trace on the two distinct surfaces
Psi_in = Op_in * Psi; Psi_in = Op_in * Psi;
Psi_out = Op_out * Psi; Psi_out = Op_out * Psi;
% Solving the adjoint problem to get the adjoint solution % Solving the adjoint problem to get the adjoint solution
% More exact way for RHS (Martin)
Rho = V\(-0.5 * g_N); Rho = V\(-0.5 * g_N);
%Rho = V\(-0.5 * M * g_N);
% Visualizing the Neumann trace % Visualizing the Neumann trace
dofs = S0_Gamma.dof; dofs = S0_Gamma.dof;
normals = mesh.nrm; normals = mesh.nrm;
sPsi = -2*Psi; 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); 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)); %quiver3(dofs(:,1),dofs(:,2),dofs(:,3),normals(:,1),normals(:,2),normals(:,3));
% Checking the solutions % Visualizing the charge density on the surface of the objects
%figure;
%plot(Psi_in);
%hold on;
%plot(Psi_out);
figure; figure;
plot(Psi); plot(mesh);
hold on; hold on;
% centers = mesh.ctr; plot(mesh,Psi);
%exact_sol = - sum(centers .* normals,2)./(sqrt(sum(centers.^2,2))).^3; title("Surface charge density");
%plot(exact_solution); colorbar;
legend(['Psi','Exact solution']);
% plot(sum(normals,2));
plot(Psi_exact);
norm(Psi-Psi_exact)
%% %%
% Classical formula for force evaluation % Classical formula for force evaluation
normals_in = mesh_in.nrm; normals_in = mesh_in.nrm;
...@@ -133,26 +115,58 @@ normals_out = mesh_out.nrm; ...@@ -133,26 +115,58 @@ normals_out = mesh_out.nrm;
classical_force_out = 0.5 * sum( mesh_out.ndv .* Psi_out.^2.* normals_out, 1) classical_force_out = 0.5 * sum( mesh_out.ndv .* Psi_out.^2.* normals_out, 1)
charge_in = sum( mesh_in.ndv .* Psi_in, 1)
charge_out = sum( mesh_out.ndv .* Psi_out, 1)
coulomb_force = charge_in * charge_out / 4 / pi / dist^2
%% %%
% Choice of velocity field for force computation % Choice of velocity field for force computation
% Input = N X 3, output = N X 3 % Input = N X 3, output = N X 3
%Nu = @(X) (sum(X.*X,2)<=R*R).*[X(:,1)==X(:,1), 0 * X(:,1) , 0 * X(:,1)]; %Nu = @(X) (sum(X.*X,2)<=R*R).*[X(:,1)==X(:,1), 0 * X(:,1) , 0 * X(:,1)];
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)]; % Velocity fields for the far away cube
Nuz = @(X) (sum(X.*X,2)>R*R).*[0* X(:,1), 0 * X(:,1) , X(:,3)==X(:,3)]; %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 cube
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)];
%Nu = @(X) (sum(X.*X,2)>R*R).*[X(:,1)==X(:,1), 0 * X(:,1) , 0 * X(:,1)]; %Nu = @(X) (sum(X.*X,2)>R*R).*[X(:,1)==X(:,1), 0 * X(:,1) , 0 * X(:,1)];
%Nu = @(X) (sum(X.*X,2)<=R*R).*[0 * X(:,1), X(:,1)==X(:,1) , 0 * X(:,1)]; %Nu = @(X) (sum(X.*X,2)<=R*R).*[0 * X(:,1), X(:,1)==X(:,1) , 0 * X(:,1)];
%Nu = @(X) (sum(X.*X,2)<=R*R).*[X(:,1)==X(:,1), 0 * X(:,1) , 0 * X(:,1)]; %Nu = @(X) (sum(X.*X,2)<=R*R).*[X(:,1)==X(:,1), 0 * X(:,1) , 0 * X(:,1)];
%Nu = @(X) (sum(X.*X,2)<=R*R).*[X(:,1)==X(:,1), 0 * X(:,1) , 0 * X(:,1)]; %Nu = @(X) (sum(X.*X,2)<=R*R).*[X(:,1)==X(:,1), 0 * X(:,1) , 0 * X(:,1)];
% 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');
% Definition of the kernel for T2 % Definition of the kernel for T2
kernelx = @(x,y,z) sum(z.*(Nu(x) - Nu(y)), 2)./(sqrt(sum(z.^2,2)).^3)/ (4*pi); kernelx = @(x,y,z) sum(z.*(Nux(x) - Nux(y)), 2)./(sqrt(sum(z.^2,2)).^3)/ (4*pi);
kernely = @(x,y,z) sum(z.*(Nuy(x) - Nuy(y)), 2)./(sqrt(sum(z.^2,2)).^3)/ (4*pi);
kernelz = @(x,y,z) sum(z.*(Nuz(x) - Nuz(y)), 2)./(sqrt(sum(z.^2,2)).^3)/ (4*pi);
t2matx = panel_oriented_assembly(mesh,kernel,S0_Gamma,S0_Gamma); t2matx = panel_oriented_assembly(mesh,kernelx,S0_Gamma,S0_Gamma);
t2maty = panel_oriented_assembly(mesh,kernely,S0_Gamma,S0_Gamma);
t2matz = panel_oriented_assembly(mesh,kernelz,S0_Gamma,S0_Gamma);
sum(sum(t2mat)) sum(sum(t2matx))
force = dot(Psi,t2mat * Rho) forcex = dot(Psi,t2matx * Rho)
forcey = dot(Psi,t2maty * Rho)
forcez = dot(Psi,t2matz * Rho)
%forcevals = [forcevals; force]; %forcevals = [forcevals; force];
str1 = "Cube_Cube_";
fname = append(str1
%end %end
\ No newline at end of file
function out = force_identical_spheres(b)
out = 1-4*b^3 - 6 * b^5 + 14 * b^6 - 8 * b^7 + 54 * b^8 - 50 * b^9 + 154 * b^10 - 264 * b^11 + 494 * b^12 - 1092 * b^13 + 1830 * b^14 - 4192 * b^15 + 7140 * b^16 - 15894 * b^17 + 28234 * b^18 - 60320 * b^19 + 112056 * b^20 - 230032 * b^21;
end
\ No newline at end of file
%function [] = force_cube_cube_ext(N)
close all; clc; clear;
addpath(genpath("../../"));
% Initializing parameters for the problem
N = 30;
% Radius of spheres