Finite element method solution data, Mathematica 13.1

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["NDSolveFEM"]

{state} =
NDSolveProcessEquations[{Laplacian[u[x, y], {x, y}] == 1,
DirichletCondition[u[x, y] == 0, True]}, u, {x, 0, 1}, {y, 0, 1}]

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

NDSolveIterate[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?

• The "Solution" method is already missing at least since v12.3. Not sure if it's a bug or compatibility issue. If I guess it right, the solution data can still be extracted with state["SolutionData"][[1, 3]] or state["SolutionVector"][[1]]. Oct 26, 2022 at 7:53
• I see! I did not see those listed when I executed state["Properties"]. I guess these are hidden. Thanks! Oct 26, 2022 at 8:08

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["NDSolveFEM"]
{state} =
NDSolveProcessEquations[{Laplacian[u[x, y], {x, y}] == 1,
DirichletCondition[u[x, y] == 0, True]},
u, {x, 0, 1}, {y, 0, 1}];
NDSolveIterate[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

NDSolveSolutionDataComponent[forwardSolution, "DependentVariables"]


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

Needs["NDSolveFEM"]
nr = ToNumericalRegion[Rectangle[{0, 0}, {1, 1/2}]];
vd = NDSolveVariableData[{"DependentVariables",
"Space"} -> {{u}, {x, y}}];
sd = NDSolveSolutionData["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:

NDSolveSolutionDataComponent[sdNew, "DependentVariables"] // Short
`

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