I'm writing my own FEM Poisson equation solver in C++ and I'm using Mathematica as a benchmark - currently it seems that my program outputs a correct solution, however I want to check whether the laplacian of the result is equal to the right-hand side. Unfortunately, the interpolation itself raises a warning "Interpolation on unstructured grids is currently only supported for InterpolationOrder->1 or InterpolationOrder->All. Order will be reduced to 1". FEM deals with mostly unstructured grids, so as a result, the order of interpolation is now 1 (piecewise linear functions). It is no surprise that when I try to compute the laplacian of the result, I get flat zero function.

For the sake of simplicty, let's use the following code

ClearAll[x, y]
c = {}; n = 10000;
For[i = 0, i <= n, i++,
 p = {x -> Random[], y -> Random[]};
 AppendTo[c, {x, y, Sin[5 x] Sinh[5 y] + 1/6 x^3 - 1/12 y^4} /. p];

The following

f = Interpolation@c;
lapf = D[f[x, y], {x, 2}] + D[f[x, y], {y, 2}];
Plot3D[lapf, {x, 0, 1}, {y, 0, 1}]

gives flat zero function, although the laplacian of Sin[5 x] Sinh[5 y] + 1/6 x^3 - 1/12 y^4 is x - y^2. What can be done in case I have only numerical data on an unstructured grid and I need to evaluate laplacian at some point?

  • $\begingroup$ What about choosing a nice, smooth function $u$, compute the Laplacian $f = \Delta u$ of it analytically, sample it (so that you obtain a discrete function $f_h$) and compare the degrees of freedom of the FEM solution $u_h$ (i.e., $A_h \, u_h = M_h \, f_h$ with stiffness matrix $A_h$ and mass matrix $M_h$) to $u$ sampled on the grid? That's what I usually do to debug my solvers... $\endgroup$ – Henrik Schumacher Jan 7 '18 at 22:41
  • $\begingroup$ I see what you're suggesting...it seems like a feasible approach I might want to try. $\endgroup$ – user16320 Jan 7 '18 at 23:04
  • $\begingroup$ Don't forget to pick the right interpolation order before taking derivatives. $\endgroup$ – Vsevolod A. Mar 15 '18 at 14:21
  • $\begingroup$ Have you tried the Obtuse package?: familydahl.se/mathematica/#Obtuse $\endgroup$ – xzczd Jul 13 '18 at 15:35
  • $\begingroup$ If noise is a concern you may want to have a look at mathematica.stackexchange.com/a/10997/1089 and mathematica.stackexchange.com/a/176375/1089 $\endgroup$ – chris Jul 13 '18 at 19:33

First, instead of procedural code to produce your sample data, how about using functional programming on a smaller sample of 200 points:

data = Apply[
  Function[{x, y}, 
    {x, y, Sin[5 x] Sinh[5 y] + (1/6) x^3 - (1/12) y^4}
  ], RandomReal[{0, 1}, {200, 2}], 1]

Then, as the warning message suggests, try

f = Interpolation[data, InterpolationOrder -> All]

A plot of f shows that it agrees well with your original function. Next compute the Laplacian

lapf = Function[{x, y}, Evaluate[Laplacian[f[x, y], {x, y}]]]

which also agrees very well with the Laplacian of f.

| improve this answer | |
  • $\begingroup$ G'day, (1) f["InterpolationOrder"] gives {9223372036854775806, 9223372036854775806} which seems remarkable; is Mathematica trying to interpolate over the whole domain? (2) In Version 11.0.1, I get a blank plot for f and lapf, please see here. This does not happen if InterpolationOrder -> All is omitted. (3) If I try it with 2000 points, and attempt to evaluate say lapf[0.5,0.5], it gives the warning "The interpolation points {<>} are not poised, so an interpolating polynomial of total degree 62 could not be found". What is going on? Cheers! $\endgroup$ – user40265 Apr 11 '18 at 4:32

You could import your higher order mesh (2nd order is possible) and construct an ElementMesh from the data like so:

mesh = ToElementMesh[
   "Coordinates" -> {{0., 0.}, {1., 0.}, {1., 1.}, {0., 1.}, {0.5, 
      0.}, {1., 0.5}, {0.5, 0.5}, {0.5, 1.}, {0., 0.5}}, 
   "MeshElements" -> {TriangleElement[{{1, 2, 3, 5, 6, 7}, {3, 4, 1, 
        8, 9, 7}}]}];

Construct some example values.

values = Function[{x, y}, x*y] @@@ mesh["Coordinates"];

Create the interpolating function.

if = ElementMeshInterpolation[{mesh}, values];

Plot the result.

Plot3D[if[x, y], {x, y} \[Element] mesh]

enter image description here

However, be careful with computing second order derivatives for that; this is more of a test of the interpolating function that it is for the PDE solver you write.

To test the solver you use an somewhat arbitrary expression and insert that in the differential equation. That will give you a right hand side. Solve the PDE with this new right hand side and check how close the solution is to the somewhat arbitrary expression you used to create the right hand side. This process is called manufactured solution.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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