# Multigrid on “not perfectly rectangular” grid

Multigrid introductions normally use a rectangular grid. Interpolation of values is then straight forward: Just interpolate linearly on the edge between two adjacent nodes of the coarse grid to find the value of the fine grid node on that edge.

For an FEM application I have a grid which is "topologically" rectangular so that node connections are as on a rectangular grid. However, the nodes are not perfectly aligned on a grid but may travel small distances in order to better fit the geometry, while still maintaing connections as in a perfect rectangular grid.

The mesh looks something like that: Mesh example. You see: Connections are "regular rectangular", but node positions not.

I can image several "reasonable" geometric interpolation schemes für such a setting.

The general question is: Does multigrid require a perfectly aligned rectangular grid, or will it also work with the situation described above, as long as the interpolation is "good"? Or is it better to use algebraic multigrid in that case? (Which I don't prefer since it's not that intuitive as geometric multigrid.)

• Im not sure I understand what you mean by the grid being topologically rectangular, but where the nodes are not aligned as in a rectangular grid. Is the grid a structured rectangular grid? Could you please clarify this, perhaps with a drawing? Is the issue that you are using rectangular elements that do not line up as would happen when using a structured cartesian grid? – James Aug 23 '15 at 23:51
• @James, I interpreted OP's question to be more like: "What happens if I have a 'cartesian-ish' grid inside a trapezoid?". – Bill Barth Aug 24 '15 at 15:47

Multigrid doesn't need a Cartesian (rectangular), uniform grid. What it needs is that you can define a fine and a coarse level (possibly recursively, if you want to go from a two-level to a multi-level scheme), and that you can define interpolation operators between these levels. The easiest way to explain this is if you indeed have a Cartesian grid, but you can in fact start with any coarse mesh, refine it once, and just like this: you have a finer mesh.

In other words, it is easiest to think of multigrid-suitable meshes not as a fine mesh and how to find the coarser meshes, but to start with a coarse mesh from which you get the finer levels by uniform refinement (i.e., every quadrilateral is subdivided into four smaller ones). Since uniform refinement is always possible, this readily gives you a hierarchy. This is opposed to coarsening, which is not always possible if you are just given an particular mesh, and so makes defining a mesh hierarchy much more complicated. (Which is why people have come up with algebraic multigrid methods to define the coarse levels based just on the matrix, without thinking of the underlying mesh from which it was created.)

• I edited an image example into the question. I would assume coarsening is extremely easy in my case by just taking only every 2nd node (per dimension) and recompute FEM for that grid to get the coarse matrix/operator. The coarse data vector can be obtained by just taking every 2nd node value. – Michael Aug 25 '15 at 8:00
• That may or may not work as well as you want (I suppose it probably will, but am not completely sure). The thing is that in your example, the coarse mesh cells do not cover the same area as the children, and so there is no nesting property of the finite element spaces: the functions that you can represent on the coarse mesh are not a subset of the functions you can represent on the fine mesh. Joe Pasciak and Jim Bramble have written papers on such cases, and I think I recall that they can be made to work. But it's not immediately obvious to me. – Wolfgang Bangerth Aug 26 '15 at 11:47

Lets say that you have the following grid composed of rectangular elements:

Now if you perform your interpolation assuming a normal structured rectangular grid then you will be introducing errors associated with this inaccurate interpolation. In other words when you restrict your residual vector and when you prolong your error vector there will be errors from the interpolation.

Now if your grid is "close" to being a normal structured cartesian grid then this may work, at least at first, but I suspect one of two things will happen depending on how far off you grid is from being rectangular:

1) You might find that the multigrid begins to converge at first. After all initially your error is large anyways and your "approximate" interpolation really just means that some nodes are slightly over represented while some are slightly under respresented. However you may find that the convergence stagnates as the solution becomes more accurate and the interpolation errors become more important.

2)Another possibility is that the multigrid does end up converging, but not as fast as it should if you had used the correct interpolation.

Basically by being off with your interpolation you are weighting the importance of certain nodes inaccurately. For example in 2D if you are weighting a given node as:

\begin{bmatrix} 0.25 & 0.5 & 0.25 \\ 0.5 & 1.0 & 0.5 \\ 0.25 & 0.5 & 0.25 \end{bmatrix}

when in truth because your grid is not exactly cartesian it should be:

\begin{bmatrix} 0.25 & 0.55 & 0.25 \\ 0.55 & 1.0 & 0.49 \\ 0.28 & 0.52 & 0.30 \end{bmatrix}

then this will result in some error. Whether this error prevents convegence will likely depend on how far off your grid is from being cartesian.

While AMG is more difficult to understand/implement it sounds like it is the correct method for your grid. Applying geometric multigrid to an "approximate" rectangular grid may work, but I would guess that it is a band-aid solution at best. Hope this helps.

Update: I think there may have been some confusion in my answer. I am not saying that geometric multigrid will only work with cartesian meshes, but rather that defining interpolation (and hence restriction) on cartesian meshes is easy whereas on non-structured meshes this may be difficult. For example consider the case of even a simple 2D domain with a triangular mesh. Refining this mesh is easy - at least conceptually - but how would you define an interpolation operator between the coarse and fine mesh? I prefer AMG simply because it performs more like a "black box" solver, i.e. doesnt need information on the underying mesh, however this is just my person bias/quirk. Geometric multigrid can work as long as you can provide accurate interpolation operators.

• I can read this as "depends how good your interpolation is", correct? So If I can come up with a "perfect" interpolation scheme, then i should be fine? – Michael Aug 25 '15 at 8:02
• I think this answer is true but misleading -- it says that if you pretend your grid is regular cartesian when it's not, you will get the wrong answer. That's true, but you could say the same thing about any numerical method. The point is that geometric multigrid is accurate on more general meshes when implemented correctly. – David Ketcheson Aug 25 '15 at 12:01
• I agree that geometric multigrid can be used with non-cartesian meshes as long as done correctly. I just took the OP question as asking whether using an "approximately" cartesian mesh, while using interpolation designed for an actual cartesian mesh, would work. – James Aug 26 '15 at 1:51
• @Michael Basically you need some way of defining interpolation operators between different grid levels. With cartesian meshes this is easy. With non-cartesian meshes this can quickly become more difficult depending on how unstructured your meshes are. As Wolfgang says in his answer you can always create a more refined mesh from a coarse one, but defining meaningful interpolation operators may be difficult. The advantage of AMG is that it behaves more like a "black box" solver in that you do not need mesh data to define interpolation operators. All you need is the matrix – James Aug 26 '15 at 2:01
• @MIchael So to answer your question. Yes, if you can come up with accurate interpolation operators then geometric multigrid will work fine. Hope this helps. – James Aug 26 '15 at 2:05