Suppose you have a region $ Omega = (0, P_1) times (0, P_2) $ which is made, for example, of two materials. A material is distributed in the form of inclusions in an integration material. We have divided the region as $ Omega = Omega_ {inc} cup Omega_ {emb} $ (inclusion and integration 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 two-dimensional periodic heat conduction problem for the unknown 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 $ 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 clarity, the above PDE 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 above problem with the FEM in Mathematica?**

**My approach so far:** I divide the solution directly into $ u (x) = g ^ T x + v (x) $ and solve for periodic $ v (x) $, that is to say solve

$$

mathrm {div} (A (x) g) + mathrm {div} (A (x) mathrm {grad} (v (x)))) = 0

$$

with the corresponding boundary conditions. Hereby, I am not very sure if I would just insert it into the GEF with `Inactivate`

is actually good for inhomogeneity $ mathrm {div} (A (x) g) $. How does Mathematica deal with this inhomogeneity, take advantage of the inactive divergence? Is there a better, simpler and more direct way to attack the original problem? Based on the code given above for generating the region, my solution code is given below:

- Mesh generation
- Coefficient depending on the region $ A (x) $
- Boundary conditions (1 Dirichlet and periodical)
- PDE with prescribed $ g $
- Solve on the mesh
- Check the periodicity along the edges and visualize the solution (seems good to me)

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

Thank you!

```
(*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)
, RegionPlot@Omegainc
)
```