8
$\begingroup$

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

Improve this question
$\endgroup$
6
  • $\begingroup$ If you're using affine (linear) elements, the unknowns are simply values of the solution at vertices - you could then evaluate the solution on the fine mesh at the coarse mesh vertices to get a representation on the coarse mesh. Alternatively, you could project the fine-mesh solution to the coarse mesh solution, but this may be much more costly. $\endgroup$
    – Jesse Chan
    Jan 15, 2015 at 5:47
  • $\begingroup$ @JesseChan I only have a solution on the initial mesh. Basically, I want my solution on the initial mesh "best represented" on a new mesh. The vertices do not necessarily have the same positions on the two meshes, so I don't understand how I could evaluate the solution of one mesh at another. Maybe I'm misunderstanding you? Also one mesh isn't necessarily coarse and one fine, one mesh might be finer in some areas, but coarser in some areas. $\endgroup$
    – Eff
    Jan 15, 2015 at 5:56
  • $\begingroup$ You had mentioned you interpolate between vertices, which lets you evaluate the solution on the fine mesh at the coarse mesh vertex positions? All you need for the solution on the coarse mesh is the value of the solution at vertex positions - though, as you said yourself, this may not give a very good representation. I can post a longer answer in more detail if that would help? $\endgroup$
    – Jesse Chan
    Jan 15, 2015 at 6:00
  • $\begingroup$ @JesseChan Could you give a hint on exactly how you would implement this? I could only come up with the following: For each vertex $v$ on the new mesh. Check if it was a vertex on the old mesh, if so, let their $u$-value be equal. If $v$ is inside an element on the old mesh, let the $u$-value be average (?) of that element's node values. It seems like there should be a better/easier implementation? $\endgroup$
    – Eff
    Jan 15, 2015 at 6:32
  • $\begingroup$ The key question here if you can make some assumption about the node positions of the new mesh with respect to the old one, or if both meshes are to be assumed as totally independent. Second you should define what a "good approximation" means to you... in which norm $\lvert \hat{u}_1 - \hat{u}_2 \rvert$ should be small? The solution proposed by Jesse simply means collocating the old solution at the nodes of the new mesh. (And may be the more sensible one, if no other info is available). $\endgroup$
    – Stefano M
    Jan 15, 2015 at 7:34

2 Answers 2

Reset to default
8
$\begingroup$

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.

Improve this answer
$\endgroup$
7
  • $\begingroup$ (+1) Thanks a lot for the answer, I will see if I can make this general implementation. Do you know what happens if we simplify the problem to the following: Let's say that, in the mesh $\mathcal{T}_2$ that we are to project onto, that all vertices also are vertices of the initial mesh $\mathcal{T}_1$ (However, there would be other vertices in the initial mesh that would not be in $\mathcal{T}_2$). Would the solution that minimizes the norm simply be given by the series expansion with coefficients of each vertex equal to the initial solution's coefficients at those vertices? If you understand. $\endgroup$
    – Eff
    Jan 15, 2015 at 17:22
  • $\begingroup$ Professor Bangerth: Is this the approach you use in deal.II? I know that deal.II supports h-adaptivity and nonlinear problems so I'm assuming this type of interpolation would be needed in that context. Thanks. $\endgroup$
    – Bill Greene
    Jan 15, 2015 at 18:12
  • $\begingroup$ @Eff: No. What you are asking is whether the interpolation (which is what you get if you just take nodal values) is the same as the projection (which is what I describe and you asked about in the original post). The answer is no. They're, in general, not the same. I also have that $\|Pu_1-u_1\| \le \|Iu_1-u_1\|$ where $Pu_1$ is the projection and $Iu_1$ is the interpolation of $u_1$ onto mesh 2. So the projection's error is better. On the other hand, I can show that $\|Pu_1-Iu_1\|\le Ch^2$, i.e., they are in a sense "close". $\endgroup$
    – Wolfgang Bangerth
    Jan 16, 2015 at 14:17
  • $\begingroup$ @BillGreene: In deal.II, we typically try very hard to avoid this kind of operation. Interpolating between unrelated meshes is very expensive because you need to evaluate $u_1$ at the quadrature points of mesh ${\cal T}_2$ and that implies searching which cell of ${\cal T}_1$ a particular quadrature point lies in. This is expensive. You also lose a significant amount of accuracy if you do it frequently. So what we try very hard is to ensure that the meshes ${\cal T}_{1,2}$ are related: e.g,. that they each result by refinement from a common mesh (i.e., they are hierarchically related). $\endgroup$
    – Wolfgang Bangerth
    Jan 16, 2015 at 14:20
  • $\begingroup$ @WolfgangBangerth Are you aware of a good reference that explains the computational difficulty of evaluating a finite element function at a large number of arbitrary points? Clearly this is true for the reasons you have explained, and I have observed it myself. But I would like a reference to cite in a paper I am writing. $\endgroup$
    – Nick Alger
    Jul 28, 2021 at 20:18
2
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ (+1) Thanks a lot. I will look into it, and see if it works "well enough." Could definitely save some work! $\endgroup$
    – Eff
    Jan 15, 2015 at 17:23

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.