Commit aa7fc941 authored by ppanchal's avatar ppanchal
Browse files

smaller steps convergence

parent ea9adf26
loglog(Nvals,abs(torque_mst(:,1)))
loglog(Nvals,abs(torques_mst(:,1)))
hold on
loglog(Nvals,abs(torque_mst(:,2)))
loglog(Nvals,abs(torque_mst(:,3)))
loglog(Nvals,abs(torques_mst(:,2)))
loglog(Nvals,abs(torques_mst(:,3)))
title('Torques vs Nvals');
figure
loglog(Nvals,abs(force_mst(:,1)))
loglog(Nvals,abs(forces_mst(:,1)))
hold on
loglog(Nvals,abs(force_mst(:,2)))
loglog(Nvals,abs(force_mst(:,3)))
\ No newline at end of file
loglog(Nvals,abs(forces_mst(:,2)))
loglog(Nvals,abs(forces_mst(:,3)))
title('Forces vs Nvals');
figure
loglog(hvals,abs(torques_mst(:,1)))
hold on
loglog(hvals,abs(torques_mst(:,2)))
loglog(hvals,abs(torques_mst(:,3)))
title('Torques vs hvals');
figure
loglog(hvals,abs(forces_mst(:,1)))
hold on
loglog(hvals,abs(forces_mst(:,2)))
loglog(hvals,abs(forces_mst(:,3)))
title('Forces vs hhvals');
......@@ -7,18 +7,20 @@ global W;
load('X3','X');
load('W3','W');
Nvals = 50:200:2000;
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);
for i = 1:sz
disp(Nvals(i));
% Get the mesh
[mesh,mesh_in,mesh_out] = tor_tor_mesh(10,3,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,1,3,'gypsi');
......@@ -41,4 +43,4 @@ for i = 1:sz
end
save('new_script_data.mat','Nvals','forces_mst','torques_bem','torques_mst','forces_bem');
save('new_script_data.mat','Nvals','forces_mst','torques_bem','torques_mst','forces_bem','hvals');
function [mesh,mesh_in,mesh_out] = sph_tor_mesh(r1,r2,R,N,seed)
rng(seed);
A = rand(3,3);
[Qrot,~] = qr(A);
% distance between the centers
dist = 3*(r1+r2+R);
% Mesh for the geometry
mesh_in = mshTorus(N,r1,r2);
mesh_out = mshSphere(N,R);
% Translate the outer torus
N_vtcs = size(mesh_out.vtx,1);
trans = ones(N_vtcs,1) * [dist 0 0];
mesh_out.vtx = mesh_out.vtx + trans;
% Rotate the inner torus
mesh_in.vtx = mesh_in.vtx * Qrot;
% Join to create the final mesh
mesh = union(mesh_in,mesh_out);
end
\ No newline at end of file
addpath(genpath("../../../../"));
clear;
clc;
format long;
global X;
global W;
load('X3','X');
load('W3','W');
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);
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,1e5,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
% torque_bem = compute_bem_torques(mesh,18,[0,0,0],Psi);
% torques_bem(i,:) = torque_bem;
%
% % Computing the forces using BEM formula
% force_bem = compute_bem_forces(mesh,18,Psi);
% forces_bem(i,:) = force_bem;
end
save('new_script_data.mat','Nvals','forces_mst','torques_bem','torques_mst','forces_bem','hvals');
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