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 :
The DampingCoefficients
and MassCoefficients
are explained beyond.