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.

  • $\begingroup$ I'm assuming these matrices arise from a coupled system of PDEs of the form $Ku+Qp=f_u$, $Q^Tu+Sp =f_p$ (it would help to state the PDEs explicitly in the question). In this case, yes, it's just a block matrix (the coupling between $u_h$ and $p_h$ at each node is taken care of by the matrices $Q$ and $Q^T$). $\endgroup$ Commented Feb 18, 2016 at 7:32
  • 1
    $\begingroup$ Hey, described system is correct. Problem is that Taylor-Hood elements ares used, so u has 9 nodes, and p has 4 nodes. So each node has varying degree of freedom. I'm not sure how to tackle this as each node has varying degree of freedom $\endgroup$
    – pozza
    Commented Feb 18, 2016 at 15:05

3 Answers 3


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.


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/

  • $\begingroup$ I don't think that's what this question is a about: Calib knows how to assemble a finite element matrix for a single PDE -- it's specifically how to deal with coupled PDEs (and in particular, saddle point problems) that they have problems with. (TL;DR: How do you assemble the coupled stiffness matrix for the system if you have already assembled the stiffness matrices for each component PDE.) $\endgroup$ Commented Feb 18, 2016 at 18:38
  • $\begingroup$ @ChristianClason exactly. $\endgroup$
    – pozza
    Commented Feb 18, 2016 at 21:21
  • $\begingroup$ Well, I was just trying to explain one way assembly of elements with varying nodal degrees-of-freedom is handled in production FE software. Is the question about "append them together in big matrix" simply a conceptual one? Certainly you wouldn't implement it that way. Beyond that, I didn't see anything at all in the question about the fact that the resulting equations are not positive definite, i.e. a saddle point problem. $\endgroup$ Commented Feb 18, 2016 at 22:28
  • $\begingroup$ Your answer helped a lot thanks Bill. In class we never went over degree of freedom map. $\endgroup$
    – pozza
    Commented Feb 19, 2016 at 3:40

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.

  • $\begingroup$ Thanks, stalli. I did have a look at PETSc, however solution needs to be computed in MATLAB :( $\endgroup$
    – pozza
    Commented Feb 18, 2016 at 15:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.