Reveal the formal PDE of FiniteElement

Coincidence had it that I needed a code to reconstruct the inactive PDE that has been parsed for a customer a few weeks back. I have then added this function to the kernel and it will be available in 12.2.

The details of the operators and their specification can be found in the documentation and @andre already added links to that documentation.

Here is the code to get the inactive PDE from the NDSolve state data:

Needs["NDSolve`FEM`"]
zeroCoefficientQ[c_] := Union[N[Flatten[c]]] === {0.}
ClearAll[GetInactivePDE]
GetInactivePDE[pdec_PDECoefficientData, vd_] := 
 Module[{lif, sif, dif, mif, hasTimeQ, tvar, vars, depVars, neqn, 
   nspace, dep, load, dload, diff, cconv, conv, react, 
   pde},
  {lif, sif, dif, mif} = pdec["All"];      

  tvar = NDSolve`SolutionDataComponent[vd, "Time"];
  If[tvar === None || tvar === {}, hasTimeQ = False;
   tvar = Sequence[];, hasTimeQ = True;];
  
  vars = NDSolve`SolutionDataComponent[vd, "Space"];
  depVars = NDSolve`SolutionDataComponent[vd, "DependentVariables"];
  neqn = Length[depVars];
  nspace = Length[vars];
  dep = (# @@ Join[{tvar}, vars]) & /@ depVars;
  
  {load, dload} = lif;
  {diff, cconv, conv, react} = sif;
  
  load = load[[All, 1]];
  dload = dload[[All, 1, All, 1]];
  conv = conv[[All, All, 1, All]];
  cconv = cconv[[All, All, All, 1]];      

  pde = If[hasTimeQ, 
    mif[[1]].D[dep, {tvar, 2}] + dif[[1]].D[dep, tvar], 
    ConstantArray[0, {Length[dep]}]];
  
  If[! zeroCoefficientQ[diff], 
   pde += (Plus @@@ 
       Table[Inactive[
          Div][-diff[[r, c]].Inactive[Grad][dep[[c]], vars], 
         vars], {r, neqn}, {c, neqn}]);];
  
  If[! zeroCoefficientQ[cconv], 
   pde += (Plus @@@ 
       Table[Inactive[Div][-cconv[[r, c]]*dep[[c]], vars], {r, 
         neqn}, {c, neqn}]);];
  
  If[! zeroCoefficientQ[dload], 
   pde += (Inactive[Div][#, vars] & /@ dload);];
  
  If[! zeroCoefficientQ[conv], 
   pde += (Plus @@@ 
       Table[conv[[r, c]].Inactive[Grad][dep[[c]], vars], {r, 
         neqn}, {c, neqn}]);];
  
  pde += react.dep;
  
  pde -= load;
  
  pde
  ]

Here is an example of it's usage:

op = -x D[u[x, y], {x, 2}] - D[u[x, y], {y, 2}] - 1;
{state} = 
  NDSolve`ProcessEquations[{op == 0, 
    DirichletCondition[u[x, y] == 0, True]}, 
   u, {x, y} ∈ Disk[]
   ];
Needs["NDSolve`FEM`"]
femd = state["FiniteElementData"];
vd = state["VariableData"];
pdec = femd["PDECoefficientData"];
pde = GetInactivePDE[pdec, vd];
pde // InputForm

{-1 + {1, 0} . Inactive[Grad][u[x, y], {x, y}] + 
  Inactive[Div][-{{x, 0}, {0, 1}} . Inactive[Grad][u[x, y], {x, y}], {x, y}]}

Note, how the x in front of the D got pulled into the Div - Grad and how that is compensated by a convection component. See for example FEMDocumentation/tutorial/FiniteElementBestPractice#588198981 that explains this behavior.


I don't know if you are aware that this is documented in details.

The problem is that the informations are dispatched over the documentation of PDECoefficentData and InitializePDECoefficients.

your code :

{state} = 
  NDSolve`ProcessEquations[
   With[{u = u[x, y]}, {-2 D[u, y, y] - 3 D[u, x, x] == 1, 
     DirichletCondition[u == 0, True]}], u, {x, 0, 1}, {y, 0, 1}];

data = state["FiniteElementData"]["PDECoefficientData"];

data["All"]
(*{{{{1}},{{{{0},{0}}}}},{{{{{3,0},{0,2}}}},{{{{0},{0}}}},{{{{0,0}}}},\
{{0}}},{{{0}}},{{{0}}}}*)  

The PDECoefficentData documentation explains this :

data["ConvectionCoefficients"]
data["DampingCoefficients"]
data["MassCoefficients"]
data["LoadCoefficients"]
(* etc ... *) 

{{{{0, 0}}}}

{{0}}

{{0}}

{{1}}

InitializePDECoefficients documentation :

enter image description here

The DampingCoefficients and MassCoefficients are explained beyond.