11
$\begingroup$

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} = 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}}}} *)

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 *)
$\endgroup$
9
$\begingroup$

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.

| improve this answer | |
$\endgroup$
7
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – xzczd Jul 11 at 8:09
  • 1
    $\begingroup$ It took me a while to notice this too (I don't rememder exactly, but it was after having read the tutorials). $\endgroup$ – andre314 Jul 11 at 8:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.