Commit 0090d7f8 authored by ppanchal's avatar ppanchal
Browse files

BEM evaluation using double layer

parent b229c718
addpath(genpath("../../../../"));
clear;
clc;
format long;
Nvals = 50:30:1000;
sz = size(Nvals,2);
torques_mst = zeros(sz,3);
torques_bem = zeros(sz,3);
forces_mst = zeros(sz,3);
forces_bem = zeros(sz,3);
hvals = zeros(sz,1);
Nux = @(X) (vecnorm(X,2,2)<18).*[X(:,1)==X(:,1), 0 * X(:,1) , 0 * X(:,1)];
Nuy = @(X) (vecnorm(X,2,2)<18).*[0 * X(:,1), 1 * (X(:,1)==X(:,1)), 0 * X(:,1)];
for i = 1:sz
disp(Nvals(i));
% Get the mesh
[mesh,mesh_in,mesh_out] = sph_tor_mesh(10,3,5,Nvals(i),10);
hvals(i) = mean(mesh.ndv,1);
% Solve the floating potential problem on mesh
[Psi,c] = solve_float_pt_ext(mesh,mesh_in,1e2,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 and parallelization
force_bem_gypsi = compute_bem_forces_gypsi(mesh,Psi,Nuy)
end
function torques = compute_bem_forces_gypsi(mesh,R,Psi)
% Function to compute BEM forces using Gypsi implementation.
function out = compute_bem_forces_gypsi(mesh,Psi,nu)
% Force computation using the double layer. One arg is Psi, other is
% Psi \nu \cdot n
S0_Gamma = fem(mesh,'P0');
normals = mesh.nrm;
dofs = S0_Gamma.dof;
Gamma = dom(mesh,3);
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];
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));
% Right vectors for forces
Psi_nu = Psi.* dot(nu(dofs),normals,2);
out = 2*dot(Psi_nu,K*Psi);
end
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment