Hello Brian,<br><br>thank you for your answer. Things are getting clearer now. I was checking the equation and specially the part that you mentioned that was making you confused: "dt*(\alpha*f_{n+1} + (1-\alpha)*f_n)". I think that the equation written in the header:<br>
<br>( M + alpha*dt* K )*U_t=(M - (1.- alpha)*dt* K)* U_{t-1} + dt*(alpha*f_{n+1} + (1-alpha)*f_n)<br><br>does not totally correspond to what it is carried out:<br><br>( M + alpha*dt* K )*U_t=(M - (1.- alpha)*dt* K)* U_{t-1} + dt*(alpha*f_{t} + (1-alpha)*f_{t-1}) <br>
<br>So the difference would be f_{t} instead of f_{n+1} and f_{t-1} instead of f_n.<br><br>However I still have a question. Where does alpha come from? In the generic PDE there is no alpha. How can I know where should an alpha or 1-alpha be in the way from the PDE to discretized equation?<br>
<br>Thank you!<br><br>Cristina<br><br>Would it make sense to calcualte U_t by using f_{t+1}<br><br><div class="gmail_quote">On Wed, Dec 15, 2010 at 5:04 PM, brian avants <span dir="ltr"><<a href="mailto:stnava@gmail.com">stnava@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">Cristina<br>
<br>
Thanks for the reminder about what's written there --- it clearly<br>
should be updated. Some additional definitions:<br>
<br>
M --- the mass matrix ( covered in any basic FEM book or webpage ).<br>
<br>
rho --- this is a scalar that determines the "density" of the mass matrix.<br>
<br>
U --- the deformation field.<br>
<br>
t --- time.<br>
<br>
K --- the stiffness matrix.<br>
<br>
alpha --- in crank-nicolson, set to 0.5 , but often used with 0 or 1.<br>
<br>
A generic PDE this would solve would read something like:<br>
<br>
$ \rho \partial U / \partial t + L U = f $<br>
<br>
where \rho and L are determined by the physics of the problem. What<br>
I am confused about is this part:<br>
<div class="im"><br>
dt*(\alpha*f_{n+1} + (1-\alpha)*f_n)<br>
<br>
</div>The n should probably be replaced with t indices that are similar to<br>
those accumulating the U solution. The code seems to verify this:<br>
<br>
void SolverCrankNicolson::RecomputeForceVector(unsigned int index)<br>
{//<br>
Float ft = m_ls->GetVectorValue(index,ForceTIndex);<br>
Float ftm1 = m_ls->GetVectorValue(index,ForceTMinus1Index);<br>
Float utm1 = m_ls->GetVectorValue(index,DiffMatrixBySolutionTMinus1Index);<br>
Float f=m_deltaT*(m_alpha*ft+(1.-m_alpha)*ftm1)+utm1;<br>
m_ls->SetVectorValue(index , f, ForceTIndex);<br>
}<br>
<br>
<br>
If I find a reference that shows this particular discretization, I<br>
will share it. However, this approach solves a common static problem<br>
if alpha = 1 and if f is static. You can verify this yourself.<br>
<br>
<br>
Brian<br>
<div class="im"><br>
<br>
<br>
<br>
> " * FEM Solver for time dependent problems; uses Crank-Nicolson implicit<br>
> discretization scheme.<br>
> *It solves the following problem:<br>
> *<br>
> * ( M + alpha*dt* K )*U_t=(M - (1.- alpha)*dt* K)* U_{t-1} +<br>
> dt*(alpha*f_{n+1} + (1-alpha)*f_n) "<br>
<br>
</div>M is the element's mass matrix<br>
<div><div></div><div class="h5"><br>
><br>
> This problem formulation differs from the one obtained for the heat equation<br>
> (for example). So, I guess that it comes from an specific PDE by applying<br>
> Crank-Nicolson discretization scheme to it. My questions then are:<br>
> 1. From which PDE does it come from? Which one is the problem that it is<br>
> trying to solve?<br>
> 2. If the PDE is different depending on the element that one chooses, should<br>
> not one get different problem formulations to solve by choosing different<br>
> elements (PDE)? And then is it the solver generic or valid for only one PDE?<br>
><br>
> Thank you!<br>
><br>
> Cristina<br>
><br>
> On Fri, Dec 10, 2010 at 8:59 PM, brian avants <<a href="mailto:stnava@gmail.com">stnava@gmail.com</a>> wrote:<br>
>><br>
>> Hi Cristina<br>
>><br>
>> The specific PDE depends upon the specific elements you choose. For<br>
>> instance, if you set up a triangular mesh via the<br>
>> 2DC0LinearTriangularMembrane element, then you will use a membrane<br>
>> energy penalty term.<br>
>><br>
>> The Crank-Nicolson solver uses the Crank-Nicolson approach to<br>
>> accumulate the solution to the registration problem over time.<br>
>><br>
>> Wikipedia provides a reasonable general explanation of the C-N<br>
>> approach and gives an example for the heat equation;<br>
>><br>
>> <a href="http://en.wikipedia.org/wiki/Crank%E2%80%93Nicolson_method" target="_blank">http://en.wikipedia.org/wiki/Crank%E2%80%93Nicolson_method</a><br>
>><br>
>> B.<br>
>><br>
>><br>
>> On Fri, Dec 10, 2010 at 10:07 AM, Cristina Oyarzun<br>
>> <<a href="mailto:coyarzunlaura@googlemail.com">coyarzunlaura@googlemail.com</a>> wrote:<br>
>> > Hello,<br>
>> ><br>
>> > I am using itk FEM registration filter. This filter makes use of<br>
>> > itkFEMSolverCrankNicolson. In the header file of the solver one finds<br>
>> > the<br>
>> > equation that the solver tries to solve. However, one can only find<br>
>> > there<br>
>> > its Crank-Nicolson formulation but not the original PDE that leads to<br>
>> > that<br>
>> > equation.<br>
>> ><br>
>> > I tried to get the discretized equation by using the equation of motion<br>
>> > that<br>
>> > is described in a book and Crank-Nicolson scheme but did not arrive to<br>
>> > the<br>
>> > same solution. Could someone tell me which on is the original PDE that<br>
>> > is<br>
>> > solved in itkFEMSolverCrankNicolson?<br>
>> ><br>
>> > Thank you!!<br>
>> ><br>
>> > Greetings,<br>
>> ><br>
>> > Cristina<br>
>> ><br>
>> > _____________________________________<br>
>> > Powered by <a href="http://www.kitware.com" target="_blank">www.kitware.com</a><br>
>> ><br>
>> > Visit other Kitware open-source projects at<br>
>> > <a href="http://www.kitware.com/opensource/opensource.html" target="_blank">http://www.kitware.com/opensource/opensource.html</a><br>
>> ><br>
>> > Kitware offers ITK Training Courses, for more information visit:<br>
>> > <a href="http://www.kitware.com/products/protraining.html" target="_blank">http://www.kitware.com/products/protraining.html</a><br>
>> ><br>
>> > Please keep messages on-topic and check the ITK FAQ at:<br>
>> > <a href="http://www.itk.org/Wiki/ITK_FAQ" target="_blank">http://www.itk.org/Wiki/ITK_FAQ</a><br>
>> ><br>
>> > Follow this link to subscribe/unsubscribe:<br>
>> > <a href="http://www.itk.org/mailman/listinfo/insight-users" target="_blank">http://www.itk.org/mailman/listinfo/insight-users</a><br>
>> ><br>
>> ><br>
>><br>
>><br>
>><br>
>> --<br>
>> ß®∫∆π<br>
><br>
><br>
> _____________________________________<br>
> Powered by <a href="http://www.kitware.com" target="_blank">www.kitware.com</a><br>
><br>
> Visit other Kitware open-source projects at<br>
> <a href="http://www.kitware.com/opensource/opensource.html" target="_blank">http://www.kitware.com/opensource/opensource.html</a><br>
><br>
> Kitware offers ITK Training Courses, for more information visit:<br>
> <a href="http://www.kitware.com/products/protraining.html" target="_blank">http://www.kitware.com/products/protraining.html</a><br>
><br>
> Please keep messages on-topic and check the ITK FAQ at:<br>
> <a href="http://www.itk.org/Wiki/ITK_FAQ" target="_blank">http://www.itk.org/Wiki/ITK_FAQ</a><br>
><br>
> Follow this link to subscribe/unsubscribe:<br>
> <a href="http://www.itk.org/mailman/listinfo/insight-users" target="_blank">http://www.itk.org/mailman/listinfo/insight-users</a><br>
><br>
><br>
<br>
<br>
<br>
--<br>
ß®∫∆π<br>
</div></div></blockquote></div><br>