Projecting Finite Element solution onto new mesh

Question

I am implementing a finite element solver in MATLAB and I have the following problem. Let's say I have a mesh $\mathcal{T}_1$ with triangular elements on a rectangular domain $\Omega\subset\mathbb{R}^2$. I know the positions of all the nodes on the mesh, and I have an approximate solution $\hat{u}_1(x,y)$ to my PDE on this mesh (this approximation is stored as a vector of real numbers for each node, and then interpolated between nodes). To be exact, I use piecewise affine basis functions, like seen in this illustration.

If I now create a new mesh $\mathcal{T}_2$ (which may or may not be more coarse/fine, and positions of nodes may be changed), is there a way I can "project" the data of $\hat{u}_1(x,y)$ onto the new mesh $\mathcal{T}_2$, call it $\hat{u}_2(x,y)$, so that $\hat{u}_1(x,y)\approx \hat{u}_2(x,y)$? In short I want $u_2(x,y)$ to be the "best representation" of $u_1(x,y)$ at the new mesh.

As an example, let's say there is a disk of radius $1$ around the origin, $B_1(0,0)\subset\Omega$, and perhaps in the initial mesh $\mathcal{T}_1$ the nodes $v_3^{(1)},v_5^{(1)},v_6^{(1)},$ and $v_8^{(1)}$ are contained in this disk. However in the new mesh $\mathcal{T}_2$, perhaps a coarser mesh, there are only $2$ nodes contained in the disk, say $v_2^{(2)}$ and $v_3^{(2)}$. I now want a solution that interpolates in the nodes, and is approximately equal to the solution on the old mesh.

I understand that if the new mesh is much more coarse, of course the solutions can't be that similar, but I just want a decent approximation. If anyone knows any existing MATLAB code to do this, or anything to help me implement it myself, I'll be happy.

It's a bit difficult to explain exactly what I mean in just words, so if it is not clear at all what I'm asking, please say so, and I'll try to clarify.

Thanks

Answer 1

When you want to say that you want $u_2$ to be the best approximation of $u_1$ on mesh ${\cal T}_2$, then you have to define what you mean by "best". Let's assume you define it as that function that has the least $L_2$ error, i.e., $$ u_2 = \arg\min_{\phi_h \in V_2} \|\phi_h-u_1\|_{L_2(\Omega)}. $$ Then indeed $u_2$ is the $L_2$ projection of $u_1$ onto the finite element space $V_2$ defined on ${\cal T}_2$, and it satisfies the condition $$ (\phi_i,u_2) = (\phi_i,u_1), \qquad i=1,\ldots,\textrm{dim}(V_2) $$ where $\phi_i$ are the shape functions defined on mesh ${\cal T}_2$.

If you want to compute this, you expand $u_2(x,y)=\sum_j U_{2,j} \phi_j(x,y)$: $$ \sum_j (\phi_i,\phi_j) U_{2,j} = (\phi_i,u_1), \qquad i=1,\ldots,\textrm{dim}(V_2), $$ so you have to solve an equation of the form $MU_2 = F$ where $M$ is the mass matrix on the ${\cal T}_2$ and $$ F_i = (\phi_i,u_1) = \int_\Omega \phi_i(x,y)u_1(x,y) \; dx \; dy. $$ The problematic part is the evaluation of this integral. Of course, you will do it using quadrature, but that will require you to evaluate $u_1$ at the quadrature points defined on mesh ${\cal T}_2$. If ${\cal T}_1$ is a uniform mesh, then evaluating $u_1$ at arbitrary points is not difficult, but if it is an unstructured mesh, then this is in general a difficult and expensive problem.

Answer 2

Since you are using MATLAB, I suggest you take a look at the scatteredInterpolant class. You construct the interpolant by giving it the x-y locations and solution values of points on the first mesh. Then you ask for interpolated values of the solution at points in the new mesh (e.g. new element integration points).

This will not provide a "best" interpolation in the minumum-L2-error sense as Wolfgang Bangerth has described. It simply does a local linear interpolation based on points in the original mesh. However it is easy to use, quite fast, and is often satisfactory.