Commit 44e48fcc authored by ppanchal's avatar ppanchal
Browse files

Parallel BEM computation

parent f0db51a0
i hval force(bem) force(MST) torque(BEM) torque (MST)
TEST messagei hval force(bem) force(MST) torque(BEM) torque (MST)01 000 -3.916580e-06-2.553163e-06-1.504464e-05 000 -5.608747e-061.607534e-05-1.960747e-052 6.814332e+00 000 5.706542e-07-2.998876e-06-1.211921e-05 000 -3.889133e-062.078131e-05-1.971177e-053 4.904223e+00 000 3.244160e-06-2.947182e-06-1.004992e-05 000 -2.559172e-062.442521e-05-2.033959e-054 3.858383e+00 000 4.762897e-06-2.914208e-06-9.123821e-06 000 -2.034297e-062.658571e-05-2.076515e-055 3.111980e+00 000 5.975181e-06-2.859049e-06-8.373631e-06 000 -1.636651e-062.827548e-05-2.115447e-056 2.866442e+00 000 6.815184e-06-2.793611e-06-7.737742e-06 000 -1.324539e-062.917643e-05-2.139922e-057 2.374319e+00 000 7.684049e-06-2.721293e-06-7.219769e-06 000 -1.086611e-063.040642e-05-2.173226e-058 2.212235e+00 000 8.302401e-06-2.649556e-06-6.767976e-06 000 -8.833717e-073.107261e-05-2.194173e-059 1.932759e+00 000 8.692780e-06-2.612316e-06-6.575881e-06 000 -8.047341e-073.170960e-05-2.212998e-0510 1.812965e+00 000 8.948203e-06-2.445856e-06-6.086192e-06 000 -9.679802e-072.983890e-05-2.110386e-0511 1.758448e+00 000 9.427427e-06-2.508557e-06-6.051231e-06 000 -5.875085e-073.250222e-05-2.240155e-0512 1.512210e+00 000 9.660693e-06-2.261545e-06-5.483985e-06 000 -9.649860e-072.954665e-05-2.073041e-0513 1.350880e+00 000 1.019104e-05-2.409082e-06-5.632338e-06 000 -4.387030e-073.361628e-05-2.276973e-0514 1.280197e+00 000 1.124030e-05-1.927124e-06-4.502836e-06 000 -9.304134e-073.173117e-05-2.137815e-0515 1.216525e+00 000 1.065969e-05-2.054842e-06-4.776711e-06 000 -7.867103e-073.024352e-05-2.084127e-0516 1.124854e+00 000 1.091814e-05-2.289899e-06-5.169507e-06 000 -2.760119e-073.449583e-05-2.309797e-0517 1.071545e+00 000 1.193030e-05-1.913061e-06-4.279273e-06 000 -5.353457e-073.406248e-05-2.253562e-0518 1.023050e+00 000 1.216031e-05-1.873271e-06-4.133686e-06 000 -4.459868e-073.436970e-05-2.268951e-0519 9.512547e-01 000 1.134262e-05-2.005473e-06-4.506044e-06 000 -5.226016e-073.205561e-05-2.165164e-0520 9.302025e-01 000 1.163570e-05-2.158892e-06-4.713133e-06 000 -1.307315e-073.534311e-05-2.342965e-0521 8.907701e-01 000 1.186587e-05-2.111805e-06-4.558492e-06 000 -8.199200e-083.558729e-05-2.353272e-0522 8.322881e-01 000 1.189978e-05-2.111206e-06-4.560676e-06 000 -8.845643e-083.568993e-05-2.356364e-0523 8.150113e-01 000 1.200728e-05-2.088358e-06-4.487976e-06 000 -6.467242e-083.580308e-05-2.361471e-0524 7.825204e-01 000 1.221074e-05-2.045428e-06-4.353182e-06 000 -2.555800e-083.601847e-05-2.370643e-0525 7.342644e-01 000 1.223856e-05-2.044878e-06-4.354861e-06 000 -2.987589e-083.610235e-05-2.373330e-0526 7.199130e-01 000 1.277026e-05-1.785069e-06-3.812307e-06 000 -2.360877e-073.522443e-05-2.308501e-0527 7.061111e-01 000 1.242608e-05-2.004085e-06-4.230650e-06 000 5.998191e-093.630017e-05-2.382111e-0528 6.800340e-01 000 1.259961e-05-1.965401e-06-4.116006e-06 000 3.861402e-083.648209e-05-2.390230e-0529 6.405029e-01 000 1.262281e-05-1.964893e-06-4.117302e-06 000 3.488255e-083.655217e-05-2.392562e-0530 6.288893e-01 000 1.331347e-05-1.967638e-06-4.020968e-06 000 2.873966e-073.997664e-05-2.572036e-0531 6.176889e-01 000 1.269668e-05-1.899011e-06-3.971648e-06 000 -3.749592e-093.591958e-05-2.355898e-0532 5.964423e-01 000 1.351139e-05-1.912526e-06-3.873351e-06 000 3.085090e-074.003976e-05-2.575312e-05
\ No newline at end of file
function torques = compute_bem_forces_gypsi(mesh,R,Psi)
S0_Gamma = fem(mesh,'P0');
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];
end
\ No newline at end of file
function forcex = compute_bem_forcex(mesh,R,Psi)
S0_Gamma = fem(mesh,'P0');
Nux = @(X) (vecnorm(X,2,2)<R).* ones(size(X,1),1)*[1 0 0];
% Kernels
kernelx = @(x,y,z) sum(z.*(Nux(x) - Nux(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);
end
\ No newline at end of file
function forcey = compute_bem_forcey(mesh,R,Psi)
S0_Gamma = fem(mesh,'P0');
Nuy = @(X) (vecnorm(X,2,2)<R).* ones(size(X,1),1)*[0 1 0];
% Kernels
kernely = @(x,y,z) sum(z.*(Nuy(x) - Nuy(y)), 2)./(vecnorm(z,2,2).^3)/ (4*pi);
t2maty = panel_oriented_assembly(mesh,kernely,S0_Gamma,S0_Gamma);
forcey = 0.5 * dot(Psi,t2maty*Psi);
end
\ No newline at end of file
function forcez = compute_bem_forcez(mesh,R,Psi)
S0_Gamma = fem(mesh,'P0');
Nuz = @(X) (vecnorm(X,2,2)<R).* ones(size(X,1),1)*[0 0 1];
% Kernels
kernelz = @(x,y,z) sum(z.*(Nuz(x) - Nuz(y)), 2)./(vecnorm(z,2,2).^3)/ (4*pi);
t2matz = panel_oriented_assembly(mesh,kernelz,S0_Gamma,S0_Gamma);
forcez = 0.5 * dot(Psi,t2matz*Psi);
end
\ No newline at end of file
function [fbem,tbem] = compute_bem_ft_par(mesh,R,Xcg,Psi)
S0_Gamma = fem(mesh,'P0');
% Translation fields
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];
% Rotational fields
Nuxr = @(X) (vecnorm(X,2,2)<R).* cross(ones(size(X,1),1)*[1 0 0],X-Xcg);
Nuyr = @(X) (vecnorm(X,2,2)<R).* cross(ones(size(X,1),1)*[0 1 0],X-Xcg);
Nuzr = @(X) (vecnorm(X,2,2)<R).* cross(ones(size(X,1),1)*[0 0 1],X-Xcg);
% Kernels translation
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);
% Kernels rotation
kernelxr = @(x,y,z) sum(z.*(Nuxr(x) - Nuxr(y)), 2)./(vecnorm(z,2,2).^3)/ (4*pi);
kernelyr = @(x,y,z) sum(z.*(Nuyr(x) - Nuyr(y)), 2)./(vecnorm(z,2,2).^3)/ (4*pi);
kernelzr = @(x,y,z) sum(z.*(Nuzr(x) - Nuzr(y)), 2)./(vecnorm(z,2,2).^3)/ (4*pi);
kernels = cell(6,1);
kernels{1} = kernelx;
kernels{2} = kernely;
kernels{3} = kernelz;
kernels{4} = kernelxr;
kernels{5} = kernelyr;
kernels{6} = kernelzr;
ft = cell(6,1);
parfor i = 1:6
t2mat = panel_oriented_assembly(mesh,kernels{i},S0_Gamma,S0_Gamma);
ft{i} = 0.5 * dot(Psi,t2mat*Psi);
end
fbem = [ft{1} ft{2} ft{3}];
tbem = [ft{4} ft{5} ft{6}];
end
\ No newline at end of file
function torquex = compute_bem_torquex(mesh,R,Xcg,Psi)
S0_Gamma = fem(mesh,'P0');
Nux = @(X) (vecnorm(X,2,2)<R).* cross(ones(size(X,1),1)*[1 0 0],X-Xcg);
% Kernels
kernelx = @(x,y,z) sum(z.*(Nux(x) - Nux(y)), 2)./(vecnorm(z,2,2).^3)/ (4*pi);
t2matx = panel_oriented_assembly(mesh,kernelx,S0_Gamma,S0_Gamma);
torquex = 0.5 * dot(Psi,t2matx*Psi);
end
\ No newline at end of file
function torquey = compute_bem_torquey(mesh,R,Xcg,Psi)
S0_Gamma = fem(mesh,'P0');
Nuy = @(X) (vecnorm(X,2,2)<R).* cross(ones(size(X,1),1)*[0 1 0],X-Xcg);
% Kernels
kernely = @(x,y,z) sum(z.*(Nuy(x) - Nuy(y)), 2)./(vecnorm(z,2,2).^3)/ (4*pi);
t2maty = panel_oriented_assembly(mesh,kernely,S0_Gamma,S0_Gamma);
torquey = 0.5 * dot(Psi,t2maty*Psi);
end
\ No newline at end of file
function torquez = compute_bem_torquez(mesh,R,Xcg,Psi)
S0_Gamma = fem(mesh,'P0');
Nuz = @(X) (vecnorm(X,2,2)<R).* cross(ones(size(X,1),1)*[0 0 1],X-Xcg);
% Kernels
kernelz = @(x,y,z) sum(z.*(Nuz(x) - Nuz(y)), 2)./(vecnorm(z,2,2).^3)/ (4*pi);
t2matz = panel_oriented_assembly(mesh,kernelz,S0_Gamma,S0_Gamma);
torquez = 0.5 * dot(Psi,t2matz*Psi);
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;
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] = 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 and parallelization
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_forcey(mesh,18,Psi);
fts(it) = force_bemy;
case 3
force_bemz = compute_bem_forcez(mesh,18,Psi);
fts(it) = force_bemz;
case 4
torque_bemx = compute_bem_torquex(mesh,18,[0 0 0],Psi);
fts(it) = torque_bemx;
case 5
torque_bemy = compute_bem_torquey(mesh,18,[0 0 0],Psi);
fts(it) = torque_bemy;
case 6
torque_bemz = compute_bem_torquez(mesh,18,[0 0 0],Psi);
fts(it) = torque_bemz;
end
end
forces_bem(i,:) = fts(1:3,1)';
torques_bem(i,:) = fts(4:6,1)';
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;
global X;
global W;
load('X3','X');
load('W3','W');
fname = "abc.txt";
out
fprintf("i hval force(bem) force(MST) torque(BEM) torque (MST)")
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');
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');
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] = 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 and parallelization
[force_bem,torque_bem] = compute_bem_ft_par(mesh,18,[0,0,0],Psi);
torques_bem(i,:) = torque_bem;
forces_bem(i,:) = force_bem;
save('tor_tor_data_par.mat','Nvals','forces_mst','torques_bem','torques_mst','forces_bem','hvals');
end
\ No newline at end of file
% Script to print the matrices in right format
% fprintf("[");
% for i = 1:81
% fprintf("%.17d,",W(i));
% end
% fprintf("\n");
fprintf("[");
for i = 1:81
for j = 1:4
if j == 4
fprintf("%.17d;",X(i,j));
else
fprintf("%.17d,",X(i,j));
end
end
end
fprintf("]\n");
\ No newline at end of file
% Function to integrate on the 4 dimensional cube, hard coded n = 5
function integral = integrate4dim(f)
%load('X','X');
%load('W','W');
global X;
global W;
% w = [0.568888888888889;
% 0.478628670499367;
% 0.478628670499367;
% 0.236926885056189;
% 0.236926885056189];
% x = [0;
% -0.538469310105683;
% 0.538469310105683;
% -0.906179845938664;
% 0.906179845938664];
% % scaling for integral from 0 to 1
% w = w * 0.5;
% x = (x+1)/2;
% integral2 = 0;
integral = dot(W,f(X(:,1),X(:,2),X(:,3),X(:,4)));
% for i = 1:5
% for j = 1:5
% for k = 1:5
% for l = 1:5
% integral2 = integral2 + w(i)*w(j)*w(k)*w(l)*f(x(i),x(j),x(k),x(l));
% end
% end
% end
% end
% disp(integral-integral2);
% Function to integrate on the 4 dimensional cube, hard coded n = 5
function integral = integrate4dim(f)
%load('X','X');
%load('W','W');
%global X;
%global W;
W = [5.95374180765122045e-03,9.52598689224197319e-03,5.95374180765122045e-03,9.52598689224197319e-03,1.52415790275871894e-02,9.52598689224197319e-03,5.95374180765122045e-03,9.52598689224197319e-03,5.95374180765122045e-03,9.52598689224197319e-03,1.52415790275871894e-02,9.52598689224197319e-03,1.52415790275871911e-02,2.43865264441395571e-02,1.52415790275871911e-02,9.52598689224197319e-03,1.52415790275871894e-02,9.52598689224197319e-03,5.95374180765122045e-03,9.52598689224197319e-03,5.95374180765122045e-03,9.52598689224197319e-03,1.52415790275871894e-02,9.52598689224197319e-03,5.95374180765122045e-03,9.52598689224197319e-03,5.95374180765122045e-03,9.52598689224197319e-03,1.52415790275871894e-02,9.52598689224197319e-03,1.52415790275871911e-02,2.43865264441395571e-02,1.52415790275871911e-02,9.52598689224197319e-03,1.52415790275871894e-02,9.52598689224197319e-03,1.52415790275871911e-02,2.43865264441395571e-02,1.52415790275871911e-02,2.43865264441395571e-02,3.90184423106233746e-02,2.43865264441395571e-02,1.52415790275871911e-02,2.43865264441395571e-02,1.52415790275871911e-02,9.52598689224197319e-03,1.52415790275871894e-02,9.52598689224197319e-03,1.52415790275871911e-02,2.43865264441395571e-02,1.52415790275871911e-02,9.52598689224197319e-03,1.52415790275871894e-02,9.52598689224197319e-03,5.95374180765122045e-03,9.52598689224197319e-03,5.95374180765122045e-03,9.52598689224197319e-03,1.52415790275871894e-02,9.52598689224197319e-03,5.95374180765122045e-03,9.52598689224197319e-03,5.95374180765122045e-03,9.52598689224197319e-03,1.52415790275871894e-02,9.52598689224197319e-03,1.52415790275871911e-02,2.43865264441395571e-02,1.52415790275871911e-02,9.52598689224197319e-03,1.52415790275871894e-02,9.52598689224197319e-03,5.95374180765122045e-03,9.52598689224197319e-03,5.95374180765122045e-03,9.52598689224197319e-03,1.52415790275871894e-02,9.52598689224197319e-03,5.95374180765122045e-03,9.52598689224197319e-03,5.95374180765122045e-03]';
X = [1.12701665379258242e-01,1.12701665379258242e-01,1.12701665379258242e-01,1.12701665379258242e-01;1.12701665379258242e-01,1.12701665379258242e-01,1.12701665379258242e-01,5.00000000000000000e-01;1.12701665379258242e-01,1.12701665379258242e-01,1.12701665379258242e-01,8.87298334620741702e-01;1.12701665379258242e-01,1.12701665379258242e-01,5.00000000000000000e-01,1.12701665379258242e-01;1.12701665379258242e-01,1.12701665379258242e-01,5.00000000000000000e-01,5.00000000000000000e-01;1.12701665379258242e-01,1.12701665379258242e-01,5.00000000000000000e-01,8.87298334620741702e-01;1.12701665379258242e-01,1.12701665379258242e-01,8.87298334620741702e-01,1.12701665379258242e-01;1.12701665379258242e-01,1.12701665379258242e-01,8.87298334620741702e-01,5.00000000000000000e-01;1.12701665379258242e-01,1.12701665379258242e-01,8.87298334620741702e-01,8.87298334620741702e-01;1.12701665379258242e-01,5.00000000000000000e-01,1.12701665379258242e-01,1.12701665379258242e-01;1.12701665379258242e-01,5.00000000000000000e-01,1.12701665379258242e-01,5.00000000000000000e-01;1.12701665379258242e-01,5.00000000000000000e-01,1.12701665379258242e-01,8.87298334620741702e-01;1.12701665379258242e-01,5.00000000000000000e-01,5.00000000000000000e-01,1.12701665379258242e-01;1.12701665379258242e-01,5.00000000000000000e-01,5.00000000000000000e-01,5.00000000000000000e-01;1.12701665379258242e-01,5.00000000000000000e-01,5.00000000000000000e-01,8.87298334620741702e-01;1.12701665379258242e-01,5.00000000000000000e-01,8.87298334620741702e-01,1.12701665379258242e-01;1.12701665379258242e-01,5.00000000000000000e-01,8.87298334620741702e-01,5.00000000000000000e-01;1.12701665379258242e-01,5.00000000000000000e-01,8.87298334620741702e-01,8.87298334620741702e-01;1.12701665379258242e-01,8.87298334620741702e-01,1.12701665379258242e-01,1.12701665379258242e-01;1.12701665379258242e-01,8.87298334620741702e-01,1.12701665379258242e-01,5.00000000000000000e-01;1.12701665379258242e-01,8.87298334620741702e-01,1.12701665379258242e-01,8.87298334620741702e-01;1.12701665379258242e-01,8.87298334620741702e-01,5.00000000000000000e-01,1.12701665379258242e-01;1.12701665379258242e-01,8.87298334620741702e-01,5.00000000000000000e-01,5.00000000000000000e-01;1.12701665379258242e-01,8.87298334620741702e-01,5.00000000000000000e-01,8.87298334620741702e-01;1.12701665379258242e-01,8.87298334620741702e-01,8.87298334620741702e-01,1.12701665379258242e-01;1.12701665379258242e-01,8.87298334620741702e-01,8.87298334620741702e-01,5.00000000000000000e-01;1.12701665379258242e-01,8.87298334620741702e-01,8.87298334620741702e-01,8.87298334620741702e-01;5.00000000000000000e-01,1.12701665379258242e-01,1.12701665379258242e-01,1.12701665379258242e-01;5.00000000000000000e-01,1.12701665379258242e-01,1.12701665379258242e-01,5.00000000000000000e-01;5.00000000000000000e-01,1.12701665379258242e-01,1.12701665379258242e-01,8.87298334620741702e-01;5.00000000000000000e-01,1.12701665379258242e-01,5.00000000000000000e-01,1.12701665379258242e-01;5.00000000000000000e-01,1.12701665379258242e-01,5.00000000000000000e-01,5.00000000000000000e-01;5.00000000000000000e-01,1.12701665379258242e-01,5.00000000000000000e-01,8.87298334620741702e-01;5.00000000000000000e-01,1.12701665379258242e-01,8.87298334620741702e-01,1.12701665379258242e-01;5.00000000000000000e-01,1.12701665379258242e-01,8.87298334620741702e-01,5.00000000000000000e-01;5.00000000000000000e-01,1.12701665379258242e-01,8.87298334620741702e-01,8.87298334620741702e-01;5.00000000000000000e-01,5.00000000000000000e-01,1.12701665379258242e-01,1.12701665379258242e-01;5.00000000000000000e-01,5.00000000000000000e-01,1.12701665379258242e-01,5.00000000000000000e-01;5.00000000000000000e-01,5.00000000000000000e-01,1.12701665379258242e-01,8.87298334620741702e-01;5.00000000000000000e-01,5.00000000000000000e-01,5.00000000000000000e-01,1.12701665379258242e-01;5.00000000000000000e-01,5.00000000000000000e-01,5.00000000000000000e-01,5.00000000000000000e-01;5.00000000000000000e-01,5.00000000000000000e-01,5.00000000000000000e-01,8.87298334620741702e-01;5.00000000000000000e-01,5.00000000000000000e-01,8.87298334620741702e-01,1.12701665379258242e-01;5.00000000000000000e-01,5.00000000000000000e-01,8.87298334620741702e-01,5.00000000000000000e-01;5.00000000000000000e-01,5.00000000000000000e-01,8.87298334620741702e-01,8.87298334620741702e-01;5.00000000000000000e-01,8.87298334620741702e-01,1.12701665379258242e-01,1.12701665379258242e-01;5.00000000000000000e-01,8.87298334620741702e-01,1.12701665379258242e-01,5.00000000000000000e-01;5.00000000000000000e-01,8.87298334620741702e-01,1.12701665379258242e-01,8.87298334620741702e-01;5.00000000000000000e-01,8.87298334620741702e-01,5.00000000000000000e-01,1.12701665379258242e-01;5.00000000000000000e-01,8.87298334620741702e-01,5.00000000000000000e-01,5.00000000000000000e-01;5.00000000000000000e-01,8.87298334620741702e-01,5.00000000000000000e-01,8.87298334620741702e-01;5.00000000000000000e-01,8.87298334620741702e-01,8.87298334620741702e-01,1.12701665379258242e-01;5.00000000000000000e-01,8.87298334620741702e-01,8.87298334620741702e-01,5.00000000000000000e-01;5.00000000000000000e-01,8.87298334620741702e-01,8.87298334620741702e-01,8.87298334620741702e-01;8.87298334620741702e-01,1.12701665379258242e-01,1.12701665379258242e-01,1.12701665379258242e-01;8.87298334620741702e-01,1.12701665379258242e-01,1.12701665379258242e-01,5.00000000000000000e-01;8.87298334620741702e-01,1.12701665379258242e-01,1.12701665379258242e-01,8.87298334620741702e-01;8.87298334620741702e-01,1.12701665379258242e-01,5.00000000000000000e-01,1.12701665379258242e-01;8.87298334620741702e-01,1.12701665379258242e-01,5.00000000000000000e-01,5.00000000000000000e-01;8.87298334620741702e-01,1.12701665379258242e-01,5.00000000000000000e-01,8.87298334620741702e-01;8.87298334620741702e-01,1.12701665379258242e-01,8.87298334620741702e-01,1.12701665379258242e-01;8.87298334620741702e-01,1.12701665379258242e-01,8.87298334620741702e-01,5.00000000000000000e-01;8.87298334620741702e-01,1.12701665379258242e-01,8.87298334620741702e-01,8.87298334620741702e-01;8.87298334620741702e-01,5.00000000000000000e-01,1.12701665379258242e-01,1.12701665379258242e-01;8.87298334620741702e-01,5.00000000000000000e-01,1.12701665379258242e-01,5.00000000000000000e-01;8.87298334620741702e-01,5.00000000000000000e-01,1.12701665379258242e-01,8.87298334620741702e-01;8.87298334620741702e-01,5.00000000000000000e-01,5.00000000000000000e-01,1.12701665379258242e-01;8.87298334620741702e-01,5.00000000000000000e-01,5.00000000000000000e-01,5.00000000000000000e-01;8.87298334620741702e-01,5.00000000000000000e-01,5.00000000000000000e-01,8.87298334620741702e-01;8.87298334620741702e-01,5.00000000000000000e-01,8.87298334620741702e-01,1.12701665379258242e-01;8.87298334620741702e-01,5.00000000000000000e-01,8.87298334620741702e-01,5.00000000000000000e-01;8.87298334620741702e-01,5.00000000000000000e-01,8.87298334620741702e-01,8.87298334620741702e-01;8.87298334620741702e-01,8.87298334620741702e-01,1.12701665379258242e-01,1.12701665379258242e-01;8.87298334620741702e-01,8.87298334620741702e-01,1.12701665379258242e-01,5.00000000000000000e-01;8.87298334620741702e-01,8.87298334620741702e-01,1.12701665379258242e-01,8.87298334620741702e-01;8.87298334620741702e-01,8.87298334620741702e-01,5.00000000000000000e-01,1.12701665379258242e-01;8.87298334620741702e-01,8.87298334620741702e-01,5.00000000000000000e-01,5.00000000000000000e-01;8.87298334620741702e-01,8.87298334620741702e-01,5.00000000000000000e-01,8.87298334620741702e-01;8.87298334620741702e-01,8.87298334620741702e-01,8.87298334620741702e-01,1.12701665379258242e-01;8.87298334620741702e-01,8.87298334620741702e-01,8.87298334620741702e-01,5.00000000000000000e-01;8.87298334620741702e-01,8.87298334620741702e-01,8.87298334620741702e-01,8.87298334620741702e-01];
% w = [0.568888888888889;
% 0.478628670499367;
% 0.478628670499367;
% 0.236926885056189;
% 0.236926885056189];
% x = [0;
% -0.538469310105683;
% 0.538469310105683;
% -0.906179845938664;
% 0.906179845938664];
% % scaling for integral from 0 to 1
% w = w * 0.5;
% x = (x+1)/2;
% integral2 = 0;
integral = dot(W,f(X(:,1),X(:,2),X(:,3),X(:,4)));
% for i = 1:5
% for j = 1:5
% for k = 1:5
% for l = 1:5
% integral2 = integral2 + w(i)*w(j)*w(k)*w(l)*f(x(i),x(j),x(k),x(l));
% end
% end
% end
% end
% disp(integral-integral2);
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