Here is a trapezoid within which I am solving a simple Laplace equation with the Dirichlet boundary conditions set at its top and bottom:

coord = {{0, 0}, {1/3, 1}, {2/3, 1}, {1, 0}};
incidents = {{1, 2}, {2, 3}, {3, 4}, {4, 1}};
bm = ToBoundaryMesh["Coordinates" -> coord2, 
   "BoundaryElements" -> {LineElement[incidents]}];
mesh = ToElementMesh[bm];
sigma = 1;
nds = NDSolveValue[{Inactive[Div][
     sigma*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} ∈ mesh]

enter image description here

Now I would like to calculate the integral of the square of the result gradient: Evaluate[Grad[nds, {x, y}].Grad[nds, {x, y}]] over x lying within the mesh, and study the dependence of the result on y.

How can I do that?

The result should probably be a nested list of with the sublists structure {yi, resi}, where yiis the point coordinate along y, and resi is the integral value in this point.

What I do understand myself, is that I could divide the initial trapezoid into several parts, and then integrate within each such a part and, thus, make a list. I am curious, if there is a more regular approach?


With two minor changes to the code

bm = ToBoundaryMesh["Coordinates" -> coord, 
     "BoundaryElements" -> {LineElement[incidents]}];
nds = NDSolveValue[{Inactive[Div][sigma*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} ∈ mesh];

the first to correct a typo, and the second for convenience, we obtain

Plot3D[nds[x, y], {x, y} ∈ mesh]

enter image description here


Plot3D[Evaluate[Grad[nds[x, y], {x, y}].Grad[nds[x, y], {x, y}]], {x, y} ∈ 
    mesh, PlotRange -> All]

enter image description here

To obtain the integral described in the question, define

f[y0_] := Quiet@NIntegrate[Evaluate[Grad[nds[x, y], {x, y}].Grad[nds[x, y], {x, y}]] 
    /. y -> y0, {x, y0/3, 1 - y0/3}]

which gives, for instance,

(* 0.672611 *)
  • $\begingroup$ I see the typo you corrected (coord instead of coord2), but cannot find, what was the second correction. Please kindly let me know. $\endgroup$ – Alexei Boulbitch Aug 5 '15 at 8:18
  • $\begingroup$ I changed u[x. y], {x, y} ∈ mesh to u, {x, y} ∈ mesh at the end of NDSolveValue. $\endgroup$ – bbgodfrey Aug 5 '15 at 13:28
  • $\begingroup$ Thank you, I appreciate your above answer very much. $\endgroup$ – Alexei Boulbitch Aug 5 '15 at 13:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.