Error using horzcat Dimensions of arrays being concatenated are not consistent. Error in MKstokes (line 52) StokesMatrix = [A, B_transpose; B, sparse_M_M];

% 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);

Unrecognized function or variable 'squaremesh'.

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);

Answer 1

[Walter Roberson]

A = [1 3 0 2; 0 2 1 0; 2 3 7 0; 9 1 0 0]

A = 4×4

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]

B = 4×2

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

  Name             Size            Bytes  Class     Attributes

  A                4x4               128  double              
  B_transpose      2x4                64  double              

whos B sparse_M_M

  Name            Size            Bytes  Class     Attributes

  B               4x2                64  double              
  sparse_M_M      2x2                32  double              

StokesMatrix = [A, B_transpose; B, sparse_M_M];

Error using horzcat
Dimensions of arrays being concatenated are not consistent.

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.

[Mostafa Kassoum]

I have changed matrix B to B = [0 2 1 0; 4 7 9 0] and I got the same mistakes. Error using horzcat

Dimensions of arrays being concatenated are not consistent.

Error in MKstokes (line 52)

StokesMatrix = [A, B_transpose; B, sparse_M_M];

[Walter Roberson]

Unfortunately we do not have your squaremesh and other functions to test with.

[Mostafa Kassoum]

I can give you a link of sequaremesh to test the code. for the functions you can take any simple examples.

For finite element package and sequaremesh I used this link

https://lyc102.github.io/ifem/projects/Stokes/#part-i-isop2-p0-element

[Mostafa Kassoum]

Example of functions in the following code:

function pde = StokesModelMustafaSin
%% Custom Stokes equation model 
% 
pde = struct('f', 0, 'exactp', @exactp, 'exactu', @exactu,'g_D',@g_D,'g_N',@g_N);

% exact velocity
    function z = exactu(p)
        x = p(:,1); 
        y = p(:,2);
        z(:,1) = sin(pi().*x).*sin(pi().*y);
        z(:,2) = sin(pi().*x).*sin(pi().*y);
    end
% exact pressure
    function z = exactp(p)
        x = p(:,1);% y = p(:,2);
        z = 1-2*x;
    end
% Dirichlet boundary condition of velocity
    function z = g_D(p)
        z = exactu(p);
    end
    function z = g_N(p)
        z(:,1) = 2*ones(size(p,1),1);
        z(:,2) = zeros(size(p,1),1);
    end
end

[Voss]

A = [1 3 0 2; 0 2 1 0; 2 3 7 0; 9 1 0 0]

A = 4×4

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]

B = 2×4

0 2 1 0 4 7 9 0

B_transpose = B.';  % Transpose of B

M = size(B,1);
sparse_M_M = eye(M, M);

StokesMatrix = [A, B_transpose; B, sparse_M_M]

StokesMatrix = 6×6

1 3 0 2 0 4 0 2 1 0 2 7 2 3 7 0 1 9 9 1 0 0 0 0 0 2 1 0 1 0 4 7 9 0 0 1

[Mostafa Kassoum]

I got the same error. I don't know what is problem while I have done the required changed.

% Define problem parameters
mu = 1.0;     % viscosity
f = [2; 0; 1; 2];   % forcing term for velocity
g = [0; 0];        % forcing term for pressure
%% Mesh creation
h = 0.25; % discretization step
[node,elem] = squaremesh([-1,1,-1,1],h);

Unrecognized function or variable 'squaremesh'.

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);

[Walter Roberson]

We still do not have squaremesh() to test with.

[Voss]

@Mostafa Kassoum: At the time sparse_M_M is created by this line:

sparse_M_M = eye(M, M);

M needs to be size(B,1), but you have M as size(B,2) when that line executes.

[Mostafa Kassoum]

I added the following lines to the code

% Define matrix M

M = [2*eye(size(A,1)) 2*(A)\B_transpose; B 2*B*(A)\B_transpose];

And I got the following error,

Error using \

Matrix dimensions must agree.

Error in MKstokestest (line 54)

M = [2*eye(size(A,1)) 2*(A)\B_transpose; B 2*B*(A)\B_transpose];

[Voss]

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];
B_transpose = B.';

This line runs without error. Is it what you intended? I have no idea.

M = [2*eye(size(A,1))  2*(A)\B_transpose; B_transpose  2*B*(A)\B]

M = 8×6

2.0000 0 0 0 -0.0644 -0.1980 0 2.0000 0 0 0.5792 1.7822 0 0 2.0000 0 -0.1584 -0.0644 0 0 0 2.0000 -0.8366 -1.5743 0 4.0000 0 0 0 0 2.0000 7.0000 0 0 0 0 1.0000 9.0000 0 0.1111 0.0556 0 0 0 0.2500 -0.5347 0.0764 0

[Mostafa Kassoum]

% 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);

Unrecognized function or variable 'squaremesh'.

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

In this code I got the following error:

Error using horzcat

Dimensions of arrays being concatenated are not consistent.

Error in MKstokestest (line 92)

disp(['Is M positive definite: ', num2str(is_pos_definite_M)]);

[Torsten]

Why do you use [ ] in the disp command ? Do you want to build a row vector of strings ?

[Mostafa Kassoum]

What can I use instead?

[Walter Roberson]

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);

Unrecognized function or variable 'ExtractBoundaryEdges'.

%% 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