Why am I "Out of memory" because of an "Error using lu"

3 views (last 30 days)
I Need to solve this PDE (the way they recommend here: Solve PDE ), but I got the following Problem:
Here´s the full Code:
%%Model parameters
% https://de.mathworks.com/help/pde/ug/c-coefficient-for-systems-for-specifycoefficients.html#bu5tna0-1
global N
N = 3;
global nu
nu = 0.34;
global E
E = 70;
global rho
rho = 2.6989;
global A
A = rho * (1- 2 * nu) * (1 + nu) / E;
global B
B = 1 - nu;
global C
C = 1/2 - nu;
global D
D = 3/2 - 2 * nu;
% time grid
tlist = 0:0.025:0.1;
model = createpde(N);
% Import geometry into the container.
importGeometry(model,'zylinder.stl');
% View the geometry with face labels.
pdegplot(model,'FaceLabels','off')
%%Specify PDE Coefficients
% Include the PDE COefficients in |model|.
specifyCoefficients(model, 'm', @mCoeffFunction, 'd', 0, 'c', @cCoeffFunction, 'a', @aCoeffFunction, 'f', @fCoeffFunction);
% set initial conditions
u0 = [0.01; 0; 0];
ut0 = [0.1; 0; 0];
setInitialConditions(model,u0, ut0);
% Create zero Dirichlet boundary conditions on all faces.
applyBoundaryCondition(model,'dirichlet','face',1:3,'u',[0.01,0,0]); % noch falsche Randbedingung, da unklar wie einzugeben (Rb ist PDE).
% Create a mesh
generateMesh(model);
pdemesh(model)
solve the PDE
result = solvepde(model, tlist);
u = result.NodalSolution;
pdeplot3D(model,'ColorMapData',u(:,1))
title('u(1) --> u')
function mMatrix = mCoeffFunction(region, state)
% https://de.mathworks.com/help/pde/ug/m-d-or-a-coefficient-for-systems.html#bu5xcw2-1
% x is r here
global N A B C D
Nr = length(region.x);
mMatrix = zeros(N^2, Nr);
mMatrix(1, :) = -region.x .* A;
mMatrix(5, :) = -region.x .* A;
mMatrix(9, :) = -region.x .* A;
end
function dMatrix = dCoeffFunction(region, state)
% https://de.mathworks.com/help/pde/ug/m-d-or-a-coefficient-for-systems.html#bu5xcw2-1
global N A B C D
Nr = length(region.x);
dMatrix = zeros(N^2, Nr);
end
function cMatrix = cCoeffFunction(region, state)
% https://de.mathworks.com/help/pde/ug/m-d-or-a-coefficient-for-systems.html#bu5xcw2-1
global N A B C D
Nr = length(region.x);
cMatrix = zeros(9*N^2, Nr);
cMatrix(getRowOfCMatrix(1, 1, 1, 1), :) = B .* region.x;
cMatrix(getRowOfCMatrix(1, 1, 2, 2), :) = C ./ region.x;
cMatrix(getRowOfCMatrix(1, 1, 3, 3), :) = C .* region.x;
cMatrix(getRowOfCMatrix(1, 2, 1, 2), :) = 1/4;
cMatrix(getRowOfCMatrix(1, 2, 2, 1), :) = 1/4;
cMatrix(getRowOfCMatrix(1, 3, 1, 3), :) = 1/4 .* region.x;
cMatrix(getRowOfCMatrix(1, 3, 3, 1), :) = 1/4 .* region.x;
cMatrix(getRowOfCMatrix(2, 1, 1, 2), :) = 1/4;
cMatrix(getRowOfCMatrix(2, 1, 2, 1), :) = 1/4;
cMatrix(getRowOfCMatrix(2, 3, 3, 2), :) = 1/4;
cMatrix(getRowOfCMatrix(2, 3, 2, 3), :) = 1/4;
cMatrix(getRowOfCMatrix(2, 2, 1, 1), :) = C .* region.x;
cMatrix(getRowOfCMatrix(2, 2, 2, 2), :) = B ./ region.x;
cMatrix(getRowOfCMatrix(2, 2, 3, 3), :) = C .* region.x;
cMatrix(getRowOfCMatrix(3, 1, 3, 1), :) = 1/4 .* region.x;
cMatrix(getRowOfCMatrix(3, 1, 1, 3), :) = 1/4 .* region.x;
cMatrix(getRowOfCMatrix(3, 3, 1, 1), :) = C .* region.x;
cMatrix(getRowOfCMatrix(3, 3, 2, 2), :) = C ./ region.x;
cMatrix(getRowOfCMatrix(3, 3, 3, 3), :) = B .* region.x;
end
function aMatrix = aCoeffFunction(region, state)
% https://de.mathworks.com/help/pde/ug/m-d-or-a-coefficient-for-systems.html#bu5xcw2-1
global N A B C D
Nr = length(region.x);
aMatrix = zeros(N^2, Nr);
aMatrix(1, :) = -B ./ region.x;
aMatrix(5, :) = -C ./ region.x;
end
function fMatrix = fCoeffFunction(region, state)
% https://de.mathworks.com/help/pde/ug/m-d-or-a-coefficient-for-systems.html#bu5xcw2-1
% state.ux(1,:) --> ur, state.uy(2, :) --> vphi, ...
global N A B C D
Nr = length(region.x);
fMatrix = zeros(N, Nr);
fMatrix(1, :) = -B .* state.ux(1, :) + D * 1 ./ region.x .* state.uy(2, :);
fMatrix(2, :) = -D ./ region.x .* state.uy(1, :) - C .* state.ux(2, :);
fMatrix(3, :) = -1/2 .* state.uz(1, :) - C .* state.ux(3, :);
end
function row = getRowOfCMatrix(i, j, k, l)
global N
row = 9 * N * (j - 1) + 9 * i + 3 * l + k - 12;
end

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!