Skip to content

GitLab

Commit 3f21861a authored by ppanchal's avatar ppanchal
Browse files

Windows commit

parent 80f55d5c
Hide whitespace changes
Inline Side-by-side
Piyush's codes/3D BEM Sauter Schwab/Exterior Problem/compute_err_ext_sph.m 0 → 100644
View file @ 3f21861a
%function l2_err = compute_err_ext_sph(N)
close all; clc; clear;
addpath(genpath("../../"));
global X;
global W;
load('X3','X');
load('W3','W');
% Initializing parameters for the problem
% Radius of spheres
Rad = 10;
% distance between the centers
% separation
s = Rad;
dist = 2*Rad + s;
% Mesh for the geometry
mesh_in = mshSphere(N,Rad);
mesh_out = mesh_in;
% Translate the outer cube
N_vtcs = size(mesh_out.vtx,1);
trans = [dist * ones(N_vtcs,1), zeros(N_vtcs,1), zeros(N_vtcs,1)];
mesh_out.vtx = mesh_out.vtx + trans;
% Slightly translating the inner sphere
mesh_in.vtx = mesh_in.vtx + 0.1*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));
%% Solving the Laplace problem using BEM
% 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
g = @(x) 1/4/pi./vecnorm(x,2,2); % Point charge at origin
% Checking the boundary condition
figure;
plot(mesh);
hold on;
PI_g = g(S0_Gamma.dof);
plot(mesh,PI_g);
title('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;
% 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;
%% Exact solution
Psi_exact = -1/4/pi./vecnorm(dofs,2,2).^3 .* (dofs.*nrms);
% Visualizing the exactness
figure;
plot(Psi./Psi_exact);
% L2 Error
l2_err = norm(Psi-Psi_exact);
close all;
%end
\ No newline at end of file
Piyush's codes/3D BEM Sauter Schwab/Exterior Problem/double_sphere_exterior_laplace.m 0 → 100644
View file @ 3f21861a
function [l2_err,av_err] = double_sphere_exterior_laplace(N)
addpath(genpath("../../"));
% Initializing parameters for the problem
%N = 50;
% Radius of sphere
Rad = 10;
% Mesh for the geometry
mesh_in = mshSphere(N,Rad);
mesh_out = mshSphere(N,Rad);
N_vtcs = size(mesh_in.vtx,1);
mesh_in.vtx = mesh_in.vtx + 3*ones(N_vtcs,3);
mesh_out.vtx = mesh_out.vtx + 25*ones(N_vtcs,3);
% 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(0,0,0);
%% Solving the Laplace problem using BEM
% Definition of FEM spaces and integration rule
S1_Gamma = fem(mesh,'P1');
S0_Gamma = fem(mesh,'P0');
order = 3;
Gamma = dom(mesh,order);
% 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
g = @(x) 1/4/pi./vecnorm(x,2,2); % Point charge at origin
% Checking the boundary condition
figure;
plot(mesh);
hold on;
PI_g = g(S0_Gamma.dof);
plot(mesh,PI_g);
title('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 = V\((0.5 * M + K)* (M\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;
%% Exact solution
Psi_exact = -1/4/pi./vecnorm(dofs,2,2).^3 .* (dot(dofs,nrms,2));
% Visualizing the exactness
figure;
plot(Psi./Psi_exact);
err_vec = Psi-Psi_exact;
% L2 Error
l2_err = norm(err_vec)
av_err = err_vec' * V * err_vec
close all;
end
\ No newline at end of file
Piyush's codes/3D BEM Sauter Schwab/Exterior Problem/error_plot.m 0 → 100644
View file @ 3f21861a
% Error computation andd plotting script
N = 20:20:160;
av_errs = 0*N;
l2_errs = 0*N;
for i = 1:size(N,2)
%[l2_errs(i),av_errs(i)] = sphere_torus_exterior_laplace(N(i));
[l2_errs(i),av_errs(i)] = double_sphere_exterior_laplace(N(i));
end
figure;
loglog(N,l2_errs,'o--');
hold on;
loglog(N,av_errs,'o--');
legend('L2 errors','Av errors');
\ No newline at end of file
Piyush's codes/3D BEM Sauter Schwab/Exterior Problem/single_sphere_exterior_laplace.m 0 → 100644
View file @ 3f21861a
%function l2_err = compute_err_ext_sph(N)
close all; clc; clear;
addpath(genpath("../../"));
for N = 20:10:500
disp(N);
% Initializing parameters for the problem
%N = 200;
% Radius of sphere
Rad = 10;
% Mesh for the geometry
mesh = mshSphere(N,Rad);
N_vtcs = size(mesh.vtx,1);
%mesh.vtx = mesh.vtx + 3*ones(N_vtcs,3);
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(0,0,0);
%% Solving the Laplace problem using BEM
% Definition of FEM spaces and integration rule
S1_Gamma = fem(mesh,'P1');
S0_Gamma = fem(mesh,'P0');
order = 3;
Gamma = dom(mesh,order);
% 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
g = @(x) 1/4/pi./vecnorm(x,2,2); % Point charge at origin
% Checking the boundary condition
figure;
plot(mesh);
hold on;
PI_g = g(S0_Gamma.dof);
plot(mesh,PI_g);
title('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 = V\((0.5 * M + K)* (M\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;
%% Exact solution
Psi_exact = -1/4/pi./vecnorm(dofs,2,2).^3 .* (dot(dofs,nrms,2));
% Visualizing the exactness
figure;
plot(Psi./Psi_exact);
err_vec = Psi-Psi_exact;
% L2 Error
l2_err = norm(err_vec)
av_err = err_vec' * V * err_vec
close all
end
\ No newline at end of file
Piyush's codes/3D BEM Sauter Schwab/Exterior Problem/sphere_torus_exterior_laplace.m 0 → 100644
View file @ 3f21861a
function [l2_err,av_err] = sphere_torus_exterior_laplace(N)
addpath(genpath("../../"));
% Initializing parameters for the problem
%N = 100;
% Radius of the sphere
Rad = 5;
% Radii for the torus
r1 = 5;
r2 = 2;
% Mesh for the geometry
mesh_out = mshSphere(N,Rad);
mesh_in = mshTorus(N,r1,r2);
% Rotate the inner torus
rng(10);
A = rand(3,3);
[Q,R] = qr(A);
mesh_in.vtx = mesh_in.vtx * Q;
% Translate the torus
tx = 15;
ty = 0;
tz = 0;
dist = norm([tx ty tz]);
N_vtcs = size(mesh_in.vtx,1);
trans = ones(N_vtcs,1) * [tx ty tz];
mesh_in.vtx = mesh_in.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(0,0,0);
%% BEM
% Definition of FEM spaces and integration rule
S1_Gamma = fem(mesh,'P1');
S0_Gamma = fem(mesh,'P0');
order = 3;
Gamma = dom(mesh,order);
% 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
g = @(x) 1/4/pi./vecnorm(x,2,2); % Point charge at origin
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));
% 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;
%% Exact solution
Psi_exact = -1/4/pi./vecnorm(dofs,2,2).^3 .* (dot(dofs,nrms,2));
% Visualizing the exactness
figure;
plot(Psi./Psi_exact);
err_vec = Psi-Psi_exact;
% L2 Error
l2_err = norm(err_vec)
av_err = err_vec' * V * err_vec
close all;
end
\ No newline at end of file
Piyush's codes/3D BEM Sauter Schwab/Torque computations/torque_cube_cube.m 0 → 100644
View file @ 3f21861a
%function [] = force_cube_cube_ext(N)
close all; clc; clear;
addpath(genpath("../../"));
global X;
global W;
load('X3','X');
load('W3','W');
% Initializing parameters for the problem
N = 50;
scale = 10;
% Dimensions of the cube
L = 10*[1,1,1];
% Mesh for the geometry
cube_vol_mesh = mshCube(N,L);
mesh_in = cube_vol_mesh.bnd;
cube_vol_mesh = mshCube(N,L);
mesh_out = cube_vol_mesh.bnd;
% Translate the outer cube
tx = scale * 2;%1;
ty = scale * 0;%3;
tz = scale * 0;%2;
dist = norm([tx ty tz]);
N_vtcs = size(mesh_out.vtx,1);
trans = [tx * ones(N_vtcs,1), ty * ones(N_vtcs,1), tz * ones(N_vtcs,1)];
mesh_out.vtx = mesh_out.vtx + trans;
% Join to create the final mesh
mesh = union(mesh_in,mesh_out);
figure;
plot(mesh);
title('Mesh');
%% BEM
% 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)