Can anyone help me to fix my 1D NonLinear General FEM Code?

6 views (last 30 days)
Federico Giai Pron
Federico Giai Pron on 16 Jul 2017
Answered: HAMZA NASSAR on 24 Jul 2018
Good Morning, I'm trying to write a simple (i.e. non optimized) non-linear general FEM code. I have written a FEM code for linear problems and now I want to extend my knowledge in the non-linear case. The code is very simple, i don't want to write the best possible code but, for the moment, only to learn the basis in FEM analysis. The code is iterative and based on the Newton's Method. It stars writing the Jacobian matrix. Each element of the matrix is given by the integrale of the Weak formulation in the domain. Obviously, the first and the N-th element are solved using the boundary conditions weak equations. Can anyone help me solving this problem? The code is:
clear all
close all
clc
%%Dominio
xi_Mesh=0;
xe_Mesh=50/1E+03;
N_Elementi_Mesh=19;
%%Formulazione debole
p=@(x,u,dudx,v,dvdx)(5);
q=@(x,u,dudx,v,dvdx)(0);
r=@(x,u,dudx,v,dvdx)(0);
f=@(x,u,dudx,v,dvdx)(5E+05);
% Non-Linear System of Equations
WF=@(x,u,dudx,v,dvdx)(p(x,u,dudx,v,dvdx)*dvdx*dudx+q(x,u,dudx,v,dvdx)*v*u+r(x,u,dudx,v,dvdx)*v*u-f(x,u,dudx,v,dvdx)*v);
%%Condizioni al contorno
ui=300;
ue=700;
WFi=@(x,u,dudx,v,dvdx)(u-ui);
WFe=@(x,u,dudx,v,dvdx)(u-ue);
%%Funzioni di base
Ordine=1;
%%Dominio di integrazione
N_Elementi_Integrale=N_Elementi_Mesh+100;
%%Equazione Non Lineare
u0=300;
h=1E-05;
Tollerance=1E-05;
n_Iterations_Max=5;
%%Soluzione
N_Nodi_Mesh=Ordine*N_Elementi_Mesh+1;
x_Mesh=linspace(xi_Mesh,xe_Mesh,N_Nodi_Mesh);
%%Dominio di integrazione
N_Nodi_Integrale=Ordine*N_Elementi_Integrale+1;
x_Integrale=linspace(xi_Mesh,xe_Mesh,N_Nodi_Integrale);
%%Inizializzazione delle variabili utilizzate all'interno del metodo iterativo
u=u0*ones(length(x_Mesh),1);
Error=1+Tollerance;
k=0;
%%Metodo di Newton
while Error>Tollerance & k<n_Iterations_Max
% Incremento del numero di iterazioni
k=k+1;
%%Interpolazione di u corrispondente all'iterazione k-1
u_Integrale=interp1(x_Mesh,u,x_Integrale);
dudx_Integrale=d(u_Integrale,x_Integrale);
%%Costruzione della matrice Jacobiana
Jacobian_Matrix=zeros(length(x_Mesh),length(x_Mesh));
for j_Mesh=1:1:length(x_Mesh)
j_Integrale=find(abs(x_Mesh(j_Mesh)-x_Integrale)==min(abs(x_Mesh(j_Mesh)-x_Integrale)));
u_down_Integrale=u_Integrale;
u_up_Integrale=u_Integrale;
u_down_Integrale(j_Integrale)=u_down_Integrale(j_Integrale)-h/2;
u_up_Integrale(j_Integrale)=u_up_Integrale(j_Integrale)+h/2;
dudx_down_Integrale=d(u_down_Integrale,x_Integrale);
dudx_up_Integrale=d(u_up_Integrale,x_Integrale);
for i=1:1:length(x_Mesh)
[phii,dphiidx]=phi_dphidx_fun(Ordine,i,N_Nodi_Mesh,x_Mesh,x_Integrale);
Integrale_down=0;
Integrale_up=0;
if i==1
for kk=1:1:length(x_Integrale)-1
%if phii(kk)~=0 || phii(kk+1)~=0
Integrale_down=Integrale_down+(WFi(x_Integrale(kk),u_down_Integrale(kk),dudx_down_Integrale(kk),phii(kk),dphiidx(kk))+WFi(x_Integrale(kk+1),u_down_Integrale(kk+1),dudx_down_Integrale(kk+1),phii(kk+1),dphiidx(kk+1)))/2*(x_Integrale(kk+1)-x_Integrale(kk));
Integrale_up=Integrale_up+(WFi(x_Integrale(kk),u_up_Integrale(kk),dudx_up_Integrale(kk),phii(kk),dphiidx(kk))+WFi(x_Integrale(kk+1),u_up_Integrale(kk+1),dudx_up_Integrale(kk+1),phii(kk+1),dphiidx(kk+1)))/2*(x_Integrale(kk+1)-x_Integrale(kk));
%end
end
elseif i==length(x_Mesh)
for kk=1:1:length(x_Integrale)-1
%if phii(kk)~=0 || phii(kk+1)~=0
Integrale_down=Integrale_down+(WFe(x_Integrale(kk),u_down_Integrale(kk),dudx_down_Integrale(kk),phii(kk),dphiidx(kk))+WFe(x_Integrale(kk+1),u_down_Integrale(kk+1),dudx_down_Integrale(kk+1),phii(kk+1),dphiidx(kk+1)))/2*(x_Integrale(kk+1)-x_Integrale(kk));
Integrale_up=Integrale_up+(WFe(x_Integrale(kk),u_up_Integrale(kk),dudx_up_Integrale(kk),phii(kk),dphiidx(kk))+WFe(x_Integrale(kk+1),u_up_Integrale(kk+1),dudx_up_Integrale(kk+1),phii(kk+1),dphiidx(kk+1)))/2*(x_Integrale(kk+1)-x_Integrale(kk));
%end
end
else
for kk=1:1:length(x_Integrale)-1
%if phii(kk)~=0 || phii(kk+1)~=0
Integrale_down=Integrale_down+(WF(x_Integrale(kk),u_down_Integrale(kk),dudx_down_Integrale(kk),phii(kk),dphiidx(kk))+WF(x_Integrale(kk+1),u_down_Integrale(kk+1),dudx_down_Integrale(kk+1),phii(kk+1),dphiidx(kk+1)))/2*(x_Integrale(kk+1)-x_Integrale(kk));
Integrale_up=Integrale_up+(WF(x_Integrale(kk),u_up_Integrale(kk),dudx_up_Integrale(kk),phii(kk),dphiidx(kk))+WF(x_Integrale(kk+1),u_up_Integrale(kk+1),dudx_up_Integrale(kk+1),phii(kk+1),dphiidx(kk+1)))/2*(x_Integrale(kk+1)-x_Integrale(kk));
%end
end
end
Jacobian_Matrix(i,j_Mesh)=(Integrale_up-Integrale_down)/h;
end
end
%%Calcolo di f(u)=Integrale(x(1),x(end))(WF)
for i=1:1:length(x_Mesh)
[phii,dphiidx]=phi_dphidx_fun(Ordine,i,N_Nodi_Mesh,x_Mesh,x_Integrale);
y(1,i)=0;
if i==1
for kk=1:1:length(x_Integrale)-1
%if phii(kk)~=0 || phii(kk+1)~=0
y(1,i)=y(1,i)+(WFi(x_Integrale(kk),u_Integrale(kk),dudx_Integrale(kk),phii(kk),dphiidx(kk))+WFi(x_Integrale(kk+1),u_Integrale(kk+1),dudx_Integrale(kk+1),phii(kk+1),dphiidx(kk+1)))/2*(x_Integrale(kk+1)-x_Integrale(kk));
%end
end
elseif i==length(x_Mesh)
for kk=1:1:length(x_Integrale)-1
%if phii(kk)~=0 || phii(kk+1)~=0
y(i)=y(1,i)+(WFe(x_Integrale(kk),u_Integrale(kk),dudx_Integrale(kk),phii(kk),dphiidx(kk))+WFe(x_Integrale(kk+1),u_Integrale(kk+1),dudx_Integrale(kk+1),phii(kk+1),dphiidx(kk+1)))/2*(x_Integrale(kk+1)-x_Integrale(kk));
%end
end
else
for kk=1:1:length(x_Integrale)-1
%if phii(kk)~=0 || phii(kk+1)~=0
y(i)=y(1,i)+(WF(x_Integrale(kk),u_Integrale(kk),dudx_Integrale(kk),phii(kk),dphiidx(kk))+WF(x_Integrale(kk+1),u_Integrale(kk+1),dudx_Integrale(kk+1),phii(kk+1),dphiidx(kk+1)))/2*(x_Integrale(kk+1)-x_Integrale(kk));
%end
end
end
end
%%Calcolo di u
u=u-Jacobian_Matrix/y;
%%Calcolo dell' errore
Error=norm(Jacobian_Matrix/y);
end
figure;
plot(x_Mesh,u);
The functions used are:
function [phi,dphidx]=phi_dphidx_fun(Ordine,k,N_Nodi_Mesh,x_Mesh,x)
% Ordine=1
switch Ordine
case 1
for i=1:1:length(x)
if k==1
if x(i)<x_Mesh(k)
phi(i)=NaN;
dphidx(i)=NaN;
elseif x(i)>=x_Mesh(k) & x(i)<=x_Mesh(k+1)
phi(i)=1-(x(i)-x_Mesh(k))/(x_Mesh(k+1)-x_Mesh(k));
dphidx(i)=0-(1-0)/(x_Mesh(k+1)-x_Mesh(k));
else
phi(i)=0;
dphidx(i)=0;
end
elseif k==N_Nodi_Mesh
if x(i)>x_Mesh(k)
phi(i)=NaN;
dphidx(i)=NaN;
elseif x(i)>=x_Mesh(k-1) & x<=x_Mesh(k)
phi(i)=(x(i)-x_Mesh(k-1))/(x_Mesh(k)-x_Mesh(k-1));
dphidx(i)=(1-0)/(x_Mesh(k)-x_Mesh(k-1));
else
phi(i)=0;
dphidx(i)=0;
end
else
if x(i)>=x_Mesh(k-1) & x(i)<x_Mesh(k)
phi(i)=(x(i)-x_Mesh(k-1))/(x_Mesh(k)-x_Mesh(k-1));
dphidx(i)=(1-0)/(x_Mesh(k)-x_Mesh(k-1));
elseif x(i)==x_Mesh(k)
phi(i)=1;
dphidx(i)=0;
elseif x(i)>x_Mesh(k) & x(i)<=x_Mesh(k+1)
phi(i)=1-(x(i)-x_Mesh(k))/(x_Mesh(k+1)-x_Mesh(k));
dphidx(i)=0-(1-0)/(x_Mesh(k+1)-x_Mesh(k));
else
phi(i)=0;
dphidx(i)=0;
end
end
end
%%Ordine=2
case 2
for i=1:1:length(x)
% k=1
if k==1
if x(i)<x_Mesh(k)
phi(i)=NaN;
dphidx(i)=NaN;
elseif x(i)>=x_Mesh(k) & x(i)<=x_Mesh(k+Ordine)
for ji=1:1:Ordine+1
for jj=1:1:Ordine+1
A(ji,jj)=x_Mesh(ji)^(Ordine+1-jj);
end
end
b(:,1)=zeros(1,Ordine+1);
b(1)=1;
c=A\b;
phi(i)=0;
dphidx(i)=0;
for ji=1:1:Ordine+1
phi(i)=phi(i)+c(ji)*x(i)^(Ordine+1-ji);
dphidx(i)=dphidx(i)+(Ordine+1-ji)*c(ji)*x(i)^(Ordine+1-ji-1);
end
else
phi(i)=0;
dphidx(i)=0;
end
% k=2*N_Elementi_Mesh+1
elseif k==2*N_Nodi_Mesh-1
if x(i)>x_Mesh(k)
phi(i)=NaN;
dphidx(i)=NaN;
elseif x(i)>=x_Mesh(k-Ordine) & x<=x_Mesh(k)
for ji=1:1:Ordine+1
for jj=1:1:Ordine+1
A(ji,jj)=x_Mesh(2*N_Nodi_Mesh-1-Ordine-1+ji)^(Ordine+1-jj);
end
end
b(:,1)=zeros(1,Ordine+1);
b(end)=1;
c=A\b;
phi(i)=0;
dphidx(i)=0;
for ji=1:1:Ordine+1
phi(i)=phi(i)+c(ji)*x(i)^(Ordine+1-ji);
dphidx(i)=dphidx(i)+(Ordine+1-ji)*c(ji)*x(i)^(Ordine+1-ji-1);
end
else
phi(i)=0;
dphidx(i)=0;
end
else
%%k!=1, k!=2*N_Nodi_Mesh+1
if mod(k,2)==1 & k~=1 & k~=2*N_Nodi_Mesh-1 % k: Dispari
if x(i)>=x_Mesh(k-Ordine) & x(i)<x_Mesh(k)
for ji=1:1:Ordine+1
for jj=1:1:Ordine+1
A(ji,jj)=x_Mesh(k-Ordine-1+ji)^(Ordine+1-jj);
end
end
b(:,1)=zeros(1,Ordine+1);
b(end)=1;
c=A\b;
phi(i)=0;
dphidx(i)=0;
for ji=1:1:Ordine+1
phi(i)=phi(i)+c(ji)*x(i)^(Ordine+1-ji);
dphidx(i)=dphidx(i)+(Ordine+1-ji)*c(ji)*x(i)^(Ordine+1-ji-1);
end
elseif x(i)==x_Mesh(k)
phi(i)=1;
dphidx(i)=0;
elseif x(i)>x_Mesh(k) & x(i)<=x_Mesh(k+Ordine)
for ji=1:1:Ordine+1
for jj=1:1:Ordine+1
A(ji,jj)=x_Mesh(k-1+ji)^(Ordine+1-jj);
end
end
b(:,1)=zeros(1,Ordine+1);
b(1)=1;
c=A\b;
phi(i)=0;
dphidx(i)=0;
for ji=1:1:Ordine+1
phi(i)=phi(i)+c(ji)*x(i)^(Ordine+1-ji);
dphidx(i)=dphidx(i)+(Ordine+1-ji)*c(ji)*x(i)^(Ordine+1-ji-1);
end
else
phi(i)=0;
dphidx(i)=0;
end
elseif mod(k,2)==0 % k: Pari
if x(i)>=x_Mesh(k-Ordine+1) & x(i)<=x_Mesh(k+Ordine-1)
for ji=1:1:Ordine+1
for jj=1:1:Ordine+1
A(ji,jj)=x_Mesh(k-Ordine+ji)^(Ordine+1-jj);
end
end
b(:,1)=zeros(1,Ordine+1);
b(Ordine)=1;
c=A\b;
phi(i)=0;
dphidx(i)=0;
for ji=1:1:Ordine+1
phi(i)=phi(i)+c(ji)*x(i)^(Ordine+1-ji);
dphidx(i)=dphidx(i)+(Ordine+1-ji)*c(ji)*x(i)^(Ordine+1-ji-1);
end
else
phi(i)=0;
dphidx(i)=0;
end
end
end
end
end
end
And a simple Finite Difference approximation of the first derivative:
function [dydx]=d(y,x)
for i=1:1:length(x)
if i==1
dydx(i)=(y(i+1)-y(i))/(x(i+1)-x(i));
elseif i==length(x)
dydx(i)=(y(i)-y(i-1))/(x(i)-x(i-1));
else
These functions works weel because I used them to solve linear problems and the codes were ok. Thanlk you very much. dydx(i)=(y(i+1)-y(i-1))/(x(i+1)-x(i-1)); end end end

Sign in to answer this question.

Answers (1)

HAMZA NASSAR
HAMZA NASSAR on 24 Jul 2018
if true
% code
function [dydx]=d(y,x)
for i=1:1:length(x)
if i==1
dydx(i)=(y(i+1)-y(i))/(x(i+1)-x(i));
elseif i==length(x)
dydx(i)=(y(i)-y(i-1))/(x(i)-x(i-1));
else
end
end
end

Categories

Find more on Data Type Identification in Help Center and File Exchange

Tags

Products

Community Treasure Hunt

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

Start Hunting!