[Insight-developers] check-in the FEM classes in itk

[Will Schroeder]

Hi Aljaz-

At 06:05 PM 11/21/2001 -0500, you wrote:

>Currently we use vnl_matrix and vnl_vector (non sparse version) for
>everything including the master stiffness matrix. Changing the master
>stiffness matrix and solver to something else has minimal impact on the
>code, because all big matrices are only present in the Solver class.
>
>The element stiffness matrix, which is returned by Element classes is
>normally small enough that can be handled by full matrix. Changing this one
>would require a bit more work, but still nothing that couldn't by done in a
>couple of hours.

Sounds good.


>> Don't worry about the regions for now. I assume that you are using
>itk::Mesh ?
>
>Actually, this is something I would like to talk about. Currently we're
>converting the input Mesh into arrays of other (non mesh) objects. The
>reasons are explained below.
>
>The whole system is currently defined only by a couple of simple arrays that
>can hold polymorphic objects --- Nodes, Elements, Loads, Boundary conditions
>and Materials. As such everything is completely independent of the Mesh
>class.

FORTRAN style has its advantages, meaning that I like the simplicity of arrays :-)


>All of the above classes are polymorphic, so we can access the basic
>required functionality of each class from a pointer to a base class (Element
>classes for example provide virtual function Ke(), which returns element's
>stiffness matrix).
>
>Linking between objects is provided by pointers (elements point to nodes for
>example). In general, we allow mixing of different types of elements and
>nodes in a single mesh. If everything is defined using pointers, we can
>enforce mesh consistency (only compatible elements can use a specific node
>class) using C++ type compatibility rules. The Mesh class however uses ID's
>and all type info is gone.

Pointers are difficult to work with when you want to save or transmit data (e.g., MPI). ID's are a fairly natural way for people to interface with a mesh, e.g., element id 100 or node 1000, etc. And the mesh design allows the ids to be a variety of types based on the PointIdentifier typedef.


>Another thing is that currently the Mesh class enforces the same PixelType
>throughout the mesh. If a PixelType is to be used to store information about
>degrees of freedom at each node, we really need to store a pointer to the
>base Node class (assuming that mixing of nodes and elements is still
>allowed). Object of specific derived Node classes will then store all the
>required info.

The pixel type could easily be a class with pointers, etc. to other information/types. That way the mesh can actually represent a variety of PixelTypes, or course at the cost of extra memory.


>If we used the Mesh class from the start, we would loose all of the above
>functionality. However, if a system is defined using Mesh class within itk,
>it can be very easily converted to our internal FEM representation.
>
>If a mesh could be expanded to provide the above functionality, that would
>be perfect. But since this expansion doesn't seem to be trivial (at least to
>my knowledge), we decided to use the conversion for now. The conversion,
>BTW, is just a simple 'for' loop that converts ID's to corresponding
>pointers and converts Cells to Elements and PixelType to Nodes.

In the short term it seems that the conversion is necessary. After the Beta, we may want to look at it again; either to modify itk::Mesh and/or modify the FEM code.

Will