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:
- Mesh generation
- Region dependent coefficient $A(x)$
- Boundary conditions (1 Dirichlet and periodic)
- PDE with prescribed $g$
- Solve on mesh
- 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["NDSolve`FEM`"]
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]]
Inactive[Div][A[x1, x2].g, {x1, x2}]. Therefore, we have a message. But FEM it is not BEM. You have to useNeumannValue[]. $\endgroup$