# Solving a PDE with a spatially piecewise-constant material parameter

I would like to solve equation for the 2D distribution of an electrical potential in the conductive medium:

$$\nabla\cdot(\sigma\nabla f) = 0$$

where $f=f(x,y)$ is the electric potential (in 2D) and sigma is the conductivity. The domain where the equation is solved is complex. For instance, like this:

 Manipulate[
t1 = Triangle[{{0, 0}, {x, y}, {1, 0}}];
t2 = Triangle[{{0, 1}, {x, y}, {1, 1}}];
t3 = Triangle[{{0, 0}, {x, y}, {0, 1}}];
t4 = Triangle[{{1, 0}, {x, y}, {1, 1}}];
Graphics[{Blue, t1, Red, t2, Green, t3, Yellow, t4},
ImageSize -> 300],
{x, 0, 1}, {y, 0, 1}]


Looking like the following:

The point is that each of the triangles should have its own value of the conductivity. It is constant within a triangle and changes by transition from one triangle to another. Say, for the sake of the trial, $\sigma=1, 2, 3$ and $4$. I would like to be able to manipulate the structure of the domain (as it is in the image shown above) and dynamically solve the equation depending upon the values $(x,y)$ fixed by the slider. The boundary conditions might be:

{ DirichletCondition[{f[x, y] == 0.}, y == 0],
DirichletCondition[{f[x, y] == 1.}, y == 1]}


and the no-flux conditions on two other sides of the square:

NeumannValue[0, x == 0 || x == 1]


Any idea of how to approach this?

Here is an approach. In essence what you are looking for is a way to detach the geometry from the equations. You can do that with ElementMarker. Let me show how such an approach might look like. We generate a boundary mesh of the region. Note the second argument to the LineElement these are integer markes. Each outer edge gets an arbitrary marker; the inner edges are not relevant and we set them to 0. We also compute the center of mass (coms) for each triangle, we'll need those just now.

Needs["NDSolveFEM"]
Manipulate[
coords = {{0, 0}, {1, 0}, {1, 1}, {0, 1}, {x, y}};
bmesh = ToBoundaryMesh["Coordinates" -> coords,
"BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4,
1}, {1, 5}, {2, 5}, {3, 5}, {4, 5}}, {1, 2, 3, 4, 0, 0, 0,
0}]}];
coms = Total[coords[[#]]]/3 & /@ {{1, 2, 5}, {2, 3, 5}, {3, 4,
5}, {4, 1, 5}};
Show[Graphics[{Red, Point[coms]}],
bmesh["Wireframe"]]
, {x, 0.1, 0.9}, {y, 0.1, 0.9}]


Next, we generate the full mesh. The center of masses mark each region. In each of these sub-regions a different to the above edge marker is attributed. During the mesh generation process each triangle that is created gets attributed a marker depending in which sub-region it resides.

rm = Thread[List[coms, {1, 2, 3, 4}]]
(* {{{0.428, 0.136}, 1}, {{0.7613333333333332, 0.46933333333333327},
2}, {{0.428, 0.8026666666666666},
3}, {{0.09466666666666668, 0.46933333333333327}, 4}} *)


Generate the mesh from the boundary:

mesh = ToElementMesh[bmesh, "RegionMarker" -> rm];


Visualize the markers of the mesh elements:

mesh["Wireframe"[
"MeshElementStyle" -> {FaceForm[Blue], FaceForm[Yellow],
FaceForm[Red], FaceForm[Green]}]]


Investigate the boundary markers:

mesh["Wireframe"["MeshElement" -> "BoundaryElements",
"MeshElementMarkerStyle" -> Purple]]


One can go further and specify a refinement for each sub-region:

rm = Thread[
List[coms, {1, 2, 3, 4}, 0.1*{0.01, 0.02, 0.005, 0.001}]];
mesh = ToElementMesh[bmesh, "RegionMarker" -> rm];
mesh["Wireframe"[
"MeshElementStyle" -> {FaceForm[Blue], FaceForm[Yellow],
FaceForm[Red], FaceForm[Green]}]]


Now, for the actual PDE solving. Sigma is now a function of ElementMarker and can have different sigma values in the respective sub-regions.

sigma = Which[ElementMarker == 1, {{3, 0}, {0, 3}},
ElementMarker == 2, {{1, 0}, {0, 1}},
ElementMarker == 3, {{10, 0}, {0, 10}},
ElementMarker == 4, {{1, 0}, {0, 1}}];


With this the call to NDSolve is decoupled from geometry. This call to NDSolve will work for any geometry that provides the needed markers. The DirichletConditions use the LineElement markers and the PDE coefficients use the mesh element markers. So even though we have an edge marker 1 and a mesh element marker 1, these are distinct because they operate on different levels.

if = NDSolveValue[{Inactive[Div][
sigma.Inactive[Grad][u[x, y], {x, y}], {x, y}] == 0,
DirichletCondition[u[x, y] == 0, ElementMarker == 1],
DirichletCondition[u[x, y] == 1, ElementMarker == 3]},
u[x, y], {x, y} \[Element] mesh];


Visualize the result:

Show[
ContourPlot[if, {x, y} \[Element] mesh,
ColorFunction -> "TemperatureMap"],
bmesh["Wireframe"]]


More information about markers can be found in the documentation for ToBoundaryMesh and ToElementMesh in the options section under RegionMarker and in the ElementMesh generation tutorial

• Great, thank you. It looks like a solution for any division of the domain into parts. Mar 30, 2015 at 11:02