Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

A lot of adaptive FEM libraries use more advanced mesh data structures to handle adding/removing nodes, edges, triangles, tetrahedra, etc. For example, the p4est library uses octree data structures for adaptive mesh refinement; you wouldn't often find octrees used for computations on a static mesh.

What changes on the linear algebra side for adaptive FEM?

The most blunt way I can conceive of would be to completely rebuild all the system matrices whenever the mesh is refined or coarsened. If mesh adaptation is a sufficiently infrequent operation, then the expense of doing so is ultimately amortized over the rest of the computation. One could easily leverage existing sparse linear algebra software (PETSc, Trilinos, etc.) with this approach.

Is this blunt method the most commonly used, or are there libraries that manage to reuse or modify the old matrix during refinement? After all, most of the mesh and the corresponding matrices are unchanged during a mesh adaptation.

share|improve this question
up vote 6 down vote accepted

Yes, the most common approach is to rebuild. Data structures that are modifiable in-place tend to be less efficient once set up, and reallocation is actually quite cheap compared to reassembly (e.g., due to nonlinearity) so it's really a fine solution. Outside of relatively rare niches with very easy solves, attempts to use dynamic data structures in the solvers will only make your application slower. It is a common perception among those who forget to measure or model performance, however.

share|improve this answer

Like Jed already said, reuse of linear algebra components such as matrices and vectors is not commonly done. It is also not necessary: setting up these components is comparably very cheap relative to the cost of solving linear systems.

If you look for things that change when going from static to adaptively refined meshes, then the biggest obstacle is dealing with hanging nodes. In deal.II, this is handled by the ConstraintMatrix class, which takes several 1000 lines of code. You can find some of the description of what this class does in the paper by myself and Oliver Kayser-Herold (linked to from my publications page). No other component (other than handling the mesh, of course) has required so much adaptation when going from fixed to adaptive meshes.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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