Error using horzcat Dimensions of arrays being concatenated are not consistent. Error in MKstokes (line 52) StokesMatrix = [A, B_transpose; B, sparse_M_M];
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
% Define problem parameters
mu = 1.0; % viscosity
f = [0; 0]; % forcing term for velocity
g = 0; % forcing term for pressure
%% Mesh creation
h = 0.25; % discretization step
[node,elem] = squaremesh([-1,1,-1,1],h);
subplot(1,2,1); showmesh(node,elem);
N = size(node,1);
NTc = size(elem,1);
% Uniform refinement to get a fine grid
[node,elem] = uniformrefine(node,elem);
% Uniform refinement to get a fine grid
subplot(1,2,2); showmesh(node,elem);
%%extract boundary conditions
[Dirichlet, bdFlag] = ExtractBoundaryEdges(elem, node);
%% Load Stokes Data Partial Differential Equations
%Different velocity and pressure equation models
pde = StokesModelMustafaSin;%StokesModelMustafa;%StokesModelMustafaNew;%Stokesdata3;%StokesModelMustafa2;%StokesModelMustafaSin;%;
option.elemType = 'isoP2P0';
option.fquadorder = 3;
tolerance = 1e-6;
max_iterations = 1000;
mesh.node = node;
mesh.elem = elem;
mesh.bdFlag = bdFlag;
%The following line of code is only run to get the eqn (equation structure) easier without defining it every time
[soln,eqn,info] = StokesisoP2P0(node,elem,bdFlag,pde, option);
% Define matrix A
A = [1 3 0 2; 0 2 1 0; 2 3 7 0; 9 1 0 0]
% Define matrix B
B = [0 2; 1 0; 4 7; 9 0]
% Assemble the finite element system
%[A, B, ~, ~] = assempde(pde,node,elem,'f',f,'g',g,'mu',mu);
% Assemble the Stokes matrix
N = size(A,1);
M = size(B,1);
%StokesMatrix = [A, B'; B, sparse(M,M)];
% Define dimensions of B'
B_transpose = B'; % Transpose of B
% Construct a sparse MxM matrix (adjust M to your desired size)
M = size(B, 2); % Number of columns in B (or rows in B')
sparse_M_M = eye(M, M);
% Concatenate A, B_transpose, B, and sparse_M_M
StokesMatrix = [A, B_transpose; B, sparse_M_M];
% Define the right-hand side
rhs = [f; g];
% Solve the Stokes system
solution = StokesMatrix \ rhs;
% Extract the velocity and pressure components from the solution
u = solution(1:N);
p = solution(N+1:end);
0 Kommentare
Antworten (1)
Walter Roberson
am 27 Nov. 2023
A = [1 3 0 2; 0 2 1 0; 2 3 7 0; 9 1 0 0]
B = [0 2; 1 0; 4 7; 9 0]
N = size(A,1);
M = size(B,1);
B_transpose = B'; % Transpose of B
M = size(B, 2); % Number of columns in B (or rows in B')
sparse_M_M = eye(M, M);
whos A B_transpose
whos B sparse_M_M
StokesMatrix = [A, B_transpose; B, sparse_M_M];
Observe that A is 4 x 4 and B_transpose is 2 x 4; you cannot "," together matrices with different numbers of rows.
Observe that B is 4 x 2 and sparse_M_M is 2 x 2; you cannot "," together matrices with different numbers of rows.
If the code had been
StokesMatrix = [A, B; B_transposse, sparse_M_M];
then that would be 4 x 4, 4 x 2 which would give 4 x 6 on the top; the bottom would be 2 x 4, 2 x 2 which would give 2 x 6 on the bottom. Then 4 x 6 ";" 2 x 6 would give a 6 x 6 output, which sounds plausible.
14 Kommentare
Walter Roberson
am 29 Nov. 2023
Bearbeitet: Walter Roberson
am 29 Nov. 2023
I found ways to extract necessary ifem files to here so that we could run the code. However, ExtractBoundaryEdges is not part of that toolbox, and the only reference I find to it online is your question here.
%load in squaremesh from third-party site
tdir = tempdir;
addpath(tdir);
smfilename = fullfile(tdir, 'squaremesh.m');
websave(smfilename, 'https://github.com/lyc102/ifem/raw/master/mesh/squaremesh.m');
shfilename = fullfile(tdir, 'showmesh.m');
websave(shfilename, 'https://raw.githubusercontent.com/lyc102/ifem/master/tool/showmesh.m');
urfilename = fullfile(tdir, 'uniformrefine.m');
websave(urfilename, 'https://github.com/lyc102/ifem/raw/master/mesh/uniformrefine.m');
mufilename = fullfile(tdir, 'myunique.m');
websave(mufilename, 'https://raw.githubusercontent.com/lyc102/ifem/master/tool/myunique.m');
%back to our program
% Define problem parameters
mu = 1.0; % viscosity
f = [1; 0; 1; 2]; % forcing term for velocity
g = [2; 1]; % forcing term for pressure
%% Mesh creation
h = 0.25; % discretization step
[node,elem] = squaremesh([-1,1,-1,1],h);
subplot(1,2,1); showmesh(node,elem);
N = size(node,1);
NTc = size(elem,1);
% Uniform refinement to get a fine grid
[node,elem] = uniformrefine(node,elem);
% Uniform refinement to get a fine grid
subplot(1,2,2); showmesh(node,elem);
%%extract boundary conditions
[Dirichlet, bdFlag] = ExtractBoundaryEdges(elem, node);
%% Load Stokes Data Partial Differential Equations
%Different velocity and pressure equation models
pde = StokesModelMustafaSin;%StokesModelMustafa;%StokesModelMustafaNew;%Stokesdata3;%StokesModelMustafa2;%StokesModelMustafaSin;%;
option.elemType = 'isoP2P0';
option.fquadorder = 3;
tolerance = 1e-6;
max_iterations = 1000;
mesh.node = node;
mesh.elem = elem;
mesh.bdFlag = bdFlag;
%The following line of code is only run to get the eqn (equation structure) easier without defining it every time
[soln,eqn,info] = StokesisoP2P0(node,elem,bdFlag,pde, option);
% Define matrix A
A = [1 0 0 2; 1 2 1 0; 0 3 1 0; 0 1 0 1];
% Define matrix B
B = [1 2 1 0; 0 1 0 0];
% Assemble the finite element system
%[A, B, ~, ~] = assempde(pde,node,elem,'f',f,'g',g,'mu',mu);
% Assemble the Stokes matrix
N = size(A,1);
M = size(B,1);
%StokesMatrix = [A, B'; B, sparse(M,M)];
% Define dimensions of B'
B_transpose = B'; % Transpose of B
% Construct a sparse MxM matrix (adjust M to your desired size)
M = size(B, 1); % Number of columns in B (or rows in B')
sparse_M_M = sparse(M, M);
% Check dimensions of matrices involved in the problematic line
size_B = size(B);
size_A = size(A);
% Perform the matrix multiplication explicitly and verify dimensions
result = 2 * B / (A)* B_transpose;
size_result = size(result);
% Concatenate A, B_transpose, B, and sparse_M_M
StokesMatrix = [A, B_transpose; B, sparse_M_M];
% Define matrix M
%M = [2*eye(size(A,1)) 2*(A)\B_transpose; B 2*B*(A)\B_transpose];
M_top = [2 * eye(size(A, 1)), 2 * (A)\B_transpose];
M_bottom = [B, 2 * B / (A) * B_transpose];
M = [M_top; M_bottom];
% Define matrix F
F = [2 * A \ f; 2 * B / (A) * f - g];
% Define the right-hand side
rhs = [f; g];
% Solve the Stokes system
solution = StokesMatrix \ rhs;
sol = M \ F;
% Extract the velocity and pressure components from the solution
u = solution(1:N);
p = solution(N+1:end);
% Calculate the residual R
R = F - M * sol;
% Define matrix C
C = [0.5 * A zeros(size(transpose(B))); zeros(size(B)) eye(size(B,1))];
% Define an appropriate non-zero vector x
% Example: Assuming x has the same dimension as the column size of B
x = [1; 0; 1; 2; 3; 0]; % You can define x as per your problem's requirements
% Calculate x^T M x using the inner product function
x_transpose_M_x = inner_x_y(M, x, x);
% Check for positive definiteness
is_pos_definite_M = x_transpose_M_x > 0;
% Display the result
disp(['Is M positive definite: ', num2str(is_pos_definite_M)]);
function [inner_pro] = inner_x_y(C, x,y)
%C = [A 0 ; 0 1]
%(x,y) = (Cx,y) = (Ax_1,y_1) + (x_2,y_2) = transpose(x_1)*A*y_1 +transpose(x_2)*y_2
%trans(y_1)*x_2
inner_pro = transpose(y) .* C .* x;
end
Siehe auch
Kategorien
Mehr zu Sparse Matrices finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!