cm-cmgui

Project Outline

Develop a test application which solves a static Laplace (e.g., steady state heat distribution) model over a two-dimensional domain with bilinear element interpolation using the CMGUI finite element data structures to store field information. The application will be written in C using the C API to CMGUI.

Project Goals

• Test the computational performance of using the CMGUI data structures/API versus the direct access arrays currently used in cm.
• Identify any deficiencies in the current CMGUI API when using the finite element data structures in this application.

Methods

Found my old afem project code (written in C) and got that working again for static 2D Laplace. Using CMISS examples 314(1,2,3,4) this code was able to obtain identical solutions as current cm. This code formed the basis of the test application.

The first step was to replace the Node and Element structures defined in my afem code with accessing field information from the CMGUI data structures. Access is primarily through the function_finite_element interface. The only field of interest in this example application is the coordinates field, or more specifically the spatial derivatives with respect to local xi at the Gaussian quadrature points.

Using the derivatives from the coordinate field at Gauss points (ds/dxi) you can invert to get dxi/ds and assemble the element stiffness matrix. The element stiffness matrix can then be added into the global stiffness matrix where we are still using node numbers/ids as indexes into the global matrix.

Once we have the global stiffness matrix, we read in and apply some boundary conditions (using the dodgy file format from my afem project) and then solve the system. This results in a vector u containing the value solution at each node and a vector fluxes containing the weighted nodal flux solutions. The solutions are then put into new fields defined at the nodes (u and fluxes) and interpolated over the elements (u only) and the full mesh written out with all fields.

Notes, Problems, and Got-ya's

Put here all the things that I find out along the way.

Node IDs as matrix/vector indices

We are currently using node IDs (numbers) as the row/column number for a given node in the global stiffness matrix and solution vectors. This is a bad thing!

A better method is to define some kind of mapping between node IDs and matrix/vector indices. Need to work out a good way to do this in a general way while maintaining efficency (speed/memory).

Boundary Conditions

So far just using a basic struct for the storage of boundary condition information and reading in from a file in the format I made up for my afem project (Have a script to generate file from the CMISS ipinit files for essential BCs only). Ideally the BCs would come from additional fields in the mesh?

Node field creator

There is some kind of problem using the same node field creator when defining multiple fields over the same elements? Need to get shane to describe the actual problem.

We got around the problem by only defining the solution value field over the elements but having both value and flux solutions at all the nodes. This seems fine for this application as it probably makes little sense interpolating nodal flux values over the elements anyway.

Global nodes for a specific element

There is currently no method in the API for accessing global node IDs for a specific field in a given element. The only method is 'calculate_FE_element_field_nodes(...)' in finite_element/finite_element.h. We need the nodes for a specific field in the element as any element may contain numerous nodes for different fields, interpolations, collapses, etc.

'calculate_FE_element_field_nodes(...)' needs a FE_field pointer to get the field it works on which is a data type hidden by the API. Luckily you can pass in a NULL FE_field and the method assumes a default of the coordinates field for that element. So this method works for us now.

The solution is to bring this method up into the API level but using function_finite_element's instead of FE_field's. The existing function is also very complex and would probably want a bit of simplification for the mode simple purpose that it would serve at the API level than what it currently does in the bowels of CMGUI.

Destroying elements

There is a method 'destroy_Cmiss_node(Cmiss_node_id *node_id_address)' in api/cmiss_finite_element.h which is simply a wrapper around the deaccessing of the actual FE_node. There is no corresponding 'destroy_Cmiss_element(Cmiss_element_id *element_id_address)' function - is there a reason for this or is it just missing??

Element dimension

The only method I can find in the API to get an element's dimension is to create a function_element for the element ('Cmiss_function_element_create(element)') and then use 'Cmiss_function_element_dimension(functionElement,&dimension)' to get the dimension. Seems like a lot of work to only get the dimension. Perhaps there could be a 'int Cmiss_element_get_dimension(Cmiss_element_id element)' wrapper function around 'int get_FE_element_dimension(FE_element *element)' added to api/cmiss_finite_element.h?

Basis function derivatives

Would be good to be able to get the derivatives of an field's basis function for a given element/xi location.

Derivative matrix

In assembling the element stiffness matrix for this model, we need the derivatives of xi with respect to global space (dxi(i)/ds(j)) at a given element/xi location. I have currently implemented this by first using the coordinates function_finite_element and the 'Cmiss_function_variable_evaluate_derivative(..)' function to independently get ds(i)/dxi(1) and ds(i)/dxi(2) and combine them into a 2x2 matrix. This matrix is then manually inverted to give the dxi(i)/ds(j) matrix.

Potentially this could be simplified and generalised by getting the full ds(i)/dxi(j) matrix from the evaluate_derivative function and using the function_matrix invert method. Currently unsure how to do this or if it is even possible. Shane thought this was getting a bit close to what David is currently working on, so since its just a simple 2d model that we'd be better off with the above hard-coded approach.

Scalar data type

Need to check the simple casting between Scalar and double. Or maybe the application should be using the Scalar data type and the API can provide things like the format specifier and functions for converting to double/float types?

Types of funnctions returned by evaluate and evaluate_derivative

It wasn't at all clear from the comments in the code how to go about getting actual numbers from the function returned from the evaluate and evaluate_derivative methods. The only method close is 'get_string_representation(..)' but thats not terribly good when you want to use the values numerically.

Shane provided the break through when he told me that you get back a function_matrix and that you could use the 'Cmiss_function_matrix_value(..)' method to get the values out. Once I knew this I found the 'Cmiss_function_get_type_id_string(..)' method which I'm guessing could be used to find out the type of a function - but I can't seem to find anything that tells me the possible values of the returned string.

I guess this falls into the documentation basket :-)

Project Outcomes

Seems to be working fine and giving the same answers as CMISS examples 314(1,2,3,4) and 312. Have yet to determine the best method for performance comparisons with back-end CMISS.

Code

Code now available "here":/opencm/wiki/2005-06.tar.gz/file_view