Finite element method solution data, Mathematica 13.1

Question

I just upgraded to Mathematica 13.1 and encountered the following issue (this is just code from the first few lines of the FEM programming documentation):

Needs["NDSolve`FEM`"]

{state} = 
 NDSolve`ProcessEquations[{Laplacian[u[x, y], {x, y}] == 1, 
   DirichletCondition[u[x, y] == 0, True]}, u, {x, 0, 1}, {y, 0, 1}]
(* {NDSolve`StateData["<" "SteadyState" ">"]} *)

femdata = state["FiniteElementData"]
(* FiniteElementData["<" 1281 ">"] *)

NDSolve`Iterate[state]; 

state["FiniteElementData"]["Solution"] // Short

which resulted in the following error

During evaluation of In[16]:= FiniteElementData::nomthd: There is no method Solution for FiniteElementData objects.

Where did the solution data go? Is this a known issue?

Answer 1

In essence this is a documentation issue. To extract the solution, there are two scenarios. Either you have a NDSolve state object (like you have shown) or you use a lower level function like PDESolve.

Consider this case:

Needs["NDSolve`FEM`"]
{state} = 
  NDSolve`ProcessEquations[{Laplacian[u[x, y], {x, y}] == 1, 
    DirichletCondition[u[x, y] == 0, True]}, 
   u, {x, 0, 1}, {y, 0, 1}];
NDSolve`Iterate[state];

On this level you should use the solution data from the state object.

{forwardSolution} = state["SolutionData"];

For more information on the solution data see the documentation of the solution data object.

The actual data extraction should then happen with SolutionDataComponent

NDSolve`SolutionDataComponent[forwardSolution, "DependentVariables"]

The second case is when you use low level functions, like in this case:

Needs["NDSolve`FEM`"]
nr = ToNumericalRegion[Rectangle[{0, 0}, {1, 1/2}]];
vd = NDSolve`VariableData[{"DependentVariables", 
     "Space"} -> {{u}, {x, y}}];
sd = NDSolve`SolutionData["Space" -> nr];
methodData = InitializePDEMethodData[vd, sd]
pdeData = InitializePDECoefficients[vd, sd, 
   "DiffusionCoefficients" -> {{-IdentityMatrix[2]}}];
bcData = InitializeBoundaryConditions[vd, 
   sd, {{DirichletCondition[u[x, y] == 0., x == 0], 
     DirichletCondition[u[x, y] == 1., x == 1]}}];

Now we call PDESolve:

sdNew = PDESolve[pdeData, bcData, vd, sd, methodData];

Note that PDESolve returns a new solution data object.

sdNew === sd
(* False *)

You can then extract the solution vector from the object in the same manner as above:

NDSolve`SolutionDataComponent[sdNew, "DependentVariables"] // Short

This is a much cleaner process than what we had before. I have updated the reference pages to reflect this.