I need to solve a PDE using FEM but it involves a function known only on the points of the mesh. I mean, I do not have a formula for such function. I want to the delay the evaluation for that function until a numeric value is need it.

Is it possible?


1 Answer 1


Yes, it is possible and it is also very reasonable to do so (for example in iterative schemes for nonlinear equations - building an interpolating function in between just wastes a lot of time).

One can use those facilities that NDSolve uses as backend to assemble mass and stiffness matrix and create a load vector by mass.f, where f is the vector of function values on vertices (this somewhat assumes that f represents a more or less smooth function). Afterwards, you have to deploy the boundary conditions into this load vector.

Here is a simple example, essential copied from here.

Creating a mesh:

R = DiscretizeRegion[
    t \[Function] (2 + Cos[5 t])/3 {Cos[t], Sin[t]}, 
    Most@Subdivide[0., 2. Pi, 2000]],
   Line[Partition[Range[2000], 2, 1, 1]]
  MaxCellMeasure -> 0.00025,
  MeshQualityGoal -> "Maximal"

enter image description here

Suppose we'd like to solve a Poisson problem with homogeneous boundary conditions:

$$\left\{\begin{array}{rcll} - \Delta u &= &f, &\text{in $\varOmega$},\\ u_{\partial \varOmega} &= &0. \end{array}\right.$$

Setting up a first order mesh, the pde, and discretize it.

Rdiscr = ToElementMesh[R, "MeshOrder" -> 1];
vd = NDSolve`VariableData[{"DependentVariables", 
     "Space"} -> {{u}, {x, y}}];
sd = NDSolve`SolutionData[{"Space"} -> {Rdiscr}];
cdata = InitializePDECoefficients[vd, sd,
   "DiffusionCoefficients" -> {{-IdentityMatrix[2]}},
   "MassCoefficients" -> {{1}}
bcdata = InitializeBoundaryConditions[vd, sd, {DirichletCondition[u[x, y] == 0., True]}];
mdata = InitializePDEMethodData[vd, sd];
dpde = DiscretizePDE[cdata, mdata, sd];
dbc = DiscretizeBoundaryConditions[bcdata, mdata, sd];
{load, stiffness, damping, mass} = dpde["All"];

Now, we can take a random vector f and create its associated load vector:

f = RandomVariate[NormalDistribution[], Length[mass]];
load = mass.Partition[f, 1];
DeployBoundaryConditions[{load, stiffness}, dbc];

Solving the equation with LinearSolve.

solution = LinearSolve[stiffness, load, Method -> "Pardiso"];

And plotting the result:

solfun = ElementMeshInterpolation[{Rdiscr}, solution];
Plot3D[solfun[x, y], {x, y} \[Element] R]

enter image description here

  • $\begingroup$ The random numbers defined on grid points are not compatible with the grid spacing requirements for a meaningful solution though. $\endgroup$ Jun 16, 2018 at 22:22
  • $\begingroup$ Thanks a lot Henrik. I am going to read every word in your message for following what you are saying. In some way I solve my problem very inefficiently: using my data I construct an interpolating function. I am assuming the "NDSolveValue" will call my function only in the points of the mesh.Thanks, again. $\endgroup$ Jun 17, 2018 at 19:30
  • $\begingroup$ Hey, Eduardo, you're welcome. Using interpolating functions does work of course. But the (current implementation of) interpolating functions on unstructured grids are somewhat slow. The main reason might be that the interpolating function has to look up in which cell its argument lives. But if the argument is a vertex position and if you know its i index in the vertex list, then it is much easier to just read the i-th component of the vector describing the function. $\endgroup$ Jun 17, 2018 at 19:45

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.