FEM: does the term $\gamma$ in "The Coefficient Form of Partial Differential Equations" has to be an explicit list?

Question

In the FEM documentation, "The Coefficient Form of Partial Differential Equations" is

I am interested in its simplification to

$$ \nabla \cdot (-c \nabla u + \gamma) = 0 $$

where $c$ should be a matrix and $\gamma$ should be a vector, see FEM documentation.

Question: does $\gamma$ have to be an explicit vector/list or can it be defined through functions returning a vector?

Example for explicit list definition in a 2D problem:

gamma[x1_,x2_]:={Exp[x1],Exp[x2]}

Example for function returning vectors in a 2D problem:

gamma[x1_,x2_]:=If[Element[{x1,x2},Disk[]],{1,2},{80,50}]

The reason for this question is that I am trying to solve a PDE with the FEM in Mathematica and I get an error which I do not understand, see example below if you are interested. Personally I am confused, since the matrix coefficient $c(x)$ in the general form can be defined through functions returning the corresponding matrices, see "Partial Differential Equations with Variable Coefficients" in the FEM documentation. Since functions returning matrices are allowed for $c(x)$, I was expecting the same for the vector $\gamma(x)$.

---

Problem 1 (standard, $\gamma = 0$)

In a rectangle $\Omega = [0,L_1] \times [0,L_2]$ with given piecewise constant $A(x)$ solve $$ \nabla \cdot (A(x) \nabla u(x)) = 0 \quad x \in \Omega $$ with boundary conditions (vector $g$ is given) $$ u(x) = g^Tx = g_1 x_1 + g_2 x_2 \quad x \in \partial \Omega \ . $$

Problem 2 ($\gamma \neq 0$)

In the very same region $\Omega$ of Problem 1 with the very same $A(x)$ and $g$ consider the linear superposition $u(x) = g^T x + v(x)$, such that $$ A(x) \nabla u(x) = A(x)g + A(x) \nabla v(x) $$ holds. Defining $$ \gamma(x) = A(x)g $$ and inserting the split into the PDE yields the equivalent problem $$ \nabla \cdot (A(x) \nabla v(x) + \gamma(x)) = 0 $$ with boundary conditions $$ v(x) = 0 \quad x \in \partial \Omega \ . $$

Code

Below you will find the complete code for the solution of Problem 1 and Problem 2, at what for Problem 2 I define 3 mathematically equivalent versions of $\gamma(x)$, but which have differences in the evaluation in Mathematica's FEM. Surprisingly,

$$ \nabla \cdot (A(x)g + A(x)\nabla v(x)) = 0 $$

is not acceptable for the FEM routine. I have to put in the PDE as follows

$$ \nabla \cdot (A(x)g) + \nabla \cdot (A(x)\nabla v(x)) = 0 $$

which computes the correct field $v(x)$ but raises the error

The error does not appear if you use the function gamma3 (see code), which is an explicit list definition. Am I doing something wrong? The function gamma2 does not work and I do not get why. What am I doing wrong?

Region, mesh and coefficient A(x)

(*Region*)
L = {5, 4};
Omega = Rectangle[{0, 0}, L];
Omegainc = Disk[{3, 2}, 1];
Omegaemb = RegionDifference[Omega, Omegainc];
RegionPlot[{Omegainc, Omegaemb}, AspectRatio -> Automatic, 
 PlotLegends -> {"\[CapitalOmega]inc", "\[CapitalOmega]emb"}]

(*Mesh*)
Needs["NDSolve`FEM`"]
mesh = ToElementMesh[Omegaemb, "RegionHoles" -> None, 
   "RegionMarker" -> {
     {{3, 2}, 1, 0.01}
     , {{0.1, 0.1}, 2, 0.5}
     }];
mesh["Wireframe"["MeshElementStyle" -> FaceForm /@ {Blue, Orange}]]

(*Region dependent coefficient A(x)*)

Ainc = DiagonalMatrix@{100, 50};
Aemb = DiagonalMatrix@{10, 20};
A[x1_, x2_] := If[Element[{x1, x2}, Omegainc], Ainc, Aemb];

Solution of Problem 1

(*Boundary conditions for u*)
g = {1, 0};
bcD = {
   DirichletCondition[u[x1, x2] == g.{x1, x2}, x1 == 0]
   , DirichletCondition[u[x1, x2] == g.{x1, x2}, x1 == L[[1]]]
   , DirichletCondition[u[x1, x2] == g.{x1, x2}, x2 == 0]
   , DirichletCondition[u[x1, x2] == g.{x1, x2}, x2 == L[[2]]]
   };
(*PDE, solve for u and visualize*)

pde = Inactive[Div][
    A[x1, x2].Inactive[Grad][u[x1, x2], {x1, x2}], {x1, x2}] == 0;
usol = NDSolveValue[{pde, bcD}, u, Element[{x1, x2}, mesh]];
Show[ContourPlot[usol[x1, x2], Element[{x1, x2}, Omega], 
  AspectRatio -> Automatic, PlotLegends -> Automatic], 
 RegionPlot@Omegainc, PlotLabel -> "u(x)"]
Plot3D[usol[x1, x2], Element[{x1, x2}, Omega], PlotLabel -> "u(x)"]

Solutions for Problem 2

(*Boundary conditions for deviation v from g.x*)
bcD = {
   DirichletCondition[v[x1, x2] == 0, x1 == 0]
   , DirichletCondition[v[x1, x2] == 0, x1 == L[[1]]]
   , DirichletCondition[v[x1, x2] == 0, x2 == 0]
   , DirichletCondition[v[x1, x2] == 0, x2 == L[[2]]]
   };
(*PDE, solution for v and visualize*)

pde = Inactive[Div][
     A[x1, x2].Inactive[Grad][v[x1, x2], {x1, x2}], {x1, x2}] + 
    Inactive[Div][A[x1, x2].g, {x1, x2}] == 0;
vsol = NDSolveValue[{pde, bcD}, v, Element[{x1, x2}, mesh]];
Show[ContourPlot[vsol[x1, x2], Element[{x1, x2}, Omega], 
  AspectRatio -> Automatic, PlotLegends -> Automatic], 
 RegionPlot@Omegainc, PlotLabel -> "v(x)"]
ContourPlot[usol[x1, x2] - (g.{x1, x2} + vsol[x1, x2]), 
 Element[{x1, x2}, Omega], PlotLegends -> Automatic, 
 AspectRatio -> Automatic, PlotLabel -> "u(x) - (g.x + v(x))"]

(*Different versions of gamma[x] for FEM*)

gamma1[x1_, x2_] := A[x1, x2].g;
gammainc = Ainc.g;
gammaemb = Aemb.g;
gamma2[x1_, x2_] := 
  If[Element[{x1, x2}, Omegainc], gammainc, gammaemb];
gamma3[x1_, x2_] := 
  If[Element[{x1, x2}, Omegainc], gammainc[[#]], gammaemb[[#]]] & /@ 
   Range@2;
(*PDE with gamma, solve for v and check*)

pde = Inactive[Div][gamma3[x1, x2], {x1, x2}] + 
    Inactive[Div][
     A[x1, x2].Inactive[Grad][v[x1, x2], {x1, x2}], {x1, x2}] == 0;
vsolgamma = NDSolveValue[{pde, bcD}, v, Element[{x1, x2}, mesh]];
ContourPlot[vsolgamma[x1, x2], Element[{x1, x2}, Omega], 
 AspectRatio -> Automatic, PlotLegends -> Automatic, 
 PlotLabel -> 
  "\!\(\*SubscriptBox[\(v\), \(\[Gamma]\)]\)(x) (based on chosen \
gamma[x])"]
ContourPlot[vsol[x1, x2] - vsolgamma[x1, x2], 
 Element[{x1, x2}, Omega], AspectRatio -> Automatic, 
 PlotLegends -> Automatic, 
 PlotLabel -> "v(x) - \!\(\*SubscriptBox[\(v\), \(\[Gamma]\)]\)(x)"]

Answer 1

I think the best way to see the exact definitions of the coefficients is on the InitializePDECoefficients ref page. $\gamma$ needs to be a vector of length n. Unfortunately, it's not as easy as that. Internally, these coefficients are converted to matrices and in this case that does not work well. It's easy to work around this. However, there are a few other issues in your question that I'd like to address also. It's probably better to use a RegionMemberFunction and to inject all values into the If statement.

Ainc = DiagonalMatrix@{100, 50};
Aemb = DiagonalMatrix@{10, 20};
(*rmf=RegionMember[Omegainc];
With[{Ainc=Ainc,Aemb=Aemb,rmf=rmf},
A[x1_,x2_]:=If[rmf[{x1,x2}],Ainc,Aemb];
]*)

It's even better (more efficient) to use ElementMarkers:

With[{Ainc = Ainc, Aemb = Aemb},
 A[x1_, x2_] := If[ElementMarker == 1, Ainc, Aemb];
 ]

For the actual computation you then need to change g to:

g2 = Partition[g, 1]
(*{{1}, {0}}*) 

With this you can then use either:

pde = Inactive[Div][
   Inactive[Plus][A[x1, x2].Inactive[Grad][v[x1, x2], {x1, x2}], 
    A[x1, x2].g2], {x1, x2}] == 0

or

pde = Inactive[Div][
    A[x1, x2].Inactive[Grad][v[x1, x2], {x1, x2}], {x1, x2}] + 
   Inactive[Div][A[x1, x2].g2, {x1, x2}] == 0

Hope this helps.