How to assemble global matrix (for a coupled) problem?

Question

I'm trying to assemble global matrices for the following system.

$$ \begin{bmatrix} K& Q\\ Q^T&S \end{bmatrix} \begin{bmatrix} u_h\\ p_h \end{bmatrix} = \begin{bmatrix} f_u\\ f_p \end{bmatrix}$$

In which, $$ K, Q, Q^T, S $$

are matrix generated from element integral, which I know how to compute. Likewise, $$f_u, f_p$$ are vectors generated from the integral over the element (right-hand side), which I'm able to compute.

The question that I have is how to compute the global matrix? I know how to compute if it was a simple $KU = F$ system. I'm not sure how to go about it in a coupled case.

So far, I have matrices all elements, but I don't know how to assemble a global matrix system. Can I just create a global matrix for $K$, $Q$, $S$, $u_h$, $p_h$, and then assume they are block matrices and append them together in a big matrix? My problem is that, in a global matrix, each node won't be coupled.

The problem is that Taylor-Hood elements are used, so $u$ has $9$ nodes, and $p$ has $4$ nodes. So each node has a varying # degrees of freedom, either 2 or 3.

Answer 1

I think you may be viewing the problem the wrong way round. The rows of your matrix correspond to your discretized equations, which depend on the degrees of freedom in the problem, not the physical locations in space. You are free (theoretically) to number the $u$ nodes and $p$ nodes however you please and you will still get $K$, a 9x9 matrix $Q$, a 9x4 matrix and $S$ a 4x4 matrix, specific to your numbering (and choice of element). The spatial coupling comes through the integrals which define these matrices, you shouldn't thing that there is anything special about the face that some nodes are spatially co-located. Indeed there are discontinuous finite element methods in which individual variables can have multiple nodes located at a single point in space.

Answer 2

If I understand your question correctly it is primarily about programming rather than a conceptual FE one.

In general purpose FE programs, it is almost always assumed that different nodes have different sets of degrees of freedom (dof). In order to assemble each element into the global matrix, a key piece of information you need is a "map" from each degree of freedom in the element matrix to the appropriate equation in the global matrix. Here is a very simple way to create that map:

1. Create an array (the map) with dimensions max_dofs_per_node x number_of_nodes in the model. (In your case 3 x n)
2. For each node in the model, set the entries in the column of the map matrix corresponding to the node to 1 if the dof is active at that node and 0 otherwise.
3. Now iterate over all entries in the map matrix and assign a global equation number if the entry is 1. Start the equation numbers at one and increment each time.

During assembly, for each element node, you can easily determine the global equation number by referring to the appropriate column in the map matrix.

The fastest way to create a sparse matrix in MATLAB is by first creating a matrix of "triplets" (row, column, value) where row and column are the appropriate equation numbers from the map matrix and then passing the triplets to the MATLAB sparse function. That approach is discussed in this note by Davis: http://blogs.mathworks.com/loren/2007/03/01/creating-sparse-finite-element-matrices-in-matlab/

Answer 3

I would highly recommend using PETSc, especially if you want to assemble in parallel, and if you are planning to use preconditioners created by combining preconditioners for individual fields.

Assembly in PETSc is as simple as it gets. Look at some of the examples (many are for coupled problems as well). Parallel assembly and solve are the most complicated part of an FE code and PETSc takes care of both.