# Reveal the formal PDE of FiniteElement

When FiniteElement method is used, the differential equations will first be transformed to certain standard form (named as formal PDE in recent FEM document), and it turns out to be critical to check what the standard form is when analyzing various issues related to FEM. Here are some examples:

Position of discontinuous coefficient influences the solution of PDE

How to input Robin boundary conditions for nonstandard Laplace equation?

Sign of conservative convection coefficient in a formal (Inactive) PDE

Stress analysis in axisymmetric bodies

The coefficient of formal PDE is available from PDECoefficientData, but its output is just hard to read. For example, with

{state} = NDSolveProcessEquations[
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}}}} *)


at hand, can you tell what's what? Can you label $$d$$, $$c$$, $$\alpha$$, etc. in the formal PDE

$$d\frac{\partial }{\partial t}u+\nabla \cdot (-c \nabla u-\alpha u+\gamma ) +\beta \cdot \nabla u+ a u -f=0$$

with corresponding values, without doubt?

Can we have a function that shows the formal PDE of FiniteElement in a easy-to-read way? A possible (but not necessary of course) input-output in my mind:

showFormalPDE@With[{u = u[x, y]}, -2 D[u, y, y] - 3 D[u, x, x] == 1]
(* -1 + Inactive[Div][(-{{3, 0}, {0, 2}}.Inactive[Grad][u[x, y], {x, y}]), {x, y}] == 0 *)


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.

In versions after 12.2 you can use:

NDSolveFEMGetInactivePDE


Else, here is the code to get the inactive PDE from the NDSolve state data:

Needs["NDSolveFEM"]
zeroCoefficientQ[c_] := Union[N[Flatten[c]]] === {0.}
ClearAll[GetInactivePDE]
GetInactivePDE[pdec_PDECoefficientData, vd_] :=
Module[{lif, sif, dif, mif, hasTimeQ, tvar, vars, depVars, neqn,
pde},
{lif, sif, dif, mif} = pdec["All"];

tvar = NDSolveSolutionDataComponent[vd, "Time"];
If[tvar === None || tvar === {}, hasTimeQ = False;
tvar = Sequence[];, hasTimeQ = True;];

vars = NDSolveSolutionDataComponent[vd, "Space"];
depVars = NDSolveSolutionDataComponent[vd, "DependentVariables"];
neqn = Length[depVars];
nspace = Length[vars];
dep = (# @@ Join[{tvar}, vars]) & /@ depVars;

{diff, cconv, conv, react} = sif;

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[
vars], {r, neqn}, {c, neqn}]);];

If[! zeroCoefficientQ[cconv],
pde += (Plus @@@
Table[Inactive[Div][-cconv[[r, c]]*dep[[c]], vars], {r,
neqn}, {c, neqn}]);];

pde += (Inactive[Div][#, vars] & /@ dload);];

If[! zeroCoefficientQ[conv],
pde += (Plus @@@
neqn}, {c, neqn}]);];

pde += react.dep;

pde
]

GetInactivePDE[state_NDSolveStateData] := Module[{femd = state["FiniteElementData"],
vd = state["VariableData"], pdec},
pdec = femd["PDECoefficientData"];
GetInactivePDE[pdec, vd]]


Here is an example of it's usage:

op = -x D[u[x, y], {x, 2}] - D[u[x, y], {y, 2}] - 1;
{state} =
NDSolveProcessEquations[{op == 0,
DirichletCondition[u[x, y] == 0, True]},
u, {x, y} ∈ Disk[]
];

pde = GetInactivePDE[state];
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.

• Side note: NDSolveFEMGetInactivePDE  is built-in in v12.2. (There's no separate document page for it though. ) Apr 29 at 6:21
• I take the liberty to add a line making GetInactivePDE in this post handle NDSolveStateData, feel free to roll back if you don't like it :) . Apr 30 at 13:25

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.

{state} =
NDSolveProcessEquations[
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"]
(* etc ... *)


{{{{0, 0}}}}

{{0}}

{{0}}

{{1}}

InitializePDECoefficients documentation :

The DampingCoefficients and MassCoefficients are explained beyond.

• Oh, I didn't notice the Details section of InitializePDECoefficients… but I think it'll be good to have a function that can generate the formal PDE automatically, checking the table in the document is still a bit painful. Jul 11 '20 at 8:09
• It took me a while to notice this too (I don't rememder exactly, but it was after having read the tutorials). Jul 11 '20 at 8:16