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["NDSolve`FEM`"];
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.
??
