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

[Cristina Oyarzun]

Hello Brian,

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:

( M + alpha*dt* K )*U_t=(M - (1.- alpha)*dt* K)* U_{t-1} + dt*(alpha*f_{n+1}
+ (1-alpha)*f_n)

does not totally correspond to what it is carried out:

( M + alpha*dt* K )*U_t=(M - (1.- alpha)*dt* K)* U_{t-1} + dt*(alpha*f_{t} +
(1-alpha)*f_{t-1})

So the difference would be f_{t} instead of f_{n+1} and f_{t-1} instead of
f_n.

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?

Thank you!

Cristina

Would it make sense to calcualte U_t by using f_{t+1}

On Wed, Dec 15, 2010 at 5:04 PM, brian avants <stnava at gmail.com> wrote:

> 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
> >
> >
>
>
>
> --
> ß®∫∆π
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20110126/5861506c/attachment-0001.htm>