Deathslayer7
2013-10-29, 01:14 PM
Hey all,
I received a code from my professor using a meshless mesh. It's a simple problem but there is a problem. It is getting close to what we need it to do but its not quite there.
http://i198.photobucket.com/albums/aa125/death_slayer7/contour.png
http://i198.photobucket.com/albums/aa125/death_slayer7/contour1.png
The picture shows a 2-D linear temperature plot. The problem being is that the contour plot should show straight lines. Even if I increase the number of points, the problem exists, though not as bad as above.
The Code: Not very long. I cut out some of it that wasn't applicable. I almost want to say it seems like a boundary error somewhere.
% Test meshless method for Steady State Heat Conduction Example
% Clear workspace:
clc
close all
clear all
% Specify random or structured grids:
random_grids = 1; % 1 = random, 2 = structured
length_x = 1.0;
length_y = 1.0;
% Provide user defined values for structured grids (for grids =2):
dx = 0.04; %1/25
dy = 0.04; %1/25
% Specify parameters for random grid generation (for grids = 1):
num_avg_spanwise_grids = 11;
temp_left_wall = 1;
temp_right_wall = 2;
if random_grids == 1
dx_avg = length_x/(num_avg_spanwise_grids-1); %=0.025
dy_avg = length_y/(num_avg_spanwise_grids-1); %=0.025
dx = dx_avg; %=0.025
dy = dy_avg; %=0.025
% Define number of nodes
num_x_nodes = length_x/dx+1; %=41 1x1 matrix
num_y_nodes = length_y/dy+1; %=41 1x1 matrix
total_nodes = num_x_nodes*num_y_nodes; %=1681 1x1 matrix
grids = 1:total_nodes; %1x1681 matrix
%Define 1681x1681 matrix and 1681x1 matrix
Z = zeros(total_nodes,total_nodes); %1681x1681 matrix
Z1 = zeros(total_nodes,1); %1681x1 matrix
% Define shape parameter value based on given grid domain:
%c = dx; % okay 0.025
%c = 0.815*sqrt(dx^2+dy^2); % moderate 0.0281
%c = 1.25*sqrt(dx^2+dy^2)/sqrt(total_nodes); % bad .0010779066786
c = 2/sqrt(total_nodes); % great 0.04878
c=0.1
% Preallocate matrices for speed:
coord_x = Z1; %1681x1 matrix
coord_y = Z1; %1681x1 matrix
K = Z; %1681x1681 matrix
F = Z1; %1681x1 matrix
delxsq = Z; %1681x1681 matrix
delysq = Z; %1681x1681 matrix
delrsq = Z; %1681x1681 matrix
total_shape_op = Z; %1681x1681 matrix
total_shape_deriv_x_op = Z; %1681x1681 matrix
total_shape_deriv_y_op = Z; %1681x1681 matrix
total_shape_laplacian_op = Z; %1681x1681 matrix
% Define all grid sets:
all_points = 1:total_nodes; %1x1681 matrix
bottom_wall = 1:num_x_nodes; %1x41 matrix
top_wall = (num_y_nodes-1)*num_x_nodes+1:total_nodes; %1x41 matrix
left_wall = 1:num_x_nodes:((num_y_nodes-1)*num_x_nodes+1); %1x41 matrix
right_wall = num_x_nodes:num_x_nodes:total_nodes; %1x41 matrix
interior_points = setdiff(all_points,[bottom_wall,top_wall,left_wall,right_wall]); %1x1521
%left_wall=left_wall' %41x1 matrix
%right_wall=right_wall' %41x1 matrix
% Randomize all points:
for L = all_points % 1:total_nodes 1x1681 matrix
coord_x(L) = rand*length_x;
coord_y(L) = rand*length_y;
end
% Specify coordinate restrictions on walls (which will result in
% correct edge and corner coordinates since points overlap properly)
% Walls:
coord_y(bottom_wall) = 0;
coord_y(top_wall) = length_y;
coord_x(left_wall) = 0;
coord_x(right_wall) = length_x;
% Plot grid coordinate locations:
figure(1)
plot(coord_x,coord_y,'bo')
title('Grid coordinate locations')
end
% Create completely sampled shape operator:
for i = 1:total_nodes
for j = 1:total_nodes
delxsq(i,j) = (coord_x(i)-coord_x(j))^2;
delysq(i,j) = (coord_y(i)-coord_y(j))^2;
delrsq(i,j) = delxsq(i,j)+delysq(i,j);
total_shape_op(i,j) = sqrt(delrsq(i,j)+c^2);
total_shape_deriv_x_op(i,j) = (coord_x(i)-coord_x(j))/(delrsq(i,j)+c^2)^(1/2);
total_shape_deriv_y_op(i,j) = (coord_y(i)-coord_y(j))/(delrsq(i,j)+c^2)^(1/2);
total_shape_laplacian_op(i,j) = (delrsq(i,j)+2*c^2)/(delrsq(i,j)+c^2)^(3/2);
end
end
% Apply interior domain conditions:
K(interior_points,:) = total_shape_laplacian_op(interior_points,:);
F(interior_points,1) = 0;
% Apply boundary domain conditions:
% Bottom wall:
K(bottom_wall,:) = total_shape_deriv_y_op(bottom_wall,:);
F(bottom_wall,1) = 0; % Value of derivative on bottom wall
% Top wall:
K(top_wall,:) = total_shape_deriv_y_op(top_wall,:);
F(top_wall,1) = 0; % Value of derivative on top wall
% Left wall:
K(left_wall,:) = total_shape_op(left_wall,:);
F(left_wall,1) = temp_left_wall; % Value of temperature on left wall
% Right wall:
K(right_wall,:) = total_shape_op(right_wall,:);
F(right_wall,1) = temp_right_wall; % Value of temperature on right wall
% Solve system of equations:
U = K\F;
% Derive actual temperature and temperature derivative values:
T = total_shape_op*U;
T_y_deriv = total_shape_deriv_y_op*U;
T_laplacian = total_shape_laplacian_op*U;
if random_grids == 1
[X,Y] = meshgrid([0:(dx_avg):length_y],[0:(dy_avg):length_y]);
T_plot = griddata(coord_x,coord_y,T,X,Y,'v4');
F_2D_interp = TriScatteredInterp(coord_x, coord_y, T);
T_plot = F_2D_interp(X,Y);
end
% Plot results and boundary conditions to test accuracy
figure(11)
contourf(X,Y,T_plot)
title('Contour plot of Temperature')
hold on
plot(coord_x,coord_y,'wo')
colorbar
hold off
figure(21)
surf(X,Y,T_plot)
title('Surface plot of Temperature')
I received a code from my professor using a meshless mesh. It's a simple problem but there is a problem. It is getting close to what we need it to do but its not quite there.
http://i198.photobucket.com/albums/aa125/death_slayer7/contour.png
http://i198.photobucket.com/albums/aa125/death_slayer7/contour1.png
The picture shows a 2-D linear temperature plot. The problem being is that the contour plot should show straight lines. Even if I increase the number of points, the problem exists, though not as bad as above.
The Code: Not very long. I cut out some of it that wasn't applicable. I almost want to say it seems like a boundary error somewhere.
% Test meshless method for Steady State Heat Conduction Example
% Clear workspace:
clc
close all
clear all
% Specify random or structured grids:
random_grids = 1; % 1 = random, 2 = structured
length_x = 1.0;
length_y = 1.0;
% Provide user defined values for structured grids (for grids =2):
dx = 0.04; %1/25
dy = 0.04; %1/25
% Specify parameters for random grid generation (for grids = 1):
num_avg_spanwise_grids = 11;
temp_left_wall = 1;
temp_right_wall = 2;
if random_grids == 1
dx_avg = length_x/(num_avg_spanwise_grids-1); %=0.025
dy_avg = length_y/(num_avg_spanwise_grids-1); %=0.025
dx = dx_avg; %=0.025
dy = dy_avg; %=0.025
% Define number of nodes
num_x_nodes = length_x/dx+1; %=41 1x1 matrix
num_y_nodes = length_y/dy+1; %=41 1x1 matrix
total_nodes = num_x_nodes*num_y_nodes; %=1681 1x1 matrix
grids = 1:total_nodes; %1x1681 matrix
%Define 1681x1681 matrix and 1681x1 matrix
Z = zeros(total_nodes,total_nodes); %1681x1681 matrix
Z1 = zeros(total_nodes,1); %1681x1 matrix
% Define shape parameter value based on given grid domain:
%c = dx; % okay 0.025
%c = 0.815*sqrt(dx^2+dy^2); % moderate 0.0281
%c = 1.25*sqrt(dx^2+dy^2)/sqrt(total_nodes); % bad .0010779066786
c = 2/sqrt(total_nodes); % great 0.04878
c=0.1
% Preallocate matrices for speed:
coord_x = Z1; %1681x1 matrix
coord_y = Z1; %1681x1 matrix
K = Z; %1681x1681 matrix
F = Z1; %1681x1 matrix
delxsq = Z; %1681x1681 matrix
delysq = Z; %1681x1681 matrix
delrsq = Z; %1681x1681 matrix
total_shape_op = Z; %1681x1681 matrix
total_shape_deriv_x_op = Z; %1681x1681 matrix
total_shape_deriv_y_op = Z; %1681x1681 matrix
total_shape_laplacian_op = Z; %1681x1681 matrix
% Define all grid sets:
all_points = 1:total_nodes; %1x1681 matrix
bottom_wall = 1:num_x_nodes; %1x41 matrix
top_wall = (num_y_nodes-1)*num_x_nodes+1:total_nodes; %1x41 matrix
left_wall = 1:num_x_nodes:((num_y_nodes-1)*num_x_nodes+1); %1x41 matrix
right_wall = num_x_nodes:num_x_nodes:total_nodes; %1x41 matrix
interior_points = setdiff(all_points,[bottom_wall,top_wall,left_wall,right_wall]); %1x1521
%left_wall=left_wall' %41x1 matrix
%right_wall=right_wall' %41x1 matrix
% Randomize all points:
for L = all_points % 1:total_nodes 1x1681 matrix
coord_x(L) = rand*length_x;
coord_y(L) = rand*length_y;
end
% Specify coordinate restrictions on walls (which will result in
% correct edge and corner coordinates since points overlap properly)
% Walls:
coord_y(bottom_wall) = 0;
coord_y(top_wall) = length_y;
coord_x(left_wall) = 0;
coord_x(right_wall) = length_x;
% Plot grid coordinate locations:
figure(1)
plot(coord_x,coord_y,'bo')
title('Grid coordinate locations')
end
% Create completely sampled shape operator:
for i = 1:total_nodes
for j = 1:total_nodes
delxsq(i,j) = (coord_x(i)-coord_x(j))^2;
delysq(i,j) = (coord_y(i)-coord_y(j))^2;
delrsq(i,j) = delxsq(i,j)+delysq(i,j);
total_shape_op(i,j) = sqrt(delrsq(i,j)+c^2);
total_shape_deriv_x_op(i,j) = (coord_x(i)-coord_x(j))/(delrsq(i,j)+c^2)^(1/2);
total_shape_deriv_y_op(i,j) = (coord_y(i)-coord_y(j))/(delrsq(i,j)+c^2)^(1/2);
total_shape_laplacian_op(i,j) = (delrsq(i,j)+2*c^2)/(delrsq(i,j)+c^2)^(3/2);
end
end
% Apply interior domain conditions:
K(interior_points,:) = total_shape_laplacian_op(interior_points,:);
F(interior_points,1) = 0;
% Apply boundary domain conditions:
% Bottom wall:
K(bottom_wall,:) = total_shape_deriv_y_op(bottom_wall,:);
F(bottom_wall,1) = 0; % Value of derivative on bottom wall
% Top wall:
K(top_wall,:) = total_shape_deriv_y_op(top_wall,:);
F(top_wall,1) = 0; % Value of derivative on top wall
% Left wall:
K(left_wall,:) = total_shape_op(left_wall,:);
F(left_wall,1) = temp_left_wall; % Value of temperature on left wall
% Right wall:
K(right_wall,:) = total_shape_op(right_wall,:);
F(right_wall,1) = temp_right_wall; % Value of temperature on right wall
% Solve system of equations:
U = K\F;
% Derive actual temperature and temperature derivative values:
T = total_shape_op*U;
T_y_deriv = total_shape_deriv_y_op*U;
T_laplacian = total_shape_laplacian_op*U;
if random_grids == 1
[X,Y] = meshgrid([0:(dx_avg):length_y],[0:(dy_avg):length_y]);
T_plot = griddata(coord_x,coord_y,T,X,Y,'v4');
F_2D_interp = TriScatteredInterp(coord_x, coord_y, T);
T_plot = F_2D_interp(X,Y);
end
% Plot results and boundary conditions to test accuracy
figure(11)
contourf(X,Y,T_plot)
title('Contour plot of Temperature')
hold on
plot(coord_x,coord_y,'wo')
colorbar
hold off
figure(21)
surf(X,Y,T_plot)
title('Surface plot of Temperature')