Commit b229c718 authored by ppanchal's avatar ppanchal
Browse files

Cube torus parallel

parent d351ae93
function [mesh,mesh_in,mesh_out] = cube_tor_mesh(r1,r2,L,N,seed)
rng(seed);
A = rand(3,3);
[Qrot,~] = qr(A);
% Mesh for the geometry
mesh_out = mshTorus(N,r1,r2);
mesh_in_vol = mshCube(N/2,L);
mesh_in = mesh_in_vol.bnd;
% Rotate the outer torus
mesh_out.vtx = mesh_out.vtx * Qrot;
% Translate the outer torus
N_vtcs = size(mesh_out.vtx,1);
trans = ones(N_vtcs,1) * [25 25 25];
mesh_out.vtx = mesh_out.vtx + trans;
% 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);
forces_bemx = zeros(sz,1);
forces_bemy = zeros(sz,1);
forces_bemz = zeros(sz,1);
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');
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
force_bemx = compute_bem_forcex(mesh,18,Psi);
forces_bemx(i) = force_bemx;
fts = zeros(6,1);
parfor it = 1:6
switch it
case 1
force_bemx = compute_bem_forcex(mesh,18,Psi);
fts(it) = force_bemx;
case 2
force_bemy = compute_bem_forcex(mesh,18,Psi);
fts(it) = force_bemx;
case 3
case 4
case 5
case 6
end
end
end
%save('tor_tor_data.mat','Nvals','forces_mst','torques_bem','torques_mst','forces_bem','hvals');
\ No newline at end of file
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);
for i = 1:sz
disp(Nvals(i));
% Get the mesh
[mesh,mesh_in,mesh_out] = cube_tor_mesh(10,3,[10 10 10],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,torque_bem] = compute_bem_ft_par(mesh,12,[0,0,0],Psi);
torques_bem(i,:) = torque_bem;
forces_bem(i,:) = force_bem;
save('cube_tor_data_par.mat','Nvals','forces_mst','torques_bem','torques_mst','forces_bem','hvals');
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