force_cube_cube_ext.m 5.02 KB
Newer Older
ppanchal's avatar
ppanchal committed
1
2
3
4
5
6
%function [] = force_cube_cube_ext(N)

close all; clc; clear;
addpath(genpath("../../"));

% Initializing parameters for the problem
ppanchal's avatar
ppanchal committed
7
8
9
N = 50;
%N  = getenv('TESTVAR');
%disp(N);
ppanchal's avatar
ppanchal committed
10
11
12
13
14
15
16
17
% Dimensions of the cube
L = [1,1,1];
%distances = [distances; offset];

% Mesh for the geometry
cube_vol_mesh = mshCube(N,L);
mesh_in = cube_vol_mesh.bnd;

ppanchal's avatar
ppanchal committed
18
cube_vol_mesh = mshCube(N,L);
ppanchal's avatar
ppanchal committed
19
20
21
mesh_out = cube_vol_mesh.bnd;

% Translate the outer cube
ppanchal's avatar
ppanchal committed
22
23
24
25
tx = 1;
ty = 3;
tz = 2;
dist = norm([tx ty tz]);
ppanchal's avatar
ppanchal committed
26
N_vtcs = size(mesh_out.vtx,1);
ppanchal's avatar
ppanchal committed
27
trans = [tx * ones(N_vtcs,1), ty * ones(N_vtcs,1), tz * ones(N_vtcs,1)];
ppanchal's avatar
ppanchal committed
28
29
30
31
32
33
mesh_out.vtx = mesh_out.vtx + trans;

% Join to create the final mesh
mesh = union(mesh_in,mesh_out);

figure;
ppanchal's avatar
ppanchal committed
34
35
plot(mesh);
title('Mesh');
ppanchal's avatar
ppanchal committed
36

ppanchal's avatar
ppanchal committed
37
%% BEM 
ppanchal's avatar
ppanchal committed
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68

% 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);
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
ppanchal's avatar
ppanchal committed
69
70
R = 1.1; % Cutoff radius, also used for the velocity field
g = @(X) (sqrt(sum(X.^2,2)) > R)*(1) + (sqrt(sum(X.^2,2)) <= R)*3;
ppanchal's avatar
ppanchal committed
71

ppanchal's avatar
ppanchal committed
72
73
74
75
76
77
figure;
plot(mesh);
hold on;
plot(mesh,g(S0_Gamma.dof));
title('Dirichlet Boundary Condition');
colorbar;
ppanchal's avatar
ppanchal committed
78
79
80
81
82
83
84

% Constructing the RHS
g_N = integral(Gamma,S0_Gamma,g);

% Exterior problem
Psi = V\((-0.5 * M + K)* (M\g_N));

ppanchal's avatar
ppanchal committed
85
% Getting Neumann trace on the two distinct surfaces
ppanchal's avatar
ppanchal committed
86
87
88
89
90
91
92
93
94
95
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;
ppanchal's avatar
ppanchal committed
96
97
98
figure;
plot(mesh);
hold on;
ppanchal's avatar
ppanchal committed
99
quiver3(dofs(:,1),dofs(:,2),dofs(:,3),normals(:,1).*sPsi,normals(:,2).*sPsi,normals(:,3).*sPsi,0);
ppanchal's avatar
ppanchal committed
100
title('Field on the surface');
ppanchal's avatar
ppanchal committed
101
102
%quiver3(dofs(:,1),dofs(:,2),dofs(:,3),normals(:,1),normals(:,2),normals(:,3));

ppanchal's avatar
ppanchal committed
103
% Visualizing the charge density on the surface of the objects
ppanchal's avatar
ppanchal committed
104
figure;
ppanchal's avatar
ppanchal committed
105
plot(mesh);
ppanchal's avatar
ppanchal committed
106
hold on;
ppanchal's avatar
ppanchal committed
107
108
109
plot(mesh,Psi);
title("Surface charge density");
colorbar;
ppanchal's avatar
ppanchal committed
110
111
112
113
114
115
116
117
118
%%
% Classical formula for force evaluation
normals_in = mesh_in.nrm;
classical_force_in = 0.5 * sum( mesh_in.ndv .* Psi_in.^2.* normals_in, 1)

normals_out = mesh_out.nrm;

classical_force_out = 0.5 * sum( mesh_out.ndv .* Psi_out.^2.* normals_out, 1)

ppanchal's avatar
ppanchal committed
119
120
121
122
123
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


ppanchal's avatar
ppanchal committed
124
125
126
127
%%
% Choice of velocity field for force computation
% 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)];
ppanchal's avatar
ppanchal committed
128
129
130
131
132
133
134
135
136
137
138
139

% Velocity fields for the far away 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)];

% 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)];


ppanchal's avatar
ppanchal committed
140
141
142
143
144
%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).*[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)];

ppanchal's avatar
ppanchal committed
145
146
147
148
149
150
151
152
153
% 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');

ppanchal's avatar
ppanchal committed
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
% 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);
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,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(t2matx))

forcex = dot(Psi,t2matx * Rho)
forcey = dot(Psi,t2maty * Rho)
forcez = dot(Psi,t2matz * Rho)
%forcevals = [forcevals; force];

ppanchal's avatar
ppanchal committed
170
171
172
173
str1 = "Cube_Cube_";
fname = append(str1,int2str(N));
save(fname);
exit;
ppanchal's avatar
ppanchal committed
174
%end