Why am I "Out of memory" because of an "Error using lu"
1 view (last 30 days)
Show older comments
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
0 Comments
Answers (0)
See Also
Categories
Find more on Eigenvalue Problems in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!