Finite Element NDSolve with input data at node points

Question

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", 
      "PointElements"}, 
     ToExpression@Import["mesh_sphere_rad=5.txt"]} // 
   ToElementMesh @@ # &;
Dimensions[mesh["Coordinates"]]
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!

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.

Needs["NDSolve`FEM`"]
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.