# FEM: Solving the heat conduction with 2D periodic condition

Suppose you have a region $$\Omega = [0,P_1] \times [0,P_2]$$ which is composed of, e.g., two materials. One material is distributed as inclusions in an embedding material. We split the region as $$\Omega = \Omega_{inc} \cup \Omega_{emb}$$ (inclusion and embedding material). Example for $$P = (10,5)$$ with 3 inclusions:

(*Periodic region with periods P*)
P = {10, 5};
Omega = Rectangle[{0, 0}, P];
centers = {{1.2, 2}, {6, 3}, {8.5, 1.5}};
Omegainc = RegionUnion[Disk[#, 1] & /@ centers];
Omegaemb = RegionDifference[Omega, Omegainc];
RegionPlot[{Omegainc, Omegaemb}, AspectRatio -> Automatic,
PlotLegends -> {"\[CapitalOmega]inc", "\[CapitalOmega]emb"}]


How would you solve the following periodic 2-dimensional heat conduction problem for the unkonwn temperature $$u(x) = u(x_1,x_2) \in \mathbb{R}$$ and periodic $$v(x) \in \mathbb{R}$$ $$\mathrm{div}(A(x) \mathrm{grad}(u(x))) = 0 \quad x \in \Omega , \quad u(x) = g^T x + v(x) \quad x \in \partial \Omega \ ,$$ $$A(x) = \begin{cases} A_{inc} & x \in \Omega_{inc} \\ A_{emb} & x \in \Omega_{emb} \end{cases}$$ The constant vector $$g \in \mathbb{R}^2$$ in the boundary conditions ($$g^Tx = g_1 x_1 + g_2 x_2$$) as well as the conductivities $$A_{inc},A_{emb} \in \mathbb{R}^{2 \times 2}$$ are given, the unknown periodic field $$v(x)$$ is to be determined (periodicity for $$v(x)$$: $$v(x_1,0) = v(x_1,P_2) \ \forall x_1 \in [0,P_1]$$ and $$v(0,x_2) = v(P_1,x2) \ \forall x_2 \in [0,P_2]$$). For complete clarity, the homogeneous PDE above can also be expressed as $$\sum_{p=1}^2 \frac{\partial}{\partial x_p} \left( \sum_{q=1}^2 A_{pq}(x) \frac{\partial u(x)}{\partial x_q} \right) = 0 \quad x \in \Omega$$ How would you solve the problem above with the FEM in Mathematica?

My approach until now: I split the solution directly into $$u(x) = g^T x + v(x)$$, insert it into the PDE and solve for the periodic $$v(x)$$, i.e., solve the inhomogeneous PDE $$\mathrm{div}(A(x)g) + \mathrm{div}(A(x)\mathrm{grad}(v(x))]) = 0$$ with corresponding periodic boundary conditions for $$v(x)$$. Hereby, I am not very sure if simply inserting this into the FEM with Inactivate is actually fine for the inhomogeneity $$\mathrm{div}(A(x)g)$$. How does Mathematica treat this inhomogeneity, does it take advantage of the inactive divergence? Based on the code given above for the region generation, my solution code is given below:

1. Mesh generation
2. Region dependent coefficient $$A(x)$$
3. Boundary conditions (1 Dirichlet and periodic)
4. PDE with prescribed $$g$$
5. Solve on mesh
6. Check periodicity along edges and visualize solution (seem fine to me)

Minor question: In step 5. Solve on mesh I get the error that $$A(x)$$ can not be transposed. What am I doing wrong there?

Thanks!

(*Mesh*)
Needs["NDSolveFEM"]
mesh = ToElementMesh[
Omegaemb
, "RegionHoles" -> None
, "RegionMarker" ->
Join[{#, 1, 0.05} & /@ centers, {{{0.1, 0.1}, 2, 0.5}}]
];
mesh["Wireframe"["MeshElementStyle" -> FaceForm /@ {Blue, Orange}]]

(*Region dependent coefficient A(x)*)
Ainc = DiagonalMatrix@{100, 50};
Aemb = DiagonalMatrix@{1, 2};
A[x1_, x2_] := If[Element[{x1, x2}, Omegainc], Ainc, Aemb];

(*Region dependent coefficient A(x)*)
Ainc = DiagonalMatrix@{100, 50};
Aemb = DiagonalMatrix@{1, 2};
A[x1_, x2_] := If[Element[{x1, x2}, Omegainc], Ainc, Aemb];

(*Boundary conditions*)
bcD = DirichletCondition[v[x1, x2] == 0, x1 == 0 && x2 == 0];
gt1 = FindGeometricTransform[{{0, 0}, {0, P[[2]]}}, {{P[[1]], 0},
P}][[2]];
gt2 = FindGeometricTransform[{{0, 0}, {P[[1]], 0}}, {{0, P[[2]]},
P}][[2]];
bcP = {
PeriodicBoundaryCondition[
v[x1, x2]
, x1 == P[[1]] && 0 <= x2 <= P[[2]]
, gt1
]
,
PeriodicBoundaryCondition[
v[x1, x2]
, x2 == P[[2]] && 0 <= x1 <= P[[1]]
, gt2
]
};

(*PDE with prescribed g*)
g = {3, 1};
pde = Inactive[Div][A[x1, x2].g, {x1, x2}] +
Inactive[Div][
A[x1, x2].Inactive[Grad][v[x1, x2], {x1, x2}], {x1, x2}] == 0;

(*Solve on mesh*)
vsol = NDSolveValue[{pde, bcD, bcP}, v, Element[{x1, x2}, mesh]];

(*Check periodictiy along edges and visualize solution*)
Plot[
vsol[x1, 0] - vsol[x1, P[[2]]], {x1, 0, P[[1]]}, PlotRange -> All,
PlotLegends -> {"v[x1,0]-v[x1,P2]"}]
Plot[vsol[0, x2] - vsol[P[[1]], x2], {x2, 0, P[[2]]},
PlotRange -> All, PlotLegends -> {"v[0,x2]-v[P1,x2]"}]
Show[ContourPlot[vsol[x1, x2], Element[{x1, x2}, Omega],
AspectRatio -> Automatic, PlotLegends -> Automatic],
RegionPlot@Omegainc]
Plot3D[vsol[x1, x2], Element[{x1, x2}, Omega]]


• I only want $v(x)$ to be periodic throughout the boundary. That is sufficient for the uniqueness of the solution (up to a constant, of course), as far as I remember, I do not need to specify anything for the flux. – Mauricio Fernández Dec 29 '19 at 23:56
• You use an undifferentiable function as differentiable in Inactive[Div][A[x1, x2].g, {x1, x2}]. Therefore, we have a message. But FEM it is not BEM. You have to use NeumannValue[]. – Alex Trounev Dec 30 '19 at 16:51
• @AlexTrounev that is not a problem. For such PDEs, you can also use, e.g., linear elements such that $A(x) \mathrm{grad}(v(x))$ would not be differentiable when passing it to $\mathrm{div}$. Using the weak formulation and integration by parts allows to get rid of the divergence operation. This should also not be a problem for the term $A(x)g$. I just found out that this case fits the term $\gamma$ in the general PDE in the FEM documentation. But I still dont get it, why I get an error there. – Mauricio Fernández Dec 30 '19 at 23:22
• Then see my answer – Alex Trounev Dec 31 '19 at 12:46

This is similar to this question and this answer also applies here.

When you use

(*PDE with prescribed g*)
g = {{3}, {1}};


Things then work as expected. The Inactive Div then does the differentiation. It's worthwhile noting that the numerical differentiation is probably less accurate then if you have a symbolic derivative.

Note the other improvements in that post also apply here.

For completeness the plot:

• Thanks, I also realized after some digging that the $\gamma$ term was the issue and decided to ask a question specifically for it in a non-periodic setting. – Mauricio Fernández Jan 1 at 11:08
• @user21 it looks like a solution, but what describes this solution and how to find out what is right? – Alex Trounev Jan 1 at 12:11
• @AlexTrounev, $\gamma$ is used for computing derivatives (divergences) of functions. This sometime is convenient to do in the FEM code, as in this case. How to check if this is correct, well like in any other case you construct an example for which you have an analytical solution: – user21 Jan 1 at 12:31
• Needs["NDSolveFEM"] test = Sin[x]; nr = ToNumericalRegion[Line[{{0}, {2 \[Pi]}}]]; vd = NDSolveVariableData[{"DependentVariables", "Space"} -> {{u}, {x}}]; sd = NDSolveSolutionData["Space" -> nr]; pdec = InitializePDECoefficients[vd, sd, "ReactionCoefficients" -> {{1}}, "LoadDerivativeCoefficients" -> {{{-test}}}]; mdata = InitializePDEMethodData[vd, sd]; dpde = DiscretizePDE[pdec, mdata, sd]; {l, s, d, m} = dpde["All"]; res = LinearSolve[s, l]; if = ElementMeshInterpolation[{nr["ElementMesh"]}, res]; Plot[Evaluate[if[x] - D[test, x]], {x, 0, 2 \[Pi]}] – user21 Jan 1 at 12:32
• @user21 I asked not at all, but in this particular case. At first glance, my solution describes the same thing, but the result is completely different. – Alex Trounev Jan 1 at 12:48

We have to use either NeumannValue[] or save the linear function in the pde as follows

(*Mesh*)Needs["NDSolveFEM"]
(*Periodic region with periods P*)

P = {10, 5};
Omega = Rectangle[{0, 0}, P];
centers = {{1.2, 2}, {6, 3}, {8.5, 1.5}};
Omegainc = RegionUnion[Disk[#, 1] & /@ centers];
Omegaemb = RegionDifference[Omega, Omegainc];
mesh = ToElementMesh[Omegaemb, "RegionHoles" -> None,
"RegionMarker" ->
Join[{#, 1, 0.05} & /@ centers, {{{0.1, 0.1}, 2, 0.5}}]];
{RegionPlot[{Omegainc, Omegaemb}, AspectRatio -> Automatic,
PlotLegends -> {"\[CapitalOmega]inc", "\[CapitalOmega]emb"}],
mesh["Wireframe"["MeshElementStyle" -> FaceForm /@ {Blue, Orange}]]}

(*Region dependent coefficient A(x)*)
Ainc = DiagonalMatrix@{100, 50};
Aemb = DiagonalMatrix@{1, 2};
A[x1_, x2_] := Boole[Element[{x1, x2}, Omegainc]] (Ainc - Aemb) + Aemb;

(*Boundary conditions*)eps = 10^-3;
bcD = DirichletCondition[v[x1, x2] == RandomReal[{-eps, eps}],
x1 == 0 && x2 == 0];
gt1 = FindGeometricTransform[{{0, 0}, {0, P[[2]]}}, {{P[[1]], 0},
P}][[2]];
gt2 = FindGeometricTransform[{{0, 0}, {P[[1]], 0}}, {{0, P[[2]]},
P}][[2]];
bcP = {PeriodicBoundaryCondition[v[x1, x2],
x1 == P[[1]] && 0 <= x2 <= P[[2]], gt1],
PeriodicBoundaryCondition[v[x1, x2],
x2 == P[[2]] && 0 <= x1 <= P[[1]], gt2]};

(*PDE with prescribed g*)
g = 3 x1 + x2;
pde = Inactive[Div][
A[x1, x2].Inactive[Grad][v[x1, x2], {x1, x2}], {x1, x2}] +
Inactive[Div][A[x1, x2].Inactive[Grad][g, {x1, x2}], {x1, x2}] ==
0;

(*Solve on mesh*)
vsol = NDSolveValue[{pde, bcP, bcD}, v, Element[{x1, x2}, mesh]]
(*Check periodictiy along edges and visualize solution*)
{Plot[vsol[x1, 0] - vsol[x1, P[[2]]], {x1, 0, P[[1]]},
PlotRange -> {-eps, eps}, PlotLegends -> {"v[x1,0]-v[x1,P2]"}],
Plot[vsol[0, x2] - vsol[P[[1]], x2], {x2, 0, P[[2]]},
PlotRange -> {-eps, eps}, PlotLegends -> {"v[0,x2]-v[P1,x2]"}],
Show[ContourPlot[vsol[x1, x2], Element[{x1, x2}, Omega],
AspectRatio -> Automatic, Contours -> 20,
ColorFunction -> "TemperatureMap", PlotLegends -> Automatic],
RegionPlot@Omegainc]}


• That does not make sense to me. The solution you get is a constant throughout the region, equal to the random value between -eps and +eps you initialized the Dirichlet condition at $(x_1,x_2) = (0,0)$. See your plot. You can put 0.3 as the Dirichlet condition and the returned solution will be 0.3 throughout the region. That seems not right. – Mauricio Fernández Dec 31 '19 at 13:43
• What solution do you want using PeriodicBoundaryCondition[]? – Alex Trounev Dec 31 '19 at 13:47
• I expect the solution $v(x)$ to do weird stuff inside the region due to the 2 different materials and to be periodic (see, e.g., the solution plot at the bottom of my question). That is all. A constant $v(x)$ trivially vanishes in the gradient such that $A(x) \mathrm{grad}(v(x))$ would vanish overall. That does not seem right to me. – Mauricio Fernández Dec 31 '19 at 13:52
• You use a linear function that is not periodic. What do you want to get? – Alex Trounev Dec 31 '19 at 14:02
• No, you are wrong, because the system found an error in your pde. On the other hand, there are no errors in my pde. And you see what happened. Try to formulate the problem as it is. What do you want to describe using this code? – Alex Trounev Dec 31 '19 at 14:33