# 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, 2021 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, 2021 at 13:25
• @user21 Looking at the inactive form of the pde: The Div-part defines the NeumannValue (-{{x, 0}, {0, 1}} . Grad [U[x, y], {x, y}]) . Is this conclusion correct? Thanks! May 16, 2022 at 14:08
• @UlrichNeumann, yes, whatever is in the Div part defines NeumannValue. May 17, 2022 at 11:22
• @user21 Thank you, in this question Understanding NeumannValue I 'm trying, in vain, to verify this result. Please read my question. May 18, 2022 at 9:47

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, 2020 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, 2020 at 8:16

The problem has been resolved by user21. I'd like to add a wrapper for NDSolveFEMGetInactivePDE to make it easier to use:

ClearAll[checkFormalPDE]

toeq[sys_] :=
Replace[Flatten@{sys},
expr : Except[_Equal | _DirichletCondition | _PeriodicBoundaryCondition] :>
expr == 0, {1}];

SetAttributes[checkFormalPDE, HoldAll];

With[{pe = NDSolveProcessEquations, gip = NDSolveFEMGetInactivePDE},
checkFormalPDE[(head : NDSolve | NDSolveValue | NDEigenvalues | NDEigensystem)[sys_,
depend_, domain__, n___Integer, rest : OptionsPattern[]]] :=
Module[{state = Sow@First@pe[sys // toeq, depend, domain, rest]}, state // gip]];


## Usage

Simply wrap it on NDSolve/NDSolveValue/NDEigenvalues/NDEigensystem:

checkFormalPDE@
NDSolve[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}]
(* {-1 + Inactive[Div][(-{{3, 0}, {0, 2}} . Inactive[Grad][u[x, y], {x, y}]), {x, y}]} *)
`
• Nice. Might be foos for the function repository? Feb 16, 2023 at 5:38
• @user21 Submitted for review :) . Feb 16, 2023 at 9:01
• Great, please share the link once it is live. Feb 16, 2023 at 15:58