Why am I "Out of memory" because of an "Error using lu"
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Antworten (0)
Siehe auch
Kategorien
Mehr zu Eigenvalue Problems 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!