How to evaluate the PDE solution dependent on the `RegionMarkers"?

Question

I am solving a classical PDE on a current through a composite conductive material:

$$\operatorname{div}(\sigma \operatorname{grad}u)=0$$

where $u=u(x,y)$ is the electrical potential and $\sigma$ is the conductivity. The problem is a 2D one. The domain consists of several sub-domains with the internal boundaries, each sub-domain has its own conductivity, sigma.

What I want is to evaluate and analyze the heat production.

$$P=\sigma\,(\operatorname{grad}u)^2$$

Say, I need to plot it along certain lines, or evaluate its value in certain points.

Below for the sake of simplicity I only consider a rectangular domain with one vertical internal boundary in the middle.

The potential in this problem can be found as follows. The code below makes the boundary and element meshes, the element mesh containing the region markers 1 and 2 indicating domains with the different conductivities:

<< NDSolve`FEM`
bm = ToBoundaryMesh[
   "Coordinates" -> {{0., 0.}, {1., 0.}, {2., 0.}, {2., 1.}, {1., 
      1.}, {0., 1.}}, 
   "BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 
        5}, {5, 6}, {6, 1}, {5, 2}}]}];
em = ToElementMesh[bm, 
   "RegionMarker" -> {{{0.5, 0.5}, 1}, {{1.5, 0.5}, 2}}];

It looks like the following

    Show[{
  em["Wireframe"["MeshElementStyle" -> Directive[EdgeForm[Red], Thin],
     "MeshElementMarkerStyle" -> Blue]],
  bm["Wireframe"["MeshElementStyle" -> Black]]
  }, ImageSize -> 300]

Now I can solve the equation assuming sigma=1 in the first domain and sigma=2 in the second:

    sigma = If[ElementMarker == 1, 1, 2];
sl1 = NDSolveValue[{D[sigma*D[u[x, y], x], x] + 
     D[sigma*D[u[x, y], y], y] == 0, DirichletCondition[u[x, y] == 0, x == 0], 
   DirichletCondition[u[x, y] == 1, x == 2]}, u[x, y], {x, y} ∈ em]

and here is the solution:

 Plot3D[sl1, {x, y} ∈ em, ColorFunction -> "Rainbow", 
 PlotRange -> {0, 1}, ImageSize -> 300]

which looks as it is expected.

Everything was OK until this point. Now a difficulty comes. I need to evaluate the heat power, sigma*(D[u[x, y], x]^2 + D[u[x, y], y]^2). However, this can only be done on the mesh, since only the mesh "knows" the values of the region markers, but without the mesh Mathematica does not know, what are they. I can, for example, integrate the heat power on the mesh:

NIntegrate[
 Evaluate[sigma*(D[sl1, x]^2 + D[sl1, y]^2)], {x, y} ∈ em]

(* 0.667 *)

I cannot, however, say, evaluate its value in a point

sigma*(D[sl1, x]^2 + D[sl1, y]^2) /. {x -> 0.5, y -> 0.5}

(* 0.444 If[ElementMarker == 1, 1, 2] *)

since it is not a mesh-related operation.

I would like to stress, that in this example there is a trivial solution:

    If[x < 1, 1, 2]*(D[sl1, x]^2 + D[sl1, y]^2) /. {x -> 0.5, y -> 0.5}
(*  0.444  *)

It cannot be applied in my real case, since the internal boundaries are multiple and have complex shapes.

I tried to add the definitions of the currents into the list of equations:

 sl2 = NDSolveValue[{D[sigma*D[u[x, y], x], x] + 
     D[sigma*D[u[x, y], y], y] == 0, jx[x, y] == sigma*D[u[x, y], x], 
   jy[x, y] == sigma*D[u[x, y], y], 
   DirichletCondition[u[x, y] == 0, x == 0], 
   DirichletCondition[u[x, y] == 1, x == 2]}, {u[x, y], jx[x, y], 
   jy[x, y]}, {x, y} ∈ em]

which gives an evident warning: "No DirichletCondition or Robin-type NeumannValue was specified for {jx,jy}; the result is not unique up to a constant", and returned an incorrect solution for the potential:

Plot3D[sl2[[1]], {x, y} ∈ em, ColorFunction -> "Rainbow", 
 PlotRange -> {0, 1}, ImageSize -> 300]

My question Have you an idea of the work around, to evaluate and print the expression for the heat power?

Answer 1

Original post

Here's a workaround that came to mind:

MeshElementToRegion[el_] := RegionUnion@Flatten@Normal@
    GraphicsComplex[em["Coordinates"], MeshElementToGraphicsPrimitives[el]]

{reg1, reg2} = RegionMember@*MeshElementToRegion /@ 
   Flatten@MeshElementSplitByMarker@em["MeshElements"];

sigma[{x_, y_}] := Which[
  reg1[{x, y}], 1,
  reg2[{x, y}], 2
  ]

Now you can evaluate your expression:

sigma[{x, y}]*(D[sl1, x]^2 + D[sl1, y]^2) /. {x -> 0.5, y -> 0.5}

0.444444

I don't think MeshElementGraphicsPrimitives does anything to take higher order meshes into account, so to ensure correctness I would use "MeshOrder" -> 1 as an option to ToElementMesh.

Explanation

Your mesh em is made up of a single TriangleElement of the form TriangleElement[incidents, markers] (see em["MeshElements"]). Incidents are indices of coordinates from the coordinate list that you supplied to ToElementMesh. For example {1, 2, 3} is a list of incidents 1, 2, and 3 corresponding to the first three coordinates in em["Coordinates"]. TriangleElement[{{1,2,3}, {4,5,6}}, {1, 2}] encodes two triangles, the first one with the element marker 1 and the second one with the element marker 2.

Again, ToElementMesh creates a single TriangleElement. Since it contains both triangles with element marker 1 and triangles with element marker 2, the last argument will be mixed, e.g. {1, 1, 2, 1, 2, 2,...}. MeshElementSplitByMarker takes a TriangleElement and splits it by the element marker.

MeshElementSplitByMarker[TriangleElement[{{1,2,3}, {4,5,6}}, {1, 2}]]

{TriangleElement[{{1, 2, 3}}, {1}], TriangleElement[{{4, 5, 6}}, {2}]}

Now, the first TriangleElement only includes triangles with element marker 1, and the second TriangleElement only includes triangles with element marker 2.

MeshElementToGraphicsPrimitives converts TriangleElements (and other mesh elements) into graphics primitives, for triangles the primitive is a polygon.

MeshElementToGraphicsPrimitives[TriangleElement[{{1, 2, 3}}, {1}]]

Polygon[{{1, 2, 3}}]

But for this to be meaningful we need to replace the incidents with their real coordinates. Is there a function for graphics which this reminds us of? Yes, this is what GraphicsComplex does.

GraphicsComplex[em["Coordinates"], Polygon[{{1, 2, 3}}]]

can be drawn, using Graphics. But GraphicsComplex is not a legitimate region. Applying Normal makes the substitution of coordinates explicit:

Normal@GraphicsComplex[em["Coordinates"], Polygon[{{1, 2, 3}}]]

{{Polygon[{{0., 0.}, {1., 0.}, {2., 0.}}]}}

and such a polygon is a legitimate region. This means that we can now use RegionMember to see if a point belongs to that polygon. Since we have previously used MeshElementSplitByMarker, the polygons each sublist will all belong to the same region and we can use RegionUnion to get a region corresponding to a particular mesh element marker.