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)"] 
