11:43 am
Hey!! Good morning! Back after indeed a very very long time! Had my 3rd semester exams; and then a long vacation of an entire week. Pampered myself for a while. Now back to the project. Pressure needs to build up for me to get back on track. Spoke to my cousin sister yesterday, and finally decided that I should start.
So, I have sort of implemented all that I have written in the previous post, only the most important one, the last block of pde solution isn't done. I tried to do it.
Firstly, followed the same steps they have given in the pdetool, parabolic part in the manual. It's missing the plotting code. So realised I need to understand the bloody code. If I don't understand it, it's senseless. To understand the code, .. I HAVE to blog! Else my mind wanders around. So, here we go..
First line in the code is:
[p,e,t]=initmesh('circleg'); % Create initial triangular mesh
As we all know, to apply the pde on every point of the geometry, we first divide the geometry into small triangles. This we call meshing. We apply the pde on every triangle and then integrate the result. This is the method of finite element I guess. Need to search a bit on what it exactly is, but I think that's it. And this function 'initmesh' initializes the creation of this mesh.
The function obviously needs the geometry as the input parameter, which we can define in different ways according to the format and flexibility of the function. The output of the function is the meshed data; p => point matrix with x and y coordinates of mesh points as the rows of the matrix, e => edge matric with starting and ending point indices, starting and ending point parameter values, edge number and left and right side subdomain number, t => triangle matrix with three corner points in anticlockwise direction and subdomain number.
Next,
[p,e,t]=refinemesh('circleg',p,e,t); % Refine triangular mesh
Refinement usually bisects every edge of the specified triangles, thus creating 4 new smaller triangles out of every specified triangle. Hence the input to the function is the geometry under consideration, and the specified initial mesh and the output of the function is the p, e, t matrices of the refined mesh for that geometry.
Next,
u0=zeros(size(p,2),1);
size(p,2) returns the size of the 2-D p matrix and zeros(size(p,2),1) returns a 2-D zero matrix of the same size. We name it u0.
Next,
ix=find(sqrt(p(1,:).^2+p(2,:).^2)<0.4); % Find indices and values of nonzero elements
This finds all the non-zero points inside the circle with radius 0.4 units, and returns its indices as a column vector.
Next,
u0(ix)=ones(size(ix));
We replace the found non-zero entries of interest by 1 in the zero matrix u0
Next,
tlist=linspace(0,0.1,20);
create an array that will have values of 0 to 0.1 in 20 steps.
Last,
u1=parabolic(u0,tlist,'circleb1',p,e,t,1,0,1,1); % Solve parabolic PDE problem
parabolic is a function that produces the solution to the FEM formulation of the scalar PDE problem. Its input parameters are,
u0 => initial value of the mesh
tlist => time for which solution to be obtained
circleb1 => Initial Boundary condition data
p, e, t => mesh point, edge, triangle data
1, 0, 1, 1 => c, a, f, d
c => viscosity
a => gravity
f => pressure
d => density
u1 => final solution at the end of simulation time.
Now I need to plot this u1.
Hey!! Good morning! Back after indeed a very very long time! Had my 3rd semester exams; and then a long vacation of an entire week. Pampered myself for a while. Now back to the project. Pressure needs to build up for me to get back on track. Spoke to my cousin sister yesterday, and finally decided that I should start.
So, I have sort of implemented all that I have written in the previous post, only the most important one, the last block of pde solution isn't done. I tried to do it.
Firstly, followed the same steps they have given in the pdetool, parabolic part in the manual. It's missing the plotting code. So realised I need to understand the bloody code. If I don't understand it, it's senseless. To understand the code, .. I HAVE to blog! Else my mind wanders around. So, here we go..
First line in the code is:
[p,e,t]=initmesh('circleg'); % Create initial triangular mesh
As we all know, to apply the pde on every point of the geometry, we first divide the geometry into small triangles. This we call meshing. We apply the pde on every triangle and then integrate the result. This is the method of finite element I guess. Need to search a bit on what it exactly is, but I think that's it. And this function 'initmesh' initializes the creation of this mesh.
The function obviously needs the geometry as the input parameter, which we can define in different ways according to the format and flexibility of the function. The output of the function is the meshed data; p => point matrix with x and y coordinates of mesh points as the rows of the matrix, e => edge matric with starting and ending point indices, starting and ending point parameter values, edge number and left and right side subdomain number, t => triangle matrix with three corner points in anticlockwise direction and subdomain number.
Next,
[p,e,t]=refinemesh('circleg',p,e,t); % Refine triangular mesh
Refinement usually bisects every edge of the specified triangles, thus creating 4 new smaller triangles out of every specified triangle. Hence the input to the function is the geometry under consideration, and the specified initial mesh and the output of the function is the p, e, t matrices of the refined mesh for that geometry.
Next,
u0=zeros(size(p,2),1);
size(p,2) returns the size of the 2-D p matrix and zeros(size(p,2),1) returns a 2-D zero matrix of the same size. We name it u0.
Next,
ix=find(sqrt(p(1,:).^2+p(2,:).^2)<0.4); % Find indices and values of nonzero elements
This finds all the non-zero points inside the circle with radius 0.4 units, and returns its indices as a column vector.
Next,
u0(ix)=ones(size(ix));
We replace the found non-zero entries of interest by 1 in the zero matrix u0
Next,
tlist=linspace(0,0.1,20);
create an array that will have values of 0 to 0.1 in 20 steps.
Last,
u1=parabolic(u0,tlist,'circleb1',p,e,t,1,0,1,1); % Solve parabolic PDE problem
parabolic is a function that produces the solution to the FEM formulation of the scalar PDE problem. Its input parameters are,
u0 => initial value of the mesh
tlist => time for which solution to be obtained
circleb1 => Initial Boundary condition data
p, e, t => mesh point, edge, triangle data
1, 0, 1, 1 => c, a, f, d
c => viscosity
a => gravity
f => pressure
d => density
u1 => final solution at the end of simulation time.
Now I need to plot this u1.
No comments:
Post a Comment