Filter löschen
Filter löschen

solving 2D Poisson equation having already generated a mesh

30 Ansichten (letzte 30 Tage)
malak
malak am 4 Feb. 2024
Bearbeitet: Shivansh am 25 Feb. 2024
So i need to resolve the 2D poisson equation using Finite Element Method, i have an L shaped domain . So first of all i generated Triangular free mesh on the domain. but i dont know how to solve the eqution using my script. The conditions at the boundary is : u(x,t)=0
thank you.
% mesh generation
N = 1;
Visualize_Plot = true; % se vero visualizza il grafico
% Input Boundary :
% mi fa scegliere il file text dalla cartella
% [file,path] = uigetfile('*.txt');
% boundary = importdata(fullfile(path,file));
% oppure importo direttamente conoscendo il nome del file
boundary = importdata(fullfile('coordinate.txt'));
% dentro al file .txt c'è una matrice con prima colonna le x e la seconda colonna le y
tic
% if Visualize_Plot
% figure
% end
[T]= CreateMesh(boundary,N,Visualize_Plot);
min = min(min(boundary(:,1)),min(boundary(:,2))); % calcolo il valore minimo delle coordinate
max = max(max(boundary(:,1)),max(boundary(:,2))); % il valore max delle coordinate
margXY = 0.1*(max - min); % margine del 10%
if Visualize_Plot
figure(1)
xlim([min-margXY , max+margXY]) % imposto il limite dell'asse delle x
ylim(xlim) % stessa cosa per l'asse delle y
xlabel('x')
ylabel('y')
title('Mesh generation')
end
toc
function [T] = CreateMesh(boundary,N,PLOT)
% matrice che ha al suo interno la descrizione del poligono
Geom = [2 % ci indica che è un poligono
size(boundary,1)-1 % Numero di segmenti
boundary(1:end-1,1) % x-coordinate
boundary(1:end-1,2)]; % y-coordinate
[dl,~] = decsg(Geom);
r = min(sqrt(sum((boundary(1:end-1,1:2) - boundary(2:end,1:2)).^2,2))); % calcolo la distanza minima tra 2 vertici consecutivi del contorno
% genero una mesh iniziale con
[p,~,t] = initmesh(dl,'Hgrad',1.3,'Hmax',r/N);
% Hmax : la dimensione massima della mesh
% Hgrad : Mesh growth rate
T = zeros(size(t,2),6);
hold on
for i=1:size(t,2) % su ciascun triangolo nella mesh (t è il num di trangoli)
T(i,:) = [p(1,t(1,i)),p(2,t(1,i)),p(1,t(2,i)),p(2,t(2,i)),p(1,t(3,i)),p(2,t(3,i))];
% memorizza i vertici del traingolo
if PLOT
plot([p(1,t(1,i)),p(1,t(2,i)),p(1,t(3,i)),p(1,t(1,i))],...
[p(2,t(1,i)),p(2,t(2,i)),p(2,t(3,i)),p(2,t(1,i))],...
'color',[0.5 0.5 0.5])
end
end
plot(boundary(:,1),boundary(:,2),'b')
hold off

Antworten (1)

Shivansh
Shivansh am 25 Feb. 2024
Bearbeitet: Shivansh am 25 Feb. 2024
Hi Malak,
It seems like you have generated a mesh without using the MATLAB inbuilt functions. To solve the Poisson equation using the Finite Element Method (FEM), you'll need to assemble the system of equations that arises from the weak form of the Poisson equation and apply the boundary conditions.
You can refer to the following steps to proceed further:
  1. Assemble the Stiffness Matrix and Load Vector:
  2. Apply Boundary Conditions.
  3. Solve the Linear System.
Here's a simplified version of the MATLAB code that you could use to perform these steps. You will need to define your own functions for computing the element stiffness matrix and load vector based on the specific form of your Poisson equation.
% Assuming you have the mesh points p and the triangulation t from your mesh generation
% p: 2xNp matrix containing the coordinates of the nodes
% t: 4xNt matrix containing the indices of the nodes for each triangle
% Number of points
Np = size(p, 2);
% Initialize the global stiffness matrix K and load vector F
K = sparse(Np, Np); % Using sparse for efficiency
F = zeros(Np, 1);
% Assemble the stiffness matrix and load vector
for i = 1:size(t, 2)
nodes = t(1:3, i); % Nodes of the i-th triangle
[Ke, Fe] = elementStiffnessMatrixAndLoadVector(p(:, nodes));
K(nodes, nodes) = K(nodes, nodes) + Ke;
F(nodes) = F(nodes) + Fe;
end
% Apply boundary conditions
% You will need to identify the indices of the nodes on the boundary
boundaryNodes = findBoundaryNodes(p, boundary);
K(boundaryNodes, :) = 0;
K(boundaryNodes, boundaryNodes) = speye(length(boundaryNodes), length(boundaryNodes));
F(boundaryNodes) = 0;
% Solve the linear system
u = K \ F;
% u now contains the solution at the mesh points
In this code, "elementStiffnessMatrixAndLoadVector" is a placeholder for a function that you need to implement which computes the local stiffness matrix "Ke" and load vector "Fe" for each element. The function "findBoundaryNodes" is also a placeholder for a function that identifies the nodes on the boundary of the domain.
You can also use to the PDE toolbox and refer to the inbuilt functions to generate mesh and solve the 2D Poisson equation. You can use the "generatemesh" function for mesh creation and "solvepde" tpo solve the eqyation. You can refer to the following documentation for more information:
  1. generatemesh: https://www.mathworks.com/help/pde/ug/generate-mesh.html.
  2. solvepde: https://www.mathworks.com/help/pde/ug/pde.pdemodel.solvepde.html.
I hope it helps!

Community Treasure Hunt

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

Start Hunting!