[Insight-users] FEM: PDE of Crank-Nicolson solver

[brian avants]

Cristina

Thanks for the reminder about what's written there --- it clearly
should be updated.   Some additional definitions:

M --- the mass matrix ( covered in any basic FEM book or webpage ).

rho --- this is a scalar that determines the "density" of the mass matrix.

U --- the deformation field.

t   --- time.

K  --- the stiffness matrix.

alpha ---  in crank-nicolson, set to 0.5 , but often used with 0 or 1.

A generic PDE this would solve would read something like:

$ \rho \partial U / \partial t  +  L U  =  f $

where \rho and  L are determined by the physics of the problem.   What
I am confused about is this part:

dt*(\alpha*f_{n+1} + (1-\alpha)*f_n)

The n should probably be replaced with t indices that are similar to
those accumulating the U solution.  The code seems to verify this:

void  SolverCrankNicolson::RecomputeForceVector(unsigned int index)
{//
  Float ft   = m_ls->GetVectorValue(index,ForceTIndex);
  Float ftm1 = m_ls->GetVectorValue(index,ForceTMinus1Index);
  Float utm1 = m_ls->GetVectorValue(index,DiffMatrixBySolutionTMinus1Index);
  Float f=m_deltaT*(m_alpha*ft+(1.-m_alpha)*ftm1)+utm1;
  m_ls->SetVectorValue(index , f, ForceTIndex);
}


If I find a reference that shows this particular discretization, I
will share it.  However, this approach solves a common static problem
if alpha = 1 and if f is static.   You can verify this yourself.


Brian




> " * FEM Solver for time dependent problems; uses Crank-Nicolson implicit
> discretization scheme.
>  *It solves the following problem:
>  *
>  *      ( M + alpha*dt* K )*U_t=(M - (1.- alpha)*dt* K)* U_{t-1} +
> dt*(alpha*f_{n+1} + (1-alpha)*f_n) "

M  is the element's mass matrix

>
> This problem formulation differs from the one obtained for the heat equation
> (for example). So, I guess that it comes from an specific PDE by applying
> Crank-Nicolson discretization scheme to it. My questions then are:
> 1. From which PDE does it come from? Which one is the problem that it is
> trying to solve?
> 2. If the PDE is different depending on the element that one chooses, should
> not one get different problem formulations to solve by choosing different
> elements (PDE)? And then is it the solver generic or valid for only one PDE?
>
> Thank you!
>
> Cristina
>
> On Fri, Dec 10, 2010 at 8:59 PM, brian avants <stnava at gmail.com> wrote:
>>
>> Hi Cristina
>>
>> The specific PDE depends upon the specific elements you choose.  For
>> instance, if you set up a triangular mesh via the
>> 2DC0LinearTriangularMembrane  element, then you will use a membrane
>> energy penalty term.
>>
>> The Crank-Nicolson solver uses the Crank-Nicolson approach to
>> accumulate the solution to the registration problem over time.
>>
>> Wikipedia provides a reasonable general explanation of the C-N
>> approach and gives an example for the heat equation;
>>
>> http://en.wikipedia.org/wiki/Crank%E2%80%93Nicolson_method
>>
>> B.
>>
>>
>> On Fri, Dec 10, 2010 at 10:07 AM, Cristina Oyarzun
>> <coyarzunlaura at googlemail.com> wrote:
>> > Hello,
>> >
>> > I am using itk FEM registration filter. This filter makes use of
>> > itkFEMSolverCrankNicolson. In the header file of the solver one finds
>> > the
>> > equation that the solver tries to solve. However, one can only find
>> > there
>> > its Crank-Nicolson formulation but not the original PDE that leads to
>> > that
>> > equation.
>> >
>> > I tried to get the discretized equation by using the equation of motion
>> > that
>> > is described in a book and Crank-Nicolson scheme but did not arrive to
>> > the
>> > same solution. Could someone tell me which on is the original PDE that
>> > is
>> > solved in itkFEMSolverCrankNicolson?
>> >
>> > Thank you!!
>> >
>> > Greetings,
>> >
>> > Cristina
>> >
>> > _____________________________________
>> > Powered by www.kitware.com
>> >
>> > Visit other Kitware open-source projects at
>> > http://www.kitware.com/opensource/opensource.html
>> >
>> > Kitware offers ITK Training Courses, for more information visit:
>> > http://www.kitware.com/products/protraining.html
>> >
>> > Please keep messages on-topic and check the ITK FAQ at:
>> > http://www.itk.org/Wiki/ITK_FAQ
>> >
>> > Follow this link to subscribe/unsubscribe:
>> > http://www.itk.org/mailman/listinfo/insight-users
>> >
>> >
>>
>>
>>
>> --
>> ß®∫∆π
>
>
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.html
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users
>
>



-- 
ß®∫∆π