double_sphere_exterior_laplace.m 2.89 KB
 ppanchal committed Jan 12, 2022 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 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 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 ``````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``````