10:30 am
So, the first problem that I am having is that the radius of the artery that I am going to be receiving will be changing for every instance of simulation. I need to pass it on to a file which is called in the embedded matlab block which inputs this changing radius.
To make it sound simple, I will create a simpler scenario, make it work on my scenario, and then test it on the actual case. Now, what we have currently is a file that takes in a pulse generated value, and a sin wave, sums it and gives the output to a scope.
I'll do this, let the pulse output be the Arterial Cross Sectiional area variations. Let the sin wave be the pressure from the heart. The embedded MATLAB block will take the CS, calculate radius, save it to the workspace as the variable r, and call another function baby. Baby will take the updated radius increment it by 10. This incremented value will be multiplied to the sin wave and output will be shown.
So, what we have to do in simulink is
1. Sin block
2. Pulse generator block
3. Embedded matlab block with the above as inputs.
4. Calculate radius and save into a workspace variable and call baby
5. Baby takes radius, multiplies by 10 and returns value.
6. This value is then multiplied to sin wave and shown as output.
In all of the above steps, steps 4 and 5 are difficult. How do we do that?
Tuesday 25 December 2012
Saturday 22 December 2012
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.
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.
Friday 21 December 2012
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. :)
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. :)
Thursday 20 December 2012
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.
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.
Wednesday 19 December 2012
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.
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.
Saturday 17 November 2012
Joining the Pieces of the Puzzle
12:50 am
So, have been delaying the real implementation so much. Why? I know what all has to be done, how it has to be done, but not sure if it will all work together. So, confirm it all, I have to implement it. But I have been delaying it for so long. Why? I know why.. It’s only because, the work that I have planned is scattered among all the blogs that I have written. Need to put it all together and sit and implement it. Was hoping I could finish all by today, but looks like that’s not happening!! So, today I enlist the things I have to do regarding the model. Here we go.
Step 1 : c1(p) block
Create a polynomial block with given values of constants, to be multiplied by pressure p which is the input from fetal heart.
Step 2 : c2(z) block
Create a waveform that approximates the waveform given in the paper. The values could be:
z = [0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80]
c2(z) = [1 3.5 7.5 10 10.2 10.3 10.4 10.5 10.6 10.7 10.8 10.9 11 11.1 11.2 11.3 11.4]
Once that is ready, approximate the input z to the closest value and give out the corresponding c2(z)
Step 3 : c2(p,z) block
multiply the outputs of step 1 and step 2. This gives a continuous output which is a variable dependent on p.
Step 4 : Location Block
Take the z input and find which arterial segment or bifurcation it belongs to (switch case), and give the distance in that segment as the output.
Step 5 : Initial Cross Section block
Find its initial cross section according to formula for that segment or bifurcation with the output from step 4.
Step 6 : Final Cross Section Block
Multiply the output from step 3 to the blood density constant continuously, find its reciprocal, add it to 1, and then multiply the whole to the output of Step 5. This gives a continuously varying cross section variable.
Step 7 : Stokes – Navier block
Write the S-function which calculates the solution to the stokes navier equation using the Thermal regulator: PDE Toolbox - Simulink application and the following:
Step 8 : Export to the workspace and store, and plot against time
ONLY THE PROBLEM IS PLOTTING AGAINST TIME!!! DISCUSS WITH SIR!!
Friday 16 November 2012
Confirming the 4th Change in Direction and Some Relief!!
11:05 pm
So, after trying very hard with the method of characteristics used in the original paper, I realized that, the mathematical part is too difficult for a child like me. So, will not be solving the stokes Navier equation using method of characteristics but will be doing it using pdetool in Simulink. This was the same idea that I had posted before. I realize that this means it is a waste of at least a week, but, I’ll call it use,.. since I know, exactly WHY I’m not using the method of characteristics. Have to enlist the reasons.
So, today, before going to sleep, the agenda is, use Simulink to create a module for the Generates the solution to the Stokes-Navier equations. I think the designing part will be easy, but what we need to verify is whether the output we are getting is right or wrong. To do that, we need to connect the entire system together, or give a temporary input that represents the practical input to the module.
1. Think how to calculate mean pressure at any point in the system which I suppose becomes the input to the Stokes-Navier module.
2. Think up if and how this module and the rest of the modules can work together.
3. Form the Stokes-Navier module, using the document mentioned before.
This document will be edited after about half an hour.
So, here I am, trying to design the Stokes – Navier module. Firstly consider this equation given in the document:
Now consider the Stokes-Navier Equation given in one of the ppts that I’m referring to which helps me understand the meaning of the equation:
Consider this Stokes-Navier equation which is given as the prototype in the mathwork manual with good explaination:
Now consider this Stokes-Navier Equation which is given as the prototype in the pdetool with notations which are somewhat different:
Now I have to create an equivalency to match the coefficients. I’m guessing that it’s as follows:
rho = d = density of blood
mu = c = viscosity of blood
The controversy is about a and f. Going by the signs of the terms, the pde toolbox help document and the ppt that I’m referring to, I guess,
p = a = pressure and
f = f = other forces
Now, hoping that these assumptions are correct, I’ve got to substitute the values of these constants. These values are, unfortunately not given in the paper, so I will have to borrow from google, possibly, wikipedia.
Taking, d = 1060 kg/m3, c = 
(from wikipedia), a will be the time varying mean arterial input pressure, f I think we’ll consider as 1(not sure).
(from wikipedia), a will be the time varying mean arterial input pressure, f I think we’ll consider as 1(not sure).
Just figured out how to apply the elasticity to the generated Cross Section of the vessel.
With the Cross section and space coordinate equations, we can find out a value S of the cross section. Now, as per the formula given in the base paper, we have,
This equation defines elasticity as the steady cross section value, divided by the product of blood density and change in the cross section due to pressure. Now consider the following:
Here, as we know, S becomes the input, we have to have a module for calculation of c2.
From here I realise that I have solved the problem of Elasticity application to vessels.
Then, there was this problem about finding the mean pressure.. well just realized how much a of a fool I am to think that pressure in an artery can come from two directions!! Haha.. So, this problem will come only at the placenta. But there too, we have a division as maternal side of the placenta, and the fetal side of the placenta, so, nowhere we will have to encounter a scenario where we will have to calculate the ‘mean’ pressure! HA! Such a fool! Shi.
So now that I have figured out almost the entire deal I am at peace. Aaah! This feels great! There’s a time for everything. Jab jab jo jo hona hai, tab tab so, so, hota hai. 
Happy!! So, Again! I’ll do all of this tomorrow!!
Happy!! So, Again! I’ll do all of this tomorrow!!
Now what remains is the thinking up that I have to do regarding the bifurcations.. There is nothing about the bifurcations that is mentioned in the paper! Only the geometrical dimentions are given! What the hell is that all about? I guess, it should be treated in the same way as the arterial segments, but the only difference is that, there are two different values of cross sections at the second end. yoyo! Very happy and relieved! Everything has its own time!
ATB!
Subscribe to:
Posts (Atom)