Convert boundary condition involving derivative to NeumannValue programmatically

Question

As discussed in e.g.

what causes the error "The dependent variable in the boundary condition needs to be linear" when using NDSolve?

When FiniteElement method is chosen for spatial discretization in NDSolve, boundary condition (b.c.) involving Derivative (include but not limited to Neumann b.c. and Robin b.c. ) cannot be directly used. We have to convert such b.c. to NeumannValue. Can the process be automated, at least for regular domain?

Answer 1

The following my (incomplete) attempt.

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

SetAttributes[allowfemdbc, HoldAll];
allowfemdbc[(head : NDSolve | NDSolveValue | NDEigenvalues | NDEigensystem)[sys_, 
   depend_, domain : Longest[{_, _, _} ..], n___Integer, rest : OptionsPattern[]], 
  wrapper_ : Identity] := 
 Module[{coordlst, boundLlst, boundRlst, dependvar, bc, bcfem, eq, dbc, bcrest, state,
    coefgrad, dim, asso, pde, neumann, dvar, bcpos, bcnew, gpu, 
   gpuvalue}, {coordlst, boundLlst, boundRlst} = {domain}\[Transpose];
  {dependvar} = Flatten@{depend} /. h_ @@ coordlst :> h;
  {bc, bcfem, eq} = 
   Internal`ProcessEquations`SeparateEquations[sys // forcetoeq, coordlst, 
     dependvar][[2 ;; 4]];
  dbc = Cases[bc, a_ /; ! FreeQ[a, Derivative]];
  bcrest = Complement[bc, dbc];
  state = Quiet[
    NDSolve`ProcessEquations[{eq, bcrest, bcfem}, depend, domain, rest][[1]], 
    NDSolve`ProcessEquations::femibcnd];
  coefgrad = 
   state["FiniteElementData"]["PDECoefficientData"]["DiffusionCoefficients"][[1]] . 
    Grad[-dependvar @@ coordlst, coordlst];
  dim = Length@coordlst;
  asso = Association@
    Flatten@{Thread[Thread[coordlst == boundLlst] -> -IdentityMatrix[dim]], 
      Thread[Thread[coordlst == boundRlst] -> IdentityMatrix[dim]]};
  pde = NDSolve`FEM`GetInactivePDE[state][[1]];
  neumann = 
   Total@Table[{dvar, bcpos} = 
      Cases[traditionalbc, 
        term : Derivative[__][dependvar][a__] :> {term, 
          And @@ Thread[coordlst == {a}]}, Infinity, 1][[1]];
     bcnew = gpu == coefgrad . asso[bcpos] /. Rule @@ bcpos;
     gpuvalue = SolveValues[{bcnew, traditionalbc}, gpu, {dvar}][[1]];
     NeumannValue[gpuvalue /. h_ @@ dvar :> dependvar @@ coordlst, 
      bcpos], {traditionalbc, dbc}];
  wrapper[head][
   Sow@Flatten@{If[MatchQ[head, NDEigenvalues | NDEigensystem], -(pde + neumann), 
       pde == neumann], bcrest, bcfem}, depend, domain, n, rest]]

To understand the code, you may want to read the following posts:

What is the best way to parse differential equations and boundary conditions to a custom function?

Reveal the formal PDE of FiniteElement

Known issues

1. Currently it cannot deal with PDE system.

2. Time dependent PDE is not supported for now, but when the PDE is time dependent in a regular domain, we can use the old good TensorProductGrid, don't we?

3. ∈ syntax for specifying domain of definition is not supported for the moment.

4. allowfemdbc is no more than a converter for NeumannValue, so if certain b.c. cannot be expressed with NeumannValue, allowfemdbc won't help, either.

Usage

The usage allowfemdbc is simple, just add it before NDSolve / NDSolveValue / NDEigenvalues NDEigensystem!

If you feel worried about the correctness of the obtained NeumannValue, set 2nd argument of allowfemdbc to Inactive to check the generated code. You can also add a Reap to extract the generated system.

Example

1

sys = {y''[x] + y[x] == 0, y[0] == 1, y'[30] == 1};

asol = DSolveValue[sys, y, {x, 0, 30}];

nsol = allowfemdbc@NDSolveValue[sys, y, {x, 0, 30}, Method -> FiniteElement];

Plot[{asol, nsol}[x] // Through // Evaluate, {x, 0, 30}, 
 PlotStyle -> {Automatic, Dashed}]

2

How to input Robin boundary conditions for nonstandard Laplace equation?

Clear[bc]
xmax = 1; ymax = 1;
epsilon = 0.5; n = 0.1;
With[{u = u[x, y]}, eq = epsilon^2 D[u, x, x] + D[u, y, y] == 1;
 {bc@x, bc@y} = {{D[u, x] == 0 /. x -> 0, 
    u == -2 epsilon n D[u, x] /. x -> xmax}, {D[u, y] == 0 /. y -> 0, 
    u == -2 n D[u, y] /. y -> ymax}};]

solfem = allowfemdbc@NDSolveValue[{eq, bc /@ {x, y}}, 
                                  u, {x, 0, xmax}, {y, 0, ymax}];

Plot3D[solfem[x, y], {x, 0, xmax}, {y, 0, ymax}]

3

How to solve 2D eigenvalue problem with robin boundary conditions

With[{u = u[x, y]}, lhs = Laplacian[u, {x, y}];
  bc = {u == 0 /. {{x -> -1}, {y -> -1}},
    {2 D[u, x] + u == 0 /. x -> 1,
     D[u, y] + u == 0 /. y -> 1}}];

tsthold = allowfemdbc[
  NDEigensystem[{lhs, bc} // Flatten, u, {x, -1, 1}, {y, -1, 1}, 4, 
   Method -> {"PDEDiscretization" -> {"FiniteElement", 
       "MeshOptions" -> {"MaxCellMeasure" -> 0.01}}}], Inactive]
(*
Inactive[NDEigensystem][{-NeumannValue[-u[x, y], y == 1] - 
   NeumannValue[-(1/2) u[x, y], x == 1] - 
   Inactive[Div][(-{{1, 0}, {0, 1}} . Inactive[Grad][u[x, y], {x, y}]), {x, y}], 
  u[-1, y] == 0, u[x, -1] == 0}, u, {x, -1, 1}, {y, -1, 1}, 4, 
 Method -> {"PDEDiscretization" -> {"FiniteElement", 
     "MeshOptions" -> {"MaxCellMeasure" -> 0.01}}}]
 *)

tst = Activate[tsthold];
Plot3D[tst[[2, 1]][x, y], {x, -1, 1}, {y, -1, 1}]