Need the complete program today

10:25 am

Referring to this for understand how to solve pde programmatically.
Really scared now.. :(
Okay, So have been going in steps.

Right now, the following program works.


[p,e,t]=initmesh('Mycircleg');
[p,e,t]=refinemesh('Mycircleg',p,e,t);
u0=zeros(size(p,2),1);
tlist=linspace(1,0.1,20);
u1=parabolic(u0,tlist,'circleb1',p,e,t,0.0035,1,f,1060);
pdesurf(p,t,u1)

What I have to do in this, are the following modifications:

1. The radius of the circle in Mycircleg.m need to be passed somehow from the main Embedded MATLAB file.
One way to do this is save the input radius variable in as a matlab variable in workspace and access it in the next file.

2. In the parabolic function, a number of replacements have to be done. We'll begin from the end. Last 4 parameters of the function are c, a, f, d.
c = viscosity of blood, 0.0035 kg/m.s. a is gravity = 1, pressure f  needs to be provided dynamically from the pressure source, d = 1060kg/m3
One way to do this is to just have another input to the block, and use that variable in place of f.

3. tlist variable!! We know that the function parabolic will evaluate the solution to the geometry for every value of the array tlist. but I need the function to be evaluated for every value of the simulation time!
I will have to try different combinations of possibilities for this.

Trying to write the s function, and serendipity: Embedded MATLAB function

10:10 am

So, I know the template to be used. I know what I want to do. Just need the courage to start doing it. Here we go.



Phusssssssssss.................... Phusss, Phusssss, ...

So fooooolish of me!! Easier way is to use Embedded MATLAB function. Already tried it out with the example that I considered in the last post. Now, to my pde function. :)

Analysis to plot the pde parabolic solution

10:15 am

Now that I have understood what the code snippette does, we've gotto think about what the output should be.. It should be the velocity profile during the simulation time.

This makes me realize that, the time tlist that we supply to the parabolic function for solution of the pde should be changed to the simulation time. => How?

 Each column in u1 is the solution at the time given by the corresponding item in tlist. For my function, parabolic, each item in tlist needs to be block.CurrentTime. What should be the size of this array? Depends on the number of itertions of the solution we want to see. I need to try this out, on a smaller scale before implementing it on my model. Need an example.

let x be a sin wave, let y be a square wave. Let us try to build a block that will create an output that is a sum of these two inputs and the simulation time taken as an input through the functions arguments. ok. Let's try that out right away, .. but for that, I need to know how,.. how EXACTLY to write a level 2 matlab s function. Here we go.

Alright, so, read through and have realized that, it's a whole new set of keywords and almost an entire language to learn, .. when writiing an s function.. Remember having a s-function writing wizard of some sort in simulink.. It will be better if I use that,.. instead of wasting my time writing the entire code, .. let's check that out. No, .. no, .. that looks even more scary, coz it uses C mex and c code. Let that be. I'll continue with writing the s function matlab code only. Studying how to write an s-function from this pdf document. Very helpful.

Writing the s-function as m-file

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.