I am trying to solve a steady-state advection-diffusion problem of the following form: $-D \nabla^2 T + \nabla \cdot (\mathbf{u} T) = -1$. I will be solving this problem in a 3D geometry (typically the surface of revolution). However, here I am phrasing the question for a disk.

The velocity field $\mathbf{u}$ is precomputed numerically from a different simulation (using an integral equation). I can compute the velocity at any given node point of the finite element mesh (with spectral accuracy). Is there any way to input $\mathbf{u}$ computed at all the mesh points (as a list) to NDSolve?

Adding more details (code)

\[CapitalOmega] = Ellipsoid[{0, 0, 0}, {5, 5, 5}];
mesh = #1 -> #2 &~
    MapThread~{{"Coordinates", "MeshElements", "BoundaryElements", 
     ToExpression@Import["mesh_sphere_rad=5.txt"]} // 
   ToElementMesh @@ # &;

velocity = Import["nodevel3D.txt", "Table"];
ux = ElementMeshInterpolation[{mesh}, velocity[[All, 1]]];
uy = ElementMeshInterpolation[{mesh}, velocity[[All, 2]]];
uz = ElementMeshInterpolation[{mesh}, velocity[[All, 3]]];

\[CapitalDelta][f_] := 
  D[f, {x, 2}] + D[f, {y, 2}] + D[f, {z, 2}]; (* Define Laplacian *)
source = -1; bval = 0; nval = 0;
Pe = 50; (* Define Peclet Number *)
sol = 
 NDSolveValue[{\[CapitalDelta][T[x, y, z]] + 
     Pe (ux[x, y, z] D[T[x, y, z], x] + 
        uy[x, y, z] D[T[x, y, z], y] + 
        uz[x, y, z] D[T[x, y, z], z]) == 
    source + NeumannValue[nval, \[CapitalGamma]N], 
   DirichletCondition[T[x, y, z] == bval, \[CapitalGamma]D]}, 
  T, {x, y, z} \[Element] \[CapitalOmega], AccuracyGoal -> 10, 
  Method -> {"FiniteElement", 
    "MeshOptions" -> {"MaxCellMeasure" -> 0.005, 
      "MeshQualityGoal" -> 1}}]

I notice that if I put {x,y,z} \[Element] mesh then their is no error thrown at me. However, if I use {x, y, z} \[Element] \[CapitalOmega] (the region on which I constructed my mesh for interpolating velocity at the first place) it says: ```InterpolatingFunction::femdmval: Input value {2.77554,3.67892,-1.91818} lies outside the range of data in the interpolating function.''' However the final solution from both these cases look almost identical. So I guess this is okay. But it will be good to know why this is happening. Thanks!

  • 1
    $\begingroup$ If there's no way to do it directly, you could probably do it indirectly via an interpolation. This might help: reference.wolfram.com/language/FEMDocumentation/ref/… $\endgroup$
    – Michael E2
    Jul 3, 2022 at 15:29
  • $\begingroup$ This may actually be helpful if I cannot do it directly. The only unfortunate part is that my velocity evaluation is spectrally accurate. However, the interpolation reduces the accuracy. It will be nice if the interpolant returns the exact values on the mesh points. $\endgroup$
    – bchakra
    Jul 3, 2022 at 16:49
  • 2
    $\begingroup$ The interpolation should not reduce the accuracy at all, if you interpolate over the same mesh as used in the FEM solution of your PDE. (NDSolve would not be evaluating the interpolation between the nodes, only at the nodes. Note that by default, NDSolve uses a quadratic mesh, and it would use values of $\bf u$ at the "quadratic" points. You would need these values to use $\bf u$ directly anyway. If you don't have them, you could use a linear mesh in NDSolve, but then you lose accuracy in the FEM solution.) $\endgroup$
    – Michael E2
    Jul 3, 2022 at 17:13
  • $\begingroup$ Let me see if I am understanding this: 1. I compute the node points: pts = mesh["Coordinates"]. 2. I evaluate precomputed $\mathbf{u}$ at these pts. 3. Use ElementMeshInterpolation to construct the interpolator. Right? Then only these values will be used and I won't lose accuracy. Is that correct? I will definitely try this! $\endgroup$
    – bchakra
    Jul 3, 2022 at 17:59
  • 1
    $\begingroup$ @MichaelE2, yes, that is correct that system matrices and solution vectors all correspond to the elements and nodes in the mesh. It's just that the evaluation of the PDE coefficients always takes place within an element, never at the nodes. Think of it as you'd do numerical integration over a rectangle, the integration points are inside the rectangle. It's during the multiplication with the integrated shape functions that these integrated coefficients get interpolated to the nodes of the elements. One of these days, I want to write a tutorial on how the FEM actually works.... $\endgroup$
    – user21
    Jul 8, 2022 at 8:17

1 Answer 1


The NDSolve when used with the finite element method never evaluates at node (mesh coordinate) points it always evaluates functions between these nodes at integration points.

The way to go about your data is to create an interpolating function with values from your function on the mesh nodes.

mesh = ToElementMesh[Disk[]];
values = Sqrt[Total[mesh["Coordinates"]^2, {2}]];
if = ElementMeshInterpolation[mesh, values]

When given to the FEM, that interpolating function if will not be evaluated at the nodes, only at the integration points. Using an interpolating function that is based on the same mesh as the FEM mesh will be reasonably fast.

One other (slower) approach would be to take the mesh, find the integration points, (like shown in the other post), evaluate your function at the integration points, create a new Delaunay mesh of the integration points and use that for the interpolating function. This will be a bit slower, but then you'd get the evaluation at the integration points use your values.

One thing to keep in mind is that FEM is a machine precision numeric method; I am not sure how much it is going to help if your integration values are supper accurate compared to the accuracy of the FEM.

For your specific problem look at the examples in the heat transfer PDE models section. Specifically, the Thermal Decomposition model has a squential PDE solution; the the solution of the Navier-Stokes equation is returned as set of interpolating functions that are an input to the convection term in the heat equation.

  • $\begingroup$ Thanks a lot for the clarification. I looked into the heat-transfer PDE section. It is not clear from the module how I can input numerical data. They have examples where the flow is solved along with the heat transfer problem (which is different from my question). Can I access the integration points of a 3D mesh and construct the interpolation function from velocity evaluations at those points? Will this be more accurate? And how to access the integration points in a 3D mesh? Thanks a lot for your help! $\endgroup$
    – bchakra
    Jul 7, 2022 at 14:12
  • $\begingroup$ @bchakra, see update. $\endgroup$
    – user21
    Jul 7, 2022 at 14:41
  • $\begingroup$ Thank you! One thing that I noticed with the interpolation function used with NDSolve is that I am getting the following message: "Input value {-4.25795,0.634379,-2.34425} lies outside the range of data in the interpolating function". Why will this happen when the interpolation is constructed on the mesh? $\endgroup$
    – bchakra
    Jul 7, 2022 at 16:14
  • $\begingroup$ Are you sure that the meshes are the same? That should not happen if the meshes are the same. To change the behavior of the interpolating function you can use the ExtrapolationHandler also look in the ref page of ElementMeshInterpolation but you first should make sure that the meshes are the same. $\endgroup$
    – user21
    Jul 8, 2022 at 8:25
  • 3
    $\begingroup$ If you want better answers you need to provide the code; everything else is guess work. $\endgroup$
    – user21
    Jul 8, 2022 at 8:26

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.