# Integration of a PDE solution over a custom-defined subdomain

Preamble: In my previous question I raised a problem of a custom-chosen integration of the solution of a PDE. In that case it was an integration along a line lying within the PDE domain. The solution given there is correct, but oversimplified, since I have oversimplified my original question. The problem, however, seems to me of general interest, since (1) it seems to push the limits of the present Mma capability and (2) it will be of a general use. Let me first formulate the question in a general form, and then go down to examples with the corresponding codes.

The question formulation in a general form: Provided one has a PDE defined in a complex (possibly multiply connected) domain, let us understand how one can integrate the result over a subdomain (that may or may not include holes).

Discussion of the question: Indeed, if one finds the PDE solution, it is not to admire it. One often needs to operate on it. Operation often means integration, evaluating derivatives, search for extrema etc. Even more, often (at least in my case) such an operation is necessary within a subdomain of the domain, or along a line lying within the domain, or so. Furthermore, all these operations can be worked around, if the domain is simple.

A) If one only needs to integrate the solution of the whole multiply connected domain, he can comfortably use the construct Integrate[expression[x, y], Element[{x, y}, mesh]] that takes care not to integrate over the holes. I assume, for instance, that the PDE in question is defined in a domain within the (x, y) plane.

B) If the domain is multiply connected, but simple enough, it may be divided into marked subdomains from the very beginning (see here: "Quad Meshes" and "Element Meshes with Subregions" and few others), and functions may be introduced that are nonzero only in the necessary subdomains and equal to the PDE solution there. This requires tedious work on the level of the mesh building, but might be a workaround.

C) If, however, the domain is complex, the task to build the mesh with the subdomains may become close to impossible. The example I give below.

So my question is related to the following problem, whether or not there are any built-in mechanisms enabling one to do something like Integrate[expression[x, y], Element[{x, y}, meshPart]], but to integrate over a part of the mesh, rather than over the whole, such that the possible holes are automatically taken care of. Or, alternatively, if there is a workaround for such a task.

This, of course, can be viewed in a more general perspective, not only in regard to integration. My focus on it is simply related to my personal tasks. I realize, of course, that such a mechanism/workaround may simply not exist.

Example: This is a trapezoid:

coord = {{0, 0}, {1/3, 1}, {2/3, 1}, {1, 0}};
trap = Polygon[coord];


This gives a list of $n$ disks with the radius $r$, such that they do not intersect the top and the bottom of the trapezoid:

distr[r_] :=
ProbabilityDistribution[(1 + y)/NIntegrate[(1 + y), {y, 2 r, 1 - 2 r}], {y, 2 r, 1 - 2 r}]
lst[r_, n_Integer] := Table[{RandomReal[{#/3, 1 - #/3}], #} &[RandomVariate[distr[r]]], {n}]


They are distributed linearly along y, but this is not important for our further discussion.

This gives the trapezoid with $n$ holes of the radius $r$:

fd[r_, n_Integer] := Fold[RegionDifference, trap, Disk[#, r] & /@ lst[r, n]];


This makes the mesh with $n$ holes of the radius $r$:

Needs["NDSolveFEM"];
mesh[r_, n_Integer] := ToElementMesh[fd[r, n], "MaxCellMeasure" -> r/20,
"MeshQualityGoal" -> Automatic, "ImproveBoundaryPosition" -> False, "MeshOrder" -> 1];


Like the mesh m1, for example:

m1 = mesh[0.05, 10]
m1["Wireframe"]


This solves the equation (for electric potential u, the Laplace equation with u equal to $1$ at the top and $0$ at the bottom):

nds1 = NDSolveValue[{Inactive[Div][Inactive[Grad][u[x, y], {x, y}], {x, y}] == 0,
DirichletCondition[u[x, y] == 0, y == 0],
DirichletCondition[u[x, y] == 1, y == 1]},
u[x, y], {x, y} \[Element] m1];


and this shows the solution for the potential:

 Show[{
ContourPlot[Evaluate[nds1], {x, y} \[Element] m1,
PlotRange -> {0, 1}, ColorFunction -> "TemperatureMap",
PlotLegends -> Automatic, ImageSize -> 300],
m1["Wireframe"]
}]


Now we come to the point: look at the two integrals below. The first integrates the solution over the mesh, the second over the whole trapezoid:

NIntegrate[Evaluate[nds1], {x, y} \[Element] m1]
NIntegrate[Evaluate[nds1], {y, 0, 1}, {x, y/3, 1 - y/3}, Method -> "LocalAdaptive"]

(*
0.197248
0.222601
*)


The first is smaller, since it takes care not to integrate over the holes, while the second does the integration there as well. Indeed, if we add the area of 10 holes to the first integral, we get:

0.197248 + Pi * 0.05^2 * 10

(*  0.275788  *)


which is pretty close to the value of the second one (this works, since u varies within the domain rather slowly).

It is quite clear that to do this correctly we should have explicitly excluded the holes from integration in the second integral. This is, however, difficult to automate, if the holes are generated randomly.

Nevertheless, the task of the integration over the whole domain is feasible. It is solved by the former approach: NIntegrate[Evaluate[nds1], {x, y} \[Element] m1].

Let us now divide the trapezoid into several parts by horizontal lines, and integrate u[x, y] within each part, such that the holes are excluded from the integration.

??

Reverting to the old strategy of using Boole seems efficient on the test case:

SeedRandom[0];         (* to give a reproducible result *)
m1 = mesh[0.05, 10];
m1["Wireframe"]

Clear[x, y, u];
nds1 = NDSolveValue[{Inactive[Div][Inactive[Grad][u[x, y], {x, y}], {x, y}] == 0,
DirichletCondition[u[x, y] == 0, y == 0],
DirichletCondition[u[x, y] == 1, y == 1]},
u[x, y], {x, y} ∈ m1];

Table[               (* slices of height 0.2 *)
With[{x1 = x1},
NIntegrate[nds1*Boole[x1 <= x <= x1 + 0.2], {x, y} ∈ m1]],
{x1, 0, 0.8, 0.2}]
Total@%              (* checks *)
NIntegrate[nds1, {x, y} ∈ m1]
(*
{0.00496725, 0.0451698, 0.0640033, 0.0425635, 0.00529073}
0.161995
0.161995
*)

• Could you please kindly comment, what does With[{x1 = x1},...  construct do? – Alexei Boulbitch Aug 10 '15 at 9:06
• @AlexeiBoulbitch It replaces every occurrence of x1 in the body with the value of x1. Since NIntegrate is HoldAll, the x1 would not automatically evaluate; I forget if I tried it without the With; it's possible (maybe likely) that NIntegrate evaluates the integrand before Table is done. See tutorial/IntroductionToDynamic#2125133640. – Michael E2 Aug 10 '15 at 10:01

One way to do it is to specify an ExtralopationHandler (see section on extrapolation) and have it return 0. for queries outside the domain. For example:

    nds1 = NDSolveValue[{Inactive[Div][
Inactive[Grad][u[x, y], {x, y}], {x, y}] == 0,
DirichletCondition[u[x, y] == 0, y == 0],
DirichletCondition[u[x, y] == 1, y == 1]},
u[x, y], {x, y} \[Element] m1, "ExtrapolationHandler" -> {(0. &),
"WarningMessage" -> False}];

NIntegrate[Evaluate[nds1], {x, y} \[Element] m1]
NIntegrate[Evaluate[nds1], {y, 0, 1}, {x, y/3, 1 - y/3},

What this does is that for every query point that is outside the domain, 0. is returned in stead of the extrapolated value, which is the default. It's a bit unfortunate that NIntegrate rejects Indeterminate like Overflow or Infinity. I would found "ExtrapolationHandler" -> {(Indeterminate &),... cleaner but the behavior of NIntegrate is not going to change; it's been that way too long.
The default behavior of Interpolation is to return an InterpolatingFunction that does extrapolation if queried outside of it's domain. For many interpolation scenarios this is a reasonable choice. For the FEM, however, this will give almost certainly a wrong solution outside the domain. The reason is we do not know anything about the region beyond it's boundary conditions. So the "ExtrapolationHandler" was introduced to allow to change the behavior of how InterpolatingFunction behave when queried outside the domain. My feeling is that Indeterminate is a good choice since it reflects the state of affairs: it is indeterminate what value should be given on extrapolation. Numerical 0. is a good choice in this specific case because 0 is invariant for integration - it does not alter the result. But arguing from a engineering perspective the value outside the region is not 0. We do not know what it is.
The first argument to the ExtrapolationHandler is a function that returns the value to be returned by the InterpolatingFunction in cases of extrapolation and it can be (if I recall correctly) any expression. The "WarningMessage" -> ... option takes True or False and controls if a message is issued on extrapolation or not.