Assembly Process for FE (Finite Element) code in C

[Babak]

Hello everyone,

I'm working on some Finite Elements(FE) codes in C and now I
encountered some problems in assembly stage. The main idea is that a
large number of 3 by 3 elemental stiffness matrices are calculated and
then they should mapped to specific elements in a very big and sparse
global stiffness matrix. I think the preferred storage format for the
sparse matrix is Compresses Row Storage (CRS) format because it is
compatible with many different sparse solver packages. Now the problem
is that I have no idea how to update and assemble "ia" (pointers to the
begining of each row), "ja" (column indices) and "a" (nonzero values)
arrays as each local stiffness matrix is calculated in the loop. Is it
a better idea to first calculate all the elemental stiffness matrices
and then assemble them in the global system once at the end or I should
update these three arrays each time a local matrix is calculated? I
think SPARSE function in MATLAB can do much of these works for me but I
can't find any similar fuction in C. Since this is my first Finite
Element program, any comments, algorithms or references on this issue
will be a great help.

[Michael Mair]

Babak schrieb:

I'm working on some Finite Elements(FE) codes in C and now I
encountered some problems in assembly stage. The main idea is that a
large number of 3 by 3 elemental stiffness matrices are calculated and
then they should mapped to specific elements in a very big and sparse
global stiffness matrix. I think the preferred storage format for the
sparse matrix is Compresses Row Storage (CRS) format because it is
compatible with many different sparse solver packages. Now the problem
is that I have no idea how to update and assemble "ia" (pointers to the
begining of each row), "ja" (column indices) and "a" (nonzero values)
arrays as each local stiffness matrix is calculated in the loop. Is it
a better idea to first calculate all the elemental stiffness matrices
and then assemble them in the global system once at the end or I should
update these three arrays each time a local matrix is calculated? I
think SPARSE function in MATLAB can do much of these works for me but I
can't find any similar fuction in C. Since this is my first Finite
Element program, any comments, algorithms or references on this issue
will be a great help.

This is not exactly a C question; there are several questions.
1) It depends on the finite element formulation (which itself
depends on the application). Let us assume that you have only
local stiffness matrices; the usual way to progress is to
calculate the local stiffness matrices elementwise and add
them to the global stiffness matrix. This part can be easily
parallelized and most formulations offer themselves.
Decide whether you want to manipulate matrix rows and columns
and the right hand side for natural and essential (Neumann or
Dirichlet, or maybe even Robin) boundary conditions or whether
you want to eliminate these rows and columns as far as possible.
Depending on the boundary direction, parametrisation, and the
boundary conditions, assembly of the boundary conditions can
become quite nontrivial.
2) The way to store the matrix depends on the intended solver
for the linear system.
If you use a multigrid method or another iterative solver,
you should not go for a too dumb compressed row or column
format.
If you intend to use a standard direct solver (say LU), then
there are enough variants for the CRS.
3) Your question about the matrix format and the "update and
assembly" may be a C question if presented correctly.
4) There is at least one finite element toolbox written in C:
ug. I do not recommend it for someone who is not deeply
immersed in Numerical Mathematics and a very good C
programmer.

There probably are newsgroups and mailing lists which can
answer questions 1)-4) better than comp.lang.c.

Back to the C part of your question: Give us C code showing
how your variety of CRS is defined; then explain your problem
in terms of C.
However, people in another newsgroup or mailing list may be
able to point you to tested and working C code for your
purposes.
Cheers
Michael
--
E-Mail: Mine is an /at/ gmx /dot/ de address.

[Rouben Rostamian]

In article <11********************@e56g2000cwe.googlegroups.c om>,
Babak <b.*******@gmail.com> wrote:

I'm working on some Finite Elements(FE) codes in C and now I
encountered some problems in assembly stage. The main idea is that a
large number of 3 by 3 elemental stiffness matrices are calculated and
then they should mapped to specific elements in a very big and sparse
global stiffness matrix. I think the preferred storage format for the
sparse matrix is Compresses Row Storage (CRS) format because it is
compatible with many different sparse solver packages. Now the problem
is that I have no idea how to update and assemble "ia" (pointers to the
begining of each row), "ja" (column indices) and "a" (nonzero values)
arrays as each local stiffness matrix is calculated in the loop. Is it
a better idea to first calculate all the elemental stiffness matrices
and then assemble them in the global system once at the end or I should
update these three arrays each time a local matrix is calculated? I
think SPARSE function in MATLAB can do much of these works for me but I
can't find any similar fuction in C. Since this is my first Finite
Element program, any comments, algorithms or references on this issue
will be a great help.

Search the web for UMFPACK. This is an open source library,
written in standard C, for storing sparse matrices and solving
sparse linear systems. It is especially well-suited for finite
element computations. Matlab's sparse matrix functions are
based on UMFPACK.

Here is how you would assemble the global stiffness matrix.
Let nelems be the number of elements. You say each element
stiffness matrix is 3x3. Here is a sketch of what you will do:

int *Ti, *Tj;
double *Tx;
int i, j, k, elemno;

Ti = malloc(3 * 3 * nelems * sizeof *Ti);
Tj = malloc(3 * 3 * nelems * sizeof *Tj);
Tx = malloc(3 * 3 * nelems * sizeof *Tx);

k = 0;
for (elemno=0; elemno<nelems; elemno++)
for (i=0; i<3; i++)
for (j=0; j<3; j++) {
Ti[k] = I;
Tj[k] = J;
Tx[k] = m[i][j];
k++;
}

Here m[][] is the stiffness matrix of the element under
consideration and I and J are the row and column numbers in
the global stiffness matrix where m[i][j] is to be stored.

After computing Ti, Tj, Tx, you would call UMFPACK's
umfpack_di_triplet_to_col() function to convert the temporary
representation Ti, Tj, Tx of the global stiffness matrix to
UMFPACK's standard representation, denoted by Ap, Ai, Ax in the
manual. Then you can free the storage allocated to Ti, Tj, Tx.

Read UMFPACK's documentation for further details.

--
Rouben Rostamian

[Martin Joergensen]

"Babak" <b.*******@gmail.com> skrev i en meddelelse
news:11********************@e56g2000cwe.googlegrou ps.com...

Hello everyone,

I'm working on some Finite Elements(FE) codes in C and now I
encountered some problems in assembly stage. The main idea is
that a
large number of 3 by 3 elemental stiffness matrices are
calculated and
then they should mapped to specific elements in a very big and
sparse
global stiffness matrix. I think the preferred storage format
for the
sparse matrix is Compresses Row Storage (CRS) format because it
is
compatible with many different sparse solver packages. Now the
problem
is that I have no idea how to update and assemble "ia"
(pointers to the
begining of each row), "ja" (column indices) and "a" (nonzero
values)
arrays as each local stiffness matrix is calculated in the
loop. Is it

I don't know about finite element, but I am trying to make a
finite difference program work using something like:

for(i=2, count = 0, no = 0 ; i< nx; i++) {
for(j=2; j< ny; j++) {
if(unknowns[i][j] == 1) {
no++;
pfvector_a[++count] = koeffB[no];
pivector_j[count] = (i-2)*(nx-2) + (j-1);
/* middle cell */
pivector_i[count] = no;
}
if(unknowns[i][j-1] == 1) {
pfvector_a[++count] = koeffA[no];
pivector_j[count] = (i-2)*(nx-2) + (j-2);
/* left */
}
if(unknowns[i][j+1] == 1) {
pfvector_a[++count] = koeffC[no];
pivector_j[count] = (i-2)*(nx-2) + (j);
/* right */
}
if(unknowns[i-1][j] == 1) {
pfvector_a[++count] = koeffE[no];
pivector_j[count] = (i-3)*(nx-2) + (j-1);
/* row up */
}
if(unknowns[i+1][j] == 1) {
pfvector_a[++count] = koeffD[no];
pivector_j[count] = (i-1)*(nx-2) + (j-1);
/* row down */
}
}
}

That is a 2D problem so for each unknown equation in the matrix I
have to look left, right, up and down in order to take boundary
effects in consideration (my first and last cell doesn't have as
many coefficients as the interior cells - actually that applies
to all cells on the domain boundaries). In my code my mesh is
defined from (1,1) to (nx, ny) - or is it (ny, nx)? It's not
completely finished yet.....

Using this method, I have a 2D array called unknowns which allows
me to do calculation on arbitrarily shapes (circle, triangle, L
or T-shaped or even more random if I like) - the routine takes
care of inserting the correct coefficients so if we're standing
on the left top boundary ofcourse there's no coefficient above
and to the left for that equation...

However, I must also say that my code doesn't work at the
moment.... Don't know if there's anything I forgot to include,
since this is also my first time working with such a solver that
requires this format (am trying to make a conjugated gradient
solver work)....

You're saying that this format with "ia" (pointers to the
begining of each row), "ja" (column indices) and "a" is called
Compressed Row Storage (CRS) format ? I didn't knew that, so if
you make any progress I would be interested in seeing how you
solve the problem - (you can also mail me, just remove the .spam
and spam. before and after the @).
Best regards / Med venlig hilsen
Martin Jørgensen

--
---------------------------------------------------------------------------
Home of Martin Jørgensen - http://www.martinjoergensen.dk